From 4665e64ecee55f9aece17ae6f283990f0029d993 Mon Sep 17 00:00:00 2001
From: Lars Hubatsch <hubatsch@pks.mpg.de>
Date: Thu, 17 Sep 2020 09:17:23 +0200
Subject: [PATCH] Radial Droplet FRAP qualitatively works.

But only up to a factor between different partitionings/diffusion coefficients outside.
Problem: Maybe need different concentration scaling for radial problem?
---
 FloryHugg_DiffUnbleached.ipynb | 57 ++++++++++++++++------------------
 1 file changed, 27 insertions(+), 30 deletions(-)

diff --git a/FloryHugg_DiffUnbleached.ipynb b/FloryHugg_DiffUnbleached.ipynb
index 87c67f8..74ecc4b 100644
--- a/FloryHugg_DiffUnbleached.ipynb
+++ b/FloryHugg_DiffUnbleached.ipynb
@@ -15,7 +15,7 @@
     "\n",
     "# domain = ms.Sphere(df.Point(0, 0, 0), 1.0)\n",
     "# mesh = ms.generate_mesh(domain, 50)\n",
-    "mesh = df.UnitIntervalMesh(100000)\n",
+    "mesh = df.UnitIntervalMesh(10000)\n",
     "dt = 0.000001\n",
     "\n",
     "F = df.FunctionSpace(mesh, 'CG', 1)"
@@ -30,17 +30,28 @@
     "def calc_sim(c0, c_tot, Ga0):\n",
     "    tc = df.TestFunction(F)\n",
     "    c = df.Function(F)\n",
-    "    # Weak form\n",
-    "    form = ((df.inner((c-c0)/dt, tc) +\n",
-    "             df.inner(df.grad(c), df.grad((1-c_tot)/Ga0*tc))) -\n",
-    "             df.inner(df.grad(c_tot), df.grad((1-c_tot)/Ga0/c_tot*c*tc))-\n",
-    "             tc*df.inner(df.grad(c), df.grad((1-c_tot)/Ga0))+\n",
-    "             tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c/Ga0))) * df.dx\n",
+    "    X = df.SpatialCoordinate(mesh)\n",
+    "#     X.interpolate(df.Expression('x[0]', degree=1))\n",
+    "    # Weak form 1D:\n",
+    "#     form = ((df.inner((c-c0)/dt, tc) +\n",
+    "#              df.inner(df.grad(c), df.grad((1-c_tot)/Ga0*tc))) -\n",
+    "#              df.inner(df.grad(c_tot), df.grad((1-c_tot)/Ga0/c_tot*c*tc))-\n",
+    "#              tc*df.inner(df.grad(c), df.grad((1-c_tot)/Ga0))+\n",
+    "#              tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c/Ga0))) * df.dx\n",
+    "    \n",
+    "    # Weak form radial symmetry:\n",
+    "    form = ((df.inner((c-c0)/dt, tc*X[0]*X[0]) +\n",
+    "             df.inner(df.grad(c), df.grad((1-c_tot)/Ga0*tc*X[0]*X[0]))) -\n",
+    "             df.inner(df.grad(c_tot), df.grad((1-c_tot)/Ga0/c_tot*c*tc*X[0]*X[0]))-\n",
+    "             tc*df.inner(df.grad(c), df.grad((1-c_tot)/Ga0*X[0]*X[0]))+\n",
+    "             tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c/Ga0*X[0]*X[0]))-\n",
+    "             (1-c_tot)/Ga0*2*X[0]*c.dx(0)*tc+\n",
+    "             (1-c_tot)/Ga0/c_tot*c*2*X[0]*c_tot.dx(0)*tc) * df.dx\n",
     "    t = 0\n",
     "    # Solve in time\n",
     "    ti = time.time()\n",
-    "    for i in range(6):\n",
-    "        print(time.time() - ti)\n",
+    "    for i in range(60):\n",
+    "#         print(time.time() - ti)\n",
     "        df.solve(form == 0, c)\n",
     "        df.assign(c0, c)\n",
     "        t += dt\n",
@@ -66,21 +77,21 @@
     "c0_1 = df.Function(F)\n",
     "c_tot_1 = df.Function(F)\n",
     "Ga0_1 = df.Function(F)\n",
-    "c_tot_1.interpolate(df.Expression('0*tanh(35000*(x[0]-0.01))+0.9', degree=1))\n",
+    "c_tot_1.interpolate(df.Expression('0*tanh(350000*(x[0]-0.01))+0.9', degree=1))\n",
     "c0_1.interpolate(df.Expression(('x[0]<0.01 ? 0 :'\n",
-    "                              '0*tanh(35000*(x[0]-0.01))+0.9'),\n",
+    "                              '0*tanh(350000*(x[0]-0.01))+0.9'),\n",
     "                             degree=1))\n",
-    "Ga0_1.interpolate(df.Expression('4.*(tanh(35000*(x[0]-0.01))+1)+1', degree=1))\n",
+    "Ga0_1.interpolate(df.Expression('4.*(tanh(350000*(x[0]-0.01))+1)+1', degree=1))\n",
     "\n",
     "# 1D, high partitioning\n",
     "c0_9 = df.Function(F)\n",
     "c_tot_9 = df.Function(F)\n",
     "Ga0_9 = df.Function(F)\n",
-    "c_tot_9.interpolate(df.Expression('0.4*tanh(-35000*(x[0]-0.01))+0.5', degree=1))\n",
+    "c_tot_9.interpolate(df.Expression('0.4*tanh(-350000*(x[0]-0.01))+0.5', degree=1))\n",
     "c0_9.interpolate(df.Expression(('x[0]<0.01 ? 0 :'\n",
-    "                              '0.4*tanh(-35000*(x[0]-0.01))+0.5'),\n",
+    "                              '0.4*tanh(-350000*(x[0]-0.01))+0.5'),\n",
     "                             degree=1))\n",
-    "Ga0_9.interpolate(df.Expression('0*(tanh(-35000*(x[0]-0.01))+1)+1', degree=1))\n",
+    "Ga0_9.interpolate(df.Expression('0*(tanh(-350000*(x[0]-0.01))+1)+1', degree=1))\n",
     "\n",
     "\n",
     "c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1)\n",
@@ -94,7 +105,7 @@
    "outputs": [],
    "source": [
     "# 1D:\n",
-    "plt.plot(np.linspace(0, 1, 10000), [c0_1([x]) for x in np.linspace(0, 1, 10000)])\n",
+    "plt.plot(np.linspace(0, 1, 10000), [1.5*c0_1([x]) for x in np.linspace(0, 1, 10000)])\n",
     "plt.plot(np.linspace(0, 1, 10000), [c0_9([x]) for x in np.linspace(0, 1, 10000)])\n",
     "plt.xlim(0, 0.1)\n",
     "# plt.ylim(0, 0.3)\n",
@@ -102,18 +113,6 @@
     "# 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": [
-    "# 1D:\n",
-    "plt.plot(np.linspace(0, 1, 10000), [c0([x]) for x in np.linspace(0, 1, 10000)])\n",
-    "# plt.xlim(0, 0.1)\n",
-    "# plt.ylim(0., 0.3)"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -146,8 +145,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import dolfin as df\n",
-    "\n",
     "mesh = df.UnitIntervalMesh(1000)\n",
     "dt = 0.001\n",
     "F = df.FunctionSpace(mesh, 'CG', 1)\n",
-- 
GitLab