Skip to content
Snippets Groups Projects
Commit d74d736a authored by Lars Hubatsch's avatar Lars Hubatsch
Browse files

Change cell order. Works OK for small times. 

Probably need proper dependence on phi_tot instead of direct dependence on x.
parent 7c35122f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
FloryHugg_DiffUnbleached.ipynb
+ 64
− 60
View file @ d74d736a
%% 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(10000)
dt = 0.0001
dt = 0.001
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)
X = df.SpatialCoordinate(mesh)
# X.interpolate(df.Expression('x[0]', degree=1))
# Weak form 1D:
# 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
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
# # Weak form radial symmetry:
form = ((df.inner((c-c0)/dt, tc*X[0]*X[0]) +
df.inner(df.grad(c), df.grad((1-c_tot)*Ga0*tc*X[0]*X[0]))) -
df.inner(df.grad(c_tot), df.grad((1-c_tot)*Ga0/c_tot*c*tc*X[0]*X[0]))-
tc*df.inner(df.grad(c), df.grad((1-c_tot)*Ga0*X[0]*X[0]))+
tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c*Ga0*X[0]*X[0]))-
(1-c_tot)*Ga0*2*X[0]*c.dx(0)*tc+
(1-c_tot)*Ga0/c_tot*c*2*X[0]*c_tot.dx(0)*tc) * df.dx
# form = ((df.inner((c-c0)/dt, tc*X[0]*X[0]) +
# df.inner(df.grad(c), df.grad((1-c_tot)*Ga0*tc*X[0]*X[0]))) -
# df.inner(df.grad(c_tot), df.grad((1-c_tot)*Ga0/c_tot*c*tc*X[0]*X[0]))-
# tc*df.inner(df.grad(c), df.grad((1-c_tot)*Ga0*X[0]*X[0]))+
# tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c*Ga0*X[0]*X[0]))-
# (1-c_tot)*Ga0*2*X[0]*c.dx(0)*tc+
# (1-c_tot)*Ga0/c_tot*c*2*X[0]*c_tot.dx(0)*tc) * df.dx
# Weak form radial symmetry:
# form = ((df.inner((c-c0)/dt, tc*X[0]*X[0]) +
# df.inner(df.grad(c), df.grad((1-c_tot+a*c_tot*c_tot)*tc*X[0]*X[0]))) -
# df.inner(df.grad(c_tot), df.grad((1-c_tot+a*c_tot*c_tot)/c_tot*c*tc*X[0]*X[0]))-
# tc*df.inner(df.grad(c), df.grad((1-c_tot+a*c_tot*c_tot)*X[0]*X[0]))+
# tc*df.inner(df.grad(c_tot), df.grad((1-c_tot+a*c_tot*c_tot)/c_tot*c*X[0]*X[0]))-
# (1-c_tot+a*c_tot*c_tot)*2*X[0]*c.dx(0)*tc+
# (1-c_tot+a*c_tot*c_tot)/c_tot*c*2*X[0]*c_tot.dx(0)*tc) * df.dx
# df.inner(df.grad(c), df.grad((1-c_tot+Ga0*c_tot*c_tot)*tc*X[0]*X[0]))) -
# df.inner(df.grad(c_tot), df.grad((1-c_tot+Ga0*c_tot*c_tot)/c_tot*c*tc*X[0]*X[0]))-
# tc*df.inner(df.grad(c), df.grad((1-c_tot+Ga0*c_tot*c_tot)*X[0]*X[0]))+
# tc*df.inner(df.grad(c_tot), df.grad((1-c_tot+Ga0*c_tot*c_tot)/c_tot*c*X[0]*X[0]))-
# (1-c_tot+Ga0*c_tot*c_tot)*2*X[0]*c.dx(0)*tc+
# (1-c_tot+Ga0*c_tot*c_tot)/c_tot*c*2*X[0]*c_tot.dx(0)*tc) * df.dx
t = 0
# Solve in time
ti = time.time()
for i in range(100):
# 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(350000*(x[0]-0.01))+0.9', degree=1))
c_tot_1.interpolate(df.Expression('0.9', degree=1))
c0_1.interpolate(df.Expression(('x[0]<0.1 ? 0 :'
'0*tanh(350000*(x[0]-0.01))+0.9'),
degree=1))
# Ga0_1.interpolate(df.Expression('4.*(tanh(350000*(x[0]-0.01))+1)+1', degree=1))
Ga0_1.interpolate(df.Expression('x[0]<0.1 ? 1:9',
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(-350000*(x[0]-0.01))+0.5', degree=1))
c_tot_9.interpolate(df.Expression('x[0]<0.1 ? 0.9 :0.1', degree=1))
# c0_9.interpolate(df.Expression(('x[0]<0.01 ? 0 :'
# '0.4*tanh(-350000*(x[0]-0.01))+0.5'),
# degree=1))
c0_9.interpolate(df.Expression('x[0]<0.1 ? 0 :0.1',
degree=1))
# Ga0_9.interpolate(df.Expression('0*(tanh(-350000*(x[0]-0.01))+1)+1', degree=1))
Ga0_9.interpolate(df.Expression('1', degree=1))
```
%% Cell type:code id: tags:
``` python
p1_i = 0.9
p1_o = 0.1
p2_i = 0.8
p2_o = 0.2
a1 = 0
def p_tot(p_i, p_o):
return str(p_i-p_o)+'*(-0.5*tanh(3500*(x[0]-0.1))+0.5)+'+str(p_o)
ct_1 = df.Function(F)
ct_2 = df.Function(F)
c0_1 = df.Function(F)
c0_2 = df.Function(F)
g_1 = df.Function(F)
g_2 = df.Function(F)
ct_1.interpolate(df.Expression(p_tot(p1_i, p1_o), degree=1))
ct_2.interpolate(df.Expression(p_tot(p2_i, p2_o), degree=1))
P1 = c_tot1(0)/c_tot1(1)
P2 = c_tot2(0)/c_tot2(1)
g_1.interpolate(df.Expression('1', degree=1))
g_2.interpolate(df.Expression(p_tot(p1_o/p2_o, 1), degree=1))
D_out1 = 1-c_tot1(1)
a = (P2*D_out1/P1-1+p2_o)/p2_o**2
c0_1.interpolate(df.Expression('x[0]<0.1 ? 0 :'+p_tot(p1_i, p1_o), degree=1))
c0_2.interpolate(df.Expression('x[0]<0.1 ? 0 :'+p_tot(p2_i, p2_o), degree=1))
```
%% Cell type:code id: tags:
``` python
# c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1)
# c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9)
# c0_1 = calc_sim(c0_1, c_tot1, 0)
# c0_2 = calc_sim(c0_2, c_tot2, a)
c0_1 = calc_sim(c0_1, ct_1, g_1)
c0_2 = calc_sim(c0_2, ct_2, g_2)
```
%% Cell type:code id: tags:
``` python
# 1D:
plt.plot(np.linspace(0, 1, 10000), [1.31*c0_1([x]) for x in np.linspace(0, 1, 10000)])
plt.plot(np.linspace(0, 1, 10000), [1.465*c0_1([x]) for x in np.linspace(0, 1, 10000)])
plt.plot(np.linspace(0, 1, 10000), [c0_2([x]) for x in np.linspace(0, 1, 10000)])
plt.xlim(0.06, 0.125)
plt.xlim(0.0, 0.125)
# 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
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:code id: tags:
``` python
p1_i = 0.9
p1_o = 0.1
p2_i = 0.8
p2_o = 0.2
a1 = 0
def p_tot(p_i, p_o):
return str(p_i-p_o)+'*(-0.5*tanh(350000*(x[0]-0.1))+0.5)+'+str(p_o)
ct_1 = df.Function(F)
ct_2 = df.Function(F)
c0_1 = df.Function(F)
c0_2 = df.Function(F)
g_1 = df.Function(F)
g_2 = df.Function(F)
ct_1.interpolate(df.Expression(p_tot(p1_i, p1_o), degree=1))
ct_2.interpolate(df.Expression(p_tot(p2_i, p2_o), degree=1))
P1 = c_tot1(0)/c_tot1(1)
P2 = c_tot2(0)/c_tot2(1)
g_1.interpolate(df.Expression('1', degree=1))
g_2.interpolate(df.Expression(p_tot(p1_o/p2_o, p1_i/p2_i), degree=1))
D_out1 = 1-c_tot1(1)
a = (P2*D_out1/P1-1+p2_o)/p2_o**2
c0_1.interpolate(df.Expression('x[0]<0.1 ? 0 :'+p_tot(p1_i, p1_o), degree=1))
c0_2.interpolate(df.Expression('x[0]<0.1 ? 0 :'+p_tot(p2_i, p2_o), degree=1))
```
%% Cell type:code id: tags:
``` python
plt.plot(np.linspace(0, 1, 1000), [ct_1([x]) for x in np.linspace(0, 1, 1000)])
plt.plot(np.linspace(0, 1, 1000), [ct_2([x]) for x in np.linspace(0, 1, 1000)])
plt.plot(np.linspace(0, 1, 1000), [(1-ct_1([x]))*ct_1([x]) for x in np.linspace(0, 1, 1000)])
plt.plot(np.linspace(0, 1, 1000), [(1-ct_2([x]))*g_2([x]) for x in np.linspace(0, 1, 1000)])
plt.show()
plt.plot(np.linspace(0, 1, 1000), [g_2([x]) for x in np.linspace(0, 1, 1000)])
# plt.plot(np.linspace(0, 1, 1000), [(1-c_tot1([x])-a*c_tot1([x])**2) for x in np.linspace(0, 1, 1000)])
```
%% Cell type:markdown id: tags:
## Radial diffusion equation
%% Cell type:code id: tags:
``` python
mesh = df.UnitIntervalMesh(1000)
dt = 0.001
F = df.FunctionSpace(mesh, 'CG', 1)
c0 = df.Function(F)
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*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(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
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment