Hi everybody!
I am working on an ESFEM problem and have the following discretization:
I am encountering an issue specifically with $\boldsymbol{T}_h^m$, which is defined as:
This is my MWE
import os
import time
import numpy as np
from mpi4py import MPI # type: ignore
from petsc4py import PETSc # type: ignore
# Matplotlib setup
import matplotlib as mpl # type: ignore
mpl.use('Agg')
import matplotlib.pyplot as plt # type: ignore
# FEniCSx and UFL imports
import dolfinx # type: ignore
import basix.ufl # type: ignore
import ufl # type: ignore
from basix.ufl import element, mixed_element # type: ignore
from dolfinx import mesh, log, default_real_type # type: ignore
from dolfinx.fem import ( # type: ignore
Function, functionspace, Expression, form, Constant # type: ignore
)
from dolfinx.fem.petsc import ( # type: ignore
NonlinearProblem, assemble_matrix, assemble_vector, create_vector # type: ignore
)
from dolfinx.nls.petsc import NewtonSolver # type: ignore
from dolfinx.io import XDMFFile, gmshio # type: ignore
from ufl import dx, grad, inner, div, dot # type: ignore
try:
import gmsh # type: ignore
except ImportError:
import sys
print("This program requires gmsh to be installed")
sys.exit(0)
start = time.perf_counter()
gmsh.initialize()
h = 0.5e-2
R = 0.504
dt = 1e-3
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
rank = MPI.COMM_WORLD.rank
if mesh_comm.rank == gmsh_model_rank:
membrane = gmsh.model.occ.addDisk(0, 0, 0, R, R) # x, y, z, x-radius, y-radius
gmsh.model.occ.synchronize()
# Make membrane a physical surface for GMSH to recognise when generating the mesh
gdim = 1 # 2D Geometric Dimension of the Mesh
gmsh.model.addPhysicalGroup(gdim, [membrane], 0) # Dimension, Entity tag, Physical tag
# Generate 2D Mesh with uniform mesh size
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", h)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", h)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(2)
mesh, cell_markers, _ = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=2)
mesh.name = "initial_mesh"
with XDMFFile(MPI.COMM_WORLD, "circle.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(cell_markers, mesh.geometry)
# Finalize GMSH
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
gmsh.finalize()
P1 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)
P2 = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,), gdim=2)
ME = functionspace(mesh, mixed_element([P1, P1], gdim=2))
MEx = functionspace(mesh, mixed_element([P2, P1], gdim=2))
V = functionspace(mesh, P2)
Q = functionspace(mesh, P1)
(Xh, Hh) = ufl.TrialFunctions(MEx)
(Xi, phi) = ufl.TestFunctions(MEx)
coordinates = mesh.geometry.x[:, :2]
q, v = ufl.TestFunctions(ME)
# +
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
Xh0 = Function(V)
Xn = Function(V)
Hh0 = Function(Q)
Hn = Function(Q)
# Split mixed functions
u1, u2 = ufl.split(u)
u10, u20 = ufl.split(u0)
k1 = 0.1
k2 = 2.0
ui = k1+k2
vi = k2/(ui)**2
# Interpolate initial condition
u.sub(0).x.array[:] = ui
u.sub(1).interpolate(lambda x: vi + np.maximum(1e-3 * x[0], 0))
# Finalize initialization of solution variables
u.x.scatter_forward()
# Compute oriented normal vectors
n = ufl.CellNormal(mesh)
xc = ufl.SpatialCoordinate(mesh)
r = xc / ufl.sqrt(ufl.dot(xc, xc)) # Radial unit vector
sign = ufl.sign(ufl.dot(n, r)) # Adjust sign for orientation consistency
n_oriented = sign * n # Corrected normal vector
t_oriented = ufl.as_vector((-n_oriented[1], n_oriented[0]))
# Define coordinate expression for interpolation
def x_exp(x):
"""Expression for coordinate interpolation."""
return np.vstack((x[0], x[1]))
# Interpolate mean curvature field
H_expr = Expression(div(sign * n), Q.element.interpolation_points())
Hh0.interpolate(H_expr)
# Define tangential projection operator
def Tan(X,nu):
"""Compute tangential projection operator."""
return (X) / dt - nu*inner((X) / dt, nu)
# Define geometric function space and functions
V_geom = functionspace(mesh, P2)
X_geom = Function(V_geom)
# Initialize displacement field
Xh0.interpolate(x_exp)
X_geom.interpolate(Xh0)
# Define simulation parameters
kp = ufl.as_vector([-2,1]) # Reaction coefficient
k = dt # Time step size
d = 10 # Diffusion coefficient
gamma = 10 # Reaction rate constant
a, b = k1, k2 # Model parameters
tol = 1e-4 # Tolerance for solver
d1, d2 = 0.5, 50 # Diffusion coefficients
# Stiffness coefficients
ks, kb = 2, 2
# Define weak formulation of the system
F = (
((u1 - u10) / k) * q * dx
+ d1 * inner(grad(u1), grad(q)) * dx
- div(u1*Tan(X_geom-Xh0, n_oriented)) * q * dx
- (gamma * (u1**2 * u2 - u1 + a)) * q * dx
+ ((u2 - u20) / k) * v * dx
- div(u2*Tan(X_geom-Xh0, n_oriented)) * v * dx
+ d2 * inner(grad(u2), grad(v)) * dx
- (gamma * (-u1**2 * u2 + b)) * v * dx
)
problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2
The specific issue arises on the line: problem = NonlinearProblem(F, u)
The rest of the implementation is working properly, but the problem appears when including the term div(u1*Tan(X_geom-Xh0, n_oriented)) * q * dx
I would really appreciate any help