Issue fem.create_form with DG0 function and ufl.derivative

dolfinx
1

I am unsure whether this is a bug or I have a mistake. Consider this MWE:

from dolfinx import fem, mesh
import ufl
from mpi4py import MPI
import numpy as np

def main():
    els_side = 4
    domain  = mesh.create_unit_square(MPI.COMM_WORLD, els_side, els_side)
    V   = fem.functionspace(domain, ("Lagrange", 1),)
    DG0 = fem.functionspace(domain, ("Discontinuous Lagrange", 0))
    solution = fem.Function(V)
    dg0_func = True
    ufl_derivative = True
    if dg0_func:
        conductivity = fem.Function(DG0)
        conductivity.x.array[:] = 1.0
    else:
        conductivity = fem.Constant(domain, 1.0)
    v = ufl.TestFunction(V)
    dx = ufl.dx
    if ufl_derivative:
        a_ufl = conductivity * ufl.dot(ufl.grad(solution), ufl.grad(v)) * dx
        form = ufl.derivative(a_ufl, solution)
    else:
        u = ufl.TrialFunction(V)
        form = conductivity * ufl.dot(ufl.grad(u), ufl.grad(v)) * dx
    j_compiled = fem.compile_form(
            domain.comm, form, form_compiler_options={"scalar_type": np.float64})
    coefficient_map = {}
    for coeff in form.coefficients():
        coefficient_map[coeff] = coeff
    constant_map = {}
    for const in form.constants():
        constant_map[const] = const
    j_instance = fem.create_form(j_compiled,
                                 [V, V],
                                 msh=domain,
                                 subdomains={},
                                 coefficient_map = coefficient_map,
                                 constant_map = constant_map)
    A = fem.assemble_matrix(j_instance)
    print(A.data)

if __name__=="__main__":
    main()

When using ufl.derivative together with a DG0 function in the form, the matrix assembled from compiling and creating the bilinear form in two steps is wrong. Here I expect to see a stiffness matrix but I obtain a null matrix.

In this MWE, both the flags ufl_derivative and dg0_func are set to True. If either of them is False, I obtain the expected stiffness matrix. Moreover, if both of them are True, the coefficient_map argument of create_form doesn’t even need the conductivity to run, when it normally complains about missing constants or coefficients.

2

The issue is the coefficient map. It seems to me that this must map to the fem.function, not back to the ufl.coefficient. This runs:

coefficient_map = {form.coefficients()[0]:conductivity}

There are two coefficients in your form; the one associated with conductivity and the one associated to solution. Conductivity happens to be the first, hence the [0]. Not sure how to do that procedurally…

1 Like
3

There is a bug in the packing of constants and coefficients in DOLFINx, see:

However, I don’t quite understand your need for compile_form and create_form in the example above, as the problem is not domain independent (similarly constants and coefficients are defined on the actual grid, not in a symbolical fashion).

Anyhow, with the proposed PR above, your code should run correctly.

2 Likes
4

@Stein many thanks for having a look. form.coefficients()[0] is actually solution so the coefficient map you propose is {solution : conductivity}, which works but doesn’t make sense.

@dokken thanks a lot as always. The example above is a MWE replicating the bug. I want to run this with 3D printing and domain decomposition examples which rely on changing the subdomain_data for both boundary and domain integrals. If I can get this to work it will avoid me recompiling my forms.

5

Ah Ok, I got weirdly lucky then… I now wonder why my modification gave the expected result :sweat_smile: