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

Simplifying weak form. Not sure whether entirely correct. 

Where does r^2 go? Inside or outside of inner product? Doesn't seem to matter?
parent 8c0d59c0
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
FloryHugg_DiffUnbleached.ipynb
+ 9
− 6
View file @ d8a3ef12
%% 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)),
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*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(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 = ((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