User specified tagging criteria (#1673)

Author: marchdf

amr-wind/utilities/tagging/CMakeLists.txt

@@ -13,4 +13,5 @@ target_sources(${amr_wind_lib_name} PRIVATE
   # Geometry based refinement options
   BoxRefiner.cpp
   CylinderRefiner.cpp
+  UDFRefiner.cpp
   )

amr-wind/utilities/tagging/UDFRefiner.H

@@ -0,0 +1,38 @@
+#ifndef UDFREFINER_H
+#define UDFREFINER_H
+
+#include "amr-wind/CFDSim.H"
+#include "amr-wind/utilities/tagging/GeometryRefinement.H"
+#include "AMReX_Parser.H"
+
+namespace amr_wind::tagging {
+
+/** AMR refinement using a user defined expression
+ *  \ingroup amr_utils
+ *
+ */
+class UDFRefiner : public GeometryType::Register<UDFRefiner>
+{
+public:
+    static std::string identifier() { return "udf"; }
+
+    explicit UDFRefiner(const CFDSim& sim, const std::string& key);
+
+    ~UDFRefiner() override = default;
+
+    void operator()(
+        const amrex::Box& /*bx*/,
+        const amrex::Geometry& geom,
+        const amrex::Array4<amrex::TagBox::TagType>& tags) const override;
+
+    const amrex::RealBox& bound_box() const override { return m_bound_box; };
+
+private:
+    const CFDSim& m_sim;
+    amrex::RealBox m_bound_box;
+    amrex::Parser m_parser;
+};
+
+} // namespace amr_wind::tagging
+
+#endif /* UDFREFINER_H */

amr-wind/utilities/tagging/UDFRefiner.cpp

@@ -0,0 +1,48 @@
+#include "amr-wind/utilities/tagging/UDFRefiner.H"
+#include "AMReX_ParmParse.H"
+
+namespace amr_wind::tagging {
+
+UDFRefiner::UDFRefiner(const CFDSim& sim, const std::string& key)
+    : m_sim(sim), m_bound_box(m_sim.repo().mesh().Geom(0).ProbDomain())
+{
+    amrex::ParmParse pp(key);
+    std::string udf_str;
+    pp.get("udf", udf_str);
+
+    amrex::Vector<amrex::Real> box_lo(AMREX_SPACEDIM, 0);
+    amrex::Vector<amrex::Real> box_hi(AMREX_SPACEDIM, 0);
+    if (pp.queryarr("box_lo", box_lo, 0, static_cast<int>(box_lo.size())) ==
+        1) {
+        m_bound_box.setLo(box_lo);
+    }
+    if (pp.queryarr("box_hi", box_hi, 0, static_cast<int>(box_hi.size())) ==
+        1) {
+        m_bound_box.setHi(box_hi);
+    }
+    m_parser.define(udf_str);
+    m_parser.registerVariables({"t", "x", "y", "z"});
+}
+
+void UDFRefiner::operator()(
+    const amrex::Box& bx,
+    const amrex::Geometry& geom,
+    const amrex::Array4<amrex::TagBox::TagType>& tag) const
+{
+    const auto& problo = geom.ProbLoArray();
+    const auto& dx = geom.CellSizeArray();
+    const auto time = m_sim.time().current_time();
+    auto udf_func = m_parser.compile<AMREX_SPACEDIM + 1>();
+
+    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+        const amrex::Real x = problo[0] + (i + 0.5) * dx[0];
+        const amrex::Real y = problo[1] + (j + 0.5) * dx[1];
+        const amrex::Real z = problo[2] + (k + 0.5) * dx[2];
+        const auto val = udf_func(time, x, y, z);
+        if (static_cast<bool>(val)) {
+            tag(i, j, k) = amrex::TagBox::SET;
+        }
+    });
+}
+
+} // namespace amr_wind::tagging

amr-wind/wind_energy/MOData.cpp

