Fenics_FloryHugg_PhiTot.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Adapted for ternary/binary mixture (FRAP project) from Standard Cahn-Hilliard Example at \n",
    "# https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/cahn-hilliard/demo_cahn-hilliard.py.rst#rst-header-id1\n",
    "# The resulting .pvd file can be opened using default settings in ParaView\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",
    "        else:\n",
    "            values[0] = 0.4                \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",
    "dt     = 5.0e-10  # time step\n",
    "theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson\n",
    "kBTG0  = 1.0e-6\n",
    "X = 3.5\n",
    "kappa = 2\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",
    "# mu_(n+theta)\n",
    "mu_mid = (1.0-theta)*mu0 + theta*mu\n",
    "c_mid =  (1.0-theta)*c0 + theta*c\n",
    "# Weak statement of the equations\n",
    "# L0 = c*q*dx - c0*q*dx + dt*dot(grad(mu_mid), grad(q))*dx\n",
    "# L0 needs to include the next time step via theta, while L1 is indendent of \n",
    "# time and therefore doesn't need discretization in time, check this is true with Chris!\n",
    "# In particular the c_mid in L0!\n",
    "L0 = (1/kBTG0*(c*q*dx - c0*q*dx) - dt*((1+2*X*c_mid)*mu_mid- 2*X*dot(grad(mu_mid), grad(mu_mid))\n",
    "                            +2*X*(2*c_mid*dot(grad(c_mid), grad(c_mid)) + c_mid**2*mu_mid) \n",
    "                            -kappa*(dot(grad(c_mid), grad(mu_mid)) - 2*c_mid*dot(grad(c_mid), grad(mu_mid))))*q*dx\n",
    "                            -kappa*(c_mid**2*dot(grad(mu_mid), grad(q))+2*c_mid*q*dot(grad(mu_mid), grad(c_mid)))*dx)\n",
    "# L1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx\n",
    "L1 = mu*v*dx + 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 = 200*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)"
   ]
  }
 ],
 "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
}