Newer
Older
Lars Hubatsch
committed
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For $\\phi_{tot}$ we obtain \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",
"We use 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)) and replace $\\mu$:\n",
"\n",
"\n",
"$\\mu = \\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}$\n",
"\n",
"\n",
"Cast into the weak form we obtain:"
]
},
Lars Hubatsch
committed
{
"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",
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 numpy as np\n",
"from scipy.optimize import curve_fit\n",
"\n",
"import random\n",
"from dolfin import *\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.45) & (x[0] < 0.55):# & (x[1] > 0.48) & (x[1] < 0.52):\n",
Lars Hubatsch
committed
"# values[0] = 0.63 + 0.02*(0.5 - random.random())\n",
" values[0] = 0.6\n",
" else:\n",
" values[0] = 0.35 \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",
"# 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",
"# Model parameters\n",
"lmbda = 1e-04 # surface parameter\n",
"dt = 5.0e-03 # 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",
Lars Hubatsch
committed
"# Form compiler options\n",
"parameters[\"form_compiler\"][\"optimize\"] = True\n",
"parameters[\"form_compiler\"][\"cpp_optimize\"] = True\n",
"# Create mesh and build function space\n",
"x_mesh = 96\n",
"y_mesh = 3\n",
"mesh = UnitSquareMesh.create(x_mesh, y_mesh, CellType.Type.quadrilateral)\n",
Lars Hubatsch
committed
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
"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",
"\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",
"# Compute the chemical potential df/dc\n",
"c = variable(c)\n",
"dfdc = ln(c/(1-c))+(1-2*X*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))+(1-2*X*c))*v*dx - lmbda*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[\"convergence_criterion\"] = \"incremental\"\n",
"solver.parameters[\"relative_tolerance\"] = 1e-6\n",
"# Output file\n",
"file = File(\"output.pvd\", \"compressed\")\n",
"\n",
"# Step in time\n",
"t = 0.0\n",
"T = 1000*dt\n",
Lars Hubatsch
committed
"while (t < T):\n",
" file << (u.split()[0], t)\n",
Lars Hubatsch
committed
" t += dt\n",
" u0.vector()[:] = u.vector()\n",
" solver.solve(problem, u.vector())"
Lars Hubatsch
committed
]
},
{
"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",
"popt, pcov = curve_fit(tanh_fit, xdata, u.compute_vertex_values()[0:x_mesh+1])\n",
"plt.plot(xdata, u.compute_vertex_values()[0:x_mesh+1])\n",
Lars Hubatsch
committed
"a, b, c, d = popt\n",
"plt.plot(xdata, tanh_fit(xdata, a, b, c, d))# np.tanh((xdata-52)/12.4)*0.36+0.5"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}