Filter Utility

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.

Utility/Filter/Filter.H

@@ -0,0 +1,168 @@
+#ifndef FILTER_H
+#define FILTER_H
+
+#include <AMReX_REAL.H>
+#include <AMReX_Array.H>
+#include <AMReX_MultiFab.H>
+
+#ifdef AMREX_USE_OMP
+#include <omp.h>
+#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 amrex::Real> const& q,
+  amrex::Array4<amrex::Real> 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<amrex::Real> _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*/