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.