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

Merge remote-tracking branch 'origin/master'

parents 93f115a6 53d9202c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Meshes/Meshes_Sphere.geo 0 → 100644
+ 30
− 0
View file @ 79b2d4bf
//+
SetFactory("OpenCASCADE");
Sphere(1) = {0, 0, 0, 1, -Pi/2, Pi/2, 2*Pi};
//+
Point(3) = {0, 0, 0, 0.015};
lc = .15;
Physical Surface("surfacedomain") = {1};
//+
Physical Volume("vilumedomain") = {1};
//+
Field[1] = Distance;
//+
Field[1].NodesList = {3};
Field[1].NNodesByEdge = 20;
Field[2] = Threshold;
Field[2].IField = 1;
Field[2].LcMin = lc / 10;
Field[2].LcMax = lc;
Field[2].DistMin = 0.3;
Field[2].DistMax = 0.5;
Field[3] = Min;
Field[3].FieldsList = {2};
Background Field = 3;
Solve_Phi_u.ipynb
+ 47
− 14
View file @ 79b2d4bf
%% Cell type:code id: tags:
``` python
from __future__ import print_function
from fenics import *
from mshr import *
from dolfin import *
import matplotlib.pyplot as plt
import time
from math import *
set_log_level(40)
def solver_disk(Radius_Drop, Radius_Disk):
dt = 0.5
mesh = generate_mesh(Circle(Point(0,0),Radius_Disk),120)
V = FunctionSpace(mesh, 'CG', 1)
Phi_u = Function(V)
v = TestFunction(V)
Phi_u0 = Function(V)
Phi_tot = Function(V)
Phi_tot = interpolate(Expression('0.5-0.4*tanh(sqrt(x[0]*x[0]+x[1]*x[1]) - 5) ', degree=2), V)
Phi_u0 = interpolate(Expression('0.0505 - (0.0495)*tanh(-sqrt(x[0]*x[0]+x[1]*x[1]) + 5) ', degree=2), V)
def boundary(x, on_boundary):
tol = 1E-14
return on_boundary
u_D = Constant(0.1)
bcC = DirichletBC(V, u_D, boundary)
Form = (inner((Phi_u-Phi_u0)/dt, v) - inner(((1-Phi_tot)/Phi_tot)*Phi_u*grad(Phi_tot) - (1-Phi_tot)*grad(Phi_u), grad(v))) *dx
t = 0
cFile = XDMFFile('Phi_u_Disk.xmdf')
cFile.write(Phi_u0, t)
ti = time.time()
for i in range(100):
solve(Form == 0, Phi_u, bcC)
assign(Phi_u0, Phi_u)
t += dt
cFile.write(Phi_u0, t)
print(time.time() - ti)
cFile.close()
return 0
def solver_sphere(Radius_Drop, Radius_Sphere):
dt = 0.15
mesh = generate_mesh(Sphere(Point(0,0,0),Radius_Sphere),30)
# mesh = generate_mesh(Sphere(Point(0,0,0),Radius_Sphere),30)
mesh = Mesh('Meshes_Sphere.xml')
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, P1)
V = FunctionSpace(mesh, 'CG', 1)
# V = FunctionSpace(mesh, 'CG', 1)
Phi_u = Function(V)
v = TestFunction(V)
Phi_u0 = Function(V)
Phi_tot = Function(V)
tol = 1E-14
Phi_tot = interpolate(Expression('0.5-0.4*tanh(sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]) - Radius_Drop) ', Radius_Drop = Radius_Drop, degree=2), V)
Phi_u0 = interpolate(Expression('0.0505 - (0.0495)*tanh(-sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]) + Radius_Drop) ',Radius_Drop = Radius_Drop , degree=2), V)
u_D = Constant(0.01)
def boundary(x, on_boundary):
tol = 1E-14
return on_boundary
bcC = DirichletBC(V, u_D, boundary)
bc = DirichletBC(V, u_D, boundary)
Form = (inner((Phi_u-Phi_u0)/dt, v) - inner(((1-Phi_tot)/Phi_tot)*Phi_u*grad(Phi_tot) - (1-Phi_tot)*grad(Phi_u), grad(v))) *dx
t = 0
cFile = XDMFFile('Phi_u_Sphere.xmdf')
cFile.write(Phi_u0, t)
ti = time.time()
for i in range(100):
solve(Form == 0, Phi_u, bcC)
for i in range(10):
#solve(Form == 0, Phi_u, bcC)
solve(Form == 0, Phi_u, bc,solver_parameters={'newton_solver': {'linear_solver': 'gmres'}})
assign(Phi_u0, Phi_u)
t += dt
cFile.write(Phi_u0, t)
print(time.time() - ti)
cFile.close()
return 0
def solver_half_frap_disk(Radius_Drop, Radius_Disk):
dt = 0.01
mesh = generate_mesh(Circle(Point(0,0),Radius_Disk),100)
V = FunctionSpace(mesh, 'CG', 1)
Phi_u = Function(V)
v = TestFunction(V)
Phi_u0 = Function(V)
Phi_tot = Function(V)
tol = 1E-14
Phi_u0_int = 0.01
Phi_u0_ext = 0.05
Phi_u0_int_notbleach = 0.35
Phi_u0 = interpolate(Expression('sqrt(x[1]*x[1]+x[0]*x[0]) <= Radius_Drop && atan2(x[1],x[0]) < 0 ? Phi_u0_int : sqrt(x[1]*x[1]+x[0]*x[0]) <= Radius_Drop && atan2(x[1],x[0]) >= 0 ? Phi_u0_int_notbleach : Phi_u0_ext' , degree=2,tol=tol,Phi_u0_int = Phi_u0_int,Phi_u0_ext=Phi_u0_ext, Phi_u0_int_notbleach = Phi_u0_int_notbleach, Radius_Drop=Radius_Drop),V)
Phi_tot_int = 0.35
Phi_tot_ext = 0.05
Phi_tot = interpolate(Expression('sqrt(x[1]*x[1]+x[0]*x[0]) <= Radius_Drop ? Phi_tot_int : Phi_tot_ext', degree=2,tol=tol,Phi_tot_int=Phi_tot_int,Phi_tot_ext=Phi_tot_ext,Radius_Drop=Radius_Drop),V)
u_D = Constant(0.05)
def boundary(x, on_boundary):
tol = 1E-14
return on_boundary
bc = DirichletBC(V, u_D, boundary)
Form = (inner((Phi_u-Phi_u0)/dt, v) - inner(((1-Phi_tot)/Phi_tot)*Phi_u*grad(Phi_tot) - (1-Phi_tot)*grad(Phi_u), grad(v))) *dx
t = 0
cFile = XDMFFile('Phi_u_half_frap_disk.xmdf')
cFile.write(Phi_u0, t)
ti = time.time()
for i in range(50):
solve(Form == 0, Phi_u, bc)
assign(Phi_u0, Phi_u)
t += dt
cFile.write(Phi_u0, t)
print(time.time() - ti)
cFile.close()
def solver_half_frap_sphere(Radius_Drop, Radius_Sphere):
dt = 0.01
mesh = generate_mesh(Sphere(Point(0,0,0),Radius_Sphere),50)
V = FunctionSpace(mesh, 'CG', 1)
#mesh = generate_mesh(Sphere(Point(0,0,0),Radius_Sphere),50)
mesh = Mesh('Meshes_Sphere.xml')
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, P1)
Phi_u = Function(V)
v = TestFunction(V)
Phi_u0 = Function(V)
Phi_tot = Function(V)
tol = 1E-14
Phi_u0_int = 0.01
Phi_u0_ext = 0.05
Phi_u0_int_notbleach = 0.35
Phi_u0 = interpolate(Expression('sqrt(x[1]*x[1]+x[0]*x[0]+x[2]*x[2]) <= Radius_Drop && atan2(x[1],x[0]) < 0 ? Phi_u0_int : sqrt(x[1]*x[1]+x[0]*x[0]+x[2]*x[2]) <= Radius_Drop && atan2(x[1],x[0]) >= 0 ? Phi_u0_int_notbleach : Phi_u0_ext' , degree=2,tol=tol,Phi_u0_int = Phi_u0_int,Phi_u0_ext=Phi_u0_ext, Phi_u0_int_notbleach = Phi_u0_int_notbleach, Radius_Drop=Radius_Drop),V)
Phi_tot_int = 0.35
Phi_tot_ext = 0.05
Phi_tot = interpolate(Expression('sqrt(x[1]*x[1]+x[0]*x[0]+x[2]*x[2]) <= Radius_Drop ? Phi_tot_int : Phi_tot_ext', degree=2,tol=tol,Phi_tot_int=Phi_tot_int,Phi_tot_ext=Phi_tot_ext,Radius_Drop=Radius_Drop),V)
u_D = Constant(0.05)
def boundary(x, on_boundary):
tol = 1E-14
return on_boundary
bc = DirichletBC(V, u_D, boundary)
Form = (inner((Phi_u-Phi_u0)/dt, v) - inner(((1-Phi_tot)/Phi_tot)*Phi_u*grad(Phi_tot) - (1-Phi_tot)*grad(Phi_u), grad(v))) *dx
t = 0
cFile = XDMFFile('Phi_u_half_frap_sphere.xmdf')
cFile.write(Phi_u0, t)
ti = time.time()
for i in range(50):
solve(Form == 0, Phi_u, bc)
#solve(Form == 0, Phi_u, bc)
solve(Form == 0, Phi_u, bc,solver_parameters={'newton_solver': {'linear_solver': 'gmres'}})
assign(Phi_u0, Phi_u)
t += dt
cFile.write(Phi_u0, t)
print(time.time() - ti)
cFile.close()
return 0
#solver_disk(3,10)
#solver_sphere(2,10)
#solver_sphere(0.2,1)
#solver_half_frap_disk(2,10)
#solver_half_frap_sphere(2,10)
solver_half_frap_sphere(0.2,1)
```
%% Output
1.384448766708374
2.8050124645233154
4.201659917831421
5.629311561584473
7.02629828453064
8.605925798416138
10.40172266960144
12.074671268463135
13.842711925506592
15.555825471878052
0
%% Cell type:code id: tags:
``` python
```
%% 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