AMRErrorTag: Add PARSER (#4534)

Author: WeiqunZhang

Src/AmrCore/AMReX_ErrorList.H

@@ -3,12 +3,13 @@
 #include <AMReX_Config.H>
 
 #include <AMReX_Array.H>
-#include <AMReX_Vector.H>
-#include <AMReX_REAL.H>
 #include <AMReX_ArrayLim.H>
+#include <AMReX_Geometry.H>
 #include <AMReX_MultiFab.H>
+#include <AMReX_Parser.H>
+#include <AMReX_REAL.H>
 #include <AMReX_TagBox.H>
-#include <AMReX_Geometry.H>
+#include <AMReX_Vector.H>
 
 #include <string>
 #include <memory>
@@ -426,7 +427,7 @@ std::ostream& operator << (std::ostream& os, const ErrorList& elst);
   {
   public:
 
-    enum TEST {GRAD=0, RELGRAD, LESS, GREATER, VORT, BOX, USER};
+    enum TEST {GRAD=0, RELGRAD, LESS, GREATER, VORT, BOX, USER, PARSER};
 
     struct UserFunc
     {
@@ -482,6 +483,12 @@ std::ostream& operator << (std::ostream& os, const ErrorList& elst);
                  const AMRErrorTagInfo& info = AMRErrorTagInfo()) noexcept
       : m_userfunc(userfunc), m_field(std::move(field)), m_info(info), m_ngrow(ngrow) {}
 
+    //! Construct AMRErrorTag from a Parser that returns Real for given
+    //! (x,y,z,t). This class's operator() can then be used to tag regions
+    //! where the Parser's return value is greater than 0.
+    explicit AMRErrorTag (Parser parser,
+                          const AMRErrorTagInfo& info = AMRErrorTagInfo());
+
     void operator() (amrex::TagBoxArray&    tb,
                      const amrex::MultiFab* mf,
                      char                   clearval,
@@ -503,9 +510,11 @@ std::ostream& operator << (std::ostream& os, const ErrorList& elst);
     Vector<Real> m_value;
     TEST m_test{BOX};
     UserFunc* m_userfunc = nullptr;
+    std::unique_ptr<Parser> m_parser;
+    ParserExecutor<4> m_parser_exe; // 4: (x,y,z,t)
     std::string m_field;
     AMRErrorTagInfo m_info;
-    int m_ngrow;
+    int m_ngrow = 0;
   };
 }
 

Src/AmrCore/AMReX_ErrorList.cpp

@@ -214,12 +214,21 @@ operator << (std::ostream&    os,
     return os;
 }
 
+AMRErrorTag::AMRErrorTag (Parser parser, const AMRErrorTagInfo& info)
+    : m_test(PARSER),
+      m_parser(std::make_unique<Parser>(std::move(parser))),
+      m_parser_exe(m_parser->compile<4>()),
+      m_info(info)
+{}
+
 int
 AMRErrorTag::SetNGrow () const noexcept
 {
-    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_test != USER, "Do not call SetNGrow with USER test");
-    static std::map<TEST,int> ng = { {GRAD,1}, {RELGRAD,1}, {LESS,0}, {GREATER,0}, {VORT,0}, {BOX,0} };
-    return ng[m_test];
+    if (m_test == GRAD || m_test == RELGRAD) {
+        return 1;
+    } else {
+        return 0;
+    }
 }
 
 void
@@ -272,6 +281,30 @@ AMRErrorTag::operator() (TagBoxArray&    tba,
                     }
                 });
             }
