Merge branch 'development' into dg/erf_particles

.gitmodules

@@ -7,3 +7,6 @@
 [submodule "Submodules/RRTMGP"]
 	path = Submodules/RRTMGP
 	url = https://github.com/E3SM-Project/rte-rrtmgp
+[submodule "Submodules/WW3"]
+	path = Submodules/WW3
+	url = https://github.com/erf-model/WW3

CMake/BuildERFExe.cmake

@@ -107,6 +107,7 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/ERF.cpp
        ${SRC_DIR}/ERF_make_new_arrays.cpp
        ${SRC_DIR}/ERF_make_new_level.cpp
+       ${SRC_DIR}/ERF_read_waves.cpp
        ${SRC_DIR}/ERF_Tagging.cpp
        ${SRC_DIR}/Advection/AdvectionSrcForMom.cpp
        ${SRC_DIR}/Advection/AdvectionSrcForState.cpp
@@ -139,6 +140,7 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/Initialization/ERF_init_from_metgrid.cpp
        ${SRC_DIR}/Initialization/ERF_init_uniform.cpp
        ${SRC_DIR}/Initialization/ERF_init1d.cpp
+       ${SRC_DIR}/Initialization/ERF_init_TurbPert.cpp
        ${SRC_DIR}/Initialization/ERF_input_sponge.cpp
        ${SRC_DIR}/IO/Checkpoint.cpp
        ${SRC_DIR}/IO/ERF_ReadBndryPlanes.cpp

Exec/ABL/GNUmakefile

@@ -19,7 +19,7 @@ USE_HIP  = FALSE
 USE_SYCL = FALSE
 
 # Debugging
-DEBUG = FALSE
+DEBUG = TRUE
 
 TEST = TRUE
 USE_ASSERTION = TRUE

Exec/ABL/inputs_mpmd

@@ -0,0 +1,68 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 5
+
+amrex.fpe_trap_invalid = 0
+
+fabarray.mfiter_tile_size = 1024 1024 1024
+
+# PROBLEM SIZE & GEOMETRY
+geometry.prob_extent =  191000     91000    1024
+amr.n_cell           =   191       91       4
+
+geometry.is_periodic = 1 1 0
+
+zlo.type = "NoSlipWall"
+zhi.type = "SlipWall"
+
+# TIME STEP CONTROL
+erf.fixed_dt           = 10  # fixed time step depending on grid resolution
+erf.max_step=1
+# DIAGNOSTICS & VERBOSITY
+erf.sum_interval   = 1       # timesteps between computing mass
+erf.v              = 1       # verbosity in ERF.cpp
+amr.v                = 1       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 0       # maximum level number allowed
+
+# CHECKPOINT FILES
+erf.check_file      = chk        # root name of checkpoint file
+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
+
+# SOLVER CHOICE
+erf.alpha_T = 0.0
+erf.alpha_C = 1.0
+erf.use_gravity = false
+
+erf.molec_diff_type = "None"
+erf.les_type        = "Smagorinsky"
+erf.Cs              = 0.1
+
+erf.init_type = "uniform"
+
+# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA)
+zlo.type      = "Most"
+# erf.most.z0   = 0.1  Don't set z0
+erf.most.zref = 128.0 # set 1/2 of the vertical dz: 1024 / 4 = 256 -> 256/2 = 128
+erf.most.roughness_type_sea = wave_coupled
+
+# PROBLEM PARAMETERS
+prob.rho_0 = 1.0
+prob.A_0 = 1.0
+
+prob.U_0 = 10.0
+prob.V_0 = 0.0
+prob.W_0 = 0.0
+prob.T_0 = 300.0
+
+# Higher values of perturbations lead to instability
+# Instability seems to be coming from BC
+prob.U_0_Pert_Mag = 0.08
+prob.V_0_Pert_Mag = 0.08 #
+prob.W_0_Pert_Mag = 0.0
+

Exec/DevTests/MovingTerrain/prob.cpp

