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

Python: problem seems more outside than inside droplet.

parent 8a392fbb
Branches flux_accuracy
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
FloryHugg_DiffUnbleached.ipynb
+ 28
− 24
View file @ 359f4464
%% 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)
F = df.FunctionSpace(mesh, 'CG', 1)
```
%% Cell type:code id: tags:
``` python
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) +
# 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 1D short:
form = ((c-c0)/dt*tc + (1-c_tot)*Ga0*
df.inner((df.grad(c)-c/c_tot*df.grad(c_tot)),
df.grad(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 = ((c-c0)/dt*tc*X[0]**2 +
# (1-c_tot)*Ga0*(c.dx(0)-c/c_tot*c_tot.dx(0))*
# tc.dx(0)*X[0]**2) * df.dx
form = ((c-c0)/dt*tc + (1-c_tot)*Ga0*
df.inner((df.grad(c)-c/c_tot*df.grad(c_tot)),
df.grad(tc)))*X[0]*X[0]*df.dx
t = 0
# Solve in time
ti = time.time()
for i in range(n_t):
# 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
print(time.time() - ti)
return c0
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)
# return '(x[0]<0.1 ? '+str(p_i)+':'+ str(p_o)+')'
def create_func(f_space, expr_str, deg):
f = df.Function(f_space)
f.interpolate(df.Expression(expr_str, degree=deg))
return f
def eval_func(func, x):
return np.array([func([x]) for x in x])
def eval_P(func, x):
return func(x[0])/func(x[1])
def eval_D(func_ga, func_tot, x):
return(func_ga(x)*(1-func_tot(x)))
```
%% Cell type:code id: tags:
``` python
bp = 0.1 # Boundary position at IC
sym = 2 # Symmetry of the problem
dt = 0.0001
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]<'+str(bp)+'? 0 :' + p_tot(0.9, 0.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)
c0_2 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.45), 1)
Ga0_2 = create_func(F,p_tot(1, 0.08), 1)
c_tot_3 = create_func(F, p_tot(0.9, 0.3), 1)
c0_3 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(0.9, 0.3), 1)
c0_3 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.3), 1)
Ga0_3 = create_func(F,p_tot(1, 0.145), 1)
c_tot_5 = create_func(F, p_tot(0.9, 0.18), 1)
c0_5 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(0.9, 0.18), 1)
c0_5 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.18), 1)
Ga0_5 = create_func(F,p_tot(1, 0.34), 1)
c_tot_8 = create_func(F, p_tot(0.9, 0.1125), 1)
c0_8 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(0.9, 0.1125), 1)
c0_8 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.1125), 1)
Ga0_8 = create_func(F,p_tot(1, 0.8), 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]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.1), 1)
Ga0_9 = create_func(F,p_tot(1, 1), 1)
sym = 0
c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1, 0.01, 10, sym)
c0_2 = calc_sim(c0_2, c_tot_2, Ga0_2, 0.01, 10, sym)
c0_3 = calc_sim(c0_3, c_tot_3, Ga0_3, 0.01, 10, sym)
c0_5 = calc_sim(c0_5, c_tot_5, Ga0_5, 0.01, 10, sym)
c0_8 = calc_sim(c0_8, c_tot_8, Ga0_8, 0.01, 10, sym)
c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9, 0.01, 10, sym)
c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1, dt, 10, sym)
c0_2 = calc_sim(c0_2, c_tot_2, Ga0_2, dt, 10, sym)
c0_3 = calc_sim(c0_3, c_tot_3, Ga0_3, dt, 10, sym)
c0_5 = calc_sim(c0_5, c_tot_5, Ga0_5, dt, 10, sym)
c0_8 = calc_sim(c0_8, c_tot_8, Ga0_8, dt, 10, sym)
c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9, dt, 10, sym)
```
%% Cell type:code id: tags:
``` python
x = np.linspace(0, 1, 10000)
plt.plot(x, eval_func(c0_1, x))
plt.plot(x, eval_func(c0_2, x))
plt.plot(x, eval_func(c0_3, x))
plt.plot(x, eval_func(c0_5, x))
plt.plot(x, eval_func(c0_8, x))
plt.plot(x, eval_func(c0_9, x))
# plt.xlim(0.08, 0.125125)
# plt.ylim(0., 0.325125)
plt.show()
# Diffusion versus partitioning
list_Ga = [Ga0_1, Ga0_2, Ga0_3, Ga0_5, Ga0_8, Ga0_9]
list_Tot = [c_tot_1, c_tot_2, c_tot_3, c_tot_5, c_tot_8, c_tot_9]
D = [eval_D(Ga, Tot, 1) for (Ga, Tot) in zip (list_Ga, list_Tot)]
P = [eval_P(Tot, [0, 1]) for Tot in list_Tot]
plt.plot(D, P)
plt.ylim(0, 11); plt.xlim(0, 1)
plt.ylabel('P'); plt.xlabel('Diffusion coefficient')
plt.plot(np.linspace(0, 1, 100), 9.5*(np.linspace(0, 1, 100))**(1/2))
```
%% Cell type:code id: tags:
``` python
# Plot outside, check invariance
x = np.linspace(0.1, 1, 1000)
x = np.linspace(bp, 1, 1000)
y1 = eval_func(c0_1, x)
y2 = eval_func(c0_2, x)
y3 = eval_func(c0_3, x)
y9 = eval_func(c0_9, x)
plt.plot((x-np.min(x)), y1)
plt.plot((x-np.min(x))/2, 2*y2)
plt.plot((x-np.min(x))/3, 3*y3)
plt.plot((x-np.min(x))/9, 9*y9)
# plt.xlim(0.0, 0.951252)
plt.plot((x-bp), y1)
plt.plot((x-bp)/2, 2*y2)
plt.plot((x-bp)/3, 3*y3)
plt.plot((x-bp)/9, 9*y9)
plt.xlim(0.0, 0.02)
plt.ylim(0.2, 1)
plt.show()
```
%% Cell type:code id: tags:
``` python
# Plot inside, check invariance
x = np.linspace(0, 0.1, 1000)
y1 = eval_func(c0_1, x)-np.min(eval_func(c0_1, x))
y2 = eval_func(c0_2, x)-np.min(eval_func(c0_2, x))
y3 = eval_func(c0_3, x)-np.min(eval_func(c0_3, x))
y9 = eval_func(c0_9, x)-np.min(eval_func(c0_9, x))
plt.plot((x-np.min(x)), y1)
plt.plot((x-np.min(x)), y2)
plt.plot((x-np.min(x)), y3)
plt.plot((x-np.min(x)), y9)
plt.xlim(0.0, 0.1251252)
plt.plot((x-np.min(x)), y1/eval_func(c0_1, [0.08]))
plt.plot((x-np.min(x)), y2/eval_func(c0_2, [0.08]))
plt.plot((x-np.min(x)), y3/eval_func(c0_3, [0.08]))
plt.plot((x-np.min(x)), y9/eval_func(c0_9, [0.08]))
plt.xlim(0.08, bp)
plt.show()
```
%% 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:
``` python
# Set B=1
pi = 0.8
po = 0.5 - np.sqrt(0.25 + (pi**2-pi))
```
%% Cell type:code id: tags:
``` python
p1_i = 0.9; p1_o = 0.1
p2_i = 0.8; p2_o = 0.2
p3_i = 0.7; p3_o = 0.3
p4_i = 0.9; p4_o = 0.9
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)
ct_4 = create_func(F, p_tot(p4_i, p4_o), 1)
g_1= create_func(F, '1', 1)
g_2= create_func(F, '0.5', 1)
g_3= create_func(F, '0.33333333333333333333', 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_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_4 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p4_i, p4_o), 1)
for i in range(10):
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_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)
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_3, x)/eval_func(c0_3, [0.099]), 'b')
plt.plot(x, eval_func(c0_4, x), 'k')
plt.xlim(0.0, 0.625)
plt.ylim(0, 1.1)
```
%% Cell type:code id: tags:
``` python
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_3, x)/0.7)
plt.xlim(0.0, 0.125)
plt.ylim(0, 1.1)
```
%% 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 = ((c-c0)/dt*q*X[0]*X[0] +
X[0]*X[0]*df.inner(df.grad(c), df.grad(q))-
c.dx(0)*2*X[0]*q) * df.dx
# Weak form 1D
# form = ((c-c0)/dt* q + df.inner(df.grad(c), df.grad(q))) * df.dx
# Weak form 1D with .dx(0) notation for derivative in 1st direction.
# form = ((c-c0)/dt*q + c.dx(0)*q.dx(0)) * 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