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

2D/1D Flory Huggins works in principle. Currently output only via paraview,... 

2D/1D Flory Huggins works in principle. Currently output only via paraview, need to figure out plotting from fenics.
parent e860bfd6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
.gitignore
+ 2
− 0
View file @ 3751a38d
.ipynb_checkpoints/
.DS_Store
*~
*.vtu
*.pvd
Fenics_FloryHugg_PhiTot.ipynb
+ 13
− 23
View file @ 3751a38d
%% 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
# Adapted for ternary/binary mixture (FRAP project) from Standard Cahn-Hilliard Example 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())
else:
values[0] = 0.4
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-14 # time step
dt = 5.0e-10 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
kBTG0 = 1
X = 2.5
kBTG0 = 1.0e-6
X = 3.5
kappa = 2
# 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, 1, 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)
# c = variable(c)
# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu
c_mid = (1.0-theta)*c0 + theta*c
# Weak statement of the equations
# L0 = c*q*dx - c0*q*dx + dt*dot(grad(mu_mid), grad(q))*dx
# L0 needs to include the next time step via theta, while L1 is indendent of
# time and therefore doesn't need discretization in time, check this is true with Chris!
# In particular the c_mid in L0!
L0 = (c*q*dx - c0*q*dx - dt*((1+2*X*c_mid)*mu_mid- 2*X*dot(grad(mu_mid), grad(mu_mid))
L0 = (1/kBTG0*(c*q*dx - c0*q*dx) - dt*((1+2*X*c_mid)*mu_mid- 2*X*dot(grad(mu_mid), grad(mu_mid))
+2*X*(2*c_mid*dot(grad(c_mid), grad(c_mid)) + c_mid**2*mu_mid)
-kappa*(dot(grad(c_mid), grad(mu_mid)) - 2*c_mid*dot(grad(c_mid), grad(mu_mid))))*q*dx
-kappa*(c_mid**2*dot(grad(mu_mid), grad(q))+2*c_mid*q*dot(grad(mu_mid), grad(c_mid)))*dx)
# dt*((1+2*X*c_mid)*mu_mid - 2*X*dot(grad(mu_mid), grad(mu_mid)) + 2*X*(2*c_mid*dot(grad(c_mid), grad(c_mid)) + c_mid**2*mu_mid)
# -kappa*(dot(grad(c_mid), grad(mu_mid)) - 2*c_mid*dot(grad(c_mid), grad(mu_mid)))*q
# -kappa*(c_mid**2*dot(grad(mu_mid), grad(q))+2*c_mid*q*dot(grad(mu_mid), grad(c_mid))))*dx)
# L1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx
L1 = mu*v*dx + 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 = 5*dt
T = 200*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