Hi everyone,
I’ve been working on implementing a code to solve a reaction-diffusion system and then move the mesh. However, I’ve been encountering some strange results, and I’m almost certain that the issue lies in the definition of the initial conditions for the mesh. I’ve set up a mixed function space for the coordinates of the mesh and the curvature. I suspect that the problem may be in this part of the code.
V01, _ = ME2.sub(0).collapse()
V02, _ = ME2.sub(1).collapse()
x_exp = Expression(x, V01.element.interpolation_points())
H_exp = Expression(div(sign*n), V02.element.interpolation_points())
normal_expr = Expression(n_oriented, V01.element.interpolation_points())
normal_vec = Function(V01)
X0.sub(0).interpolate(x_exp)
X0.sub(1).interpolate(H_exp)
Specifically, I suspect the issue lies in X0.sub(0).interpolate(x_exp) . I’ve already attempted to interpolate the mesh coordinate into X0.sub(0) , but the mesh always seems to move in a strange manner.
What is the correct way to define this function in the initial step?
Here is my MWE
# +
import os
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import math
import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, mesh
from dolfinx.fem import Function, functionspace, Expression, form
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.io import XDMFFile, gmshio
from dolfinx.mesh import CellType, GhostMode
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, div,dot
# Save all logging to file
log.set_output_file("log.txt")
# -
try:
import gmsh # type: ignore
except ImportError:
import sys
print("This demo requires gmsh to be installed")
sys.exit(0)
def gmsh_circle(model: gmsh.model, name: str) -> gmsh.model:
"""Create a Gmsh model of a circle.
"""
model.add(name)
model.setCurrent(name)
gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.8)
circle = model.occ.addCircle(0, 0, 0, 1, tag=1)
# Synchronize OpenCascade representation with gmsh model
model.occ.synchronize()
# Add physical marker for cells. It is important to call this
# function after OpenCascade synchronization
model.add_physical_group(dim=1, tags=[circle])
# Generate the mesh
model.mesh.generate(dim=1)
model.mesh.setOrder(2)
return model
def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
"""Create a DOLFINx from a Gmsh model and output to file.
"""
msh, ct, ft = gmshio.model_to_mesh(model, comm, rank=0,gdim=2)
msh.name = name
ct.name = f"{msh.name}_cells"
ft.name = f"{msh.name}_facets"
with XDMFFile(msh.comm, filename, mode) as file:
msh.topology.create_connectivity(0, 1)
file.write_mesh(msh)
file.write_meshtags(ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")
file.write_meshtags(ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
# Create model
model = gmsh.model()
model = gmsh_circle(model, "circle")
model.setCurrent("circle")
create_mesh(MPI.COMM_SELF, model, "circle", f"out_gmsh/mesh_rankCircle_{MPI.COMM_WORLD.rank}.xdmf", "w")
with XDMFFile(MPI.COMM_WORLD, "out_gmsh/mesh_rankCircle_0.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh(name="circle")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
gmsh.finalize
P1 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)
ME = functionspace(mesh, mixed_element([P1, P1], gdim=2))
n = ufl.CellNormal(mesh)
x = ufl.SpatialCoordinate(mesh)
r = x/ufl.sqrt(ufl.dot(x, x))
sign = ufl.sign(ufl.dot(n, r))
n_oriented = sign*n
V1 = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,), gdim=2)
V2 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)
ME2 = functionspace(mesh, mixed_element([V1, V2],gdim=2))
# Trial and test functions of the space `ME` are now defined:
q, v = ufl.TestFunctions(ME)
Xi, phi = ufl.TestFunctions(ME2)
# ```{index} split functions
# ```
# +
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
# X = Function(ME2) # current solution
X0 = Function(ME2) # solution from previous converged step
X0.name = "X0"
(Xh, H) = ufl.TrialFunctions(ME2) # current solution
#(Xh0, H0) = ufl.TrialFunctions(ME2) # solution from previous converged step
# Split mixed functions
u1, u2 = ufl.split(u)
u10, u20 = ufl.split(u0)
# # Split mixed functions
# Xh, H = ufl.split(X)
Xh0, H0 = ufl.split(X0)
au = 0.9
bu = 1.1
av = 0.8
bv = 0.95
# Interpolate initial condition
u.sub(0).interpolate(lambda x: (bu-au)*np.random.rand(x.shape[1]) +au)
u.sub(1).interpolate(lambda x: (bv-av)*np.random.rand(x.shape[1]) +av)
u.x.scatter_forward()
V01, _ = ME2.sub(0).collapse()
V02, _ = ME2.sub(1).collapse()
x_exp = Expression(x, V01.element.interpolation_points())
H_exp = Expression(div(sign*n), V02.element.interpolation_points())
normal_expr = Expression(n_oriented, V01.element.interpolation_points())
normal_vec = Function(V01)
X0.sub(0).interpolate(x_exp)
X0.sub(1).interpolate(H_exp)
X0.x.scatter_forward()
print (f"X array {X0.sub(0).x.array[:]}")
# -
#X0.x.array[:] = uh.x.array
u_vec = ufl.as_vector([u1, u2])
kp = ufl.as_vector([1.5 , 0])
# Next, various model parameters are defined:
dt = 5.0e-05 # time step
# Weak statement of the equations
k = dt
d = 10
gamma = 29
a = 0.1
b = 0.9
tol = 1e-4
d1 = 1.0
d2 = d
ks = 25
kb = 2
F = ((u1 - u10) / k)*q*dx + d1*inner(grad(u1), grad(q))*dx\
-(gamma*(u1*u1*u2-u1+a))*q*dx \
+ ((u2 - u20) / k)*v*dx + d2*inner(grad(u2), grad(v))*dx\
-(gamma*(-u1*u1*u2+b))*v*dx
ax = dot(((Xh) / k),n_oriented)*phi*dx + kb*inner(grad(H),grad(phi))*dx - 1/2 *kb* H0**2 * H *phi*dx + ks*H*phi*dx\
+ dot(H*n_oriented,Xi)*dx - inner(grad(Xh),grad(Xi))*dx
#+ kb*inner(grad(H),grad(phi))*dx + ks*H*phi*dx\
Lx = dot(((Xh0) / k),n_oriented)*phi*dx + dot(kp,u_vec)*phi*dx
# Create nonlinear problem and Newton solver
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
# We can customize the linear solver
ksp = solver.krylov_solver
opts = PETSc.Options() # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()
X_h = Function(ME2)
X_h.name = "X_h"
problemX = LinearProblem(ax, Lx, u = X_h)
# Step in time
t = 0.0
# Output file
name = "mesh_at_t"+str(t)
mesh.name = name
# Reduce run time if on test (CI) server
if "CI" in os.environ.keys() or "GITHUB_ACTIONS" in os.environ.keys():
T = 3 * dt
else:
T = 5 * dt
# Get the sub-space for c and the corresponding dofs
normal_vec.interpolate(normal_expr)
V0, dofs = ME.sub(0).collapse()
W02, dofs2 = ME2.sub(0).collapse()
u1 = u.sub(0)
u2 = u.sub(1)
Xh0 = X0.sub(0)
H0 = X0.sub(1)
file1 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/outputu1.xdmf", "w")
file3 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/normal.xdmf", "w")
file4 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/H.xdmf", "w")
file1.write_mesh(mesh)
file3.write_mesh(mesh)
file4.write_mesh(mesh)
file1.write_function(u1, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file3.write_function(normal_vec, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file4.write_function(X0.sub(1), t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
norm_stop = 0.1e-6
fileL2 = open('L2.dat', 'w')
z = 0
inc = 0
impresion = 1
#print (f"MESH array {mesh.geometry.x[:,]}")
while (t < T):
t += dt
problemX.solve()
r = solver.solve(u)
inc += 1
print(f"Step {int(t/dt)}: num iterations: {r[0]}")
l2_norm = np.linalg.norm(u.x.array[dofs]-u0.x.array[dofs])/dt
#print (f"L2_norm = {l2_norm}")
fileL2.write(f"{t}\t{l2_norm}\n")
u0.x.array[:] = u.x.array
X0.x.array[:] = X_h.x.array
num=inc/impresion-1
# print (f"MESH array shape {mesh.geometry.x[:,].shape}")
# print (f"X array shape {X_h.x.array[dofs2].shape}")
# print (f"MESH array {mesh.geometry.x[:,]}")
# print (f"X array {X_h.x.array[dofs2]}")
#print (f"X array {X0.sub(0).x.array[:]}")
X_coord = X0.x.array[dofs2].reshape(len(mesh.geometry.x[:,1]),2)
DeltaX = mesh.geometry.x[:,[0,1]]-X_coord
#print (f"Delta X array {DeltaX}")
mesh.geometry.x[:,0] = X_coord[:,0]
mesh.geometry.x[:,1] = X_coord[:,1]
normal_vec.interpolate(normal_expr)
if (int(num)==z) or (int(num-z)==0):
us, vs = ufl.split(u)
# Save solution to file (VTK)
name = "mesh_at_t"+str(t)
mesh.name = name
file1.write_mesh(mesh)
file1.write_function(u1, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file3.write_mesh(mesh)
file3.write_function(normal_vec, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file4.write_mesh(mesh)
file4.write_function(X0.sub(0), t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
z += 1
# if l2_norm < norm_stop:
# print("l2_norm convergence")
# break
file1.close()
fileL2.close()
print(l2_norm)
#print(H.x.array)
Anyway, if you could find out if the error lies in another part, I would appreciate it,
Thank you