From c1ed77876337fe59069ade2bd98682f01d334dfe Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Fri, 2 Feb 2024 15:37:38 -0700 Subject: [PATCH] Fix flat top with wrfinput and allow for top extrapolation with metgrid. (#1423) --- Exec/DevTests/MetGrid/inputs_flowmas_metgrid | 71 +++++++++++++++++++ ..._fitch_flowmas => inputs_flowmas_wrfinput} | 17 ++--- .../Initialization/ERF_init_from_wrfinput.cpp | 60 +++++++++------- Source/Initialization/Metgrid_utils.H | 42 +++++++++-- 4 files changed, 152 insertions(+), 38 deletions(-) create mode 100644 Exec/DevTests/MetGrid/inputs_flowmas_metgrid rename Exec/DevTests/MetGrid/{inputs_fitch_flowmas => inputs_flowmas_wrfinput} (76%) diff --git a/Exec/DevTests/MetGrid/inputs_flowmas_metgrid b/Exec/DevTests/MetGrid/inputs_flowmas_metgrid new file mode 100644 index 000000000..889dbfa5a --- /dev/null +++ b/Exec/DevTests/MetGrid/inputs_flowmas_metgrid @@ -0,0 +1,71 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 1000 + +amrex.fpe_trap_invalid = 1 +amrex.fpe_trap_zero = 1 +amrex.fpe_trap_overflow = 1 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 447000 387000 20500 +amr.n_cell = 149 129 40 + +geometry.is_periodic = 0 0 0 + +xlo.type = "Outflow" +xhi.type = "Outflow" +ylo.type = "Outflow" +yhi.type = "Outflow" +zlo.type = "NoSlipWall" +zhi.type = "SlipWall" + +# TIME STEP CONTROL +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 +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 = 25 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 25 # number of timesteps between plotfiles +erf.plot_vars_1 = density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse dens_hse pert_pres pert_dens rhoQ1 rhoQ2 rhoQ3 qv qc + +# SOLVER CHOICE +erf.alpha_T = 1.0 +erf.alpha_C = 1.0 +erf.use_gravity = true + +erf.molec_diff_type = "None" +erf.les_type = "Smagorinsky" +erf.Cs = 0.1 + +erf.moisture_model = "Kessler" + +erf.use_terrain = true +erf.terrain_smoothing = 0 +erf.terrain_z_levels = 0 53.25373067 116.3299934 190.8583428 278.6665766 381.7196783 502.1312033 642.3454199 \ + 805.2193582 992.889535 1207.388303 1450.761988 1724.488474 2029.43012 2365.936136 \ + 2733.669456 3131.532666 3557.60783 4009.253851 4483.532727 4978.420557 5494.911024 \ + 6034.232547 6597.033978 7183.139415 7792.547066 8425.074263 9080.62066 9761.028885 \ + 10467.85208 11200.42746 11958.47509 12753.28089 13605.36124 14516.05106 15467.59286 \ + 16453.44037 17477.73597 18503.25481 19533.31103 20500.0 + +# INITIALIZATION WITH ATM DATA +erf.metgrid_bdy_width = 5 +erf.metgrid_bdy_set_width = 1 +erf.init_type = "metgrid" +erf.nc_init_file_0 = "met_em.d01.2015-04-01_00_00_00.nc" "met_em.d01.2015-04-01_08_00_00.nc" + +#There will be no OpenMP tiling +fabarray.mfiter_tile_size = 1024 1024 1024 diff --git a/Exec/DevTests/MetGrid/inputs_fitch_flowmas b/Exec/DevTests/MetGrid/inputs_flowmas_wrfinput similarity index 76% rename from Exec/DevTests/MetGrid/inputs_fitch_flowmas rename to Exec/DevTests/MetGrid/inputs_flowmas_wrfinput index 7e4ebc6e5..c46c95764 100644 --- a/Exec/DevTests/MetGrid/inputs_fitch_flowmas +++ b/Exec/DevTests/MetGrid/inputs_flowmas_wrfinput @@ -6,7 +6,7 @@ amrex.fpe_trap_zero = 1 amrex.fpe_trap_overflow = 1 # PROBLEM SIZE & GEOMETRY -geometry.prob_extent = 447000 387000 15000 +geometry.prob_extent = 447000 387000 20500 amr.n_cell = 149 129 40 geometry.is_periodic = 0 0 0 @@ -38,8 +38,8 @@ erf.check_int = 25 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name -erf.plot_int_1 = 25 # number of timesteps between plotfiles -erf.plot_vars_1 = density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres +erf.plot_int_1 = 1 # number of timesteps between plotfiles +erf.plot_vars_1 = density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse dens_hse pert_pres pert_dens rhoQ1 rhoQ2 rhoQ3 qv qc # SOLVER CHOICE erf.alpha_T = 1.0 @@ -50,13 +50,14 @@ erf.molec_diff_type = "None" erf.les_type = "Smagorinsky" erf.Cs = 0.1 -erf.moisture_model = "NullMoist" +erf.moisture_model = "Kessler" # INITIALIZATION WITH ATM DATA -erf.metgrid_bdy_width = 7 -erf.metgrid_bdy_set_width = 1 -erf.init_type = "metgrid" -erf.nc_init_file_0 = "met_em.d01.2015-04-01_00_00_00.nc" "met_em.d01.2015-04-01_08_00_00.nc" +erf.wrfbdy_width = 5 +erf.wrfbdy_set_width = 1 +erf.init_type = "real" +erf.nc_init_file_0 = "wrfinput_d01" +erf.nc_bdy_file = "wrfbdy_d01" #There will be no OpenMP tiling fabarray.mfiter_tile_size = 1024 1024 1024 diff --git a/Source/Initialization/ERF_init_from_wrfinput.cpp b/Source/Initialization/ERF_init_from_wrfinput.cpp index 780d4633f..0da111d0d 100644 --- a/Source/Initialization/ERF_init_from_wrfinput.cpp +++ b/Source/Initialization/ERF_init_from_wrfinput.cpp @@ -64,7 +64,8 @@ init_state_from_wrfinput (int lev, FArrayBox& state_fab, const Vector& NC_zvel_fab, const Vector& NC_rho_fab, const Vector& NC_rhotheta_fab, - MoistureType moisture_type); + MoistureType moisture_type, + const int& n_qstate); void init_msfs_from_wrfinput (int lev, FArrayBox& msfu_fab, @@ -73,7 +74,8 @@ init_msfs_from_wrfinput (int lev, FArrayBox& msfu_fab, const Vector& NC_MSFV_fab, const Vector& NC_MSFM_fab); void -init_terrain_from_wrfinput (int lev, const Box& domain, FArrayBox& z_phys, +init_terrain_from_wrfinput (int lev, const Real& z_top, + const Box& domain, FArrayBox& z_phys, const Vector& NC_PH_fab, const Vector& NC_PHB_fab); @@ -132,7 +134,7 @@ ERF::init_from_wrfinput (int lev) } auto& lev_new = vars_new[lev]; - + int n_qstate = micro.Get_Qstate_Size(); #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif @@ -149,7 +151,8 @@ ERF::init_from_wrfinput (int lev) init_state_from_wrfinput(lev, cons_fab, xvel_fab, yvel_fab, zvel_fab, NC_QVAPOR_fab, NC_QCLOUD_fab, NC_QRAIN_fab, NC_xvel_fab, NC_yvel_fab, NC_zvel_fab, - NC_rho_fab, NC_rhoth_fab, solverChoice.moisture_type); + NC_rho_fab, NC_rhoth_fab, + solverChoice.moisture_type, n_qstate); } // mf #ifdef _OPENMP @@ -168,14 +171,14 @@ ERF::init_from_wrfinput (int lev) } // mf const Box& domain = geom[lev].Domain(); - + const Real& z_top = geom[lev].ProbHi(2); if (solverChoice.use_terrain) { std::unique_ptr& z_phys = z_phys_nd[lev]; for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { FArrayBox& z_phys_nd_fab = (*z_phys)[mfi]; - init_terrain_from_wrfinput(lev, domain, z_phys_nd_fab, NC_PH_fab, NC_PHB_fab); + init_terrain_from_wrfinput(lev, z_top, domain, z_phys_nd_fab, NC_PH_fab, NC_PHB_fab); } // mf make_J (geom[lev],*z_phys_nd[lev],* detJ_cc[lev]); @@ -265,7 +268,8 @@ init_state_from_wrfinput (int lev, const Vector& NC_zvel_fab, const Vector& NC_rho_fab, const Vector& NC_rhotheta_fab, - MoistureType moisture_type) + MoistureType moisture_type, + const int& n_qstate) { int nboxes = NC_xvel_fab.size(); for (int idx = 0; idx < nboxes; idx++) @@ -299,8 +303,10 @@ init_state_from_wrfinput (int lev, state_fab.template copy(NC_QCLOUD_fab[idx], 0, RhoQ2_comp, 1); state_fab.template mult(NC_rho_fab[idx] , 0, RhoQ2_comp, 1); - state_fab.template copy(NC_QRAIN_fab[idx], 0, RhoQ3_comp, 1); - state_fab.template mult(NC_rho_fab[idx] , 0, RhoQ3_comp, 1); + if (n_qstate >= 3) { + state_fab.template copy(NC_QRAIN_fab[idx], 0, RhoQ3_comp, 1); + state_fab.template mult(NC_rho_fab[idx] , 0, RhoQ3_comp, 1); + } } } // idx } @@ -389,7 +395,8 @@ init_base_state_from_wrfinput (int lev, const Box& valid_bx, const Real l_rdOcp, * @param NC_PHB_fab Vector of FArrayBox objects storing WRF terrain coordinate data (PHB) */ void -init_terrain_from_wrfinput (int lev, const Box& domain, FArrayBox& z_phys, +init_terrain_from_wrfinput (int lev, const Real& z_top, + const Box& domain, FArrayBox& z_phys, const Vector& NC_PH_fab, const Vector& NC_PHB_fab) { @@ -440,27 +447,28 @@ init_terrain_from_wrfinput (int lev, const Box& domain, FArrayBox& z_phys, nc_ph_arr (ii,jj-1,khi-1) + nc_ph_arr (ii-1,jj-1,khi-1) + nc_phb_arr(ii,jj ,khi-1) + nc_phb_arr(ii-1,jj ,khi-1) + nc_phb_arr(ii,jj-1,khi-1) + nc_phb_arr(ii-1,jj-1,khi-1) ) / CONST_GRAV; - z_arr(i, j, k) = 2.0 * z_khi - z_khim1; - } else { + z_arr(i, j, k) = 2.0 * z_top - z_khim1; + } else if (k == khi) { + z_arr(i, j, k) = z_top; + } else { z_arr(i, j, k) = 0.25 * ( nc_ph_arr (ii,jj ,k) + nc_ph_arr (ii-1,jj ,k) + nc_ph_arr (ii,jj-1,k) + nc_ph_arr (ii-1,jj-1,k) + nc_phb_arr(ii,jj ,k) + nc_phb_arr(ii-1,jj ,k) + nc_phb_arr(ii,jj-1,k) + nc_phb_arr(ii-1,jj-1,k) ) / CONST_GRAV; - - // - // Fill values outside the domain on the fly (we will need these to make detJ in ghost cells) - // - if (i == domlo.x) z_arr(i-1,j,k) = z_arr(i,j,k); - if (i == domhi.x) z_arr(i+1,j,k) = z_arr(i,j,k); - if (j == domlo.y) z_arr(i,j-1,k) = z_arr(i,j,k); - if (j == domhi.y) z_arr(i,j+1,k) = z_arr(i,j,k); - - if (i == domlo.x && j == domlo.y) z_arr(i-1,j-1,k) = z_arr(i,j,k); - if (i == domlo.x && j == domhi.y) z_arr(i-1,j+1,k) = z_arr(i,j,k); - if (i == domhi.x && j == domlo.y) z_arr(i+1,j-1,k) = z_arr(i,j,k); - if (i == domhi.x && j == domhi.y) z_arr(i+1,j+1,k) = z_arr(i,j,k); - } // k + + // + // Fill values outside the domain on the fly (we will need these to make detJ in ghost cells) + // + if (i == domlo.x) z_arr(i-1,j,k) = z_arr(i,j,k); + if (i == domhi.x) z_arr(i+1,j,k) = z_arr(i,j,k); + if (j == domlo.y) z_arr(i,j-1,k) = z_arr(i,j,k); + if (j == domhi.y) z_arr(i,j+1,k) = z_arr(i,j,k); + + if (i == domlo.x && j == domlo.y) z_arr(i-1,j-1,k) = z_arr(i,j,k); + if (i == domlo.x && j == domhi.y) z_arr(i-1,j+1,k) = z_arr(i,j,k); + if (i == domhi.x && j == domlo.y) z_arr(i+1,j-1,k) = z_arr(i,j,k); + if (i == domhi.x && j == domhi.y) z_arr(i+1,j+1,k) = z_arr(i,j,k); }); } // idx } diff --git a/Source/Initialization/Metgrid_utils.H b/Source/Initialization/Metgrid_utils.H index ad96e9b09..3c19125e4 100644 --- a/Source/Initialization/Metgrid_utils.H +++ b/Source/Initialization/Metgrid_utils.H @@ -189,8 +189,8 @@ interpolate_column_metgrid (const int& i, amrex::Real z0, z1; int klow = -1; int khi0 = -1; - amrex::Real dzlow = 10000.0; - amrex::Real dzhi0 = -10000.0; + amrex::Real dzlow = 1.0e12; + amrex::Real dzhi0 = -1.0e12; for (int kk = 0; kk < kmax_orig; kk++) { amrex::Real orig_z_stag = 0.0; if (stag == 'M') { @@ -218,6 +218,7 @@ interpolate_column_metgrid (const int& i, orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i,j-1,kk)); } } + amrex::Real dz = z - orig_z_stag; if ((dz < 0.0) && (dz > dzhi0)) { dzhi0 = dz; @@ -231,10 +232,10 @@ interpolate_column_metgrid (const int& i, } } // kk + // extrapolate below the bottom surface if (klow == -1) { - // extrapolate int khi1 = -1; - amrex::Real dzhi1 = -10000.0; + amrex::Real dzhi1 = -1.0e12; for (int kk = 0; kk < kmax_orig; kk++) { amrex::Real orig_z_stag = 0.0; if (stag == 'M') { @@ -272,6 +273,39 @@ interpolate_column_metgrid (const int& i, amrex::Real y0 = orig_data(i,j,khi0,src_comp); amrex::Real y1 = orig_data(i,j,khi1,src_comp); return ( y0-(y1-y0)/(z1-z0)*(z0-z) ); + + // Extrapolate above the top surface + } else if (khi0 == -1) { + khi0 = klow - 1; + int khi1 = klow; + if (stag == 'M') { + z0 = orig_z(i,j,khi0); + } + else if (stag == 'X') { + if (i == 0) { + z0 = orig_z(i,j,khi0); + } + else if (i == imax_orig) { + z0 = orig_z(imax_orig-1,j,khi0); + } + else { + z0 = 0.5*(orig_z(i,j,khi0)+orig_z(i-1,j,khi0)); + } + } + else if (stag == 'Y') { + if (j == 0) { + z0 = orig_z(i,j,khi0); + } + else if (j == jmax_orig) { + z0 = orig_z(i,jmax_orig-1,khi0); + } + else { + z0 = 0.5*(orig_z(i,j,khi0)+orig_z(i,j-1,khi0)); + } + } + amrex::Real y0 = orig_data(i,j,khi0,src_comp); + amrex::Real y1 = orig_data(i,j,khi1,src_comp); + return ( y0+(y1-y0)/(z1-z0)*(z-z0) ); } else { // interpolate amrex::Real y0 = orig_data(i,j,klow,src_comp);