firedrakeproject / firedrake

Firedrake is an automated system for the portable solution of partial differential equations using the finite element method (FEM)
https://firedrakeproject.org
Other
482 stars 157 forks source link

BUG: incorrect assembly of constant zero forms #3598

Open APaganini opened 4 weeks ago

APaganini commented 4 weeks ago

Describe the bug firedrake.assemble does not work correctly on zero forms with constant scalar values (and in one case it even throws an IndexError

Steps to Reproduce Run the following

from firedrake import *
from pyadjoint.adjfloat import AdjFloat
import numpy as np
mesh = UnitSquareMesh(3,3)
print("Expect 1, got ", assemble(1*dx(domain=mesh)))
print("Expect nan, got ", assemble(Constant(np.nan)*dx(domain=mesh)))
print("Expect nan, got ", assemble(np.nan*dx(domain=mesh)))
print("Expect nan, got ", assemble(AdjFloat(np.nan)*dx(domain=mesh)))
print("Expect 0, got ", assemble(Constant(0)*dx(domain=mesh)))
print("Expect 0, got ", assemble(0*dx(domain=mesh)))

Error message Expect 1, got 1.0000000000000002 Expect nan, got nan Expect nan, got 0.0 Expect nan, got 0.0 Expect 0, got 0.0 Traceback (most recent call last): File "/Users/admp1/Documents/bug_assemble.py", line 10, in print("Expect 0, got ", assemble(0dx(domain=mesh))) File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func File "/Users/admp1/Documents/FIREDRAKE/firedrake/src/firedrake/firedrake/adjoint_utils/assembly.py", line 30, in wrapper output = assemble(form, args, *kwargs) File "/Users/admp1/Documents/FIREDRAKE/firedrake/src/firedrake/firedrake/assemble.py", line 132, in assemble return get_assembler(expr, args, **kwargs).assemble(tensor=tensor) File "/Users/admp1/Documents/FIREDRAKE/firedrake/src/firedrake/firedrake/assemble.py", line 969, in assemble tensor = self.allocate() File "/Users/admp1/Documents/FIREDRAKE/firedrake/src/firedrake/firedrake/assemble.py", line 1115, in allocate comm=self._form.ufl_domains()[0]._comm IndexError: tuple index out of range

Environment:

Firedrake Configuration:
    package_manager: True
    minimal_petsc: False
    mpicc: None
    mpicxx: None
    mpif90: None
    mpiexec: None
    disable_ssh: True
    honour_petsc_dir: False
    with_parmetis: False
    slepc: False
    packages: []
    honour_pythonpath: False
    opencascade: False
    tinyasm: False
    torch: False
    petsc_int_type: int32
    cache_dir: /Users/admp1/Documents/FIREDRAKE/firedrake/.cache
    complex: False
    remove_build_files: False
    with_blas: download
    netgen: False
Additions:
    None
Environment:
    PYTHONPATH: None
    PETSC_ARCH: None
    PETSC_DIR: None
Status of components:
---------------------------------------------------------------------------
|Package             |Branch                        |Revision  |Modified  |
---------------------------------------------------------------------------
|FInAT               |master                        |f352325   |False     |
|PyOP2               |master                        |7bef38fa  |False     |
|fiat                |master                        |d0bea63   |False     |
|firedrake           |master                        |149f8fda6 |False     |
|h5py                |firedrake                     |6b512e5e  |False     |
|libsupermesh        |master                        |84becef   |False     |
|loopy               |main                          |967461ba  |False     |
|petsc               |firedrake                     |272f13d92c7|False     |
|pyadjoint           |master                        |908b636   |False     |
|pytest-mpi          |main                          |f2566a1   |False     |
|tsfc                |master                        |021589c   |False     |
|ufl                 |master                        |fbd288e6  |False     |
---------------------------------------------------------------------------
Ig-dolci commented 4 weeks ago

The error with assemble(0*dx(domain=mesh)) occurs because the ufl.algorithms.map_integrands does not map integrals for Zero as coded at this line. I am unsure if the correct approach is also to map the Zero integrand. @dham any thoughts?

connorjward commented 4 weeks ago

I don't necessarily think that this is a UFL problem. It looks like we basically just preprocess the form to give us an "empty" form so things like form.ufl_domains() gives us back an empty tuple, hence the IndexError. I think a reasonable fix would be to modify assemble.py so it can deal with empty forms.