diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml index 77a39486..2f1cf238 100644 --- a/.github/workflows/linux.yml +++ b/.github/workflows/linux.yml @@ -162,16 +162,22 @@ jobs: steps: - uses: actions/checkout@v3 - name: System Dependencies - run: .github/workflows/dependencies/dependencies_gcc10.sh + run: | + .github/workflows/dependencies/dependencies_gcc10.sh + sudo apt-get install -y python3-setuptools + python3 -m pip install --user scipy + python3 -m pip install --user matplotlib + python3 -m pip install --user numpy - name: Repo Dependencies run: Utils/CloneDeps.sh - name: GenerateTurbFile - env: - AMREX_HOME: ${GITHUB_WORKSPACE}/Submodules/PelePhysics/Submodules/amrex - working-directory: ./Exec/RegTests/TurbInflow/TurbFileHIT + working-directory: ./Submodules/PelePhysics/Support/TurbFileHIT run: | + ./gen_hit_ic.py -k0 4 -N 32 -Nk 256 make -j 2 COMP=gnu - ./PeleTurb3d.gnu.ex input hit_file=../../HITDecay/hit_ic_4_32.dat input_ncell=32 amrex.abort_on_unused_inputs=1 + ./PeleTurb3d.gnu.ex input hit_file=hit_ic_4_32.dat input_ncell=32 amrex.abort_on_unused_inputs=1 + mkdir ${GITHUB_WORKSPACE}/Exec/RegTests/TurbInflow/TurbFileHIT + cp -r TurbTEST ${GITHUB_WORKSPACE}/Exec/RegTests/TurbInflow/TurbFileHIT - name: Build working-directory: ./Exec/RegTests/TurbInflow/ run: | diff --git a/Docs/sphinx/manual/LMeXControls.rst b/Docs/sphinx/manual/LMeXControls.rst index 4c5202d3..0c70e5bd 100644 --- a/Docs/sphinx/manual/LMeXControls.rst +++ b/Docs/sphinx/manual/LMeXControls.rst @@ -36,6 +36,15 @@ Computational domain definition peleLM.lo_bc = Interior Interior Inflow peleLM.hi_bc = Interior Interior Inflow +If specifying boundaries as ``Inflow``, the bcnormal function must be defined +in the ``pelelmex_prob.H`` file for the case to define the inflow conditions. +``Inflow`` boundaries may also be augmented with spatially and temporally +varying turbulent fluctuations using the ``TurbInflow`` utility from +PelePhysics. See the ``Exec/RegTests/TurbInflow`` test for an example of how +to use this capability and the +`PelePhysics documentation `_ +for the relevant input file flags. + Grid/AMR parameters ------------------- @@ -239,7 +248,7 @@ PeleLMeX algorithm peleLM.gravity = 0.0 0.0 -9.81 # [OPT, DEF=Vec{0.0}] Gravity vector [MKS] peleLM.gradP0 = 0.0 0.0 10.0 # [OPT, DEF=Vec{0.0}] Average background pressure gradient [Pa/m] peleLM.do_periodic_channel = 0 # [OPT, DEF= 0] Add an automatic pressure gradient to maintain initial condition mass flow rate in periodic channel - peleLM.periodic_channel_dir = 2 # [OPT, DEF= -1] Required if do_periodic_channel != 0. Direction to apply pressure gradient. + peleLM.periodic_channel_dir = 2 # [OPT, DEF= -1] Required if do_periodic_channel != 0. Direction to apply pressure gradient. peleLM.closed_chamber = 0 # [OPT] Override the automatic detection of closed chamber (based on Outflow(s)) peleLM.floor_species = 0 # [OPT, DEF=0] Crudely enforce mass fraction positivity peleLM.deltaT_verbose = 0 # [OPT, DEF=0] Verbose of the deltaT iterative solve algorithm @@ -258,7 +267,7 @@ Transport coeffs and LES peleLM.Prandtl = 0.7 # [OPT, DEF=0.7] If fixed_Pr or doing LES, specifies the Prandtl number peleLM.Schmidt = 0.7 # [OPT, DEF=0.7] If doing LES, specifies the Schmidt number peleLM.Lewis = 1.0 # [OPT, DEF=1.0] If fixed_Le, specifies the Lewis number - + peleLM.les_model = "None" # [OPT, DEF="None"] Model to compute turbulent viscosity: None, Smagorinsky, WALE, Sigma peleLM.les_cs_smag = 0.18 # [OPT, DEF=0.18] If using Smagorinsky LES model, provides model coefficient peleLM.les_cm_wale = 0.60 # [OPT, DEF=0.60] If using WALE LES model, provides model coefficient @@ -280,7 +289,7 @@ Chemistry integrator cvode.solve_type = denseAJ_direct # [OPT, DEF=GMRES] Linear solver employed for CVODE Newton direction cvode.max_order = 4 # [OPT, DEF=2] Maximum order of the BDF method in CVODE cvode.max_substeps = 10000 # [OPT, DEF=10000] Maximum number of substeps for the linear solver in CVODE - + Note that the last five parameters belong to the Reactor class of PelePhysics but are specified here for completeness. In particular, CVODE is the adequate choice of integrator to tackle PeleLMeX large time step sizes. Several linear solvers are available depending on whether or not GPU are employed: on CPU, `dense_direct` is a finite-difference direct solver, `denseAJ_direct` is an analytical-jacobian direct solver (preferred choice), `sparse_direct` is an analytical-jacobian sparse direct solver based on the KLU library and `GMRES` is a matrix-free iterative solver; on GPU `GMRES` is a matrix-free iterative solver (available on all the platforms), `sparse_direct` is a batched block-sparse direct solve based on NVIDIA's cuSparse (only with CUDA), `magma_direct` is a batched block-dense direct solve based on the MAGMA library (available with CUDA and HIP. Different `cvode.solve_type` should be tried before increasing the `cvode.max_substeps`. Embedded Geometry diff --git a/Exec/RegTests/HITDecay/README b/Exec/RegTests/HITDecay/README index e25d1e91..8f706cbd 100644 --- a/Exec/RegTests/HITDecay/README +++ b/Exec/RegTests/HITDecay/README @@ -3,25 +3,8 @@ originally implemented in PeleC and used to study compressible turbulence in: E. Motheau and J. Wakefield, "Investigation of finite-volume methods to capture shocks and turbulence spectra in compressible flows", Commun. in Appl. Math. and Comput. Sci, 15-1 (2020), 1--36. https://arxiv.org/abs/1902.06665 -The script gen_hit_ic.py can be used to generate an initial solution. +The script gen_hit_ic.py from PelePhysics/Support/TurbFileHIT can be used to generate an initial solution. Here, there is no scaling by the turbulent Mach number and the density is kept constant, so that the Taylor-scale Reynolds number is defined by Re=0.5 urms / nu where umrs can be set in the input file (prob.urms0=...) and nu is computed from the mixture composition. -At the beginning of a simulation, a file initialConditions.txt is created to hold the values of urms0, lambda0 and tau, which will be used for post-processing. - -The post-processing chain is composed of two major parts: computing the temporal evolution of the turbulent kinetic energy, the square of magnitude of vorticity as well as the square of the divergence. Then the turbulent spectra is computed for the last plotfile found. - -IMPORTANT: -The following tools must be pre-compiled and the executable located at the root of this directory, where the post.sh script is located: -AmrDeriveSpectrum (https://github.com/AMReX-Astro/AmrDeriveSpectrum) -AugmentPlotfile (located in amrex/Tools/C_util/AugmentPlotfile) - -To proceed to an analysis, please create a directory and copy inside all the plotfiles as well as the initialConditions.txt file. Then execute the post.sh script as in the following example: - -./PeleLMeX3d.gnu.MPI.ex inputs.3d -mkdir ANALYSIS -cp -r plt_0000* ANALYSIS/ -cp initialConditions.txt ANALYSIS/ -./post.sh ANALYSIS/ 1 - -The output csv files will be located in the ANALYSIS directory. +At the beginning of a simulation, a file initialConditions.txt is created to hold the values of urms0, lambda0 and tau, which will be used for post-processing. \ No newline at end of file diff --git a/Exec/RegTests/HITDecay/gen_hit_ic.py b/Exec/RegTests/HITDecay/gen_hit_ic.py deleted file mode 100755 index a3fc013a..00000000 --- a/Exec/RegTests/HITDecay/gen_hit_ic.py +++ /dev/null @@ -1,463 +0,0 @@ -#!/usr/bin/env python -# -# Generate a table of the velocity fluctuations for the homogeneous -# isotropic turbulence case at a specific k0 (default to 4) -# -# Order of operations: -# 1. velocity fluctuations generated on a 512^3 grid in wavenumber space -# 2. Coefficients associated to wavenumbers that cannot be represented on -# the desired grid are set to 0 (sharp wavenumber cutoff) -# 3. inverse Fourier transform of the velocity fluctuations (512^3 grid) -# 4. velocity fluctuations resampled on the desired grid (N^3) -# -# The velocity fluctuations are normalized by urms0 so to get the -# actual velocity fluctuations, one must multiply these velocities by -# the appropriate urms0. -# -# - -# ======================================================================== -# -# Imports -# -# ======================================================================== -import argparse -import sys -import time -from datetime import timedelta -import numpy as np -import scipy.interpolate as spi -import matplotlib as mpl - -mpl.use("Agg") -import matplotlib.pyplot as plt - - -# ======================================================================== -# -# Parse arguments -# -# ======================================================================== -parser = argparse.ArgumentParser( - description="Generate the velocity fluctuations for the HIT IC" -) -parser.add_argument( - "-k0", help="Wave number containing highest energy", type=float, default=4.0 -) -parser.add_argument("-N", help="Resolution", type=int, default=16) -parser.add_argument( - "-s", "--seed", help="Random number generator seed", type=int, default=42 -) -parser.add_argument( - "-p", "--plot", help="Save a plot of the x-velocity", action="store_true" -) -args = parser.parse_args() - -# =============================================================================== -# -# Some defaults variables -# -# =============================================================================== -plt.rc("text", usetex=True) -plt.rc("font", family="serif", serif="Times") -cmap_med = [ - "#F15A60", - "#7AC36A", - "#5A9BD4", - "#FAA75B", - "#9E67AB", - "#CE7058", - "#D77FB4", - "#737373", -] -cmap = [ - "#EE2E2F", - "#008C48", - "#185AA9", - "#F47D23", - "#662C91", - "#A21D21", - "#B43894", - "#010202", -] -dashseq = [ - (None, None), - [10, 5], - [10, 4, 3, 4], - [3, 3], - [10, 4, 3, 4, 3, 4], - [3, 3], - [3, 3], -] -markertype = ["s", "d", "o", "p", "h"] - -# ======================================================================== -# -# Function definitions -# -# ======================================================================== -def div0(a, b): - """ Ignore division by 0, just replace it by 0, - - From: http://stackoverflow.com/questions/26248654/numpy-return-0-with-divide-by-zero - e.g. div0( [-1, 0, 1], 0 ) -> [0, 0, 0] - """ - with np.errstate(divide="ignore", invalid="ignore"): - c = np.true_divide(a, b) - c[~np.isfinite(c)] = 0 # -inf inf NaN - return c - - -# ======================================================================== -def abs2(x): - """This is equivalent to np.abs(x)**2 or x*np.conj(x) - - To make it faster, add this right before the function definition - import numba - @numba.vectorize([numba.float64(numba.complex128),numba.float32(numba.complex64)]) - """ - return x.real ** 2 + x.imag ** 2 - - -# ======================================================================== -# -# Main -# -# ======================================================================== - -# Timer -start = time.time() - -# ======================================================================== -# 1. velocity fluctuations generated on a 512^3 grid in wavenumber space - -# Dimension of the large cube -N = 512 -halfN = int(round(0.5 * N)) -xs = 0 -xe = 2.0 * np.pi * 0.01 -L = xe - xs -dx = L / N - -# Only work if N and args.N are even -if not ((args.N % 2 == 0) and N % 2 == 0): - print("N or args.N is not even. Exiting") - sys.exit(1) - -# Get cell centered values and meshed grid -x = np.linspace(xs, xe, N + 1) -xc = (x[1:] + x[:-1]) / 2 # get cell center coordinates -X, Y, Z = np.meshgrid(xc, xc, xc, indexing="ij") - -# Get the wave numbers and associated quantities -k = np.concatenate((np.arange(halfN), np.arange(-halfN, 0, 1)), axis=0) -khalf = np.arange(halfN + 1) -k1, k2, k3 = np.meshgrid(k, k, khalf, indexing="ij") -kmag = np.sqrt(k1 ** 2 + k2 ** 2 + k3 ** 2) -k12 = np.sqrt(k1 ** 2 + k2 ** 2) -k1k12 = div0(k1, k12) -k2k12 = div0(k2, k12) -k3kmag = div0(k3, kmag) -k12kmag = div0(k12, kmag) - -# Generate data - -# # Toy Fourier data corresponding to uo = cos(X) * cos(2*Y) * cos(3*Z) -# uo = np.cos(X) * np.cos(2*Y) * np.cos(3*Z) -# uf = np.fft.rfftn(uo) -# vf = np.copy(uf) -# wf = np.copy(uf) - -# Energy spectrum -Ek = ( - 16.0 - * np.sqrt(2.0 / np.pi) - * (kmag ** 4) - / (args.k0 ** 5) - * np.exp(-2.0 * (kmag ** 2) / (args.k0 ** 2)) -) - -# Draw random numbers -np.random.seed(args.seed) -phi1 = np.random.uniform(0, 2 * np.pi, np.shape(kmag)) -phi2 = np.random.uniform(0, 2 * np.pi, np.shape(kmag)) -phi3 = np.random.uniform(0, 2 * np.pi, np.shape(kmag)) - -# the random quantities -prefix = np.sqrt(2.0 * div0(Ek, 4.0 * np.pi * (kmag ** 2))) -a = prefix * np.exp(1j * phi1) * np.cos(phi3) -b = prefix * np.exp(1j * phi2) * np.sin(phi3) - -# the random velocities -uf = k2k12 * a + k1k12 * k3kmag * b -vf = k2k12 * k3kmag * b - k1k12 * a -wf = -k12kmag * b - -# Impose the 3D spherical symmetry (to ensure we have a real signal) -# equiv: uf[-l,-m,0] = np.conj(uf[ l, m,0]) for l=0..N/2 and m=0..N/2 -uf[N:halfN:-1, N:halfN:-1, 0] = np.conj(uf[1:halfN, 1:halfN, 0]) -# symmetry on first column -uf[N:halfN:-1, 0, 0] = np.conj(uf[1:halfN, 0, 0]) -# symmetry on first row -uf[0, N:halfN:-1, 0] = np.conj(uf[0, 1:halfN, 0]) -# symmetry about the (halfN,halfN) element -uf[halfN - 1 : 0 : -1, N : halfN - 1 : -1, 0] = np.conj( - uf[halfN + 1 : N, 1 : halfN + 1, 0] -) - -vf[N:halfN:-1, N:halfN:-1, 0] = np.conj(vf[1:halfN, 1:halfN, 0]) -vf[halfN - 1 : 0 : -1, N : halfN - 1 : -1, 0] = np.conj( - vf[halfN + 1 : N, 1 : halfN + 1, 0] -) -vf[N:halfN:-1, 0, 0] = np.conj(vf[1:halfN, 0, 0]) -vf[0, N:halfN:-1, 0] = np.conj(vf[0, 1:halfN, 0]) - -wf[N:halfN:-1, N:halfN:-1, 0] = np.conj(wf[1:halfN, 1:halfN, 0]) -wf[halfN - 1 : 0 : -1, N : halfN - 1 : -1, 0] = np.conj( - wf[halfN + 1 : N, 1 : halfN + 1, 0] -) -wf[N:halfN:-1, 0, 0] = np.conj(wf[1:halfN, 0, 0]) -wf[0, N:halfN:-1, 0] = np.conj(wf[0, 1:halfN, 0]) - -# Normalize. Because we are generating the data in wavenumber space, -# we have to multiply by N**3 because in the definition of the numpy -# ifftn there is a 1/N**n. -uf = uf * N ** 3 -vf = vf * N ** 3 -wf = wf * N ** 3 - -# # Quick check on energy content (make sure you add both the current -# # contribution and the one we are neglecting because we are assuming -# # real input data) -# print('Energy = int E(k) dk = 0.5 * int (uf**2 + vf**2 wf**2) dk1 dk2 dk3 = {0:.10f} ~= 3/2'.format( -# (np.sum(abs2(uf ) + abs2(vf ) + abs2(wf )) + -# np.sum(abs2(uf[:,:,1:-1]) + abs2(vf[:,:,1:-1]) + abs2(wf[:,:,1:-1]))) -# * 0.5 / N**6)) - -# if plotting, save the original field (before filtering) -if args.plot: - uo = np.fft.irfftn(uf) - Eko = ( - 16.0 - * np.sqrt(2.0 / np.pi) - * (khalf ** 4) - / (args.k0 ** 5) - * np.exp(-2.0 * (khalf ** 2) / (args.k0 ** 2)) - ) - - # Get the spectrum from 3D velocity field - kbins = np.arange(1, halfN + 1) - Nbins = len(kbins) - whichbin = np.digitize(kmag.flat, kbins) - ncount = np.bincount(whichbin) - - KI = (abs2(uf) + abs2(vf) + abs2(wf)) * 0.5 / N ** 6 - KI[:, :, 1:-1] += ( - (abs2(uf[:, :, 1:-1]) + abs2(vf[:, :, 1:-1]) + abs2(wf[:, :, 1:-1])) - * 0.5 - / N ** 6 - ) - - Eku = np.zeros(len(ncount) - 1) - for n in range(1, len(ncount)): - Eku[n - 1] = np.sum(KI.flat[whichbin == n]) - - ku = 0.5 * (kbins[0 : Nbins - 1] + kbins[1:Nbins]) + 1 - Eku = Eku[1:Nbins] - - -# ======================================================================== -# 2. Coefficients associated to wavenumbers that cannot be represented -# on the desired grid are set to 0 (sharp wavenumber cutoff) -kmagc = 0.5 * args.N -uf[kmag > kmagc] = 0.0 -vf[kmag > kmagc] = 0.0 -wf[kmag > kmagc] = 0.0 - - -# ======================================================================== -# 3. inverse Fourier transform of the velocity fluctuations (512^3 grid) -u = np.fft.irfftn(uf, s=(N, N, N)) -v = np.fft.irfftn(vf, s=(N, N, N)) -w = np.fft.irfftn(wf, s=(N, N, N)) - -# Another energy content check -print( - "Energy = 1/V * int E(x,y,z) dV = 0.5/V * int (u**2 + v**2 + w**2) dx dy dz = {0:.10f} ~= 3/2".format( - np.sum(u ** 2 + v ** 2 + w ** 2) * 0.5 * (dx / L) ** 3 - ) -) - -# # Enstrophy check -# _, dudy, dudz = np.gradient(u, dx) -# dvdx, _, dvdz = np.gradient(v, dx) -# dwdx, dwdy, _ = np.gradient(w, dx) -# wx = dwdy-dvdz -# wy = dudz-dwdx -# wz = dvdx-dudy -# lambda0 = 2.0/args.k0 -# print('Enstrophy = 0.5/V * int (wx**2 + wy**2 + wz**2) dx dy dz= -# {0:.10f} ~= '.format(np.sum(wx**2+wy**2+wz**2) * 0.5 * (dx/L)**3 * -# lambda0**2)) - -# ======================================================================== -# 4. velocity fluctuations re-sampled on the desired grid (N^3) -xr = np.linspace(xs, xe, args.N + 1) -xrc = (xr[1:] + xr[:-1]) / 2 -Xr, Yr, Zr = np.meshgrid(xrc, xrc, xrc, indexing="ij") - -Xr = Xr.reshape(-1, order="F") -Yr = Yr.reshape(-1, order="F") -Zr = Zr.reshape(-1, order="F") - -ur = spi.interpn((xc, xc, xc), u, (Xr, Yr, Zr), method="linear") -vr = spi.interpn((xc, xc, xc), v, (Xr, Yr, Zr), method="linear") -wr = spi.interpn((xc, xc, xc), w, (Xr, Yr, Zr), method="linear") - - -# ======================================================================== -# Save the data in Fortran ordering -fname = "hit_ic_{0:d}_{1:d}.dat".format(int(args.k0), args.N) -data = np.vstack((Xr, Yr, Zr, ur, vr, wr)).T -np.savetxt(fname, data, fmt="%.18e", delimiter=",", header="x, y, z, u, v, w") - - -# ======================================================================== -# plot (only u fluctuations) -if args.plot: - import matplotlib as mpl - - mpl.use("Agg") - import matplotlib.pyplot as plt - - datmin = u.min() - datmax = u.max() - # print("min/max u =",datmin,datmax) - - # Original data - # transpose and origin change bc I used meshgrid with ij and not xy - fig, ax = plt.subplots(nrows=3, ncols=3, figsize=(14, 14)) - ax[0, 0].imshow( - uo[:, :, 0].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[0, 0].set_title("Original data (x,y)") - ax[0, 1].imshow( - uo[:, 0, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[0, 1].set_title("Original data (x,z)") - ax[0, 2].imshow( - uo[0, :, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[0, 2].set_title("Original data (y,z)") - - # Filtered original data - ax[1, 0].imshow( - u[:, :, 0].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[1, 0].set_title("Filtered original data (x,y)") - ax[1, 1].imshow( - u[:, 0, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[1, 1].set_title("Filtered original data (x,z)") - ax[1, 2].imshow( - u[0, :, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[1, 2].set_title("Filtered original data (y,z)") - - # Downsampled filtered data - ur = ur.reshape(args.N, args.N, args.N, order="F") - ax[2, 0].imshow( - ur[:, :, 0].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[2, 0].set_title("Downsampled data (x,y)") - ax[2, 1].imshow( - ur[:, 0, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[2, 1].set_title("Downsampled data (x,z)") - ax[2, 2].imshow( - ur[0, :, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[2, 2].set_title("Downsampled data (y,z)") - - plt.savefig("hit_ic_u_{0:d}_{1:d}.png".format(int(args.k0), args.N), format="png") - - # Fourier coefficients of original data - fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(14, 8)) - ax[0, 0].imshow(np.real(uf[:, :, 0].T), origin="lower", cmap="RdBu_r") - ax[0, 0].set_title("Real Fourier coefficients (x,y)") - ax[0, 1].imshow(np.real(uf[:, 0, :].T), origin="lower", cmap="RdBu_r") - ax[0, 1].set_title("Real Fourier coefficients (x,z)") - ax[0, 2].imshow(np.real(uf[0, :, :].T), origin="lower", cmap="RdBu_r") - ax[0, 2].set_title("Real Fourier coefficients (y,z)") - ax[1, 0].imshow(np.imag(uf[:, :, 0].T), origin="lower", cmap="RdBu_r") - ax[1, 0].set_title("Imag Fourier coefficients (x,y)") - ax[1, 1].imshow(np.imag(uf[:, 0, :].T), origin="lower", cmap="RdBu_r") - ax[1, 1].set_title("Imag Fourier coefficients (x,z)") - ax[1, 2].imshow(np.imag(uf[0, :, :].T), origin="lower", cmap="RdBu_r") - ax[1, 2].set_title("Imag Fourier coefficients (y,z)") - plt.savefig("hit_ic_uf_{0:d}_{1:d}.png".format(int(args.k0), args.N), format="png") - - # Spectrum - plt.figure(20) - ax = plt.gca() - p = plt.loglog(khalf, Eko, color=cmap[-1], lw=2) - p[0].set_dashes(dashseq[0]) - p = plt.loglog(ku, Eku, color=cmap[0], lw=2) - p[0].set_dashes(dashseq[1]) - plt.ylim([1e-16, 10]) - plt.xlabel(r"$k$", fontsize=22, fontweight="bold") - plt.ylabel(r"$E(k)$", fontsize=22, fontweight="bold") - plt.setp(ax.get_xmajorticklabels(), fontsize=18, fontweight="bold") - plt.setp(ax.get_ymajorticklabels(), fontsize=18, fontweight="bold") - plt.savefig( - "hit_ic_spectrum_{0:d}_{1:d}.png".format(int(args.k0), args.N), format="png" - ) - -# output timer -end = time.time() - start -print("Elapsed time " + str(timedelta(seconds=end)) + " (or {0:f} seconds)".format(end)) diff --git a/Exec/RegTests/TurbInflow/README.md b/Exec/RegTests/TurbInflow/README.md index d01cbf2d..c4918b13 100644 --- a/Exec/RegTests/TurbInflow/README.md +++ b/Exec/RegTests/TurbInflow/README.md @@ -1,3 +1,5 @@ ## TurbInflow -Injection of a 3D precomputed turbulent velocity field at the boundary of a PeleLMeX simulation. Testing the -turbulence injection functionality. +Injection of a 3D precomputed turbulent velocity field at the boundary of a +PeleLMeX simulation. Testing the turbulence injection functionality. Use the +TurbFileHIT support code from PelePhysics to generate the turbulent +boundary data. diff --git a/Exec/RegTests/TurbInflow/TurbFileHIT/GNUmakefile b/Exec/RegTests/TurbInflow/TurbFileHIT/GNUmakefile deleted file mode 100644 index 57cc1e91..00000000 --- a/Exec/RegTests/TurbInflow/TurbFileHIT/GNUmakefile +++ /dev/null @@ -1,27 +0,0 @@ -TOP = ../../../.. -AMREX_HOME ?= ${TOP}/amrex - -# Local -EBASE = PeleTurb - -# AMReX -DIM = 3 -DEBUG = FALSE -PRECISION = DOUBLE -TINY_PROFILE = FALSE - -# Compilation -COMP = gnu -USE_MPI = FALSE - -include $(AMREX_HOME)/Tools/GNUMake/Make.defs - -Bpack += ./Make.package - -Pdirs := Base Boundary -Bpack += $(foreach dir, $(Pdirs), $(AMREX_HOME)/Src/$(dir)/Make.package) -Blocs += $(foreach dir, $(Pdirs), $(AMREX_HOME)/Src/$(dir)) -include $(Bpack) -INCLUDE_LOCATIONS += $(Blocs) - -include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/Exec/RegTests/TurbInflow/TurbFileHIT/HITData.H b/Exec/RegTests/TurbInflow/TurbFileHIT/HITData.H deleted file mode 100644 index 73857782..00000000 --- a/Exec/RegTests/TurbInflow/TurbFileHIT/HITData.H +++ /dev/null @@ -1,21 +0,0 @@ -#ifndef HITDATA_H -#define HITDATA_H - -#include -#include - -struct HITData -{ - int input_ncell = 0; - amrex::Real urms0 = 1.0; - amrex::Real uin_norm = 1.0; - - amrex::Real Linput = 0.0; - - amrex::Real* d_uinput = nullptr; - amrex::Real* d_vinput = nullptr; - amrex::Real* d_winput = nullptr; - amrex::Real* d_xarray = nullptr; - amrex::Real* d_xdiff = nullptr; -}; -#endif diff --git a/Exec/RegTests/TurbInflow/TurbFileHIT/Make.package b/Exec/RegTests/TurbInflow/TurbFileHIT/Make.package deleted file mode 100644 index 04d5df0e..00000000 --- a/Exec/RegTests/TurbInflow/TurbFileHIT/Make.package +++ /dev/null @@ -1,2 +0,0 @@ -CEXE_headers += HITData.H main_K.H Utilities.H -CEXE_sources += main.cpp Utilities.cpp diff --git a/Exec/RegTests/TurbInflow/TurbFileHIT/README.md b/Exec/RegTests/TurbInflow/TurbFileHIT/README.md deleted file mode 100644 index c61c73f7..00000000 --- a/Exec/RegTests/TurbInflow/TurbFileHIT/README.md +++ /dev/null @@ -1,14 +0,0 @@ -This code enable generating synthetic HIT files usable with PelePhysics -TurbInflow capabilities. - -First generate the data using the python script: -./gen_hit_ic.py -k0 4 -N 128 - -To generate a synthetic HIT field discretized with 128 cells and most energetic eddies -at a wave number of 4. - -Then compile the C++ executable (AMReX needed): -make - -And the executable to generature the turbfile (adapt the input file to your needs): -./PeleTurb3d.gnu.ex input diff --git a/Exec/RegTests/TurbInflow/TurbFileHIT/Utilities.H b/Exec/RegTests/TurbInflow/TurbFileHIT/Utilities.H deleted file mode 100644 index 1306fc86..00000000 --- a/Exec/RegTests/TurbInflow/TurbFileHIT/Utilities.H +++ /dev/null @@ -1,72 +0,0 @@ -#ifndef _UTILITIES_H -#define _UTILITIES_H - -#include - -AMREX_FORCE_INLINE -std::string -read_file(std::ifstream& in) -{ - return static_cast( - std::stringstream() << in.rdbuf()) - .str(); -} - -void read_binary( - const std::string& iname, - const size_t nx, - const size_t ny, - const size_t nz, - const size_t ncol, - amrex::Vector& data); - -void read_csv( - const std::string& iname, - const size_t nx, - const size_t ny, - const size_t nz, - amrex::Vector& data); - -// ----------------------------------------------------------- -// Search for the closest index in an array to a given value -// using the bisection technique. -// INPUTS/OUTPUTS: -// xtable(0:n-1) => array to search in (ascending order) -// n => array size -// x => x location -// idxlo <=> output st. xtable(idxlo) <= x < xtable(idxlo+1) -// ----------------------------------------------------------- -AMREX_GPU_DEVICE -AMREX_FORCE_INLINE -void -locate(const amrex::Real* xtable, const int n, const amrex::Real& x, int& idxlo) -{ - // If x is out of bounds, return boundary index - if (x >= xtable[n - 1]) { - idxlo = n - 1; - return; - } - if (x <= xtable[0]) { - idxlo = 0; - return; - } - - // Do the bisection - idxlo = 0; - int idxhi = n - 1; - bool notdone = true; - while (notdone) { - if (idxhi - idxlo <= 1) { - notdone = false; - } else { - const int idxmid = (idxhi + idxlo) / 2; - if (x >= xtable[idxmid]) { - idxlo = idxmid; - } else { - idxhi = idxmid; - } - } - } -} - -#endif diff --git a/Exec/RegTests/TurbInflow/TurbFileHIT/Utilities.cpp b/Exec/RegTests/TurbInflow/TurbFileHIT/Utilities.cpp deleted file mode 100644 index ea77b032..00000000 --- a/Exec/RegTests/TurbInflow/TurbFileHIT/Utilities.cpp +++ /dev/null @@ -1,87 +0,0 @@ -#include "Utilities.H" - -// ----------------------------------------------------------- -// Read a binary file -// INPUTS/OUTPUTS: -// iname => filename -// nx => input resolution -// ny => input resolution -// nz => input resolution -// data <= output data -// ----------------------------------------------------------- -void -read_binary( - const std::string& iname, - const size_t nx, - const size_t ny, - const size_t nz, - const size_t ncol, - amrex::Vector& data /*needs to be double*/) -{ - std::ifstream infile(iname, std::ios::in | std::ios::binary); - if (not infile.is_open()) { - amrex::Abort("Unable to open input file " + iname); - } - - for (size_t i = 0; i < nx * ny * nz * ncol; i++) { - infile.read(reinterpret_cast(&data[i]), sizeof(data[i])); - } - infile.close(); -} - -// ----------------------------------------------------------- -// Read a csv file -// INPUTS/OUTPUTS: -// iname => filename -// nx => input resolution -// ny => input resolution -// nz => input resolution -// data <= output data -// ----------------------------------------------------------- -void -read_csv( - const std::string& iname, - const size_t nx, - const size_t ny, - const size_t nz, - amrex::Vector& data) -{ - std::ifstream infile(iname, std::ios::in); - const std::string memfile = read_file(infile); - if (not infile.is_open()) { - amrex::Abort("Unable to open input file " + iname); - } - infile.close(); - std::istringstream iss(memfile); - - // Read the file - size_t nlines = 0; - std::string firstline; - std::string line; - std::getline(iss, firstline); // skip header - while (getline(iss, line)) { - ++nlines; - } - - // Quick sanity check - if (nlines != nx * ny * nz) { - amrex::Abort( - "Number of lines in the input file (= " + std::to_string(nlines) + - ") does not match the input resolution (=" + std::to_string(nx) + ")"); - } - - // Read the data from the file - iss.clear(); - iss.seekg(0, std::ios::beg); - std::getline(iss, firstline); // skip header - int cnt = 0; - while (std::getline(iss, line)) { - std::istringstream linestream(line); - std::string value; - while (getline(linestream, value, ',')) { - std::istringstream sinput(value); - sinput >> data[cnt]; - cnt++; - } - } -} diff --git a/Exec/RegTests/TurbInflow/TurbFileHIT/gen_hit_ic.py b/Exec/RegTests/TurbInflow/TurbFileHIT/gen_hit_ic.py deleted file mode 100755 index 3d847237..00000000 --- a/Exec/RegTests/TurbInflow/TurbFileHIT/gen_hit_ic.py +++ /dev/null @@ -1,463 +0,0 @@ -#!/usr/bin/env python -# -# Generate a table of the velocity fluctuations for the homogeneous -# isotropic turbulence case at a specific k0 (default to 4) -# -# Order of operations: -# 1. velocity fluctuations generated on a 512^3 grid in wavenumber space -# 2. Coefficients associated to wavenumbers that cannot be represented on -# the desired grid are set to 0 (sharp wavenumber cutoff) -# 3. inverse Fourier transform of the velocity fluctuations (512^3 grid) -# 4. velocity fluctuations resampled on the desired grid (N^3) -# -# The velocity fluctuations are normalized by urms0 so to get the -# actual velocity fluctuations, one must multiply these velocities by -# the appropriate urms0. -# -# - -# ======================================================================== -# -# Imports -# -# ======================================================================== -import argparse -import sys -import time -from datetime import timedelta -import numpy as np -import scipy.interpolate as spi -import matplotlib as mpl - -mpl.use("Agg") -import matplotlib.pyplot as plt - - -# ======================================================================== -# -# Parse arguments -# -# ======================================================================== -parser = argparse.ArgumentParser( - description="Generate the velocity fluctuations for the HIT IC" -) -parser.add_argument( - "-k0", help="Wave number containing highest energy", type=float, default=4.0 -) -parser.add_argument("-N", help="Resolution", type=int, default=16) -parser.add_argument( - "-s", "--seed", help="Random number generator seed", type=int, default=42 -) -parser.add_argument( - "-p", "--plot", help="Save a plot of the x-velocity", action="store_true" -) -args = parser.parse_args() - -# =============================================================================== -# -# Some defaults variables -# -# =============================================================================== -plt.rc("text", usetex=True) -plt.rc("font", family="serif", serif="Times") -cmap_med = [ - "#F15A60", - "#7AC36A", - "#5A9BD4", - "#FAA75B", - "#9E67AB", - "#CE7058", - "#D77FB4", - "#737373", -] -cmap = [ - "#EE2E2F", - "#008C48", - "#185AA9", - "#F47D23", - "#662C91", - "#A21D21", - "#B43894", - "#010202", -] -dashseq = [ - (None, None), - [10, 5], - [10, 4, 3, 4], - [3, 3], - [10, 4, 3, 4, 3, 4], - [3, 3], - [3, 3], -] -markertype = ["s", "d", "o", "p", "h"] - -# ======================================================================== -# -# Function definitions -# -# ======================================================================== -def div0(a, b): - """ Ignore division by 0, just replace it by 0, - - From: http://stackoverflow.com/questions/26248654/numpy-return-0-with-divide-by-zero - e.g. div0( [-1, 0, 1], 0 ) -> [0, 0, 0] - """ - with np.errstate(divide="ignore", invalid="ignore"): - c = np.true_divide(a, b) - c[~np.isfinite(c)] = 0 # -inf inf NaN - return c - - -# ======================================================================== -def abs2(x): - """This is equivalent to np.abs(x)**2 or x*np.conj(x) - - To make it faster, add this right before the function definition - import numba - @numba.vectorize([numba.float64(numba.complex128),numba.float32(numba.complex64)]) - """ - return x.real ** 2 + x.imag ** 2 - - -# ======================================================================== -# -# Main -# -# ======================================================================== - -# Timer -start = time.time() - -# ======================================================================== -# 1. velocity fluctuations generated on a 512^3 grid in wavenumber space - -# Dimension of the large cube -N = 512 -halfN = int(round(0.5 * N)) -xs = 0 -xe = 2.0 * np.pi -L = xe - xs -dx = L / N - -# Only work if N and args.N are even -if not ((args.N % 2 == 0) and N % 2 == 0): - print("N or args.N is not even. Exiting") - sys.exit(1) - -# Get cell centered values and meshed grid -x = np.linspace(xs, xe, N + 1) -xc = (x[1:] + x[:-1]) / 2 # get cell center coordinates -X, Y, Z = np.meshgrid(xc, xc, xc, indexing="ij") - -# Get the wave numbers and associated quantities -k = np.concatenate((np.arange(halfN), np.arange(-halfN, 0, 1)), axis=0) -khalf = np.arange(halfN + 1) -k1, k2, k3 = np.meshgrid(k, k, khalf, indexing="ij") -kmag = np.sqrt(k1 ** 2 + k2 ** 2 + k3 ** 2) -k12 = np.sqrt(k1 ** 2 + k2 ** 2) -k1k12 = div0(k1, k12) -k2k12 = div0(k2, k12) -k3kmag = div0(k3, kmag) -k12kmag = div0(k12, kmag) - -# Generate data - -# # Toy Fourier data corresponding to uo = cos(X) * cos(2*Y) * cos(3*Z) -# uo = np.cos(X) * np.cos(2*Y) * np.cos(3*Z) -# uf = np.fft.rfftn(uo) -# vf = np.copy(uf) -# wf = np.copy(uf) - -# Energy spectrum -Ek = ( - 16.0 - * np.sqrt(2.0 / np.pi) - * (kmag ** 4) - / (args.k0 ** 5) - * np.exp(-2.0 * (kmag ** 2) / (args.k0 ** 2)) -) - -# Draw random numbers -np.random.seed(args.seed) -phi1 = np.random.uniform(0, 2 * np.pi, np.shape(kmag)) -phi2 = np.random.uniform(0, 2 * np.pi, np.shape(kmag)) -phi3 = np.random.uniform(0, 2 * np.pi, np.shape(kmag)) - -# the random quantities -prefix = np.sqrt(2.0 * div0(Ek, 4.0 * np.pi * (kmag ** 2))) -a = prefix * np.exp(1j * phi1) * np.cos(phi3) -b = prefix * np.exp(1j * phi2) * np.sin(phi3) - -# the random velocities -uf = k2k12 * a + k1k12 * k3kmag * b -vf = k2k12 * k3kmag * b - k1k12 * a -wf = -k12kmag * b - -# Impose the 3D spherical symmetry (to ensure we have a real signal) -# equiv: uf[-l,-m,0] = np.conj(uf[ l, m,0]) for l=0..N/2 and m=0..N/2 -uf[N:halfN:-1, N:halfN:-1, 0] = np.conj(uf[1:halfN, 1:halfN, 0]) -# symmetry on first column -uf[N:halfN:-1, 0, 0] = np.conj(uf[1:halfN, 0, 0]) -# symmetry on first row -uf[0, N:halfN:-1, 0] = np.conj(uf[0, 1:halfN, 0]) -# symmetry about the (halfN,halfN) element -uf[halfN - 1 : 0 : -1, N : halfN - 1 : -1, 0] = np.conj( - uf[halfN + 1 : N, 1 : halfN + 1, 0] -) - -vf[N:halfN:-1, N:halfN:-1, 0] = np.conj(vf[1:halfN, 1:halfN, 0]) -vf[halfN - 1 : 0 : -1, N : halfN - 1 : -1, 0] = np.conj( - vf[halfN + 1 : N, 1 : halfN + 1, 0] -) -vf[N:halfN:-1, 0, 0] = np.conj(vf[1:halfN, 0, 0]) -vf[0, N:halfN:-1, 0] = np.conj(vf[0, 1:halfN, 0]) - -wf[N:halfN:-1, N:halfN:-1, 0] = np.conj(wf[1:halfN, 1:halfN, 0]) -wf[halfN - 1 : 0 : -1, N : halfN - 1 : -1, 0] = np.conj( - wf[halfN + 1 : N, 1 : halfN + 1, 0] -) -wf[N:halfN:-1, 0, 0] = np.conj(wf[1:halfN, 0, 0]) -wf[0, N:halfN:-1, 0] = np.conj(wf[0, 1:halfN, 0]) - -# Normalize. Because we are generating the data in wavenumber space, -# we have to multiply by N**3 because in the definition of the numpy -# ifftn there is a 1/N**n. -uf = uf * N ** 3 -vf = vf * N ** 3 -wf = wf * N ** 3 - -# # Quick check on energy content (make sure you add both the current -# # contribution and the one we are neglecting because we are assuming -# # real input data) -# print('Energy = int E(k) dk = 0.5 * int (uf**2 + vf**2 wf**2) dk1 dk2 dk3 = {0:.10f} ~= 3/2'.format( -# (np.sum(abs2(uf ) + abs2(vf ) + abs2(wf )) + -# np.sum(abs2(uf[:,:,1:-1]) + abs2(vf[:,:,1:-1]) + abs2(wf[:,:,1:-1]))) -# * 0.5 / N**6)) - -# if plotting, save the original field (before filtering) -if args.plot: - uo = np.fft.irfftn(uf) - Eko = ( - 16.0 - * np.sqrt(2.0 / np.pi) - * (khalf ** 4) - / (args.k0 ** 5) - * np.exp(-2.0 * (khalf ** 2) / (args.k0 ** 2)) - ) - - # Get the spectrum from 3D velocity field - kbins = np.arange(1, halfN + 1) - Nbins = len(kbins) - whichbin = np.digitize(kmag.flat, kbins) - ncount = np.bincount(whichbin) - - KI = (abs2(uf) + abs2(vf) + abs2(wf)) * 0.5 / N ** 6 - KI[:, :, 1:-1] += ( - (abs2(uf[:, :, 1:-1]) + abs2(vf[:, :, 1:-1]) + abs2(wf[:, :, 1:-1])) - * 0.5 - / N ** 6 - ) - - Eku = np.zeros(len(ncount) - 1) - for n in range(1, len(ncount)): - Eku[n - 1] = np.sum(KI.flat[whichbin == n]) - - ku = 0.5 * (kbins[0 : Nbins - 1] + kbins[1:Nbins]) + 1 - Eku = Eku[1:Nbins] - - -# ======================================================================== -# 2. Coefficients associated to wavenumbers that cannot be represented -# on the desired grid are set to 0 (sharp wavenumber cutoff) -kmagc = 0.5 * args.N -uf[kmag > kmagc] = 0.0 -vf[kmag > kmagc] = 0.0 -wf[kmag > kmagc] = 0.0 - - -# ======================================================================== -# 3. inverse Fourier transform of the velocity fluctuations (512^3 grid) -u = np.fft.irfftn(uf, s=(N, N, N)) -v = np.fft.irfftn(vf, s=(N, N, N)) -w = np.fft.irfftn(wf, s=(N, N, N)) - -# Another energy content check -print( - "Energy = 1/V * int E(x,y,z) dV = 0.5/V * int (u**2 + v**2 + w**2) dx dy dz = {0:.10f} ~= 3/2".format( - np.sum(u ** 2 + v ** 2 + w ** 2) * 0.5 * (dx / L) ** 3 - ) -) - -# # Enstrophy check -# _, dudy, dudz = np.gradient(u, dx) -# dvdx, _, dvdz = np.gradient(v, dx) -# dwdx, dwdy, _ = np.gradient(w, dx) -# wx = dwdy-dvdz -# wy = dudz-dwdx -# wz = dvdx-dudy -# lambda0 = 2.0/args.k0 -# print('Enstrophy = 0.5/V * int (wx**2 + wy**2 + wz**2) dx dy dz= -# {0:.10f} ~= '.format(np.sum(wx**2+wy**2+wz**2) * 0.5 * (dx/L)**3 * -# lambda0**2)) - -# ======================================================================== -# 4. velocity fluctuations re-sampled on the desired grid (N^3) -xr = np.linspace(xs, xe, args.N + 1) -xrc = (xr[1:] + xr[:-1]) / 2 -Xr, Yr, Zr = np.meshgrid(xrc, xrc, xrc, indexing="ij") - -Xr = Xr.reshape(-1, order="F") -Yr = Yr.reshape(-1, order="F") -Zr = Zr.reshape(-1, order="F") - -ur = spi.interpn((xc, xc, xc), u, (Xr, Yr, Zr), method="linear") -vr = spi.interpn((xc, xc, xc), v, (Xr, Yr, Zr), method="linear") -wr = spi.interpn((xc, xc, xc), w, (Xr, Yr, Zr), method="linear") - - -# ======================================================================== -# Save the data in Fortran ordering -fname = "hit_ic_{0:d}_{1:d}.dat".format(int(args.k0), args.N) -data = np.vstack((Xr, Yr, Zr, ur, vr, wr)).T -np.savetxt(fname, data, fmt="%.18e", delimiter=",", header="x, y, z, u, v, w") - - -# ======================================================================== -# plot (only u fluctuations) -if args.plot: - import matplotlib as mpl - - mpl.use("Agg") - import matplotlib.pyplot as plt - - datmin = u.min() - datmax = u.max() - # print("min/max u =",datmin,datmax) - - # Original data - # transpose and origin change bc I used meshgrid with ij and not xy - fig, ax = plt.subplots(nrows=3, ncols=3, figsize=(14, 14)) - ax[0, 0].imshow( - uo[:, :, 0].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[0, 0].set_title("Original data (x,y)") - ax[0, 1].imshow( - uo[:, 0, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[0, 1].set_title("Original data (x,z)") - ax[0, 2].imshow( - uo[0, :, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[0, 2].set_title("Original data (y,z)") - - # Filtered original data - ax[1, 0].imshow( - u[:, :, 0].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[1, 0].set_title("Filtered original data (x,y)") - ax[1, 1].imshow( - u[:, 0, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[1, 1].set_title("Filtered original data (x,z)") - ax[1, 2].imshow( - u[0, :, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[1, 2].set_title("Filtered original data (y,z)") - - # Downsampled filtered data - ur = ur.reshape(args.N, args.N, args.N, order="F") - ax[2, 0].imshow( - ur[:, :, 0].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[2, 0].set_title("Downsampled data (x,y)") - ax[2, 1].imshow( - ur[:, 0, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[2, 1].set_title("Downsampled data (x,z)") - ax[2, 2].imshow( - ur[0, :, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[2, 2].set_title("Downsampled data (y,z)") - - plt.savefig("hit_ic_u_{0:d}_{1:d}.png".format(int(args.k0), args.N), format="png") - - # Fourier coefficients of original data - fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(14, 8)) - ax[0, 0].imshow(np.real(uf[:, :, 0].T), origin="lower", cmap="RdBu_r") - ax[0, 0].set_title("Real Fourier coefficients (x,y)") - ax[0, 1].imshow(np.real(uf[:, 0, :].T), origin="lower", cmap="RdBu_r") - ax[0, 1].set_title("Real Fourier coefficients (x,z)") - ax[0, 2].imshow(np.real(uf[0, :, :].T), origin="lower", cmap="RdBu_r") - ax[0, 2].set_title("Real Fourier coefficients (y,z)") - ax[1, 0].imshow(np.imag(uf[:, :, 0].T), origin="lower", cmap="RdBu_r") - ax[1, 0].set_title("Imag Fourier coefficients (x,y)") - ax[1, 1].imshow(np.imag(uf[:, 0, :].T), origin="lower", cmap="RdBu_r") - ax[1, 1].set_title("Imag Fourier coefficients (x,z)") - ax[1, 2].imshow(np.imag(uf[0, :, :].T), origin="lower", cmap="RdBu_r") - ax[1, 2].set_title("Imag Fourier coefficients (y,z)") - plt.savefig("hit_ic_uf_{0:d}_{1:d}.png".format(int(args.k0), args.N), format="png") - - # Spectrum - plt.figure(20) - ax = plt.gca() - p = plt.loglog(khalf, Eko, color=cmap[-1], lw=2) - p[0].set_dashes(dashseq[0]) - p = plt.loglog(ku, Eku, color=cmap[0], lw=2) - p[0].set_dashes(dashseq[1]) - plt.ylim([1e-16, 10]) - plt.xlabel(r"$k$", fontsize=22, fontweight="bold") - plt.ylabel(r"$E(k)$", fontsize=22, fontweight="bold") - plt.setp(ax.get_xmajorticklabels(), fontsize=18, fontweight="bold") - plt.setp(ax.get_ymajorticklabels(), fontsize=18, fontweight="bold") - plt.savefig( - "hit_ic_spectrum_{0:d}_{1:d}.png".format(int(args.k0), args.N), format="png" - ) - -# output timer -end = time.time() - start -print("Elapsed time " + str(timedelta(seconds=end)) + " (or {0:f} seconds)".format(end)) diff --git a/Exec/RegTests/TurbInflow/TurbFileHIT/input b/Exec/RegTests/TurbInflow/TurbFileHIT/input deleted file mode 100644 index 5a59a576..00000000 --- a/Exec/RegTests/TurbInflow/TurbFileHIT/input +++ /dev/null @@ -1,4 +0,0 @@ -hit_file = hit_ic_2_64.dat -input_ncell = 64 -urms0 = 1.0 -TurbFile = TurbTEST diff --git a/Exec/RegTests/TurbInflow/TurbFileHIT/main.cpp b/Exec/RegTests/TurbInflow/TurbFileHIT/main.cpp deleted file mode 100644 index 2d21413e..00000000 --- a/Exec/RegTests/TurbInflow/TurbFileHIT/main.cpp +++ /dev/null @@ -1,262 +0,0 @@ -#include - -#include "AMReX_ParmParse.H" -#include -#include -#include - -using namespace amrex; - -static void -Extend(FArrayBox& xfab, FArrayBox& vfab, const Box& domain) -{ - Box tbx = vfab.box(); - - tbx.setBig(0, domain.bigEnd(0) + 3); - - const int ygrow = AMREX_SPACEDIM == 3 ? 3 : 1; - - tbx.setBig(1, domain.bigEnd(1) + ygrow); - - xfab.resize(tbx, 1); - - Box orig_box = vfab.box(); - vfab.shift(0, 1); - vfab.shift(1, 1); - xfab.copy(vfab); // (0,0) - - vfab.shift(0, domain.length(0)); // (1,0) - xfab.copy(vfab); - vfab.shift(1, domain.length(1)); // (1,1) - xfab.copy(vfab); - vfab.shift(0, -domain.length(0)); // (0,1) - xfab.copy(vfab); - vfab.shift(0, -domain.length(0)); // (-1,1) - xfab.copy(vfab); - vfab.shift(1, -domain.length(1)); // (-1,0) - xfab.copy(vfab); - vfab.shift(1, -domain.length(1)); // (-1,-1) - xfab.copy(vfab); - vfab.shift(0, domain.length(0)); // (0,-1) - xfab.copy(vfab); - vfab.shift(0, domain.length(0)); // (1,-1) - xfab.copy(vfab); - vfab.shift(0, -domain.length(0) - 1); - vfab.shift(1, domain.length(1) - 1); - - if (vfab.box() != orig_box) - Abort("Oops, something bad happened"); -} - -void -readHIT(HITData* a_data) -{ - ParmParse pp; - std::string hit_file("IC"); - pp.query("hit_file", hit_file); - - pp.query("input_ncell", a_data->input_ncell); - if (a_data->input_ncell == 0) { - amrex::Abort("for HIT, the input ncell cannot be 0 !"); - } - - int binfmt = 0; // Default is ASCII format - pp.query("input_binaryformat", binfmt); - pp.query("urms0", a_data->urms0); - - amrex::Real lambda0 = 0.5; - amrex::Real tau = lambda0 / a_data->urms0; - - // Output initial conditions - std::ofstream ofs("initialConditions.txt", std::ofstream::out); - amrex::Print(ofs) << "lambda0, urms0, tau \n"; - amrex::Print(ofs).SetPrecision(17) - << lambda0 << "," << a_data->urms0 << "," << tau << std::endl; - ofs.close(); - - // Read initial velocity field - const size_t nx = a_data->input_ncell; - const size_t ny = a_data->input_ncell; - const size_t nz = a_data->input_ncell; - amrex::Vector data( - nx * ny * nz * 6); /* this needs to be double */ - if (binfmt) { - read_binary(hit_file, nx, ny, nz, 6, data); - } else { - read_csv(hit_file, nx, ny, nz, data); - } - - // Extract position and velocities - amrex::Vector xinput; - amrex::Vector uinput; - amrex::Vector vinput; - amrex::Vector winput; - amrex::Vector xdiff; - amrex::Vector xarray; - - xinput.resize(nx * ny * nz); - uinput.resize(nx * ny * nz); - vinput.resize(nx * ny * nz); - winput.resize(nx * ny * nz); - - for (long i = 0; i < xinput.size(); i++) { - xinput[i] = data[0 + i * 6]; - uinput[i] = data[3 + i * 6] * a_data->urms0 / a_data->uin_norm; - vinput[i] = data[4 + i * 6] * a_data->urms0 / a_data->uin_norm; - winput[i] = data[5 + i * 6] * a_data->urms0 / a_data->uin_norm; - } - - // Get the xarray table and the differences. - xarray.resize(nx); - for (long i = 0; i < xarray.size(); i++) { - xarray[i] = xinput[i]; - } - xdiff.resize(nx); - std::adjacent_difference(xarray.begin(), xarray.end(), xdiff.begin()); - xdiff[0] = xdiff[1]; - - // Make sure the search array is increasing - if (not std::is_sorted(xarray.begin(), xarray.end())) { - amrex::Abort("Error: non ascending x-coordinate array."); - } - - // Pass data to the prob_parm - a_data->Linput = xarray[nx - 1] + 0.5 * xdiff[nx - 1]; - - a_data->d_xarray = - (amrex::Real*)amrex::The_Arena()->alloc(nx * sizeof(amrex::Real)); - a_data->d_xdiff = - (amrex::Real*)amrex::The_Arena()->alloc(nx * sizeof(amrex::Real)); - a_data->d_uinput = - (amrex::Real*)amrex::The_Arena()->alloc(nx * ny * nz * sizeof(amrex::Real)); - a_data->d_vinput = - (amrex::Real*)amrex::The_Arena()->alloc(nx * ny * nz * sizeof(amrex::Real)); - a_data->d_winput = - (amrex::Real*)amrex::The_Arena()->alloc(nx * ny * nz * sizeof(amrex::Real)); - - for (int i = 0; i < nx; i++) { - a_data->d_xarray[i] = xarray[i]; - a_data->d_xdiff[i] = xdiff[i]; - } - for (int i = 0; i < nx * ny * nz; i++) { - a_data->d_uinput[i] = uinput[i]; - a_data->d_vinput[i] = vinput[i]; - a_data->d_winput[i] = winput[i]; - } -} - -int -main(int argc, char* argv[]) -{ - Initialize(argc, argv); - { - ParmParse pp; - - HITData data; - - readHIT(&data); - - int ncell = data.input_ncell; - Real xlo = data.d_xarray[0]; - Real xhi = data.d_xarray[ncell - 1]; - - Box box_turb( - IntVect(AMREX_D_DECL(0, 0, 0)), - IntVect(AMREX_D_DECL(ncell - 1, ncell - 1, ncell - 1))); - RealBox rb_turb( - {AMREX_D_DECL(xlo, xlo, xlo)}, {AMREX_D_DECL(xhi, xhi, xhi)}); - int coord_turb(0); - Array per_turb = {AMREX_D_DECL(1, 1, 1)}; - Geometry geom_turb(box_turb, rb_turb, coord_turb, per_turb); - const Real* dx_turb = geom_turb.CellSize(); - - // Fill the velocity FAB with HIT - FArrayBox vel_turb(box_turb, AMREX_SPACEDIM); - Array4 const& fab = vel_turb.array(); - auto geomdata = geom_turb.data(); - AMREX_PARALLEL_FOR_3D( - box_turb, i, j, k, { fillVelFab(i, j, k, fab, geomdata, data); }); - - std::ofstream fabOut; - fabOut.open("Turb_D", std::ios::out | std::ios::trunc); - vel_turb.writeOn(fabOut); - // Write turbulence file for TurbInflow - std::string TurbDir("Turb"); - pp.query("TurbFile", TurbDir); - - if (ParallelDescriptor::IOProcessor()) - if (!UtilCreateDirectory(TurbDir, 0755)) - CreateDirectoryFailed(TurbDir); - - std::string Hdr = TurbDir; - Hdr += "/HDR"; - std::string Dat = TurbDir; - Dat += "/DAT"; - - std::ofstream ifsd, ifsh; - - ifsh.open(Hdr.c_str(), std::ios::out | std::ios::trunc); - if (!ifsh.good()) - FileOpenFailed(Hdr); - - ifsd.open(Dat.c_str(), std::ios::out | std::ios::trunc); - if (!ifsd.good()) - FileOpenFailed(Dat); - - // - // Write the first part of the Turb header. - // Note that this is solely for periodic style inflow files. - // - Box box_turb_io(box_turb); - box_turb_io.setBig(0, box_turb.bigEnd(0) + 3); - box_turb_io.setBig(1, box_turb.bigEnd(1) + 3); - box_turb_io.setBig(2, box_turb.bigEnd(2) + 1); - - ifsh << box_turb_io.length(0) << ' ' << box_turb_io.length(1) << ' ' - << box_turb_io.length(2) << '\n'; - - ifsh << rb_turb.length(0) + 2 * dx_turb[0] << ' ' - << rb_turb.length(1) + 2 * dx_turb[1] << ' ' << rb_turb.length(2) - << '\n'; - - ifsh << per_turb[0] << ' ' << per_turb[1] << ' ' << per_turb[2] << '\n'; - - // Dump field as a "turbulence file" - IntVect sm = box_turb.smallEnd(); - IntVect bg = box_turb.bigEnd(); - int dir = AMREX_SPACEDIM - 1; - FArrayBox xfab, TMP; - // - // We work on one cell wide Z-planes. - // We first do the lo AMREX_SPACEDIM plane. - // And then all the other planes in xhi -> xlo order. - // - for (int d = 0; d < AMREX_SPACEDIM; ++d) { - bg[dir] = sm[dir]; - { - Box bx(sm, bg); - TMP.resize(bx, 1); - TMP.copy(vel_turb, bx, d, bx, 0, 1); - Extend(xfab, TMP, box_turb); - ifsh << ifsd.tellp() << std::endl; - xfab.writeOn(ifsd); - } - for (int i = box_turb.bigEnd(dir); i >= box_turb.smallEnd(dir); i--) { - sm[dir] = i; - bg[dir] = i; - Box bx(sm, bg); - TMP.resize(bx, 1); - TMP.copy(vel_turb, bx, d, bx, 0, 1); - Extend(xfab, TMP, box_turb); - ifsh << ifsd.tellp() << std::endl; - xfab.writeOn(ifsd); - } - } - amrex::The_Arena()->free(data.d_xarray); - amrex::The_Arena()->free(data.d_xdiff); - amrex::The_Arena()->free(data.d_uinput); - amrex::The_Arena()->free(data.d_vinput); - amrex::The_Arena()->free(data.d_winput); - } - Finalize(); -} diff --git a/Exec/RegTests/TurbInflow/TurbFileHIT/main_K.H b/Exec/RegTests/TurbInflow/TurbFileHIT/main_K.H deleted file mode 100644 index 3a8ff369..00000000 --- a/Exec/RegTests/TurbInflow/TurbFileHIT/main_K.H +++ /dev/null @@ -1,96 +0,0 @@ -#ifndef MAIN_K_H -#define MAIN_K_H - -#include -#include -#include - -#include -#include - -AMREX_GPU_DEVICE -AMREX_FORCE_INLINE -void -fillVelFab( - int i, - int j, - int k, - amrex::Array4 const& vfab, - amrex::GeometryData const& geomdata, - HITData const& a_data) -{ - - const amrex::Real* prob_lo = geomdata.ProbLo(); - const amrex::Real* prob_hi = geomdata.ProbHi(); - const amrex::Real* dx = geomdata.CellSize(); - - amrex::Real x[3] = { - prob_lo[0] + static_cast(i + 0.5) * dx[0], - prob_lo[1] + static_cast(j + 0.5) * dx[1], - prob_lo[2] + static_cast(k + 0.5) * dx[2]}; - - // Fill in the velocities - amrex::Real uinterp[3] = {0.0}; - - // Interpolation factors - amrex::Real mod[3] = {0.0}; - int idx[3] = {0}; - int idxp1[3] = {0}; - amrex::Real slp[3] = {0.0}; - for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { - mod[idim] = std::fmod(x[idim], a_data.Linput); - locate(a_data.d_xarray, a_data.input_ncell, mod[idim], idx[idim]); - idxp1[idim] = (idx[idim] + 1) % a_data.input_ncell; - slp[idim] = - (mod[idim] - a_data.d_xarray[idx[idim]]) / a_data.d_xdiff[idx[idim]]; - } - - int inSize = a_data.input_ncell; - - const amrex::Real f0 = (1 - slp[0]) * (1 - slp[1]) * (1 - slp[2]); - const amrex::Real f1 = slp[0] * (1 - slp[1]) * (1 - slp[2]); - const amrex::Real f2 = (1 - slp[0]) * slp[1] * (1 - slp[2]); - const amrex::Real f3 = (1 - slp[0]) * (1 - slp[1]) * slp[2]; - const amrex::Real f4 = slp[0] * (1 - slp[1]) * slp[2]; - const amrex::Real f5 = (1 - slp[0]) * slp[1] * slp[2]; - const amrex::Real f6 = slp[0] * slp[1] * (1 - slp[2]); - const amrex::Real f7 = slp[0] * slp[1] * slp[2]; - - uinterp[0] = - a_data.d_uinput[idx[0] + inSize * (idx[1] + inSize * idx[2])] * f0 + - a_data.d_uinput[idxp1[0] + inSize * (idx[1] + inSize * idx[2])] * f1 + - a_data.d_uinput[idx[0] + inSize * (idxp1[1] + inSize * idx[2])] * f2 + - a_data.d_uinput[idx[0] + inSize * (idx[1] + inSize * idxp1[2])] * f3 + - a_data.d_uinput[idxp1[0] + inSize * (idx[1] + inSize * idxp1[2])] * f4 + - a_data.d_uinput[idx[0] + inSize * (idxp1[1] + inSize * idxp1[2])] * f5 + - a_data.d_uinput[idxp1[0] + inSize * (idxp1[1] + inSize * idx[2])] * f6 + - a_data.d_uinput[idxp1[0] + inSize * (idxp1[1] + inSize * idxp1[2])] * f7; - - uinterp[1] = - a_data.d_vinput[idx[0] + inSize * (idx[1] + inSize * idx[2])] * f0 + - a_data.d_vinput[idxp1[0] + inSize * (idx[1] + inSize * idx[2])] * f1 + - a_data.d_vinput[idx[0] + inSize * (idxp1[1] + inSize * idx[2])] * f2 + - a_data.d_vinput[idx[0] + inSize * (idx[1] + inSize * idxp1[2])] * f3 + - a_data.d_vinput[idxp1[0] + inSize * (idx[1] + inSize * idxp1[2])] * f4 + - a_data.d_vinput[idx[0] + inSize * (idxp1[1] + inSize * idxp1[2])] * f5 + - a_data.d_vinput[idxp1[0] + inSize * (idxp1[1] + inSize * idx[2])] * f6 + - a_data.d_vinput[idxp1[0] + inSize * (idxp1[1] + inSize * idxp1[2])] * f7; - - uinterp[2] = - a_data.d_winput[idx[0] + inSize * (idx[1] + inSize * idx[2])] * f0 + - a_data.d_winput[idxp1[0] + inSize * (idx[1] + inSize * idx[2])] * f1 + - a_data.d_winput[idx[0] + inSize * (idxp1[1] + inSize * idx[2])] * f2 + - a_data.d_winput[idx[0] + inSize * (idx[1] + inSize * idxp1[2])] * f3 + - a_data.d_winput[idxp1[0] + inSize * (idx[1] + inSize * idxp1[2])] * f4 + - a_data.d_winput[idx[0] + inSize * (idxp1[1] + inSize * idxp1[2])] * f5 + - a_data.d_winput[idxp1[0] + inSize * (idxp1[1] + inSize * idx[2])] * f6 + - a_data.d_winput[idxp1[0] + inSize * (idxp1[1] + inSize * idxp1[2])] * f7; - - // - // Fill Velocity - // - vfab(i, j, k, 0) = uinterp[0]; - vfab(i, j, k, 1) = uinterp[1]; - vfab(i, j, k, 2) = uinterp[2]; -} -#endif diff --git a/Submodules/PelePhysics b/Submodules/PelePhysics index c0f30b91..80f2cfa0 160000 --- a/Submodules/PelePhysics +++ b/Submodules/PelePhysics @@ -1 +1 @@ -Subproject commit c0f30b91b1d21f98c2b13bb1c37f5c86f29f7a0a +Subproject commit 80f2cfa0edd2e928a93bb9dff45e68d32d8699f7 diff --git a/Utils/CloneDeps.sh b/Utils/CloneDeps.sh index 170df4cd..daa28d2f 100755 --- a/Utils/CloneDeps.sh +++ b/Utils/CloneDeps.sh @@ -1,4 +1,4 @@ #!/usr/bin/env bash echo "Getting PeleLMeX dependencies - tests ... " -git submodule update --init --recursive +git submodule update --init --recursive --recommend-shallow