Solve time dependent problem with Error :Integral is missing an integral domain

General
1

I try to solve a time dependent problem of this kind : u_t -\epsilon div(grad(u))+dot(b,grad(u))+cu=f
with u=0 on the boundary.I use the SUPG and Euler Backwards scheme to solve this problem.
I use grids with 20,40,80,160 on each direction for space discretization while I choose number of time steps as 5 from t=0 to 1.
With this number of time steps, I encounter the following error
This integral is missing an integration domain. ERROR:UFL:This integral is missing an integration domain. Traceback (most recent call last): File "/home/fenics/FSI_Solver/SUPG_TEST.py", line 56, in <module> a=fem.form(u*v*dx+dt*epsilon*inner(grad(u),grad(v))*dx+dt*dot(b,grad(u))*v*dx+dt*c*u*v*dx+dt*delta*(-epsilon*div(grad(u))+dot(b,grad(u))+c*u)*(dot(b,grad(v)))*dx+(u*delta*dot(b,grad(v))*dx)) File "/usr/local/lib/python3.10/dist-packages/ufl/measure.py", line 413, in __rmul__ error("This integral is missing an integration domain.") File "/usr/local/lib/python3.10/dist-packages/ufl/log.py", line 135, in error raise self._exception_type(self._format_raw(*message)) ufl.log.UFLException: This integral is missing an integration domain.
Is this because the time step is so small?Since when I use 10 as number of time steps, it works well but it breaks down with 5,6,7,8.
The complete code is as follows:

from dolfinx import mesh,fem
import numpy as np
from ufl import div,grad,diff,dot,inner,dx,ds,CellDiameter,TrialFunction,TestFunction,sin,pi,exp
import ufl
from mpi4py import MPI
from petsc4py import PETSc

t=0
T=1


for nx in [20, 40, 80, 160,320]:
    
    ny=nx
    num_steps=8
    dt=(T-t)/num_steps
    domain=mesh.create_unit_square(MPI.COMM_WORLD,nx,ny,mesh.CellType.triangle)
    V=fem.FunctionSpace(domain,("CG",1))

    epsilon=fem.Constant(domain,PETSc.ScalarType(1e-8))
    b=fem.Constant(domain,np.array([1,-1],dtype=PETSc.ScalarType))
    c=fem.Constant(domain,PETSc.ScalarType(1))

    #define exact solution
    def u_exact(mod, t):
        return lambda x: mod.exp(mod.sin(2*mod.pi*t))*mod.sin(2*mod.pi*x[0])*mod.sin(2*mod.pi*x[1])


    x=ufl.SpatialCoordinate(domain)
    t_var=ufl.variable(t)

    exact_soln_ufl=u_exact(ufl,t_var)  #used for linear form
    exact_soln_np=u_exact(np,t)    #used for interpolation


    #Apply homogeneous Dirichlet boundary condition
    boundary_fun=fem.Function(V)
    boundary_fun.interpolate(exact_soln_np)

    tdim=domain.topology.dim
    fdim=tdim-1
    domain.topology.create_connectivity(fdim,tdim)
    boundary_facets=mesh.exterior_facet_indices(domain.topology)
    bc=fem.dirichletbc(boundary_fun,fem.locate_dofs_topological(V,fdim,boundary_facets))

    #Initial condition
    u_old=fem.Function(V)
    u_old.interpolate(exact_soln_np)

    #Define right hand side f
    f=diff(exact_soln_ufl(x),t_var)-epsilon*div(grad(exact_soln_ufl(x)))+dot(b,grad(exact_soln_ufl(x)))+c*exact_soln_ufl(x)

    u,v=TrialFunction(V),TestFunction(V)

    h=CellDiameter(domain)

    delta=h/4

    a=fem.form(u*v*dx+dt*epsilon*inner(grad(u),grad(v))*dx+dt*dot(b,grad(u))*v*dx+dt*c*u*v*dx+dt*delta*(-epsilon*div(grad(u))+dot(b,grad(u))+c*u)*(dot(b,grad(v)))*dx+(u*delta*dot(b,grad(v))*dx))
    L=fem.form(u_old*v*dx+dt*f*(v+delta*dot(b,grad(v)))*dx+u_old*delta*dot(b,grad(v))*dx)

    A=fem.petsc.assemble_matrix(a,bcs=[bc])
    A.assemble()
    b=fem.petsc.create_vector(L)
    uh=fem.Function(V)

    solver=PETSc.KSP().create(domain.comm)
    solver.setOperators(A)
    solver.setType(PETSc.KSP.Type.PREONLY)
    solver.getPC().setType(PETSc.PC.Type.LU)

    V_ex = fem.FunctionSpace(domain, ("CG", 2))
    u_ex = fem.Function(V_ex)


#Use Euler Backwards for time discretization
    for n in range(num_steps):
        t+=dt
        exact_soln_np=u_exact(np,t)
        boundary_fun.interpolate(exact_soln_np)

        with b.localForm() as loc_b:
            loc_b.set(0)
        fem.petsc.assemble_vector(b,L)

        fem.petsc.apply_lifting(b,[a],[[bc]])
        b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,mode=PETSc.ScatterMode.REVERSE)
        fem.petsc.set_bc(b,[bc])

        solver.solve(b,uh.vector)
        uh.x.scatter_forward()

        u_old.x.array[:] = uh.x.array

 # Compute L2 norm of the discrete solution at final time
L2 = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(fem.form((uh)**2 * dx)),op=MPI.SUM))
if domain.comm.rank == 0:
        print(f"L2: {L2:.2e}")

