{
 "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
}