Merge branch 'development' into xyuan/rad_develop

ewquon · Nov 9, 2023 · 6989ae2 · 6989ae2
2 parents 138124c + ac2e922
commit 6989ae2
Show file tree

Hide file tree

Showing 42 changed files with 2,270 additions and 938 deletions.
diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake
@@ -45,7 +45,6 @@ function(build_erf_lib erf_lib_name)
 
   if(ERF_ENABLE_NETCDF)
     target_sources(${erf_lib_name} PRIVATE
-                   ${SRC_DIR}/IO/NCBuildFABs.cpp
                    ${SRC_DIR}/IO/NCInterface.cpp
                    ${SRC_DIR}/IO/NCPlotFile.cpp
                    ${SRC_DIR}/IO/NCCheckpoint.cpp
@@ -114,6 +113,7 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/BoundaryConditions/BoundaryConditions_zvel.cpp
        ${SRC_DIR}/BoundaryConditions/BoundaryConditions_bndryreg.cpp
        ${SRC_DIR}/BoundaryConditions/BoundaryConditions_wrfbdy.cpp
+       ${SRC_DIR}/BoundaryConditions/BoundaryConditions_metgrid.cpp
        ${SRC_DIR}/BoundaryConditions/ERF_FillPatch.cpp
        ${SRC_DIR}/BoundaryConditions/ERF_FillPatcher.cpp
        ${SRC_DIR}/BoundaryConditions/ERF_PhysBCFunct.cpp

diff --git a/Exec/CMakeLists.txt b/Exec/CMakeLists.txt
@@ -24,4 +24,5 @@ else ()
   add_subdirectory(DevTests/MovingTerrain)
   add_subdirectory(DevTests/ParticlesOverWoA)
   add_subdirectory(DevTests/MiguelDev)
+  add_subdirectory(DevTests/MetGrid)
 endif()
diff --git a/Exec/DevTests/MetGrid/CMakeLists.txt b/Exec/DevTests/MetGrid/CMakeLists.txt
@@ -0,0 +1,12 @@
+set(erf_exe_name erf_metgrid_dev)
+
+add_executable(${erf_exe_name} "")
+target_sources(${erf_exe_name}
+   PRIVATE
+     prob.cpp
+)
+
+target_include_directories(${erf_exe_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
+
+include(${CMAKE_SOURCE_DIR}/CMake/BuildERFExe.cmake)
+build_erf_exe(${erf_exe_name})
diff --git a/Exec/DevTests/MetGrid/GNUmakefile b/Exec/DevTests/MetGrid/GNUmakefile
@@ -19,14 +19,19 @@ USE_HIP  = FALSE
 USE_SYCL = FALSE
 
 # Debugging
-DEBUG = FALSE
+#DEBUG = FALSE
 DEBUG = TRUE
 
 TEST = TRUE
 USE_ASSERTION = TRUE
 
 USE_NETCDF = TRUE
 
+#USE_MOISTURE = FALSE
+USE_MOISTURE = TRUE
+
+#USE_WARM_NO_PRECIP = TRUE
+
 # GNU Make
 Bpack := ./Make.package
 Blocs := .

diff --git a/Exec/DevTests/MetGrid/inputs b/Exec/DevTests/MetGrid/inputs
@@ -1,11 +1,13 @@
 # ------------------  INPUTS TO MAIN PROGRAM  -------------------
-max_step = 1
+max_step = 100
 
 amrex.fpe_trap_invalid = 1
+amrex.fpe_trap_zero = 1
+amrex.fpe_trap_overflow = 1
 
 # PROBLEM SIZE & GEOMETRY
-geometry.prob_extent =  10600 5000 5000
-amr.n_cell           =   140 80 60
+geometry.prob_extent =  28000 16000 8000
+amr.n_cell           =   140 80 100
 
 geometry.is_periodic = 0 0 0
 
@@ -17,9 +19,10 @@ zlo.type = "NoSlipWall"
 zhi.type = "SlipWall"
 
 # TIME STEP CONTROL
-erf.fixed_dt           = 0.1  # fixed time step depending on grid resolution
+erf.fixed_dt           = 1.0  # fixed time step depending on grid resolution
 
 erf.use_terrain = true
+erf.terrain_smoothing = 2
 
 # DIAGNOSTICS & VERBOSITY
 erf.sum_interval   = -1       # timesteps between computing mass
@@ -36,22 +39,26 @@ erf.check_int       = 100        # number of timesteps between checkpoints
 
 # PLOTFILES
 erf.plot_file_1     = plt        # prefix of plotfile name
-erf.plot_int_1      = 10         # number of timesteps between plotfiles
-erf.plot_vars_1     = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys
+erf.plot_int_1      = 1          # number of timesteps between plotfiles
+erf.plot_vars_1     = qv Rhoqv density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres KE QKE
 
 # SOLVER CHOICE
-erf.alpha_T = 0.0
+erf.alpha_T = 1.0
 erf.alpha_C = 1.0
 erf.use_gravity = true
 
 erf.molec_diff_type = "None"
 erf.les_type        = "Smagorinsky"
+#erf.les_type        = "Deardorff"
 erf.Cs              = 0.1
 
+#erf.terrain_z_levels = 0 130 354 583 816 1054 1549 2068 2615 3193 3803 4450 5142 5892 6709 7603 8591 9702 10967 12442 14230 16610 18711 20752 22133 23960 26579 28493 31236 33699 36068 40000
+
 # INITIALIZATION WITH ATM DATA
-erf.init_type    = "metgrid"
-erf.nc_init_file_0 = "met_em_d01.nc"
-#erf.nc_bdy_file = ""
+erf.metgrid_bdy_width = 5
+erf.metgrid_bdy_set_width = 1
+erf.init_type      = "metgrid"
+erf.nc_init_file_0 = "met_em.d01.2022-06-18_00:00:00.nc" "met_em.d01.2022-06-18_06:00:00.nc" "met_em.d01.2022-06-18_12:00:00.nc" "met_em.d01.2022-06-18_18:00:00.nc"
 
 #There will be no OpenMP tiling
 fabarray.mfiter_tile_size = 1024 1024 1024
diff --git a/Exec/DevTests/MetGrid/prob.H b/Exec/DevTests/MetGrid/prob.H
@@ -12,7 +12,7 @@ struct ProbParm : ProbParmDefaults {
 class Problem : public ProblemBase
 {
 public:
-    Problem();
+    Problem(const amrex::Real* problo, const amrex::Real* probhi);
 
 protected:
     std::string name() override { return "Metgrid"; }

diff --git a/Exec/DevTests/MetGrid/prob.cpp b/Exec/DevTests/MetGrid/prob.cpp
@@ -8,5 +8,5 @@ amrex_probinit(const amrex_real* problo, const amrex_real* probhi)
     return std::make_unique<Problem>(problo, probhi);
 }
 
-Problem::Problem()
+Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi)
 {}
diff --git a/Exec/LLJ/GNUmakefile b/Exec/LLJ/GNUmakefile
@@ -25,7 +25,7 @@ TEST = TRUE
 USE_ASSERTION = TRUE
 
 USE_NETCDF  = TRUE
-USE_HDF5    = TRUE
+#USE_HDF5    = TRUE
 
 # GNU Make
 Bpack := ./Make.package

diff --git a/Exec/RegTests/WPS_Test/GNUmakefile b/Exec/RegTests/WPS_Test/GNUmakefile
@@ -27,7 +27,7 @@ USE_ASSERTION = TRUE
 USE_MOISTURE = TRUE
 
 USE_NETCDF  = TRUE
-USE_HDF5    = TRUE
+#USE_HDF5    = TRUE
 
 # GNU Make
 Bpack := ./Make.package

diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H
@@ -33,8 +33,14 @@ public:
                       amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_old,
                       amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Theta_prim,
                       amrex::Vector<std::unique_ptr<amrex::MultiFab>>& z_phys_nd,
-                      int ng_in)
-      : m_geom(geom), m_ma(geom,vars_old,Theta_prim,z_phys_nd)
+                      amrex::Vector<amrex::Vector<std::unique_ptr<amrex::MultiFab>>>& sst_lev,
+                      amrex::Vector<amrex::Vector<std::unique_ptr<amrex::iMultiFab>>>& lmask_lev,
+                      amrex::Real start_bdy_time = 0.0,
+                      amrex::Real bdy_time_interval = 0.0)
+    : m_start_bdy_time(start_bdy_time),
+      m_bdy_time_interval(bdy_time_interval),
+      m_geom(geom),
+      m_ma(geom,vars_old,Theta_prim,z_phys_nd)
     {
         amrex::ParmParse pp("erf");
         pp.query("most.z0", z0_const);
@@ -86,116 +92,137 @@ public:
             amrex::Abort("Undefined MOST roughness type!");
         }
 
+        // Size the MOST params for all levels
         int nlevs = m_geom.size();
         z_0.resize(nlevs);
         u_star.resize(nlevs);
         t_star.resize(nlevs);
         t_surf.resize(nlevs);
         olen.resize(nlevs);
 
-        for (int lev = 0; lev < nlevs; lev++)
-        {
-            // Z0 heights
+        // Get pointers to SST and LANDMASK data
+        m_sst_lev.resize(nlevs);
+        m_lmask_lev.resize(nlevs);
+        for (int lev(0); lev<nlevs; ++lev) {
+            int nt_tot = sst_lev[lev].size();
+            m_sst_lev[lev].resize(nt_tot);
+            m_lmask_lev[lev].resize(nt_tot);
+            for (int nt(0); nt<nt_tot; ++nt) {
+                m_sst_lev[lev][nt]   = sst_lev[lev][nt].get();
+                m_lmask_lev[lev][nt] = lmask_lev[lev][nt].get();
+            }
+        }
+
+        for (int lev = 0; lev < nlevs; lev++) {
+            // Attributes for MFs and FABs
+            //--------------------------------------------------------
+            auto& mf = vars_old[lev][Vars::cons];
+            // Create a 2D ba, dm, & ghost cells
+            amrex::BoxArray ba  = mf.boxArray();
+            amrex::BoxList bl2d = ba.boxList();
+            for (auto& b : bl2d) {
+                b.setRange(2,0);
+            }
+            amrex::BoxArray ba2d(std::move(bl2d));
+            const amrex::DistributionMapping& dm = mf.DistributionMap();
+            const int ncomp   = 1;
+            amrex::IntVect ng = mf.nGrowVect(); ng[2]=0;
+
+            // Z0 heights FAB
             //--------------------------------------------------------
-            amrex::Box bx = amrex::grow(m_geom[lev].Domain(),ng_in);
+            amrex::Box bx = amrex::grow(m_geom[lev].Domain(),ng);
             bx.setSmall(2,0);
             bx.setBig(2,0);
             z_0[lev].resize(bx,1);
             z_0[lev].setVal<amrex::RunOn::Device>(z0_const);
 
             // 2D MFs for U*, T*, T_surf
             //--------------------------------------------------------
-            { // CC vars
-                auto& mf = vars_old[lev][Vars::cons];
-                // Create a 2D ba, dm, & ghost cells
-                amrex::BoxArray ba  = mf.boxArray();
-                amrex::BoxList bl2d = ba.boxList();
-                for (auto& b : bl2d) {
-                    b.setRange(2,0);
-                }
-                amrex::BoxArray ba2d(std::move(bl2d));
-                const amrex::DistributionMapping& dm = mf.DistributionMap();
-                const int ncomp   = 1;
-                amrex::IntVect ng = mf.nGrowVect(); ng[2]=0;
-
-                u_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
-                u_star[lev]->setVal(1.E34);
-
-                t_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
-                t_star[lev]->setVal(1.E34);
-
-                olen[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
-                olen[lev]->setVal(1.E34);
-
-                t_surf[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
-                if (erf_st) {
-                   t_surf[lev]->setVal(surf_temp);
-                } else {
-                    t_surf[lev]->setVal(0.0);
-                }
+            u_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
+            u_star[lev]->setVal(1.E34);
+
+            t_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
+            t_star[lev]->setVal(1.E34);
+
+            olen[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
+            olen[lev]->setVal(1.E34);
+
+            t_surf[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
+
+            if (m_sst_lev[lev][0]) { // Valid SST data at t==0
+                theta_type = ThetaCalcType::SURFACE_TEMPERATURE;
+                amrex::MultiFab::Copy(*(t_surf[lev]), *(m_sst_lev[lev][0]), 0, 0, 1, ng);
+            } else if (erf_st) {     // Constant temp
+                t_surf[lev]->setVal(surf_temp);
+            } else {
+                t_surf[lev]->setVal(0.0);
             }
         }// lev
     }
 
     void
-    update_surf_temp (amrex::Real cur_time)
-    {
-        if (surf_heating_rate != 0) {
-            int nlevs = m_geom.size();
-            for (int lev = 0; lev < nlevs; lev++)
-            {
-               t_surf[lev]->setVal(surf_temp + surf_heating_rate*cur_time);
-               amrex::Print() << "Surface temp at t=" << cur_time
-                   << ": "
-                   << surf_temp + surf_heating_rate*cur_time
-                   << std::endl;
-            }
-        }
-    }
+    update_fluxes (const int& lev,
+                   const amrex::Real& time,
+                   int max_iters = 25);
+
+    template <typename FluxIter>
+    void
+    compute_fluxes (const int& lev,
+                    const int& max_iters,
+                    const FluxIter& most_flux);
 
     void
-    impose_most_bcs (const int lev,
+    impose_most_bcs (const int& lev,
                      const amrex::Vector<amrex::MultiFab*>& mfs,
                      amrex::MultiFab* eddyDiffs,
                      amrex::MultiFab* z_phys);
 
     template<typename FluxCalc>
     void
-    compute_most_bcs (const int lev,
+    compute_most_bcs (const int& lev,
                       const amrex::Vector<amrex::MultiFab*>& mfs,
                       amrex::MultiFab* eddyDiffs,
                       amrex::MultiFab* z_phys,
                       const amrex::Real& dz_no_terrain,
                       const FluxCalc& flux_comp);
 
     void
-    update_fluxes (int lev,
-                   amrex::Real cur_time,
-                   int max_iters = 25);
+    time_interp_tsurf(const int& lev,
+                      const amrex::Real& time);
 
-    template <typename FluxIter>
     void
-    compute_fluxes (const int& lev,
-                    const int& max_iters,
-                    const FluxIter& most_flux);
+    update_surf_temp (const amrex::Real& time)
+    {
+        if (surf_heating_rate != 0) {
+            int nlevs = m_geom.size();
+            for (int lev = 0; lev < nlevs; lev++)
+            {
+               t_surf[lev]->setVal(surf_temp + surf_heating_rate*time);
+               amrex::Print() << "Surface temp at t=" << time
+                   << ": "
+                   << surf_temp + surf_heating_rate*time
+                   << std::endl;
+            }
+        }
+    }
 
     void
-    update_mac_ptrs (int lev,
+    update_mac_ptrs (const int& lev,
                      amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_old,
                      amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Theta_prim)
     { m_ma.update_field_ptrs(lev,vars_old,Theta_prim); }
 
     const amrex::MultiFab*
-    get_u_star (int lev) { return u_star[lev].get(); }
+    get_u_star (const int& lev) { return u_star[lev].get(); }
 
     const amrex::MultiFab*
-    get_t_star (int lev) { return t_star[lev].get(); }
+    get_t_star (const int& lev) { return t_star[lev].get(); }
 
     const amrex::MultiFab*
-    get_olen (int lev) { return olen[lev].get(); }
+    get_olen (const int& lev) { return olen[lev].get(); }
 
     const amrex::MultiFab*
-    get_mac_avg (int lev, int comp) { return m_ma.get_average(lev,comp); }
+    get_mac_avg (const int& lev, int comp) { return m_ma.get_average(lev,comp); }
 
     enum struct FluxCalcType {
         MOENG = 0,      ///< Moeng functional form
@@ -225,6 +252,8 @@ private:
     amrex::Real surf_temp_flux{0};
     amrex::Real cnk_a{0.0185};
     amrex::Real depth{30.0};
+    amrex::Real m_start_bdy_time;
+    amrex::Real m_bdy_time_interval;
     amrex::Vector<amrex::Geometry>  m_geom;
     amrex::Vector<amrex::FArrayBox> z_0;
 
@@ -233,6 +262,9 @@ private:
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> t_star;
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> olen;
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> t_surf;
+
+    amrex::Vector<amrex::Vector<amrex::MultiFab*>>  m_sst_lev;
+    amrex::Vector<amrex::Vector<amrex::iMultiFab*>> m_lmask_lev;
 };
 
 #endif /* ABLMOST_H */