Discrepancy in fem.Expression interpolation between DOLFINx 0.9 and 0.11

General
1

Hi everyone,

I’ve recently upgraded my solver from DOLFINx 0.9 to 0.11 and noticed a small discrepancy in my non-linear solver’s convergence history. After investigating, the only difference I could find is in the L2 norm of my local time-step vector k, which is stored in a CG_1 function space. In the definition of the time step I compute the cell size h_cube (cell_volume^1/3) in DG_0 and then interpolate it into CG_1 using dolfinx.fem.Expression. While the DG_0 norm of h_cube is identical between both versions, the resulting CG_1 norm of k differs slightly:

=== DOLFINx 0.9 ===
norm h_cube : 3.5573675499613735
norm k      : 1.4503997478694577

=== DOLFINx 0.11 ===
norm h_cube : 3.5573675499613735
norm k      : 1.4521587135232983

Although the difference in k is small,I was wondering if that was normal?
Here is the code that allows you to reproduce this in version 0.11:

import dolfinx
import dolfinx.fem.petsc
import ufl
import numpy as np
from dolfinx.io.gmsh import read_from_msh
from mpi4py import MPI
from petsc4py import PETSc

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

# Load mesh
Meshdata = read_from_msh("/path/to/mesh.msh", comm, gdim=3)
mesh = Meshdata.mesh

DG0 = dolfinx.fem.functionspace(mesh, ("DG", 0))
CG1 = dolfinx.fem.functionspace(mesh, ("CG", 1))

# Compute cell volumes in DG0
cell_vol = dolfinx.fem.Function(DG0)
v_dg = ufl.TestFunction(DG0)
dolfinx.fem.petsc.assemble_vector(
    cell_vol.x.petsc_vec, 
    dolfinx.fem.form(v_dg * ufl.dx)
)
cell_vol.x.petsc_vec.ghostUpdate(
    addv=PETSc.InsertMode.ADD, 
    mode=PETSc.ScatterMode.REVERSE
)

# Compute h_cube = volume^(1/3) in DG0
h_cube = dolfinx.fem.Function(DG0)
h_cube.x.array[:] = cell_vol.x.array**(1/3)
h_cube.x.scatter_forward()

# Interpolate DG0 expression into CG1 space
k = dolfinx.fem.Function(CG1)
dt_expr = dolfinx.fem.Expression(h_cube, CG1.element.interpolation_points)
k.interpolate(dt_expr)
k.x.scatter_forward()

# PETSc L2 Norms
global_norm_h = h_cube.x.petsc_vec.norm(PETSc.NormType.NORM_2)
global_norm_k = k.x.petsc_vec.norm(PETSc.NormType.NORM_2)

if rank == 0:
    print(f"=== DOLFINx {dolfinx.__version__} ===")
    print(f"norm h_cube : {global_norm_h:.16f}")
    print(f"norm k      : {global_norm_k:.16f}")

Thank you in advance!

2

Hi @Nathan_Langlet,

the problematic line here is k.interpolate(dt_expr). The expression dt_expr is a DG_0 function evaluated at CG_1 DOFs. These are for this setup exactly the vertices of the element, in particular they lie on the element boundary where dt_expr is discontinuous. Therefore this is not a well-defined operation.

Different values are obtained, most likely (did not test) because one of the element orders changed between the versions and thus a different DG value is used in some of the DOFs.

What would you like to achieve with the k function? One possible way to produce a ‘CG_1 representative of h_\text{cube}’ would be by projection:

u = ufl.TrialFunction(CG1)
v = ufl.TestFunction(CG1)

problem = dolfinx.fem.petsc.LinearProblem(
    ufl.inner(u, v) * ufl.dx,
    ufl.inner(h_cube, v) * ufl.dx,
    petsc_options_prefix="hcube_projection_", # comment out for v0.9.0
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
    },
)

k = problem.solve()
k.x.scatter_forward()

this yields for different versions:

=== DOLFINx 0.9.0 ===
norm h_cube : 4.2627697484524711
norm k      : 2.0077298577662206

=== DOLFINx 0.10.0 ===
norm h_cube : 4.2627697484524711
norm k      : 2.0077298577662206

=== DOLFINx 0.11.0.post0 ===
norm h_cube : 4.2627697484524711
norm k      : 2.0077298577662206