Tangential gradient issue

General
1

Hi everybody!

I am working on an ESFEM problem and have the following discretization:

image
image1048×168 21.7 KB

I am encountering an issue specifically with $\boldsymbol{T}_h^m$​, which is defined as:

image

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