@@ -155,47 +155,57 @@ Problem::init_custom_terrain (const Geometry& geom,
                               MultiFab& z_phys_nd,
                               const Real& time)
 {
-    // Domain cell size and real bounds
-    auto dx = geom.CellSizeArray();
-
-    // Domain valid box (z_nd is nodal)
-    const amrex::Box& domain = geom.Domain();
-    int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
-    int domlo_z = domain.smallEnd(2);
-
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    // Populate bottom plane
-    int k0 = domlo_z;
-
-    Real Ampl        = parms.Ampl;
-    Real wavelength  = 100.;
-    Real kp          = 2.0 * PI / wavelength;
-    Real g           = CONST_GRAV;
-    Real omega       = std::sqrt(g * kp);
-
-    for (MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Clip indices for ghost-cells
-            int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
-
-            // Location of nodes
-            Real x = ii  * dx[0];
-
-            // Wave height
-            Real height = Ampl * std::sin(kp * x - omega * time);
-
-            // Populate terrain height
-            z_arr(i,j,k0) = height;
-        });
+    // Check if a valid csv file exists for the terrain
+    std::string fname;
+    amrex::ParmParse pp("erf");
+    auto valid_fname = pp.query("terrain_file_name",fname);
+    if (valid_fname) {
+        this->read_custom_terrain(fname,geom,z_phys_nd,time);
+    } else {
+
+        // Domain cell size and real bounds
+        auto dx = geom.CellSizeArray();
+
+        // Domain valid box (z_nd is nodal)
+        const amrex::Box& domain = geom.Domain();
+        int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
+        int domlo_z = domain.smallEnd(2);
+
+        // Number of ghost cells
+        int ngrow = z_phys_nd.nGrow();
+
+        // Populate bottom plane
+        int k0 = domlo_z;
+
+        Real Ampl        = parms.Ampl;
+        Real wavelength  = 100.;
+        Real kp          = 2.0 * PI / wavelength;
+        Real g           = CONST_GRAV;
+        Real omega       = std::sqrt(g * kp);
+
+        for (MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
+        {
+            // Grown box with no z range
+            amrex::Box xybx = mfi.growntilebox(ngrow);
+            xybx.setRange(2,0);
+
+            amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
+
+            ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int)
+            {
+
+                // Clip indices for ghost-cells
+                int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
+
+                // Location of nodes
+                Real x = ii  * dx[0];
+
+                // Wave height
+                Real height = Ampl * std::sin(kp * x - omega * time);
+
+                // Populate terrain height
+                z_arr(i,j,k0) = height;
+            });
+        }
     }
 }

Exec/Make.ERF

@@ -155,6 +155,10 @@ ifeq ($(USE_WINDFARM), TRUE)
 	INCLUDE_LOCATIONS += $(ERF_WINDFARM_EWP_DIR)
 endif
 
+ifeq ($(USE_WW3_COUPLING), TRUE)
+	DEFINES += -DERF_USE_WW3_COUPLING
+endif
+
 ERF_LSM_DIR = $(ERF_SOURCE_DIR)/LandSurfaceModel
 include $(ERF_LSM_DIR)/Make.package
 VPATH_LOCATIONS   += $(ERF_LSM_DIR)

Exec/RegTests/ParticlesOverWoA/prob.cpp

@@ -126,54 +126,64 @@ void
 Problem::init_custom_terrain(
     const Geometry& geom,
     MultiFab& z_phys_nd,
-    const Real& /*time*/)
+    const Real& time)
 {
-    // Domain cell size and real bounds
-    auto dx = geom.CellSizeArray();
-    auto ProbLoArr = geom.ProbLoArray();
-    auto ProbHiArr = geom.ProbHiArray();
-
-    // Domain valid box (z_nd is nodal)
-    const amrex::Box& domain = geom.Domain();
-    int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
-    // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
-    int domlo_z = domain.smallEnd(2);
-
-    // User function parameters
-    Real a    = 0.5;
-    Real num  = 8 * a * a * a;
-    Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]);
-    // Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]);
-
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    // Populate bottom plane
-    int k0 = domlo_z;
-
-    for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Clip indices for ghost-cells
-            int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
-            // int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
-
-            // Location of nodes
-            Real x = (ii  * dx[0] - xcen);
-            // Real y = (jj  * dx[1] - ycen);
-
-            // WoA Hill in x-direction
-            Real height = num / (x*x + 4 * a * a);
-
-            // Populate terrain height
-            z_arr(i,j,k0) = height;
-        });
+    // Check if a valid csv file exists for the terrain
+    std::string fname;
+    amrex::ParmParse pp("erf");
+    auto valid_fname = pp.query("terrain_file_name",fname);
+    if (valid_fname) {
+        this->read_custom_terrain(fname,geom,z_phys_nd,time);
+    } else {
+
+        // Domain cell size and real bounds
+        auto dx = geom.CellSizeArray();
+        auto ProbLoArr = geom.ProbLoArray();
+        auto ProbHiArr = geom.ProbHiArray();
+
+        // Domain valid box (z_nd is nodal)
+        const amrex::Box& domain = geom.Domain();
+        int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
+        // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
+        int domlo_z = domain.smallEnd(2);
+
+        // User function parameters
+        Real a    = 0.5;
+        Real num  = 8 * a * a * a;
+        Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]);
+        // Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]);
+
+        // Number of ghost cells
+        int ngrow = z_phys_nd.nGrow();
+
+        // Populate bottom plane
+        int k0 = domlo_z;
+
+        for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
+        {
+            // Grown box with no z range
+            amrex::Box xybx = mfi.growntilebox(ngrow);
+            xybx.setRange(2,0);
+
+            amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
+
+            ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int)
+            {
+
+                // Clip indices for ghost-cells
+                int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
+                // int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
+
+                // Location of nodes
+                Real x = (ii  * dx[0] - xcen);
+                // Real y = (jj  * dx[1] - ycen);
+
+                // WoA Hill in x-direction
+                Real height = num / (x*x + 4 * a * a);
+
+                // Populate terrain height
+                z_arr(i,j,k0) = height;
+            });
+        }
     }
 }