Parameter changes to make cool movie. Now again 2D.

5d0ec8c1 · Lars Hubatsch · 367f8eb5 · 5d0ec8c1
Commit 5d0ec8c1 authored 5 years ago by Lars Hubatsch
--- a/Fenics_FloryHugg_Binary.ipynb
+++ b/Fenics_FloryHugg_Binary.ipynb
@@ -23,7 +23,7 @@
    "        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",
+    "        if (x[0] > 0.48) & (x[0] < 0.52) & (x[1] > 0.48) & (x[1] < 0.52):\n",
    "#             values[0] = 0.63 + 0.02*(0.5 - random.random())\n",
    "            values[0] = 0.6\n",
    "        else:\n",
@@ -45,13 +45,13 @@
    "        assemble(self.a, tensor=A)\n",
    "# Model parameters\n",
    "lmbda  = 1.0e-02  # surface parameter\n",
-    "dt     = 5.0e-03  # time step\n",
+    "dt     = 5.0e-02  # 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",
+    "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",
@@ -94,12 +94,12 @@
    "\n",
    "# Step in time\n",
    "t = 0.0\n",
-    "T = 3000*dt\n",
+    "T = 100*dt\n",
    "while (t < T):\n",
+    "    file << (u.split()[0], t)\n",
    "    t += dt\n",
    "    u0.vector()[:] = u.vector()\n",
-    "    solver.solve(problem, u.vector())\n",
-    "    file << (u.split()[0], t)"
+    "    solver.solve(problem, u.vector())"
   ]
  },
  {
@@ -116,13 +116,6 @@
    "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": {

 %% Cell type:code id: tags:

 ``` python
 # Extended Standard Cahn-Hilliard Example to Binary Flory Huggins as discussed on 28/11/2019
 # 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
 # Runs with fenics 2019.01
 # The resulting .pvd file can be opened using default settings in ParaView

 import matplotlib.pyplot as plt
 import numpy as np
 from scipy.optimize import curve_fit

 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):
-        if x[0] > 0.5:
+        if (x[0] > 0.48) & (x[0] < 0.52) & (x[1] > 0.48) & (x[1] < 0.52):
 #             values[0] = 0.63 + 0.02*(0.5 - random.random())
            values[0] = 0.6
        else:
            values[0] = 0.35
        values[1] = 0.0
 #         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-03  # time step
+dt     = 5.0e-02  # time step
 theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
 # Form compiler options
 parameters["form_compiler"]["optimize"]     = True
 parameters["form_compiler"]["cpp_optimize"] = True
 # Create mesh and build function space
-mesh = UnitSquareMesh.create(96, 1, CellType.Type.quadrilateral)
+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)
 X = 2.5
 dfdc = ln(c/(1-c))+(1-2*X*c)
 # mu_(n+theta)
 mu_mid = (1.0-theta)*mu0 + theta*mu
 # Weak statement of the equations
 L0 = c*q*dx - c0*q*dx + dt*c*(1-c)*dot(grad(mu_mid), grad(q))*dx
 L1 = mu*v*dx - (ln(c/(1-c))+(1-2*X*c))*v*dx - lmbda*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 = 3000*dt
+T = 100*dt
 while (t < T):
+    file << (u.split()[0], t)
    t += dt
    u0.vector()[:] = u.vector()
    solver.solve(problem, u.vector())
-    file << (u.split()[0], t)
 ```

 %% Cell type:code id: tags:

 ``` python
 def tanh_fit(x, a, b, c, d):
     return np.tanh((x-a)/b)*c+d
 xdata = np.linspace(0, 96, 97)
 popt, pcov = curve_fit(tanh_fit, xdata, u.compute_vertex_values()[0:97])
 plt.plot(xdata, u.compute_vertex_values()[0:97])
 a, b, c, d = popt
 plt.plot(xdata, tanh_fit(xdata, a, b, c, d))# np.tanh((xdata-52)/12.4)*0.36+0.5
 ```
-
-%% Cell type:code id: tags:
-
-``` python
-```