@@ -88,11 +88,6 @@ void MOData::update_fluxes(int max_iters)
         ++iter;
     } while ((std::abs(utau_iter - utau) > 1e-5) && iter <= max_iters);
 
-    amrex::Print() << "MOData::update_fluxes: output\n"
-                   << "  utau = " << utau << "\n"
-                   << "  q = " << surf_temp_flux << "\n"
-                   << "  L = " << obukhov_len << std::endl;
-
     if (iter >= max_iters) {
         amrex::Print()
             << "MOData::update_fluxes: Convergence criteria not met after "

docs/sphinx/user/inputs_tagging.rst

@@ -122,7 +122,8 @@ Refinement using geometry
 
 This section controls refinement using pre-defined geometric shapes. Currently,
 two options are supported: 1. ``box`` -- refines the region inside a hexahedral
-block, and 2. ``cylinder`` -- refines the region inside a cylindrical block.
+block, 2. ``cylinder`` -- refines the region inside a cylindrical block, and
+3. ``udf`` -- refines along an analytical function depending on time and spatial coordinate.
 
 .. input_param:: tagging.GeometryRefinement.shapes
 
@@ -181,11 +182,21 @@ Example::
   tagging.g2.c1.outer_radius = 300.0
   tagging.g2.c1.inner_radius = 275.0
 
+  tagging.g3.type = GeometryRefinement
+  tagging.g3.shapes = udf0
+  tagging.g3.level = 0
+  tagging.g3.udf0.type = udf
+  tagging.g3.udf0.udf = "if(x > 500, 1.0, 0.0)"
+  tagging.g3.udf0.box_lo = 0.0 0.0 0.0
+  tagging.g3.udf0.box_hi = 1000.0 1000.0 500.0
 
-This example defines two different refinement definitions acting on level 0 and
-1 respectively. The refinement at level 0 (``g1``) contains two box regions,
-whereas the refinement at level 1 (``g2``) only contains one cylinder
-definition.
+
+This example defines three different refinement definitions acting on
+levels 0, 1 and 0 respectively. The first refinement at level 0
+(``g1``) contains two box regions, whereas the refinement at level 1
+(``g2``) only contains one cylinder definition. The second refinement
+at level 0 (``g3``) uses a user defined function (UDF) to specify a
+geometrical refinement region.
 
 **Refinement using hexahedral block definitions**
 
@@ -212,6 +223,41 @@ The axis and the extents along the axis are defined by two position vectors
 optional ``inner_radius`` can be specified to restrict tagging to an annulus
 between the inner and outer radii.
 
+**Refinement using a user specified function**
+
+The UDF uses AMReX's Parser class to parse a string in the input file,
+``udf`` key in the input file. For more information on the AMReX
+Parser and supported functions, consult the `AMReX Parser
+documentation
+<https://amrex-codes.github.io/amrex/docs_html/Basics.html#parser>`_.
+The UDF can be (almost) arbitrarily complex, leveraging standard
+functions, boolean operators, local variables, etc. The following
+assumptions are imposed on the UDF: 1. it must evaluate to 0 (tagging
+off) or 1 (tagging on) (the return value of the UDF is cast to a
+boolean and that is used to set the tags), 2. the allowed variables
+are ``t``, ``x``, ``y``, and ``z``.
+
+.. input_param:: tagging.GeometryRefinement.UDFRefiner.udf
+
+   **type:** String, required
+
+   String specifying the user defined function.
+
+.. input_param:: tagging.GeometryRefinement.UDFRefiner.box_lo
+
+   **type:** Vector<Real>, optional
+
+   List of the low corner values for a bounding box where the tagging
+   will be active. By default the bounding box will span the entire domain.
+
+.. input_param:: tagging.GeometryRefinement.UDFRefiner.box_hi
+
+   **type:** Vector<Real>, optional
+
+   List of the high corner values for a bounding box where the tagging
+   will be active. By default the bounding box will span the entire domain.
+
+
 Refinement using Q-Criterion
 `````````````````````````````````````
 

test/CMakeLists.txt

@@ -287,6 +287,7 @@ add_test_re(ctv_mol_mesh_map_explicit)
 add_test_re(abl_specifiedmol_neutral)
 add_test_re(abl_specifiedmol_stable)
 add_test_re(abl_specifiedmol_unstable)
+add_test_re(udf_refinement)
 
 if(AMR_WIND_ENABLE_NETCDF)
   add_test_re(abl_bndry_output)

test/test_files/udf_refinement/udf_refinement.inp

@@ -0,0 +1,84 @@
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            SIMULATION STOP            #
+#.......................................#
+time.stop_time               =   22000.0     # Max (simulated) time to evolve
+time.max_step                =   10          # Max number of time steps
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#         TIME STEP COMPUTATION         #
+#.......................................#
+time.fixed_dt         =   0.5        # Use this constant dt if > 0
+time.cfl              =   0.95         # CFL factor
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            INPUT AND OUTPUT           #
+#.......................................#
+time.plot_interval            =  10       # Steps between plot files
+time.checkpoint_interval      =  5       # Steps between checkpoint files
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#               PHYSICS                 #
+#.......................................#
+incflo.gravity          =   0.  0. -9.81  # Gravitational force (3D)
+incflo.density             = 1.0          # Reference density
+
+incflo.use_godunov = 1
+# default godunov_type is weno_z
+incflo.diffusion_type = 1
+transport.viscosity = 1.0e-5
+transport.laminar_prandtl = 0.7
+transport.turbulent_prandtl = 0.3333
+transport.reference_temperature = 300.0
+turbulence.model = Smagorinsky
+Smagorinsky_coeffs.Cs = 0.135
+
+
+incflo.physics = ABL
+ICNS.source_terms = BoussinesqBuoyancy CoriolisForcing ABLForcing
+CoriolisForcing.latitude = 41.3
+ABLForcing.abl_forcing_height = 90
+
+incflo.velocity = 6.128355544951824  5.142300877492314 0.0
+
+ABL.temperature_heights = 650.0 750.0 1000.0
+ABL.temperature_values = 300.0 308.0 308.75
+
+ABL.kappa = .41
+ABL.surface_roughness_z0 = 0.15
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#        ADAPTIVE MESH REFINEMENT       #
+#.......................................#
+# amr.n_cell              = 48 48 48    # Grid cells at coarsest AMRlevel
+amr.n_cell              = 128 128 128    # Grid cells at coarsest AMRlevel
+amr.max_level           = 1           # Max AMR level in hierarchy
+amr.max_grid_size       = 8
+amr.n_error_buf         = 0
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#              GEOMETRY                 #
+#.......................................#
+geometry.prob_lo        =   0.       0.     0.  # Lo corner coordinates
+geometry.prob_hi        =   1000.  1000.  1000.  # Hi corner coordinates
+geometry.is_periodic    =   1   1   0   # Periodicity x y z (0/1)
+
+# Boundary conditions
+zlo.type =   "wall_model"
+
+zhi.type =   "slip_wall"
+zhi.temperature_type = "fixed_gradient"
+zhi.temperature = 0.003 # tracer is used to specify potential temperature gradient
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#              VERBOSITY                #
+#.......................................#
+incflo.verbose          =   0          # incflo_level
+
+tagging.labels = udf
+tagging.udf.type = GeometryRefinement
+tagging.udf.shapes = udf0
+tagging.udf.level  = 0
+tagging.udf.udf0.type = udf
+tagging.udf.udf0.udf = "x0 = 500.0; y0 = 500.0; z0 = 500.0; r1 = 200.0; r2 = 250.0; xp = x - x0; yp = y - y0; zp = z - z0; d2 = (xp^2 + yp^2 + zp^2); if((r1**2 <= d2) and (d2 <= r2^2), 1.0, 0.0)"
+tagging.udf.udf0.box_lo = 0.0 0.0 0.0
+tagging.udf.udf0.box_hi = 1000.0 1000.0 500.0