Hi, so I am trying to create some simulation data for a Cahn-Hilliard system with periodic boundary conditions using fenicsx. I used the code from Cahn-Hilliard equation — DOLFINx 0.10.0.0 documentation to get started and have adjusted it for the boundary conditions.
My problem is that my boundary conditions don’t appear to be implemented since the simulation is not periodic. Here is my code, I will highlight the region of interest and if anyone with more knowledge than myself could point me in the right direction, it would be much appreciated.
import os
from petsc4py import PETSc
import numpy as np
import matplotlib.pyplot as plt
import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, plot, default_scalar_type
from dolfinx import mesh, fem, log
from dolfinx.fem import Function, functionspace, dirichletbc, locate_dofs_geometrical
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_unit_square, locate_entities_boundary
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.common import Timer
from dolfinx_mpc import MultiPointConstraint
from ufl import dx, grad, inner
from mpi4py import MPI
OFF_SCREEN = False
try:
import pyvista as pv
import pyvistaqt as pvqt
have_pyvista = True
if OFF_SCREEN:
pv.start_xvfb(wait=0.5)
except ModuleNotFoundError:
print("pyvista and pyvistaqt are required to visualise the solution")
have_pyvista = False
# Save all logging to file
log.set_output_file("log.txt")
# Next, various model parameters are defined:
lmbda = 0.0025 # surface parameter
dt = 0.0001 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicholson
for seed in range(11):
# A unit square mesh with 96 cells edges in each direction is created,
# and on this mesh a
# {py:class}`FunctionSpace <dolfinx.fem.FunctionSpace>` `ME` is built
# using a pair of linear Lagrange elements.
**# I THINK THE PROBLEM LIES BETWEEN HERE...**
msh = create_unit_square(MPI.COMM_WORLD, 128, 128, CellType.triangle)
P1 = element("Lagrange", msh.basix_cell(), 1)
ME = functionspace(msh, mixed_element([P1, P1]))
# Define periodic boundary conditions
tol = 250 * np.finfo(default_scalar_type).resolution
bcs = []
def periodic_boundary_x(x):
return np.isclose(x[0], 1, tol) & ~np.isclose(x[1], 0, tol) & ~np.isclose(x[1], 1, tol)
def periodic_relation_x(x):
out_x = np.zeros_like(x)
out_x[0] = x[0] - 1
out_x[1] = x[1]
return out_x
def periodic_boundary_y(x):
return np.isclose(x[1], 1, tol) & ~np.isclose(x[0], 0, tol) & ~np.isclose(x[0], 1, tol)
def periodic_relation_y(x):
out_x = np.zeros_like(x)
out_x[0] = x[0]
out_x[1] = x[1] - 1
return out_x
with Timer("~PERIODIC: Initialize MPC"):
mpc = MultiPointConstraint(ME)
for i in range(ME.num_sub_spaces):
subspace = ME.sub(i)
mpc.create_periodic_constraint_geometrical(subspace,periodic_boundary_x, periodic_relation_x, bcs)
mpc.create_periodic_constraint_geometrical(subspace, periodic_boundary_y, periodic_relation_y, bcs)
mpc.finalize()
# Trial and test functions of the space `ME` are now defined:
q, v = ufl.TestFunctions(mpc.function_space)
u = Function(mpc.function_space) # current solution
u0 = Function(mpc.function_space) # solution from previous converged step
**# ...AND HERE**
# Split mixed functions
c, mu = ufl.split(u)
c0, mu0 = ufl.split(u0)
# -
# Zero u
u.x.array[:] = 0.0
# Seed of random number generator
rng = np.random.default_rng(seed)
# Interpolate initial condition
u.sub(0).interpolate(lambda x: 0.6 + 0.001*(rng.random(x.shape[1]) - 0.5))
u.x.scatter_forward()
# -
# The first line creates an object of type `InitialConditions`. The
# following two lines make `u` and `u0` interpolants of `u_init` (since
# `u` and `u0` are finite element functions, they may not be able to
# represent a given function exactly, but the function can be
# approximated by interpolating it in a finite element space).
#
# ```{index} automatic differentiation
# ```
#
# The chemical potential $df/dc$ is computed using UFL automatic
# differentiation:
# Compute the chemical potential df/dc
c = ufl.variable(c)
f = c**2 * (1-c)**2
dfdc = ufl.diff(f, c)
# The first line declares that `c` is a variable that some function can
# be differentiated with respect to. The next line is the function $f$
# defined in the problem statement, and the third line performs the
# differentiation of `f` with respect to the variable `c`.
#
# It is convenient to introduce an expression for $\mu_{n+\theta}$:
# mu_(n+theta)
mu_mid = (1.0 - theta) * mu0 + theta * mu
# which is then used in the definition of the variational forms:
# Weak statement of the equations
F0 = inner(c, q) * dx - inner(c0, q) * dx + dt * inner(grad(mu_mid), grad(q)) * dx
F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
F = F0 + F1
# This is a statement of the time-discrete equations presented as part
# of the problem statement, using UFL syntax.
#
# ```{index} single: Newton solver; (in Cahn-Hilliard demo)
# ```
#
# The DOLFINx Newton solver requires a
# {py:class}`NonlinearProblem<dolfinx.fem.NonlinearProblem>` object to
# solve a system of nonlinear equations
# +
# 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"
sys = PETSc.Sys() # type: ignore
# For factorisation prefer MUMPS, then superlu_dist, then default.
if sys.hasExternalPackage("mumps"):
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
elif sys.hasExternalPackage("superlu_dist"):
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "superlu_dist"
ksp.setFromOptions()
# -
# The setting of `convergence_criterion` to `"incremental"` specifies
# that the Newton solver should compute a norm of the solution increment
# to check for convergence (the other possibility is to use
# `"residual"`, or to provide a user-defined check). The tolerance for
# convergence is specified by `rtol`
# To run the solver and save the output to a VTK file for later
# visualization, the solver is advanced in time from $t_{n}$ to
# $t_{n+1}$ until a terminal time $T$ is reached:
# +
# Output file
file = XDMFFile(MPI.COMM_WORLD, "demo_ch/output.xdmf", "w")
file.write_mesh(msh)
# Step in time
t = 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 = 1.01
# Get the sub-space for c and the corresponding dofs in the mixed space
# vector
V0, dofs = mpc.function_space.sub(0).collapse()
# Prepare viewer for plotting the solution during the computation
if have_pyvista:
# Create a VTK 'mesh' with 'nodes' at the function dofs
topology, cell_types, x = plot.vtk_mesh(V0)
grid = pv.UnstructuredGrid(topology, cell_types, x)
# Set output data
grid.point_data["c"] = u.x.array[dofs].real
grid.set_active_scalars("c")
p = pvqt.BackgroundPlotter(title="concentration", auto_update=True)
p.add_mesh(grid, clim=[0, 1])
p.view_xy(True)
p.add_text(f"time: {t}", font_size=12, name="timelabel")
c = u.sub(0)
u0.x.array[:] = u.x.array
while t < T:
#t += dt
r = solver.solve(u)
u_array = u.vector.array
u0.x.array[:] = u.x.array
file.write_function(c, t)
# Update the plot window
if have_pyvista:
p.add_text(f"time: {t:.2e}", font_size=12, name="timelabel")
grid.point_data["c"] = u.x.array[dofs].real
p.app.processEvents()
# Save data and screenshot at select timesteps
data_path = f""
np.save(data_path, u.x.array[dofs].real)
screenshot_path = f""
p.screenshot(screenshot_path)
print(f"Data saved at t = {t} for seed = {seed}")
t += dt
file.close()