Hi everyone, I’m encountering an issue with incorrect temporal convergence order results while solving the Cahn-Hilliard-Navier-Stokes equations using FEniCS. I wonder if anyone can offer some guidance or share insights.
To provide more details, here’s the variational form (weak form) of my coupled Cahn-Hilliard-Navier-Stokes system:
P1 = FiniteElement(“Lagrange”, mesh.ufl_cell(), 1)
B = FiniteElement(“Bubble”, mesh.ufl_cell(), 3)
V_mini = VectorElement(NodalEnrichedElement(P1, B))
P1_p = FiniteElement(“Lagrange”, mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, MixedElement([P1, P1, P1, V_mini, P1_p]))
# Create MixedElement with all 5 spaces
u = Function(ME)
u0 = Function(ME)
ME = u.function_space()
mesh = u.function_space().mesh()
du = TrialFunction(ME)
q, v, w, r, s = TestFunctions(ME)
dc, dmu, dsig, dvel, dp = split(du)
c, mu, sig, vel, p = split(u)
c0, mu0, sig0, vel0, p0 = split(u0)
''' nonlinear functions '''
c = variable(c)
DPsi = (c**3-c0) # derivative of potential, dpsi/dphi
m_phi = 1.0 #0.5*(1+c0)**2 + 1e-12 # mobility function m(phi)
''' first equation, testing with q '''
L01 = c*q*dx - c0*q*dx
L02 = dt*inner(dot(grad(c0), vel),q)*dx
L03 = dt*m_phi*dot(grad(mu), grad(q))*dx
L0 = L01 + L02 + L03
''' second equation, testing with v '''
L11 = mu*v*dx - self.beta/self.epsilon*DPsi*v*dx
L12 = - self.beta*self.epsilon*dot(grad(c), grad(v))*dx
L13 = self.chi*sig*v*dx
L1 = L11 + L12 + L13
''' third equation, testing with w '''
L21 = dt * dot(grad(sig), grad(w)) * dx
if self.dynamic:
L21 += (sig * w * dx - sig0 * w * dx)
L22 = - dt * self.eta * dot(grad(c), grad(w)) * dx
L24 = dt * inner(dot(grad(sig0), vel), w) * dx
L2 = L21 + L22 + L24
''' Equation 4 (vel), testing with r (placeholder) '''
L31 = inner(vel - vel0, r) * dx
L32 = dt * self.eta1 * inner(nabla_grad(vel), nabla_grad(r)) * dx
L33 = 0.5 * dt * inner(dot(vel0, nabla_grad(vel)), r) * dx
L34 = - 0.5 * dt * inner(dot(vel0, nabla_grad(r)), vel) * dx
L35 = - dt * nabla_div(r) * p * dx
L36 = - dt * inner(dot(grad(c0), r), mu) * dx
L37 = - 50 * dt * inner(dot(grad(sig0), r), sig) * dx
L38 = dt * inner(dot(grad(sig0), r), c) * dx
L3 = L31 + L32 + L33 + L34 + L35 + L36 + L37 + L38
''' Equation 5 (p - Continuity Equation), testing with s '''
L4 = div(vel)*s*dx
''' Compute directional derivative about u in the direction of du (Jacobian) '''
L = L0 + L1 + L2 + L3 + L4
a = derivative(L, u, du)
''' Create nonlinear problem and Newton solver '''
problem = MyNonlinearProblem(a,L)
=======================================
— Summary and Temporal Convergence Order —
— Component: c —
L-infinity(L2) error (dt=1/256): 2.362466e-01
L-infinity(L2) error (dt=1/512): 1.565832e-01
Estimated temporal convergence order (p): 0.5934
— Component: mu —
L-infinity(L2) error (dt=1/256): 4.969451e-01
L-infinity(L2) error (dt=1/512): 3.096248e-01
Estimated temporal convergence order (p): 0.6826
— Component: sig —
L-infinity(L2) error (dt=1/256): 4.436670e-03
L-infinity(L2) error (dt=1/512): 2.669192e-03
Estimated temporal convergence order (p): 0.7331
— Component: vel —
L-infinity(L2) error (dt=1/256): 2.965885e-02
L-infinity(L2) error (dt=1/512): 1.777768e-02
Estimated temporal convergence order (p): 0.7384
— Component: p —
L-infinity(L2) error (dt=1/256): 3.088319e-01
L-infinity(L2) error (dt=1/512): 2.721607e-01
Estimated temporal convergence order (p): 0.1824
I’m not sure if the issue stems from my code (e.g., variational form implementation, time discretization scheme, or coupling strategy) or other numerical settings. I’d really appreciate any help or insights from those familiar with temporal convergence analysis for coupled systems in FEniCS! If needed, I can provide more details about my weak formulations, time-stepping loop, or solver configurations."