From 80f2cfa0edd2e928a93bb9dff45e68d32d8699f7 Mon Sep 17 00:00:00 2001 From: Bruce Perry <53018946+baperry2@users.noreply.github.com> Date: Mon, 15 Apr 2024 14:52:47 -0600 Subject: [PATCH] Move turb inflow helpers from PeleLMeX and add documentation (#489) * Bring turbulent inflow generation stuff from PeleLMeX * some progress on documentation * more progress on turbinflow documentation * add option for size of wavenumber space in turbulence generation * make multiple turbinflows additive rather than replacing each other * fix typo and math format --- Docs/sphinx/Support.rst | 132 +++++++ Docs/sphinx/Utility.rst | 37 +- Docs/sphinx/index.rst | 5 +- Source/Utility/TurbInflow/turbinflow.cpp | 4 +- Support/TurbFileHIT/GNUmakefile | 26 ++ Support/TurbFileHIT/HITData.H | 21 + Support/TurbFileHIT/Make.package | 2 + Support/TurbFileHIT/README.md | 14 + Support/TurbFileHIT/Utilities.H | 72 ++++ Support/TurbFileHIT/Utilities.cpp | 87 +++++ Support/TurbFileHIT/gen_hit_ic.py | 467 +++++++++++++++++++++++ Support/TurbFileHIT/input | 4 + Support/TurbFileHIT/main.cpp | 262 +++++++++++++ Support/TurbFileHIT/main_K.H | 96 +++++ 14 files changed, 1223 insertions(+), 6 deletions(-) create mode 100644 Docs/sphinx/Support.rst create mode 100644 Support/TurbFileHIT/GNUmakefile create mode 100644 Support/TurbFileHIT/HITData.H create mode 100644 Support/TurbFileHIT/Make.package create mode 100644 Support/TurbFileHIT/README.md create mode 100644 Support/TurbFileHIT/Utilities.H create mode 100644 Support/TurbFileHIT/Utilities.cpp create mode 100755 Support/TurbFileHIT/gen_hit_ic.py create mode 100644 Support/TurbFileHIT/input create mode 100644 Support/TurbFileHIT/main.cpp create mode 100644 Support/TurbFileHIT/main_K.H diff --git a/Docs/sphinx/Support.rst b/Docs/sphinx/Support.rst new file mode 100644 index 000000000..ce740705b --- /dev/null +++ b/Docs/sphinx/Support.rst @@ -0,0 +1,132 @@ +.. highlight:: rst + +.. _sec:Support: + +******* +Support +******* + +The main purpose of PelePhysics is to contain a library of routines for shared use by the other Pele codes, +which are contained in the ``Source`` directory. In the ``Support`` directory, PelePhysics also contains +several stand-alone tools and scripts that are used to generate data and source code that support these purposes. +These tools include: + +* ``ceptr``: Chemistry Evaluation for Pele Through Recasting, machine generation of mechanism code +* ``TurbFileHIT``: Isotropic turbulence initial conditions and inflows +* ``MechanismPAH``: Append PAH module to chemical mechanism data +* ``liquidProp``: Fits to NIST data for condensed phase species + +ceptr +===== + +CEPTR is used to generate the C++ mechanism files from Cantera-format mechanism data for all the chemical +mechanisms in the ``Mechanisms/`` directory. For a full description of CEPTR and how to use it, see the +:ref:`CEPTR subsection ` of the Chemistry section of this documentation. + +.. _sec_turbfile: + +TurbFileHIT +=========== + +This support code contains two separate pieces: a python script that generates a synthetic 3D isotropic turbulence velocity field and +an AMReX-based C++ code that converts that data into a format that can be used by the PelePhysics :ref:`TurbInflow Utility ` +to provide turbulent inflow boundary conditions in either PeleC or PeleLMeX. + +Generating an HIT File +~~~~~~~~~~~~~~~~~~~~~~ + +The ``gen_hit_ic.py`` python script creates a synthetic isotropic turbulence velocity field based on the methodology described +in the appendix of `Johnsen et al. (2010) J. Comp Phys. `_. The resulting velocity +field is stored in files with names following the pattern ``hit_ic_k0_N.dat``, where ``k0`` is the most +energetic wave number and ``N`` is the size of the grid. The initial +condition for a grid of size :math:`N^3` is generated as follows: + +1. velocity fluctuations generated on a :math:`N_k^3` grid in wavenumber space (default :math:`N_k = 512`) +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 (:math:`N_k^3` grid) +4. velocity fluctuations resampled on the desired grid (:math:`N^3`) + +This script accepts the following options: :: + + ./gen_hit_ic.py --help + usage: gen_hit_ic.py [-h] [-k0 K0] [-N N] [-Nk NK] [-s SEED] [-p] + + Generate the velocity fluctuations for the HIT IC + + options: + -h, --help show this help message and exit + -k0 K0 Wave number containing highest energy + -N N Resolution + -Nk NK Resolution in wavenumber space for intermediate step + -s SEED, --seed SEED Random number generator seed + -p, --plot Save a plot of the x-velocity + +Generating an initial condition file is as easy as: :: + + ./gen_hit_ic.py -N 16 + +.. Note:: + + This script utilizes standard python libraries such as SciPy and Matplotlib. If you'd like to ensure that you are using appropriate + versions of the dependencies, follow the :ref:`CEPTR instructions ` for setting up ``poetry`` to manage + dependencies, then run the script using ``poetry run -C ../ceptr/ python gen_hit_ic.py -N 16``. + +Generating Inflow Files +~~~~~~~~~~~~~~~~~~~~~~~ + +The ``TurbFileHIT`` directory also contains C++ source code for a small program that reads the ``.dat`` files created using the python script +above and saves them as AMReX data structures that can be read by the TurbInflow utility, which uses Taylor's hypothesis to provide +temporally varying boundary data by marching through the 3rd dimension of the HIT files at a constant velocity. + +Build options can be specified in the ``GNUmakefile``. Notably, if you did not recursively clone PelePhysics with the AMReX submodule, you +will need to specify a path to AMReX (``AMREX_HOME``). But in most cases, you should be able to directly build the utility simply by running: :: + + make -j + +Then run the resulting executable using the provided input file: :: + + ./PeleTurb3d.gnu.ex input + +Adjust the ``hit_file`` and ``input_ncell`` parameters in the input file to match the file name and number of cells from the previous step. +The ``urms0`` parameter is a scale factor to rescale the velocity fluctuations. ``TurbFile`` specified the name of the directory where the +output files will be saved. After successful execution, the output directory should contain two files: ``HDR`` and ``DAT``. + + +MechanismPAH +============ + +This functionality takes an existing set of mechanism, transport, and thermodynamic data files in Chemkin format and attempts to add a PAH module to it. It checks if any species are duplicated by comparing the atomic composition and the enthalpy curves. It check if any reactions are duplicated. If any species or reactions are duplicated, these are skipped. Once the new yaml file is created, ignition delay time at 1 and 20 atm is compared between the original and new mechanism to see what impact the PAH module has on the mechanism. Plot files of these values are created for you to decide if the differences are significant. + +Usage +~~~~~ + +You need to have the NumPy, Cantera, and MatPlotLib python modules. In order to run, use the following :: + + python addPAHmech.py --mech origmech.inp --transport origtrans.dat --thermo origthermo.dat --fuelname NC10H22 + +where ``origmech.inp``, ``origthermo.dat``, and ``origtrans.dat`` are the initial mechanism, thermodynamic, +and transport files to have the PAH module amended to. + +Disclaimer +~~~~~~~~~~ + +The resulting IDT should be studied to determine if the new mechanism is now compromised with the addition of the PAH module. This is left up to the user's discretion. This has only been tested on a few mechanisms and might have bugs. + +liquidProp +========== + +This is a python script that reads in an NIST property file for a condensed or saturated phase of a species. +Files for rho, mu, and lambda (thermal conductivity) should be provided in a single directory and +named ``rho.dat``, ``mu.dat``, and ``lambda.dat``. The usage for this script is:: + + $ python propcoeff.py -h + + usage: propcoeff.py [-h] --species NC10H22 [--file_loc FILE_LOC] [--units UNITS] [--vars VARS [VARS ...]] + + options: + -h, --help show this help message and exit + --species NC10H22 Species name + --file_loc FILE_LOC Location of data files. Files should be called rho.dat, mu.dat, and/or lambda.dat + --units UNITS Units, either MKS or CGS + --vars VARS [VARS ...] + Which variables to fit, ex. mu lambda rho diff --git a/Docs/sphinx/Utility.rst b/Docs/sphinx/Utility.rst index d2b7fbc88..5d2e93a1b 100644 --- a/Docs/sphinx/Utility.rst +++ b/Docs/sphinx/Utility.rst @@ -44,16 +44,49 @@ An additional script is provided to allow plotting of the PMF solutions. This sc poetry -C ../../Support/ceptr/ run python plotPMF.py +.. _sec_turbinflow: Turbulent Inflows ================= -Placeholder. PelePhysics supports the capability of the flow solvers to have spatially and temporally varying inflow conditions based on precomputed turbulence data. +PelePhysics supports the capability of the flow solvers to have spatially and temporally varying inflow conditions based on precomputed turbulence data (3 components of velocity fluctuations). +A three-dimensional synthetic istropic turbulence is generated and Taylor's hypothesis is used to convert this data into two-dimensional planar data by moving through the third dimension at fixed velocity (the TurbInflow capability is currently not supported for 2D simulations). To reduce memory requirements, a fixed number of planes in the third dimension are read at a time and then replaced when exhausted. This number of planes can be specified at runtime to balance I/O and memory requirements. Multiple turbulent inflow patches can be applied on any domain boundary or on multiple domain boundaries. If multiple patches overlap on the same face, the data from each overlapping patch are superimposed. + +This PelePhysics utility provides the machinery to load the data files and interpolate in space and time onto boundary patches, which may or may not cover an entire boundary. Additional code within PeleC and PeleLMeX is required to drive this functionality, and the documentation for the relevant code should be consulted to fully understand the necessary steps. Typically, the TurbInflow capability provides veloicity fluctuations, and the mean inflow velocity must be provided through another means. For each code, an example test case for the capability is provided in `Exec/RegTests/TurbInflow`. Note that the turbulence data is always a square of size 2pi and has a fluctuation velocity scale of unity. Inputs are available as part of this utility to rescale the data as needed. If differently shaped inlet patches are required, this must be done by masking undesired parts of the patch on the PeleC or PeleLMeX side of the implementation. Generating a turbulence file ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -A python script is used to generate a synthetic turbulence spectrum. +The relevant data files are generated using the tools in ``Support/TurbFileHIT`` +and usage instructions are available in the :ref:`documentation ` on these tools. + +Input file options +~~~~~~~~~~~~~~~~~~ + +The input file options that drive this utility (and thus are inherited by PeleC and PeleLMeX) are +provided below. :: + + turbinflows=low high # Names of injections (can provide any number) + + turbinflow.low.turb_file = TurbFileHIT/TurbTEST # Path to directory created in previous step + turbinflow.low.dir = 1 # Boundary normal direction (0,1, or 2) for patch + turbinflow.low.side = "low" # Boundary side (low or high) for patch + turbinflow.low.turb_scale_loc = 633.151 # Factor by which to scale the spatial coordinate between the data file and simulation + turbinflow.low.turb_scale_vel = 1.0 # Factor by which to scale the velcoity between the data file and simulation + turbinflow.low.turb_center = 0.005 0.005 # Center point where turbulence patch will be applied + turbinflow.low.turb_conv_vel = 5. # Velocity to move through the 3rd dimension to simulate time evolution + turbinflow.low.turb_nplane = 32 # Number of planes to read and store at a time + turbinflow.low.time_offset = 0.0 # Offset in time for reading through the 3rd dimension + + turbinflow.high.turb_file = TurbFileHIT/TurbTEST # All same as above, but for second injection patch + turbinflow.high.dir = 1 + turbinflow.high.side = "high" + turbinflow.high.turb_scale_loc = 633.151 + turbinflow.high.turb_scale_vel = 1.0 + turbinflow.high.turb_center = 0.005 0.005 + turbinflow.high.turb_conv_vel = 5. + turbinflow.high.turb_nplane = 32 + turbinflow.high.time_offset = 0.0006 Plt File Management =================== diff --git a/Docs/sphinx/index.rst b/Docs/sphinx/index.rst index f81e1c72d..225690454 100644 --- a/Docs/sphinx/index.rst +++ b/Docs/sphinx/index.rst @@ -8,7 +8,7 @@ PelePhysics =========== `PelePhysics` is a repository of physics databases and implementation code for use within the other `Pele` codes. In particular, the choice of chemistry and transport models as well as associated functions and capabilities are managed in `PelePhysics`. `PelePhysics` has an official project `homepage `_, and can be obtained via -`GitHub `_. The documentation pages appearing here are distributed with the code in the ``Docs/sphinx`` folder as "restructured text" files. +`GitHub `_. The documentation pages appearing here are distributed with the code in the ``Docs/sphinx`` folder as "restructured text" files. The html is built automatically with certain pushes to the `PelePhysics` GibHub repository. A local version can also be built as follows :: cd ${PELE_PHYSICS_DIR}/build @@ -19,7 +19,7 @@ point your web browser at the file ``${PELE_PHYSICS_DIR}/build/html/index.html`` .. |a| image:: ./Visualization/PeleSuite.png -.. table:: +.. table:: :align: center +-----+ @@ -40,6 +40,7 @@ point your web browser at the file ``${PELE_PHYSICS_DIR}/build/html/index.html`` Spray.rst Soot.rst Utility.rst + Support.rst Tutorials.rst DeveloperGuide.rst diff --git a/Source/Utility/TurbInflow/turbinflow.cpp b/Source/Utility/TurbInflow/turbinflow.cpp index a944de9f5..266badc7a 100644 --- a/Source/Utility/TurbInflow/turbinflow.cpp +++ b/Source/Utility/TurbInflow/turbinflow.cpp @@ -382,8 +382,8 @@ TurbInflow::fill_turb_plane( ydata[ii] = cy[0] * zdata[ii][0] + cy[1] * zdata[ii][1] + cy[2] * zdata[ii][2]; } - vd(i, j, k, n) = cx[0] * ydata[0] + cx[1] * ydata[1] + cx[2] * ydata[2]; - vd(i, j, k, n) *= velScale; + vd(i, j, k, n) += + velScale * (cx[0] * ydata[0] + cx[1] * ydata[1] + cx[2] * ydata[2]); } } }); diff --git a/Support/TurbFileHIT/GNUmakefile b/Support/TurbFileHIT/GNUmakefile new file mode 100644 index 000000000..1a326427f --- /dev/null +++ b/Support/TurbFileHIT/GNUmakefile @@ -0,0 +1,26 @@ +AMREX_HOME ?= ../../Submodules/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/Support/TurbFileHIT/HITData.H b/Support/TurbFileHIT/HITData.H new file mode 100644 index 000000000..738577829 --- /dev/null +++ b/Support/TurbFileHIT/HITData.H @@ -0,0 +1,21 @@ +#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/Support/TurbFileHIT/Make.package b/Support/TurbFileHIT/Make.package new file mode 100644 index 000000000..04d5df0ea --- /dev/null +++ b/Support/TurbFileHIT/Make.package @@ -0,0 +1,2 @@ +CEXE_headers += HITData.H main_K.H Utilities.H +CEXE_sources += main.cpp Utilities.cpp diff --git a/Support/TurbFileHIT/README.md b/Support/TurbFileHIT/README.md new file mode 100644 index 000000000..6c5a2073c --- /dev/null +++ b/Support/TurbFileHIT/README.md @@ -0,0 +1,14 @@ +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 generate the turbfile (adapt the input file to your needs): +./PeleTurb3d.gnu.ex input diff --git a/Support/TurbFileHIT/Utilities.H b/Support/TurbFileHIT/Utilities.H new file mode 100644 index 000000000..1306fc866 --- /dev/null +++ b/Support/TurbFileHIT/Utilities.H @@ -0,0 +1,72 @@ +#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/Support/TurbFileHIT/Utilities.cpp b/Support/TurbFileHIT/Utilities.cpp new file mode 100644 index 000000000..ea77b0320 --- /dev/null +++ b/Support/TurbFileHIT/Utilities.cpp @@ -0,0 +1,87 @@ +#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/Support/TurbFileHIT/gen_hit_ic.py b/Support/TurbFileHIT/gen_hit_ic.py new file mode 100755 index 000000000..b74c84cc7 --- /dev/null +++ b/Support/TurbFileHIT/gen_hit_ic.py @@ -0,0 +1,467 @@ +#!/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 Nk^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 (Nk^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( + "-Nk", help="Resolution in wavenumber space for intermediate step", type=int, + default=512 +) +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 Nk^3 grid in wavenumber space + +# Dimension of the large cube +N = args.Nk +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 (Nk^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/Support/TurbFileHIT/input b/Support/TurbFileHIT/input new file mode 100644 index 000000000..5a59a5767 --- /dev/null +++ b/Support/TurbFileHIT/input @@ -0,0 +1,4 @@ +hit_file = hit_ic_2_64.dat +input_ncell = 64 +urms0 = 1.0 +TurbFile = TurbTEST diff --git a/Support/TurbFileHIT/main.cpp b/Support/TurbFileHIT/main.cpp new file mode 100644 index 000000000..2d21413e5 --- /dev/null +++ b/Support/TurbFileHIT/main.cpp @@ -0,0 +1,262 @@ +#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/Support/TurbFileHIT/main_K.H b/Support/TurbFileHIT/main_K.H new file mode 100644 index 000000000..3a8ff369e --- /dev/null +++ b/Support/TurbFileHIT/main_K.H @@ -0,0 +1,96 @@ +#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