Hi everyone, I’m trying to simulate a simple system of coupled PDEs
\begin{align}
u_t &=u_{xx}+u\left(1-u-v\right) \\
v_t &=v_{xx}+v\left(1-u-v\right) \\
u_x(t, 0) &=0, \quad u(t, 1)=1 \\
v(t, 0) &=0, \quad v_x(t, 1)=0 \\
u(0, x) &=x^2 \\
v(0, x) &=x(x-2)
\end{align}
I modified the script after the Cahn–Hilliard example, but the solution seems to get u and v swapped.
The script I use is as follows. Any help is appreciated.
"""
Coupled time-dependent 1D convection-diffusion equation with constant coefficients on [0, 1]
u_t=u_{xx}+fu, v_t=v_{xx}+fv, where fu=u(1-u-v), fv=v(1-u-v)
u_x(t,0)=0, u(t,1)=1,
v(t,0)=0, v_x(t,1)=0
u(0,x)=x^2, v(0,x)=x^2-2x
"""
# %%
import matplotlib.pylab as plt
import numpy as np
from basix.ufl import element, mixed_element
from dolfinx import mesh
from dolfinx.fem import Function, dirichletbc, functionspace, locate_dofs_topological
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from ufl import TestFunctions, derivative, dx, grad, inner, split
# Define temporal parameters
t = 0
T = 1 # final time
num_steps = 100 # time steps
dt = T / num_steps
# Create mesh and define function space
# Define the function spaces
nx = 200 # number of intervals in space
domain = mesh.create_interval(MPI.COMM_WORLD, nx, [0, 1])
P1 = element("Lagrange", domain.basix_cell(), 1)
ME = functionspace(domain, mixed_element([P1, P1])) # Total space for all DOFs
# Trial and test functions of the space `ME` are now defined:
q, v = TestFunctions(ME)
u = Function(ME) # current solution
u_n = Function(ME) # solution from previous converged step
# Split mixed functions
u0, u1 = split(u)
u0_n, u1_n = split(u_n)
# Create initial condition
u.x.array[:] = 0
initial_condition_u0 = lambda x: x[0] ** 2
u.sub(0).interpolate(initial_condition_u0)
initial_condition_u1 = lambda x: x[0] * (x[0] - 2)
u.sub(1).interpolate(initial_condition_u1)
# Define boundary conditions
def left_boundary(x):
return np.isclose(x[0], 0)
def right_boundary(x):
return np.isclose(x[0], 1)
fdim = domain.topology.dim - 1
# left boundary
boundary_facets_left = mesh.locate_entities_boundary(domain, fdim, left_boundary)
boundary_dofs_left = locate_dofs_topological(ME.sub(0), fdim, boundary_facets_left)
# right boundary
boundary_facets_right = mesh.locate_entities_boundary(domain, fdim, right_boundary)
boundary_dofs_right = locate_dofs_topological(ME.sub(1), fdim, boundary_facets_right)
bc_u0_right = dirichletbc(PETSc.ScalarType(1), boundary_dofs_right, ME.sub(0))
bc_u1_left = dirichletbc(PETSc.ScalarType(0), boundary_dofs_left, ME.sub(1))
# dirichlet boundary conditions
bcs = [bc_u0_right, bc_u1_left]
# neumann boundary conditions are all zero, making no contributions
# Variational formulation of the system
f0 = u0 * (1 - u0 - u1)
f1 = u1 * (1 - u0 - u1)
# Define the nonlinear form
F = (
inner(u0 - u0_n, q) * dx
+ dt * inner(grad(u0), grad(q)) * dx
+ dt * f0 * q * dx
+ inner(u1 - u1_n, v) * dx
+ dt * inner(grad(u1), grad(v)) * dx
+ dt * f1 * v * dx
)
# Compute the Jacobian (derivative of F with respect to u)
J = derivative(F, u)
# Create the nonlinear problem
problem = NonlinearProblem(F, u, bcs=bcs, J=J)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6 # Set a relative tolerance for convergence
# Time-dependent output
# xdmf = XDMFFile(domain.comm, "result/test.xdmf", "w")
# xdmf.write_mesh(domain)
times = []
solutions_u0 = []
solutions_u1 = []
times.append(t)
solutions_u0.append(u.sub(0).collapse().x.array.copy())
solutions_u1.append(u.sub(1).collapse().x.array.copy())
# Updating the solution and right hand side per time step
for i in range(num_steps):
t += dt
# Solve linear problem
solver.solve(u)
u.x.scatter_forward()
u_n.x.array[:] = u.x.array
# Extract the solution components
u0_out, u1_out = u.sub(0).collapse(), u.sub(1).collapse()
# Save the solution every ten steps
if (i + 1) % 10 == 0:
times.append(t)
solutions_u0.append(u0_out.x.array.copy())
solutions_u1.append(u1_out.x.array.copy())
# xdmf.write_function(uh, t)
# Plotting the evolution of u and v at every tenth time step
x = domain.geometry.x[:, 0]
# Create two subplots for u and v
_, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 10), sharex=True)
# Plotting the evolution of u
for i, u0_data in enumerate(solutions_u0):
ax1.plot(x, u0_data, label=rf"$t={{{times[i]:.2f}}}$", alpha=0.7)
ax1.set_ylabel(r"$u$")
ax1.set_title("Evolution of u with respect to t")
ax1.legend()
ax1.grid(True)
# Plotting the evolution of v
for i, u1_data in enumerate(solutions_u1):
ax2.plot(x, u1_data, label=rf"$t={{{times[i]:.2f}}}$", alpha=0.7)
ax2.set_xlabel(r"$x$")
ax2.set_ylabel(r"$v$")
ax2.set_title("Evolution of v with respect to t")
ax2.legend()
ax2.grid(True)
# Show the plots
plt.tight_layout()
plt.show()
# Write solution to file
# xdmf.write_function(u, t)
# xdmf.close()
# %%