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

@@ -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 <sec:ceptr>` 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 <sec_turbinflow>`
+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. <http://dx.doi.org/10.1016/j.jcp.2009.10.028>`_. 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 <sec_ceptr_software>` 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

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 <pmf-file> <variable-index>
 
+.. _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 <sec_turbfile>` 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
 ===================

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 <https://amrex-combustion.github.io/>`_, and can be obtained via
-`GitHub <https://github.com/AMReX-Combustion/PelePhysics>`_.  The documentation pages appearing here are distributed with the code in the ``Docs/sphinx`` folder as "restructured text" files.  
+`GitHub <https://github.com/AMReX-Combustion/PelePhysics>`_.  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
 

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]);
       }
     }
   });

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

Support/TurbFileHIT/HITData.H

@@ -0,0 +1,21 @@
+#ifndef HITDATA_H
+#define HITDATA_H
+
+#include <AMReX_REAL.H>
+#include <AMReX_GpuMemory.H>
+
+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

Support/TurbFileHIT/Make.package

@@ -0,0 +1,2 @@
+CEXE_headers += HITData.H main_K.H Utilities.H
+CEXE_sources += main.cpp Utilities.cpp

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

Support/TurbFileHIT/Utilities.H

@@ -0,0 +1,72 @@
+#ifndef _UTILITIES_H
+#define _UTILITIES_H
+
+#include <AMReX_FArrayBox.H>
+
+AMREX_FORCE_INLINE
+std::string
+read_file(std::ifstream& in)
+{
+  return static_cast<std::stringstream const&>(
+           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<amrex::Real>& data);
+
+void read_csv(
+  const std::string& iname,
+  const size_t nx,
+  const size_t ny,
+  const size_t nz,
+  amrex::Vector<amrex::Real>& 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