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

1D works again. Much more general functions. 

3D can be rescaled by maximum. 1D seems perfect up until material conservation kicks in at boundary.
parent d74d736a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
FloryHugg_DiffUnbleached.ipynb
+ 56
− 123
View file @ 0cea5929
%% 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.001
F = df.FunctionSpace(mesh, 'CG', 1)
```
%% Cell type:code id: tags:
``` python
def calc_sim(c0, c_tot, Ga0):
def calc_sim(c0, c_tot, Ga0, dt, n_t, sym):
tc = df.TestFunction(F)
c = df.Function(F)
X = df.SpatialCoordinate(mesh)
# X.interpolate(df.Expression('x[0]', degree=1))
if sym == 0:
# Weak form 1D:
form = ((df.inner((c-c0)/dt, tc) +
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
elif sym == 2:
# 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
# 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*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)
for i in range(n_t):
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)
def create_func(f_space, expr_str, deg):
f = df.Function(f_space)
f.interpolate(df.Expression(expr_str, degree=deg))
return f
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))
def eval_func(func, x):
return [func([x]) for x in x]
```
%% 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)
c_tot_1 = create_func(F, p_tot(0.9, 0.9), 1)
c0_1 = create_func(F, 'x[0]<0.1 ? 0 :' + p_tot(0.9, 0.9), 1)
Ga0_1 = create_func(F, p_tot(1, 1/9), 1)
# c0_1 = calc_sim(c0_1, c_tot1, 0)
# c0_2 = calc_sim(c0_2, c_tot2, a)
c_tot_9 = create_func(F, p_tot(0.9, 0.1), 1)
c0_9 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(0.9, 0.1), 1)
Ga0_9 = create_func(F,p_tot(1, 1), 1)
c0_1 = calc_sim(c0_1, ct_1, g_1)
c0_2 = calc_sim(c0_2, ct_2, g_2)
c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1, 0.001, 10, 0)
c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9, 0.001, 10, 0)
```
%% Cell type:code id: tags:
``` python
# 1D:
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.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)])
x = np.linspace(0, 1, 10000)
plt.plot(x, eval_func(c0_1, x))
plt.plot(x, eval_func(c0_9, x))
plt.xlim(0, 0.2)
```
%% 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)
```
p1_i = 0.9; p1_o = 0.1
p2_i = 0.8; p2_o = 0.2
p3_i = 0.9; p3_o = 0.9
%% Cell type:code id: tags:
ct_1 = create_func(F, p_tot(p1_i, p1_o), 1)
ct_2 = create_func(F, p_tot(p2_i, p2_o), 1)
ct_3 = create_func(F, p_tot(p3_i, p3_o), 1)
``` python
plt.plot(np.linspace(0, 1, 1000), [c_tot_9([x]) for x in np.linspace(0, 1, 1000)])
g_1= create_func(F, '1', 1)
g_2= create_func(F, p_tot(p1_o/p2_o, 1), 1)
g_3= create_func(F, p_tot(1, 1), 1)
c0_1 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(p1_i, p1_o), 1)
c0_2 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(p2_i, p2_o), 1)
c0_3 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(p3_i, p3_o), 1)
c0_1 = calc_sim(c0_1, ct_1, g_1, 0.001, 50, 0)
c0_2 = calc_sim(c0_2, ct_2, g_2, 0.001, 50, 0)
c0_3 = calc_sim(c0_3, ct_3, g_3, 0.001, 50, 0)
```
%% Cell type:code id: tags:
``` python
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)])
plt.plot(x, eval_func(c0_1, x)/np.max(eval_func(c0_1, x)))
plt.plot(x, eval_func(c0_2, x)/np.max(eval_func(c0_2, x)))
plt.plot(x, eval_func(c0_3, x)/np.array(0.47))
plt.xlim(0.0, 0.125)
plt.ylim(0, 1.2)
```
%% 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