remove NMOIST from usage (#1779)

* remove NMOIST from usage

* remove unused

* fix for incompressible

Source/BoundaryConditions/ABLMost.H

@@ -47,14 +47,16 @@ public:
                       amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Hwave,
                       amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Lwave,
                       amrex::Vector<std::unique_ptr<amrex::MultiFab>>& eddyDiffs,
+                      int n_qstate,
                       amrex::Real start_bdy_time = 0.0,
                       amrex::Real bdy_time_interval = 0.0)
     : m_exp_most(use_exp_most),
       m_rotate(use_rot_most),
       m_start_bdy_time(start_bdy_time),
       m_bdy_time_interval(bdy_time_interval),
       m_geom(geom),
-      m_ma(geom,vars_old,Theta_prim,Qv_prim,Qr_prim,z_phys_nd)
+      m_ma(geom,vars_old,Theta_prim,Qv_prim,Qr_prim,z_phys_nd),
+      m_n_qstate(n_qstate)
     {
         // We have a moisture model if Qv_prim is a valid pointer
         use_moisture = (Qv_prim[0].get());
@@ -527,6 +529,8 @@ private:
     amrex::Vector<amrex::MultiFab*>  m_Hwave_lev;
     amrex::Vector<amrex::MultiFab*>  m_Lwave_lev;
     amrex::Vector<amrex::MultiFab*>  m_eddyDiffs_lev;
+
+    int m_n_qstate;
 };
 
 #endif /* ABLMOST_H */

Source/BoundaryConditions/ABLMost.cpp

@@ -577,12 +577,9 @@ ABLMost::compute_pblh (const int& lev,
                        MultiFab* z_phys_cc,
                        const PBLHeightEstimator& est)
 {
-    int moist_flag = 0;
-    int n_qstate = vars[lev][Vars::cons].nComp() - (NVAR_max - NMOIST_max);
-    if (n_qstate > 3) {
-        moist_flag = (n_qstate > 3) ? RhoQ4_comp : RhoQ3_comp;
-    } else if (n_qstate > 0) {
-        moist_flag = 1;
+    int moist_flag = (m_n_qstate > 0) ? 1 : 0;
+    if (m_n_qstate > 3) {
+        moist_flag = RhoQ4_comp;
     }
 
     est.compute_pblh(m_geom[lev],z_phys_cc, pblh[lev].get(),

Source/Diffusion/Diffusion.H

@@ -66,7 +66,8 @@ void DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain,
                              const amrex::Array4<const amrex::Real>& tm_arr,
                              const amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> grav_gpu,
                              const amrex::BCRec* bc_ptr,
-                             const bool use_most);
+                             const bool use_most,
+                             const int n_qstate);
 
 void DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain,
                              int start_comp, int num_comp,
@@ -104,7 +105,8 @@ void DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain,
                              const amrex::Array4<const amrex::Real>& tm_arr,
                              const amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> grav_gpu,
                              const amrex::BCRec* bc_ptr,
-                             const bool use_most);
+                             const bool use_most,
+                             const int n_qstate);
 
 
 

Source/Diffusion/DiffusionSrcForState_N.cpp

@@ -63,7 +63,8 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
                         const Array4<const Real>& tm_arr,
                         const GpuArray<Real,AMREX_SPACEDIM> grav_gpu,
                         const BCRec* bc_ptr,
-                        const bool use_most)
+                        const bool use_most,
+                        const int n_qstate)
 {
     BL_PROFILE_VAR("DiffusionSrcForState_N()",DiffusionSrcForState_N);
 
@@ -666,12 +667,9 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
         int qty_index = RhoQKE_comp;
         auto pbl_mynn_B1_l = turbChoice.pbl_mynn.B1;
 
-        int moist_flag = 0;
-        int n_qstate = cell_data.nComp() - (NVAR_max - NMOIST_max);
+        int moist_comp = (n_qstate > 0) ? 1 : 0;
         if (n_qstate > 3) {
-            moist_flag = (n_qstate > 3) ? RhoQ4_comp : RhoQ3_comp;
-        } else if (n_qstate > 0) {
-            moist_flag = 1;
+            moist_comp = RhoQ4_comp;
         }
 
         ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -687,7 +685,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
             cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim,
                                                                   mu_turb,cellSizeInv,domain,
                                                                   pbl_mynn_B1_l,tm_arr(i,j,0),
-                                                                  moist_flag,
+                                                                  moist_comp,
                                                                   c_ext_dir_on_zlo, c_ext_dir_on_zhi,
                                                                   u_ext_dir_on_zlo, u_ext_dir_on_zhi,
                                                                   v_ext_dir_on_zlo, v_ext_dir_on_zhi);

