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

Works also for longer time courses.

parent c14c647e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
FloryHugg_DiffUnbleached.ipynb
+ 3
− 3
View file @ 3c661086
%% 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.0000002
dt = 0.000001
F = df.FunctionSpace(mesh, 'CG', 1)
```
%% 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
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'),
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))
# 1D, no 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'),
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))
```
%% 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) + 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
t = 0
# Solve in time
ti = time.time()
for i in range(10):
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
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.ylim(0, 0.3)
# 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.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