Nuclear Channel#
This notebook implements the ISF algorithm to simulate inviscid fluids. The mesh is generated with gmsh and loaded in Python within the FEniCSx framework.
[1]:
import numpy as np
from dolfinx.io.gmshio import model_to_mesh, read_from_msh
from mpi4py import MPI
from dolfinx.io import XDMFFile
from dolfinx import mesh
from IPython.display import clear_output
import gmsh
from petsc4py import PETSc
mesh_comm = MPI.COMM_WORLD
model_rank = 0
mesh_factor = 0.035
use_msh = True
gdim = 3
if use_msh:
domain, ct, ft = read_from_msh("mesh/nucchannel_gmsh.msh", comm=mesh_comm, rank=model_rank, gdim = gdim)
else:
# Initialize the gmsh module
gmsh.initialize()
# Load the .geo file
gmsh.merge('mesh/nucchannel_gmsh.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)
# Linear Finite Element
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.optimize("Netgen")
clear_output()
# gmsh.write('mesh/nucchannel_gmsh.vtk')
gmsh.write('mesh/nucchannel_gmsh.msh')
domain, _, ft = model_to_mesh(gmsh.model, MPI.COMM_WORLD, model_rank, gdim = gdim )
gmsh.finalize()
Info : Reading 'mesh/nucchannel_gmsh.msh'...
Info : 88 entities
Info : 40475 nodes
Info : 219289 elements
Info : Done reading 'mesh/nucchannel_gmsh.msh'
Let us define the boundary markers
[2]:
ft.name = "Facet markers"
boundary_markers = dict()
boundary_markers['inlet'] = 100
boundary_markers['outlet'] = 110
boundary_markers['lateral_walls'] = 120
boundary_markers['control_rod'] = 130
boundary_markers['top_right_nuc'] = 140
boundary_markers['top_left_nuc'] = 150
boundary_markers['bottom_left_nuc'] = 160
boundary_markers['bottom_right_nuc'] = 170
domain_marker = 1000
print(domain.geometry.x.shape)
(40475, 3)
The ISF algorithm is based on the analogy between hydrodynamics and quantum mechinacs: in fact, it can be shown thta the Euler equations \begin{equation} \partial_t\mathbf{u}=-\left(\mathbf{u}\cdot \nabla\right)\mathbf{u} -\nabla p +\mathbf{g} \qquad \text{ with } \nabla \cdot \mathbf{u} = 0 \end{equation} are equivalent to a 2-component wave function \(\Psi=[\psi_1,\psi_2]^T\in\mathbb{C}^2\), solution to Schrodinger equation \begin{equation} i\hbar \partial_t \Psi =- 0.5 \cdot \hbar^2 \Delta \Psi + p \Psi + \mathbb{G}\cdot \Psi \end{equation} and constrained to \begin{equation} \mathfrak{Re}\left\{\left(\Psi^T\right)^*i\,\Delta\Psi\right\} = 0 \end{equation}
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
Solution of other physics (e.g., temperature advection-diffusion)
[3]:
import tqdm
import os
import sys
sys.path.append('/home/paolo/Github/ISF-dev/.')
from src.prediction import schrodinger
from src.ISFvelocity import computeVelocity
from src.projection import poisson, phaseShift
from src.scalar_adv_diff import adv_diffusion
Time Loop#
An inlet velocity \(\mathbf{u}_{in}\) is used and the total simulation time is \(T\).
[4]:
# Inlet velocity
u_in = np.array([0, 0., 1.])
# Time parameters
T = 4.
dt = 1e-2
num_steps = int(T/dt)
# Equivalent viscosity
hbar = 0.05
# Assigning physical paramemeters
phys_parameters = {'uin': u_in, 'hbar': hbar}
gravity_dict = {'gravity': np.array([0, 0., -9.81]), 'beta': 0.5, 'Tref': 0}
# 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([1 + 0 * 1j, .1 + 0 * 1j]))
step1.assemble(phys_parameters, dt=dt, gravity_dict=gravity_dict)
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, 'fully_Neumann': False, 'is_dirichlet': True}
step2.assemble(bc_poisson_values, gauss = True)
print('Assembling Schrodinger phase-shift')
step3 = phaseShift(domain, degree_psi = degree_psi)
step3.assemble(hbar)
##### Store solution #####
path = 'Results/NuclearRod/'
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")
psi1_xdmf.write_mesh(domain)
psi2_xdmf.write_mesh(domain)
p_xdmf.write_mesh(domain)
u_xdmf.write_mesh(domain)
save_tilde = False
if save_tilde:
psi1_Tilde_xdmf = XDMFFile(domain.comm, path+"psi1_tilde.xdmf", "w")
psi2_Tilde_xdmf = XDMFFile(domain.comm, path+"psi2_tilde.xdmf", "w")
psi1_Tilde_xdmf.write_mesh(domain)
psi2_Tilde_xdmf.write_mesh(domain)
Assembling Schrodinger predict
Assembling Schrodinger velocity
Assembling Poisson
Assembling Schrodinger phase-shift
The advection diffusion equation reads \begin{equation} \frac{\partial T}{\partial t} + \mathbf{u}\cdot \nabla T - \alpha \Delta T =0 \end{equation}
[5]:
temperature_bound_marks = {
'dirichlet' : dict(),
'neumann' : dict(),
'robin' : None
}
# Define class for passive scalar transport
_names = {'dirichlet': ['inlet', 'control_rod'],
'neumann': ['top_right_nuc', 'top_left_nuc', 'bottom_left_nuc', 'bottom_right_nuc']}
for key in list(_names.keys()):
for name in _names[key]:
temperature_bound_marks[key][name] = boundary_markers[name]
particles = adv_diffusion(domain, ft, temperature_bound_marks)
bc_values = {'dirichlet': dict(), 'neumann': dict()}
# Define boundaries values
assigned_temperatures = [1, 0.5]
for ii in range(len(assigned_temperatures)):
bc_values['dirichlet'][_names['dirichlet'][ii]] = PETSc.ScalarType(assigned_temperatures[ii])
assigned_heat_flux = [-0.05, -0.025, 0.03, 0.015]
for ii in range(len(assigned_heat_flux)):
bc_values['neumann'][_names['neumann'][ii]] = PETSc.ScalarType(assigned_heat_flux[ii])
# Assembling the differential problem
particles.assemble(bc_values, dt = dt, alpha = 5e-4, T_0=gravity_dict['Tref'])
# Create file for storage
T_xdmf = XDMFFile(domain.comm, path+"T.xdmf", "w")
T_xdmf.write_mesh(domain)
Let us solve the transient.
[7]:
from dolfinx import common
# Sampling rate to store the data
LL = 10
steps = ["1a", "1b", "1c", "2 ", "3a", "3b", "4"]
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"):
if np.isclose(gravity_dict['beta'], 0):
psi1_adv, psi2_adv = step1.advect()
else:
psi1_adv, psi2_adv = step1.advect(temperature=particles.Told)
## 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)
## Step 4: Temperature advection-diffusion
with common.Timer("~Step 4"):
T_new = particles.solve(u_new)
if np.isclose((nn+1) % LL, 0):
pressure.name = "p"
Psi_1.name = "psi_1"
Psi_2.name = "psi_2"
u_new.name = "U"
T_new.name = "T"
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)
T_xdmf.write_function(T_new, (nn+1) * dt)
if save_tilde:
psi1_Tilde_xdmf.write_function(PsiTilde_1, (nn+1) * dt)
psi2_Tilde_xdmf.write_function(PsiTilde_2, (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()
T_xdmf.close()
if save_tilde:
psi1_Tilde_xdmf.close()
psi2_Tilde_xdmf.close()
del progress
Advancing in time: 5%|▍ | 19/400 [04:53<6:49:34, 64.50s/it]
In the end, let us print the time required by each step during the simulation
[ ]:
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: 6.168e-01, Max time: 6.399e-01
Step 1b: Min time: 2.409e-03, Max time: 2.770e-03
Step 1c: Min time: 7.203e-01, Max time: 8.996e-01
Step 2 : Min time: 8.623e-01, Max time: 9.457e-01
Step 3a: Min time: 6.385e-01, Max time: 7.601e-01
Step 3a: Min time: 6.385e-01, Max time: 7.601e-01