diff --git a/Fenics_Cahn_Hilliard.ipynb b/Fenics_Cahn_Hilliard.ipynb
deleted file mode 100644
index cb709ff760d514d53a8b81202b7057c4093b47a8..0000000000000000000000000000000000000000
--- a/Fenics_Cahn_Hilliard.ipynb
+++ /dev/null
@@ -1,122 +0,0 @@
-{
- "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",
- " 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 = 2.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, 3, 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.2\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 - 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 = 5000*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
-}
diff --git a/Fenics_FloryHugg_PhiTot.ipynb b/Fenics_FloryHugg_PhiTot.ipynb
index ab2a465dcd091f7b65b3fe40550cc64e8eaf6578..9c3710e4041ed0f2bf3834314f4c30908dee00ea 100644
--- a/Fenics_FloryHugg_PhiTot.ipynb
+++ b/Fenics_FloryHugg_PhiTot.ipynb
@@ -19,9 +19,10 @@
" 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.63 + 0.02*(0.5 - random.random())\n",
+ " values[0] = 0.6\n",
" else:\n",
- " values[0] = 0.4 \n",
+ " values[0] = 0.35 \n",
" values[1] = 0.0\n",
" def value_shape(self):\n",
" return (2,)\n",
@@ -36,16 +37,16 @@
" def J(self, A, x):\n",
" assemble(self.a, tensor=A)\n",
"# Model parameters\n",
- "dt = 5.0e-10 # time step\n",
+ "dt = 5.0e-11 # 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",
+ "X = 2.2\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",
+ "mesh = UnitSquareMesh.create(96, 3, CellType.Type.quadrilateral)\n",
"P1 = FiniteElement(\"Lagrange\", mesh.ufl_cell(), 1)\n",
"ME = FunctionSpace(mesh, P1*P1)\n",
"# Define trial and test functions\n",
@@ -74,9 +75,9 @@
"# 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",
+ " +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",
+ " -dt*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",
@@ -93,12 +94,12 @@
"\n",
"# Step in time\n",
"t = 0.0\n",
- "T = 200*dt\n",
+ "T = 2000*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)"
+ " file << (u.split()[0], t)\n",
+ " t += dt"
]
}
],