Python bleaching works for sphere and slab for standard case.

29995129 · Lars Hubatsch · da695776 · 29995129
Commit 29995129 authored 4 years ago by Lars Hubatsch
--- a/FloryHugg_DiffUnbleached.ipynb
+++ b/FloryHugg_DiffUnbleached.ipynb
@@ -13,9 +13,9 @@
    "import time\n",
    "df.set_log_level(40)\n",
    "\n",
-    "domain = ms.Sphere(df.Point(0, 0, 0), 1.0)\n",
-    "# mesh = ms.generate_mesh(domain, 35)\n",
-    "mesh = df.UnitSquareMesh(50,50)\n",
+    "# domain = ms.Sphere(df.Point(0, 0, 0), 1.0)\n",
+    "# mesh = ms.generate_mesh(domain, 50)\n",
+    "mesh = df.UnitIntervalMesh(1000)\n",
    "dt = 0.01\n",
    "\n",
    "F = df.FunctionSpace(mesh, 'CG', 1)\n",
@@ -26,16 +26,21 @@
    "tc = df.TestFunction(F)\n",
    "\n",
    "# Interpolate c_tot and initial conditions\n",
-    "# c_tot.interpolate(df.Expression('-0.4*tanh((sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/0.1-0.5)) ', degree=1))\n",
-    "c_tot.interpolate(df.Expression('0.4*tanh(-35*(sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))-0.2)) +0.5', degree=1))\n",
+    "# 3D:\n",
+    "# c_tot.interpolate(df.Expression('0.4*tanh(-350*(sqrt((x[0])*(x[0])+(x[1])*(x[1])+(x[2])*(x[2]))-0.2))+0.5', degree=1))\n",
+    "# c0.interpolate(df.Expression(('(x[0]<0.5) && sqrt((x[0])*(x[0])+(x[1])*(x[1])+(x[2])*(x[2]))<0.2 ? 0 :'\n",
+    "#                               '0.4*tanh(-350*(sqrt((x[0])*(x[0])+(x[1])*(x[1])+(x[2])*(x[2]))-0.2)) + 0.5'),\n",
+    "#                              degree=1))\n",
+    "# 1D:\n",
+    "c_tot.interpolate(df.Expression('0.4*tanh(-350*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.2))+0.5', degree=1))\n",
+    "c0.interpolate(df.Expression(('sqrt((x[0]-0.5)*(x[0]-0.5))<0.2 ? 0 :'\n",
+    "                              '0.4*tanh(-350*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.2)) +0.5'),\n",
+    "                             degree=1))\n",
    "# tanh from left to right initial conditions\n",
    "# c_tot.interpolate(df.Expression('0.4*tanh(5*(x[0]-0.5)) - 1', degree=1))\n",
    "# Step function left to right initial condition\n",
    "# c0.interpolate(df.Expression('(x[0]<0.5) ? 0 : 1.0', degree=1))\n",
    "# Somewhat realistic half FRAP\n",
-    "c0.interpolate(df.Expression(('(x[0]<0.5) && sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))<0.2 ? 0 :'\n",
-    "                              '0.4*tanh(-35*(sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))-0.2)) +0.5'),\n",
-    "                             degree=1))\n",
    "\n",
    "# Weak form\n",
    "form = ((df.inner((c-c0)/dt, tc) + (1-c_tot)* df.inner(df.grad(c), df.grad(tc))) -\n",
@@ -48,8 +53,8 @@
    "\n",
    "# Solve in time\n",
    "ti = time.time()\n",
-    "for i in range(80):\n",
-    "    print(t)\n",
+    "for i in range(8):\n",
+    "    print(time.time() - ti)\n",
    "    df.solve(form == 0, c)\n",
    "    df.assign(c0, c)\n",
    "    t += dt\n",
