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

Correcting weak form spherical diffusion equation.

parent 880c691b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
FloryHugg_DiffUnbleached.ipynb
+ 5
− 4
View file @ 8c0d59c0
%% 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 = (df.inner((c-c0)/dt, tc) +
df.inner((1-c_tot)*Ga0*(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
t = 0
# Solve in time
ti = time.time()
for i in range(n_t):
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(35000*(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])
```
%% Cell type:code id: tags:
``` python
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)
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_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)
Ga0_3 = create_func(F,p_tot(1, 0.145), 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)
Ga0_9 = create_func(F,p_tot(1, 1), 1)
c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1, 0.001, 10, 0)
c0_2 = calc_sim(c0_2, c_tot_2, Ga0_2, 0.001, 10, 0)
c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9, 0.001, 10, 0)
c0_3 = calc_sim(c0_3, c_tot_3, Ga0_3, 0.001, 10, 0)
```
%% 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_9, x))
plt.xlim(0.0, 0.125)
plt.show()
# Diffusion versus partitioning
plt.plot([Ga0_1(1)*(1-c_tot_1(1)), Ga0_2(1)*(1-c_tot_2(1)), Ga0_3(1)*(1-c_tot_3(1)), Ga0_9(1)*(1-c_tot_9(1))],
[c_tot_1(0)/c_tot_1(1), c_tot_2(0)/c_tot_2(1), c_tot_3(0)/c_tot_3(1), c_tot_9(0)/c_tot_9(1)])
plt.ylim(0, 9); plt.xlim(0, 1)
plt.plot(np.linspace(0, 1, 100), 9.5*np.sqrt(np.linspace(0, 1, 100)))
```
%% 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 = (df.inner((c-c0)/dt, q*X[0]*X[0]) +
df.inner(df.grad(c), df.grad(X[0]*X[0]*q))-
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 = (df.inner((c-c0)/dt, q) +
# df.inner(df.grad(c), df.grad(q))) * df.dx
# 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