3D - Vortical Flow
Click to view the program’s source code
Compilation and Run
#! bash
np=$(getconf _NPROCESSORS_ONLN 2>/dev/null || sysctl -n hw.ncpu)
make vortical_flow_3d
mpirun -n $np ./build/bin/vortical_flow_3d
# If using IntelMPI
# mpiexec -n $np ./build/bin/vortical_flow_3d
Description
The sample program in [root_dir]/src/apps/vortical_flow_3d.f90 simulates the temporal evolution of a vortex with an axial jet. The base vortex model is the \( q \)-vortex, with a velocity field defined as:
where \( q \) is the swirling parameter, determining the relative strength of the swirl. \( r’ \) and \( \phi’ \) represent translated polar variables centered at the \( q \)-vortex’s center. In this example, \( q = 1 \). Initially, two \( q \)-vortices form a co-rotating vortex pair, with a distance of 4 between their centers \( (x,y) = (\pm 2, 0) \). The program assumes a uniform background flow \( \mathbf{U} = U_z \hat{\mathbf{e}}_z \), resulting in a total velocity field \( \mathbf{U} + \mathbf{V} \). Here, \( U_z \) is set to \( -0.5 \).
The evolution of the vortex field \( \mathbf{V}(r,\phi,z,t) \) is governed by:
\[\frac{\partial \mathbf{V}}{\partial t} = - \nabla \varphi - \left( \nabla \times \mathbf{V} \right) \times (\mathbf{U} + \mathbf{V}) + {\rm{VISC}} \cdot \nabla^2 \mathbf{V},\] \[\nabla \cdot \mathbf{V} = 0,\]where \( \varphi \) is the specific energy term. The viscosity constant \( {\rm{VISC}} (\equiv {\rm{Re}}^{-1}) \) represents the inverse Reynolds number of the flow. Using the Toroidal-Poloidal (TP) decomposition:
\[\mathbf{V} = \nabla \times \left\{\psi \hat{\mathbf{e}}_z \right\} + \nabla \times \left[ \nabla \times \left\{ \chi \hat{\mathbf{e}}_z \right\} \right],\]the equations reduce to:
\[\frac{\partial}{\partial t} \begin{bmatrix} \psi \\ \chi \end{bmatrix} = - \mathbb{P} \left[ \left( \nabla \times \mathbf{V} \right) \times (\mathbf{U} + \mathbf{V}) \right] + {\rm{VISC}} \cdot \nabla^2 \begin{bmatrix} \psi \\ \chi \end{bmatrix},\]where \( \mathbb{P} \) is the toroidal-poloidal projection operator. This operator extracts the toroidal and poloidal scalars from the vector field \( \left( \nabla \times \mathbf{V} \right) \times (\mathbf{U} + \mathbf{V}) \).
The default parameters are specified in [root_dir]/input.params:
!!! COMPUTATIONAL DOMAIN INFO !!!
# ---------- NR ----------- NP ----------- NZ ----------------------------------
200 128 128
# ------ NRCHOP ------- NPCHOP ------- NZCHOP ----------------------------------
200 65 65
# --------- ELL --------- ZLEN ------ ZLENxPI ---(IF ZLENxPI==T, ZLEN=ZLEN*PI)--
4.D0 2.D0 T
#
!!! TIME STEPPING INFO !!!
# ---------- DT ----------- TI --------- TOTT ----------- NI --------- TOTN ----
1.D-2 0.D0 5.D1 0 5000
#
!!! FIELD PROPERTY INFO !!!
# -------- VISC ----- HYPERPOW ---- HYPERVISC ----------------------------------
1.D-4 8 5.D-7
#
!!! FIELD SAVE INFO !!!
# --------------------- FLDDIR -------------------------------------------------
./output/fld
# ---- ISFLDSAV -- FLDSAVINTVL ---(IF ISFLDSAV!=T, FIELDS ARE NOT SAVED)--------
T 100
#
!!! DATA LOGGING INFO !!!
# --------------------- LOGDIR -------------------------------------------------
./output/dat
# ---- ISLOGSAV -- LOGSAVINTVL ---(IF ISLOGSAV!=T, LOGS ARE NOT GENERATED)------
F 1
!!! STABILITY CONTROL INFO (OPTIONAL) !!!
# ------- IS_SVV ------- CUTOFF ------- TARGET ------ STRENGTH ------- RELAX ---
T 0.75 2.D-2 0.12 0.25
/* ------------------------------ END OF INPUT ------------------------------ */
The final optional input block controls the bounded spectral vanishing-viscosity filter used by the vortex example. It applies a 2/3 de-alias mask to the periodic directions and increases damping only when the modal tail exceeds svv_target; a mode enters that tail when any normalized radial, azimuthal, or axial mode coordinate exceeds svv_cutoff. Set the first value to F to disable it. The filter is a tail-energy safeguard, not a controller that forces a -5/3 spectrum. With the periodic 2/3 mask, a cutoff of 0.75 primarily guards the mapped-radial tail; choose a cutoff below about 2/3 if the top retained azimuthal and axial bands must also participate. The vortex example removes the mapped far-field (x=1) mode after each step and aborts collectively if a non-finite modal value is produced.
On macOS with Homebrew GNU Fortran, install gfortran, open-mpi, and cmake, build the external libraries with FC=gfortran CC=cc BUILD_JOBS=2 ./external/CMake_build.sh, then run the same make command.
From these parameters, the simulation is set as follows:
- Collocation points in physical space: \( NR = 200 \), \( NP = 128 \), \( NZ = 128 \).
- Spectral elements in spectral space: \( NRCHOP = 200 \), \( NPCHOP = 65 \), \(NZCHOP = 65\).
- Mapping parameter: \( L = 4.0 \) (one-half of \(NR \) are in \( 0 \le r \le 4.0 \)).
- Domain Length in \( z \): \( 2 \pi \).
- Time integration parameters:
- Time step: \( dt = 0.01 \).
- Initial time: \( t_0 = 0 \) (0th step).
- Final time: \( t_f = 50 \) (5,000th step).
- Scalar field properties: \( \rm{VISC} = 10^{-4} \) (i.e., \( {\rm{Re}} = 10^4 \)) and \( \rm{HYPERVISC} = 5 \times 10^{-7} \) with \( \rm{HYPERPOW} = 8 \).
- Field information is stored at every 100 steps in
[root_dir]/output/fld/. - No log information is stored.
Note that a non-zero hyperviscous dissipation term \( - {\rm{HYPERVISC}} \cdot (-\nabla^2)^{({\rm{HYPERPOW}}/2)} \), along with the physical dissipation term \( {\rm{VISC}} \cdot \nabla^2 \), is additionally taken into account. The hyperviscosity term does not represent physical dissipation but prevents numerical instability caused by high-frequency oscillations at small scales. It should remain minimal to avoid affecting the physical accuracy of the results.
For a stability study, add a velocity-spectrum diagnostic and check that the resolved inertial range is unchanged as the tail filter activates. The saved sample output is vorticity magnitude, so it can expose cutoff pile-up but is not by itself a kinetic-energy spectrum or a -5/3-slope validation. The default controls are intentionally conservative and should be retuned with resolution and Reynolds number rather than used as a universal turbulence model.
Animated Results
The program outputs 51 field sets of vorticity magnitude (i.e., \( \lVert \nabla \times \mathbf{V} \rVert \)) from \( t = t_0 = 0 \) to \( t = t_f = 50 \) at intervals of \( 1 \). The following Python matplotlib code helps visualize the vortex merging process, using a 3-dimensional iso-contour (i.e., iso-surface) of vorticity magnitude.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.patches as mpatches
from skimage.measure import marching_cubes as mc
from scipy.interpolate import RegularGridInterpolator
# Load static data outside the animate function
NR = 200; NP = 128; NZ = 128
dt = 1.E-2; dn = 1E2
coords_r = np.loadtxt('../output/dat/coords_r.dat', skiprows=1, max_rows=NR)
coords_p = np.loadtxt('../output/dat/coords_p.dat', skiprows=1, max_rows=NP+1)
coords_z = np.loadtxt('../output/dat/coords_z.dat', skiprows=1, max_rows=NZ+1)
# Initialize the figure and 3D axis once
fig = plt.figure(figsize=(10, 8), facecolor='none')
ax = fig.add_subplot(projection='3d', facecolor='none')
xlim = 5; ylim = 5; zlim = 2*np.pi
ax.set_xlabel("Horizontal (x)"); ax.set_ylabel("Vertical (y)"); ax.set_zlabel("Axial (z)")
ax.set_xlim(-xlim, xlim); ax.set_ylim(-ylim, ylim); ax.set_zlim(0, zlim)
# Animation function
def animate(frame):
nn = frame * dn
f_val = []
for z in range(NZ):
f_val.append(np.loadtxt('../output/fld/vormagfld' + f"{nn:06n}" + '_PPP', skiprows=z+1+(NR+8)*z, max_rows=NR))
f_val.append(f_val[0])
f_val = np.swapaxes(np.swapaxes(np.array(f_val), 0, 1), 1, 2)
f_val = f_val[0:NR, 0:NP+1, 0:NZ+1]
f_val[:,NP,:] = f_val[:,0,:] # phi of 2pi == phi of 0
f_val[:,:,NZ] = f_val[:,:,0] # z of ZLEN == z of 0
# Interpolation grid in Cartesian coordinates
cartesian_x = np.linspace(-xlim, xlim, 100) # Number of points can be adjusted
cartesian_y = np.linspace(-ylim, ylim, 100)
cartesian_z = np.linspace(0, zlim, 100)
X, Y, Z = np.meshgrid(cartesian_x, cartesian_y, cartesian_z, indexing="ij")
# Interpolation to create a smooth isosurface
R = np.sqrt(X**2 + Y**2)
P = np.arctan2(Y, X) % (2 * np.pi) # Wrap angles to [0, 2*pi]
interpolator = RegularGridInterpolator( (coords_r, coords_p, coords_z),
f_val, bounds_error=False, fill_value=0 )
cylindrical_points = np.stack((R.ravel(), P.ravel(), Z.ravel()), axis=-1)
f_cartesian = interpolator(cylindrical_points).reshape(X.shape)
# Isosurface extraction
iso_value = 1.0
verts, faces, normals, values = mc(f_cartesian, level=iso_value)
# Transform vertices directly in Cartesian coordinates
verts_cartesian = verts * np.array(
[2 * xlim / (len(cartesian_x) - 1), 2 * ylim / (len(cartesian_y) - 1), zlim / (len(cartesian_z) - 1)]
) + np.array([-xlim, -ylim, 0])
# Add labels and adjust view
ax.clear()
ax.set_xlabel("Horizontal (x)")
ax.set_ylabel("Vertical (y)")
ax.set_zlabel("Axial (z)")
ax.set_xlim(-xlim, xlim)
ax.set_ylim(-ylim, ylim)
ax.set_zlim(0, zlim)
ax.set_title(f"Isosurface (Vorticity Magnitude={iso_value}) at t = {dt*nn}")
# Plot the surface with shading and transparency
ax.plot_trisurf( verts_cartesian[:, 0], verts_cartesian[:, 1],
faces, verts_cartesian[:, 2],
lw=.5, alpha=.9, shade=True )
plt.tight_layout()
# Run the animation
anim = animation.FuncAnimation(fig, animate, frames=range(0, 51), interval=50)
anim.save('animation.mp4', writer='ffmpeg', fps=15, dpi=150)
Bonus
Adding a slight perturbation in the axial direction can significantly alter the dynamics of vortex evolution. Consider superposing a sinusoidally perturbed axial jet initially located at the center, with a velocity field defined as:
\[\mathbf{V} (r, \phi, z) = \exp(-r^2) \left( \frac{1 + \cos z}{2} \right) \hat{\mathbf{e}}_{z}.\]As a result, the vortex system undergoes more drastic evolution, with the emergence of pronounced axially varying flow structures. This scenario demonstrates that slight perturbations can significantly affect the flow dynamics, leading to complex vortical flow patterns.