Source/Diffusion/DiffusionSrcForState_T.cpp

@@ -76,7 +76,8 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
                         const Array4<const Real>& tm_arr,
                         const GpuArray<Real,AMREX_SPACEDIM> grav_gpu,
                         const BCRec* bc_ptr,
-                        const bool use_most)
+                        const bool use_most,
+                        const int n_qstate)
 {
     BL_PROFILE_VAR("DiffusionSrcForState_T()",DiffusionSrcForState_T);
 
@@ -732,12 +733,9 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
         int qty_index = RhoQKE_comp;
         auto pbl_mynn_B1_l = turbChoice.pbl_mynn.B1;
 
-        int moist_flag = 0;
-        int n_qstate = cell_data.nComp() - (NVAR_max - NMOIST_max);
+        int moist_comp = (n_qstate > 0) ? 1 : 0;
         if (n_qstate > 3) {
-            moist_flag = (n_qstate > 3) ? RhoQ4_comp : RhoQ3_comp;
-        } else if (n_qstate > 0) {
-            moist_flag = 1;
+            moist_comp = RhoQ4_comp;
         }
 
         ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -754,7 +752,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
             cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim,
                                                                   mu_turb,cellSizeInv,domain,
                                                                   pbl_mynn_B1_l,tm_arr(i,j,0),
-                                                                  moist_flag,
+                                                                  moist_comp,
                                                                   c_ext_dir_on_zlo, c_ext_dir_on_zhi,
                                                                   u_ext_dir_on_zlo, u_ext_dir_on_zhi,
                                                                   v_ext_dir_on_zlo, v_ext_dir_on_zhi,

Source/ERF.cpp

@@ -536,7 +536,8 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0)
       if (is_it_time_for_action(istep[0], time, dt_lev0, bndry_output_planes_interval, bndry_output_planes_per) &&
           time >= bndry_output_planes_start_time)
       {
-         m_w2d->write_planes(istep[0], time, vars_new);
+         bool is_moist = (micro->Get_Qstate_Size() > 0);
+         m_w2d->write_planes(istep[0], time, vars_new, is_moist);
       }
     }
 
@@ -853,7 +854,8 @@ ERF::InitData ()
 
         Real time = 0.;
         if (time >= bndry_output_planes_start_time) {
-            m_w2d->write_planes(0, time, vars_new);
+            bool is_moist = (micro->Get_Qstate_Size() > 0);
+            m_w2d->write_planes(0, time, vars_new, is_moist);
         }
     }
 
@@ -978,9 +980,12 @@ ERF::InitData ()
             }
         }
 
+        const int n_qstate = micro->Get_Qstate_Size();
+
         m_most = std::make_unique<ABLMost>(geom, use_exp_most, use_rot_most,
                                            vars_old, Theta_prim, Qv_prim, Qr_prim, z_phys_nd,
-                                           sst_lev, lmask_lev, lsm_data, lsm_flux, Hwave, Lwave, eddyDiffs_lev
+                                           sst_lev, lmask_lev, lsm_data, lsm_flux,
+                                           Hwave, Lwave, eddyDiffs_lev, n_qstate
 #ifdef ERF_USE_NETCDF
                                            ,start_bdy_time, bdy_time_interval
 #endif

Source/ERF_make_new_arrays.cpp

@@ -120,7 +120,7 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
         ncomp = vars_new[lev-1][Vars::cons].nComp();
     } else {
         int n_qstate   = micro->Get_Qstate_Size();
-        ncomp = NVAR_max - (NMOIST_max - n_qstate);
+        ncomp = NDRY + NSCALARS + n_qstate;
     }
 
     // ********************************************************************************************

Source/IO/ERF_WriteBndryPlanes.H

@@ -18,7 +18,8 @@ public:
                                amrex::Vector<amrex::Geometry>& geom);
 
     void write_planes (int t_step, amrex::Real time,
-                       amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_new);
+                       amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_new,
+                       bool is_moist);
 
 private:
 

