From 2efcf4c0af16f1a730b97965489b2ff82549f653 Mon Sep 17 00:00:00 2001
From: Lars Hubatsch <hubatsch@pks.mpg.de>
Date: Wed, 16 Sep 2020 15:48:01 +0200
Subject: [PATCH] Radially symmetric diffusion works.
Take care of volume elements everywhere! Needs to be in all integrals (also t) and mass conservation check.
---
FloryHugg_DiffUnbleached.ipynb | 22 ++++++++++++++--------
1 file changed, 14 insertions(+), 8 deletions(-)
diff --git a/FloryHugg_DiffUnbleached.ipynb b/FloryHugg_DiffUnbleached.ipynb
index ad6d433..87c67f8 100644
--- a/FloryHugg_DiffUnbleached.ipynb
+++ b/FloryHugg_DiffUnbleached.ipynb
@@ -149,28 +149,34 @@
"import dolfin as df\n",
"\n",
"mesh = df.UnitIntervalMesh(1000)\n",
- "dt = 0.01\n",
+ "dt = 0.001\n",
"F = df.FunctionSpace(mesh, 'CG', 1)\n",
"c0 = df.Function(F)\n",
- "c0.interpolate(df.Expression('x[0]<0.5 ? 0:1', degree=1))\n",
+ "c0.interpolate(df.Expression('x[0]<0.5 && x[0]>0.2 ? 1:0', degree=1))\n",
"q = df.TestFunction(F)\n",
"c = df.Function(F)\n",
"X = df.SpatialCoordinate(mesh)\n",
- "\n",
+ "g = df.Expression('.00', degree=1)\n",
+ "u_D = df.Expression('1', degree=1)\n",
+ "def boundary(x, on_boundary):\n",
+ " return on_boundary\n",
+ "bc = df.DirichletBC(F, u_D, boundary)\n",
"# Weak form spherical symmetry\n",
- "form = (df.inner((c-c0)/dt, q) +\n",
- " df.inner(df.grad(c), df.grad(4*3.1415*X[0]*X[0]*q))-\n",
- " c.dx(0)*8*3.1415*X[0]*q) * df.dx\n",
+ "form = (df.inner((c-c0)/dt, q*X[0]*X[0]) +\n",
+ " df.inner(df.grad(c), df.grad(X[0]*X[0]*q))-\n",
+ " c.dx(0)*2*X[0]*q) * df.dx\n",
"# Weak form 1D\n",
"# form = (df.inner((c-c0)/dt, q) +\n",
"# df.inner(df.grad(c), df.grad(q))) * df.dx\n",
"t = 0\n",
"\n",
"# Solve in time\n",
- "for i in range(6):\n",
+ "for i in range(60):\n",
+ " print(np.sum([x*x*c0([x]) for x in np.linspace(0, 1, 1000)]))\n",
" df.solve(form == 0, c)\n",
" df.assign(c0, c)\n",
- " t += dt"
+ " t += dt\n",
+ " plt.plot(np.linspace(0, 1, 1000), [c0([x]) for x in np.linspace(0, 1, 1000)])"
]
},
{
--
GitLab