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

Correcting weak form, now all terms. Seems to work.

parent 610d21fd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
FloryHugg_DiffUnbleached.ipynb
+ 10
− 52
View file @ ba597e58
%% 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) +
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
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
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))
c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1)
c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9)
```
%% Output
9.5367431640625e-07
0.2761721611022949
0.5432901382446289
0.8246011734008789
1.096081256866455
1.3647711277008057
1.6441941261291504
1.6689300537109375e-06
0.272655725479126
0.5432279109954834
0.8257148265838623
1.1101319789886475
1.3879718780517578
1.6780588626861572
%% 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)
# 3D:
# plt.plot(np.linspace(0, 0.5, 1000), [c0([x, 0, 0]) for x in np.linspace(0, 0.5, 1000)])
```
%% Output
(0.495, 0.505)
%% 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