Source/IO/ERF_WriteBndryPlanes.cpp

@@ -52,7 +52,7 @@ int WriteBndryPlanes::bndry_lev = 0;
  * @param geom Vector of Geometry containing the geometry at each level in the AMR
  */
 WriteBndryPlanes::WriteBndryPlanes (Vector<BoxArray>& grids,
-                                    Vector<Geometry>& geom): m_geom(geom)
+                                    Vector<Geometry>& geom) : m_geom(geom)
 {
     ParmParse pp("erf");
 
@@ -117,7 +117,8 @@ WriteBndryPlanes::WriteBndryPlanes (Vector<BoxArray>& grids,
  * @param vars_new Grid data for all variables across the AMR hierarchy
  */
 void WriteBndryPlanes::write_planes (const int t_step, const Real time,
-                                     Vector<Vector<MultiFab>>& vars_new)
+                                     Vector<Vector<MultiFab>>& vars_new,
+                                     bool is_moist)
 {
     BL_PROFILE("ERF::WriteBndryPlanes::write_planes");
 
@@ -143,9 +144,6 @@ void WriteBndryPlanes::write_planes (const int t_step, const Real time,
     Box target_box_shifted(IntVect(0,0,0),new_hi);
     BoxArray ba_shifted(target_box_shifted);
 
-    int n_moist_var = NMOIST_max - (S.nComp() - NVAR_max);
-    bool ismoist = (n_moist_var >= 1);
-
     for (int i = 0; i < m_var_names.size(); i++)
     {
         std::string var_name = m_var_names[i];
@@ -172,7 +170,7 @@ void WriteBndryPlanes::write_planes (const int t_step, const Real time,
             for (MFIter mfi(Temp, TilingIfNotGPU()); mfi.isValid(); ++mfi)
             {
                 const Box& bx = mfi.tilebox();
-                if (ismoist) {
+                if (is_moist) {
                     derived::erf_dermoisttemp(bx, Temp[mfi], 0, 1, S[mfi], m_geom[bndry_lev], time, nullptr, bndry_lev);
                 } else {
                     derived::erf_dertemp(bx, Temp[mfi], 0, 1, S[mfi], m_geom[bndry_lev], time, nullptr, bndry_lev);

Source/IO/Plotfile.cpp

@@ -46,14 +46,11 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>
     // since they may be in any order in the input list
     Vector<std::string> tmp_plot_names;
 
-    int n_qstate   = micro->Get_Qstate_Size();
-    int ncomp_cons = NVAR_max - (NMOIST_max - n_qstate);
-
-    for (int i = 0; i < ncomp_cons; ++i) {
+    for (int i = 0; i < cons_names.size(); ++i) {
         if ( containerHasElement(plot_var_names, cons_names[i]) ) {
-            if ( (solverChoice.moisture_type == MoistureType::SAM) || (derived_names[i] != "rhoQ4" &&
-                                                                       derived_names[i] != "rhoQ5" &&
-                                                                       derived_names[i] != "rhoQ6") )
+            if ( (solverChoice.moisture_type == MoistureType::SAM) || (cons_names[i] != "rhoQ4" &&
+                                                                       cons_names[i] != "rhoQ5" &&
+                                                                       cons_names[i] != "rhoQ6") )
             {
                 tmp_plot_names.push_back(cons_names[i]);
             } // moisture_type

Source/IndexDefines.H

@@ -8,26 +8,35 @@
  * Definition of indexing parameters
 */
 
-// Cell-centered state variables
-#define Rho_comp       0
-#define RhoTheta_comp  1
-#define RhoKE_comp     2 // for Deardorff LES Model
-#define RhoQKE_comp    3 // for MYNN or YSU PBL Model
-#define RhoScalar_comp 4
-#define RhoQ1_comp     5
-#define RhoQ2_comp     6
-#define RhoQ3_comp     7
-#define RhoQ4_comp     8
-#define RhoQ5_comp     9
-#define RhoQ6_comp     10
+// This defines the ACTUAL number of non-moisture vars =
+//      rho, rhotheta, rhoKE, rhoQKE
+#define NDRY       4
+
+// This defines the ACTUAL number of non-moisture scalar vars
+#define NSCALARS   1
+
+// This defines the MAXIMUM number of moisture vars
+#define NMOIST_max 6
 
 // This is the number of components if using moisture
 // We use this to allocate the 1d arrays of boundary condition types,
 //    but not to allocate actual solution data
-#define NVAR_max 11
+#define NVAR_max (NDRY + NSCALARS + NMOIST_max)
 
-// This defines the maximum number of moisture vars
-#define NMOIST_max 6
+// Cell-centered state variables
+#define Rho_comp       0
+#define RhoTheta_comp  (Rho_comp+1)
+#define RhoKE_comp     (Rho_comp+2) // for Deardorff LES Model
+#define RhoQKE_comp    (Rho_comp+3) // for MYNN or YSU PBL Model
+
+#define RhoScalar_comp (RhoQKE_comp+1)
+
+#define RhoQ1_comp     (RhoScalar_comp+NSCALARS)
+#define RhoQ2_comp     (RhoQ1_comp+1)
+#define RhoQ3_comp     (RhoQ1_comp+2)
+#define RhoQ4_comp     (RhoQ1_comp+3)
+#define RhoQ5_comp     (RhoQ1_comp+4)
+#define RhoQ6_comp     (RhoQ1_comp+5)
 
 // Cell-centered primitive variables
 #define PrimTheta_comp   (RhoTheta_comp -1)

Source/PBL/PBLModels.H

@@ -92,19 +92,19 @@ ComputeVerticalDerivativesPBL (int i, int j, int k,
                                amrex::Real& dthetadz,
                                amrex::Real& dudz,
                                amrex::Real& dvdz,
-                               const int moist_flag = 0)
+                               const int moist_comp = 0)
 {
     if ( k==izmax && c_ext_dir_on_zhi ) {
-        dthetadz = (1.0/3.0)*(       -Thetav(i,j,k-1,cell_data,moist_flag)
-                              - 3.0 * Thetav(i,j,k  ,cell_data,moist_flag)
-                              + 4.0 * Thetav(i,j,k+1,cell_data,moist_flag) )*dz_inv;
+        dthetadz = (1.0/3.0)*(       -Thetav(i,j,k-1,cell_data,moist_comp)
+                              - 3.0 * Thetav(i,j,k  ,cell_data,moist_comp)
+                              + 4.0 * Thetav(i,j,k+1,cell_data,moist_comp) )*dz_inv;
     } else if ( k==izmin && c_ext_dir_on_zlo ) {
-        dthetadz = (1.0/3.0)*(        Thetav(i,j,k+1,cell_data,moist_flag)
-                              + 3.0 * Thetav(i,j,k  ,cell_data,moist_flag)
-                              - 4.0 * Thetav(i,j,k-1,cell_data,moist_flag) )*dz_inv;
+        dthetadz = (1.0/3.0)*(        Thetav(i,j,k+1,cell_data,moist_comp)
+                              + 3.0 * Thetav(i,j,k  ,cell_data,moist_comp)
+                              - 4.0 * Thetav(i,j,k-1,cell_data,moist_comp) )*dz_inv;
     } else {
-        dthetadz = 0.5*(  Thetav(i,j,k+1,cell_data,moist_flag)
-                        - Thetav(i,j,k-1,cell_data,moist_flag) )*dz_inv;
+        dthetadz = 0.5*(  Thetav(i,j,k+1,cell_data,moist_comp)
+                        - Thetav(i,j,k-1,cell_data,moist_comp) )*dz_inv;
     }
 
     if ( k==izmax && u_ext_dir_on_zhi ) {
@@ -154,7 +154,7 @@ ComputeQKESourceTerms (int i, int j, int k,
                        const amrex::Box& domain,
                        amrex::Real pbl_mynn_B1_l,
                        const amrex::Real theta_mean,
-                       const int moist_flag,
+                       const int moist_comp,
                        bool c_ext_dir_on_zlo,
                        bool c_ext_dir_on_zhi,
                        bool u_ext_dir_on_zlo,
@@ -177,7 +177,7 @@ ComputeQKESourceTerms (int i, int j, int k,
                                   u_ext_dir_on_zlo, u_ext_dir_on_zhi,
                                   v_ext_dir_on_zlo, v_ext_dir_on_zhi,
                                   dthetadz, dudz, dvdz,
-                                  moist_flag);
+                                  moist_comp);
 
     // Note: Transport terms due to turbulence and pressure are included when
     // DiffusionSrcForState_* is called from ERF_slow_rhs_post.