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

Radial Droplet FRAP qualitatively works. 

But only up to a factor between different partitionings/diffusion coefficients outside.
Problem: Maybe need different concentration scaling for radial problem?
parent 2efcf4c0
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
FloryHugg_DiffUnbleached.ipynb
+ 27
− 30
View file @ 4665e64e
%% 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(100000)
mesh = df.UnitIntervalMesh(10000)
dt = 0.000001
F = df.FunctionSpace(mesh, 'CG', 1)
```
%% Cell type:code id: tags:
``` python
def calc_sim(c0, c_tot, Ga0):
tc = df.TestFunction(F)
c = df.Function(F)
# Weak form
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
X = df.SpatialCoordinate(mesh)
# X.interpolate(df.Expression('x[0]', degree=1))
# 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 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
t = 0
# Solve in time
ti = time.time()
for i in range(6):
print(time.time() - ti)
for i in range(60):
# print(time.time() - ti)
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(35000*(x[0]-0.01))+0.9', degree=1))
c_tot_1.interpolate(df.Expression('0*tanh(350000*(x[0]-0.01))+0.9', degree=1))
c0_1.interpolate(df.Expression(('x[0]<0.01 ? 0 :'
'0*tanh(35000*(x[0]-0.01))+0.9'),
'0*tanh(350000*(x[0]-0.01))+0.9'),
degree=1))
Ga0_1.interpolate(df.Expression('4.*(tanh(35000*(x[0]-0.01))+1)+1', degree=1))
Ga0_1.interpolate(df.Expression('4.*(tanh(350000*(x[0]-0.01))+1)+1', 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(-35000*(x[0]-0.01))+0.5', degree=1))
c_tot_9.interpolate(df.Expression('0.4*tanh(-350000*(x[0]-0.01))+0.5', degree=1))
c0_9.interpolate(df.Expression(('x[0]<0.01 ? 0 :'
'0.4*tanh(-35000*(x[0]-0.01))+0.5'),
'0.4*tanh(-350000*(x[0]-0.01))+0.5'),
degree=1))
Ga0_9.interpolate(df.Expression('0*(tanh(-35000*(x[0]-0.01))+1)+1', degree=1))
Ga0_9.interpolate(df.Expression('0*(tanh(-350000*(x[0]-0.01))+1)+1', degree=1))
c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1)
c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9)
```
%% Cell type:code id: tags:
``` python
# 1D:
plt.plot(np.linspace(0, 1, 10000), [c0_1([x]) for x in np.linspace(0, 1, 10000)])
plt.plot(np.linspace(0, 1, 10000), [1.5*c0_1([x]) for x in np.linspace(0, 1, 10000)])
plt.plot(np.linspace(0, 1, 10000), [c0_9([x]) for x in np.linspace(0, 1, 10000)])
plt.xlim(0, 0.1)
# 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)])
```
%% Cell type:code id: tags:
``` python
# 1D:
plt.plot(np.linspace(0, 1, 10000), [c0([x]) for x in np.linspace(0, 1, 10000)])
# plt.xlim(0, 0.1)
# plt.ylim(0., 0.3)
```
%% 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)
```
%% Cell type:code id: tags:
``` python
plt.plot(np.linspace(0, 1, 1000), [c_tot_9([x]) for x in np.linspace(0, 1, 1000)])
```
%% Cell type:markdown id: tags:
## Radial diffusion equation
%% Cell type:code id: tags:
``` python
import dolfin as df
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