Tube Bundle#
This notebook implements the ISF algorithm to simulate inviscid fluids, governed by the incompressible Euler equations. The mesh is generated using the gmsh library for Python and the ISF differential problems are solved using FEniCSx, a Finite Element library.
[1]:
import numpy as np
from dolfinx.io.gmshio import model_to_mesh
from mpi4py import MPI
from dolfinx.io import XDMFFile
from dolfinx import mesh
import gmsh
mesh_comm = MPI.COMM_WORLD
model_rank = 0
mesh_factor = 0.035
# Initialize the gmsh module
gmsh.initialize()
# Load the .geo file
gmsh.merge('mesh/tube_bundle2D.geo')
gmsh.model.geo.synchronize()
# Set algorithm (adaptive = 1, Frontal-Delaunay = 6)
gmsh.option.setNumber("Mesh.Algorithm", 6)
gmsh.option.setNumber("Mesh.MeshSizeFactor", mesh_factor)
gdim = 2
# Linear Finite Element
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.optimize("Netgen")
Info : Reading 'mesh/tube_bundle2D.geo'...
Info : Done reading 'mesh/tube_bundle2D.geo'
Info : Meshing 1D...
Info : [ 0%] Meshing curve 1 (Line)
Info : [ 10%] Meshing curve 2 (Line)
Info : [ 20%] Meshing curve 3 (Line)
Info : [ 20%] Meshing curve 4 (Line)
Info : [ 30%] Meshing curve 5 (Line)
Info : [ 30%] Meshing curve 6 (Line)
Info : [ 40%] Meshing curve 7 (Line)
Info : [ 40%] Meshing curve 8 (Line)
Info : [ 50%] Meshing curve 9 (Line)
Info : [ 50%] Meshing curve 10 (Line)
Info : [ 60%] Meshing curve 11 (Circle)
Info : [ 70%] Meshing curve 12 (Circle)
Info : [ 70%] Meshing curve 13 (Circle)
Info : [ 80%] Meshing curve 14 (Circle)
Info : [ 80%] Meshing curve 15 (Circle)
Info : [ 90%] Meshing curve 16 (Circle)
Info : [ 90%] Meshing curve 17 (Circle)
Info : [100%] Meshing curve 18 (Circle)
Info : Done meshing 1D (Wall 0.026944s, CPU 0.045533s)
Info : Meshing 2D...
Info : Meshing surface 1 (Plane, Frontal-Delaunay)
Info : Done meshing 2D (Wall 0.429316s, CPU 0.497053s)
Info : 4154 nodes 8320 elements
Info : Optimizing mesh (Netgen)...
Info : Done optimizing mesh (Wall 2.79024e-06s, CPU 4e-06s)
Loading the mesh
[2]:
# Import into dolfinx
domain, _, ft = model_to_mesh(gmsh.model, MPI.COMM_WORLD, model_rank, gdim = gdim )
ft.name = "Facet markers"
boundary_markers = dict()
boundary_markers['inlet'] = 10
boundary_markers['adiab_walls'] = 20
boundary_markers['coldest_walls'] = 30
boundary_markers['hottest_walls'] = 40
boundary_markers['hot_centre'] = 50
boundary_markers['cold_centre'] = 60
boundary_markers['outlet'] = 70
domain_marker = 100
# Finalize the gmsh module
gmsh.finalize()
The ISF algorithm is based on the analogy between hydrodynamics and quantum mechinacs: in fact, it can be shown thta the Euler equations
\[\partial_t\mathbf{u}=-\left(\mathbf{u}\cdot \nabla\right)\mathbf{u} -\nabla p \qquad \text{ with } \nabla \cdot \mathbf{u} = 0\]
are equivalent to a 2-component wave function \(\Psi=[\psi_1,\psi_2]^T\in\mathbb{C}^2\), solution to Schrodinger equation
\[i\hbar \partial_t \Psi =- 0.5 \cdot \hbar^2 \Delta \Psi + p \Psi\]
and constrained to
\[\mathfrak{Re}\left\{\left(\Psi^T\right)^*i\,\Delta\Psi\right\} = 0\]
The solution algorithm is divided into three main steps, at each time iteration:
Advection: through the free-particle Schrodinger
Pressure projection: Poisson problem to have a divergence free velocity
Phase shift
[3]:
import tqdm
import time
import os
# Mesh generation
from dolfinx.fem import form, assemble_scalar
import ufl
from ufl import dot, FacetNormal
import matplotlib.pyplot as plt
from matplotlib import rcParams
from matplotlib import cm
plt.rcParams.update({
"text.usetex": True,
"font.family": "serif"
})
rcParams['text.latex.preamble'] = r'\usepackage{amssymb} \usepackage{amsmath} \usepackage{amsthm} \usepackage{mathtools}'
The steps of the algorithm are imported from the src folder.
[4]:
import sys
sys.path.append('../../.')
from src.prediction import schrodinger
from src.ISFvelocity import computeVelocity
from src.projection import poisson, phaseShift
Time Loop#
An inlet velocity \(\mathbf{u}_{in} = [1, 0]^T\) is used and the total simulation time is \(T\), using a time step \(\Delta t = 10^{-2}\).
[6]:
# Inlet velocity
u_in = np.array([1., 0., 0.])
# Time parameters
T = 5.
dt = 1e-2
num_steps = int(T/dt)
# Equivalent viscosity
hbar = 0.1
phys_parameters = {'uin': u_in,
'hbar': hbar}
# Degree of the polynomial for FEM - the Neumann BC at the outlet for step 2 either requires a finer mesh or degree_psi = 2
degree_psi = 1
### Assembling forms #####
print('Assembling Schrodinger predict')
step1 = schrodinger(domain, ft, boundary_marks=boundary_markers, degree_psi = degree_psi, c = np.array([.5 + 0 * 1j, .5 + 0 * 1j]))
step1.assemble(phys_parameters, dt=dt, gravity_dict={'gravity': np.array([0, -9.81, 0])})
print('Assembling Schrodinger velocity')
ISF_vel = computeVelocity(domain, degree_psi=degree_psi)
ISF_vel.assemble(hbar)
print('Assembling Poisson')
step2 = poisson(domain, ft, boundary_marks=boundary_markers)
bc_poisson_values = {'uin': u_in, 'is_dirichlet': True, 'fully_Neumann': False}
step2.assemble(bc_poisson_values, gauss = False)
print('Assembling Schrodinger phase-shift')
step3 = phaseShift(domain, degree_psi = degree_psi)
step3.assemble(hbar)
##### Store solution #####
path = './Results/'
if not os.path.exists(path):
os.system("mkdir "+path)
path = './Results/TubeBundle/'
if not os.path.exists(path):
os.system("mkdir "+path)
psi1_xdmf = XDMFFile(domain.comm, path+"psi1.xdmf", "w")
psi2_xdmf = XDMFFile(domain.comm, path+"psi2.xdmf", "w")
p_xdmf = XDMFFile(domain.comm, path+"p.xdmf", "w")
u_xdmf = XDMFFile(domain.comm, path+"U.xdmf", "w")
p_xdmf.write_mesh(domain)
psi1_xdmf.write_mesh(domain)
psi2_xdmf.write_mesh(domain)
u_xdmf.write_mesh(domain)
Assembling Schrodinger predict
Assembling Schrodinger velocity
Assembling Poisson
Assembling Schrodinger phase-shift
Let us solve the transient.
[7]:
from dolfinx import common
# Sampling rate to store the data
LL = 10
steps = ["1a", "1b", "1c", "2 ", "3a", "3a"]
t_steps = list()
progress = tqdm.tqdm(desc="Advancing in time", total=num_steps)
for nn in range(num_steps):
## Step 1a: advection
with common.Timer("~Step 1a"):
# psi1_adv, psi2_adv = step1.advect(temperature=particles.Told)
psi1_adv, psi2_adv = step1.advect(temperature=None)
## Step 1b: normalisation
with common.Timer("~Step 1b"):
PsiTilde_1 = psi1_adv.copy()
PsiTilde_2 = psi2_adv.copy()
norm = np.sqrt(np.conj(psi1_adv.x.array[:]) * psi1_adv.x.array[:] + np.conj(psi2_adv.x.array[:]) * psi2_adv.x.array[:])
PsiTilde_1.x.array[:] = psi1_adv.x.array[:] / norm.real
PsiTilde_2.x.array[:] = psi2_adv.x.array[:] / norm.real
## Step 1c: compute non-solenoidal velocity
with common.Timer("~Step 1c"):
uTilde = ISF_vel.solve(PsiTilde_1, PsiTilde_2)
## Step 2: pressure projection
with common.Timer("~Step 2 "):
pressure = step2.solve(uTilde)
## Step 3a: phase shift
with common.Timer("~Step 3a"):
Psi_1, Psi_2 = step3.solve(PsiTilde_1, PsiTilde_2, pressure)
## Step 3b: compute new velocity
with common.Timer("~Step 3b"):
u_new = ISF_vel.solve(Psi_1, Psi_2)
if np.isclose((nn+1) % LL, 0):
pressure.name = "p"
Psi_1.name = "psi_1"
Psi_2.name = "psi_2"
u_new.name = "U"
p_xdmf.write_function(pressure, (nn+1) * dt)
psi1_xdmf.write_function(Psi_1, (nn+1) * dt)
psi2_xdmf.write_function(Psi_2, (nn+1) * dt)
u_xdmf.write_function(u_new, (nn+1) * dt)
# Update Old
with Psi_1.vector.localForm() as loc_, step1.psi1_old.vector.localForm() as loc_n:
loc_.copy(loc_n)
with Psi_2.vector.localForm() as loc_, step1.psi2_old.vector.localForm() as loc_n:
loc_.copy(loc_n)
t_steps.append(list())
for step in steps:
t_steps[nn].append(mesh_comm.gather(common.timing("~Step "+step), root=0))
progress.update(1)
p_xdmf.close()
psi1_xdmf.close()
psi2_xdmf.close()
u_xdmf.close()
Advancing in time: 100%|██████████| 500/500 [03:01<00:00, 2.82it/s]
In the end, let us print the time required by each step during the simulation
[8]:
if mesh_comm.rank == 0:
print("Time-step breakdown")
step_arr = np.asarray(t_steps)
time_per_run = step_arr[:, :, 0, 1] / step_arr[:, :, 0, 0]
for jj, step in enumerate(steps):
print(f"Step "+step+f": Min time: {np.min(time_per_run[:, jj]):.3e}, Max time: {np.max(time_per_run[:, jj]):.3e}")
Time-step breakdown
Step 1a: Min time: 9.562e-02, Max time: 1.051e-01
Step 1b: Min time: 5.757e-04, Max time: 6.342e-04
Step 1c: Min time: 1.251e-02, Max time: 1.380e-02
Step 2 : Min time: 2.094e-01, Max time: 2.726e-01
Step 3a: Min time: 1.771e-02, Max time: 1.992e-02
Step 3a: Min time: 1.771e-02, Max time: 1.992e-02