How to use DG solution in CG variational form?


I have a challenging variational form to implement in fenics. My CG variational form depends on a term that is the discrete solution of an evolution equation at each time step, for example \int \frac{u_{n+1}-u_n}{\Delta t} z^h w, where z^h is the discrete solution of an ODE with element-wise values. To solve for z^h, because there is not a standard way to solve ODEs with fenics, I am solving the problem as a PDE with piecewise constant DG functions. And z^h MUST be solved with the same time integration method as the PDE.

The problem now is that I need z^h to be one value per element to multiply times my PDE. How can I do this? Is there a way to get the value at the quadrature point(s)? Perhaps there is an easier way that I do not know of…

I have put together a sample of the DG solve, but I do no know if the code is even helpful for this question:

import numpy as np
from dolfinx import mesh, fem, nls, plot, cpp
from mpi4py import MPI
from petsc4py import PETSc
import ufl

dt = 0.1
# create domains
nx, ny = 20,20
domainRef = mesh.create_rectangle(comm=MPI.COMM_WORLD, points=((0., 0.), (1., 1.)),n=(nx, ny), cell_type=mesh.CellType.triangle)
x_ref = domainRef.geometry.x

domainCur = mesh.create_rectangle(comm=MPI.COMM_WORLD, points=((0., 0.), (1., 1.)),n=(nx, ny), cell_type=mesh.CellType.triangle)
x_cur = domainCur.geometry.x

# corresponding function spaces
CG = fem.FunctionSpace(domainCur, ("CG",1))
CG_vec = fem.VectorFunctionSpace(domainRef, ("CG",1),dim=2)

delX = fem.Function(CG)
delX.interpolate(lambda x: 0.2*x[1])
x_cur[:,0] += delX.x.array[:]

# parameters to integrate
DG = fem.FunctionSpace(domainRef, ("DG",0))
w = ufl.TestFunction(DG)
z_np1 = fem.Function(DG)
z_n = fem.Function(DG)

# functions for ODE
vel = fem.Function(CG_vec)
vel.sub(0).interpolate(lambda x: (0.2*x[1])/dt)
vel.sub(1).interpolate(lambda x: (x[1] - x[1])/dt)
det = ufl.JacobianDeterminant(domainRef)

# ODE residual
Res = ufl.inner((z_np1 - z_n)/dt, w)*ufl.dx - ufl.inner(det*ufl.div(vel), w)*ufl.dx

# solver setup
J = ufl.derivative(Res, z_np1)
residual = fem.form(Res)
jacobian = fem.form(J)

# Iteration dependent structure
A = fem.petsc.create_matrix(jacobian)
L = fem.petsc.create_vector(residual)

solver = PETSc.KSP().create(domainRef.comm)
DeltaA = fem.Function(DG)

# Newton loop. Overkill, I know...
while i < 5:
    # Assemble Jacobian and residual
    with L.localForm() as loc_L:
    fem.petsc.assemble_matrix(A, jacobian)
    fem.petsc.assemble_vector(L, residual)
    L.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)

    # Scale residual by -1

    # Solve linearized problem
    solver.solve(L, DeltaA.vector)

    # Update du_{i+1} = du_i + delta x_i
    z_np1.x.array[:] += DeltaA.x.array

    i += 1
    # Compute norm of update
    correction_norm = DeltaA.vector.norm(0)
    print(f"Iteration {i}: Correction norm {correction_norm}")
    if correction_norm < 1e-10:


Thank you!

If you have solved the ODE for a given time step, say in variable z_np1 you can have that in your variational formulation by standard multiplication?

1./dt*inner((u_np1-u_n)*z_np1, v)*dx 

or am I missing something about your question?

The shapes do not match. In your example above,

1./dt*inner((u_np1-u_n)*z_np1, v)*dx

The function u_np1 would be a CG function, but after a quick follow on test:

CG = fem.FunctionSpace(domainRef, ("CG",1))
u = ufl.TrialFunction(CG)
v = ufl.TestFunction(CG)
Res = ufl.inner(u*z_np1, v)*ufl.dx

I see that there is no error with setting the residual in this way. Is there a place where I can learn more about the finite element assembly procedure in DOLFINx?

Please show the code that produced that error. If v and u are in scalar spaces it should work.

What do you want to know? There are plenty of demos,
and a tutorial at:

Garth Wells also gave a very good talk a few months back: FEM@LLNL | FEniCSx: Design of the Next Generation FEniCS Libraries for Finite Element Methods - YouTube

Thank you! I’ve found all the demos and tutorials great, and perhaps I’ve overlooked some helpful information within those. In this case, puzzeled why I can multiply a DG function by a CG function. They are defined on the same mesh, sure, but the number of DOFs don’t match. So I’m wondering how DOLFINx does the multiplication between non-matchin data structures. Is it element-by-element, and the multiplicaiton happens after the function values are taken to the quadrature points?

DOLFIN takes a variational form written in the unified form language, and generates integration kernels (Where things are evaluated at quadrature points, and summed accordingly).

If you want to understand the general framework, having a look at the talk by Garth that I linked above, or reading
(The code generation optimizations has changed with DOLFINx, but the general concept is similar).

Thank you! I really appreciate it!