Thanks in advance!

2

I would start by trying:

from ufl import Measure
dx = ufl.Measure("dx", domain=domain)

and see if that works.

Note that if you wrap dt as a k = dolfinx.fem.Constant(mesh, dt) and use k in your variational form, it should give you some performance benefit, as the form does not have to be recompiled if you ever need to change T, t, or num_steps

4

Thank you very much.
Sorry for another question. Since I define f as
f=diff(exact_soln_ufl(x),t_var)-epsilon*div(grad(exact_soln_ufl(x)))+dot(b,grad(exact_soln_ufl(x)))+c*exact_soln_ufl(x)
it depends on time t . So when I solve the linear system at every time step, I need to update the time.If I try to use

 for n in range(num_steps):
        t+=dt

to update the time, will the f(t) be automatically updated to f(t+dt)?
Further, I use t_var=ufl.variable(t) to define the time t as a variable for differentiating with time.So if t is a variable now or we create a new variable t_var which has the value of t?I dont know how this works when defining bilinear form and when I compute the right hand side.

5

You should wrap t as a dolfinx.fem.Constant to make sure that it is automatically update when you call
t.value += dt
When you call dolfinx.fem.form you generate C code for assembly, and the only parameters that are modifiable after this are dolfinx.fem.Constant, dolfinx.fem.Function and the mesh geometry.

6

Thank you.
Now I define the paarameters and function f as follows

domain=mesh.create_unit_square(MPI.COMM_WORLD,nx,ny,mesh.CellType.triangle)
V=fem.FunctionSpace(domain,("CG",1))

k=T/num_steps
t=fem.Constant(domain,PETSc.ScalarType(0))
dt=fem.Constant(domain,PETSc.ScalarType(k))
def u_exact(mod, t):
            return lambda x: mod.exp(mod.sin(2*mod.pi*t))*mod.sin(2*mod.pi*x[0])*mod.sin(2*mod.pi*x[1])
x=ufl.SpatialCoordinate(domain)
t_var=ufl.variable(t)
exact_soln_ufl=u_exact(ufl,t_var)  #used for linear form
exact_soln_np=u_exact(np,t)    #used for interpolation
f=diff(exact_soln_ufl(x),t_var)-epsilon*div(grad(exact_soln_ufl(x)))+dot(b,grad(exact_soln_ufl(x)))+c*exact_soln_ufl(x)

and in iteration I use

for n in range(num_steps):
        t.value+=k

Will it automatically update the function f? And t_var is related with t so they change values in the same time as t changes?

7

Yes, as you are not changing the object t, but the internal value it represents, t_var is keeping track of the change.

8

Very helpful!Thank you very much.

9

Hi @dokken, I’ve just hit the error reported here, “Integral is missing an integral domain”. Your suggestion to explicitly define the dx measure resolves the problem.

What is strange is that my script runs fine with exactly the same code (without Measure), only with a different value for one of my parameters. Not sure I can reduce to a short MWE to demonstrate, but on the face of it the parameter that is different should be completely unrelated to the domain or the dx Measure.

Is it possible to provide some background on the solution given here? Why would the default ufl.dx object work fine with one value of the parameter, but fail, requiring an explicit Measure, for another value of the parameter, where the code is otherwise the same?

10

UFL consists several important classes, here I’ll refer to two of them, namely ufl.core.expr.Expr and ufl.form.Form.

A ufl.core.expr.Expr is anything that can be an integrand in ufl. This means it could be 1, u*v, ufl.inner(ufl.div(ufl.grad(u)), v) etc.

These expressions are usually associated with a mesh through its test or trial functions, coefficients or constants.

A form is whenever you want to integrate this expression over a domain, i.e. right multiply by ufl.dx or ufl.Measure("dx").
For this right multiplication to work (and create a ufl.integral.Integral, it requires an integration domain (a mesh, or ufl_domain in ufl-speak): https://github.com/FEniCS/ufl/blob/77ae57c4b9594851f6bff4a63b7eecb81477c486/ufl/measure.py#L380-L388

Some forms do not have this property, for instance
1*ufl.dx
This can be solved by either making 1 a dolfinx.fem.Constant(mesh, 1) (which is also good for other reasons), or attaching the ufl_domain to the measure.

However, ufl does some clever graph simplification and matching of terms, so that whenever you multiply by 0, it reduced the ufl-form:
i.e.

v = ufl.TestFunction(V)
L = v * ufl.dx

is fine, while

0*v*ufl.dx

produces the error with missing integration domain, as

0 * v

produces

Zero((), (), ())

which does not have a ufl domain.

I guess this is what is happening for on of the terms of your variational form.

The best way to counteract this is to wrap material parameters as dolfinx.fem.Constant (which is also beneficial for code generation), or simply always define the integration domain in the ufl.Measure, as it is the first place it looks: 77ae57c4b9594851f6bff4a63b7eecb81477c486/ufl/measure.py#L380-L388

without an example I cannot really pinpoint what goes wrong in your code. This is just my best guess from my understanding of UFL (which might not be 100 % correct, so anyone feel free to jump in and add more context or corrections).

1 Like
11

That makes sense about how const*ufl.dx is handled, it’s it’s certainly true I need to do through my code and replace some with dolfinx.fem.Constant. It’s strange though, in my instance the parameter that triggers the error is not a constant (not even 0), it’s a python function used as an argument inside other ufl maths functions (counterpart to the old Expression from legacy dolfin). In lieu of an available MWE maybe I should send you the code and make you a coauthor :slight_smile:

Setting an explicit Measure is easy enough to do in any case.