search

FEniCS / dolfinx

Next generation FEniCS problem solving environment
https://fenicsproject.org
GNU Lesser General Public License v3.0
727 stars 177 forks source link

Add tests of quadrature elements #2872

Open mscroggs opened 10 months ago

mscroggs commented 10 months ago

Describe new/missing feature

Quadrature elements are not thoroughly tested. It would be good to add more tests. Open to suggestions of test that could be added.

A first test is added in #2870.

Suggestion user interface

No response

bleyerj commented 9 months ago

Hi, on dolfinx:nightly, the following test works for scalar elements shape=() but fails for vector elements shape=(1,) or shape=(2,)

def test_vector_element(shape):
    msh = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

    dx_m = ufl.Measure(
        "dx",
        domain=msh,
        metadata={"quadrature_degree": 1, "quadrature_scheme": "default"},
    )

    Qe = quadrature_element(
        msh.topology.cell_name(), value_shape=shape, scheme="default", degree=1
    )
    Quad = fem.functionspace(msh, Qe)
    q_ = ufl.TestFunction(Quad)
    dq = ufl.TrialFunction(Quad)
    one = fem.Function(Quad)
    one.vector.set(1.0)
    mass_L_form = fem.form(ufl.inner(one, q_) * dx_m)
    mass_v = fem.assemble_vector(mass_L_form)
    mass_a_form = fem.form(ufl.inner(dq, q_) * dx_m)
    mass_A = fem.assemble_matrix(mass_a_form)
    assert np.isclose(sum(mass_v.array), 1.0)
    assert np.allclose(np.diag(mass_A.to_dense()), mass_v.array)

The error is:

    mass_L_form = fem.form(ufl.inner(one, q_) * dx_m)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 189, in form
    return _create_form(form)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 184, in _create_form
    return _form(form)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 141, in _form
    ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 56, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 204, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 199, in compile_forms
    raise e
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 190, in compile_forms
    impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 260, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/compiler.py", line 97, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, options)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 86, in analyze_ufl_objects
    form_data = tuple(_analyze_form(form, options) for form in forms)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 86, in <genexpr>
    form_data = tuple(_analyze_form(form, options) for form in forms)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 193, in _analyze_form
    assert np.allclose(e._points, custom_q[0])
AttributeError: '_BlockedElement' object has no attribute '_points'
mscroggs commented 9 months ago

Hi, on dolfinx:nightly, the following test works for scalar elements shape=() but fails for vector elements shape=(1,) or shape=(2,)

def test_vector_element(shape):
    msh = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

    dx_m = ufl.Measure(
        "dx",
        domain=msh,
        metadata={"quadrature_degree": 1, "quadrature_scheme": "default"},
    )

    Qe = quadrature_element(
        msh.topology.cell_name(), value_shape=shape, scheme="default", degree=1
    )
    Quad = fem.functionspace(msh, Qe)
    q_ = ufl.TestFunction(Quad)
    dq = ufl.TrialFunction(Quad)
    one = fem.Function(Quad)
    one.vector.set(1.0)
    mass_L_form = fem.form(ufl.inner(one, q_) * dx_m)
    mass_v = fem.assemble_vector(mass_L_form)
    mass_a_form = fem.form(ufl.inner(dq, q_) * dx_m)
    mass_A = fem.assemble_matrix(mass_a_form)
    assert np.isclose(sum(mass_v.array), 1.0)
    assert np.allclose(np.diag(mass_A.to_dense()), mass_v.array)

The error is:

    mass_L_form = fem.form(ufl.inner(one, q_) * dx_m)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 189, in form
    return _create_form(form)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 184, in _create_form
    return _form(form)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 141, in _form
    ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 56, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 204, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 199, in compile_forms
    raise e
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 190, in compile_forms
    impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 260, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/compiler.py", line 97, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, options)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 86, in analyze_ufl_objects
    form_data = tuple(_analyze_form(form, options) for form in forms)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 86, in <genexpr>
    form_data = tuple(_analyze_form(form, options) for form in forms)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 193, in _analyze_form
    assert np.allclose(e._points, custom_q[0])
AttributeError: '_BlockedElement' object has no attribute '_points'

I think I know how to fix this, PR to FFCx hopefully incoming soon. I've opened https://github.com/FEniCS/dolfinx/pull/2907 so we can add this test once I've got it working

bleyerj commented 9 months ago

Thanks @mscroggs, maybe you can avoid or change the last check 

assert np.allclose(np.diag(mass_A.to_dense()), mass_v.array)

as it uses a dense format, I did that very quickly