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

Found third value for partitioning/gamma that works, beyond 0.9/0.1, 0.9/0.9.

parent 2d39e5e2
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
FloryHugg_DiffUnbleached.ipynb
+ 20
− 3
View file @ 44eadcae
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import dolfin as df import dolfin as df
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import mshr as ms import mshr as ms
import numpy as np import numpy as np
import time import time
df.set_log_level(40) df.set_log_level(40)
# domain = ms.Sphere(df.Point(0, 0, 0), 1.0) # domain = ms.Sphere(df.Point(0, 0, 0), 1.0)
# mesh = ms.generate_mesh(domain, 50) # mesh = ms.generate_mesh(domain, 50)
mesh = df.UnitIntervalMesh(10000) mesh = df.UnitIntervalMesh(10000)
F = df.FunctionSpace(mesh, 'CG', 1) F = df.FunctionSpace(mesh, 'CG', 1)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def calc_sim(c0, c_tot, Ga0, dt, n_t, sym): def calc_sim(c0, c_tot, Ga0, dt, n_t, sym):
tc = df.TestFunction(F) tc = df.TestFunction(F)
c = df.Function(F) c = df.Function(F)
X = df.SpatialCoordinate(mesh) X = df.SpatialCoordinate(mesh)
# X.interpolate(df.Expression('x[0]', degree=1)) # X.interpolate(df.Expression('x[0]', degree=1))
if sym == 0: if sym == 0:
# Weak form 1D: # 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), df.grad((1-c_tot)*Ga0*tc)) -
# df.inner(df.grad(c_tot), df.grad((1-c_tot)*Ga0/c_tot*c*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), df.grad((1-c_tot)*Ga0))+
# tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c*Ga0))) * df.dx # tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c*Ga0))) * df.dx
# # Weak form 1D short: # # Weak form 1D short:
form = (df.inner((c-c0)/dt, tc) + form = (df.inner((c-c0)/dt, tc) +
df.inner((1-c_tot)*Ga0*(df.grad(c)-c/c_tot*df.grad(c_tot)), df.inner((1-c_tot)*Ga0*(df.grad(c)-c/c_tot*df.grad(c_tot)),
df.grad(tc))) * df.dx df.grad(tc))) * df.dx
elif sym == 2: elif sym == 2:
# Weak form radial symmetry: # Weak form radial symmetry:
# form = (df.inner((c-c0)/dt, tc*X[0]*X[0]) + # 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), 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]))- # 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), 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]))- # 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*2*X[0]*c.dx(0)*tc+
# (1-c_tot)*Ga0/c_tot*c*2*X[0]*c_tot.dx(0)*tc) * df.dx # (1-c_tot)*Ga0/c_tot*c*2*X[0]*c_tot.dx(0)*tc) * df.dx
# # Weak form radial symmetry: # # Weak form radial symmetry:
form = (df.inner((c-c0)/dt, tc*X[0]**2) + form = (df.inner((c-c0)/dt, tc*X[0]**2) +
df.inner((1-c_tot)*Ga0*(df.grad(c)-c/c_tot*df.grad(c_tot)), df.inner((1-c_tot)*Ga0*(df.grad(c)-c/c_tot*df.grad(c_tot)),
df.grad(tc)*X[0]**2)) * df.dx df.grad(tc)*X[0]**2)) * df.dx
t = 0 t = 0
# Solve in time # Solve in time
ti = time.time() ti = time.time()
for i in range(n_t): for i in range(n_t):
df.solve(form == 0, c) df.solve(form == 0, c)
df.assign(c0, c) df.assign(c0, c)
t += dt t += dt
print(time.time() - ti) print(time.time() - ti)
return c0 return c0
def p_tot(p_i, p_o): def p_tot(p_i, p_o):
return str(p_i-p_o)+'*(-0.5*tanh(35000*(x[0]-0.1))+0.5)+'+str(p_o) return str(p_i-p_o)+'*(-0.5*tanh(35000*(x[0]-0.1))+0.5)+'+str(p_o)
# return '(x[0]<0.1 ? '+str(p_i)+':'+ str(p_o)+')' # return '(x[0]<0.1 ? '+str(p_i)+':'+ str(p_o)+')'
def create_func(f_space, expr_str, deg): def create_func(f_space, expr_str, deg):
f = df.Function(f_space) f = df.Function(f_space)
f.interpolate(df.Expression(expr_str, degree=deg)) f.interpolate(df.Expression(expr_str, degree=deg))
return f return f
def eval_func(func, x): def eval_func(func, x):
return np.array([func([x]) for x in x]) return np.array([func([x]) for x in x])
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
c_tot_1 = create_func(F, p_tot(0.9, 0.9), 1) 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) 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) Ga0_1 = create_func(F, p_tot(1, 1/9), 1)
c_tot_2 = create_func(F, p_tot(0.9, 0.45), 1)
c0_2 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(0.9, 0.45), 1)
Ga0_2 = create_func(F,p_tot(1, 0.08), 1)
c_tot_9 = create_func(F, p_tot(0.9, 0.1), 1) 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) 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) Ga0_9 = create_func(F,p_tot(1, 1), 1)
c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1, 0.001, 10, 2) c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1, 0.001, 50, 0)
c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9, 0.001, 10, 2) c0_2 = calc_sim(c0_2, c_tot_2, Ga0_2, 0.001, 50, 0)
c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9, 0.001, 50, 0)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
x = np.linspace(0, 1, 10000) x = np.linspace(0, 1, 10000)
plt.plot(x, eval_func(c0_1, x)) plt.plot(x, eval_func(c0_1, x))
plt.plot(x, eval_func(c0_2, x))
plt.plot(x, eval_func(c0_9, x)) plt.plot(x, eval_func(c0_9, x))
plt.xlim(0, 0.2) plt.xlim(0.0, 0.125)
```
%% Cell type:code id: tags:
``` python
# Check partitioning
cp = eval_func(c0_9, x)
print('Partitioning: ' + str(np.max(cp[1:1050])/np.min(cp[999:1050])))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Set B=1 # Set B=1
pi = 0.8 pi = 0.8
po = 0.5 - np.sqrt(0.25 + (pi**2-pi)) po = 0.5 - np.sqrt(0.25 + (pi**2-pi))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
p1_i = 0.9; p1_o = 0.1 p1_i = 0.9; p1_o = 0.1
p2_i = 0.8; p2_o = 0.2 p2_i = 0.8; p2_o = 0.2
p3_i = 0.7; p3_o = 0.3 p3_i = 0.7; p3_o = 0.3
p4_i = 0.9; p4_o = 0.9 p4_i = 0.9; p4_o = 0.9
ct_1 = create_func(F, p_tot(p1_i, p1_o), 1) 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_2 = create_func(F, p_tot(p2_i, p2_o), 1)
ct_3 = create_func(F, p_tot(p3_i, p3_o), 1) ct_3 = create_func(F, p_tot(p3_i, p3_o), 1)
ct_4 = create_func(F, p_tot(p4_i, p4_o), 1) ct_4 = create_func(F, p_tot(p4_i, p4_o), 1)
g_1= create_func(F, '1', 1) g_1= create_func(F, '1', 1)
g_2= create_func(F, '0.5', 1) g_2= create_func(F, '0.5', 1)
g_3= create_func(F, '0.33333333333333333333', 1) g_3= create_func(F, '0.33333333333333333333', 1)
g_4= create_func(F, '1', 1) g_4= create_func(F, '1', 1)
c0_1 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p1_i, p1_o), 1) c0_1 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p1_i, p1_o), 1)
c0_2 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p2_i, p2_o), 1) c0_2 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p2_i, p2_o), 1)
c0_3 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p3_i, p3_o), 1) c0_3 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p3_i, p3_o), 1)
c0_4 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p4_i, p4_o), 1) c0_4 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p4_i, p4_o), 1)
for i in range(10): for i in range(10):
c0_1 = calc_sim(c0_1, ct_1, g_1, 0.001, 2, 2) c0_1 = calc_sim(c0_1, ct_1, g_1, 0.001, 2, 2)
# c0_2 = calc_sim(c0_2, ct_2, g_2, 0.01, 2, 2) # c0_2 = calc_sim(c0_2, ct_2, g_2, 0.01, 2, 2)
# c0_3 = calc_sim(c0_3, ct_3, g_3, 0.01, 2, 2) # c0_3 = calc_sim(c0_3, ct_3, g_3, 0.01, 2, 2)
c0_4 = calc_sim(c0_4, ct_4, g_4, 0.001, 2, 2) c0_4 = calc_sim(c0_4, ct_4, g_4, 0.001, 2, 2)
plt.plot(x, eval_func(c0_1, x), 'r')#/eval_func(c0_1, [0.099]) plt.plot(x, eval_func(c0_1, x), 'r')#/eval_func(c0_1, [0.099])
# plt.plot(x, eval_func(c0_2, x), 'g') # plt.plot(x, eval_func(c0_2, x), 'g')
# plt.plot(x, eval_func(c0_3, x)/eval_func(c0_3, [0.099]), 'b') # plt.plot(x, eval_func(c0_3, x)/eval_func(c0_3, [0.099]), 'b')
plt.plot(x, eval_func(c0_4, x), 'k') plt.plot(x, eval_func(c0_4, x), 'k')
plt.xlim(0.0, 0.625) plt.xlim(0.0, 0.625)
plt.ylim(0, 1.1) plt.ylim(0, 1.1)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
plt.plot(x, eval_func(c0_1, x)/0.9) plt.plot(x, eval_func(c0_1, x)/0.9)
plt.plot(x, eval_func(c0_2, x)/0.8) plt.plot(x, eval_func(c0_2, x)/0.8)
plt.plot(x, eval_func(c0_3, x)/0.7) plt.plot(x, eval_func(c0_3, x)/0.7)
plt.xlim(0.0, 0.125) plt.xlim(0.0, 0.125)
plt.ylim(0, 1.1) plt.ylim(0, 1.1)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Radial diffusion equation ## Radial diffusion equation
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
mesh = df.UnitIntervalMesh(1000) mesh = df.UnitIntervalMesh(1000)
dt = 0.001 dt = 0.001
F = df.FunctionSpace(mesh, 'CG', 1) F = df.FunctionSpace(mesh, 'CG', 1)
c0 = df.Function(F) c0 = df.Function(F)
c0.interpolate(df.Expression('x[0]<0.5 && x[0]>0.2 ? 1:0', degree=1)) c0.interpolate(df.Expression('x[0]<0.5 && x[0]>0.2 ? 1:0', degree=1))
q = df.TestFunction(F) q = df.TestFunction(F)
c = df.Function(F) c = df.Function(F)
X = df.SpatialCoordinate(mesh) X = df.SpatialCoordinate(mesh)
g = df.Expression('.00', degree=1) g = df.Expression('.00', degree=1)
u_D = df.Expression('1', degree=1) u_D = df.Expression('1', degree=1)
def boundary(x, on_boundary): def boundary(x, on_boundary):
return on_boundary return on_boundary
bc = df.DirichletBC(F, u_D, boundary) bc = df.DirichletBC(F, u_D, boundary)
# Weak form spherical symmetry # Weak form spherical symmetry
form = (df.inner((c-c0)/dt, q*X[0]*X[0]) + form = (df.inner((c-c0)/dt, q*X[0]*X[0]) +
df.inner(df.grad(c), df.grad(X[0]*X[0]*q))- df.inner(df.grad(c), df.grad(X[0]*X[0]*q))-
c.dx(0)*2*X[0]*q) * df.dx c.dx(0)*2*X[0]*q) * df.dx
# Weak form 1D # Weak form 1D
# form = (df.inner((c-c0)/dt, q) + # form = (df.inner((c-c0)/dt, q) +
# df.inner(df.grad(c), df.grad(q))) * df.dx # df.inner(df.grad(c), df.grad(q))) * df.dx
t = 0 t = 0
# Solve in time # Solve in time
for i in range(60): for i in range(60):
print(np.sum([x*x*c0([x]) for x in np.linspace(0, 1, 1000)])) print(np.sum([x*x*c0([x]) for x in np.linspace(0, 1, 1000)]))
df.solve(form == 0, c) df.solve(form == 0, c)
df.assign(c0, c) df.assign(c0, c)
t += dt t += dt
plt.plot(np.linspace(0, 1, 1000), [c0([x]) for x in np.linspace(0, 1, 1000)]) plt.plot(np.linspace(0, 1, 1000), [c0([x]) for x in np.linspace(0, 1, 1000)])
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` 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