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

Swapping Gamma and D works in principle for 1D.

parent 29995129
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
FloryHugg_DiffUnbleached.ipynb
+ 44
− 12
View file @ 21570d0d
%% 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(1000)
dt = 0.01
mesh = df.UnitIntervalMesh(10000)
dt = 0.0000002
F = df.FunctionSpace(mesh, 'CG', 1)
c0 = df.Function(F)
c_tot = df.Function(F)
c = df.Function(F)
tc = df.TestFunction(F)
Ga0 = df.Function(F)
# 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:
c_tot.interpolate(df.Expression('0.4*tanh(-350*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.2))+0.5', degree=1))
c0.interpolate(df.Expression(('sqrt((x[0]-0.5)*(x[0]-0.5))<0.2 ? 0 :'
'0.4*tanh(-350*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.2)) +0.5'),
# 1D, high partitioning
# c_tot.interpolate(df.Expression('0.4*tanh(-35000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001))+0.5', degree=1))
# c0.interpolate(df.Expression(('sqrt((x[0]-0.5)*(x[0]-0.5))<0.001 ? 0 :'
# '0.4*tanh(-35000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001)) +0.5'),
# degree=1))
# Ga0.interpolate(df.Expression('0*(tanh(-35000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001))+1)+1', degree=1))
# 1D, no partitioning
c_tot.interpolate(df.Expression('0*tanh(-35000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001))+0.9', degree=1))
c0.interpolate(df.Expression(('sqrt((x[0]-0.5)*(x[0]-0.5))<0.001 ? 0 :'
'0*tanh(-35000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001)) +0.9'),
degree=1))
Ga0.interpolate(df.Expression('4.*(tanh(35000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001))+1)+1', degree=1))
# tanh from left to right initial conditions
# c_tot.interpolate(df.Expression('0.4*tanh(5*(x[0]-0.5)) - 1', degree=1))
# Step function left to right initial condition
# c0.interpolate(df.Expression('(x[0]<0.5) ? 0 : 1.0', degree=1))
# Somewhat realistic half FRAP
# Weak form
form = ((df.inner((c-c0)/dt, tc) + (1-c_tot)* df.inner(df.grad(c), df.grad(tc))) -
(1-c_tot)/c_tot* c*df.inner(df.grad(c_tot), df.grad(tc))) * df.dx
form = ((df.inner((c-c0)/dt, tc) + 1/Ga0*(1-c_tot)* df.inner(df.grad(c), df.grad(tc))) -
1/Ga0*(1-c_tot)/c_tot* c*df.inner(df.grad(c_tot), df.grad(tc))) * df.dx
# Open file for writing
cFile = df.XDMFFile('conc.xdmf')
t = 0
cFile.write(c0, t)
# Solve in time
ti = time.time()
for i in range(8):
for i in range(10):
print(time.time() - ti)
df.solve(form == 0, c)
df.assign(c0, c)
t += dt
cFile.write(c0, t)
cFile.close()
print(time.time() - ti)
```
%% Cell type:code id: tags:
``` python
# 1D:
plt.plot(np.linspace(0, 1, 1000), [c0([x]) for x in np.linspace(0, 1, 1000)])
plt.plot(np.linspace(0, 1, 10000), [c0([x]) for x in np.linspace(0, 1, 10000)])
plt.xlim(0.495, 0.505)
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
c0([0, 0, 0])
# 1D:
plt.plot(np.linspace(0, 1, 10000), [c0([x]) for x in np.linspace(0, 1, 10000)])
plt.xlim(0.495, 0.505)
plt.ylim(0., 0.3)
```
%% Cell type:code id: tags:
``` python
plt.plot(np.linspace(0, 1, 2000), [Ga0([x]) for x in np.linspace(0, 1, 2000)])
```
%% Cell type:code id: tags:
``` python
plt.plot(np.linspace(0, 1, 1000), [c_tot([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