Skip to content
Snippets Groups Projects
Commit f3ba88ea authored by Lars Hubatsch's avatar Lars Hubatsch
Browse files

Amended Cahn Hilliard code in a minimal way to account for Flory Huggins... 

Amended Cahn Hilliard code in a minimal way to account for Flory Huggins binary mixture with the corresponding mobility scaling. Seems to work, but not sure about the volume fraction to concentration conversion.
parent 3751a38d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Fenics_Cahn_Hilliard.ipynb
+ 16
− 15
View file @ f3ba88ea
%% Cell type:code id: tags:
``` python
# Standard Cahn-Hilliard Example that runs with fenics 2019.1
# Examples can be found at https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/cahn-hilliard/demo_cahn-hilliard.py.rst#rst-header-id1
# The resulting .pvd file can be opened using default settings in ParaView
import random
from dolfin import *
# Class representing the intial conditions
class InitialConditions(UserExpression):
def __init__(self, **kwargs):
random.seed(2 + MPI.rank(MPI.comm_world))
super().__init__(**kwargs)
def eval(self, values, x):
values[0] = 0.63 + 0.02*(0.5 - random.random())
if x[0] > 0.5:
# values[0] = 0.63 + 0.02*(0.5 - random.random())
values[0] = 0.6
else:
values[0] = 0.35
values[1] = 0.0
# values[0] = 0.63 + 0.02*(0.5 - random.random())
values[1] = 0.0
def value_shape(self):
return (2,)
# Class for interfacing with the Newton solver
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)
# Model parameters
lmbda = 1.0e-02 # surface parameter
dt = 5.0e-06 # time step
lmbda = 2.0e-02 # surface parameter
dt = 5.0e-03 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
# Form compiler options
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
# Create mesh and build function space
mesh = UnitSquareMesh.create(96, 96, CellType.Type.quadrilateral)
mesh = UnitSquareMesh.create(96, 3, CellType.Type.quadrilateral)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)
# Define trial and test functions
du = TrialFunction(ME)
q, v = TestFunctions(ME)
# Define functions
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
# Split mixed functions
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)
# Compute the chemical potential df/dc
c = variable(c)
f = 100*c**2*(1-c)**2
dfdc = diff(f, c)
# f = 100*c**2*(1-c)**2
# dfdc = diff(f, c)
X = 2.2
dfdc = ln(c/(1-c))+(1-2*X*c)
# mu_(n+theta)
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
L0 = c*q*dx - c0*q*dx + dt*c*(1-c)*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 = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
# Output file
file = File("output.pvd", "compressed")
# Step in time
t = 0.0
T = 50*dt
T = 5000*dt
while (t < T):
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
file << (u.split()[0], t)
```
%% Cell type:code id: tags:
``` python
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment