Solving Cahn-Hilliard equation with a multi-grid method.
Hi, I am a beginner at FEniCS.
I have tried solving the Chan-Hilliard equation with a multi-grid solver.
But It is not stable, I must set up a very smaller time-step than the standard Krylov-solver.
I think my implementation is wrong. But I have no idea. Please help me, thanks.
Following is my source code.
import random
from dolfin import *
from mpi4py import MPI
comm = MPI.COMM_WORLD
N = 60
mesh_1 = UnitSquareMesh(N, N)
mesh_2 = refine(mesh_1)
mesh_3 = refine(mesh_2)
meshes = [mesh_1, mesh_2, mesh_3]
P1 = FiniteElement("Lagrange", mesh_1.ufl_cell(), 1)
P2 = FiniteElement("Lagrange", mesh_2.ufl_cell(), 1)
P3 = FiniteElement("Lagrange", mesh_3.ufl_cell(), 1)
MEs = [FunctionSpace(mesh_1, P1*P1), FunctionSpace(mesh_2, P2*P2), FunctionSpace(mesh_3, P3*P3)]
print(mesh_3.num_vertices())
# Create collection of PETSc DM objects
dm_collection = PETScDMCollection(MEs)
class InitialConditions(UserExpression):
def __init__(self, **kwargs):
random.seed(2)
super().__init__(**kwargs)
def eval(self, values, x):
values[0] = 0.63 + 0.02*(0.5 - random.random())
values[1] = 0.0
def value_shape(self):
return (2,)
class CahnHilliardEquation(NonlinearProblem):
def __init__(self, a, L):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
def F(self, b, x):
assemble(self.L, tensor=b)
def J(self, A, x):
assemble(self.a, tensor=A)
class KrylovSolver(NewtonSolver):
def __init__(self):
NewtonSolver.__init__(self, comm, PETScKrylovSolver(), PETScFactory.instance())
def solver_setup(self, A, P, problem, iteration):
self.linear_solver().set_operator(A)
# Set PETSc solver type
PETScOptions.set("ksp_type", "gmres")
PETScOptions.set("pc_type", "mg")
# Set PETSc MG type and levels
PETScOptions.set("pc_mg_levels", 3)
PETScOptions.set("pc_mg_galerkin")
# Set smoother
PETScOptions.set("mg_levels_ksp_type", "chebyshev")
PETScOptions.set("mg_levels_pc_type", "jacobi")
# Set tolerance and monitor residual
PETScOptions.set("ksp_monitor_true_residual")
PETScOptions.set("ksp_atol", 1.0e-12)
PETScOptions.set("ksp_rtol", 1.0e-12)
self.linear_solver().set_from_options()
self.linear_solver().set_dm(dm_collection.get_dm(-1))
self.linear_solver().set_dm_active(False)
# Model parameters
lmbda = 1.0e-02 # surface parameter
dt = 5.0e-10 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
parameters["linear_algebra_backend"] = "PETSc"
du = TrialFunction(MEs[-1])
q, v = TestFunctions(MEs[-1])
u = Function(MEs[-1]) # current solution
u0 = Function(MEs[-1]) # solution from previous converged step
dc, dmu = split(du)
c, mu = split(u)
c0, mu0 = split(u0)
# Create intial conditions and interpolate
u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)
c = variable(c)
f = 100*c**2*(1-c)**2
dfdc = diff(f, c)
mu_mid = (1.0-theta)*mu0 + theta*mu
# Weak statement of the equations
L0 = c*q*dx - c0*q*dx + dt*dot(grad(mu_mid), grad(q))*dx
L1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx
L = L0 + L1
# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)
# Create nonlinear problem and Newton solver
problem = CahnHilliardEquation(a, L)
solver = KrylovSolver()
# Output file
file = File("output.pvd", "compressed")
# Step in time
t = 0.0
T = 50*dt
while (t < T):
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
file << (u.split()[0], t)