Update mesh over time

mesh
1

Hello everyone,

I have implemented a model for the solution of a PDE system where the mesh is changing over time. However, when I visualize the solution, the mesh appears that it is not changing, although when I print the coordinates, they have changed. Is this a mistake in how I am exporting the solution each time? Or is there an error in the solution?

Here is my code:

# # Schnackeberg model 1D
#
import os

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, plot, mesh
from dolfinx.fem import Function, functionspace, Expression, FunctionSpaceBase
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile, VTKFile
from dolfinx.mesh import CellType, create_unit_square, GhostMode
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, SpatialCoordinate


# Save all logging to file
log.set_output_file("log.txt")
# -

dt = 5.0e-04  # time step

# A unit interval

msh = mesh.create_unit_interval(comm=MPI.COMM_WORLD,
                            nx=50,
                            ghost_mode=GhostMode.shared_facet)
P1 = element("Lagrange", msh.basix_cell(), 1)
ME = functionspace(msh, mixed_element([P1, P1]))


# Trial and test functions of the space `ME` are now defined:

q, v = ufl.TestFunctions(ME)

# +
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step

# Split mixed functions
u1, u2 = ufl.split(u)
u10, u20 = ufl.split(u0)

au = 0.9
bu = 1.1
av = 0.8
bv = 0.95
# u_n = InitialConditions(au,bu,av,bv)

# 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()

#dt = 1e-3 # time step size
k = dt
d = 10
gamma = 250 
a = 0.1
b = 0.9
tol = 1e-4

d1 = 1.0
d2 = d

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 

# +
# 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 used inside the NewtonSolver by
# modifying the PETSc options
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()
# -

# +
# Output file
file1 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg1D/outputu1.xdmf", "w")
#file2 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg1D/outputu2.xdmf", "w")
file1.write_mesh(msh)
#file2.write_mesh(msh)

# Step in time
t = 0.0

#  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 = 1000 * dt

# Get the sub-space for c and the corresponding dofs in the mixed space
# vector
V0, dofs = ME.sub(0).collapse()

# Prepare viewer for plotting the solution during the computation

u1 = u.sub(0)
u2 = u.sub(1)

u0.x.array[:] = u.x.array

x_cord = msh.geometry.x[:,0]

x_2 = x_cord[-1]
x_1 = x_cord[0]

xi = 2*(x_cord - x_1)/(x_2-x_1)-1

v = 0.01
mov = xi**2

print(mov)

while (t < T):
    
    t += dt
    r = solver.solve(u)
    print(f"Step {int(t/dt)}: num iterations: {r[0]}")
    u0.x.array[:] = u.x.array

    file1.write_function(u1, t)
    #file2.write_function(u2, t)
    
    x_cord = msh.geometry.x[:,0]

    x_2 = x_cord[-1]
    x_1 = x_cord[0]

    xi = 2*(x_cord - x_1)/(x_2-x_1)-1

    v = 0.01
    mov = xi

    msh.geometry.x[:,0] = msh.geometry.x[:,0]+v*mov
 


file1.close()
#file2.close()

print(msh.geometry.x[:,0])

Thank you.

2

The XDMFFile does not store a new mesh geometry at each time-step if not explicitly told to do so (as this is very expensive).

Here is a minimal reproducible example:

from mpi4py import MPI
import dolfinx

mesh = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, 10, 10, dolfinx.cpp.mesh.CellType.triangle)
mesh.name = "initial_mesh"

xdmf = dolfinx.io.XDMFFile(mesh.comm, "output.xdmf", "w")
xdmf.write_mesh(mesh)


def f(x):
    return x[0]**2, x[1]**2


V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (2, )))
u = dolfinx.fem.Function(V)
u.interpolate(f)

xdmf.write_function(
    u, 0.0, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")

# Perturb mesh
u.x.array[:] -= 1
mesh.geometry.x[:, :] += 0.1
mesh.name = "perturbed_mesh"
xdmf.write_mesh(mesh)
xdmf.write_function(
    u, 1.0, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
xdmf.close()
2 Likes