Flow over cylinder (2D)#
This notebook implements the ISF algorithm to simulate inviscid fluids, governed by the incompressible Euler equations in a flow over cylinder case. The mesh is generated using the gmsh library for Python and the ISF differential problems are solved using FEniCSx, a Finite Element library.
The geometry is saved into a .geo file which is loaded and meshed using the gmsh package in Python.
[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.5
# Initialize the gmsh module
gmsh.initialize()
# Load the .geo file
gmsh.merge('mesh/cyl_dfg2D.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/cyl_dfg2D.geo'...
Info : Done reading 'mesh/cyl_dfg2D.geo'
Info : Meshing 1D...
Info : [ 0%] Meshing curve 1 (Circle)
Info : [ 20%] Meshing curve 2 (Circle)
Info : [ 30%] Meshing curve 3 (Circle)
Info : [ 40%] Meshing curve 4 (Circle)
Info : [ 50%] Meshing curve 5 (Line)
Info : [ 70%] Meshing curve 6 (Line)
Info : [ 80%] Meshing curve 7 (Line)
Info : [ 90%] Meshing curve 8 (Line)
Info : Done meshing 1D (Wall 0.00224929s, CPU 0.009821s)
Info : Meshing 2D...
Info : Meshing surface 11 (Plane, Frontal-Delaunay)
Info : Done meshing 2D (Wall 0.410715s, CPU 0.586401s)
Info : 7998 nodes 16003 elements
Info : Optimizing mesh (Netgen)...
Info : Done optimizing mesh (Wall 3.05474e-06s, CPU 5e-06s)
The mesh is loaded into dolfinx.
[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['outlet'] = 30
boundary_markers['walls'] = 20
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 os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
from matplotlib import cm
# If latex is not installed on your machine, comment these lines
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}\).
[5]:
# Inlet velocity
u_in = np.array([1., 0., 0.])
# Time parameters
T = 5.
dt = 2e-2
num_steps = int(T/dt)
# Equivalent viscosity
hbar = 0.025
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)
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/FlowCyl/'
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.
[6]:
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: 0%| | 0/250 [00:00<?, ?it/s]Advancing in time: 100%|██████████| 250/250 [03:02<00:00, 1.30it/s]
In the end, let us print the time required by each step during the simulation
[7]:
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: 2.148e-01, Max time: 2.444e-01
Step 1b: Min time: 9.480e-04, Max time: 1.192e-03
Step 1c: Min time: 3.106e-02, Max time: 4.047e-02
Step 2 : Min time: 3.861e-01, Max time: 4.876e-01
Step 3a: Min time: 4.272e-02, Max time: 5.080e-02
Step 3a: Min time: 4.272e-02, Max time: 5.080e-02