From 769ae634f5841212fe6ffd8fd19b36e63cecf185 Mon Sep 17 00:00:00 2001 From: Bruce Perry <53018946+baperry2@users.noreply.github.com> Date: Fri, 1 Dec 2023 10:48:41 -0700 Subject: [PATCH] Filter Utility (#445) * Explicit filtering of hydro fluxes and source term. (#112) This commit implements explicit filtering of the hydro fluxes and source term (for both MOL and PPM code paths). This is necessary to perform explicitly filtered Large Eddy Simulations. Contains: - regressions tests - documentation - CMake implementation * Remove legacy Fortran code. (#252) * Use default for trivial constructors and only build HIP in CI for gfx908 (#434) * Use default for trivial constructors. * Only build HIP in CI for gfx908. * Update PelePhysics. * Fixes to header files (#439) * Fixes to header files. * Update PelePhysics submodule. * More fixes. * Initialize xdiff in HIT. * Invert some mfiter loops (#446) Co-authored-by: Jon Rood * Add GPU syncs because MF ParallelFor is non-blocking (#471) * Rename header guards (#506) * Replace _OPENMP with AMREX_USE_OMP to support Ascent. (#519) * Address some warnings (#554) * remove dependence of Filter on Utilities (#710) * move files * docs for filter * use amrex Pi instead of PeleC Pi --------- Co-authored-by: Marc T. Henry de Frahan Co-authored-by: Jon Rood --- Docs/sphinx/Utility.rst | 53 +++- Utility/Filter/Filter.H | 168 ++++++++++++ Utility/Filter/Filter.cpp | 491 ++++++++++++++++++++++++++++++++++++ Utility/Filter/Make.package | 5 + 4 files changed, 714 insertions(+), 3 deletions(-) create mode 100644 Utility/Filter/Filter.H create mode 100644 Utility/Filter/Filter.cpp create mode 100644 Utility/Filter/Make.package diff --git a/Docs/sphinx/Utility.rst b/Docs/sphinx/Utility.rst index e26e811ec..d2b7fbc88 100644 --- a/Docs/sphinx/Utility.rst +++ b/Docs/sphinx/Utility.rst @@ -14,7 +14,7 @@ In addition to routines for evaluating chemical reactions, transport properties, * Output of runtime ``Diagnostics`` This section provides relevant notes on using these utilities across the Pele family of codes. There may be additional subtleties based on the code-specific implementation of these capabilities, in which case it will be necessary to refer to the documentation of the code of interest. - + Premixed Flame Initialization ============================= @@ -26,7 +26,7 @@ This code has two runtime parameters that may be set in the input file: :: pmf.do_cellAverage = 1 The first parameter specifies the path to the PMF data file. This file contains a two-line header followed by whitespace-delimited data columns in the order: position (cm), temperature (K), velocity (cm/s), density(g/cm3), species mole fractions. Sample files are provided in the relevant PeleLMeX (and PeleC) examples, and the procedure to generate these files with a provided script is described below. The second parameter specifies whether the PMF code does a finite volume-style integral over the queried cell (``pmf.do_cellAverage = 1``) or whether the code finds an interpolated value at the midpoint of the queried cell (``pmf.do_cellAverage = 0``) - + Generating a PMF file ~~~~~~~~~~~~~~~~~~~~~ @@ -48,7 +48,7 @@ An additional script is provided to allow plotting of the PMF solutions. This sc Turbulent Inflows ================= -Placeholder. PelePhysics supports the capability of the flow solvers to have spatially and temporally varying inflow conditions based on precomputed turbulence data. +Placeholder. PelePhysics supports the capability of the flow solvers to have spatially and temporally varying inflow conditions based on precomputed turbulence data. Generating a turbulence file ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -64,3 +64,50 @@ Diagnostics =========== Placeholder. Once the porting of diagnostics from PeleLMeX to PelePhysics/PeleC is complete, documentation can be added here. + +Filter +====== + +A utility for filtering data stored in AMReX data structures. When initializing the ``Filter`` class, the filter type +and filter width to grid ratio are specified. A variety of filter types are supported: + +* ``type = 0``: no filtering +* ``type = 1``: standard box filter +* ``type = 2``: standard Gaussian filter + +We have also implemented a set of filters defined in Sagaut & Grohens (1999) Int. J. Num. Meth. Fluids: + +* ``type = 3``: 3 point box filter approximation (Eq. 26) +* ``type = 4``: 5 point box filter approximation (Eq. 27) +* ``type = 5``: 3 point box filter optimized approximation (Table 1) +* ``type = 6``: 5 point box filter optimized approximation (Table 1) +* ``type = 7``: 3 point Gaussian filter approximation +* ``type = 8``: 5 point Gaussian filter approximation (Eq. 29) +* ``type = 9``: 3 point Gaussian filter optimized approximation (Table 1) +* ``type = 10``: 5 point Gaussian filter optimized approximation (Table 1) + +.. warning:: This utility is not aware of EB or domain boundaries. If the filter stencil extends across these boundaries, + the boundary cells are treated as if they are fluid cells. + It is up to the user to ensure an adequate number of ghost cells in the arrays are appropriately populated, + using the ``get_filter_ngrow()`` member function of the class to determine the required number of ghost cells. + + +Developing +~~~~~~~~~~ + +The weights for these filters are set in ``Filter.cpp``. To add a +filter type, one needs to add an enum to the ``filter_types`` and +define a corresponding ``set_NAME_weights`` function to be called at +initialization. + +The application of a filter can be done on a Fab or MultiFab. The loop nesting +ordering was chosen to be performant on existing HPC architectures and +discussed in PeleC milestone reports. An example call to the filtering operation is + +:: + + les_filter = Filter(les_filter_type, les_filter_fgr); + ... + les_filter.apply_filter(bxtmp, flux[i], filtered_flux[i], Density, NUM_STATE); + +The user must ensure that the correct number of grow cells is present in the Fab or MultiFab. diff --git a/Utility/Filter/Filter.H b/Utility/Filter/Filter.H new file mode 100644 index 000000000..9ec2c674e --- /dev/null +++ b/Utility/Filter/Filter.H @@ -0,0 +1,168 @@ +#ifndef FILTER_H +#define FILTER_H + +#include +#include +#include + +#ifdef AMREX_USE_OMP +#include +#endif + +// Filter types +enum filter_types { + no_filter = 0, + box, // 1 + gaussian, // 2 + box_3pt_approx, // 3 + box_5pt_approx, // 4 + box_3pt_optimized_approx, // 5 + box_5pt_optimized_approx, // 6 + gaussian_3pt_approx, // 7 + gaussian_5pt_approx, // 8 + gaussian_3pt_optimized_approx, // 9 + gaussian_5pt_optimized_approx, // 10 + num_filter_types +}; + +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +void +run_filter( + const int i, + const int j, + const int k, + const int nc, + const int ng, + const int nstart, + const amrex::Real* w, + amrex::Array4 const& q, + amrex::Array4 const& qh) +{ + for (int n = -ng; n <= ng; n++) { + for (int m = -ng; m <= ng; m++) { + for (int l = -ng; l <= ng; l++) { + qh(i, j, k, nc + nstart) += w[l + ng] * w[m + ng] * w[n + ng] * + q(i + l, j + m, k + n, nc + nstart); + } + } + } +} + +class Filter +{ + +public: + explicit Filter(const int type = box, const int fgr = 2) + : _type(type), _fgr(fgr) + { + + switch (_type) { + + case box: + set_box_weights(); + break; + + case gaussian: + set_gaussian_weights(); + break; + + case box_3pt_approx: + case gaussian_3pt_approx: // same as box_3pt_approx + set_box_3pt_approx_weights(); + break; + + case box_5pt_approx: + set_box_5pt_approx_weights(); + break; + + case box_3pt_optimized_approx: + set_box_3pt_optimized_approx_weights(); + break; + + case box_5pt_optimized_approx: + set_box_5pt_optimized_approx_weights(); + break; + + case gaussian_5pt_approx: + set_gaussian_5pt_approx_weights(); + break; + + case gaussian_3pt_optimized_approx: + set_gaussian_3pt_optimized_approx_weights(); + break; + + case gaussian_5pt_optimized_approx: + set_gaussian_5pt_optimized_approx_weights(); + break; + + case no_filter: + default: + _fgr = 1; + _ngrow = 0; + _nweights = 2 * _ngrow + 1; + _weights.resize(_nweights); + _weights[0] = 1.; + break; + + } // end switch + } + + // Default destructor + ~Filter() = default; + + int get_filter_ngrow() const { return _ngrow; } + + void apply_filter(const amrex::MultiFab& in, amrex::MultiFab& out); + + void apply_filter( + const amrex::MultiFab& in, + amrex::MultiFab& out, + const int nstart, + const int ncnt); + + void apply_filter( + const amrex::Box& cbox, const amrex::FArrayBox& in, amrex::FArrayBox& out); + + void apply_filter( + const amrex::Box& cbox, + const amrex::FArrayBox& in, + amrex::FArrayBox& out, + const int nstart, + const int ncnt); + + void apply_filter( + const amrex::Box& box, + const amrex::FArrayBox& in, + amrex::FArrayBox& out, + const int nstart, + const int ncnt, + const int ncomp); + +private: + int _type; + int _fgr; + int _ngrow; + int _nweights; + amrex::Vector _weights; + + void set_box_weights(); + + void set_gaussian_weights(); + + void set_box_3pt_approx_weights(); + + void set_box_5pt_approx_weights(); + + void set_box_3pt_optimized_approx_weights(); + + void set_box_5pt_optimized_approx_weights(); + + void set_gaussian_5pt_approx_weights(); + + void set_gaussian_3pt_optimized_approx_weights(); + + void set_gaussian_5pt_optimized_approx_weights(); +}; + +#endif /*_FILTER_H*/ diff --git a/Utility/Filter/Filter.cpp b/Utility/Filter/Filter.cpp new file mode 100644 index 000000000..a2dd940e9 --- /dev/null +++ b/Utility/Filter/Filter.cpp @@ -0,0 +1,491 @@ +#include "Filter.H" + +// Set the filter weights for the standard box filter +void +Filter::set_box_weights() +{ + + AMREX_ASSERT(_fgr % 2 == 0 || _fgr == 1); + _ngrow = _fgr / 2; + _nweights = 2 * _ngrow + 1; + _weights.resize(_nweights); + + // Set the weights + for (int i = 0; i < _nweights; i++) { + _weights[i] = 1.0 / _fgr; + } + + // Only half the cell is used at the ends + if (_fgr > 1) { + _weights[0] = 0.5 * _weights[0]; + _weights[_nweights - 1] = _weights[0]; + } +} + +// Set the filter weights for the standard gaussian filter +void +Filter::set_gaussian_weights() +{ + + AMREX_ASSERT(_fgr % 2 == 0); + _ngrow = _fgr / 2; + _nweights = 2 * _ngrow + 1; + _weights.resize(_nweights); + + // Set the weights + amrex::Real gamma = 6.0; + amrex::Real sigma = std::sqrt(1.0 / (2.0 * gamma)) * _fgr; + for (int i = 0; i < _nweights; i++) { + // equivalent to: std::sqrt(gamma / (PI * _fgr * _fgr)) * std::exp((-gamma + // * (i-_ngrow) * (i-_ngrow)) / (_fgr * _fgr)); + _weights[i] = + 1.0 / (std::sqrt(2.0 * amrex::Math::pi()) * sigma) * + std::exp((-(i - _ngrow) * (i - _ngrow)) / (2 * sigma * sigma)); + } + // normalize to ensure it all sums to one + amrex::Real sum = std::accumulate(_weights.begin(), _weights.end(), 0.0); + for (int i = 0; i < _nweights; i++) { + _weights[i] /= sum; + } +} + +// Set the filter weights for the 3pt polynomial truncation +// approximation of the box filter. See Eq. 26 in Sagaut & Grohens +// (1999) Int. J. Num. Meth. Fluids. +void +Filter::set_box_3pt_approx_weights() +{ + + _ngrow = 1; + _nweights = 2 * _ngrow + 1; + _weights.resize(_nweights); + + // Set the weights + _weights[0] = _fgr * _fgr / 24.0; + _weights[1] = (12.0 - _fgr * _fgr) / 12.0; + _weights[2] = _weights[0]; +} + +// Set the filter weights for the 5pt polynomial truncation +// approximation of the box filter. See Eq. 27 in Sagaut & Grohens +// (1999) Int. J. Num. Meth. Fluids (though there are typos). +void +Filter::set_box_5pt_approx_weights() +{ + + _ngrow = 2; + _nweights = 2 * _ngrow + 1; + _weights.resize(_nweights); + + // Set the weights + int _fgr2 = _fgr * _fgr; + int _fgr4 = _fgr2 * _fgr2; + _weights[0] = (3.0 * _fgr4 - 20.0 * _fgr2) / 5760.0; + _weights[1] = (80.0 * _fgr2 - 3.0 * _fgr4) / 1440.0; + _weights[2] = (3.0 * _fgr4 - 100.0 * _fgr2 + 960.0) / 960.0; + _weights[3] = _weights[1]; + _weights[4] = _weights[0]; +} + +// Set the filter weights for the 3pt optimized approximation of the +// box filter. See Table I in Sagaut & Grohens (1999) +// Int. J. Num. Meth. Fluids. +void +Filter::set_box_3pt_optimized_approx_weights() +{ + + _ngrow = 1; + _nweights = 2 * _ngrow + 1; + _weights.resize(_nweights); + + amrex::Real ratio; + switch (_fgr) { + + case 1: + ratio = 0.079; + break; + + case 2: + ratio = 0.274; + break; + + case 3: + ratio = 1.377; + break; + + case 4: + ratio = -2.375; + break; + + case 5: + ratio = -1.000; + break; + + case 6: + ratio = -0.779; + break; + + case 7: + ratio = -0.680; + break; + + case 8: + ratio = -0.627; + break; + + case 9: + ratio = -0.596; + break; + + case 10: + ratio = -0.575; + break; + + default: // default to standard box filter + set_box_weights(); + return; + break; + + } // end switch + + // Set the weights + _weights[0] = ratio / (1 + 2.0 * ratio); + _weights[1] = 1.0 - 2.0 * _weights[0]; + _weights[2] = _weights[0]; +} + +// Set the filter weights for the 5pt optimized approximation of the +// box filter. See Table I in Sagaut & Grohens (1999) +// Int. J. Num. Meth. Fluids. +void +Filter::set_box_5pt_optimized_approx_weights() +{ + + _ngrow = 2; + _nweights = 2 * _ngrow + 1; + _weights.resize(_nweights); + + amrex::Real ratio1; + amrex::Real ratio2; + switch (_fgr) { + + case 1: + ratio1 = 0.0886; + ratio2 = -0.0169; + break; + + case 2: + ratio1 = 0.3178; + ratio2 = -0.0130; + break; + + case 3: + ratio1 = 1.0237; + ratio2 = 0.0368; + break; + + case 4: + ratio1 = 2.4414; + ratio2 = 0.5559; + break; + + case 5: + ratio1 = 0.2949; + ratio2 = 0.7096; + break; + + case 6: + ratio1 = -0.5276; + ratio2 = 0.4437; + break; + + case 7: + ratio1 = -0.6708; + ratio2 = 0.3302; + break; + + case 8: + ratio1 = -0.7003; + ratio2 = 0.2767; + break; + + case 9: + ratio1 = -0.7077; + ratio2 = 0.2532; + break; + + case 10: + ratio1 = -0.6996; + ratio2 = 0.2222; + break; + + default: // default to standard box filter + set_box_weights(); + return; + break; + + } // end switch + + // Set the weights + _weights[0] = ratio2 / (1 + 2.0 * ratio1 + 2.0 * ratio2); + _weights[1] = ratio1 / ratio2 * _weights[0]; + _weights[2] = 1.0 - 2.0 * _weights[0] - 2.0 * _weights[1]; + _weights[3] = _weights[1]; + _weights[4] = _weights[0]; +} + +// Set the filter weights for the 5pt polynomial truncation +// approximation of the gaussian filter. See Eq. 29 in Sagaut & +// Grohens (1999) Int. J. Num. Meth. Fluids. +void +Filter::set_gaussian_5pt_approx_weights() +{ + + _ngrow = 2; + _nweights = 2 * _ngrow + 1; + _weights.resize(_nweights); + + // Set the weights + int _fgr2 = _fgr * _fgr; + int _fgr4 = _fgr2 * _fgr2; + _weights[0] = (_fgr4 - 4.0 * _fgr2) / 1152.0; + _weights[1] = (16.0 * _fgr2 - _fgr4) / 288.0; + _weights[2] = (_fgr4 - 20.0 * _fgr2 + 192.0) / 192.0; + _weights[3] = _weights[1]; + _weights[4] = _weights[0]; +} + +// Set the filter weights for the 3pt optimized approximation of the +// gaussian filter. See Table I in Sagaut & Grohens (1999) +// Int. J. Num. Meth. Fluids. +void +Filter::set_gaussian_3pt_optimized_approx_weights() +{ + + _ngrow = 1; + _nweights = 2 * _ngrow + 1; + _weights.resize(_nweights); + + amrex::Real ratio; + switch (_fgr) { + + case 1: + ratio = 0.0763; + break; + + case 2: + ratio = 0.2527; + break; + + case 3: + ratio = 1.1160; + break; + + case 4: + ratio = -3.144; + break; + + case 5: + ratio = -1.102; + break; + + case 6: + ratio = -0.809; + break; + + case 7: + ratio = -0.696; + break; + + case 8: + ratio = -0.638; + break; + + case 9: + ratio = -0.604; + break; + + case 10: + ratio = -0.581; + break; + + default: // default to the 3pt gaussian filter + set_box_3pt_approx_weights(); + return; + break; + + } // end switch + + // Set the weights + _weights[0] = ratio / (1 + 2.0 * ratio); + _weights[1] = 1.0 - 2.0 * _weights[0]; + _weights[2] = _weights[0]; +} + +// Set the filter weights for the 5pt optimized approximation of the +// gaussian filter. See Table I in Sagaut & Grohens (1999) +// Int. J. Num. Meth. Fluids. +void +Filter::set_gaussian_5pt_optimized_approx_weights() +{ + + _ngrow = 2; + _nweights = 2 * _ngrow + 1; + _weights.resize(_nweights); + + amrex::Real ratio1; + amrex::Real ratio2; + switch (_fgr) { + + case 1: + ratio1 = 0.0871; + ratio2 = -0.0175; + break; + + case 2: + ratio1 = 0.2596; + ratio2 = -0.0021; + break; + + case 3: + ratio1 = 0.4740; + ratio2 = 0.0785; + break; + + case 4: + ratio1 = 0.1036; + ratio2 = 0.2611; + break; + + case 5: + ratio1 = -0.4252; + ratio2 = 0.3007; + break; + + case 6: + ratio1 = -0.6134; + ratio2 = 0.2696; + break; + + case 7: + ratio1 = -0.6679; + ratio2 = 0.2419; + break; + + case 8: + ratio1 = -0.6836; + ratio2 = 0.2231; + break; + + case 9: + ratio1 = -0.6873; + ratio2 = 0.2103; + break; + + case 10: + ratio1 = -0.6870; + ratio2 = 0.2014; + break; + + default: // default to 5pt gaussian filter + set_gaussian_5pt_approx_weights(); + return; + break; + + } // end switch + + // Set the weights + _weights[0] = ratio2 / (1 + 2.0 * ratio1 + 2.0 * ratio2); + _weights[1] = ratio1 / ratio2 * _weights[0]; + _weights[2] = 1.0 - 2.0 * _weights[0] - 2.0 * _weights[1]; + _weights[3] = _weights[1]; + _weights[4] = _weights[0]; +} + +// Run the filtering operation on a MultiFab +void +Filter::apply_filter(const amrex::MultiFab& in, amrex::MultiFab& out) +{ + apply_filter(in, out, 0, out.nComp()); +} + +void +Filter::apply_filter( + const amrex::MultiFab& in, + amrex::MultiFab& out, + const int nstart, + const int ncnt) +{ + + // Ensure enough grow cells + AMREX_ASSERT(in.nGrow() >= out.nGrow() + _ngrow); + + const auto& ins = in.const_arrays(); + const auto& outs = out.arrays(); + out.setVal(0, nstart, ncnt); + + amrex::Gpu::DeviceVector weights(_weights.size()); + amrex::Real* w = weights.data(); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, _weights.begin(), _weights.end(), + weights.begin()); + const int captured_ngrow = _ngrow; + + const amrex::IntVect ngs(out.nGrow()); + amrex::ParallelFor( + in, ngs, ncnt - nstart, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k, int nc) noexcept { + run_filter(i, j, k, nc, captured_ngrow, nstart, w, ins[nbx], outs[nbx]); + }); + amrex::Gpu::synchronize(); +} + +// Run the filtering operation on a FAB +void +Filter::apply_filter( + const amrex::Box& cbox, const amrex::FArrayBox& in, amrex::FArrayBox& out) +{ + apply_filter(cbox, in, out, 0, out.nComp()); +} + +void +Filter::apply_filter( + const amrex::Box& cbox, + const amrex::FArrayBox& in, + amrex::FArrayBox& out, + const int nstart, + const int ncnt) +{ + BL_PROFILE("Filter::apply_filter()"); + AMREX_ASSERT(in.nComp() == out.nComp()); + + const int ncomp = in.nComp(); + apply_filter(cbox, in, out, nstart, ncnt, ncomp); +} + +void +Filter::apply_filter( + const amrex::Box& box, + const amrex::FArrayBox& in, + amrex::FArrayBox& out, + const int nstart, + const int ncnt, + const int /*ncomp*/) +{ + const auto q = in.const_array(); + out.setVal(0.0, box, nstart, ncnt); + auto qh = out.array(); + + amrex::Gpu::DeviceVector weights(_weights.size()); + amrex::Real* w = weights.data(); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, _weights.begin(), _weights.end(), + weights.begin()); + + const int captured_ngrow = _ngrow; + amrex::ParallelFor( + box, ncnt - nstart, + [=] AMREX_GPU_DEVICE(int i, int j, int k, int nc) noexcept { + run_filter(i, j, k, nc, captured_ngrow, nstart, w, q, qh); + }); +} diff --git a/Utility/Filter/Make.package b/Utility/Filter/Make.package new file mode 100644 index 000000000..147f76e85 --- /dev/null +++ b/Utility/Filter/Make.package @@ -0,0 +1,5 @@ +CEXE_sources += Filter.cpp +CEXE_headers += Filter.H + +VPATH_LOCATIONS += $(PELE_PHYSICS_HOME)/Utility/Filter +INCLUDE_LOCATIONS += $(PELE_PHYSICS_HOME)/Utility/Filter