At least it does something but solution oscillates like crazy, maybe system too small?

b1d89e4d · Lars Hubatsch · de8e7f7b · b1d89e4d
Commit b1d89e4d authored 5 years ago by Lars Hubatsch
--- a/Fenics_FloryHugg_PhiTot.ipynb
+++ b/Fenics_FloryHugg_PhiTot.ipynb
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Standard Cahn-Hilliard Example that runs with fenics 2019.1\n",
+    "# Examples 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",
+    "# 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",
+    "        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-14  # time step\n",
+    "theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson\n",
+    "kBTG0  = 1\n",
+    "X = 2.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, 96, 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",
+    "# 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 = (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",
+    "#     dt*((1+2*X*c_mid)*mu_mid - 2*X*dot(grad(mu_mid), grad(mu_mid)) + 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\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 = 5*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": []
+  }
+ ],
+ "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
+# Standard Cahn-Hilliard Example that runs with fenics 2019.1
+# Examples 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
+# The resulting .pvd file can be opened using default settings in ParaView
+
+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):
+        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-14  # time step
+theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
+kBTG0  = 1
+X = 2.5
+kappa = 2
+# Form compiler options
+parameters["form_compiler"]["optimize"]     = True
+parameters["form_compiler"]["cpp_optimize"] = True
+# Create mesh and build function space
+mesh = UnitSquareMesh.create(96, 96, 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)
+# mu_(n+theta)
+mu_mid = (1.0-theta)*mu0 + theta*mu
+c_mid =  (1.0-theta)*c0 + theta*c
+# Weak statement of the equations
+# L0 = c*q*dx - c0*q*dx + dt*dot(grad(mu_mid), grad(q))*dx
+# L0 needs to include the next time step via theta, while L1 is indendent of
+# time and therefore doesn't need discretization in time, check this is true with Chris!
+# In particular the c_mid in L0!
+L0 = (c*q*dx - c0*q*dx - dt*((1+2*X*c_mid)*mu_mid- 2*X*dot(grad(mu_mid), grad(mu_mid))
+                            +2*X*(2*c_mid*dot(grad(c_mid), grad(c_mid)) + c_mid**2*mu_mid)
+                            -kappa*(dot(grad(c_mid), grad(mu_mid)) - 2*c_mid*dot(grad(c_mid), grad(mu_mid))))*q*dx
+                            -kappa*(c_mid**2*dot(grad(mu_mid), grad(q))+2*c_mid*q*dot(grad(mu_mid), grad(c_mid)))*dx)
+#     dt*((1+2*X*c_mid)*mu_mid - 2*X*dot(grad(mu_mid), grad(mu_mid)) + 2*X*(2*c_mid*dot(grad(c_mid), grad(c_mid)) + c_mid**2*mu_mid)
+#      -kappa*(dot(grad(c_mid), grad(mu_mid)) - 2*c_mid*dot(grad(c_mid), grad(mu_mid)))*q
+#      -kappa*(c_mid**2*dot(grad(mu_mid), grad(q))+2*c_mid*q*dot(grad(mu_mid), grad(c_mid))))*dx)
+# L1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx
+L1 = mu*v*dx + 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 = 5*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
+```