Not able to evolve the solution with time for mixed elements for a nonlinear problem

Hello,
I am solving a mixed-element problem. In the first part of the code, I solve a Laplace problem with no source term, and pure Dirichlet boundary conditions. The solution works as expected and I do not get any errors in the solution. The code is as follows:

from dolfinx import default_scalar_type as dst
from dolfinx import mesh, fem, io, log
from dolfinx.io import VTXWriter
from dolfinx.fem import Function, Constant, locate_dofs_topological, functionspace, dirichletbc
from dolfinx.mesh import locate_entities_boundary
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from ufl import dot, grad, inner, dx, TestFunction, split
from basix.ufl import element, mixed_element
from pathlib import Path

import numpy as np

cache_dir = Path('cache_Fick')

DOLFIN_EPS = 3.0e-16
c_0        = 1.0
c_L        = 0.0

L          = 20
domain     = mesh.create_interval(MPI.COMM_WORLD, 50, [0.0, L])
P1         = element("Lagrange", domain.basix_cell(), 2)
P2         = element("Lagrange", domain.basix_cell(), 2)
P1P2       = mixed_element([P1, P2])
VQ         = functionspace(domain, P1P2)
fdim       = domain.topology.dim - 1
v1v2       = TestFunction(VQ)
v1, v2     = split(v1v2)
c1c2_n     = Function(VQ) 
c1_n, c2_n = split(c1c2_n)

def left_boundary(x):
    return np.isclose(x[0], 0.0)
left_boundary_facets = locate_entities_boundary(domain, fdim, left_boundary)
V10                  = VQ.sub(0)
V1                   = V10.collapse()[0]
c1left               = Function(V1)
c1left.x.array[:]    = c_0
left_boundary_dofs1  = locate_dofs_topological((V10, V1), fdim, left_boundary_facets)
bc_01                = dirichletbc(c1left, left_boundary_dofs1, V10)
V20                  = VQ.sub(1)
V2                   = V20.collapse()[0]
c2left               = Function(V2)
c2left.x.array[:]    = c_0
left_boundary_dofs2  = locate_dofs_topological((V20, V2), fdim, left_boundary_facets)
bc_02                = dirichletbc(c2left, left_boundary_dofs2, V20)

def right_boundary(x):
    return np.isclose(x[0], L)
right_boundary_facets = locate_entities_boundary(domain, fdim, right_boundary)
c1right               = Function(V1)
c1right.x.array[:]    = c_L
right_boundary_dofs1  = locate_dofs_topological((V10,V1), fdim, right_boundary_facets)
bc_L1                 = dirichletbc(c1right, right_boundary_dofs1, V10)
c2right               = Function(V2)
c2right.x.array[:]    = c_L
right_boundary_dofs2  = locate_dofs_topological((V20, V2), fdim, right_boundary_facets)
bc_L2                 = dirichletbc(c2right, right_boundary_dofs2, V20)

bc      = [bc_01, bc_02, bc_L1, bc_L2]
F       = - dot(grad(c1_n), grad(v1)) * dx - dot(grad(c2_n), grad(v2)) * dx
problem = NonlinearProblem(F, c1c2_n, bc, jit_options={'cache_dir': cache_dir})
solver  = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol                  = 1e-6
solver.max_it                = 100
solver.report                = True
solver.solve(c1c2_n) 

Next, I want to use this steady-state solution as initial condition to solve a transient problem. However, the boundary conditions for the transient problem are the same as the steady state problem, hence I use the same boundary conditions defined for the steady state problem. However, my residual is different, so I implement a different residual. In addition, I define a new function that stores the solution at every time step. The part of the code that implements these things are as follows:

Ti        = 0
Tf        = 200
num_steps = 200
dt        = (Tf-Ti)/num_steps

c1c2            = Function(VQ)  
c1c2.interpolate(c1c2_n)
c1c2.x.scatter_forward()
c1, c2          = split(c1c2)
dt              = Constant(domain, dst(dt))

F_tr         = -dot(grad(c1), grad(v1)) * dx - inner(1, v1) * dx - ((c1 - c1_n) / dt) * v1 * dx \
-dot(grad(c2), grad(v2)) * dx - inner(1, v2) * dx - ((c2 - c2_n) / dt) * v2 * dx
problem_tr   = NonlinearProblem(F_tr, c1c2, bc, jit_options={'cache_dir': cache_dir})
solver_tr    = NewtonSolver(MPI.COMM_WORLD, problem_tr)
solver_tr.convergence_criterion = "incremental"
solver_tr.rtol                  = 1e-6
solver_tr.max_it                = 100
solver_tr.report                = True

folder = Path("resultsDEMO")
folder.mkdir(exist_ok=True, parents=True)
vtx_c1 = VTXWriter(domain.comm, folder / "ins_c1.bp", c1c2.sub(0).collapse(), engine="BP4")
vtx_c2 = VTXWriter(domain.comm, folder / "ins_c2.bp", c1c2.sub(1).collapse(), engine="BP4")
vtx_c1.write(Ti)
vtx_c2.write(Ti)

Next, I implement the time loop here:

t = Ti
while t < Tf + DOLFIN_EPS:
    t += dst(dt)
    print(f'Time {t}')
    solver_tr.solve(c1c2)
    c1c2.x.scatter_forward()
    c1c2_n.interpolate(c1c2)
    c1c2_n.x.scatter_forward()
    vtx_c1.write(t)
    vtx_c2.write(t)

When I do this for the case with no mixed elements, it works fine. However, when I have mixed elements the solution does not update inside the time loop. I have seen several posts before regarding this topic, but none of the proposed solutions worked for me specifically. Can you specify what specific mistake I am making here.

I tried several approaches today to solve this problem, but none of them have worked so far. Some of the ideas I have tried are as follows:

  1. Use the split function inside the loop, although there was a suggestion against that in one of the past posts. That did not work.
  2. Use the syntax of prevSoln.interpolate(newSoln) in the time loop, as some people have done that. However, even that does not seem to work for this case.
  3. I also suspected that the correct weak form is not getting implemented in the time loop, so I changed the residual for the nonlinear solver to something nonsensical just to check. The nonlinear solver diverged, which makes me believe that I am solving the correct equation and plotting something wrong. I may be wrong. In the interim if anybody has some suggestions then that would be very welcome.

The issue is that you have collapsed the function within the initialization of the VTXWriter.
What you should do for efficiency is:

V, V_to_VQ = c1c2.sub(0).function_space.collapse()
Q, Q_to_VQ = c1c2.sub(1).function_space.collapse()
c1_out = dolfinx.fem.Function(V)
c2_out = dolfinx.fem.Function(Q)

c1_out.x.array[:] = c1c2.x.array[V_to_VQ]
c2_out.x.array[:] = c1c2.x.array[Q_to_VQ]

vtx_c1 = VTXWriter(domain.comm, folder / "ins_c1.bp", c1_out, engine="BP4")
vtx_c2 = VTXWriter(domain.comm, folder / "ins_c2.bp", c2_out, engine="BP4")

and then within the loop do:

t = Ti
while t < Tf + DOLFIN_EPS:
    t += dst(dt)
    print(f'Time {t}')
    solver_tr.solve(c1c2)
    c1c2.x.scatter_forward()
    c1c2_n.interpolate(c1c2)
    c1c2_n.x.scatter_forward()
    # Update output variables and write to disk
    c1_out.x.array[:] = c1c2.x.array[V_to_VQ]
    c2_out.x.array[:] = c1c2.x.array[Q_to_VQ]
    vtx_c1.write(t)
    vtx_c2.write(t)
1 Like