{
"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",
"\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",
"\n",
"\n",
"$\\partial_t \\phi_{\\mathrm{tot}} = \\nabla \\cdot \\Gamma_{\\mathrm{tot}}\\nabla\\mu_\\mathrm{tot}$\n",
"\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",
"\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} }$"
]
},
{
"cell_type": "code",
"execution_count": null,
"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",
"# 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",
"import numpy as np\n",
"from scipy.optimize import curve_fit\n",
"import time\n",
"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",
"# Class representing the intial conditions\n",
"# 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",
"# 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",
"# Model parameters\n",
"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",
"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",
"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",
"# 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",
"# 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",
"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",
"solver.parameters[\"convergence_criterion\"] = \"incremental\"\n",
"solver.parameters[\"relative_tolerance\"] = 1e-6"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Output file\n",
"# file = File(\"output.pvd\", \"compressed\")\n",
"file_c = XDMFFile('c_long.xdmf')\n",
"\n",
"# Step in time\n",
"t = 0.0\n",
"T = 3*dt\n",
"ti = time.time()\n",
"while (t < T):\n",
"# file << (u.split()[0], t)\n",
" file_c.write(u.split()[0], t)\n",
" t += dt\n",
" u0.vector()[:] = u.vector()\n",
" solver.solve(problem, u.vector())\n",
" print(time.time() - ti)\n",
"file_c.close()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"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",
"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",
"plt.plot(xdata[65:], ydata[65:])\n",
"a, b, c, d = popt\n",
"plt.plot(xdata[65:], tanh_fit(xdata[65:], a, b, c, d))# np.tanh((xdata-52)/12.4)*0.36+0.5"
]
},
{
"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"
]
}
],
"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"
}
},
"nbformat": 4,
"nbformat_minor": 4
}