Newer
Older
Lars Hubatsch
committed
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Binary mixture based on Flory Huggins free energy\n",
"For $\\phi_{tot}$ in a binary mixture, assuming a Flory Huggins free energy density, we obtain (Overleaf):\n",
"Reference for FEM: https://www.comsol.com/multiphysics/finite-element-method\n",
"\n",
"$\\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]$\n",
"\n",
"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).\n",
"$\\partial_t \\phi_{\\mathrm{tot}} = \\nabla \\cdot \\Gamma_{\\mathrm{tot}}\\nabla\\mu_\\mathrm{tot}$\n",
"$\\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]$\n",
"The **boundary conditions** for the flux of the chemical potential and concentration read:\n",
"\n",
"$\\nabla\\mu_\\mathrm{tot}\\cdot n=0 $ on $\\partial \\Omega\\;,$ \n",
"\n",
"$\\nabla\\phi_\\mathrm{tot}\\cdot n = 0$ on $\\partial \\Omega$\n",
"\n",
"Cast into the **weak form** and integrating by parts we obtain:\n",
"\n",
"$\\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\\;, $\n",
"\n",
"$\\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$\n",
"\n",
"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).\n",
"\n",
"According to the review (Weber et al. 2019, eqn. 2.29) the boundary length scale is \n",
"\n",
"$w = \\sqrt{2\\kappa/b} = \\left( \\frac{\\kappa}{k_\\mathrm{B}T\\chi}\\right) ^{3/5} \\sqrt{ \\frac{\\chi}{\\chi-2} }$"
Lars Hubatsch
committed
{
"cell_type": "code",
Lars Hubatsch
committed
"metadata": {},
"outputs": [],
"source": [
"# Extended Standard Cahn-Hilliard Example to Binary Flory Huggins \n",
"# as discussed on 28/11/2019\n",
"# Example can be found at \n",
"# https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/\n",
"# documented/cahn-hilliard/demo_cahn-hilliard.py.rst#rst-header-id1\n",
Lars Hubatsch
committed
"# Runs with fenics 2019.01\n",
"# The resulting .pvd file can be opened using default settings in ParaView\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import mshr as ms\n",
Lars Hubatsch
committed
"import numpy as np\n",
"from scipy.optimize import curve_fit\n",
"import time\n",
Lars Hubatsch
committed
"import random\n",
"from dolfin import *\n",
"def create_expr(p_list):\n",
" string = ''\n",
" for p in p_list:\n",
" string = (string+'(x[0]-'+str(p[0])+')*(x[0]-'+str(p[0])+')+(x[1]-'+str(p[1])+\n",
" ')*(x[1]-'+str(p[1])+')+(x[2]-'+str(p[2])+')*(x[2]-'+str(p[2])+\n",
" ')<=.05*.05 ? .76 :' )\n",
" string = string + '.268'\n",
" return string\n",
"\n",
"f = Expression((create_expr([[0.5, 0.5, 0.5], [0.2, 0.2, 0.2], [0.8, 0.8, 0.8],\n",
" [0.2, 0.8, 0.8], [0.2, 0.8, 0.2], [0.2, 0.2, 0.8],\n",
" [0.8, 0.2, 0.2], [0.8, 0.8, 0.2], [0.8, 0.2, 0.8],\n",
" [0.65, 0.65, 0.65], [0.35, 0.35, 0.35], [0.65, 0.35, 0.65],\n",
" [0.65, 0.65, 0.35], [0.35, 0.65, 0.65], [0.65, 0.35, 0.35],\n",
" [0.35, 0.35, 0.65], [0.35, 0.65, 0.35]]), 1), degree=1)\n",
"\n",
Lars Hubatsch
committed
"# Class representing the intial conditions\n",
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
"# class InitialConditions(UserExpression):\n",
"# def __init__(self, **kwargs):\n",
"# random.seed(2 + MPI.rank(MPI.comm_world))\n",
"# super().__init__(**kwargs)\n",
"# def eval(self, values, x):\n",
"# # if (x[0] > -0.1) & (x[0] < 0.1) & (x[1] > -0.1) & (x[1] < 0.1) & (x[2] > -0.1) & (x[2] < 0.1):\n",
"# # # values[0] = 0.63 + 0.02*(0.5 - random.random())\n",
"# # values[0] = 0.82\n",
"# # if (x[0]**2+x[1]**2+x[2]**2 < 0.25):\n",
"# # values[0] = 0.73\n",
"# # if (self.in_rad(x, [0.5, 0.5, 0.5], 0.05)):\n",
"# # values[0] = 0.75\n",
"# if (self.in_rad(x, [0.5, 0.5, 0.5], 0.05) or\n",
"# self.in_rad(x, [0.2, 0.2, 0.2], 0.05) or\n",
"# self.in_rad(x, [0.8, 0.8, 0.8], 0.05) or\n",
"# self.in_rad(x, [0.2, 0.8, 0.8], 0.05) or\n",
"# self.in_rad(x, [0.2, 0.8, 0.2], 0.05) or\n",
"# self.in_rad(x, [0.2, 0.2, 0.8], 0.05) or\n",
"# self.in_rad(x, [0.8, 0.2, 0.2], 0.05) or\n",
"# self.in_rad(x, [0.8, 0.8, 0.2], 0.05) or\n",
"# self.in_rad(x, [0.8, 0.2, 0.8], 0.05) or\n",
"# self.in_rad(x, [0.65, 0.65, 0.65], 0.05) or\n",
"# self.in_rad(x, [0.35, 0.35, 0.35], 0.05) or\n",
"# self.in_rad(x, [0.65, 0.35, 0.65], 0.05) or\n",
"# self.in_rad(x, [0.65, 0.65, 0.35], 0.05) or\n",
"# self.in_rad(x, [0.35, 0.65, 0.65], 0.05) or\n",
"# self.in_rad(x, [0.65, 0.35, 0.35], 0.05) or\n",
"# self.in_rad(x, [0.35, 0.35, 0.65], 0.05) or\n",
"# self.in_rad(x, [0.35, 0.65, 0.35], 0.05)):\n",
"# values[0] = 0.76\n",
"# else:\n",
"# values[0] = 0.268 \n",
"# values[1] = 0.0\n",
"# # values[0] = 0.63 + 0.02*(0.5 - random.random())\n",
"# values[1] = 0.0\n",
"# def value_shape(self):\n",
"# return (2,)\n",
"# def in_rad(self, x, drop, rad):\n",
"# if np.sqrt((x[0]-drop[0])**2+(x[1]-drop[1])**2+(x[2]-drop[2])**2)<rad:\n",
"# return True\n",
"# else:\n",
"# return False\n",
"\n",
Lars Hubatsch
committed
"# Class for interfacing with the Newton solver\n",
"class CahnHilliardEquation(NonlinearProblem):\n",
" def __init__(self, a, L):\n",
" NonlinearProblem.__init__(self)\n",
" self.L = L\n",
" self.a = a\n",
" def F(self, b, x):\n",
" assemble(self.L, tensor=b)\n",
" def J(self, A, x):\n",
" assemble(self.a, tensor=A)\n",
"\n",
Lars Hubatsch
committed
"# Model parameters\n",
Lars Hubatsch
committed
"kappa = 0.25*10e-05 # surface parameter\n",
"dt = 0.5e-02 # time step\n",
"X = 2.2 # Chi Flory Huggins parameter\n",
"# time stepping family, e.g.: \n",
"# theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson\n",
"theta = 0.5 \n",
"# Form compiler optionsmai\n",
Lars Hubatsch
committed
"parameters[\"form_compiler\"][\"optimize\"] = True\n",
"parameters[\"form_compiler\"][\"cpp_optimize\"] = True\n",
"mesh = Mesh('Meshes/box_spheres.xml')\n",
"# domain = ms.Sphere(Point(0, 0, 0), 1.0)\n",
"# mesh = ms.generate_mesh(domain, 100)\n",
Lars Hubatsch
committed
"P1 = FiniteElement(\"Lagrange\", mesh.ufl_cell(), 1)\n",
"ME = FunctionSpace(mesh, P1*P1)\n",
"# Define trial and test functions\n",
"du = TrialFunction(ME)\n",
"q, v = TestFunctions(ME)\n",
"# Define functions\n",
"u = Function(ME) # current solution\n",
"u0 = Function(ME) # solution from previous converged step\n",
"tis = time.time()\n",
Lars Hubatsch
committed
"# Split mixed functions\n",
"dc, dmu = split(du)\n",
"c, mu = split(u)\n",
"c0, mu0 = split(u0)\n",
"# Create intial conditions and interpolate\n",
"# u_init = InitialConditions(degree=1)\n",
"# u.interpolate(u_init)\n",
"# u0.interpolate(u_init)\n",
"u.interpolate(f)\n",
"u0.interpolate(f)\n",
"print(time.time()-tis)\n",
Lars Hubatsch
committed
"# Compute the chemical potential df/dc\n",
"c = variable(c)\n",
"# mu_(n+theta)\n",
"mu_mid = (1.0-theta)*mu0 + theta*mu\n",
"# Weak statement of the equations\n",
"L0 = c*q*dx - c0*q*dx + dt*c*(1-c)*dot(grad(mu_mid), grad(q))*dx\n",
"L1 = mu*v*dx - (ln(c/(1-c))+X*(1-2*c))*v*dx - kappa*dot(grad(c), grad(v))*dx \n",
Lars Hubatsch
committed
"L = L0 + L1\n",
"# Compute directional derivative about u in the direction of du (Jacobian)\n",
"a = derivative(L, u, du)\n",
"# Create nonlinear problem and Newton solver\n",
"problem = CahnHilliardEquation(a, L)\n",
"solver = NewtonSolver()\n",
"# solver.parameters[\"linear_solver\"] = \"lu\"\n",
"solver.parameters[\"linear_solver\"] = 'gmres'\n",
"solver.parameters[\"preconditioner\"] = 'ilu'\n",
Lars Hubatsch
committed
"solver.parameters[\"convergence_criterion\"] = \"incremental\"\n",
Lars Hubatsch
committed
"solver.parameters[\"relative_tolerance\"] = 1e-6"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
Lars Hubatsch
committed
"# Output file\n",
"# file = File(\"output.pvd\", \"compressed\")\n",
"file_c = XDMFFile('c_long.xdmf')\n",
Lars Hubatsch
committed
"\n",
"# Step in time\n",
"t = 0.0\n",
"T = 3*dt\n",
"ti = time.time()\n",
Lars Hubatsch
committed
"while (t < T):\n",
"# file << (u.split()[0], t)\n",
" file_c.write(u.split()[0], t)\n",
Lars Hubatsch
committed
" t += dt\n",
" u0.vector()[:] = u.vector()\n",
" solver.solve(problem, u.vector())\n",
" print(time.time() - ti)\n",
"file_c.close()"
Lars Hubatsch
committed
]
},
{
"cell_type": "code",
Lars Hubatsch
committed
"metadata": {},
Lars Hubatsch
committed
"source": [
"def tanh_fit(x, a, b, c, d):\n",
" return np.tanh((x-a)/b)*c+d\n",
"# xdata = np.linspace(0, x_mesh, x_mesh+1)\n",
"# ydata = u.compute_vertex_values()[0:x_mesh+1]\n",
Lars Hubatsch
committed
"xdata = np.linspace(0.1, 0.5, 100)\n",
"ydata = np.array([u0([x, 0.5, 0.5])[0] for x in xdata])\n",
"popt, pcov = curve_fit(tanh_fit, xdata, ydata)\n",
Lars Hubatsch
committed
"plt.plot(xdata[65:], ydata[65:])\n",
Lars Hubatsch
committed
"a, b, c, d = popt\n",
Lars Hubatsch
committed
"plt.plot(xdata[65:], tanh_fit(xdata[65:], a, b, c, d))# np.tanh((xdata-52)/12.4)*0.36+0.5"
Lars Hubatsch
committed
]
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mesh.num_cells()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Thoughts on stability of algorithm\n",
"- relatively robust to grid size changes\n",
"- initial conditions seem to have strong influence\n",
"- time step also plays a role for Newton solver\n",
"- best results with small time step and initial conditions close to final concentrations"
]
Lars Hubatsch
committed
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.2"
Lars Hubatsch
committed
}
},
"nbformat": 4,
"nbformat_minor": 4
}