From 68b5393e8954d6264df6a0fcc7249a9a64a34aee Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Thu, 7 Dec 2023 14:36:32 -0800 Subject: [PATCH 1/5] WIP metgrid initialization debugging lateral boundaries. --- Exec/DevTests/MetGrid/GNUmakefile | 2 - Exec/DevTests/MetGrid/inputs | 10 +++-- Source/IO/ReadFromMetgrid.cpp | 42 +++++++++++++++---- .../Initialization/ERF_init_from_metgrid.cpp | 5 ++- Source/Initialization/Metgrid_utils.H | 4 +- Submodules/AMReX | 2 +- 6 files changed, 48 insertions(+), 17 deletions(-) diff --git a/Exec/DevTests/MetGrid/GNUmakefile b/Exec/DevTests/MetGrid/GNUmakefile index 0b412a6bb..15f8bda96 100644 --- a/Exec/DevTests/MetGrid/GNUmakefile +++ b/Exec/DevTests/MetGrid/GNUmakefile @@ -27,9 +27,7 @@ USE_ASSERTION = TRUE USE_NETCDF = TRUE -#USE_MOISTURE = FALSE USE_MOISTURE = TRUE - #USE_WARM_NO_PRECIP = TRUE # GNU Make diff --git a/Exec/DevTests/MetGrid/inputs b/Exec/DevTests/MetGrid/inputs index 00870c0c9..b6997b59c 100644 --- a/Exec/DevTests/MetGrid/inputs +++ b/Exec/DevTests/MetGrid/inputs @@ -19,7 +19,7 @@ zlo.type = "NoSlipWall" zhi.type = "SlipWall" # TIME STEP CONTROL -erf.fixed_dt = 1.0 # fixed time step depending on grid resolution +erf.fixed_dt = 2.0 # fixed time step depending on grid resolution erf.use_terrain = true erf.terrain_smoothing = 2 @@ -40,7 +40,7 @@ erf.check_int = 100 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name 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 +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 pert_dens KE QKE # SOLVER CHOICE erf.alpha_T = 1.0 @@ -52,12 +52,14 @@ erf.les_type = "Smagorinsky" #erf.les_type = "Deardorff" erf.Cs = 0.1 -erf.moisture_model = "SAM" +erf.moisture_model = "Kessler" +#erf.do_cloud = false +#erf.do_precip = false #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.metgrid_bdy_width = 5 +erf.metgrid_bdy_width = 1 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" diff --git a/Source/IO/ReadFromMetgrid.cpp b/Source/IO/ReadFromMetgrid.cpp index 43a82fb7d..4d4058bda 100644 --- a/Source/IO/ReadFromMetgrid.cpp +++ b/Source/IO/ReadFromMetgrid.cpp @@ -27,13 +27,41 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname, auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4); { // Global Attributes (int) std::vector attr; - ncf.get_attr("FLAG_PSFC", attr); flag_psfc = attr[0]; - ncf.get_attr("FLAG_MAPFAC_U", attr); flag_msfu = attr[0]; - ncf.get_attr("FLAG_MAPFAC_V", attr); flag_msfv = attr[0]; - ncf.get_attr("FLAG_MAPFAC_M", attr); flag_msfm = attr[0]; - ncf.get_attr("FLAG_HGT_M", attr); flag_hgt = attr[0]; - ncf.get_attr("FLAG_SST", attr); flag_sst = attr[0]; - ncf.get_attr("FLAG_LANDMASK", attr); flag_lmask = attr[0]; + if (ncf.has_attr("FLAG_PSFC")) { + ncf.get_attr("FLAG_PSFC", attr); flag_psfc = attr[0]; + } else { + flag_psfc = 0; + } + if (ncf.has_attr("FLAG_MAPFAC_U")) { + ncf.get_attr("FLAG_MAPFAC_U", attr); flag_msfu = attr[0]; + } else { + flag_msfu = 0; + } + if (ncf.has_attr("FLAG_MAPFAC_V")) { + ncf.get_attr("FLAG_MAPFAC_V", attr); flag_msfv = attr[0]; + } else { + flag_msfv = 0; + } + if (ncf.has_attr("FLAG_MAPFAC_M")) { + ncf.get_attr("FLAG_MAPFAC_M", attr); flag_msfm = attr[0]; + } else { + flag_msfm = 0; + } + if (ncf.has_attr("FLAG_HGT_M")) { + ncf.get_attr("FLAG_HGT_M", attr); flag_hgt = attr[0]; + } else { + flag_hgt = 0; + } + if (ncf.has_attr("FLAG_SST")) { + ncf.get_attr("FLAG_SST", attr); flag_sst = attr[0]; + } else { + flag_sst = 0; + } + if (ncf.has_attr("FLAG_LANDMASK")) { + ncf.get_attr("FLAG_LANDMASK", attr); flag_lmask = attr[0]; + } else { + flag_lmask = 0; + } ncf.get_attr("WEST-EAST_GRID_DIMENSION", attr); NC_nx = attr[0]; ncf.get_attr("SOUTH-NORTH_GRID_DIMENSION", attr); NC_ny = attr[0]; } diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index 6037c0de7..bc2d51d1c 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -351,6 +351,9 @@ ERF::init_from_metgrid (int lev) fabs_for_bcs, mask_c_arr); } // mf + // Set up boxes for lateral boundary arrays. + if (metgrid_bdy_width-1 <= metgrid_bdy_set_width) metgrid_bdy_set_width = metgrid_bdy_width; + if (metgrid_bdy_width == metgrid_bdy_set_width) metgrid_bdy_width += 1; // NOTE: fabs_for_bcs is defined over the whole domain on each rank. // However, the operations needed to define the data on the ERF @@ -500,7 +503,7 @@ ERF::init_from_metgrid (int lev) xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag; ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag; } else if (ivar == MetGridBdyVars::QV) { - multiply_rho = true; + multiply_rho = false; xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag; ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag; } // MetGridBdyVars::QV diff --git a/Source/Initialization/Metgrid_utils.H b/Source/Initialization/Metgrid_utils.H index af06ecd5d..43dd92bf9 100644 --- a/Source/Initialization/Metgrid_utils.H +++ b/Source/Initialization/Metgrid_utils.H @@ -153,14 +153,14 @@ calc_rho_p (const int& kmax, // integrate from the top back down to get dry pressure and density. Pd_vec[kmax] = Pm_vec[kmax]; - Rhod_vec[kmax] = 1.0/(R_d/p_0*Thetam_vec[kmax]*std::pow(Pd_vec[kmax]/p_0, -iGamma)); + Rhod_vec[kmax] = 1.0/(R_d/p_0*Thetad_vec[kmax]*std::pow(Pd_vec[kmax]/p_0, -iGamma)); for (int k=kmax-1; k>=0; k--) { amrex::Real dz = z_vec[k+1]-z_vec[k]; Rhod_vec[k] = Rhod_vec[k+1]; // an initial guess. for (int it=0; it < maxiter; it++) { Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*CONST_GRAV; if (Pd_vec[k] < 0.0) Pd_vec[k] = 0.0; - Rhod_vec[k] = 1.0/(R_d/p_0*Thetam_vec[k]*std::pow(Pd_vec[k]/p_0, -iGamma)); + Rhod_vec[k] = 1.0/(R_d/p_0*Thetad_vec[k]*std::pow(Pd_vec[k]/p_0, -iGamma)); } // it } // k } diff --git a/Submodules/AMReX b/Submodules/AMReX index b2052f2d2..7ee29121e 160000 --- a/Submodules/AMReX +++ b/Submodules/AMReX @@ -1 +1 @@ -Subproject commit b2052f2d2e5ce44317450ba13de705a3e01ef0ea +Subproject commit 7ee29121ed70d7e255ad98a8b1690d345cb4fb33 From a290beb5c1818203f9968582d6b2db9b58a1f0d9 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Fri, 15 Dec 2023 12:49:55 -0800 Subject: [PATCH 2/5] WIP debugging lateral boundary condition updates for metgrid initialization and forcing. --- Exec/DevTests/MetGrid/inputs | 12 +- Exec/RegTests/Bubble/inputs_squall2d_x | 14 +- .../BoundaryConditions_metgrid.cpp | 44 +++--- .../Initialization/ERF_init_from_metgrid.cpp | 148 ++++++++++-------- Source/Initialization/Metgrid_utils.H | 23 +-- 5 files changed, 120 insertions(+), 121 deletions(-) diff --git a/Exec/DevTests/MetGrid/inputs b/Exec/DevTests/MetGrid/inputs index b6997b59c..0fb942e26 100644 --- a/Exec/DevTests/MetGrid/inputs +++ b/Exec/DevTests/MetGrid/inputs @@ -1,5 +1,5 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 100 +max_step = 10 amrex.fpe_trap_invalid = 1 amrex.fpe_trap_zero = 1 @@ -11,10 +11,10 @@ amr.n_cell = 140 80 100 geometry.is_periodic = 0 0 0 -xlo.type = "Outflow" -xhi.type = "Outflow" -ylo.type = "Outflow" -yhi.type = "Outflow" +xlo.type = "outflow" +xhi.type = "outflow" +ylo.type = "outflow" +yhi.type = "outflow" zlo.type = "NoSlipWall" zhi.type = "SlipWall" @@ -40,7 +40,7 @@ erf.check_int = 100 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name 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 pert_dens KE QKE +erf.plot_vars_1 = qt qp qv qc qi qrain qsnow qgraup density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres pert_dens KE QKE # SOLVER CHOICE erf.alpha_T = 1.0 diff --git a/Exec/RegTests/Bubble/inputs_squall2d_x b/Exec/RegTests/Bubble/inputs_squall2d_x index d561668c4..9d8c67dce 100644 --- a/Exec/RegTests/Bubble/inputs_squall2d_x +++ b/Exec/RegTests/Bubble/inputs_squall2d_x @@ -34,7 +34,7 @@ erf.check_int = 9000 # 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 x_velocity y_velocity z_velocity pressure theta pres_hse +erf.plot_vars_1 = qt qp qv qc qi qrain qsnow qgraup density x_velocity y_velocity z_velocity pressure theta pres_hse # SOLVER CHOICES erf.use_gravity = true @@ -44,14 +44,16 @@ erf.rotational_time_period = 86400.0 erf.use_rayleigh_damping = false # TODO: enable # PHYSICS OPTIONS -erf.les_type = "Deardorff" +#erf.les_type = "Deardorff" +erf.les_type = "None" erf.pbl_type = "None" -erf.moisture_model = "SAM" +#erf.moisture_model = "SAM" +erf.moisture_model = "Kessler" -erf.molec_diff_type = "Constant" +erf.molec_diff_type = "ConstantAlpha" erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities -erf.dynamicViscosity= 0.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s -erf.alpha_T = 0.0 # [m^2/s] +erf.dynamicViscosity= 75.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s +erf.alpha_T = 75.0 # [m^2/s] # INITIAL CONDITIONS erf.init_type = "input_sounding" diff --git a/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp b/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp index 37f4dc04f..0bf84d0c0 100644 --- a/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp @@ -79,9 +79,10 @@ ERF::fill_from_metgrid (const Vector& mfs, int offset = (var_idx == Vars::cons) ? icomp_cons : 0; // Loop over each component - for (int comp_idx(0); comp_idx < comp_var[var_idx]; ++comp_idx) + for (int comp_idx(offset); comp_idx < (comp_var[var_idx]+offset); ++comp_idx) { - int width = metgrid_bdy_set_width; + int width = metgrid_bdy_width; + int set_width = metgrid_bdy_set_width; // Variable can be read from met_em files //------------------------------------ @@ -90,19 +91,19 @@ ERF::fill_from_metgrid (const Vector& mfs, int ivar = ind_map[var_idx][comp_idx]; IntVect ng_vect = mf.nGrowVect(); ng_vect[2] = 0; -// if (ivar == MetGridBdyVars::U) { -// amrex::Print() << "fill_from_metgrid U var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl; -// } else if (ivar == MetGridBdyVars::V) { -// amrex::Print() << "fill_from_metgrid V var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl; -// } else if (ivar == MetGridBdyVars::R) { -// amrex::Print() << "fill_from_metgrid R var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl; -// } else if (ivar == MetGridBdyVars::T) { -// amrex::Print() << "fill_from_metgrid T var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl; -// } else if (ivar == MetGridBdyVars::QV) { -// amrex::Print() << "fill_from_metgrid QV var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl; -// } else { -// amrex::Print() << "fill_from_metgrid UNKNOWN" << std::endl; -// } + if (ivar == MetGridBdyVars::U) { + amrex::Print() << "fill_from_metgrid U var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << " ng_vect=" << ng_vect << std::endl; + } else if (ivar == MetGridBdyVars::V) { + amrex::Print() << "fill_from_metgrid V var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << " ng_vect=" << ng_vect << std::endl; + } else if (ivar == MetGridBdyVars::R) { + amrex::Print() << "fill_from_metgrid R var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << " ng_vect=" << ng_vect << std::endl; + } else if (ivar == MetGridBdyVars::T) { + amrex::Print() << "fill_from_metgrid T var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << " ng_vect=" << ng_vect << std::endl; + } else if (ivar == MetGridBdyVars::QV) { + amrex::Print() << "fill_from_metgrid QV var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << " ng_vect=" << ng_vect << std::endl; + } else { + amrex::Print() << "fill_from_metgrid UNKNOWN" << std::endl; + } // We have data at fixed time intervals we will call dT // Then to interpolate, given time, we can define n = (time/dT) @@ -125,10 +126,8 @@ ERF::fill_from_metgrid (const Vector& mfs, // Grown tilebox so we fill exterior ghost cells as well Box gbx = mfi.growntilebox(ng_vect); const Array4& dest_arr = mf.array(mfi); - - // Call w/o interior ghost cells Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; - compute_interior_ghost_bxs_xy(gbx, domain, width, 0, + compute_interior_ghost_bxs_xy(gbx, domain, width, set_width, bx_xlo, bx_xhi, bx_ylo, bx_yhi, ng_vect); @@ -170,7 +169,6 @@ ERF::fill_from_metgrid (const Vector& mfs, // Variable not read or computed from met_em files //------------------------------------ } else { - width = metgrid_bdy_width; IntVect ng_vect = mf.nGrowVect(); ng_vect[2] = 0; #ifdef AMREX_USE_OMP @@ -190,14 +188,14 @@ ERF::fill_from_metgrid (const Vector& mfs, ParallelFor(bx_xlo, bx_xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - int jj = std::max(j , dom_lo.y); - jj = std::min(jj, dom_hi.y); + int jj = std::max(j , dom_lo.y+width); + jj = std::min(jj, dom_hi.y-width); dest_arr(i,j,k,comp_idx) = dest_arr(dom_lo.x+width,jj,k,comp_idx); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - int jj = std::max(j , dom_lo.y); - jj = std::min(jj, dom_hi.y); + int jj = std::max(j , dom_lo.y+width); + jj = std::min(jj, dom_hi.y-width); dest_arr(i,j,k,comp_idx) = dest_arr(dom_hi.x-width,jj,k,comp_idx); }); diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index bc2d51d1c..f2170fd4f 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -235,7 +235,6 @@ ERF::init_from_metgrid (int lev) #else int MetGridBdyEnd = MetGridBdyVars::NumTypes-1; #endif - //amrex::Vector > fabs_for_bcs; amrex::Vector> fabs_for_bcs; fabs_for_bcs.resize(ntimes); for (int it(0); it < ntimes; it++) { @@ -476,7 +475,8 @@ ERF::init_from_metgrid (int lev) amrex::Box xlo_plane, xhi_plane, ylo_plane, yhi_plane; for (int it(0); it < ntimes; it++) { - const Array4& R_bcs_arr = fabs_for_bcs[it][MetGridBdyVars::R].const_array(); +// const Array4& R_bcs_arr = fabs_for_bcs[it][MetGridBdyVars::R].const_array(); +// const Array4& R_bcs_arr = fabs_for_bcs[0][MetGridBdyVars::R].const_array(); for (int ivar(MetGridBdyVars::U); ivar < MetGridBdyEnd; ivar++) { @@ -484,7 +484,11 @@ ERF::init_from_metgrid (int lev) auto xhi_arr = bdy_data_xhi[it][ivar].array(); auto ylo_arr = bdy_data_ylo[it][ivar].array(); auto yhi_arr = bdy_data_yhi[it][ivar].array(); - const Array4& fabs_for_bcs_arr = fabs_for_bcs[it][ivar].const_array(); +// const Array4& fabs_for_bcs_arr = fabs_for_bcs[it][ivar].const_array(); +// --------------------------------------------------------------------------------------- +// DJW: remove this after debugging is complete. This enforces steady boundary conditions. + const Array4& fabs_for_bcs_arr = fabs_for_bcs[0][ivar].const_array(); +// --------------------------------------------------------------------------------------- if (ivar == MetGridBdyVars::U) { multiply_rho = false; @@ -511,26 +515,30 @@ ERF::init_from_metgrid (int lev) // west boundary amrex::ParallelFor(xlo_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; - xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; +// amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; +// xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; + xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); }); // xvel at east boundary amrex::ParallelFor(xhi_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; - xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; +// amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; +// xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; + xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); }); // xvel at south boundary amrex::ParallelFor(ylo_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; - ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; +// amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; +// ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; + ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); }); // xvel at north boundary amrex::ParallelFor(yhi_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; - yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; +// amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; +// yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; + yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); }); } // ivar @@ -686,6 +694,15 @@ init_state_from_metgrid (const Real l_rdOcp, }); } +// --------------------------------------------------------------------------------------- +// DJW: remove this after debugging is complete. This enforces quiescent conditions.. + if (it == 0) { + x_vel_fab.template setVal(0.0); + y_vel_fab.template setVal(0.0); + } + fabs_for_bcs[it][MetGridBdyVars::U].template setVal(0.0); + fabs_for_bcs[it][MetGridBdyVars::V].template setVal(0.0); +// --------------------------------------------------------------------------------------- // ******************************************************** // W @@ -701,45 +718,6 @@ init_state_from_metgrid (const Real l_rdOcp, state_fab.template setVal(0.0); } - - // ******************************************************** - // theta - // ******************************************************** - { // calculate potential temperature. - Box bx = NC_rhum_fab[it].box() & tbxc; - auto const temp = NC_temp_fab[it].const_array(); - auto const pres = NC_pres_fab[it].const_array(); - auto theta = theta_fab[it].array(); - - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - theta(i,j,k) = getThgivenPandT(temp(i,j,k),pres(i,j,k),l_rdOcp); - //theta(i,j,k) = 300.0; // TODO: Remove when not needed. Force an isothermal atmosphere for debugging. - }); - } - - // vertical interpolation of potential temperature. - { - Box bx2d = NC_temp_fab[it].box() & tbxc; - bx2d.setRange(2,0); - auto const orig_data = theta_fab[it].const_array(); - auto const orig_z = NC_ght_fab[it].const_array(); - auto new_data = state_fab.array(); - auto bc_data = fabs_for_bcs[it][MetGridBdyVars::T].array(); - auto const new_z = z_phys_nd_fab.const_array(); - - int kmax = amrex::ubound(tbxc).z; - - ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept - { - for (int k = 0; k<=kmax; k++) { - Real Interp_Val = interpolate_column_metgrid(i,j,k,'M',0,orig_z,orig_data,new_z); - if (mask_c_arr(i,j,k)) bc_data(i,j,k,0) = Interp_Val; - if (it==0) new_data(i,j,k,RhoTheta_comp) = Interp_Val; - } - }); - } - #if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) // ******************************************************** // specific humidity / relative humidity / mixing ratio @@ -790,13 +768,46 @@ init_state_from_metgrid (const Real l_rdOcp, } #endif - // TODO: TEMPORARY CODE TO RUN QUIESCENT, REMOVE WHEN NOT NEEDED. -// if (it == 0) { -// x_vel_fab.template setVal(0.0); // TODO: temporary code to initialize with quiescent atmosphere. -// y_vel_fab.template setVal(0.0); // TODO: temporary code to initialize with quiescent atmosphere. -// } -// fabs_for_bcs[it][MetGridBdyVars::U].template setVal(0.0); // TODO: temporary code to force with quiescent atmosphere. -// fabs_for_bcs[it][MetGridBdyVars::V].template setVal(0.0); // TODO: temporary code to force with quiescent atmosphere. + // ******************************************************** + // theta + // ******************************************************** + { // calculate potential temperature. + Box bx = NC_rhum_fab[it].box() & tbxc; + auto const temp = NC_temp_fab[it].const_array(); + auto const pres = NC_pres_fab[it].const_array(); + auto theta = theta_fab[it].array(); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + theta(i,j,k) = getThgivenPandT(temp(i,j,k),pres(i,j,k),l_rdOcp); +// --------------------------------------------------------------------------------------- +// DJW: remove this after debugging is complete. This enforces isothermal conditions.. + //theta(i,j,k) = 300.0; +// --------------------------------------------------------------------------------------- + }); + } + + // vertical interpolation of potential temperature. + { + Box bx2d = NC_temp_fab[it].box() & tbxc; + bx2d.setRange(2,0); + auto const orig_data = theta_fab[it].const_array(); + auto const orig_z = NC_ght_fab[it].const_array(); + auto new_data = state_fab.array(); + auto bc_data = fabs_for_bcs[it][MetGridBdyVars::T].array(); + auto const new_z = z_phys_nd_fab.const_array(); + + int kmax = amrex::ubound(tbxc).z; + + ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept + { + for (int k = 0; k<=kmax; k++) { + Real Interp_Val = interpolate_column_metgrid(i,j,k,'M',0,orig_z,orig_data,new_z); + if (mask_c_arr(i,j,k)) bc_data(i,j,k,0) = Interp_Val; + if (it==0) new_data(i,j,k,RhoTheta_comp) = Interp_Val; + } + }); + } } // it } @@ -840,7 +851,6 @@ init_base_state_from_metgrid (const Real l_rdOcp, // Device vectors for columnwise operations Gpu::DeviceVector z_vec_d(kmax+2,0); Real* z_vec = z_vec_d.data(); - Gpu::DeviceVector Thetad_vec_d(kmax+1,0); Real* Thetad_vec = Thetad_vec_d.data(); Gpu::DeviceVector Thetam_vec_d(kmax+1,0); Real* Thetam_vec = Thetam_vec_d.data(); Gpu::DeviceVector Rhod_vec_d(kmax+1,0); Real* Rhod_vec = Rhod_vec_d.data(); Gpu::DeviceVector Rhom_vec_d(kmax+1,0); Real* Rhom_vec = Rhom_vec_d.data(); @@ -874,14 +884,17 @@ init_base_state_from_metgrid (const Real l_rdOcp, { for (int k=0; k<=kmax; k++) { z_vec[k] = new_z(i,j,k); - Thetad_vec[k] = new_data(i,j,k,RhoTheta_comp); #if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) Q_vec[k] = new_data(i,j,k,RhoQ_comp); + amrex::Real qvf = 1.0+(R_v/R_d+1.0)*Q_vec[0]; +#else + amrex::Real qvf = 1.0; #endif + Thetam_vec[k] = new_data(i,j,k,RhoTheta_comp)*qvf; } z_vec[kmax+1] = new_z(i,j,kmax+1); - calc_rho_p(kmax,flag_psfc_vec[0],orig_psfc(i,j,0),Thetad_vec,Thetam_vec, + calc_rho_p(kmax,flag_psfc_vec[0],orig_psfc(i,j,0),Thetam_vec, #if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) Q_vec, #endif @@ -937,15 +950,18 @@ init_base_state_from_metgrid (const Real l_rdOcp, amrex::ParallelFor(valid_bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { for (int k=0; k<=kmax; k++) { - z_vec[k] = new_z(i,j,k); - Thetad_vec[k] = Theta_arr(i,j,k); + z_vec[k] = new_z(i,j,k); #if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) - Q_vec[k] = Q_arr(i,j,k); + Q_vec[k] = Q_arr(i,j,k); + amrex::Real qvf = 1.0+(R_v/R_d+1.0)*Q_vec[0]; +#else + amrex::Real qvf = 1.0; #endif + Thetam_vec[k] = Theta_arr(i,j,k)*qvf; } - z_vec[kmax+1] = new_z(i,j,kmax+1); + z_vec[kmax+1] = new_z(i,j,kmax+1); - calc_rho_p(kmax,flag_psfc_vec[it],orig_psfc(i,j,0),Thetad_vec,Thetam_vec, + calc_rho_p(kmax,flag_psfc_vec[it],orig_psfc(i,j,0),Thetam_vec, #if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) Q_vec, #endif @@ -958,7 +974,7 @@ init_base_state_from_metgrid (const Real l_rdOcp, #if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) Q_arr(i,j,k) = Rhod_vec[k]*Q_vec[k]; #endif - Theta_arr(i,j,k) = Rhod_vec[k]*Thetad_vec[k]; + Theta_arr(i,j,k) = Rhod_vec[k]*Theta_arr(i,j,k); } } // k }); diff --git a/Source/Initialization/Metgrid_utils.H b/Source/Initialization/Metgrid_utils.H index 43dd92bf9..bb7f830c1 100644 --- a/Source/Initialization/Metgrid_utils.H +++ b/Source/Initialization/Metgrid_utils.H @@ -99,7 +99,6 @@ void calc_rho_p (const int& kmax, const int& flag_psfc, const amrex::Real& psfc, - amrex::Real* Thetad_vec, amrex::Real* Thetam_vec, #if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) amrex::Real* Q_vec, @@ -121,29 +120,13 @@ calc_rho_p (const int& kmax, Pm_vec[0] = p_0*exp(-t_0/a+std::pow((std::pow(t_0/a, 2)-2.0*CONST_GRAV*z_vec[0]/(a*R_d)), 0.5)); } - // calculate virtual potential temperature at the surface. - { -#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) - amrex::Real qvf = 1.0+(R_v/R_d+1.0)*Q_vec[0]; -#else - amrex::Real qvf = 1.0; -#endif - Thetam_vec[0] = Thetad_vec[0]*qvf; - } - // calculate moist density at the surface. Rhom_vec[0] = 1.0/(R_d/p_0*Thetam_vec[0]*std::pow(Pm_vec[0]/p_0, -iGamma)); // integrate from the surface to the top boundary. for (int k=1; k<=kmax; k++) { amrex::Real dz = z_vec[k]-z_vec[k-1]; -#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) - amrex::Real qvf = 1.0+(R_v/R_d+1.0)*Q_vec[k]; -#else - amrex::Real qvf = 1.0; -#endif - Thetam_vec[k] = Thetad_vec[k]*qvf; - Rhom_vec[k] = Rhom_vec[k-1]; // an initial guess. + Rhom_vec[k] = Rhom_vec[k-1]; // an initial guess. for (int it=0; it=0; k--) { amrex::Real dz = z_vec[k+1]-z_vec[k]; Rhod_vec[k] = Rhod_vec[k+1]; // an initial guess. for (int it=0; it < maxiter; it++) { Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*CONST_GRAV; if (Pd_vec[k] < 0.0) Pd_vec[k] = 0.0; - Rhod_vec[k] = 1.0/(R_d/p_0*Thetad_vec[k]*std::pow(Pd_vec[k]/p_0, -iGamma)); + Rhod_vec[k] = 1.0/(R_d/p_0*Thetam_vec[k]*std::pow(Pd_vec[k]/p_0, -iGamma)); } // it } // k } From cbf463b8d49c6000331e722edf1a024f618acc52 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Thu, 18 Jan 2024 11:37:15 -0800 Subject: [PATCH 3/5] WIP debugging metgrid initialization and forcing --- Exec/DevTests/MetGrid/inputs | 6 +- .../BoundaryConditions_metgrid.cpp | 45 ++++++- .../Initialization/ERF_init_from_metgrid.cpp | 115 +++++++++++------- Source/Initialization/Metgrid_utils.H | 29 +++-- 4 files changed, 140 insertions(+), 55 deletions(-) diff --git a/Exec/DevTests/MetGrid/inputs b/Exec/DevTests/MetGrid/inputs index 0fb942e26..cc128fdb0 100644 --- a/Exec/DevTests/MetGrid/inputs +++ b/Exec/DevTests/MetGrid/inputs @@ -1,5 +1,5 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 10 +max_step = 1 amrex.fpe_trap_invalid = 1 amrex.fpe_trap_zero = 1 @@ -59,8 +59,8 @@ erf.moisture_model = "Kessler" #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.metgrid_bdy_width = 1 -erf.metgrid_bdy_set_width = 1 +erf.metgrid_bdy_width = 3 +erf.metgrid_bdy_set_width = 3 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" diff --git a/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp b/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp index 0bf84d0c0..07ba21bee 100644 --- a/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp @@ -104,6 +104,7 @@ ERF::fill_from_metgrid (const Vector& mfs, } else { amrex::Print() << "fill_from_metgrid UNKNOWN" << std::endl; } + amrex::Print() << " n_time=" << n_time << " \toma=" << oma << " \talpha=" << alpha << std::endl; // We have data at fixed time intervals we will call dT // Then to interpolate, given time, we can define n = (time/dT) @@ -127,10 +128,39 @@ ERF::fill_from_metgrid (const Vector& mfs, Box gbx = mfi.growntilebox(ng_vect); const Array4& dest_arr = mf.array(mfi); Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; - compute_interior_ghost_bxs_xy(gbx, domain, width, set_width, + compute_interior_ghost_bxs_xy(gbx, domain, width, 0, bx_xlo, bx_xhi, bx_ylo, bx_yhi, ng_vect); + amrex::AllPrint() << "fill_from_metgrid bx_xlo=" << bx_xlo << std::endl; + amrex::AllPrint() << "fill_from_metgrid bx_xhi=" << bx_xhi << std::endl; + amrex::AllPrint() << "fill_from_metgrid bx_ylo=" << bx_ylo << std::endl; + amrex::AllPrint() << "fill_from_metgrid bx_yhi=" << bx_yhi << std::endl; + if (ivar == MetGridBdyVars::QV) { + amrex::Print() << " QV BEFORE" << std::endl; + int jj = (dom_hi.y - dom_lo.y)/2+dom_lo.y; + for (int kk = 9; kk >= 0; kk--) { + for (int ii = 0; ii < 10; ii++) { + amrex::Print() << " " << dest_arr(ii,jj,kk,comp_idx); + } + amrex::Print() << std::endl; + } + amrex::Print() << " bdatxlo_n" << std::endl; + for (int k = 9; k >= 0; k--) { + for (int ii = std::max(dom_lo.x, amrex::lbound(bx_xlo).x); ii <= std::min(dom_hi.x, amrex::ubound(bx_xlo).x); ii++) { + amrex::Print() << " " << bdatxlo_n(ii,jj,k,0); + } + amrex::Print() << std::endl; + } + amrex::Print() << " bdatxlo_np1" << std::endl; + for (int k = 9; k >= 0; k--) { + for (int ii = std::max(dom_lo.x, amrex::lbound(bx_xlo).x); ii <= std::min(dom_hi.x, amrex::ubound(bx_xlo).x); ii++) { + amrex::Print() << " " << bdatxlo_np1(ii,jj,k,0); + } + amrex::Print() << std::endl; + } + } + // x-faces (includes exterior y ghost cells) ParallelFor(bx_xlo, bx_xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) @@ -164,6 +194,17 @@ ERF::fill_from_metgrid (const Vector& mfs, dest_arr(i,j,k,comp_idx) = oma * bdatyhi_n (i,jj,k,0) + alpha * bdatyhi_np1(i,jj,k,0); }); + + if (ivar == MetGridBdyVars::QV) { + amrex::Print() << " QV AFTER" << std::endl; + int jj = (dom_hi.y - dom_lo.y)/2+dom_lo.y; + for (int kk = 9; kk >= 0; kk--) { + for (int ii = 0; ii < 10; ii++) { + amrex::Print() << " " << dest_arr(ii,jj,kk,comp_idx); + } + amrex::Print() << std::endl; + } + } } // mfi // Variable not read or computed from met_em files @@ -182,7 +223,7 @@ ERF::fill_from_metgrid (const Vector& mfs, Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; compute_interior_ghost_bxs_xy(gbx, domain, width, 0, bx_xlo, bx_xhi, - bx_ylo, bx_yhi, ng_vect); + bx_ylo, bx_yhi, ng_vect, true); // x-faces (includes y ghost cells) ParallelFor(bx_xlo, bx_xhi, diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index f2170fd4f..caac0e6dc 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -368,6 +368,16 @@ ERF::init_from_metgrid (int lev) } } + for (int it(0); it < ntimes; it++) { + const Array4& fabs_for_bcs_arr = fabs_for_bcs[it][MetGridBdyVars::QV].const_array(); + amrex::Print() << " fabs_for_bcs[" << it << "][" << MetGridBdyVars::QV << "] after ParallelAllReduce::Sum" << std::endl; + for (int k=9; k >= 0; k--) { + for (int ii=0; ii <= 9; ii++) { + amrex::Print() << " " << fabs_for_bcs_arr(ii, 40, k); + } + amrex::Print() << std::endl; + } + } // Set up boxes for lateral boundary arrays. bdy_data_xlo.resize(ntimes); @@ -434,26 +444,31 @@ ERF::init_from_metgrid (int lev) for (int ivar(MetGridBdyVars::U); ivar < MetGridBdyEnd; ivar++) { for (int it(0); it < ntimes; it++) { if (ivar == MetGridBdyVars::U) { + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tbdy_data_*[" << it << "] for U" << std::endl; bdy_data_xlo[it].push_back(FArrayBox(xlo_plane_x_stag, 1)); bdy_data_xhi[it].push_back(FArrayBox(xhi_plane_x_stag, 1)); bdy_data_ylo[it].push_back(FArrayBox(ylo_plane_x_stag, 1)); bdy_data_yhi[it].push_back(FArrayBox(yhi_plane_x_stag, 1)); } else if (ivar == MetGridBdyVars::V) { + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tbdy_data_*[" << it << "] for V" << std::endl; bdy_data_xlo[it].push_back(FArrayBox(xlo_plane_y_stag, 1)); bdy_data_xhi[it].push_back(FArrayBox(xhi_plane_y_stag, 1)); bdy_data_ylo[it].push_back(FArrayBox(ylo_plane_y_stag, 1)); bdy_data_yhi[it].push_back(FArrayBox(yhi_plane_y_stag, 1)); } else if (ivar == MetGridBdyVars::R) { + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tbdy_data_*[" << it << "] for R" << std::endl; bdy_data_xlo[it].push_back(FArrayBox(xlo_plane_no_stag, 1)); bdy_data_xhi[it].push_back(FArrayBox(xhi_plane_no_stag, 1)); bdy_data_ylo[it].push_back(FArrayBox(ylo_plane_no_stag, 1)); bdy_data_yhi[it].push_back(FArrayBox(yhi_plane_no_stag, 1)); } else if (ivar == MetGridBdyVars::T) { + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tbdy_data_*[" << it << "] for T" << std::endl; bdy_data_xlo[it].push_back(FArrayBox(xlo_plane_no_stag, 1)); bdy_data_xhi[it].push_back(FArrayBox(xhi_plane_no_stag, 1)); bdy_data_ylo[it].push_back(FArrayBox(ylo_plane_no_stag, 1)); bdy_data_yhi[it].push_back(FArrayBox(yhi_plane_no_stag, 1)); } else if (ivar == MetGridBdyVars::QV) { + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tbdy_data_*[" << it << "] for QV" << std::endl; bdy_data_xlo[it].push_back(FArrayBox(xlo_plane_no_stag, 1)); bdy_data_xhi[it].push_back(FArrayBox(xhi_plane_no_stag, 1)); bdy_data_ylo[it].push_back(FArrayBox(ylo_plane_no_stag, 1)); @@ -475,69 +490,78 @@ ERF::init_from_metgrid (int lev) amrex::Box xlo_plane, xhi_plane, ylo_plane, yhi_plane; for (int it(0); it < ntimes; it++) { -// const Array4& R_bcs_arr = fabs_for_bcs[it][MetGridBdyVars::R].const_array(); -// const Array4& R_bcs_arr = fabs_for_bcs[0][MetGridBdyVars::R].const_array(); - for (int ivar(MetGridBdyVars::U); ivar < MetGridBdyEnd; ivar++) { auto xlo_arr = bdy_data_xlo[it][ivar].array(); auto xhi_arr = bdy_data_xhi[it][ivar].array(); auto ylo_arr = bdy_data_ylo[it][ivar].array(); auto yhi_arr = bdy_data_yhi[it][ivar].array(); -// const Array4& fabs_for_bcs_arr = fabs_for_bcs[it][ivar].const_array(); + const Array4& fabs_for_bcs_arr = fabs_for_bcs[it][ivar].const_array(); // --------------------------------------------------------------------------------------- // DJW: remove this after debugging is complete. This enforces steady boundary conditions. - const Array4& fabs_for_bcs_arr = fabs_for_bcs[0][ivar].const_array(); +// const Array4& fabs_for_bcs_arr = fabs_for_bcs[0][ivar].const_array(); // --------------------------------------------------------------------------------------- if (ivar == MetGridBdyVars::U) { multiply_rho = false; xlo_plane = xlo_plane_x_stag; xhi_plane = xhi_plane_x_stag; ylo_plane = ylo_plane_x_stag; yhi_plane = yhi_plane_x_stag; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane " << xlo_plane << " \tfor U" << std::endl; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane " << xhi_plane << " \tfor U" << std::endl; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane " << ylo_plane << " \tfor U" << std::endl; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane " << yhi_plane << " \tfor U" << std::endl; } else if (ivar == MetGridBdyVars::V) { multiply_rho = false; xlo_plane = xlo_plane_y_stag; xhi_plane = xhi_plane_y_stag; ylo_plane = ylo_plane_y_stag; yhi_plane = yhi_plane_y_stag; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane " << xlo_plane << " \tfor V" << std::endl; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane " << xhi_plane << " \tfor V" << std::endl; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane " << ylo_plane << " \tfor V" << std::endl; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane " << yhi_plane << " \tfor V" << std::endl; } else if (ivar == MetGridBdyVars::R) { multiply_rho = false; xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag; ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane " << xlo_plane << " \tfor R" << std::endl; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane " << xhi_plane << " \tfor R" << std::endl; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane " << ylo_plane << " \tfor R" << std::endl; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane " << yhi_plane << " \tfor R" << std::endl; } else if (ivar == MetGridBdyVars::T) { multiply_rho = false; xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag; ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane " << xlo_plane << " \tfor T" << std::endl; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane " << xhi_plane << " \tfor T" << std::endl; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane " << ylo_plane << " \tfor T" << std::endl; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane " << yhi_plane << " \tfor T" << std::endl; } else if (ivar == MetGridBdyVars::QV) { multiply_rho = false; xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag; ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane " << xlo_plane << " \tfor QV" << std::endl; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane " << xhi_plane << " \tfor QV" << std::endl; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane " << ylo_plane << " \tfor QV" << std::endl; + amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane " << yhi_plane << " \tfor QV" << std::endl; } // MetGridBdyVars::QV // west boundary amrex::ParallelFor(xlo_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { -// amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; -// xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); }); // xvel at east boundary amrex::ParallelFor(xhi_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { -// amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; -// xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); }); // xvel at south boundary amrex::ParallelFor(ylo_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { -// amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; -// ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); }); // xvel at north boundary amrex::ParallelFor(yhi_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { -// amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; -// yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); }); @@ -667,6 +691,13 @@ init_state_from_metgrid (const Real l_rdOcp, if (it==0) new_data(i,j,k,0) = Interp_Val; } }); + amrex::Print() << " mask_u_arr" << std::endl; + for (int k=9; k >= 0; k--) { + for (int ii=0; ii >= 3; ii++) { + amrex::Print() << " " << mask_u_arr(ii,40,k); + } + amrex::Print() << std::endl; + } } @@ -757,6 +788,7 @@ init_state_from_metgrid (const Real l_rdOcp, #elif defined(ERF_USE_WARM_NO_PRECIP) int state_indx = RhoQv_comp; #endif + amrex::Print() << " bx2d for QV in init_state_from_metgrid " << bx2d << std::endl; ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { for (int k = 0; k<=kmax; k++) { @@ -765,7 +797,16 @@ init_state_from_metgrid (const Real l_rdOcp, if (it==0) new_data(i,j,k,state_indx) = Interp_Val; } }); + amrex::Print() << " mask_c_arr for QV right" << std::endl; + amrex::Print() << mask_c_arr << std::endl; + // for (int k=9; k >= 0; k--) { + // for (int ii=0; ii >= 3; ii++) { + // amrex::Print() << " " << mask_c_arr(ii,40,k); + // } + // amrex::Print() << std::endl; + // } } + #endif // ******************************************************** @@ -850,14 +891,14 @@ init_base_state_from_metgrid (const Real l_rdOcp, int kmax = amrex::ubound(valid_bx).z; // Device vectors for columnwise operations - Gpu::DeviceVector z_vec_d(kmax+2,0); Real* z_vec = z_vec_d.data(); - Gpu::DeviceVector Thetam_vec_d(kmax+1,0); Real* Thetam_vec = Thetam_vec_d.data(); - Gpu::DeviceVector Rhod_vec_d(kmax+1,0); Real* Rhod_vec = Rhod_vec_d.data(); - Gpu::DeviceVector Rhom_vec_d(kmax+1,0); Real* Rhom_vec = Rhom_vec_d.data(); - Gpu::DeviceVector Pd_vec_d(kmax+1,0); Real* Pd_vec = Pd_vec_d.data(); - Gpu::DeviceVector Pm_vec_d(kmax+1,0); Real* Pm_vec = Pm_vec_d.data(); + Gpu::DeviceVector z_vec_d(kmax+2,0); Real* z_vec = z_vec_d.data(); + Gpu::DeviceVector Theta_vec_d(kmax+1,0); Real* Theta_vec = Theta_vec_d.data(); + Gpu::DeviceVector Rhod_vec_d(kmax+1,0); Real* Rhod_vec = Rhod_vec_d.data(); + Gpu::DeviceVector Rhom_vec_d(kmax+1,0); Real* Rhom_vec = Rhom_vec_d.data(); + Gpu::DeviceVector Pd_vec_d(kmax+1,0); Real* Pd_vec = Pd_vec_d.data(); + Gpu::DeviceVector Pm_vec_d(kmax+1,0); Real* Pm_vec = Pm_vec_d.data(); #if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) - Gpu::DeviceVector Q_vec_d(kmax+1,0); Real* Q_vec = Q_vec_d.data(); + Gpu::DeviceVector Q_vec_d(kmax+1,0); Real* Q_vec = Q_vec_d.data(); #endif // Device vectors for psfc flags @@ -883,26 +924,20 @@ init_base_state_from_metgrid (const Real l_rdOcp, amrex::ParallelFor(valid_bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { for (int k=0; k<=kmax; k++) { - z_vec[k] = new_z(i,j,k); -#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) - Q_vec[k] = new_data(i,j,k,RhoQ_comp); - amrex::Real qvf = 1.0+(R_v/R_d+1.0)*Q_vec[0]; -#else - amrex::Real qvf = 1.0; -#endif - Thetam_vec[k] = new_data(i,j,k,RhoTheta_comp)*qvf; + z_vec[k] = new_z(i,j,k); + Theta_vec[k] = new_data(i,j,k,RhoTheta_comp); } z_vec[kmax+1] = new_z(i,j,kmax+1); - calc_rho_p(kmax,flag_psfc_vec[0],orig_psfc(i,j,0),Thetam_vec, + calc_rho_p(l_rdOcp, kmax,flag_psfc_vec[0],orig_psfc(i,j,0),Theta_vec, #if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) Q_vec, #endif z_vec,Rhod_vec,Rhom_vec,Pd_vec,Pm_vec); for (int k=0; k<=kmax; k++) { - p_hse_arr(i,j,k) = Pd_vec[k]; - r_hse_arr(i,j,k) = Rhod_vec[k]; + p_hse_arr(i,j,k) = Pd_vec[k]; + r_hse_arr(i,j,k) = Rhom_vec[k]; } }); @@ -951,30 +986,24 @@ init_base_state_from_metgrid (const Real l_rdOcp, { for (int k=0; k<=kmax; k++) { z_vec[k] = new_z(i,j,k); -#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) - Q_vec[k] = Q_arr(i,j,k); - amrex::Real qvf = 1.0+(R_v/R_d+1.0)*Q_vec[0]; -#else - amrex::Real qvf = 1.0; -#endif - Thetam_vec[k] = Theta_arr(i,j,k)*qvf; + Theta_vec[k] = Theta_arr(i,j,k); } z_vec[kmax+1] = new_z(i,j,kmax+1); - calc_rho_p(kmax,flag_psfc_vec[it],orig_psfc(i,j,0),Thetam_vec, + calc_rho_p(l_rdOcp, kmax,flag_psfc_vec[it],orig_psfc(i,j,0),Theta_vec, #if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) Q_vec, #endif z_vec,Rhod_vec,Rhom_vec,Pd_vec,Pm_vec); for (int k=0; k<=kmax; k++) { - p_hse_arr(i,j,k) = Pd_vec[k]; + p_hse_arr(i,j,k) = Pd_vec[k]; if (mask_c_arr(i,j,k)) { - r_hse_arr(i,j,k) = Rhod_vec[k]; #if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) - Q_arr(i,j,k) = Rhod_vec[k]*Q_vec[k]; + Q_arr(i,j,k) = Rhom_vec[k]*Q_vec[k]; #endif - Theta_arr(i,j,k) = Rhod_vec[k]*Theta_arr(i,j,k); + r_hse_arr(i,j,k) = Rhom_vec[k]; + Theta_arr(i,j,k) = Rhom_vec[k]*Theta_arr(i,j,k); } } // k }); diff --git a/Source/Initialization/Metgrid_utils.H b/Source/Initialization/Metgrid_utils.H index bb7f830c1..3929522f8 100644 --- a/Source/Initialization/Metgrid_utils.H +++ b/Source/Initialization/Metgrid_utils.H @@ -96,10 +96,11 @@ init_base_state_from_metgrid (const amrex::Real l_rdOcp, AMREX_FORCE_INLINE AMREX_GPU_DEVICE void -calc_rho_p (const int& kmax, +calc_rho_p (const amrex::Real& rdOcp, + const int& kmax, const int& flag_psfc, const amrex::Real& psfc, - amrex::Real* Thetam_vec, + amrex::Real* Theta_vec, #if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) amrex::Real* Q_vec, #endif @@ -121,7 +122,12 @@ calc_rho_p (const int& kmax, } // calculate moist density at the surface. - Rhom_vec[0] = 1.0/(R_d/p_0*Thetam_vec[0]*std::pow(Pm_vec[0]/p_0, -iGamma)); +#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) + amrex::Real Q = Q_vec[0]; +#else + amrex::Real Q = 0.0; +#endif + Rhom_vec[0] = getRhogivenThetaPress(Theta_vec[0], Pm_vec[0], rdOcp, Q)*(1.0+Q); // integrate from the surface to the top boundary. for (int k=1; k<=kmax; k++) { @@ -130,20 +136,29 @@ calc_rho_p (const int& kmax, for (int it=0; it=0; k--) { amrex::Real dz = z_vec[k+1]-z_vec[k]; Rhod_vec[k] = Rhod_vec[k+1]; // an initial guess. for (int it=0; it < maxiter; it++) { Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*CONST_GRAV; if (Pd_vec[k] < 0.0) Pd_vec[k] = 0.0; - Rhod_vec[k] = 1.0/(R_d/p_0*Thetam_vec[k]*std::pow(Pd_vec[k]/p_0, -iGamma)); +#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) + Q = Q_vec[k]; +#endif + Rhod_vec[k] = getRhogivenThetaPress(Theta_vec[k], Pd_vec[k], rdOcp, Q); } // it } // k } From 1533a2ddf6dcef4f35b7a246cb5f193774f6b47f Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Fri, 19 Jan 2024 13:11:54 -0800 Subject: [PATCH 4/5] WIP debugging metgrid initialization and forcing --- Exec/DevTests/MetGrid/inputs | 4 +- .../BoundaryConditions_metgrid.cpp | 60 +---- .../Initialization/ERF_init_from_metgrid.cpp | 226 +++++++----------- Source/Initialization/Metgrid_utils.H | 25 +- Source/Utils/InteriorGhostCells.cpp | 2 - 5 files changed, 101 insertions(+), 216 deletions(-) diff --git a/Exec/DevTests/MetGrid/inputs b/Exec/DevTests/MetGrid/inputs index cc128fdb0..5de9460c5 100644 --- a/Exec/DevTests/MetGrid/inputs +++ b/Exec/DevTests/MetGrid/inputs @@ -40,7 +40,7 @@ erf.check_int = 100 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name erf.plot_int_1 = 1 # number of timesteps between plotfiles -erf.plot_vars_1 = qt qp qv qc qi qrain qsnow qgraup density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres pert_dens KE QKE +erf.plot_vars_1 = qt qv qc density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres pert_dens KE QKE # SOLVER CHOICE erf.alpha_T = 1.0 @@ -60,7 +60,7 @@ erf.moisture_model = "Kessler" # INITIALIZATION WITH ATM DATA erf.metgrid_bdy_width = 3 -erf.metgrid_bdy_set_width = 3 +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" diff --git a/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp b/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp index a8138d393..d1dc4c369 100644 --- a/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp @@ -15,7 +15,7 @@ void ERF::fill_from_metgrid (const Vector& mfs, const Real time, bool cons_only, - int /*icomp_cons*/, + int icomp_cons, int ncomp_cons) { int lev = 0; @@ -64,6 +64,9 @@ ERF::fill_from_metgrid (const Vector& mfs, const auto& dom_lo = amrex::lbound(domain); const auto& dom_hi = amrex::ubound(domain); + // Offset only applys to cons (we may fill a subset of these vars) + int offset = (var_idx == Vars::cons) ? icomp_cons : 0; + // Loop over each component for (int comp_idx(offset); comp_idx < (comp_var[var_idx]+offset); ++comp_idx) { @@ -77,21 +80,6 @@ ERF::fill_from_metgrid (const Vector& mfs, int ivar = ind_map[var_idx][comp_idx]; IntVect ng_vect = mf.nGrowVect(); ng_vect[2] = 0; - if (ivar == MetGridBdyVars::U) { - amrex::Print() << "fill_from_metgrid U var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << " ng_vect=" << ng_vect << std::endl; - } else if (ivar == MetGridBdyVars::V) { - amrex::Print() << "fill_from_metgrid V var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << " ng_vect=" << ng_vect << std::endl; - } else if (ivar == MetGridBdyVars::R) { - amrex::Print() << "fill_from_metgrid R var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << " ng_vect=" << ng_vect << std::endl; - } else if (ivar == MetGridBdyVars::T) { - amrex::Print() << "fill_from_metgrid T var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << " ng_vect=" << ng_vect << std::endl; - } else if (ivar == MetGridBdyVars::QV) { - amrex::Print() << "fill_from_metgrid QV var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << " ng_vect=" << ng_vect << std::endl; - } else { - amrex::Print() << "fill_from_metgrid UNKNOWN" << std::endl; - } - amrex::Print() << " n_time=" << n_time << " \toma=" << oma << " \talpha=" << alpha << std::endl; - // We have data at fixed time intervals we will call dT // Then to interpolate, given time, we can define n = (time/dT) // and alpha = (time - n*dT) / dT, then we define the data at time @@ -118,35 +106,6 @@ ERF::fill_from_metgrid (const Vector& mfs, bx_xlo, bx_xhi, bx_ylo, bx_yhi, ng_vect); - amrex::AllPrint() << "fill_from_metgrid bx_xlo=" << bx_xlo << std::endl; - amrex::AllPrint() << "fill_from_metgrid bx_xhi=" << bx_xhi << std::endl; - amrex::AllPrint() << "fill_from_metgrid bx_ylo=" << bx_ylo << std::endl; - amrex::AllPrint() << "fill_from_metgrid bx_yhi=" << bx_yhi << std::endl; - if (ivar == MetGridBdyVars::QV) { - amrex::Print() << " QV BEFORE" << std::endl; - int jj = (dom_hi.y - dom_lo.y)/2+dom_lo.y; - for (int kk = 9; kk >= 0; kk--) { - for (int ii = 0; ii < 10; ii++) { - amrex::Print() << " " << dest_arr(ii,jj,kk,comp_idx); - } - amrex::Print() << std::endl; - } - amrex::Print() << " bdatxlo_n" << std::endl; - for (int k = 9; k >= 0; k--) { - for (int ii = std::max(dom_lo.x, amrex::lbound(bx_xlo).x); ii <= std::min(dom_hi.x, amrex::ubound(bx_xlo).x); ii++) { - amrex::Print() << " " << bdatxlo_n(ii,jj,k,0); - } - amrex::Print() << std::endl; - } - amrex::Print() << " bdatxlo_np1" << std::endl; - for (int k = 9; k >= 0; k--) { - for (int ii = std::max(dom_lo.x, amrex::lbound(bx_xlo).x); ii <= std::min(dom_hi.x, amrex::ubound(bx_xlo).x); ii++) { - amrex::Print() << " " << bdatxlo_np1(ii,jj,k,0); - } - amrex::Print() << std::endl; - } - } - // x-faces (includes exterior y ghost cells) ParallelFor(bx_xlo, bx_xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) @@ -180,17 +139,6 @@ ERF::fill_from_metgrid (const Vector& mfs, dest_arr(i,j,k,comp_idx) = oma * bdatyhi_n (i,jj,k,0) + alpha * bdatyhi_np1(i,jj,k,0); }); - - if (ivar == MetGridBdyVars::QV) { - amrex::Print() << " QV AFTER" << std::endl; - int jj = (dom_hi.y - dom_lo.y)/2+dom_lo.y; - for (int kk = 9; kk >= 0; kk--) { - for (int ii = 0; ii < 10; ii++) { - amrex::Print() << " " << dest_arr(ii,jj,kk,comp_idx); - } - amrex::Print() << std::endl; - } - } } // mfi // Variable not read or computed from met_em files diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index faae8a85f..f387b6de5 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -17,19 +17,19 @@ ERF::init_from_metgrid (int lev) { bool use_moisture = (solverChoice.moisture_type != MoistureType::None); if (use_moisture) { - amrex::Print() << "Init with met_em with valid moisture model." << std::endl; + Print() << "Init with met_em with valid moisture model." << std::endl; } else { - amrex::Print() << "Init with met_em without moisture model." << std::endl; + Print() << "Init with met_em without moisture model." << std::endl; } int nboxes = num_boxes_at_level[lev]; int ntimes = num_files_at_level[lev]; if (nc_init_file.empty()) - amrex::Error("NetCDF initialization file name must be provided via input"); + Error("NetCDF initialization file name must be provided via input"); if (nc_init_file[lev].empty()) - amrex::Error("NetCDF initialization file name must be provided via input"); + Error("NetCDF initialization file name must be provided via input"); // At least two met_em files are necessary to calculate tendency terms. AMREX_ALWAYS_ASSERT(ntimes >= 2); @@ -72,7 +72,7 @@ ERF::init_from_metgrid (int lev) for (int it = 0; it < ntimes; it++) { #ifndef AMREX_USE_GPU - amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << ": nc_init_file[" << lev << "][" << it << "]\t" << nc_init_file[lev][it] << std::endl; + AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << ": nc_init_file[" << lev << "][" << it << "]\t" << nc_init_file[lev][it] << std::endl; #endif read_from_metgrid(lev, boxes_at_level[lev][0], nc_init_file[lev][it], NC_dateTime[it], NC_epochTime[it], @@ -85,18 +85,18 @@ ERF::init_from_metgrid (int lev) NC_MSFU_fab[it], NC_MSFV_fab[it], NC_MSFM_fab[it], NC_sst_fab[it], NC_lmask_iab[it]); #ifndef AMREX_USE_GPU - amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": flag_psfc \t" << flag_psfc[it] << std::endl; - amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": flag_msfu \t" << flag_msfu[it] << std::endl; - amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": flag_msfv \t" << flag_msfv[it] << std::endl; - amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": flag_msfm \t" << flag_msfm[it] << std::endl; - amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": flag_hgt \t" << flag_hgt[it] << std::endl; - amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": NC_nx \t" << NC_nx[it] << std::endl; - amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": NC_ny \t" << NC_ny[it] << std::endl; - amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": NC_dx \t" << NC_dx[it] << std::endl; - amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": NC_dy \t" << NC_dy[it] << std::endl; - amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": NC_epochTime\t" << NC_epochTime[it] << std::endl; + AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": flag_psfc \t" << flag_psfc[it] << std::endl; + AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": flag_msfu \t" << flag_msfu[it] << std::endl; + AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": flag_msfv \t" << flag_msfv[it] << std::endl; + AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": flag_msfm \t" << flag_msfm[it] << std::endl; + AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": flag_hgt \t" << flag_hgt[it] << std::endl; + AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": NC_nx \t" << NC_nx[it] << std::endl; + AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": NC_ny \t" << NC_ny[it] << std::endl; + AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": NC_dx \t" << NC_dx[it] << std::endl; + AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": NC_dy \t" << NC_dy[it] << std::endl; + AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": NC_epochTime\t" << NC_epochTime[it] << std::endl; // NC_dateTime is only on the IOProcessor. - amrex::Print() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << ": NC_dateTime \t" << NC_dateTime[it] << std::endl; + Print() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << ": NC_dateTime \t" << NC_dateTime[it] << std::endl; #endif } // it @@ -115,9 +115,9 @@ ERF::init_from_metgrid (int lev) for (int it = 1; it < ntimes; it++) { Real NC_dt = NC_epochTime[it]-NC_epochTime[it-1]; #ifndef AMREX_USE_GPU - amrex::Print() << " " << nc_init_file[lev][it-1] << " / " << nc_init_file[lev][it] << " are " << NC_dt << " seconds apart" << std::endl; + Print() << " " << nc_init_file[lev][it-1] << " / " << nc_init_file[lev][it] << " are " << NC_dt << " seconds apart" << std::endl; #endif - if (NC_dt != bdy_time_interval) amrex::Error("Time interval between consecutive met_em files must be consistent."); + if (NC_dt != bdy_time_interval) Error("Time interval between consecutive met_em files must be consistent."); } // Set up a FAB for mixing ratio and another for potential temperature. @@ -178,8 +178,8 @@ ERF::init_from_metgrid (int lev) const Array4& src_arr = src.const_array(); ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { - int li = amrex::min(amrex::max(i, i_lo), i_hi); - int lj = amrex::min(amrex::max(j, j_lo), j_hi); + int li = min(max(i, i_lo), i_hi); + int lj = min(max(j, j_lo), j_hi); dst_arr(i,j,0) = src_arr(li,lj,0); }); } @@ -199,8 +199,8 @@ ERF::init_from_metgrid (int lev) const Array4& src_arr = src.const_array(); ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { - int li = amrex::min(amrex::max(i, i_lo), i_hi); - int lj = amrex::min(amrex::max(j, j_lo), j_hi); + int li = min(max(i, i_lo), i_hi); + int lj = min(max(j, j_lo), j_hi); dst_arr(i,j,0) = src_arr(li,lj,0); }); } @@ -231,7 +231,7 @@ ERF::init_from_metgrid (int lev) if (use_moisture) MetGridBdyEnd = MetGridBdyVars::NumTypes; // Zero out fabs_for_bcs on the global domain - amrex::Vector> fabs_for_bcs; + Vector> fabs_for_bcs; fabs_for_bcs.resize(ntimes); for (int it(0); it < ntimes; it++) { fabs_for_bcs[it].resize(MetGridBdyEnd); @@ -243,16 +243,16 @@ ERF::init_from_metgrid (int lev) ldomain = geom[lev].Domain(); ldomain.convert(IntVect(1,0,0)); auto ng = lev_new[Vars::xvel].nGrowVect(); - gdomain = amrex::grow(ldomain,ng); + gdomain = grow(ldomain,ng); } else if (nvar==MetGridBdyVars::V) { ldomain = geom[lev].Domain(); ldomain.convert(IntVect(0,1,0)); auto ng = lev_new[Vars::yvel].nGrowVect(); - gdomain = amrex::grow(ldomain,ng); + gdomain = grow(ldomain,ng); } else { ldomain = geom[lev].Domain(); auto ng = lev_new[Vars::cons].nGrowVect(); - gdomain = amrex::grow(ldomain,ng); + gdomain = grow(ldomain,ng); } #ifdef AMREX_USE_GPU fabs_for_bcs[it][nvar].resize(gdomain, 1, The_Pinned_Arena()); @@ -357,20 +357,9 @@ ERF::init_from_metgrid (int lev) // complete data set; initialized to 0 above. for (int it(0); it < ntimes; it++) { for (int nvar(0); nvar& fabs_for_bcs_arr = fabs_for_bcs[it][MetGridBdyVars::QV].const_array(); - amrex::Print() << " fabs_for_bcs[" << it << "][" << MetGridBdyVars::QV << "] after ParallelAllReduce::Sum" << std::endl; - for (int k=9; k >= 0; k--) { - for (int ii=0; ii <= 9; ii++) { - amrex::Print() << " " << fabs_for_bcs_arr(ii, 40, k); - } - amrex::Print() << std::endl; + ParallelAllReduce::Sum(fabs_for_bcs[it][nvar].dataPtr(), + fabs_for_bcs[it][nvar].size(), + ParallelContext::CommunicatorAll()); } } @@ -382,8 +371,8 @@ ERF::init_from_metgrid (int lev) const auto& lo = geom[lev].Domain().loVect(); const auto& hi = geom[lev].Domain().hiVect(); - amrex::IntVect plo(lo); - amrex::IntVect phi(hi); + IntVect plo(lo); + IntVect phi(hi); plo[0] = lo[0]; plo[1] = lo[1]; plo[2] = lo[2]; phi[0] = lo[0]+metgrid_bdy_width-1; phi[1] = hi[1]; phi[2] = hi[2]; @@ -414,65 +403,60 @@ ERF::init_from_metgrid (int lev) Box yhi_plane_y_stag = pbx_yhi; yhi_plane_y_stag.shiftHalf(1,1); #ifndef AMREX_USE_GPU - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tplo \t" << plo << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tphi \t" << phi << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane_no_stag \t" << xlo_plane_no_stag << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane_x_stag \t" << xlo_plane_x_stag << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane_y_stag \t" << xlo_plane_y_stag << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tplo \t" << plo << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tphi \t" << phi << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane_no_stag \t" << xhi_plane_no_stag << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane_x_stag \t" << xhi_plane_x_stag << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane_y_stag \t" << xhi_plane_y_stag << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tplo \t" << plo << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tphi \t" << phi << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane_no_stag \t" << ylo_plane_no_stag << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane_x_stag \t" << ylo_plane_x_stag << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane_y_stag \t" << ylo_plane_y_stag << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tplo \t" << plo << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tphi \t" << phi << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane_no_stag \t" << yhi_plane_no_stag << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane_x_stag \t" << yhi_plane_x_stag << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane_y_stag \t" << yhi_plane_y_stag << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tplo \t" << plo << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tphi \t" << phi << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane_no_stag \t" << xlo_plane_no_stag << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane_x_stag \t" << xlo_plane_x_stag << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane_y_stag \t" << xlo_plane_y_stag << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tplo \t" << plo << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tphi \t" << phi << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane_no_stag \t" << xhi_plane_no_stag << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane_x_stag \t" << xhi_plane_x_stag << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane_y_stag \t" << xhi_plane_y_stag << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tplo \t" << plo << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tphi \t" << phi << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane_no_stag \t" << ylo_plane_no_stag << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane_x_stag \t" << ylo_plane_x_stag << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane_y_stag \t" << ylo_plane_y_stag << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tplo \t" << plo << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tphi \t" << phi << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane_no_stag \t" << yhi_plane_no_stag << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane_x_stag \t" << yhi_plane_x_stag << std::endl; + AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane_y_stag \t" << yhi_plane_y_stag << std::endl; #endif for (int ivar(MetGridBdyVars::U); ivar < MetGridBdyEnd; ivar++) { for (int it(0); it < ntimes; it++) { if (ivar == MetGridBdyVars::U) { - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tbdy_data_*[" << it << "] for U" << std::endl; bdy_data_xlo[it].push_back(FArrayBox(xlo_plane_x_stag, 1)); bdy_data_xhi[it].push_back(FArrayBox(xhi_plane_x_stag, 1)); bdy_data_ylo[it].push_back(FArrayBox(ylo_plane_x_stag, 1)); bdy_data_yhi[it].push_back(FArrayBox(yhi_plane_x_stag, 1)); } else if (ivar == MetGridBdyVars::V) { - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tbdy_data_*[" << it << "] for V" << std::endl; bdy_data_xlo[it].push_back(FArrayBox(xlo_plane_y_stag, 1)); bdy_data_xhi[it].push_back(FArrayBox(xhi_plane_y_stag, 1)); bdy_data_ylo[it].push_back(FArrayBox(ylo_plane_y_stag, 1)); bdy_data_yhi[it].push_back(FArrayBox(yhi_plane_y_stag, 1)); } else if (ivar == MetGridBdyVars::R) { - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tbdy_data_*[" << it << "] for R" << std::endl; bdy_data_xlo[it].push_back(FArrayBox(xlo_plane_no_stag, 1)); bdy_data_xhi[it].push_back(FArrayBox(xhi_plane_no_stag, 1)); bdy_data_ylo[it].push_back(FArrayBox(ylo_plane_no_stag, 1)); bdy_data_yhi[it].push_back(FArrayBox(yhi_plane_no_stag, 1)); } else if (ivar == MetGridBdyVars::T) { - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tbdy_data_*[" << it << "] for T" << std::endl; bdy_data_xlo[it].push_back(FArrayBox(xlo_plane_no_stag, 1)); bdy_data_xhi[it].push_back(FArrayBox(xhi_plane_no_stag, 1)); bdy_data_ylo[it].push_back(FArrayBox(ylo_plane_no_stag, 1)); bdy_data_yhi[it].push_back(FArrayBox(yhi_plane_no_stag, 1)); } else if (ivar == MetGridBdyVars::QV) { - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tbdy_data_*[" << it << "] for QV" << std::endl; bdy_data_xlo[it].push_back(FArrayBox(xlo_plane_no_stag, 1)); bdy_data_xhi[it].push_back(FArrayBox(xhi_plane_no_stag, 1)); bdy_data_ylo[it].push_back(FArrayBox(ylo_plane_no_stag, 1)); bdy_data_yhi[it].push_back(FArrayBox(yhi_plane_no_stag, 1)); } else { #ifndef AMREX_USE_GPU - amrex::Print() << "Unexpected ivar " << ivar << std::endl; + Print() << "Unexpected ivar " << ivar << std::endl; #endif - amrex::Abort("See Initialization/ERF_init_from_metgrid.cpp"); + Abort("See Initialization/ERF_init_from_metgrid.cpp"); } } // it } // ivar @@ -482,7 +466,7 @@ ERF::init_from_metgrid (int lev) // at subsequent times. We can optimize this later if needed. For now, we need to fill // the lateral boundary arrays using the info set aside earlier. bool multiply_rho = false; - amrex::Box xlo_plane, xhi_plane, ylo_plane, yhi_plane; + Box xlo_plane, xhi_plane, ylo_plane, yhi_plane; for (int it(0); it < ntimes; it++) { for (int ivar(MetGridBdyVars::U); ivar < MetGridBdyEnd; ivar++) { @@ -501,42 +485,22 @@ ERF::init_from_metgrid (int lev) multiply_rho = false; xlo_plane = xlo_plane_x_stag; xhi_plane = xhi_plane_x_stag; ylo_plane = ylo_plane_x_stag; yhi_plane = yhi_plane_x_stag; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane " << xlo_plane << " \tfor U" << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane " << xhi_plane << " \tfor U" << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane " << ylo_plane << " \tfor U" << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane " << yhi_plane << " \tfor U" << std::endl; } else if (ivar == MetGridBdyVars::V) { multiply_rho = false; xlo_plane = xlo_plane_y_stag; xhi_plane = xhi_plane_y_stag; ylo_plane = ylo_plane_y_stag; yhi_plane = yhi_plane_y_stag; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane " << xlo_plane << " \tfor V" << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane " << xhi_plane << " \tfor V" << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane " << ylo_plane << " \tfor V" << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane " << yhi_plane << " \tfor V" << std::endl; } else if (ivar == MetGridBdyVars::R) { multiply_rho = false; xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag; ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane " << xlo_plane << " \tfor R" << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane " << xhi_plane << " \tfor R" << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane " << ylo_plane << " \tfor R" << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane " << yhi_plane << " \tfor R" << std::endl; } else if (ivar == MetGridBdyVars::T) { multiply_rho = false; xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag; ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane " << xlo_plane << " \tfor T" << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane " << xhi_plane << " \tfor T" << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane " << ylo_plane << " \tfor T" << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane " << yhi_plane << " \tfor T" << std::endl; } else if (ivar == MetGridBdyVars::QV) { multiply_rho = false; xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag; ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txlo_plane " << xlo_plane << " \tfor QV" << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \txhi_plane " << xhi_plane << " \tfor QV" << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tylo_plane " << ylo_plane << " \tfor QV" << std::endl; - amrex::AllPrint() << " proc-" << ParallelDescriptor::MyProc() << " \tyhi_plane " << yhi_plane << " \tfor QV" << std::endl; } // MetGridBdyVars::QV // west boundary @@ -578,8 +542,8 @@ init_terrain_from_metgrid (FArrayBox& z_phys_nd_fab, for (int it = 0; it < ntimes; it++) { #ifndef AMREX_USE_GPU - amrex::Print() << " SIZE OF HGT FAB " << NC_hgt_fab[it].box() << std::endl; - amrex::Print() << " SIZE OF ZP FAB " << z_phys_nd_fab.box() << std::endl; + Print() << " SIZE OF HGT FAB " << NC_hgt_fab[it].box() << std::endl; + Print() << " SIZE OF ZP FAB " << z_phys_nd_fab.box() << std::endl; #endif // This copies from NC_zphys on z-faces to z_phys_nd on nodes @@ -602,8 +566,8 @@ init_terrain_from_metgrid (FArrayBox& z_phys_nd_fab, Box bxv = bx; bxv.growHi(1,1); #ifndef AMREX_USE_GPU - amrex::Print() << "FROM BOX " << from_box << std::endl; - amrex::Print() << "BX " << bx << std::endl; + Print() << "FROM BOX " << from_box << std::endl; + Print() << "BX " << bx << std::endl; #endif ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -657,10 +621,10 @@ init_state_from_metgrid (const bool use_moisture, const Vector& NC_pres_fab, Vector& theta_fab, Vector& mxrat_fab, - amrex::Vector>& fabs_for_bcs, - const amrex::Array4& mask_c_arr, - const amrex::Array4& mask_u_arr, - const amrex::Array4& mask_v_arr) + Vector>& fabs_for_bcs, + const Array4& mask_c_arr, + const Array4& mask_u_arr, + const Array4& mask_v_arr) { int ntimes = NC_hgt_fab.size(); for (int it = 0; it < ntimes; it++) @@ -677,7 +641,7 @@ init_state_from_metgrid (const bool use_moisture, auto bc_data = fabs_for_bcs[it][MetGridBdyVars::U].array(); auto const new_z = z_phys_nd_fab.const_array(); - int kmax = amrex::ubound(tbxu).z; + int kmax = ubound(tbxu).z; ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { @@ -687,13 +651,6 @@ init_state_from_metgrid (const bool use_moisture, if (it==0) new_data(i,j,k,0) = Interp_Val; } }); - amrex::Print() << " mask_u_arr" << std::endl; - for (int k=9; k >= 0; k--) { - for (int ii=0; ii >= 3; ii++) { - amrex::Print() << " " << mask_u_arr(ii,40,k); - } - amrex::Print() << std::endl; - } } @@ -709,7 +666,7 @@ init_state_from_metgrid (const bool use_moisture, auto bc_data = fabs_for_bcs[it][MetGridBdyVars::V].array(); auto const new_z = z_phys_nd_fab.const_array(); - int kmax = amrex::ubound(tbxv).z; + int kmax = ubound(tbxv).z; ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { @@ -771,7 +728,7 @@ init_state_from_metgrid (const bool use_moisture, auto bc_data = fabs_for_bcs[it][MetGridBdyVars::T].array(); auto const new_z = z_phys_nd_fab.const_array(); - int kmax = amrex::ubound(tbxc).z; + int kmax = ubound(tbxc).z; ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { @@ -815,7 +772,7 @@ init_state_from_metgrid (const bool use_moisture, auto bc_data = fabs_for_bcs[it][MetGridBdyVars::QV].array(); auto const new_z = z_phys_nd_fab.const_array(); - int kmax = amrex::ubound(tbxc).z; + int kmax = ubound(tbxc).z; int state_indx = RhoQ1_comp; ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept @@ -858,7 +815,7 @@ init_state_from_metgrid (const bool use_moisture, auto bc_data = fabs_for_bcs[it][MetGridBdyVars::T].array(); auto const new_z = z_phys_nd_fab.const_array(); - int kmax = amrex::ubound(tbxc).z; + int kmax = ubound(tbxc).z; ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { @@ -902,20 +859,19 @@ init_base_state_from_metgrid (const bool use_moisture, const Vector& NC_ght_fab, const Vector& NC_psfc_fab, Vector>& fabs_for_bcs, - const amrex::Array4& mask_c_arr) + const Array4& mask_c_arr) { int RhoQ_comp = RhoQ1_comp; - int kmax = amrex::ubound(valid_bx).z; + int kmax = ubound(valid_bx).z; // Device vectors for columnwise operations - Gpu::DeviceVector z_vec_d(kmax+2,0); Real* z_vec = z_vec_d.data(); - Gpu::DeviceVector Thetad_vec_d(kmax+1,0); Real* Thetad_vec = Thetad_vec_d.data(); - Gpu::DeviceVector Thetam_vec_d(kmax+1,0); Real* Thetam_vec = Thetam_vec_d.data(); - Gpu::DeviceVector Rhod_vec_d(kmax+1,0); Real* Rhod_vec = Rhod_vec_d.data(); - Gpu::DeviceVector Rhom_vec_d(kmax+1,0); Real* Rhom_vec = Rhom_vec_d.data(); - Gpu::DeviceVector Pd_vec_d(kmax+1,0); Real* Pd_vec = Pd_vec_d.data(); - Gpu::DeviceVector Pm_vec_d(kmax+1,0); Real* Pm_vec = Pm_vec_d.data(); - Gpu::DeviceVector Q_vec_d(kmax+1,0); Real* Q_vec = Q_vec_d.data(); + Gpu::DeviceVector z_vec_d(kmax+2,0); Real* z_vec = z_vec_d.data(); + Gpu::DeviceVector Theta_vec_d(kmax+1,0); Real* Theta_vec = Theta_vec_d.data(); + Gpu::DeviceVector Rhod_vec_d(kmax+1,0); Real* Rhod_vec = Rhod_vec_d.data(); + Gpu::DeviceVector Rhom_vec_d(kmax+1,0); Real* Rhom_vec = Rhom_vec_d.data(); + Gpu::DeviceVector Pd_vec_d(kmax+1,0); Real* Pd_vec = Pd_vec_d.data(); + Gpu::DeviceVector Pm_vec_d(kmax+1,0); Real* Pm_vec = Pm_vec_d.data(); + Gpu::DeviceVector Q_vec_d(kmax+1,0); Real* Q_vec = Q_vec_d.data(); // Device vectors for psfc flags Gpu::DeviceVectorflag_psfc_d(flag_psfc.size()); @@ -940,15 +896,14 @@ init_base_state_from_metgrid (const bool use_moisture, ParallelFor(valid_bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { for (int k=0; k<=kmax; k++) { - z_vec[k] = new_z(i,j,k); - Thetad_vec[k] = new_data(i,j,k,RhoTheta_comp); - Q_vec[k] = (use_moisture) ? new_data(i,j,k,RhoQ_comp) : 0.0; + z_vec[k] = new_z(i,j,k); + Theta_vec[k] = new_data(i,j,k,RhoTheta_comp); + Q_vec[k] = (use_moisture) ? new_data(i,j,k,RhoQ_comp) : 0.0; } z_vec[kmax+1] = new_z(i,j,kmax+1); - calc_rho_p(kmax,flag_psfc_vec[0],orig_psfc(i,j,0),Thetad_vec,Thetam_vec, - Q_vec, - z_vec,Rhod_vec,Rhom_vec,Pd_vec,Pm_vec); + calc_rho_p(l_rdOcp,kmax,flag_psfc_vec[0],orig_psfc(i,j,0),Theta_vec, + Q_vec,z_vec,Rhod_vec,Rhom_vec,Pd_vec,Pm_vec); for (int k=0; k<=kmax; k++) { p_hse_arr(i,j,k) = Pd_vec[k]; @@ -998,22 +953,21 @@ init_base_state_from_metgrid (const bool use_moisture, ParallelFor(valid_bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { for (int k=0; k<=kmax; k++) { - z_vec[k] = new_z(i,j,k); - Thetad_vec[k] = Theta_arr(i,j,k); + z_vec[k] = new_z(i,j,k); + Theta_vec[k] = Theta_arr(i,j,k); Q_vec[k] = (use_moisture) ? Q_arr(i,j,k) : 0.0; } z_vec[kmax+1] = new_z(i,j,kmax+1); - calc_rho_p(kmax,flag_psfc_vec[it],orig_psfc(i,j,0),Thetad_vec,Thetam_vec, - Q_vec, - z_vec,Rhod_vec,Rhom_vec,Pd_vec,Pm_vec); + calc_rho_p(l_rdOcp,kmax,flag_psfc_vec[it],orig_psfc(i,j,0),Theta_vec, + Q_vec,z_vec,Rhod_vec,Rhom_vec,Pd_vec,Pm_vec); for (int k=0; k<=kmax; k++) { p_hse_arr(i,j,k) = Pd_vec[k]; if (mask_c_arr(i,j,k)) { r_hse_arr(i,j,k) = Rhod_vec[k]; - Q_arr(i,j,k) = (use_moisture) ? Rhod_vec[k]*Q_vec[k] : 0.0; - Theta_arr(i,j,k) = Rhod_vec[k]*Thetad_vec[k]; + Q_arr(i,j,k) = (use_moisture) ? Rhod_vec[k]*Q_vec[k] : 0.0; + Theta_arr(i,j,k) = Rhod_vec[k]*Theta_vec[k]; } } // k }); @@ -1058,7 +1012,7 @@ init_msfs_from_metgrid (FArrayBox& msfu_fab, msfm_fab.template copy(NC_MSFM_fab[it]); } else { #ifndef AMREX_USE_GPU - amrex::Print() << " MAPFAC_M not present in met_em files. Setting to 1.0" << std::endl; + Print() << " MAPFAC_M not present in met_em files. Setting to 1.0" << std::endl; #endif msfm_fab.template setVal(1.0); } @@ -1068,7 +1022,7 @@ init_msfs_from_metgrid (FArrayBox& msfu_fab, msfu_fab.template copy(NC_MSFU_fab[it]); } else { #ifndef AMREX_USE_GPU - amrex::Print() << " MAPFAC_U not present in met_em files. Setting to 1.0" << std::endl; + Print() << " MAPFAC_U not present in met_em files. Setting to 1.0" << std::endl; #endif msfu_fab.template setVal(1.0); } @@ -1078,7 +1032,7 @@ init_msfs_from_metgrid (FArrayBox& msfu_fab, msfv_fab.template copy(NC_MSFV_fab[it]); } else { #ifndef AMREX_USE_GPU - amrex::Print() << " MAPFAC_V not present in met_em files. Setting to 1.0" << std::endl; + Print() << " MAPFAC_V not present in met_em files. Setting to 1.0" << std::endl; #endif msfv_fab.template setVal(1.0); } diff --git a/Source/Initialization/Metgrid_utils.H b/Source/Initialization/Metgrid_utils.H index 4e2cc6ac3..c646b2209 100644 --- a/Source/Initialization/Metgrid_utils.H +++ b/Source/Initialization/Metgrid_utils.H @@ -98,7 +98,7 @@ init_base_state_from_metgrid (const bool use_moisture, AMREX_FORCE_INLINE AMREX_GPU_DEVICE void -calc_rho_p (const amrex::Real& rdOcp, +calc_rho_p (const amrex::Real& l_rdOcp, const int& kmax, const int& flag_psfc, const amrex::Real& psfc, @@ -122,7 +122,7 @@ calc_rho_p (const amrex::Real& rdOcp, } // calculate moist density at the surface. - Rhom_vec[0] = getRhogivenThetaPress(Theta_vec[0], Pm_vec[0], rdOcp, Q)*(1.0+Q); + Rhom_vec[0] = getRhogivenThetaPress(Theta_vec[0], Pm_vec[0], l_rdOcp, Q_vec[0])*(1.0+Q_vec[0]); // integrate from the surface to the top boundary. for (int k=1; k<=kmax; k++) { @@ -131,35 +131,20 @@ calc_rho_p (const amrex::Real& rdOcp, for (int it=0; it=0; k--) { amrex::Real dz = z_vec[k+1]-z_vec[k]; Rhod_vec[k] = Rhod_vec[k+1]; // an initial guess. for (int it=0; it < maxiter; it++) { Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*CONST_GRAV; if (Pd_vec[k] < 0.0) Pd_vec[k] = 0.0; -#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) - Q = Q_vec[k]; -#else - Q = 0.0 -#endif - Rhod_vec[k] = getRhogivenThetaPress(Theta_vec[k], Pd_vec[k], rdOcp, Q); + Rhod_vec[k] = getRhogivenThetaPress(Theta_vec[k], Pd_vec[k], l_rdOcp, Q_vec[k]); } // it } // k } diff --git a/Source/Utils/InteriorGhostCells.cpp b/Source/Utils/InteriorGhostCells.cpp index f1d0148cb..d88affd04 100644 --- a/Source/Utils/InteriorGhostCells.cpp +++ b/Source/Utils/InteriorGhostCells.cpp @@ -437,8 +437,6 @@ wrfbdy_compute_interior_ghost_rhs (const std::string& init_type, } } // mfi -// return; // DJW debugging return statement to shut off relaxation zone. - // Compute RHS in relaxation region //========================================================== for (int ivar(ivarU); ivar < BdyEnd; ivar++) From a51636c941a44568a55f09aa0b677ca0ae48cee6 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Tue, 23 Jan 2024 14:59:08 -0800 Subject: [PATCH 5/5] WIP debugging metgrid initialization and forcing --- Exec/DevTests/MetGrid/inputs | 8 +-- .../Initialization/ERF_init_from_metgrid.cpp | 50 +++++++++++-------- Source/Initialization/Metgrid_utils.H | 14 ++++-- 3 files changed, 45 insertions(+), 27 deletions(-) diff --git a/Exec/DevTests/MetGrid/inputs b/Exec/DevTests/MetGrid/inputs index 31653e86b..89708d118 100644 --- a/Exec/DevTests/MetGrid/inputs +++ b/Exec/DevTests/MetGrid/inputs @@ -1,5 +1,5 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 1 +max_step = 10 amrex.fpe_trap_invalid = 1 amrex.fpe_trap_zero = 1 @@ -41,6 +41,7 @@ erf.check_int = 100 # number of timesteps between checkpoints erf.plot_file_1 = plt # prefix of plotfile name erf.plot_int_1 = 1 # number of timesteps between plotfiles erf.plot_vars_1 = qt qv qc density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres pert_dens KE QKE +#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 pert_dens KE QKE # SOLVER CHOICE erf.alpha_T = 1.0 @@ -53,14 +54,15 @@ erf.les_type = "Deardorff" erf.Cs = 0.1 erf.moisture_model = "Kessler" +#erf.moisture_model = "Null" #erf.do_cloud = false #erf.do_precip = false #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.metgrid_bdy_width = 5 -erf.metgrid_bdy_set_width = 5 +erf.metgrid_bdy_width = 10 +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" diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index c20e9607b..36a195a55 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -479,7 +479,10 @@ ERF::init_from_metgrid (int lev) auto xhi_arr = bdy_data_xhi[it][ivar].array(); auto ylo_arr = bdy_data_ylo[it][ivar].array(); auto yhi_arr = bdy_data_yhi[it][ivar].array(); - const Array4& fabs_for_bcs_arr = fabs_for_bcs[it][ivar].const_array(); +// const Array4& fabs_for_bcs_arr = fabs_for_bcs[it][ivar].const_array(); + + // TODO: stationary boundary conditions, remove or comment when not needed for debugging. + const Array4& fabs_for_bcs_arr = fabs_for_bcs[0][ivar].const_array(); if (ivar == MetGridBdyVars::U) { xlo_plane = xlo_plane_x_stag; xhi_plane = xhi_plane_x_stag; @@ -675,12 +678,12 @@ init_state_from_metgrid (const bool use_moisture, // --------------------------------------------------------------------------------------- // DJW: remove this after debugging is complete. This enforces quiescent conditions.. - if (it == 0) { - x_vel_fab.template setVal(0.0); - y_vel_fab.template setVal(0.0); - } - fabs_for_bcs[it][MetGridBdyVars::U].template setVal(0.0); - fabs_for_bcs[it][MetGridBdyVars::V].template setVal(0.0); +// if (it == 0) { +// x_vel_fab.template setVal(0.0); +// y_vel_fab.template setVal(0.0); +// } +// fabs_for_bcs[it][MetGridBdyVars::U].template setVal(0.0); +// fabs_for_bcs[it][MetGridBdyVars::V].template setVal(0.0); // --------------------------------------------------------------------------------------- // ******************************************************** @@ -782,12 +785,12 @@ init_state_from_metgrid (const bool use_moisture, } // TODO: TEMPORARY CODE TO RUN QUIESCENT, REMOVE WHEN NOT NEEDED. -// if (it == 0) { -// x_vel_fab.template setVal(0.0); // TODO: temporary code to initialize with quiescent atmosphere. -// y_vel_fab.template setVal(0.0); // TODO: temporary code to initialize with quiescent atmosphere. -// } -// fabs_for_bcs[it][MetGridBdyVars::U].template setVal(0.0); // TODO: temporary code to force with quiescent atmosphere. -// fabs_for_bcs[it][MetGridBdyVars::V].template setVal(0.0); // TODO: temporary code to force with quiescent atmosphere. + if (it == 0) { + x_vel_fab.template setVal(0.0); // TODO: temporary code to initialize with quiescent atmosphere. + y_vel_fab.template setVal(0.0); // TODO: temporary code to initialize with quiescent atmosphere. + } + fabs_for_bcs[it][MetGridBdyVars::U].template setVal(0.0); // TODO: temporary code to force with quiescent atmosphere. + fabs_for_bcs[it][MetGridBdyVars::V].template setVal(0.0); // TODO: temporary code to force with quiescent atmosphere. } // it } @@ -863,29 +866,35 @@ init_base_state_from_metgrid (const bool use_moisture, Q_vec[k] = (use_moisture) ? new_data(i,j,k,RhoQ_comp) : 0.0; } z_vec[kmax+1] = new_z(i,j,kmax+1); + if ((i == 40) && (j == 40)) { + for (int k=0; k<=kmax; k++) { + Print() << " z_vec[" << k << "] = " << z_vec[k] << std::endl; + } + } calc_rho_p(l_rdOcp,kmax,flag_psfc_vec[0],orig_psfc(i,j,0),Theta_vec, Q_vec,z_vec,Rhod_vec,Rhom_vec,Pd_vec,Pm_vec); for (int k=0; k<=kmax; k++) { p_hse_arr(i,j,k) = Pd_vec[k]; - r_hse_arr(i,j,k) = Rhom_vec[k]; + r_hse_arr(i,j,k) = Rhod_vec[k]; } }); ParallelFor(valid_bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - new_data(i,j,k,Rho_comp) = r_hse_arr(i,j,k); new_data(i,j,k,RhoScalar_comp) = 0.0; // RhoTheta and RhoQ1 or RhoQv currently hold Theta and Q1 or Qv. Multiply by Rho. + new_data(i,j,k,Rho_comp) = r_hse_arr(i,j,k); Real RhoTheta = r_hse_arr(i,j,k)*new_data(i,j,k,RhoTheta_comp); new_data(i,j,k,RhoTheta_comp) = RhoTheta; if (use_moisture){ Real RhoQ = r_hse_arr(i,j,k)*new_data(i,j,k,RhoQ_comp); new_data(i,j,k,RhoQ_comp) = RhoQ; } - pi_hse_arr(i,j,k) = getExnergivenP(p_hse_arr(i,j,k), l_rdOcp); - //pi_hse_arr(i,j,k) = getExnergivenRTh(RhoTheta, l_rdOcp); +// p_hse_arr(i,j,k) = getPgivenRTh(RhoTheta); + //pi_hse_arr(i,j,k) = getExnergivenP(p_hse_arr(i,j,k), l_rdOcp); + pi_hse_arr(i,j,k) = getExnergivenRTh(RhoTheta, l_rdOcp); }); } @@ -926,10 +935,11 @@ init_base_state_from_metgrid (const bool use_moisture, for (int k=0; k<=kmax; k++) { p_hse_arr(i,j,k) = Pd_vec[k]; +// p_hse_arr(i,j,k) = getPgivenRTh(Rhod_vec[k]*Theta_vec[k]); if (mask_c_arr(i,j,k)) { - r_hse_arr(i,j,k) = Rhom_vec[k]; - if (use_moisture) Q_arr(i,j,k) = Rhom_vec[k]*Q_vec[k]; - Theta_arr(i,j,k) = Rhom_vec[k]*Theta_vec[k]; + r_hse_arr(i,j,k) = Rhod_vec[k]; + if (use_moisture) Q_arr(i,j,k) = Rhod_vec[k]*Q_vec[k]; + Theta_arr(i,j,k) = Rhod_vec[k]*Theta_vec[k]; } } // k }); diff --git a/Source/Initialization/Metgrid_utils.H b/Source/Initialization/Metgrid_utils.H index c646b2209..e4cc1410a 100644 --- a/Source/Initialization/Metgrid_utils.H +++ b/Source/Initialization/Metgrid_utils.H @@ -122,30 +122,36 @@ calc_rho_p (const amrex::Real& l_rdOcp, } // calculate moist density at the surface. - Rhom_vec[0] = getRhogivenThetaPress(Theta_vec[0], Pm_vec[0], l_rdOcp, Q_vec[0])*(1.0+Q_vec[0]); +// Rhom_vec[0] = getRhogivenThetaPress(Theta_vec[0], Pm_vec[0], l_rdOcp, Q_vec[0])*(1.0+Q_vec[0]); + Rhom_vec[0] = 1.0/(R_d/p_0*Theta_vec[0]*(1.0+(R_v/R_d)*Q_vec[0])*std::pow(Pm_vec[0]/p_0, -iGamma)); // integrate from the surface to the top boundary. for (int k=1; k<=kmax; k++) { + amrex::Real qvf = 1.0+(R_v/R_d)*Q_vec[k]; amrex::Real dz = z_vec[k]-z_vec[k-1]; Rhom_vec[k] = Rhom_vec[k-1]; // an initial guess. for (int it=0; it=0; k--) { amrex::Real dz = z_vec[k+1]-z_vec[k]; Rhod_vec[k] = Rhod_vec[k+1]; // an initial guess. for (int it=0; it < maxiter; it++) { Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*CONST_GRAV; if (Pd_vec[k] < 0.0) Pd_vec[k] = 0.0; - Rhod_vec[k] = getRhogivenThetaPress(Theta_vec[k], Pd_vec[k], l_rdOcp, Q_vec[k]); +// Rhod_vec[k] = getRhogivenThetaPress(Theta_vec[k], Pd_vec[k], l_rdOcp, Q_vec[k]); + Rhod_vec[k] = 1.0/(R_d/p_0*Theta_vec[k]*std::pow(Pd_vec[k]/p_0, -iGamma)); } // it + Pd_vec[k] = getPgivenRTh(Rhod_vec[k]*Theta_vec[k]); } // k }