@@ -64,10 +69,27 @@
   "metadata": {},
   "outputs": [],
   "source": [
-    "# c_tot.interpolate(df.Expression('-0.4*tanh((sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/0.1-0.1)) ', degree=1))\n",
-    "c_tot.interpolate(df.Expression('0.4*tanh(-35*(sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))-0.2)) +0.5', degree=1))\n",
-    "plt.plot(np.linspace(0, 1, 100), [c_tot([x, 0.5]) for x in np.linspace(0, 1, 100)])"
+    "# 1D:\n",
+    "plt.plot(np.linspace(0, 1, 1000), [c0([x]) for x in np.linspace(0, 1, 1000)])\n",
+    "# 3D:\n",
+    "# plt.plot(np.linspace(0, 0.5, 1000), [c0([x, 0, 0]) for x in np.linspace(0, 0.5, 1000)])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "c0([0, 0, 0])"
   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
  }
 ],
 "metadata": {

 %% Cell type:code id: tags:

 ``` python
 import dolfin as df
 import matplotlib.pyplot as plt
 import mshr as ms
 import numpy as np
 import time
 df.set_log_level(40)

-domain = ms.Sphere(df.Point(0, 0, 0), 1.0)
-# mesh = ms.generate_mesh(domain, 35)
-mesh = df.UnitSquareMesh(50,50)
+# domain = ms.Sphere(df.Point(0, 0, 0), 1.0)
+# mesh = ms.generate_mesh(domain, 50)
+mesh = df.UnitIntervalMesh(1000)
 dt = 0.01

 F = df.FunctionSpace(mesh, 'CG', 1)

 c0 = df.Function(F)
 c_tot = df.Function(F)
 c = df.Function(F)
 tc = df.TestFunction(F)

 # Interpolate c_tot and initial conditions
-# c_tot.interpolate(df.Expression('-0.4*tanh((sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/0.1-0.5)) ', degree=1))
-c_tot.interpolate(df.Expression('0.4*tanh(-35*(sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))-0.2)) +0.5', degree=1))
+# 3D:
+# c_tot.interpolate(df.Expression('0.4*tanh(-350*(sqrt((x[0])*(x[0])+(x[1])*(x[1])+(x[2])*(x[2]))-0.2))+0.5', degree=1))
+# c0.interpolate(df.Expression(('(x[0]<0.5) && sqrt((x[0])*(x[0])+(x[1])*(x[1])+(x[2])*(x[2]))<0.2 ? 0 :'
+#                               '0.4*tanh(-350*(sqrt((x[0])*(x[0])+(x[1])*(x[1])+(x[2])*(x[2]))-0.2)) + 0.5'),
+#                              degree=1))
+# 1D:
+c_tot.interpolate(df.Expression('0.4*tanh(-350*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.2))+0.5', degree=1))
+c0.interpolate(df.Expression(('sqrt((x[0]-0.5)*(x[0]-0.5))<0.2 ? 0 :'
+                              '0.4*tanh(-350*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.2)) +0.5'),
+                             degree=1))
 # tanh from left to right initial conditions
 # c_tot.interpolate(df.Expression('0.4*tanh(5*(x[0]-0.5)) - 1', degree=1))
 # Step function left to right initial condition
 # c0.interpolate(df.Expression('(x[0]<0.5) ? 0 : 1.0', degree=1))
 # Somewhat realistic half FRAP
-c0.interpolate(df.Expression(('(x[0]<0.5) && sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))<0.2 ? 0 :'
-                              '0.4*tanh(-35*(sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))-0.2)) +0.5'),
-                             degree=1))

 # Weak form
 form = ((df.inner((c-c0)/dt, tc) + (1-c_tot)* df.inner(df.grad(c), df.grad(tc))) -
        (1-c_tot)/c_tot* c*df.inner(df.grad(c_tot), df.grad(tc))) * df.dx

 # Open file for writing
 cFile = df.XDMFFile('conc.xdmf')
 t = 0
 cFile.write(c0, t)

 # Solve in time
 ti = time.time()
-for i in range(80):
-    print(t)
+for i in range(8):
+    print(time.time() - ti)
    df.solve(form == 0, c)
    df.assign(c0, c)
    t += dt
    cFile.write(c0, t)
 cFile.close()
 print(time.time() - ti)
 ```

 %% Cell type:code id: tags:

 ``` python
-# c_tot.interpolate(df.Expression('-0.4*tanh((sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/0.1-0.1)) ', degree=1))
-c_tot.interpolate(df.Expression('0.4*tanh(-35*(sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))-0.2)) +0.5', degree=1))
-plt.plot(np.linspace(0, 1, 100), [c_tot([x, 0.5]) for x in np.linspace(0, 1, 100)])
+# 1D:
+plt.plot(np.linspace(0, 1, 1000), [c0([x]) for x in np.linspace(0, 1, 1000)])
+# 3D:
+# plt.plot(np.linspace(0, 0.5, 1000), [c0([x, 0, 0]) for x in np.linspace(0, 0.5, 1000)])
+```
+
+%% Cell type:code id: tags:
+
+``` python
+c0([0, 0, 0])
+```
+
+%% Cell type:code id: tags:
+
+``` python
 ```