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

Deleting simple Cahn Hilliard example. First trial Flory Huggins with many... 

Deleting simple Cahn Hilliard example. First trial Flory Huggins with many terms (simple substitution of y=laplace(phi_tot) instead of the entire chemical potential) still doesn't work properply. For now continuing with amended Cahn Hilliard which accounts for Flory Huggins free energy. The other Flory Huggins code probably has a coding/fenics issue I don't see.
parent 5d0ec8c1
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Fenics_Cahn_Hilliard.ipynb deleted 100644 → 0
+ 0
− 122
View file @ 5d0ec8c1
%% 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):
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 = 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, 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)
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*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 = 5000*dt
while (t < T):
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
file << (u.split()[0], t)
```
Fenics_FloryHugg_PhiTot.ipynb
+ 12
− 11
View file @ ec39d8f8
%% Cell type:code id: tags:
``` python
# 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):
if x[0] > 0.5:
values[0] = 0.63 + 0.02*(0.5 - random.random())
# values[0] = 0.63 + 0.02*(0.5 - random.random())
values[0] = 0.6
else:
values[0] = 0.4
values[0] = 0.35
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
dt = 5.0e-10 # time step
dt = 5.0e-11 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
kBTG0 = 1.0e-6
X = 3.5
X = 2.2
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, 1, 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)
# 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 = (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)
+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
-dt*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 = 200*dt
T = 2000*dt
while (t < T):
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
file << (u.split()[0], t)
t += dt
```
......
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