+            else if (m_test == PARSER)
+            {
+                const auto& plo = geom.ProbLoArray();
+                const auto& dx  = geom.CellSizeArray();
+                auto tag_update = m_info.m_derefine ? clearval : tagval;
+                auto fn = m_parser_exe;
+                ParallelFor(tba, [=] AMREX_GPU_DEVICE (int bi, int i, int j, int k) noexcept
+                {
+                    auto x = plo[0]+(Real(i)+Real(0.5))*dx[0];
+#if (AMREX_SPACEDIM > 1)
+                    auto y = plo[1]+(Real(j)+Real(0.5))*dx[1];
+#else
+                    auto y = Real(0);
+#endif
+#if (AMREX_SPACEDIM == 3)
+                    auto z = plo[2]+(Real(k)+Real(0.5))*dx[2];
+#else
+                    auto z = Real(0);
+#endif
+                    if (fn(x,y,z,time) > 0) {
+                        tagma[bi](i,j,k) = tag_update;
+                    }
+                });
+            }
             else
             {
                 auto const& datma   = mf->const_arrays();

Tests/Amr/Advection_AmrCore/Exec/inputs-errorfn

@@ -0,0 +1,85 @@
+# *****************************************************************
+# Run until nsteps == max_step or time == stop_time,
+#     whichever comes first
+# *****************************************************************
+max_step  = 1000
+stop_time = 2.0
+
+# *****************************************************************
+# Are we restarting from an existing checkpoint file?
+# *****************************************************************
+#amr.restart  = chk00060 # restart from this checkpoint file
+
+# *****************************************************************
+# Problem size and geometry
+# *****************************************************************
+geometry.prob_lo     =  0.0  0.0  0.0
+geometry.prob_hi     =  1.0  1.0  0.125
+geometry.is_periodic =  1    1    1
+
+# *****************************************************************
+# VERBOSITY
+# *****************************************************************
+amr.v              = 1       # verbosity in Amr
+
+# *****************************************************************
+# Resolution and refinement
+# *****************************************************************
+amr.n_cell          = 64 64 8
+amr.max_level       = 2       # maximum level number allowed --
+                              # number of levels = max_level + 1
+
+amr.ref_ratio       = 2 2 2 2 # refinement ratio between levels
+
+# *****************************************************************
+# Control of grid creation
+# *****************************************************************
+# Blocking factor for grid creation in each dimension --
+#   this ensures that every grid is coarsenable by a factor of 8 --
+#   this is mostly relevant for multigrid performance
+amr.blocking_factor_x = 8
+amr.blocking_factor_y = 8
+amr.blocking_factor_z = 8
+
+amr.max_grid_size   = 16
+
+amr.regrid_int      = 2       # how often to regrid
+
+# *****************************************************************
+# Time step control
+# *****************************************************************
+adv.cfl            = 0.7     # CFL constraint for explicit advection
+
+adv.do_subcycle    = 1       # Do we subcycle in time?
+
+# *****************************************************************
+# Should we reflux at coarse-fine boundaries?
+# *****************************************************************
+adv.do_reflux = 1
+
+# *****************************************************************
+# Tagging -  if phi > 1.01 at level 0, then refine
+#            if phi > 1.1  at level 1, then refine
+#            if phi > 1.5  at level 2, then refine
+# *****************************************************************
+adv.phierr = 1.01  1.1  1.5
+adv.errfn = "1.e-3-((x-0.2*t)**2+(y-0.3*t)**2+z**2)" # tag moving sphere
+
+# *****************************************************************
+# Plotfile name and frequency
+# *****************************************************************
+amr.plot_file  = plt    # root name of plot file
+amr.plot_int   =  10    # number of timesteps between plot files
+                        # if negative then no plot files will be written
+
+# *****************************************************************
+# Checkpoint name and frequency
+# *****************************************************************
+amr.chk_file = chk      # root name of checkpoint file
+amr.chk_int  = -1       # number of timesteps between checkpoint files
+                        # if negative then no checkpoint files will be written
+
+# *****************************************************************
+# Particles
+# *****************************************************************
+amr.do_tracers = 1      # Turn tracer particles on or off

Tests/Amr/Advection_AmrCore/Source/AmrCoreAdv.H

@@ -9,6 +9,7 @@
 
 #include <AMReX_AmrCore.H>
 #include <AMReX_BCRec.H>
+#include <AMReX_ErrorList.H>
 #include <AMReX_FluxRegister.H>
 #include <AMReX_FillPatcher.H>
 
@@ -215,6 +216,9 @@ private:
     //number of ghost cells on facevel
     int nghost = 2;
 
+    std::string error_fn;
+    amrex::Vector<amrex::AMRErrorTag> error_tag;
+
 #ifdef AMREX_PARTICLES
     void init_particles ();
     int       do_tracers = 0;

Tests/Amr/Advection_AmrCore/Source/AmrCoreAdv.cpp

@@ -90,6 +90,14 @@ AmrCoreAdv::AmrCoreAdv ()
     // fillpatcher[lev] is for filling data on level lev using the data on
     // lev-1 and lev.
     fillpatcher.resize(nlevs_max+1);
+
+    // amrex::AMRErrorTag supports a number of common tagging
+    // approaches. Here we set it up to use amrex::Parser.
+    if (! error_fn.empty()) {
+        amrex::Parser parser(error_fn);
+        parser.registerVariables({"x","y","z","t"});
+        error_tag.emplace_back(std::move(parser));
+    }
 }
 
 AmrCoreAdv::~AmrCoreAdv () = default;
@@ -349,7 +357,7 @@ void AmrCoreAdv::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba
 // tag all cells for refinement
 // overrides the pure virtual function in AmrCore
 void
-AmrCoreAdv::ErrorEst (int lev, TagBoxArray& tags, Real /*time*/, int /*ngrow*/)
+AmrCoreAdv::ErrorEst (int lev, TagBoxArray& tags, Real time, int /*ngrow*/)
 {
     static bool first = true;
     static Vector<Real> phierr;
@@ -370,18 +378,16 @@ AmrCoreAdv::ErrorEst (int lev, TagBoxArray& tags, Real /*time*/, int /*ngrow*/)
         }
     }
 
-    if (lev >= phierr.size()) { return; }
-
-//    const int clearval = TagBox::CLEAR;
+    const int clearval = TagBox::CLEAR;
     const int   tagval = TagBox::SET;
 
-    const MultiFab& state = phi_new[lev];
+    if (lev < phierr.size())
+    {
+        const MultiFab& state = phi_new[lev];
 
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if(Gpu::notInLaunchRegion())
 #endif
-    {
-
         for (MFIter mfi(state,TilingIfNotGPU()); mfi.isValid(); ++mfi)
         {
             const Box& bx  = mfi.tilebox();
@@ -396,6 +402,10 @@ AmrCoreAdv::ErrorEst (int lev, TagBoxArray& tags, Real /*time*/, int /*ngrow*/)
             });
         }
     }
+
+    for (auto const& etag : error_tag) {
+        etag(tags, nullptr, clearval, tagval, time, lev, geom[lev]);
+    }
 }
 
 // read in some parameters from inputs file
@@ -425,6 +435,7 @@ AmrCoreAdv::ReadParameters ()
         pp.query("cfl", cfl);
         pp.query("do_reflux", do_reflux);
         pp.query("do_subcycle", do_subcycle);
+        pp.query("errfn", error_fn);
     }
 
 #ifdef AMREX_PARTICLES