From 21570d0dd2815aa05290498163d1553cb52e11e5 Mon Sep 17 00:00:00 2001
From: Lars Hubatsch <hubatsch@pks.mpg.de>
Date: Tue, 15 Sep 2020 11:53:44 +0200
Subject: [PATCH] Swapping Gamma and D works in principle for 1D.
---
FloryHugg_DiffUnbleached.ipynb | 56 ++++++++++++++++++++++++++--------
1 file changed, 44 insertions(+), 12 deletions(-)
diff --git a/FloryHugg_DiffUnbleached.ipynb b/FloryHugg_DiffUnbleached.ipynb
index 0ab22ba..3b5c4c2 100644
--- a/FloryHugg_DiffUnbleached.ipynb
+++ b/FloryHugg_DiffUnbleached.ipynb
@@ -15,8 +15,8 @@
"\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",
+ "mesh = df.UnitIntervalMesh(10000)\n",
+ "dt = 0.0000002\n",
"\n",
"F = df.FunctionSpace(mesh, 'CG', 1)\n",
"\n",
@@ -24,18 +24,27 @@
"c_tot = df.Function(F)\n",
"c = df.Function(F)\n",
"tc = df.TestFunction(F)\n",
- "\n",
+ "Ga0 = df.Function(F)\n",
"# Interpolate c_tot and initial conditions\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",
+ "# 1D, high partitioning\n",
+ "# c_tot.interpolate(df.Expression('0.4*tanh(-35000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001))+0.5', degree=1))\n",
+ "# c0.interpolate(df.Expression(('sqrt((x[0]-0.5)*(x[0]-0.5))<0.001 ? 0 :'\n",
+ "# '0.4*tanh(-35000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001)) +0.5'),\n",
+ "# degree=1))\n",
+ "# Ga0.interpolate(df.Expression('0*(tanh(-35000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001))+1)+1', degree=1))\n",
+ "\n",
+ "# 1D, no partitioning\n",
+ "c_tot.interpolate(df.Expression('0*tanh(-35000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001))+0.9', degree=1))\n",
+ "c0.interpolate(df.Expression(('sqrt((x[0]-0.5)*(x[0]-0.5))<0.001 ? 0 :'\n",
+ " '0*tanh(-35000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001)) +0.9'),\n",
" degree=1))\n",
+ "Ga0.interpolate(df.Expression('4.*(tanh(35000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001))+1)+1', degree=1))\n",
+ "\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",
@@ -43,8 +52,8 @@
"# Somewhat realistic half FRAP\n",
"\n",
"# Weak form\n",
- "form = ((df.inner((c-c0)/dt, tc) + (1-c_tot)* df.inner(df.grad(c), df.grad(tc))) -\n",
- " (1-c_tot)/c_tot* c*df.inner(df.grad(c_tot), df.grad(tc))) * df.dx\n",
+ "form = ((df.inner((c-c0)/dt, tc) + 1/Ga0*(1-c_tot)* df.inner(df.grad(c), df.grad(tc))) -\n",
+ " 1/Ga0*(1-c_tot)/c_tot* c*df.inner(df.grad(c_tot), df.grad(tc))) * df.dx\n",
"\n",
"# Open file for writing\n",
"cFile = df.XDMFFile('conc.xdmf')\n",
@@ -53,7 +62,7 @@
"\n",
"# Solve in time\n",
"ti = time.time()\n",
- "for i in range(8):\n",
+ "for i in range(10):\n",
" print(time.time() - ti)\n",
" df.solve(form == 0, c)\n",
" df.assign(c0, c)\n",
@@ -70,7 +79,9 @@
"outputs": [],
"source": [
"# 1D:\n",
- "plt.plot(np.linspace(0, 1, 1000), [c0([x]) for x in np.linspace(0, 1, 1000)])\n",
+ "plt.plot(np.linspace(0, 1, 10000), [c0([x]) for x in np.linspace(0, 1, 10000)])\n",
+ "plt.xlim(0.495, 0.505)\n",
+ "plt.ylim(0, 0.3)\n",
"# 3D:\n",
"# plt.plot(np.linspace(0, 0.5, 1000), [c0([x, 0, 0]) for x in np.linspace(0, 0.5, 1000)])"
]
@@ -81,7 +92,28 @@
"metadata": {},
"outputs": [],
"source": [
- "c0([0, 0, 0])"
+ "# 1D:\n",
+ "plt.plot(np.linspace(0, 1, 10000), [c0([x]) for x in np.linspace(0, 1, 10000)])\n",
+ "plt.xlim(0.495, 0.505)\n",
+ "plt.ylim(0., 0.3)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(np.linspace(0, 1, 2000), [Ga0([x]) for x in np.linspace(0, 1, 2000)])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(np.linspace(0, 1, 1000), [c_tot([x]) for x in np.linspace(0, 1, 1000)])"
]
},
{
--
GitLab