Radially symmetric diffusion works.

Take care of volume elements everywhere! Needs to be in all integrals (also t) and mass conservation check.

Radially symmetric diffusion works.
2efcf4c0 · Lars Hubatsch · 3de621c6 · 2efcf4c0
Commit 2efcf4c0 authored 4 years ago by Lars Hubatsch
--- 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)])"
   ]
  },
  {

 %% 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, 50)
 mesh = df.UnitIntervalMesh(100000)
 dt = 0.000001

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

 %% Cell type:code id: tags:

 ``` python
 def calc_sim(c0, c_tot, Ga0):
    tc = df.TestFunction(F)
    c = df.Function(F)
    # Weak form
    form = ((df.inner((c-c0)/dt, tc) +
             df.inner(df.grad(c), df.grad((1-c_tot)/Ga0*tc))) -
             df.inner(df.grad(c_tot), df.grad((1-c_tot)/Ga0/c_tot*c*tc))-
             tc*df.inner(df.grad(c), df.grad((1-c_tot)/Ga0))+
             tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c/Ga0))) * df.dx
    t = 0
    # Solve in time
    ti = time.time()
    for i in range(6):
        print(time.time() - ti)
        df.solve(form == 0, c)
        df.assign(c0, c)
        t += dt
    print(time.time() - ti)
    return c0
 ```

 %% Cell type:code id: tags:

 ``` python
 # Interpolate c_tot and initial conditions
 # 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, no partitioning
 c0_1 = df.Function(F)
 c_tot_1 = df.Function(F)
 Ga0_1 = df.Function(F)
 c_tot_1.interpolate(df.Expression('0*tanh(35000*(x[0]-0.01))+0.9', degree=1))
 c0_1.interpolate(df.Expression(('x[0]<0.01 ? 0 :'
                              '0*tanh(35000*(x[0]-0.01))+0.9'),
                             degree=1))
 Ga0_1.interpolate(df.Expression('4.*(tanh(35000*(x[0]-0.01))+1)+1', degree=1))

 # 1D, high partitioning
 c0_9 = df.Function(F)
 c_tot_9 = df.Function(F)
 Ga0_9 = df.Function(F)
 c_tot_9.interpolate(df.Expression('0.4*tanh(-35000*(x[0]-0.01))+0.5', degree=1))
 c0_9.interpolate(df.Expression(('x[0]<0.01 ? 0 :'
                              '0.4*tanh(-35000*(x[0]-0.01))+0.5'),
                             degree=1))
 Ga0_9.interpolate(df.Expression('0*(tanh(-35000*(x[0]-0.01))+1)+1', degree=1))


 c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1)
 c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9)
 ```

 %% Cell type:code id: tags:

 ``` python
 # 1D:
 plt.plot(np.linspace(0, 1, 10000), [c0_1([x]) for x in np.linspace(0, 1, 10000)])
 plt.plot(np.linspace(0, 1, 10000), [c0_9([x]) for x in np.linspace(0, 1, 10000)])
 plt.xlim(0, 0.1)
 # plt.ylim(0, 0.3)
 # 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
 # 1D:
 plt.plot(np.linspace(0, 1, 10000), [c0([x]) for x in np.linspace(0, 1, 10000)])
 # plt.xlim(0, 0.1)
 # plt.ylim(0., 0.3)
 ```

 %% Cell type:code id: tags:

 ``` python
 plt.plot(np.linspace(0, 1, 2000), [Ga0_1([x]) for x in np.linspace(0, 1, 2000)])
 plt.xlim(0, 0.1)
 ```

 %% Cell type:code id: tags:

 ``` python
 plt.plot(np.linspace(0, 1, 1000), [c_tot_9([x]) for x in np.linspace(0, 1, 1000)])
 ```

 %% Cell type:markdown id: tags:

 ## Radial diffusion equation

 %% Cell type:code id: tags:

 ``` python
 import dolfin as df

 mesh = df.UnitIntervalMesh(1000)
-dt = 0.01
+dt = 0.001
 F = df.FunctionSpace(mesh, 'CG', 1)
 c0 = df.Function(F)
-c0.interpolate(df.Expression('x[0]<0.5 ? 0:1', degree=1))
+c0.interpolate(df.Expression('x[0]<0.5 && x[0]>0.2 ? 1:0', degree=1))
 q = df.TestFunction(F)
 c = df.Function(F)
 X = df.SpatialCoordinate(mesh)
-
+g = df.Expression('.00', degree=1)
+u_D = df.Expression('1', degree=1)
+def boundary(x, on_boundary):
+    return on_boundary
+bc = df.DirichletBC(F, u_D, boundary)
 # Weak form spherical symmetry
-form = (df.inner((c-c0)/dt, q) +
-        df.inner(df.grad(c), df.grad(4*3.1415*X[0]*X[0]*q))-
-        c.dx(0)*8*3.1415*X[0]*q) * df.dx
+form = (df.inner((c-c0)/dt, q*X[0]*X[0]) +
+        df.inner(df.grad(c), df.grad(X[0]*X[0]*q))-
+        c.dx(0)*2*X[0]*q) * df.dx
 # Weak form 1D
 # form = (df.inner((c-c0)/dt, q) +
 #         df.inner(df.grad(c), df.grad(q))) * df.dx
 t = 0

 # Solve in time
-for i in range(6):
+for i in range(60):
+    print(np.sum([x*x*c0([x]) for x in np.linspace(0, 1, 1000)]))
    df.solve(form == 0, c)
    df.assign(c0, c)
    t += dt
+    plt.plot(np.linspace(0, 1, 1000), [c0([x]) for x in np.linspace(0, 1, 1000)])
 ```

 %% Cell type:code id: tags:

 ``` python
 ```