Metgrid with SST and LANDMASK (#1280)

* WIP initialization from met_em files * WIP metgrid initialization now timesteps in serial. * WIP metgrid initialization now checks met_em global attr * WIP started writing procedures for setting lateral boundary conditions from met_em files. * WIP progress on methods for setting lateral boundary conditions from met_em files. * WIP. Initialization & forcing from met_em files runs in serial but not yet in parallel. Currently limited to met_em generated from GFS because met_em variables differ slightly depending on the origin model processed. Runs with mositure hard-coded to 0 look great. Runs with moisture do not look great and there is either a bug or error in the methodology for calculating dry density and dry pressure. * WIP cleaned up errors from prior merge with development branch updates * WIP debugging initialization and forcing from met_em files. * WIP merging and fixing merge errors from recent commits from development. * WIP more debugging for WPS met_em*.nc initialization and forcing. * Add SST and LANDMASK read and copy. * This seems to be working but dies with MOST due to 0 eddy viscosity in ghost cell. * Amrex submodule. * Unify how we initialize time and do temporal interpolation with init from wrf and metgrid. * Cmake compile with metgrid. * Remove debug print. * Template build FABs from NETCDF to allow different data types. * Fix parallel issue on CPU. Need to implement GPU option. * Minor fixes for GPU. * Parallel run and GPU compile successful. * Cmake build with metgrid. --------- Co-authored-by: wiersema1 <[email protected]> Co-authored-by: Ann Almgren <[email protected]> Co-authored-by: Aaron Lattanzi <[email protected]>
erf-model · Nov 8, 2023 · 3395bc0 · 3395bc0
1 parent 80397f6
commit 3395bc0
Show file tree

Hide file tree

Showing 41 changed files with 2,148 additions and 797 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
@@ -113,6 +112,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 */