Hello,
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)
solver.setOperators(A)
DeltaA = fem.Function(DG)
z_n.x.set(1)
z_np1.x.set(1)
i=0
# Newton loop. Overkill, I know...
while i < 5:
# Assemble Jacobian and residual
with L.localForm() as loc_L:
loc_L.set(0)
A.zeroEntries()
fem.petsc.assemble_matrix(A, jacobian)
A.assemble()
fem.petsc.assemble_vector(L, residual)
L.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
# Scale residual by -1
L.scale(-1)
# Solve linearized problem
solver.solve(L, DeltaA.vector)
DeltaA.x.scatter_forward()
# 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:
break
# NOW SETUP AND SOLVE FULL PDE SYSTEM....
Thank you!