From de8e7f7b32c4999c8a4aeb23484f822865d67b69 Mon Sep 17 00:00:00 2001
From: Lars Hubatsch <hubatsch@pks.mpg.de>
Date: Wed, 27 Nov 2019 10:11:07 +0100
Subject: [PATCH] Adding working exampel code from fenics docs for Cahn
Hilliard.
---
Fenics_Cahn_Hilliard.ipynb | 121 +++++++++++++++++++++++++++++++++++++
1 file changed, 121 insertions(+)
create mode 100644 Fenics_Cahn_Hilliard.ipynb
diff --git a/Fenics_Cahn_Hilliard.ipynb b/Fenics_Cahn_Hilliard.ipynb
new file mode 100644
index 0000000..b38153e
--- /dev/null
+++ b/Fenics_Cahn_Hilliard.ipynb
@@ -0,0 +1,121 @@
+{
+ "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-06 # 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, 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",
+ "# Weak statement of the equations\n",
+ "L0 = c*q*dx - c0*q*dx + dt*dot(grad(mu_mid), grad(q))*dx\n",
+ "L1 = mu*v*dx - dfdc*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 = 50*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
+}
--
GitLab