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

Updating Fenics_FloryHugg_Binary examples for Kafa.

parent e1a73c16
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Fenics_FloryHugg_Binary.ipynb
+ 15
− 13
View file @ 89a9ec55
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Binary mixture based on Flory Huggins free energy # Binary mixture based on Flory Huggins free energy
For $\phi_{tot}$ in a binary mixture, assuming a Flory Huggins free energy density, we obtain (Overleaf): For $\phi_{tot}$ in a binary mixture, assuming a Flory Huggins free energy density, we obtain (Overleaf):
Reference for FEM: https://www.comsol.com/multiphysics/finite-element-method Reference for FEM: https://www.comsol.com/multiphysics/finite-element-method
$\partial_\mathrm{t}\phi_\mathrm{tot}=k_\mathrm{B}T\Gamma_\mathrm{0} \,\nabla\cdot\left[\left[1-2\chi\phi_\mathrm{tot}(1-\phi_\mathrm{tot})\right]\nabla \phi_\mathrm{tot}-\phi_\mathrm{tot}(1-\phi_\mathrm{tot})\kappa\nabla\nabla^2\phi_\mathrm{tot}\right]$ $\partial_\mathrm{t}\phi_\mathrm{tot}=k_\mathrm{B}T\Gamma_\mathrm{0} \,\nabla\cdot\left[\left[1-2\chi\phi_\mathrm{tot}(1-\phi_\mathrm{tot})\right]\nabla \phi_\mathrm{tot}-\phi_\mathrm{tot}(1-\phi_\mathrm{tot})\kappa\nabla\nabla^2\phi_\mathrm{tot}\right]$
To solve this equation we follow the [Fenics Cahn-Hilliard example](https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/cahn-hilliard/demo_cahn-hilliard.py.rst#rst-header-id1) (the maths is explained [here](https://fenicsproject.org/docs/dolfinx/dev/python/demos/cahn-hilliard/demo_cahn-hilliard.py.html), for our free energy see overleaf). To solve this equation we follow the [Fenics Cahn-Hilliard example](https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/cahn-hilliard/demo_cahn-hilliard.py.rst#rst-header-id1) (the maths is explained [here](https://fenicsproject.org/docs/dolfinx/dev/python/demos/cahn-hilliard/demo_cahn-hilliard.py.html), for our free energy see overleaf).
$\partial_t \phi_{\mathrm{tot}} = \nabla \cdot \Gamma_{\mathrm{tot}}\nabla\mu_\mathrm{tot}$ $\partial_t \phi_{\mathrm{tot}} = \nabla \cdot \Gamma_{\mathrm{tot}}\nabla\mu_\mathrm{tot}$
$\mu_\mathrm{tot}=\nu\frac{\partial f}{\partial \phi_\mathrm{tot}}=\,k_\mathrm{B}T\left[\ln\phi_\mathrm{tot} - \ln\left(1-\phi_\mathrm{tot}\right)+\chi(1-2\phi_\mathrm{tot})-\kappa\nabla^2 \phi_\mathrm{tot} +1-n \right]$ $\mu_\mathrm{tot}=\nu\frac{\partial f}{\partial \phi_\mathrm{tot}}=\,k_\mathrm{B}T\left[\ln\phi_\mathrm{tot} - \ln\left(1-\phi_\mathrm{tot}\right)+\chi(1-2\phi_\mathrm{tot})-\kappa\nabla^2 \phi_\mathrm{tot} +1-n \right]$
The **boundary conditions** for the flux of the chemical potential and concentration read: The **boundary conditions** for the flux of the chemical potential and concentration read:
$\nabla\mu_\mathrm{tot}\cdot n=0 $ on $\partial \Omega\;,$ $\nabla\mu_\mathrm{tot}\cdot n=0 $ on $\partial \Omega\;,$
$\nabla\phi_\mathrm{tot}\cdot n = 0$ on $\partial \Omega$ $\nabla\phi_\mathrm{tot}\cdot n = 0$ on $\partial \Omega$
Cast into the **weak form** and integrating by parts we obtain: Cast into the **weak form** and integrating by parts we obtain:
$\int_{\Omega} \partial_t c \;q \; dx = \int_{\partial \Omega} \Gamma_\mathrm{tot}\nabla\mu_\mathrm{tot}q \; ds - \int_\Omega \Gamma_\mathrm{tot} \nabla \mu_\mathrm{tot} \nabla q\;dx = - \int_\Omega \Gamma_\mathrm{tot} \nabla \mu_\mathrm{tot} \nabla q\;dx\;, $ $\int_{\Omega} \partial_t c \;q \; dx = \int_{\partial \Omega} \Gamma_\mathrm{tot}\nabla\mu_\mathrm{tot}q \; ds - \int_\Omega \Gamma_\mathrm{tot} \nabla \mu_\mathrm{tot} \nabla q\;dx = - \int_\Omega \Gamma_\mathrm{tot} \nabla \mu_\mathrm{tot} \nabla q\;dx\;, $
$\int_\Omega \mu_\mathrm{tot} q \; dx = k_\mathrm{B} T \int_\Omega\left(\ln\left(\frac{\phi_\mathrm{tot}}{1-\phi_\mathrm{tot}}\right)+\chi(1-2\phi_\mathrm{tot})\right)q + \kappa\nabla\phi_\mathrm{tot}\nabla q \; dx$ $\int_\Omega \mu_\mathrm{tot} q \; dx = k_\mathrm{B} T \int_\Omega\left(\ln\left(\frac{\phi_\mathrm{tot}}{1-\phi_\mathrm{tot}}\right)+\chi(1-2\phi_\mathrm{tot})\right)q + \kappa\nabla\phi_\mathrm{tot}\nabla q \; dx$
where we dropped the boundary terms (not explicitly written in second equation) due to the boundary conditions/weak form. We also disregard 1-n, because constants don't change the flux in the absence of chemical reactions (an offset doesn't change the solution in this case, try out by adding $-6.0*v*dx$ to L1, for this check out fenics tutorial Poisson equation). where we dropped the boundary terms (not explicitly written in second equation) due to the boundary conditions/weak form. We also disregard 1-n, because constants don't change the flux in the absence of chemical reactions (an offset doesn't change the solution in this case, try out by adding $-6.0*v*dx$ to L1, for this check out fenics tutorial Poisson equation).
According to the review (Weber et al. 2019, eqn. 2.29) the boundary length scale is According to the review (Weber et al. 2019, eqn. 2.29) the boundary length scale is
$w = \sqrt{2\kappa/b} = \left( \frac{\kappa}{k_\mathrm{B}T\chi}\right) ^{3/5} \sqrt{ \frac{\chi}{\chi-2} }$ $w = \sqrt{2\kappa/b} = \left( \frac{\kappa}{k_\mathrm{B}T\chi}\right) ^{3/5} \sqrt{ \frac{\chi}{\chi-2} }$
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Extended Standard Cahn-Hilliard Example to Binary Flory Huggins # Extended Standard Cahn-Hilliard Example to Binary Flory Huggins
# as discussed on 28/11/2019 # as discussed on 28/11/2019
# Example can be found at # Example can be found at
# https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/ # https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/
# documented/cahn-hilliard/demo_cahn-hilliard.py.rst#rst-header-id1 # documented/cahn-hilliard/demo_cahn-hilliard.py.rst#rst-header-id1
# Runs with fenics 2019.01 # Runs with fenics 2019.01
# The resulting .pvd file can be opened using default settings in ParaView # The resulting .pvd file can be opened using default settings in ParaView
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import mshr as ms import mshr as ms
import numpy as np import numpy as np
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
import time import time
import random import random
from dolfin import * from dolfin import *
fol = '/Users/hubatsch/ownCloud/Dropbox_Lars_Christoph/frap_theory/' fol = '/Users/hubatsch/Desktop/frap_theory/'
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def create_expr(p_list): def create_expr(p_list):
string = '' string = ''
for p in p_list: for p in p_list:
string = (string+'(x[0]-'+str(p[0])+')*(x[0]-'+str(p[0])+ string = (string+'(x[0]-'+str(p[0])+')*(x[0]-'+str(p[0])+
')+(x[1]-'+str(p[1])+ ')+(x[1]-'+str(p[1])+
')*(x[1]-'+str(p[1])+')+(x[2]-'+str(p[2])+')*(x[2]-'+ ')*(x[1]-'+str(p[1])+')+(x[2]-'+str(p[2])+')*(x[2]-'+
str(p[2])+')<=.05*.05 ? .76 :' ) str(p[2])+')<=.05*.05 ? .76 :' )
string = string + '.268' string = string + '.268'
return string return string
# f = Expression((create_expr([[0.5, 0.5, 0.5]]), 1), degree=1) # f = Expression((create_expr([[0.5, 0.5, 0.5]]), 1), degree=1)
file = fol+'Meshes/points_clean.txt' file = fol+'Meshes/points_clean.txt'
points = np.loadtxt(file, delimiter=',') points = np.loadtxt(file, delimiter=',')
f = Expression((create_expr(points), 1), degree=1) f = Expression((create_expr(points), 1), degree=1)
# Class for interfacing with the Newton solver # Class for interfacing with the Newton solver
class CahnHilliardEquation(NonlinearProblem): class CahnHilliardEquation(NonlinearProblem):
def __init__(self, a, L): def __init__(self, a, L):
NonlinearProblem.__init__(self) NonlinearProblem.__init__(self)
self.L = L self.L = L
self.a = a self.a = a
def F(self, b, x): def F(self, b, x):
assemble(self.L, tensor=b) assemble(self.L, tensor=b)
def J(self, A, x): def J(self, A, x):
assemble(self.a, tensor=A) assemble(self.a, tensor=A)
# Model parameters # Model parameters
kappa = 0.25*10e-05 # surface parameter kappa = 0.25*10e-05 # surface parameter
dt = 0.5e-02 # time step dt = 0.5e-02 # time step
X = 2.2 # Chi Flory Huggins parameter X = 2.2 # Chi Flory Huggins parameter
# time stepping family, e.g.: # time stepping family, e.g.:
# theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson # theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
theta = 0.5 theta = 0.5
# Form compiler optionsmai # Form compiler optionsmai
parameters["form_compiler"]["optimize"] = True parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True parameters["form_compiler"]["cpp_optimize"] = True
# mesh = Mesh('/Users/hubatsch/Desktop/test_par/Meshes/spheres.xml') # mesh = Mesh('/Users/hubatsch/Desktop/test_par/Meshes/spheres.xml')
# domain = ms.Sphere(Point(0, 0, 0), 1.0) # domain = ms.Sphere(Point(0, 0, 0), 1.0)
# mesh = ms.generate_mesh(domain, 100) # mesh = ms.generate_mesh(domain, 100)
mesh = Mesh() mesh = Mesh()
with XDMFFile("Meshes/mesh1.xdmf") as infile: with XDMFFile("Meshes/mesh1.xdmf") as infile:
infile.read(mesh) infile.read(mesh)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1) P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1) ME = FunctionSpace(mesh, P1*P1)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1) P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1) ME = FunctionSpace(mesh, P1*P1)
# Define trial and test functions # Define trial and test functions
du = TrialFunction(ME) du = TrialFunction(ME)
q, v = TestFunctions(ME) q, v = TestFunctions(ME)
# Define functions # Define functions
u = Function(ME) # current solution u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step u0 = Function(ME) # solution from previous converged step
tis = time.time() tis = time.time()
# Split mixed functions # Split mixed functions
dc, dmu = split(du) dc, dmu = split(du)
c, mu = split(u) c, mu = split(u)
c0, mu0 = split(u0) c0, mu0 = split(u0)
# Create intial conditions and interpolate # Create intial conditions and interpolate
# u_init = InitialConditions(degree=1) # u_init = InitialConditions(degree=1)
# u.interpolate(u_init) # u.interpolate(u_init)
# u0.interpolate(u_init) # u0.interpolate(u_init)
u.interpolate(f) u.interpolate(f)
u0.interpolate(f) u0.interpolate(f)
print(time.time()-tis) print(time.time()-tis)
# Compute the chemical potential df/dc # Compute the chemical potential df/dc
c = variable(c) c = variable(c)
# mu_(n+theta) # mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu mu_mid = (1.0-theta)*mu0 + theta*mu
# Weak statement of the equations # Weak statement of the equations
L0 = c*q*dx - c0*q*dx + dt*c*(1-c)*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 - (ln(c/(1-c))+X*(1-2*c))*v*dx - kappa*dot(grad(c), grad(v))*dx L1 = mu*v*dx - (ln(c/(1-c))+X*(1-2*c))*v*dx - kappa*dot(grad(c), grad(v))*dx
L = L0 + L1 L = L0 + L1
# Compute directional derivative about u in the direction of du (Jacobian) # Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du) a = derivative(L, u, du)
# Create nonlinear problem and Newton solver # Create nonlinear problem and Newton solver
problem = CahnHilliardEquation(a, L) problem = CahnHilliardEquation(a, L)
solver = NewtonSolver() solver = NewtonSolver()
# solver.parameters["linear_solver"] = "lu" # solver.parameters["linear_solver"] = "lu"
solver.parameters["linear_solver"] = 'gmres' solver.parameters["linear_solver"] = 'gmres'
solver.parameters["preconditioner"] = 'ilu' solver.parameters["preconditioner"] = 'ilu'
# solver.parameters["convergence_criterion"] = "incremental" # solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6 solver.parameters["relative_tolerance"] = 1e-6
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
mesh.num_cells() mesh.num_cells()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Output file # Output file
# file = File("output.pvd", "compressed") # file = File("output.pvd", "compressed")
file_c = XDMFFile('c_long.xdmf') file_c = XDMFFile('c_long.xdmf')
# Step in time # Step in time
t = 0.0 t = 0.0
T = 70*dt T = 5*dt
ti = time.time() ti = time.time()
while (t < T): while (t < T):
# file << (u.split()[0], t) # file << (u.split()[0], t)
file_c.write(u.split()[0], t) file_c.write(u.split()[0], t)
t += dt t += dt
u0.vector()[:] = u.vector() u0.vector()[:] = u.vector()
solver.solve(problem, u.vector()) solver.solve(problem, u.vector())
print(time.time() - ti) print(time.time() - ti)
file_c.close() file_c.close()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def tanh_fit(x, a, b, c, d): def tanh_fit(x, a, b, c, d):
return np.tanh((x-a)/b)*c+d return np.tanh((x-a)/b)*c+d
# xdata = np.linspace(0, x_mesh, x_mesh+1) # xdata = np.linspace(0, x_mesh, x_mesh+1)
# ydata = u.compute_vertex_values()[0:x_mesh+1] # ydata = u.compute_vertex_values()[0:x_mesh+1]
xdata = np.linspace(0.1, 0.5, 100) xdata = np.linspace(0.1, 0.5, 100)
ydata = np.array([u0([x, 0.5, 0.5])[0] for x in xdata]) ydata = np.array([u0([x, 0.5, 0.5])[0] for x in xdata])
popt, pcov = curve_fit(tanh_fit, xdata, ydata) popt, pcov = curve_fit(tanh_fit, xdata, ydata)
plt.plot(xdata[65:], ydata[65:], lw=1) plt.plot(xdata[65:], ydata[65:], lw=1)
a, b, c, d = popt a, b, c, d = popt
plt.plot(xdata[65:], tanh_fit(xdata[65:], a, b, c, d), lw=1) plt.plot(xdata[65:], tanh_fit(xdata[65:], a, b, c, d), lw=1)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Thoughts on stability of algorithm ### Thoughts on stability of algorithm
- relatively robust to grid size changes - relatively robust to grid size changes
- initial conditions seem to have strong influence - initial conditions seem to have strong influence
- time step also plays a role for Newton solver - time step also plays a role for Newton solver
- best results with small time step and initial conditions close to final concentrations - best results with small time step and initial conditions close to final concentrations
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
**Example for how to use Spheres class to create n positions of non-overlapping spheres inside a box of dimensions ls.** ### Example for use of 'Spheres' class
...to create n positions of non-overlapping spheres inside a box of dimensions ls.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from fem_utils import Spheres from fem_utils import Spheres
a = Spheres(n=17, ls = [2, 2, 2]) a = Spheres(n=17, ls = [2, 2, 2])
a.create_list() a.create_list()
a.point_list a.point_list
a.write_to_txt() a.write_to_txt()
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
**Example on how to use convert_msh_to_xdmf to convert a mesh created in gmsh via 'gmsh -3 file -format msh2' to xdmf, which can be parallelised, in contrast to .xml meshes.** The resulting points in points.txt can now be used in a .geo file, e.g. box_with_spheres. Simply copy them in, then create the mesh (see below).
%% Cell type:code id: tags: %% Cell type:markdown id: tags:
``` python ### Example on how to use convert_msh_to_xdmf
from fem_utils import convert_msh_to_xdmf ...to convert a mesh created in gmsh via 'gmsh -3 file -format msh2' to xdmf, which can be parallelised, in contrast to .xml meshes.
path_to_msh = "/Users/hubatsch/ownCloud/Dropbox_Lars_Christoph/frap_theory/Meshes/box_with_spheres.msh"
path_to_xdmf = "mesh1.xdmf"
convert_msh_to_xdmf(path_to_msh, path_to_xdmf)
```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from fem_utils import convert_msh_to_xdmf
path_to_msh = "/Users/hubatsch/Desktop/frap_theory/Meshes/box_with_spheres.msh"
path_to_xdmf = "Meshes/mesh1.xdmf"
convert_msh_to_xdmf(path_to_msh, path_to_xdmf)
``` ```
......
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