Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import dolfin as df\n",
"import matplotlib.pyplot as plt\n",
"import mshr as ms\n",
"import numpy as np\n",
"import time\n",
"df.set_log_level(40)\n",
"\n",
"domain = ms.Sphere(df.Point(0, 0, 0), 1.0)\n",
"# mesh = ms.generate_mesh(domain, 35)\n",
"mesh = df.UnitSquareMesh(50,50)\n",
"dt = 0.01\n",
"\n",
"F = df.FunctionSpace(mesh, 'CG', 1)\n",
"\n",
"c0 = df.Function(F)\n",
"c_tot = df.Function(F)\n",
"c = df.Function(F)\n",
"tc = df.TestFunction(F)\n",
"\n",
"# Interpolate c_tot and initial conditions\n",
"# c_tot.interpolate(df.Expression('-0.4*tanh((sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/0.1-0.5)) ', degree=1))\n",
"c_tot.interpolate(df.Expression('0.4*tanh(-35*(sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))-0.2)) +0.5', degree=1))\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",
"# c0.interpolate(df.Expression('(x[0]<0.5) ? 0 : 1.0', degree=1))\n",
"# Somewhat realistic half FRAP\n",
"c0.interpolate(df.Expression(('(x[0]<0.5) && sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))<0.2 ? 0 :'\n",
" '0.4*tanh(-35*(sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))-0.2)) +0.5'),\n",
" degree=1))\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",
"\n",
"# Open file for writing\n",
"cFile = df.XDMFFile('conc.xdmf')\n",
"t = 0\n",
"cFile.write(c0, t)\n",
"\n",
"# Solve in time\n",
"ti = time.time()\n",
"for i in range(80):\n",
" print(t)\n",
" df.solve(form == 0, c)\n",
" df.assign(c0, c)\n",
" t += dt\n",
" cFile.write(c0, t)\n",
"cFile.close()\n",
"print(time.time() - ti)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# c_tot.interpolate(df.Expression('-0.4*tanh((sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/0.1-0.1)) ', degree=1))\n",
"c_tot.interpolate(df.Expression('0.4*tanh(-35*(sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))-0.2)) +0.5', degree=1))\n",
"plt.plot(np.linspace(0, 1, 100), [c_tot([x, 0.5]) for x in np.linspace(0, 1, 100)])"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}