Difference assembly dolfin vs dolfinx

dolfinx
1

Hello all!

I’m trying to update some legacy fenics code to fenicsx and I can’t fix the following problem:

When I try to assemble a mass matrix in a simple DG space, the resulting matrices coincide (more or less) for polynomial degrees 1 and 2, but they have different entries for polynomial degrees 3 and 4.

I create the corresponding spaces with

import dolfin as df

import scipy.sparse as scp
import numpy as np

import dolfinx as dfx
from mpi4py import MPI
import ufl

deg = 4
no_cells = 5

# mesh and function spaces
mesh_df = df.UnitIntervalMesh(no_cells)
mesh_dfx = dfx.mesh.create_unit_interval(MPI.COMM_SELF, no_cells)

V_df = df.FunctionSpace(mesh_df, 'DG', deg)
V_dfx = dfx.fem.FunctionSpace(mesh_dfx, ('DG', deg))

and the assembly of the matrices with

# Test and trial functions
u_df = df.TestFunction(V_df)
v_df = df.TrialFunction(V_df)

v_dfx = ufl.TestFunction(V_dfx)
u_dfx = ufl.TrialFunction(V_dfx)

# Setting up bilinear forms and assembling matrices
m_df_form = u_df*v_df*df.dx
M_df = df.PETScMatrix()
df.assemble(m_df_form, tensor=M_df)

m_dfx = u_dfx*v_dfx*ufl.dx
m_dfx_form = dfx.fem.form(m_dfx)
M_dfx = dfx.fem.petsc.assemble_matrix(m_dfx_form)
M_dfx.assemble()

Since in fenics I worked with scipy matrices I’d first like to get the program running with fenicsx and scipy matrices and switch to petsc later, so I converted the matrices to a csr format:

# convert PETSc to scipy.sparse csr matrix
indptr, indices, data = M_df.mat().getValuesCSR()
M_csr_df = scp.csr_matrix((data, indices, indptr))

indptr, indices, data = M_dfx.getValuesCSR()
M_csr_dfx = scp.csr_matrix((data, indices, indptr))

To compare the entries in the matrices (they have a different order inside a cell in dolfin and dolfinx for higher polynomial degrees), I cut off any values smaller than a tolerance, then sort them by value and subtract the results.

tol = 1e-14
df_vals = np.sort( M_csr_df.A[np.abs(M_csr_df.A) > tol] )
dfx_vals = np.sort( M_csr_dfx.A[np.abs(M_csr_dfx.A) > tol] )

diff = df_vals - dfx_vals
max_val = np.max(np.abs(diff))

For max_val I get the following results:
deg = 1 → max_val = 2.7755575615628914e-17
deg = 2 → max_val = 3.122502256758253e-17
deg = 3 → max_val = 0.006666666666666662
deg = 4 → max_val = 0.014814814814814788

I don’t know what I should do differently, so that it also works for polynomial degrees higher than 2.

My dolfinx version is 0.6.0 if that helps.

Thanks for your help!

2

This is because you are using a different lagrange variant in DOLFINx, see:
https://docs.fenicsproject.org/dolfinx/v0.6.0/python/demos/demo_lagrange_variants.html?highlight=ufl_wrapper
To illustrate this, consider the following one element problem

import dolfinx as dfx
from mpi4py import MPI
import ufl
import scipy.sparse as scp
import basix.ufl_wrapper
import basix

deg = 3
no_cells = 1

# mesh and function spaces
def compute_matrix(element):
    mesh_dfx = dfx.mesh.create_unit_interval(MPI.COMM_SELF, no_cells)
    ufl_element = basix.ufl_wrapper.BasixElement(element)

    V_dfx = dfx.fem.FunctionSpace(mesh_dfx, ufl_element)

    v_dfx = ufl.TestFunction(V_dfx)
    u_dfx = ufl.TrialFunction(V_dfx)

    # Setting up bilinear forms and assembling matrices

    m_dfx = u_dfx*v_dfx*ufl.dx
    m_dfx_form = dfx.fem.form(m_dfx)
    M_dfx = dfx.fem.petsc.assemble_matrix(m_dfx_form)
    M_dfx.assemble()

    indptr, indices, data = M_dfx.getValuesCSR()
    M_csr_dfx = scp.csr_matrix((data, indices, indptr))
    return M_csr_dfx

element = basix.create_element(basix.ElementFamily.P, basix.CellType.interval, deg,
                               basix.LagrangeVariant.gll_warped, True)

element2 = basix.create_element(basix.ElementFamily.P, basix.CellType.interval, deg,
                               basix.LagrangeVariant.equispaced, True)

M1 = compute_matrix(element)
M2 = compute_matrix(element2)
print(M1.todense())
print(M2.todense())

Giving:

[[ 0.07142857  0.01190476  0.02661986 -0.02661986]
 [ 0.01190476  0.07142857 -0.02661986  0.02661986]
 [ 0.02661986 -0.02661986  0.35714286  0.05952381]
 [-0.02661986  0.02661986  0.05952381  0.35714286]]
[[ 0.07619048  0.01130952  0.05892857 -0.02142857]
 [ 0.01130952  0.07619048 -0.02142857  0.05892857]
 [ 0.05892857 -0.02142857  0.38571429 -0.04821429]
 [-0.02142857  0.05892857 -0.04821429  0.38571429]]

while they are equal for deg=2, as there are only one dof, which makes gll_warped == equispaced

3 Likes
3

Is there a recommended lagrange variant? / Which Lagrange variant coincides with the standard one from fenics for DG elements?

4

As you can see in

GLL fixes Runge’s phenomenon that is experienced with equispaced elements, thus that is what is used under the hood in dolfinx as default:

Old dolfin uses equispaced elements.

1 Like
5

How can I create a a VectorFunctionSpace with equispaced nodes, as in your element2?

When I try

element = basix.create_element(basix.ElementFamily.P, basix.CellType.interval, deg, basix.LagrangeVariant.equispaced, True)
ufl_element = basix.ufl_wrapper.BasixElement(element)
V = dfx.fem.VectorFunctionSpace(mesh_dfx, ufl_element)
//
V = dfx.fem.VectorFunctionSpace(mesh_dfx, ufl_element, dim)

it doesn’t work. I don’t get an error message, but my program just keeps running and running…

6

Use: basix.ufl_wrapper — Basix 0.6.0 documentation
and dfx.fem.FunctionSpace(mesh_dfx, ufl_element)

1 Like