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

Profiles amended. Bleaching close to zero instead of middle. Phi_tot accordingly.

parent ba597e58
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
FloryHugg_DiffUnbleached.ipynb
+ 18
− 15
View file @ 02a451b9
%% 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)
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
t = 0
# Solve in time
ti = time.time()
for i in range(6):
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, high partitioning
# 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.4*tanh(-350000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001))+0.5', degree=1))
c0_1.interpolate(df.Expression(('sqrt((x[0]-0.5)*(x[0]-0.5))<0.001 ? 0 :'
'0.4*tanh(-350000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001)) +0.5'),
c_tot_1.interpolate(df.Expression('0*tanh(35000*(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'),
degree=1))
Ga0_1.interpolate(df.Expression('0*(tanh(-350000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001))+1)+1', degree=1))
Ga0_1.interpolate(df.Expression('4.*(tanh(35000*(x[0]-0.01))+1)+1', degree=1))
# 1D, no partitioning
# 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*tanh(-350000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001))+0.9', degree=1))
c0_9.interpolate(df.Expression(('sqrt((x[0]-0.5)*(x[0]-0.5))<0.001 ? 0 :'
'0*tanh(-350000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001)) +0.9'),
c_tot_9.interpolate(df.Expression('0.4*tanh(-35000*(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'),
degree=1))
Ga0_9.interpolate(df.Expression('4.*(tanh(350000*(sqrt((x[0]-0.5)*(x[0]-0.5))-0.001))+1)+1', degree=1))
Ga0_9.interpolate(df.Expression('0*(tanh(-35000*(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), [c0_9([x]) for x in np.linspace(0, 1, 10000)])
plt.xlim(0.495, 0.505)
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.495, 0.505)
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)
```
%% Cell type:code id: tags:
``` python
plt.plot(np.linspace(0, 1, 2000), [Ga0([x]) for x in np.linspace(0, 1, 2000)])
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([x]) for x in np.linspace(0, 1, 1000)])
plt.plot(np.linspace(0, 1, 1000), [c_tot_9([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