From 0cea5929464a1a7c40260e0d3ab30f4d52f9ae65 Mon Sep 17 00:00:00 2001
From: Lars Hubatsch <hubatsch@pks.mpg.de>
Date: Tue, 22 Sep 2020 15:22:57 +0200
Subject: [PATCH] 1D works again. Much more general functions.
3D can be rescaled by maximum. 1D seems perfect up until material conservation kicks in at boundary.
---
FloryHugg_DiffUnbleached.ipynb | 179 +++++++++++----------------------
1 file changed, 56 insertions(+), 123 deletions(-)
diff --git a/FloryHugg_DiffUnbleached.ipynb b/FloryHugg_DiffUnbleached.ipynb
index d631829..021188e 100644
--- a/FloryHugg_DiffUnbleached.ipynb
+++ b/FloryHugg_DiffUnbleached.ipynb
@@ -16,8 +16,6 @@
"# domain = ms.Sphere(df.Point(0, 0, 0), 1.0)\n",
"# mesh = ms.generate_mesh(domain, 50)\n",
"mesh = df.UnitIntervalMesh(10000)\n",
- "dt = 0.001\n",
- "\n",
"F = df.FunctionSpace(mesh, 'CG', 1)"
]
},
@@ -27,26 +25,27 @@
"metadata": {},
"outputs": [],
"source": [
- "def calc_sim(c0, c_tot, Ga0):\n",
+ "def calc_sim(c0, c_tot, Ga0, dt, n_t, sym):\n",
" tc = df.TestFunction(F)\n",
" c = df.Function(F)\n",
" X = df.SpatialCoordinate(mesh)\n",
"# X.interpolate(df.Expression('x[0]', degree=1))\n",
+ " if sym == 0:\n",
" # Weak form 1D:\n",
- " form = ((df.inner((c-c0)/dt, tc) +\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",
+ " elif sym == 2:\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",
" # 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*c_tot*c_tot)*tc*X[0]*X[0]))) -\n",
@@ -59,95 +58,23 @@
" t = 0\n",
" # Solve in time\n",
" ti = time.time()\n",
- " for i in range(100):\n",
- "# print(time.time() - ti)\n",
+ " for i in range(n_t):\n",
" df.solve(form == 0, c)\n",
" df.assign(c0, c)\n",
" t += dt\n",
" print(time.time() - ti)\n",
- " return c0"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "jupyter": {
- "source_hidden": true
- }
- },
- "outputs": [],
- "source": [
- "# 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",
- "\n",
- "\n",
- "# 1D, no partitioning\n",
- "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(350000*(x[0]-0.01))+0.9', degree=1))\n",
- "c_tot_1.interpolate(df.Expression('0.9', degree=1))\n",
- "c0_1.interpolate(df.Expression(('x[0]<0.1 ? 0 :'\n",
- " '0*tanh(350000*(x[0]-0.01))+0.9'),\n",
- " degree=1))\n",
- "# Ga0_1.interpolate(df.Expression('4.*(tanh(350000*(x[0]-0.01))+1)+1', degree=1))\n",
- "Ga0_1.interpolate(df.Expression('x[0]<0.1 ? 1:9',\n",
- " 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(-350000*(x[0]-0.01))+0.5', degree=1))\n",
- "c_tot_9.interpolate(df.Expression('x[0]<0.1 ? 0.9 :0.1', degree=1))\n",
- "# c0_9.interpolate(df.Expression(('x[0]<0.01 ? 0 :'\n",
- "# '0.4*tanh(-350000*(x[0]-0.01))+0.5'),\n",
- "# degree=1))\n",
- "c0_9.interpolate(df.Expression('x[0]<0.1 ? 0 :0.1',\n",
- " degree=1))\n",
- "# Ga0_9.interpolate(df.Expression('0*(tanh(-350000*(x[0]-0.01))+1)+1', degree=1))\n",
- "Ga0_9.interpolate(df.Expression('1', degree=1))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "p1_i = 0.9\n",
- "p1_o = 0.1\n",
- "p2_i = 0.8\n",
- "p2_o = 0.2\n",
- "a1 = 0\n",
+ " return c0\n",
"\n",
"def p_tot(p_i, p_o):\n",
" return str(p_i-p_o)+'*(-0.5*tanh(3500*(x[0]-0.1))+0.5)+'+str(p_o)\n",
"\n",
- "ct_1 = df.Function(F)\n",
- "ct_2 = df.Function(F)\n",
- "c0_1 = df.Function(F)\n",
- "c0_2 = df.Function(F)\n",
- "g_1 = df.Function(F)\n",
- "g_2 = df.Function(F)\n",
- "\n",
- "ct_1.interpolate(df.Expression(p_tot(p1_i, p1_o), degree=1))\n",
- "ct_2.interpolate(df.Expression(p_tot(p2_i, p2_o), degree=1))\n",
- "P1 = c_tot1(0)/c_tot1(1)\n",
- "P2 = c_tot2(0)/c_tot2(1)\n",
- "\n",
- "g_1.interpolate(df.Expression('1', degree=1))\n",
- "g_2.interpolate(df.Expression(p_tot(p1_o/p2_o, 1), degree=1))\n",
- "D_out1 = 1-c_tot1(1)\n",
- "a = (P2*D_out1/P1-1+p2_o)/p2_o**2\n",
+ "def create_func(f_space, expr_str, deg):\n",
+ " f = df.Function(f_space)\n",
+ " f.interpolate(df.Expression(expr_str, degree=deg))\n",
+ " return f\n",
"\n",
- "c0_1.interpolate(df.Expression('x[0]<0.1 ? 0 :'+p_tot(p1_i, p1_o), degree=1))\n",
- "c0_2.interpolate(df.Expression('x[0]<0.1 ? 0 :'+p_tot(p2_i, p2_o), degree=1))"
+ "def eval_func(func, x):\n",
+ " return [func([x]) for x in x]"
]
},
{
@@ -156,29 +83,16 @@
"metadata": {},
"outputs": [],
"source": [
- "# c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1)\n",
- "# c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9)\n",
+ "c_tot_1 = create_func(F, p_tot(0.9, 0.9), 1)\n",
+ "c0_1 = create_func(F, 'x[0]<0.1 ? 0 :' + p_tot(0.9, 0.9), 1)\n",
+ "Ga0_1 = create_func(F, p_tot(1, 1/9), 1)\n",
"\n",
- "# c0_1 = calc_sim(c0_1, c_tot1, 0)\n",
- "# c0_2 = calc_sim(c0_2, c_tot2, a)\n",
+ "c_tot_9 = create_func(F, p_tot(0.9, 0.1), 1)\n",
+ "c0_9 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(0.9, 0.1), 1)\n",
+ "Ga0_9 = create_func(F,p_tot(1, 1), 1)\n",
"\n",
- "c0_1 = calc_sim(c0_1, ct_1, g_1)\n",
- "c0_2 = calc_sim(c0_2, ct_2, g_2)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# 1D:\n",
- "plt.plot(np.linspace(0, 1, 10000), [1.465*c0_1([x]) for x in np.linspace(0, 1, 10000)])\n",
- "plt.plot(np.linspace(0, 1, 10000), [c0_2([x]) for x in np.linspace(0, 1, 10000)])\n",
- "plt.xlim(0.0, 0.125)\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)])"
+ "c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1, 0.001, 10, 0)\n",
+ "c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9, 0.001, 10, 0)"
]
},
{
@@ -187,8 +101,10 @@
"metadata": {},
"outputs": [],
"source": [
- "plt.plot(np.linspace(0, 1, 2000), [Ga0_1([x]) for x in np.linspace(0, 1, 2000)])\n",
- "plt.xlim(0, 0.1)"
+ "x = np.linspace(0, 1, 10000)\n",
+ "plt.plot(x, eval_func(c0_1, x))\n",
+ "plt.plot(x, eval_func(c0_9, x))\n",
+ "plt.xlim(0, 0.2)"
]
},
{
@@ -197,7 +113,25 @@
"metadata": {},
"outputs": [],
"source": [
- "plt.plot(np.linspace(0, 1, 1000), [c_tot_9([x]) for x in np.linspace(0, 1, 1000)])"
+ "p1_i = 0.9; p1_o = 0.1\n",
+ "p2_i = 0.8; p2_o = 0.2\n",
+ "p3_i = 0.9; p3_o = 0.9\n",
+ "\n",
+ "ct_1 = create_func(F, p_tot(p1_i, p1_o), 1)\n",
+ "ct_2 = create_func(F, p_tot(p2_i, p2_o), 1)\n",
+ "ct_3 = create_func(F, p_tot(p3_i, p3_o), 1)\n",
+ "\n",
+ "g_1= create_func(F, '1', 1)\n",
+ "g_2= create_func(F, p_tot(p1_o/p2_o, 1), 1)\n",
+ "g_3= create_func(F, p_tot(1, 1), 1)\n",
+ "\n",
+ "c0_1 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(p1_i, p1_o), 1)\n",
+ "c0_2 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(p2_i, p2_o), 1)\n",
+ "c0_3 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(p3_i, p3_o), 1)\n",
+ "\n",
+ "c0_1 = calc_sim(c0_1, ct_1, g_1, 0.001, 50, 0)\n",
+ "c0_2 = calc_sim(c0_2, ct_2, g_2, 0.001, 50, 0)\n",
+ "c0_3 = calc_sim(c0_3, ct_3, g_3, 0.001, 50, 0)"
]
},
{
@@ -206,12 +140,11 @@
"metadata": {},
"outputs": [],
"source": [
- "plt.plot(np.linspace(0, 1, 1000), [(1-ct_1([x]))*ct_1([x]) for x in np.linspace(0, 1, 1000)])\n",
- "plt.plot(np.linspace(0, 1, 1000), [(1-ct_2([x]))*g_2([x]) for x in np.linspace(0, 1, 1000)])\n",
- "plt.show()\n",
- "\n",
- "plt.plot(np.linspace(0, 1, 1000), [g_2([x]) for x in np.linspace(0, 1, 1000)])\n",
- "# plt.plot(np.linspace(0, 1, 1000), [(1-c_tot1([x])-a*c_tot1([x])**2) for x in np.linspace(0, 1, 1000)])"
+ "plt.plot(x, eval_func(c0_1, x)/np.max(eval_func(c0_1, x)))\n",
+ "plt.plot(x, eval_func(c0_2, x)/np.max(eval_func(c0_2, x)))\n",
+ "plt.plot(x, eval_func(c0_3, x)/np.array(0.47))\n",
+ "plt.xlim(0.0, 0.125)\n",
+ "plt.ylim(0, 1.2)"
]
},
{
--
GitLab