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

Temporarily adding file for project with Stefano/Frank/Christoph.

parent 801c65bb
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Fenics_FloryHugg_Binary_Single_Mol.ipynb 0 → 100644
+ 291
− 0
View file @ 82342de8
%% Cell type:markdown id: tags:
# 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):
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]$
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}$
$\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:
$\nabla\mu_\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:
$\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$
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
$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:
``` python
## Extended Standard Cahn-Hilliard Example to Binary Flory Huggins
# as discussed on 28/11/2019
# Example 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
# Runs with fenics 2019.01
# The resulting .pvd file can be opened using default settings in ParaView
import matplotlib.pyplot as plt
import mshr as ms
import numpy as np
from scipy.optimize import curve_fit
import time
import random
from dolfin import *
def create_expr(p_list):
string = ''
for p in p_list:
string = (string+'(x[0]-'+str(p[0])+')*(x[0]-'+str(p[0])+')+(x[1]-'+str(p[1])+
')*(x[1]-'+str(p[1])+')+(x[2]-'+str(p[2])+')*(x[2]-'+str(p[2])+
')<=0.15*0.15 ? 0.8102 :' )
string = string + '.196'
return string
f = Expression((create_expr([[0.5, 0.5, 0.5], [0.2, 0.2, 0.2], [0.8, 0.8, 0.8],
[0.2, 0.8, 0.8], [0.2, 0.8, 0.2], [0.2, 0.2, 0.8],
[0.8, 0.2, 0.2], [0.8, 0.8, 0.2], [0.8, 0.2, 0.8],
[0.65, 0.65, 0.65], [0.35, 0.35, 0.35], [0.65, 0.35, 0.65],
[0.65, 0.65, 0.35], [0.35, 0.65, 0.65], [0.65, 0.35, 0.35],
[0.35, 0.35, 0.65], [0.35, 0.65, 0.35]]), 1), degree=1)
f = Expression((create_expr([[0., 0., 0.]]), 1), degree=1)
# 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
kappa = 2.92*10e-06 # surface parameter
dt = 0.5e-01 # time step
X = 7/3 # Chi Flory Huggins parameter
# time stepping family, e.g.:
# theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
theta = 0.5
# Form compiler optionsmai
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
mesh = Mesh('Meshes/balls.xml')
# domain = ms.Sphere(Point(0, 0, 0), 1.0)
# mesh = ms.generate_mesh(domain, 100)
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
tis = time.time()
# Split mixed functions
dc, dmu = split(du)
c, mu = split(u)
c0, mu0 = split(u0)
# Create intial conditions and interpolate
u.interpolate(f)
u0.interpolate(f)
print(time.time()-tis)
# Compute the chemical potential df/dc
c = variable(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 - (ln(c/(1-c))+X*(1-2*c))*v*dx - kappa*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["linear_solver"] = 'gmres'
solver.parameters["preconditioner"] = 'ilu'
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
```
%% Cell type:code id: tags:
``` python
# Output file
file_c = XDMFFile('c_long.xdmf')
# Step in time
t = 0.0
T = 1000*dt
ti = time.time()
while (t < T):
file_c.write(u.split()[0], t)
t += 200*dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
# xdata = np.linspace(0.1, 1, 1000)
# ydata = np.array([u0([x, 1, 1])[0] for x in xdata])
# np.savetxt('xdata'+str(t)+'.txt', xdata, delimiter=',', fmt='%.4f')
# np.savetxt('ydata'+str(t)+'.txt', ydata, delimiter=',', fmt='%.4f')
print(time.time() - ti)
file_c.close()
```
%% Cell type:code id: tags:
``` python
mesh.num_cells()
```
%% Cell type:code id: tags:
``` python
def tanh_fit(x, a, b, c, d):
return np.tanh((x-a)/b)*c+d
# xdata = np.linspace(0, x_mesh, x_mesh+1)
# ydata = u.compute_vertex_values()[0:x_mesh+1]
xdata = np.linspace(0.6, 0, 1000)
ydata = np.array([u0([0, x, 0])[0] for x in xdata])
# xdata = np.linspace(0., 0.49, 1000)
popt, pcov = curve_fit(tanh_fit, xdata, ydata, [0.2, 0.01, 1, 1])
fig, (ax1, ax2) = plt.subplots(1, 2)
fig.set_figwidth(16)
fig.set_figheight(5)
a, b, c, d = popt
from matplotlib import rc
rc('font',**{'family':'Helvetica','sans-serif':['Helvetica'],'size':16})
rc('text', usetex=True)
ax1.plot(xdata[:], ydata[:], lw=2, label='FEM model')
ax1.plot(xdata[:], tanh_fit(xdata[:], a, b, c, d), lw=3, ls='--'
, label='tanh fit')
ax1.legend()
ax1.set(xlabel=r'x', ylabel=r'$\phi_\mathrm{tot}$')
ax2.plot(xdata[:], (ydata[:]-tanh_fit(xdata[:], a, b, c, d))/(2*c))
ax2.set(xlabel=r'x', ylabel=r'relative error')
plt.show()
# fig.savefig('test.pdf')
```
%% Cell type:code id: tags:
``` python
print(b)
print(c+d)
print(1-c-d)
```
%% Cell type:code id: tags:
``` python
def tanh_fit(x, a, b, c, d):
return np.tanh((x-a)/b)*c+d
# xdata = np.linspace(0, x_mesh, x_mesh+1)
# ydata = u.compute_vertex_values()[0:x_mesh+1]
xdata = np.linspace(0.49, 0, 1000)
ydata = np.array([u0([0, x, 0])[0] for x in xdata])
xdata = np.linspace(0., 0.49, 1000)
popt, pcov = curve_fit(tanh_fit, xdata, ydata)
fig, (ax1, ax2) = plt.subplots(1, 2)
fig.set_figwidth(16)
fig.set_figheight(5)
a, b, c, d = popt
from matplotlib import rc
rc('font',**{'family':'Helvetica','sans-serif':['Helvetica'],'size':16})
rc('text', usetex=True)
ax1.plot(xdata[:], ydata[:], lw=2, label='FEM model')
ax1.plot(xdata[:], tanh_fit(xdata[:], a, b, c, d), lw=3, ls='--'
, label='tanh fit')# np.tanh((xdata-52)/12.4)*0.36+0.5
ax1.legend()
ax1.set(xlabel=r'x', ylabel=r'$\phi_\mathrm{tot}$')
ax2.plot(xdata[:], (ydata[:]-tanh_fit(xdata[:], a, b, c, d))/(2*c))
ax2.set(xlabel=r'x', ylabel=r'relative error')
plt.show()
fig.savefig('test.pdf')
```
%% Cell type:code id: tags:
``` python
mesh.num_cells()
```
%% Cell type:markdown id: tags:
## Thoughts on stability of algorithm
- relatively robust to grid size changes
- initial conditions seem to have strong influence
- time step also plays a role for Newton solver
- best results with small time step and initial conditions close to final concentrations
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