Flow over NACA0012-Hydrofoil#
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.05
# Initialize the gmsh module
gmsh.initialize()
# Load the .geo file
gmsh.merge('mesh/naca0012_2d.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/naca0012_2d.geo'...
Info : Done reading 'mesh/naca0012_2d.geo'
Info : Meshing 1D...
Info : [ 0%] Meshing curve 1 (Nurb)
Info : [ 20%] Meshing curve 2 (Nurb)
Info : [ 30%] Meshing curve 3 (Nurb)
Info : [ 40%] Meshing curve 4 (Nurb)
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.00190617s, CPU 0.008755s)
Info : Meshing 2D...
Info : Meshing surface 11 (Plane, Frontal-Delaunay)
Info : Done meshing 2D (Wall 0.894243s, CPU 1.03443s)
Info : 15891 nodes 31594 elements
Info : Optimizing mesh (Netgen)...
Info : Done optimizing mesh (Wall 3.23728e-06s, CPU 1.2e-05s)
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['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 = 1e-2
num_steps = int(T/dt)
# Equivalent viscosity
hbar = 0.075
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/NACA/'
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/50 [00:00<?, ?it/s]Advancing in time: 100%|██████████| 50/50 [02:33<00:00, 3.05s/it]
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: 5.026e-01, Max time: 5.298e-01
Step 1b: Min time: 1.628e-03, Max time: 1.942e-03
Step 1c: Min time: 9.038e-02, Max time: 1.071e-01
Step 2 : Min time: 2.246e+00, Max time: 2.332e+00
Step 3a: Min time: 9.520e-02, Max time: 1.049e-01
Step 3a: Min time: 9.520e-02, Max time: 1.049e-01