Adding new file with Cahn Hilliard extended to Flory Huggins binary mixture....

Adding new file with Cahn Hilliard extended to Flory Huggins binary mixture. This was used in email to Chris.

Adding new file with Cahn Hilliard extended to Flory Huggins binary mixture....
367f8eb5 · Lars Hubatsch · f3ba88ea · 367f8eb5
Commit 367f8eb5 authored 5 years ago by Lars Hubatsch
--- a/Fenics_FloryHugg_Binary.ipynb
+++ b/Fenics_FloryHugg_Binary.ipynb
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Extended Standard Cahn-Hilliard Example to Binary Flory Huggins as discussed on 28/11/2019\n",
+    "# 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\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 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.5:\n",
+    "#             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  = 1.0e-02  # surface parameter\n",
+    "dt     = 5.0e-03  # time step\n",
+    "theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson\n",
+    "# Form compiler options\n",
+    "parameters[\"form_compiler\"][\"optimize\"]     = True\n",
+    "parameters[\"form_compiler\"][\"cpp_optimize\"] = True\n",
+    "# Create mesh and build function space\n",
+    "mesh = UnitSquareMesh.create(96, 1, CellType.Type.quadrilateral)\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",
+    "\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",
+    "# f    = 100*c**2*(1-c)**2\n",
+    "# dfdc = diff(f, c)\n",
+    "X = 2.5\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 = 3000*dt\n",
+    "while (t < T):\n",
+    "    t += dt\n",
+    "    u0.vector()[:] = u.vector()\n",
+    "    solver.solve(problem, u.vector())\n",
+    "    file << (u.split()[0], t)"
+   ]
+  },
+  {
+   "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, 96, 97)\n",
+    "popt, pcov = curve_fit(tanh_fit, xdata, u.compute_vertex_values()[0:97])\n",
+    "plt.plot(xdata, u.compute_vertex_values()[0:97])\n",
+    "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": []
+  }
+ ],
+ "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
+}
+%% 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 numpy as np
+from scipy.optimize import curve_fit
+
+import random
+from dolfin import *
+# Class representing the intial conditions
+class InitialConditions(UserExpression):
+    def __init__(self, **kwargs):
+        random.seed(2 + MPI.rank(MPI.comm_world))
+        super().__init__(**kwargs)
+    def eval(self, values, x):
+        if x[0] > 0.5:
+#             values[0] = 0.63 + 0.02*(0.5 - random.random())
+            values[0] = 0.6
+        else:
+            values[0] = 0.35
+        values[1] = 0.0
+#         values[0] = 0.63 + 0.02*(0.5 - random.random())
+        values[1] = 0.0
+    def value_shape(self):
+        return (2,)
+# 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
+lmbda  = 1.0e-02  # surface parameter
+dt     = 5.0e-03  # time step
+theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
+# Form compiler options
+parameters["form_compiler"]["optimize"]     = True
+parameters["form_compiler"]["cpp_optimize"] = True
+# Create mesh and build function space
+mesh = UnitSquareMesh.create(96, 1, CellType.Type.quadrilateral)
+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
+
+# Split mixed functions
+dc, dmu = split(du)
+c,  mu  = split(u)
+c0, mu0 = split(u0)
+# Create intial conditions and interpolate
+u_init = InitialConditions(degree=1)
+u.interpolate(u_init)
+u0.interpolate(u_init)
+# Compute the chemical potential df/dc
+c = variable(c)
+# f    = 100*c**2*(1-c)**2
+# dfdc = diff(f, c)
+X = 2.5
+dfdc = ln(c/(1-c))+(1-2*X*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))+(1-2*X*c))*v*dx - lmbda*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["convergence_criterion"] = "incremental"
+solver.parameters["relative_tolerance"] = 1e-6
+# Output file
+file = File("output.pvd", "compressed")
+
+# Step in time
+t = 0.0
+T = 3000*dt
+while (t < T):
+    t += dt
+    u0.vector()[:] = u.vector()
+    solver.solve(problem, u.vector())
+    file << (u.split()[0], t)
+```
+
+%% 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, 96, 97)
+popt, pcov = curve_fit(tanh_fit, xdata, u.compute_vertex_values()[0:97])
+plt.plot(xdata, u.compute_vertex_values()[0:97])
+a, b, c, d = popt
+plt.plot(xdata, tanh_fit(xdata, a, b, c, d))# np.tanh((xdata-52)/12.4)*0.36+0.5
+```
+
+%% Cell type:code id: tags:
+
+``` python
+```