My refering the tutorial available at “Extension for periodic conditons (dolfinx_mpc):” in github by Jorgen, I tried to implement periodic boundary conditions in cahn hilliard equation in both top_bottom and left_right boundaries. Below I am attaching my code for the reference and the screenshort of the result when I ran the simulation for 100 timesteps.
Code
comm = MPI.COMM_WORLD
msh = create_unit_square(MPI.COMM_WORLD, 96,96 , CellType.triangle)
ndim = msh.geometry.dim
lmbda = 1.0e-02 # surface parameter
dt = 5.0e-06 # time step
theta = 1 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicholson
P1 = element("Lagrange", msh.basix_cell(), 1, dtype=default_real_type)
ME = functionspace(msh, mixed_element([P1, P1]))
q, v = ufl.TestFunctions(ME)
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
# Split mixed functions
c, mu = ufl.split(u)
c0, mu0 = ufl.split(u0)
# Initialize boundary conditions list
bcs = []
# Define the tolerance for boundary detection
tol = 1e-8 # Adjust this as needed based on mesh size and precision
# Define periodic boundary functions for x and y directions
def left_right_boundary(x):
"""Identify points on the left (x=0) and right (x=1) boundaries for periodicity."""
return np.isclose(x[0], 1, atol=tol)
def top_bottom_boundary(x):
"""Identify points on the bottom (y=0) and top (y=1) boundaries for periodicity."""
return np.isclose(x[1], 1, atol=tol)
# Define the relation for periodic boundaries for both x and y directions
def periodic_relation_x(x):
"""Map points from right boundary (x=1) to the left boundary (x=0) and vice versa."""
out_x = np.copy(x)
out_x[0] = 1 - x[0] # Reverse x-coordinate to apply left-right periodicity
out_x[1] = x[1] # y-coordinate remains the same
return out_x
def periodic_relation_y(x):
"""Map points from top boundary (y=1) to the bottom boundary (y=0) and vice versa."""
out_x = np.copy(x)
out_x[0] = x[0] # x-coordinate remains the same
out_x[1] = 1 - x[1] # Reverse y-coordinate to apply top-bottom periodicity
return out_x
# Apply MultiPointConstraint to each subspace individually
with Timer("~PERIODIC: Initialize MPC for Cahn-Hilliard"):
# Initialize the MPC constraint object for each subspace
mpc_c = MultiPointConstraint(ME.sub(0)) # Subspace for concentration
mpc_mu = MultiPointConstraint(ME.sub(1)) # Subspace for chemical potential
# Apply periodic constraints on the concentration subspace
mpc_c.create_periodic_constraint_geometrical(ME.sub(0), left_right_boundary, periodic_relation_x, bcs)
mpc_c.create_periodic_constraint_geometrical(ME.sub(0), top_bottom_boundary, periodic_relation_y, bcs)
mpc_c.finalize()
# Apply periodic constraints on the chemical potential subspace
mpc_mu.create_periodic_constraint_geometrical(ME.sub(1), left_right_boundary, periodic_relation_x, bcs)
mpc_mu.create_periodic_constraint_geometrical(ME.sub(1), top_bottom_boundary, periodic_relation_y, bcs)
mpc_mu.finalize()
# Zero u
u.x.array[:] = 0.0
# Interpolate initial condition
rng = np.random.default_rng(42)
u.sub(0).interpolate(lambda x: 0.63 + 0.02 * (0.5 - rng.random(x.shape[1])))
u.x.scatter_forward()
c = ufl.variable(c)
f = 100 * c**2 * (1 - c) ** 2
dfdc = ufl.diff(f, c)
# mu_(n+theta)
mu_mid = (1.0 - theta) * mu0 + theta * mu
V0, dofs_c = ME.sub(0).collapse()
V1, dofs_mu = ME.sub(1).collapse()
# 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
# 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-5
solver.max_it = 200 # Maximum number of iterations
# We can customize the linear solver used inside the NewtonSolver
# 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 superlu_dist, then MUMPS, then default
if sys.hasExternalPackage("superlu_dist"):
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "superlu_dist"
elif sys.hasExternalPackage("mumps"):
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
t = 0.0
if "CI" in os.environ.keys() or "GITHUB_ACTIONS" in os.environ.keys():
T = 30 * dt
else:
T = 300* dt
xdmf_file = XDMFFile(MPI.COMM_WORLD, "cahn_hilliard_results_test_server_pbc_0.51.xdmf", "w")
xdmf_file.write_mesh(msh)
# Simulation loop
c = u.sub(0)
mu = u.sub(1)
u0.x.array[:] = u.x.array
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
xdmf_file.write_function(c, t)
xdmf_file.close()
I DON’T THING THAT I HAVE APPLIED THE PERIODIC BOUNDARY CONDITIONS CORRECTLY, AS THE CONCENTRATION IS NOT same AT THE LEFT AND RIGHT EDGE.
I guess, I have to somehow incorporate mpc_c and mpc_mu in the newton solver
Any help is appreciated.