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.