From 4157f68b73ddea75613bd10d06dc76a23475333b Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Sun, 11 Feb 2024 11:47:41 -0800 Subject: [PATCH] Balance pressure. (#1427) * Balance pressure. * Overwrite the bdy vars with init state, coalesce the met and wrf operations, fix FArrayBox on GPU without managed memory. * Missed input change. * Add commented out unit test and change inputs to use MOST for bottom BC. * Fill halo cells of HSE arrays for safety. --------- Co-authored-by: Aaron Lattanzi Co-authored-by: Ann Almgren --- CMake/BuildERFExe.cmake | 3 +- Docs/sphinx_doc/BoundaryConditions.rst | 4 +- Exec/DevTests/MetGrid/inputs | 10 +- Exec/DevTests/MetGrid/inputs_flowmas_metgrid | 16 +- .../MetGrid/inputs_flowmas_metgrid_NoZlev | 8 +- Exec/DevTests/MetGrid/inputs_flowmas_wrfinput | 16 +- Exec/LLJ/inputs | 15 +- .../WPS_Test/inputs_real_ChisholmView | 10 +- .../BoundaryConditions_metgrid.cpp | 193 ------------ ...bdy.cpp => BoundaryConditions_realbdy.cpp} | 22 +- Source/BoundaryConditions/ERF_FillPatch.cpp | 13 +- Source/BoundaryConditions/Make.package | 3 +- Source/ERF.H | 25 +- Source/ERF.cpp | 17 +- Source/IO/Checkpoint.cpp | 6 +- Source/IO/NCWpsFile.H | 6 +- Source/IO/ReadFromMetgrid.cpp | 6 +- Source/IO/ReadFromWRFBdy.cpp | 279 +++++++----------- Source/IO/ReadFromWRFInput.cpp | 32 +- Source/IndexDefines.H | 11 + .../Initialization/ERF_init_from_metgrid.cpp | 102 ++----- .../Initialization/ERF_init_from_wrfinput.cpp | 157 ++++++---- Source/TimeIntegration/ERF_slow_rhs_post.cpp | 40 +++ Source/TimeIntegration/TI_slow_rhs_fun.H | 73 ++--- Source/Utils/InteriorGhostCells.cpp | 90 ++++-- Source/Utils/Utils.H | 30 +- 26 files changed, 497 insertions(+), 690 deletions(-) delete mode 100644 Source/BoundaryConditions/BoundaryConditions_metgrid.cpp rename Source/BoundaryConditions/{BoundaryConditions_wrfbdy.cpp => BoundaryConditions_realbdy.cpp} (92%) diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 9e18901e1..a877d1c67 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -100,8 +100,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/BoundaryConditions/BoundaryConditions_yvel.cpp ${SRC_DIR}/BoundaryConditions/BoundaryConditions_zvel.cpp ${SRC_DIR}/BoundaryConditions/BoundaryConditions_bndryreg.cpp - ${SRC_DIR}/BoundaryConditions/BoundaryConditions_wrfbdy.cpp - ${SRC_DIR}/BoundaryConditions/BoundaryConditions_metgrid.cpp + ${SRC_DIR}/BoundaryConditions/BoundaryConditions_realbdy.cpp ${SRC_DIR}/BoundaryConditions/ERF_FillPatch.cpp ${SRC_DIR}/BoundaryConditions/ERF_FillPatcher.cpp ${SRC_DIR}/BoundaryConditions/ERF_PhysBCFunct.cpp diff --git a/Docs/sphinx_doc/BoundaryConditions.rst b/Docs/sphinx_doc/BoundaryConditions.rst index 8d24483b1..d1057b548 100644 --- a/Docs/sphinx_doc/BoundaryConditions.rst +++ b/Docs/sphinx_doc/BoundaryConditions.rst @@ -145,9 +145,9 @@ Specified Domain BCs When we use specified lateral boundary conditions, we read time-dependent Dirichlet data from a file. The user may specify (in the inputs file) the total width of the interior Dirichlet and relaxation region with -``erf.wrfbdy_width = `` (yellow + blue) +``erf.real_width = `` (yellow + blue) and analogously the width of the interior Dirichlet region may be specified with -``erf.wrfbdy_set_width = `` (yellow). +``erf.real_set_width = `` (yellow). We note that all dycore variables are set and relaxed, but moisture and other scalars are only set in the yellow region if present in the boundary file. diff --git a/Exec/DevTests/MetGrid/inputs b/Exec/DevTests/MetGrid/inputs index db926600c..ccf0ee7dc 100644 --- a/Exec/DevTests/MetGrid/inputs +++ b/Exec/DevTests/MetGrid/inputs @@ -15,9 +15,13 @@ xlo.type = "Outflow" xhi.type = "Outflow" ylo.type = "Outflow" yhi.type = "Outflow" -zlo.type = "NoSlipWall" +#zlo.type = "NoSlipWall" zhi.type = "SlipWall" +zlo.type = "Most" +erf.most.z0 = 0.1 +erf.most.zref = 40.0 + # TIME STEP CONTROL erf.fixed_dt = 0.5 # fixed time step depending on grid resolution @@ -56,8 +60,8 @@ erf.moisture_model = "NullMoist" #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 = 7 -erf.metgrid_bdy_set_width = 1 +erf.real_width = 7 +erf.real_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" diff --git a/Exec/DevTests/MetGrid/inputs_flowmas_metgrid b/Exec/DevTests/MetGrid/inputs_flowmas_metgrid index 828d69fba..16059ed5d 100644 --- a/Exec/DevTests/MetGrid/inputs_flowmas_metgrid +++ b/Exec/DevTests/MetGrid/inputs_flowmas_metgrid @@ -15,9 +15,13 @@ xlo.type = "Outflow" xhi.type = "Outflow" ylo.type = "Outflow" yhi.type = "Outflow" -zlo.type = "NoSlipWall" +#zlo.type = "NoSlipWall" zhi.type = "SlipWall" +zlo.type = "Most" +erf.most.z0 = 0.1 +erf.most.zref = 26.63 + # TIME STEP CONTROL erf.fixed_dt = 1.0 # fixed time step depending on grid resolution @@ -31,7 +35,7 @@ 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 +erf.check_int = 1000 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name @@ -60,10 +64,10 @@ erf.terrain_z_levels = 0 53.25373067 116.3299934 190.8583428 278.6665766 381.71 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" +erf.real_width = 5 +erf.real_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_flowmas_metgrid_NoZlev b/Exec/DevTests/MetGrid/inputs_flowmas_metgrid_NoZlev index fcfc67c85..3c00ada13 100644 --- a/Exec/DevTests/MetGrid/inputs_flowmas_metgrid_NoZlev +++ b/Exec/DevTests/MetGrid/inputs_flowmas_metgrid_NoZlev @@ -54,10 +54,10 @@ erf.terrain_smoothing = 2 erf.buoyancy_type = 1 # 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.real_width = 7 +erf.real_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_flowmas_wrfinput b/Exec/DevTests/MetGrid/inputs_flowmas_wrfinput index 19ff3dbd4..b40247a84 100644 --- a/Exec/DevTests/MetGrid/inputs_flowmas_wrfinput +++ b/Exec/DevTests/MetGrid/inputs_flowmas_wrfinput @@ -15,9 +15,13 @@ xlo.type = "Outflow" xhi.type = "Outflow" ylo.type = "Outflow" yhi.type = "Outflow" -zlo.type = "NoSlipWall" +#zlo.type = "NoSlipWall" zhi.type = "SlipWall" +zlo.type = "Most" +erf.most.z0 = 0.1 +erf.most.zref = 26.63 + # TIME STEP CONTROL erf.fixed_dt = 1.0 # fixed time step depending on grid resolution @@ -54,11 +58,11 @@ erf.terrain_smoothing = 0 erf.buoyancy_type = 1 # INITIALIZATION WITH ATM DATA -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" +erf.real_width = 5 +erf.real_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/Exec/LLJ/inputs b/Exec/LLJ/inputs index fc9868508..4973d93ee 100644 --- a/Exec/LLJ/inputs +++ b/Exec/LLJ/inputs @@ -25,12 +25,6 @@ yhi.type = "Outflow" zlo.type = "NoSlipWall" zhi.type = "SlipWall" -# RELAXATION WIDTHS -erf.wrfbdy_width = 5 -erf.wrfbdy_set_width = 1 -erf.cf_width = 5 -erf.cf_set_width = 1 - # TIME STEP CONTROL erf.fixed_dt = 0.001 # fixed time step depending on grid resolution erf.no_substepping = 1 @@ -63,17 +57,20 @@ prob.rho_0 = 1.0 prob.T_0 = 300.0 # MULTILEVEL FLAGS -#amr.max_level = 0 -amr.max_level = 1 +erf.cf_width = 5 +erf.cf_set_width = 1 +amr.max_level = 1 amr.ref_ratio_vect = 3 3 1 erf.refinement_indicators = box1 erf.box1.max_level = 1 erf.box1.in_box_lo = 60. 51. erf.box1.in_box_hi = 525. 309. -#erf.coupling_type = "TwoWay" +#erf.coupling_type = "TwoWay" erf.nc_init_file_1 = "wrfinput_d02" # INITIALIZATION WITH ATM DATA +erf.real_width = 5 +erf.real_set_width = 1 erf.init_type = "real" erf.nc_init_file_0 = "wrfinput_d01" erf.nc_bdy_file = "wrfbdy_d01" diff --git a/Exec/RegTests/WPS_Test/inputs_real_ChisholmView b/Exec/RegTests/WPS_Test/inputs_real_ChisholmView index 097433d67..9f9fba31a 100644 --- a/Exec/RegTests/WPS_Test/inputs_real_ChisholmView +++ b/Exec/RegTests/WPS_Test/inputs_real_ChisholmView @@ -21,10 +21,6 @@ yhi.type = "Outflow" zlo.type = "NoSlipWall" zhi.type = "SlipWall" -# SET/RELAXATION REGIONS -erf.wrfbdy_width = 1 -erf.wrfbdy_set_width = 1 - # TIME STEP CONTROL erf.fixed_dt = 0.5 # fixed time step depending on grid resolution @@ -70,6 +66,8 @@ prob.rho_0 = 1.0 prob.T_0 = 300.0 # INITIALIZATION WITH ATM DATA -erf.init_type = "real" +erf.real_width = 1 +erf.real_set_width = 1 +erf.init_type = "real" erf.nc_init_file_0 = "wrfinput_chisholmview_d01" -erf.nc_bdy_file = "wrfbdy_chisholmview_d01_width1" +erf.nc_bdy_file = "wrfbdy_chisholmview_d01_width1" diff --git a/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp b/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp deleted file mode 100644 index 7c7be9a4c..000000000 --- a/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp +++ /dev/null @@ -1,193 +0,0 @@ -#include "ERF.H" -#include "Utils.H" - -using namespace amrex; - -#ifdef ERF_USE_NETCDF -/* - * Impose boundary conditions using data read in from met_em files - * - * @param[out] mfs Vector of MultiFabs to be filled - * @param[in] time time at which the data should be filled - */ - -void -ERF::fill_from_metgrid (const Vector& mfs, - const Real time, - bool cons_only, - int icomp_cons, - int ncomp_cons) -{ - int lev = 0; - - // Time interpolation - Real dT = bdy_time_interval; - Real time_since_start = time - start_bdy_time; - int n_time = static_cast( time_since_start / dT); - amrex::Real alpha = (time_since_start - n_time * dT) / dT; - AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0); - amrex::Real oma = 1.0 - alpha; - - // Flags for read vars and index mapping - // Cons includes [Rho RhoTheta RhoKE RhoQKE RhoScalar RhoQ1 RhoQ2 RhoQ3] - Vector cons_read = {1, 1, 0, 0, 0, 1, 0, 0, 0}; - Vector cons_map = {MetGridBdyVars::R, MetGridBdyVars::T, 0, 0, 0, MetGridBdyVars::QV, 0, 0, 0}; - - Vector> is_read; - is_read.push_back( cons_read ); - is_read.push_back( {1} ); // xvel - is_read.push_back( {1} ); // yvel - is_read.push_back( {0} ); // zvel - - Vector> ind_map; - ind_map.push_back( cons_map ); - ind_map.push_back( {MetGridBdyVars::U} ); // xvel - ind_map.push_back( {MetGridBdyVars::V} ); // yvel - ind_map.push_back( {0} ); // zvel - - // Nvars to loop over - Vector comp_var = {ncomp_cons, 1, 1, 1}; - - // End of vars loop - int var_idx_end = (cons_only) ? Vars::cons + 1 : Vars::NumTypes; - - // Loop over all variable types - for (int var_idx = Vars::cons; var_idx < var_idx_end; ++var_idx) - { - MultiFab& mf = *mfs[var_idx]; - - // - // Note that "domain" is mapped onto the type of box the data is in - // - Box domain = geom[lev].Domain(); - domain.convert(mf.boxArray().ixType()); - 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) - { - int width = metgrid_bdy_set_width; - - // Variable can be read from met_em files - //------------------------------------ - if (is_read[var_idx][comp_idx]) - { - int ivar = ind_map[var_idx][comp_idx]; - IntVect ng_vect = mf.nGrowVect(); ng_vect[2] = 0; - - // 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 - // as alpha * (data at time n+1) + (1 - alpha) * (data at time n) - const auto& bdatxlo_n = bdy_data_xlo[n_time ][ivar].const_array(); - const auto& bdatxlo_np1 = bdy_data_xlo[n_time+1][ivar].const_array(); - const auto& bdatxhi_n = bdy_data_xhi[n_time ][ivar].const_array(); - const auto& bdatxhi_np1 = bdy_data_xhi[n_time+1][ivar].const_array(); - const auto& bdatylo_n = bdy_data_ylo[n_time ][ivar].const_array(); - const auto& bdatylo_np1 = bdy_data_ylo[n_time+1][ivar].const_array(); - const auto& bdatyhi_n = bdy_data_yhi[n_time ][ivar].const_array(); - const auto& bdatyhi_np1 = bdy_data_yhi[n_time+1][ivar].const_array(); - -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(mf,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - // 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, - bx_xlo, bx_xhi, - bx_ylo, bx_yhi, ng_vect); - - // x-faces (includes exterior y ghost cells) - ParallelFor(bx_xlo, bx_xhi, - [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - int ii = std::max(i , dom_lo.x); - int jj = std::max(j , dom_lo.y); - jj = std::min(jj, dom_hi.y); - dest_arr(i,j,k,comp_idx) = oma * bdatxlo_n (ii,jj,k,0) - + alpha * bdatxlo_np1(ii,jj,k,0); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - int ii = std::min(i , dom_hi.x); - int jj = std::max(j , dom_lo.y); - jj = std::min(jj, dom_hi.y); - dest_arr(i,j,k,comp_idx) = oma * bdatxhi_n (ii,jj,k,0) - + alpha * bdatxhi_np1(ii,jj,k,0); - }); - - // y-faces (do not include exterior x ghost cells) - ParallelFor(bx_ylo, bx_yhi, - [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - int jj = std::max(j , dom_lo.y); - dest_arr(i,j,k,comp_idx) = oma * bdatylo_n (i,jj,k,0) - + alpha * bdatylo_np1(i,jj,k,0); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - int jj = std::min(j , dom_hi.y); - dest_arr(i,j,k,comp_idx) = oma * bdatyhi_n (i,jj,k,0) - + alpha * bdatyhi_np1(i,jj,k,0); - }); - } // mfi - - // Variable not read or computed from met_em files - //------------------------------------ - } else { - IntVect ng_vect = mf.nGrowVect(); ng_vect[2] = 0; - -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(mf,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - // Grown tilebox so we fill exterior ghost cells as well - 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, 0, - bx_xlo, bx_xhi, - bx_ylo, bx_yhi, ng_vect); - - // x-faces (includes y ghost cells) - 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); - 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); - dest_arr(i,j,k,comp_idx) = dest_arr(dom_hi.x-width,jj,k,comp_idx); - }); - - // y-faces (does not include x ghost cells) - ParallelFor(bx_ylo, bx_yhi, - [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - dest_arr(i,j,k,comp_idx) = dest_arr(i,dom_lo.y+width,k,comp_idx); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - dest_arr(i,j,k,comp_idx) = dest_arr(i,dom_hi.y-width,k,comp_idx); - }); - } // mfi - } // is_read - } // comp - } // var -} -#endif diff --git a/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp b/Source/BoundaryConditions/BoundaryConditions_realbdy.cpp similarity index 92% rename from Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp rename to Source/BoundaryConditions/BoundaryConditions_realbdy.cpp index 620cea160..970dacc35 100644 --- a/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_realbdy.cpp @@ -12,11 +12,11 @@ using namespace amrex; */ void -ERF::fill_from_wrfbdy (const Vector& mfs, - const Real time, - bool cons_only, - int icomp_cons, - int ncomp_cons) +ERF::fill_from_realbdy (const Vector& mfs, + const Real time, + bool cons_only, + int icomp_cons, + int ncomp_cons) { int lev = 0; @@ -30,7 +30,7 @@ ERF::fill_from_wrfbdy (const Vector& mfs, // Flags for read vars and index mapping Vector cons_read = {1, 1, 0, 0, 0, 1, 0, 0, 0}; - Vector cons_map = {WRFBdyVars::R, WRFBdyVars::T, 0, 0, 0, WRFBdyVars::QV, 0, 0, 0}; + Vector cons_map = {RealBdyVars::R, RealBdyVars::T, 0, 0, 0, RealBdyVars::QV, 0, 0, 0}; Vector> is_read; is_read.push_back( cons_read ); @@ -40,8 +40,8 @@ ERF::fill_from_wrfbdy (const Vector& mfs, Vector> ind_map; ind_map.push_back( cons_map ); - ind_map.push_back( {WRFBdyVars::U} ); // xvel - ind_map.push_back( {WRFBdyVars::V} ); // yvel + ind_map.push_back( {RealBdyVars::U} ); // xvel + ind_map.push_back( {RealBdyVars::V} ); // yvel ind_map.push_back( {0} ); // zvel // Nvars to loop over @@ -69,7 +69,7 @@ ERF::fill_from_wrfbdy (const Vector& mfs, // Loop over each component for (int comp_idx(offset); comp_idx < (comp_var[var_idx]+offset); ++comp_idx) { - int width = wrfbdy_set_width; + int width = real_set_width; // Variable can be read from wrf bdy //------------------------------------ @@ -163,13 +163,13 @@ ERF::fill_from_wrfbdy (const Vector& mfs, { 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); + 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+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); + dest_arr(i,j,k,comp_idx) = dest_arr(dom_hi.x-width,jj,k,comp_idx); }); // y-faces (does not include x ghost cells) diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index a7edfc689..e5907331b 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -91,9 +91,8 @@ ERF::FillPatch (int lev, Real time, const Vector& mfs, bool fillset) #ifdef ERF_USE_NETCDF // We call this here because it is an ERF routine - if (use_real_bcs) { - if (init_type == "real" && lev==0) fill_from_wrfbdy (mfs,time,false,0,ncomp_cons); - if (init_type == "metgrid" && lev==0) fill_from_metgrid(mfs,time,false,0,ncomp_cons); + if (use_real_bcs && (lev==0)) { + fill_from_realbdy(mfs,time,false,0,ncomp_cons); } #endif @@ -243,11 +242,9 @@ ERF::FillIntermediatePatch (int lev, Real time, IntVect ngvect_vels = IntVect(ng_vel ,ng_vel ,ng_vel); #ifdef ERF_USE_NETCDF - // NOTE: This routine needs to be aware of what FillIntermediatePatch is operating on - // --- i.e., cons_only and which cons indices (icomp_cons & ncomp_cons) - if (use_real_bcs) { - if (init_type == "real" && lev==0) fill_from_wrfbdy(mfs, time, cons_only, icomp_cons, ncomp_cons); - if (init_type == "metgrid" && lev==0) fill_from_metgrid(mfs, time, cons_only, icomp_cons, ncomp_cons); + // We call this here because it is an ERF routine + if (use_real_bcs && (lev==0)) { + fill_from_realbdy(mfs,time,false,0,ncomp_cons); } #endif diff --git a/Source/BoundaryConditions/Make.package b/Source/BoundaryConditions/Make.package index b9bb5be57..41886ad34 100644 --- a/Source/BoundaryConditions/Make.package +++ b/Source/BoundaryConditions/Make.package @@ -3,8 +3,7 @@ CEXE_sources += BoundaryConditions_yvel.cpp CEXE_sources += BoundaryConditions_zvel.cpp CEXE_sources += BoundaryConditions_cons.cpp CEXE_sources += BoundaryConditions_bndryreg.cpp -CEXE_sources += BoundaryConditions_wrfbdy.cpp -CEXE_sources += BoundaryConditions_metgrid.cpp +CEXE_sources += BoundaryConditions_realbdy.cpp CEXE_headers += ABLMost.H CEXE_sources += ABLMost.cpp diff --git a/Source/ERF.H b/Source/ERF.H index dee2b187b..4a9d771da 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -310,16 +310,11 @@ public: amrex::Real time); #ifdef ERF_USE_NETCDF - void fill_from_wrfbdy (const amrex::Vector& mfs, - amrex::Real time, - bool cons_only, - int icomp_cons, - int ncomp_cons); - void fill_from_metgrid (const amrex::Vector& mfs, - amrex::Real time, - bool cons_only, - int icomp_cons, - int ncomp_cons); + void fill_from_realbdy (const amrex::Vector& mfs, + amrex::Real time, + bool cons_only, + int icomp_cons, + int ncomp_cons); #endif #ifdef ERF_USE_NETCDF @@ -729,14 +724,10 @@ private: // NetCDF initialization (wrfinput) file static amrex::Vector> nc_init_file; - // NetCDF initialization (wrfbdy) file + // NetCDF initialization (wrfbdy/met_em) file static std::string nc_bdy_file; - int wrfbdy_width{0}; - int wrfbdy_set_width{0}; - - // NetCDF initialization (met_em) file - int metgrid_bdy_width{0}; - int metgrid_bdy_set_width{0}; + int real_width{0}; + int real_set_width{0}; // Text input_sounding file static std::string input_sounding_file; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 25b53aeca..435055a4e 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1131,18 +1131,11 @@ ERF::ReadParameters () pp.query("input_bndry_planes", input_bndry_planes); // Query the set and total widths for wrfbdy interior ghost cells - pp.query("wrfbdy_width", wrfbdy_width); - pp.query("wrfbdy_set_width", wrfbdy_set_width); - AMREX_ALWAYS_ASSERT(wrfbdy_width >= 0); - AMREX_ALWAYS_ASSERT(wrfbdy_set_width >= 0); - AMREX_ALWAYS_ASSERT(wrfbdy_width >= wrfbdy_set_width); - - // Query the set and total widths for metgrid_bdy interior ghost cells - pp.query("metgrid_bdy_width", metgrid_bdy_width); - pp.query("metgrid_bdy_set_width", metgrid_bdy_set_width); - AMREX_ALWAYS_ASSERT(metgrid_bdy_width >= 0); - AMREX_ALWAYS_ASSERT(metgrid_bdy_set_width >= 0); - AMREX_ALWAYS_ASSERT(metgrid_bdy_width >= metgrid_bdy_set_width); + pp.query("real_width", real_width); + pp.query("real_set_width", real_set_width); + AMREX_ALWAYS_ASSERT(real_width >= 0); + AMREX_ALWAYS_ASSERT(real_set_width >= 0); + AMREX_ALWAYS_ASSERT(real_width >= real_set_width); // Query the set and total widths for crse-fine interior ghost cells pp.query("cf_width", cf_width); diff --git a/Source/IO/Checkpoint.cpp b/Source/IO/Checkpoint.cpp index b37f7ce01..46689c94c 100644 --- a/Source/IO/Checkpoint.cpp +++ b/Source/IO/Checkpoint.cpp @@ -196,7 +196,7 @@ ERF::WriteCheckpointFile () const bdy_h_file << num_var << "\n"; bdy_h_file << start_bdy_time << "\n"; bdy_h_file << bdy_time_interval << "\n"; - bdy_h_file << wrfbdy_width << "\n"; + bdy_h_file << real_width << "\n"; for (int ivar(0); ivar> num_var; bdy_h_file >> start_bdy_time; bdy_h_file >> bdy_time_interval; - bdy_h_file >> wrfbdy_width; + bdy_h_file >> real_width; bx_v.resize(4*num_var); for (int ivar(0); ivar> bx_v[4*ivar ]; @@ -452,7 +452,7 @@ ERF::ReadCheckpointFile () ParallelDescriptor::Barrier(); ParallelDescriptor::Bcast(&start_bdy_time,1,ioproc); ParallelDescriptor::Bcast(&bdy_time_interval,1,ioproc); - ParallelDescriptor::Bcast(&wrfbdy_width,1,ioproc); + ParallelDescriptor::Bcast(&real_width,1,ioproc); ParallelDescriptor::Bcast(&num_time,1,ioproc); ParallelDescriptor::Bcast(&num_var,1,ioproc); diff --git a/Source/IO/NCWpsFile.H b/Source/IO/NCWpsFile.H index 5e4994d64..fd949ca9e 100644 --- a/Source/IO/NCWpsFile.H +++ b/Source/IO/NCWpsFile.H @@ -258,12 +258,12 @@ fill_fab_from_arrays (int iv, var_name == "MAPFAC_V" || var_name == "MAPFAC_VY") my_box.setType(amrex::IndexType(amrex::IntVect(0,1,0))); if (var_name == "W" || var_name == "WW") my_box.setType(amrex::IndexType(amrex::IntVect(0,0,1))); + amrex::Arena* Arena_Used = amrex::The_Arena(); #ifdef AMREX_USE_GPU // Make sure temp lives on CPU since nc_arrays lives on CPU only - temp.resize(my_box,1,amrex::The_Pinned_Arena()); -#else - temp.resize(my_box,1); + Arena_Used = amrex::The_Pinned_Arena(); #endif + temp.resize(my_box,1, Arena_Used); amrex::Array4 fab_arr = temp.array(); int ioff = temp.box().smallEnd()[0]; diff --git a/Source/IO/ReadFromMetgrid.cpp b/Source/IO/ReadFromMetgrid.cpp index b3dd3f1b4..5cea6c73d 100644 --- a/Source/IO/ReadFromMetgrid.cpp +++ b/Source/IO/ReadFromMetgrid.cpp @@ -29,22 +29,22 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname, std::vector attr; ncf.get_attr("FLAG_PSFC", attr); flag_psfc = attr[0]; - /* UNCOMMENT FOR FLOWMAS */ + /* UNCOMMENT FOR FLOWMAS flag_msfu = 0; //ncf.get_attr("FLAG_MAPFAC_U", attr); flag_msfu = attr[0]; flag_msfv = 0; //ncf.get_attr("FLAG_MAPFAC_V", attr); flag_msfv = attr[0]; flag_msfm = 0; //ncf.get_attr("FLAG_MAPFAC_M", attr); flag_msfm = attr[0]; flag_hgt = 1; //ncf.get_attr("FLAG_HGT_M", attr); flag_hgt = attr[0]; flag_sst = 0; //ncf.get_attr("FLAG_SST", attr); flag_sst = attr[0]; flag_lmask = 0; //ncf.get_attr("FLAG_LANDMASK", attr); flag_lmask = 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]; - */ + 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/IO/ReadFromWRFBdy.cpp b/Source/IO/ReadFromWRFBdy.cpp index 97de503ea..52c0f59ed 100644 --- a/Source/IO/ReadFromWRFBdy.cpp +++ b/Source/IO/ReadFromWRFBdy.cpp @@ -175,6 +175,11 @@ read_from_wrfbdy (const std::string& nc_bdy_file, const Box& domain, Box y_line_no_stag(IntVect(lo[0], 0, 0), IntVect(hi[0], 0, 0)); + Arena* Arena_Used = The_Arena(); +#ifdef AMREX_USE_GPU + Arena_Used = The_Pinned_Arena(); +#endif + if (bdyType == WRFBdyTypes::x_lo) { // ******************************************************************************* @@ -192,28 +197,28 @@ read_from_wrfbdy (const std::string& nc_bdy_file, const Box& domain, if (bdyVarType == WRFBdyVars::U) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_xlo[nt].push_back(FArrayBox(xlo_plane_x_stag, 1)); // U + bdy_data_xlo[nt].push_back(FArrayBox(xlo_plane_x_stag, 1, Arena_Used)); // U } } else if (bdyVarType == WRFBdyVars::V) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_xlo[nt].push_back(FArrayBox(xlo_plane_y_stag , 1)); // V + bdy_data_xlo[nt].push_back(FArrayBox(xlo_plane_y_stag , 1, Arena_Used)); // V } } else if (bdyVarType == WRFBdyVars::R) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_xlo[nt].push_back(FArrayBox(xlo_plane_no_stag, 1)); // R + bdy_data_xlo[nt].push_back(FArrayBox(xlo_plane_no_stag, 1, Arena_Used)); // R } } else if (bdyVarType == WRFBdyVars::T) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_xlo[nt].push_back(FArrayBox(xlo_plane_no_stag, 1)); // T + bdy_data_xlo[nt].push_back(FArrayBox(xlo_plane_no_stag, 1, Arena_Used)); // T } } else if (bdyVarType == WRFBdyVars::QV) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_xlo[nt].push_back(FArrayBox(xlo_plane_no_stag, 1)); // QV + bdy_data_xlo[nt].push_back(FArrayBox(xlo_plane_no_stag, 1, Arena_Used)); // QV } } else if (bdyVarType == WRFBdyVars::MU || bdyVarType == WRFBdyVars::PC) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_xlo[nt].push_back(FArrayBox(xlo_line, 1)); + bdy_data_xlo[nt].push_back(FArrayBox(xlo_line, 1, Arena_Used)); } } @@ -238,28 +243,28 @@ read_from_wrfbdy (const std::string& nc_bdy_file, const Box& domain, if (bdyVarType == WRFBdyVars::U) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_xhi[nt].push_back(FArrayBox(xhi_plane_x_stag, 1)); // U + bdy_data_xhi[nt].push_back(FArrayBox(xhi_plane_x_stag, 1, Arena_Used)); // U } } else if (bdyVarType == WRFBdyVars::V) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_xhi[nt].push_back(FArrayBox(xhi_plane_y_stag , 1)); // V + bdy_data_xhi[nt].push_back(FArrayBox(xhi_plane_y_stag , 1, Arena_Used)); // V } } else if (bdyVarType == WRFBdyVars::R) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_xhi[nt].push_back(FArrayBox(xhi_plane_no_stag, 1)); // R + bdy_data_xhi[nt].push_back(FArrayBox(xhi_plane_no_stag, 1, Arena_Used)); // R } } else if (bdyVarType == WRFBdyVars::T) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_xhi[nt].push_back(FArrayBox(xhi_plane_no_stag, 1)); // T + bdy_data_xhi[nt].push_back(FArrayBox(xhi_plane_no_stag, 1, Arena_Used)); // T } } else if (bdyVarType == WRFBdyVars::QV) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_xhi[nt].push_back(FArrayBox(xhi_plane_no_stag, 1)); // QV + bdy_data_xhi[nt].push_back(FArrayBox(xhi_plane_no_stag, 1, Arena_Used)); // QV } } else if (bdyVarType == WRFBdyVars::MU || bdyVarType == WRFBdyVars::PC) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_xhi[nt].push_back(FArrayBox(xhi_line, 1)); // MU + bdy_data_xhi[nt].push_back(FArrayBox(xhi_line, 1, Arena_Used)); // MU } } @@ -280,28 +285,28 @@ read_from_wrfbdy (const std::string& nc_bdy_file, const Box& domain, if (bdyVarType == WRFBdyVars::U) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_ylo[nt].push_back(FArrayBox(ylo_plane_x_stag , 1)); // U + bdy_data_ylo[nt].push_back(FArrayBox(ylo_plane_x_stag , 1, Arena_Used)); // U } } else if (bdyVarType == WRFBdyVars::V) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_ylo[nt].push_back(FArrayBox(ylo_plane_y_stag, 1)); // V + bdy_data_ylo[nt].push_back(FArrayBox(ylo_plane_y_stag, 1, Arena_Used)); // V } } else if (bdyVarType == WRFBdyVars::R) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_ylo[nt].push_back(FArrayBox(ylo_plane_no_stag, 1)); // R + bdy_data_ylo[nt].push_back(FArrayBox(ylo_plane_no_stag, 1, Arena_Used)); // R } } else if (bdyVarType == WRFBdyVars::T) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_ylo[nt].push_back(FArrayBox(ylo_plane_no_stag, 1)); // T + bdy_data_ylo[nt].push_back(FArrayBox(ylo_plane_no_stag, 1, Arena_Used)); // T } } else if (bdyVarType == WRFBdyVars::QV) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_ylo[nt].push_back(FArrayBox(ylo_plane_no_stag, 1)); // QV + bdy_data_ylo[nt].push_back(FArrayBox(ylo_plane_no_stag, 1, Arena_Used)); // QV } } else if (bdyVarType == WRFBdyVars::MU || bdyVarType == WRFBdyVars::PC) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_ylo[nt].push_back(FArrayBox(ylo_line, 1)); // PC + bdy_data_ylo[nt].push_back(FArrayBox(ylo_line, 1, Arena_Used)); // PC } } @@ -326,28 +331,28 @@ read_from_wrfbdy (const std::string& nc_bdy_file, const Box& domain, if (bdyVarType == WRFBdyVars::U) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_yhi[nt].push_back(FArrayBox(yhi_plane_x_stag , 1)); // U + bdy_data_yhi[nt].push_back(FArrayBox(yhi_plane_x_stag , 1, Arena_Used)); // U } } else if (bdyVarType == WRFBdyVars::V) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_yhi[nt].push_back(FArrayBox(yhi_plane_y_stag, 1)); // V + bdy_data_yhi[nt].push_back(FArrayBox(yhi_plane_y_stag, 1, Arena_Used)); // V } } else if (bdyVarType == WRFBdyVars::R) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_yhi[nt].push_back(FArrayBox(yhi_plane_no_stag, 1)); // R + bdy_data_yhi[nt].push_back(FArrayBox(yhi_plane_no_stag, 1, Arena_Used)); // R } } else if (bdyVarType == WRFBdyVars::T) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_yhi[nt].push_back(FArrayBox(yhi_plane_no_stag, 1)); // T + bdy_data_yhi[nt].push_back(FArrayBox(yhi_plane_no_stag, 1, Arena_Used)); // T } } else if (bdyVarType == WRFBdyVars::QV) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_yhi[nt].push_back(FArrayBox(yhi_plane_no_stag, 1)); // QV + bdy_data_yhi[nt].push_back(FArrayBox(yhi_plane_no_stag, 1, Arena_Used)); // QV } } else if (bdyVarType == WRFBdyVars::MU || bdyVarType == WRFBdyVars::PC) { for (int nt(0); nt < ntimes; ++nt) { - bdy_data_yhi[nt].push_back(FArrayBox(yhi_line, 1)); // PC + bdy_data_yhi[nt].push_back(FArrayBox(yhi_line, 1, Arena_Used)); // PC } } } @@ -530,7 +535,8 @@ convert_wrfbdy_data (int which, const Box& domain, Vector>& bd const FArrayBox& NC_C1H_fab, const FArrayBox& NC_C2H_fab, const FArrayBox& NC_RDNW_fab, const FArrayBox& NC_xvel_fab, const FArrayBox& NC_yvel_fab, - const FArrayBox& NC_rho_fab, const FArrayBox& NC_rhotheta_fab) + const FArrayBox& NC_rho_fab, const FArrayBox& NC_rhotheta_fab, + const FArrayBox& NC_QVAPOR_fab) { // These were filled from wrfinput Array4 c1h_arr = NC_C1H_fab.const_array(); @@ -556,147 +562,90 @@ convert_wrfbdy_data (int which, const Box& domain, Vector>& bd int jlo = domain.smallEnd()[1]; int jhi = domain.bigEnd()[1]; - const auto & bx_u = bdy_data[0][WRFBdyVars::U].box(); - ParallelFor(bx_u, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real xmu; - if (i == ilo) { - xmu = mu_arr(i,j,0) + mub_arr(i,j,0); - } else if (i > ihi) { - xmu = mu_arr(i-1,j,0) + mub_arr(i-1,j,0); - } else { - xmu = ( mu_arr(i,j,0) + mu_arr(i-1,j,0) - +mub_arr(i,j,0) + mub_arr(i-1,j,0)) * 0.5; - } - Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); - Real new_bdy = bdy_u_arr(i,j,k) / xmu_mult; - bdy_u_arr(i,j,k) = new_bdy; - }); - -#ifndef AMREX_USE_GPU - if (nt == 0) { - FArrayBox diff(bx_u,1); - diff.template copy(bdy_data[0][WRFBdyVars::U]); - diff.template minus(NC_xvel_fab); - if (which == 0) - amrex::Print() << "Max norm of diff between initial U and bdy U on lo x face: " << diff.norm(0) << std::endl; - if (which == 1) - amrex::Print() << "Max norm of diff between initial U and bdy U on hi x face: " << diff.norm(0) << std::endl; - if (which == 2) - amrex::Print() << "Max norm of diff between initial U and bdy U on lo y face: " << diff.norm(0) << std::endl; - if (which == 3) - amrex::Print() << "Max norm of diff between initial U and bdy U on hi y face: " << diff.norm(0) << std::endl; - } -#endif - - const auto & bx_v = bdy_data[0][WRFBdyVars::V].box(); - ParallelFor(bx_v, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real xmu; - if (j == jlo) { - xmu = mu_arr(i,j,0) + mub_arr(i,j,0); - } else if (j > jhi) { - xmu = mu_arr(i,j-1,0) + mub_arr(i,j-1,0); - } else { - xmu = ( mu_arr(i,j,0) + mu_arr(i,j-1,0) - +mub_arr(i,j,0) + mub_arr(i,j-1,0) ) * 0.5; - } - Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); - Real new_bdy = bdy_v_arr(i,j,k) / xmu_mult; - bdy_v_arr(i,j,k) = new_bdy; - }); - -#ifndef AMREX_USE_GPU - if (nt == 0) { - FArrayBox diff(bx_v,1); - diff.template copy(bdy_data[0][WRFBdyVars::V]); - diff.template minus(NC_yvel_fab); - if (which == 0) - amrex::Print() << "Max norm of diff between initial V and bdy V on lo x face: " << diff.norm(0) << std::endl; - if (which == 1) - amrex::Print() << "Max norm of diff between initial V and bdy V on hi x face: " << diff.norm(0) << std::endl; - if (which == 2) - amrex::Print() << "Max norm of diff between initial V and bdy V on lo y face: " << diff.norm(0) << std::endl; - if (which == 3) - amrex::Print() << "Max norm of diff between initial V and bdy V on hi y face: " << diff.norm(0) << std::endl; - } -#endif - - const auto & bx_t = bdy_data[0][WRFBdyVars::T].box(); // Note this is currently "THM" aka the perturbational moist pot. temp. - - // Define density - ParallelFor(bx_t, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - - Real xmu = c1h_arr(0,0,k) * (mu_arr(i,j,0) + mub_arr(i,j,0)) + c2h_arr(0,0,k); - - Real dpht = (ph_arr(i,j,k+1) + phb_arr(i,j,k+1)) - (ph_arr(i,j,k) + phb_arr(i,j,k)); - - bdy_r_arr(i,j,k) = -xmu / ( dpht * rdnw_arr(0,0,k) ); - - //if (nt == 0 and std::abs(r_arr(i,j,k) - bdy_r_arr(i,j,k)) > 0.) { - // amrex::Print() << "INIT VS BDY DEN " << IntVect(i,j,k) << " " << r_arr(i,j,k) << " " << bdy_r_arr(i,j,k) << - // " " << std::abs(r_arr(i,j,k) - bdy_r_arr(i,j,k)) << std::endl; - //} - }); - - // Define theta - amrex::Real theta_ref = 300.; - ParallelFor(bx_t, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - - Real xmu = (mu_arr(i,j,0) + mub_arr(i,j,0)); - Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); - - Real new_bdy_Th = bdy_t_arr(i,j,k) / xmu_mult + theta_ref; - - Real qv_fac = (1. + bdy_qv_arr(i,j,k) / 0.622 / xmu_mult); - - new_bdy_Th /= qv_fac; - - bdy_t_arr(i,j,k) = new_bdy_Th * bdy_r_arr(i,j,k); - - //if (nt == 0 and std::abs(rth_arr(i,j,k) - bdy_t_arr(i,j,k)) > 0.) { - // amrex::Print() << "INIT VS BDY TH " << IntVect(i,j,k) << " " << rth_arr(i,j,k) << " " << bdy_t_arr(i,j,k) << - // " " << std::abs(th_arr(i,j,k) - bdy_t_arr(i,j,k)) << std::endl; - //} - }); - -#ifndef AMREX_USE_GPU - if (nt == 0) { - FArrayBox diff(bx_t,1); - diff.template copy(bdy_data[0][WRFBdyVars::R]); - //diff.template mult(NC_rho_fab); - diff.template minus(NC_rho_fab); - if (which == 0) - amrex::Print() << "Max norm of diff between initial r and bdy r on lo x face: " << diff.norm(0) << std::endl; - if (which == 1) - amrex::Print() << "Max norm of diff between initial r and bdy r on hi x face: " << diff.norm(0) << std::endl; - if (which == 2) - amrex::Print() << "Max norm of diff between initial r and bdy r on lo y face: " << diff.norm(0) << std::endl; - if (which == 3) - amrex::Print() << "Max norm of diff between initial r and bdy r on hi y face: " << diff.norm(0) << std::endl; - - diff.template copy(bdy_data[0][WRFBdyVars::T]); - diff.template minus(NC_rhotheta_fab); - if (which == 0) - amrex::Print() << "Max norm of diff between initial rTh and bdy rTh on lo x face: " << diff.norm(0) << std::endl; - if (which == 1) - amrex::Print() << "Max norm of diff between initial rTh and bdy rTh on hi x face: " << diff.norm(0) << std::endl; - if (which == 2) - amrex::Print() << "Max norm of diff between initial rTh and bdy rTh on lo y face: " << diff.norm(0) << std::endl; - if (which == 3) - amrex::Print() << "Max norm of diff between initial rTh and bdy rTh on hi y face: " << diff.norm(0) << std::endl; - } -#endif - - // Define Qv - const auto & bx_qv = bdy_data[0][WRFBdyVars::QV].box(); - ParallelFor(bx_qv, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - - Real xmu = (mu_arr(i,j,0) + mub_arr(i,j,0)); - Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); - - Real new_bdy_QV = bdy_qv_arr(i,j,k) / xmu_mult; + if (nt==0) { + bdy_data[0][WRFBdyVars::U].template copy(NC_xvel_fab); + bdy_data[0][WRFBdyVars::V].template copy(NC_yvel_fab); + bdy_data[0][WRFBdyVars::R].template copy(NC_rho_fab); + bdy_data[0][WRFBdyVars::T].template copy(NC_rhotheta_fab); + bdy_data[0][WRFBdyVars::QV].template copy(NC_QVAPOR_fab); + bdy_data[0][WRFBdyVars::QV].template mult(NC_rho_fab); + } else { + // Define u velocity + const auto & bx_u = bdy_data[0][WRFBdyVars::U].box(); + ParallelFor(bx_u, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real xmu; + if (i == ilo) { + xmu = mu_arr(i,j,0) + mub_arr(i,j,0); + } else if (i > ihi) { + xmu = mu_arr(i-1,j,0) + mub_arr(i-1,j,0); + } else { + xmu = ( mu_arr(i,j,0) + mu_arr(i-1,j,0) + +mub_arr(i,j,0) + mub_arr(i-1,j,0)) * 0.5; + } + Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); + Real new_bdy = bdy_u_arr(i,j,k) / xmu_mult; + bdy_u_arr(i,j,k) = new_bdy; + }); + + // Define v velocity + const auto & bx_v = bdy_data[0][WRFBdyVars::V].box(); + ParallelFor(bx_v, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real xmu; + if (j == jlo) { + xmu = mu_arr(i,j,0) + mub_arr(i,j,0); + } else if (j > jhi) { + xmu = mu_arr(i,j-1,0) + mub_arr(i,j-1,0); + } else { + xmu = ( mu_arr(i,j,0) + mu_arr(i,j-1,0) + +mub_arr(i,j,0) + mub_arr(i,j-1,0) ) * 0.5; + } + Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); + Real new_bdy = bdy_v_arr(i,j,k) / xmu_mult; + bdy_v_arr(i,j,k) = new_bdy; + }); + + // Define density + const auto & bx_t = bdy_data[0][WRFBdyVars::T].box(); // Note this is currently "THM" aka the perturbational moist pot. temp. + ParallelFor(bx_t, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real xmu = c1h_arr(0,0,k) * (mu_arr(i,j,0) + mub_arr(i,j,0)) + c2h_arr(0,0,k); + Real dpht = (ph_arr(i,j,k+1) + phb_arr(i,j,k+1)) - (ph_arr(i,j,k) + phb_arr(i,j,k)); + bdy_r_arr(i,j,k) = -xmu / ( dpht * rdnw_arr(0,0,k) ); + //if (nt == 0 and std::abs(r_arr(i,j,k) - bdy_r_arr(i,j,k)) > 0.) { + // amrex::Print() << "INIT VS BDY DEN " << IntVect(i,j,k) << " " << r_arr(i,j,k) << " " << bdy_r_arr(i,j,k) << + // " " << std::abs(r_arr(i,j,k) - bdy_r_arr(i,j,k)) << std::endl; + //} + }); + + // Define theta + amrex::Real theta_ref = 300.; + ParallelFor(bx_t, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real xmu = (mu_arr(i,j,0) + mub_arr(i,j,0)); + Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); + Real new_bdy_Th = bdy_t_arr(i,j,k) / xmu_mult + theta_ref; + Real qv_fac = (1. + bdy_qv_arr(i,j,k) / 0.622 / xmu_mult); + new_bdy_Th /= qv_fac; + bdy_t_arr(i,j,k) = new_bdy_Th * bdy_r_arr(i,j,k); + //if (nt == 0 and std::abs(rth_arr(i,j,k) - bdy_t_arr(i,j,k)) > 0.) { + // amrex::Print() << "INIT VS BDY TH " << IntVect(i,j,k) << " " << rth_arr(i,j,k) << " " << bdy_t_arr(i,j,k) << + // " " << std::abs(th_arr(i,j,k) - bdy_t_arr(i,j,k)) << std::endl; + //} + }); + + // Define Qv + const auto & bx_qv = bdy_data[0][WRFBdyVars::QV].box(); + ParallelFor(bx_qv, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real xmu = (mu_arr(i,j,0) + mub_arr(i,j,0)); + Real xmu_mult = c1h_arr(0,0,k) * xmu + c2h_arr(0,0,k); + Real new_bdy_QV = bdy_qv_arr(i,j,k) / xmu_mult; + bdy_qv_arr(i,j,k) = new_bdy_QV * bdy_r_arr(i,j,k); + }); - bdy_qv_arr(i,j,k) = new_bdy_QV * bdy_r_arr(i,j,k); - }); + } // nt ==0 } // ntimes } #endif // ERF_USE_NETCDF diff --git a/Source/IO/ReadFromWRFInput.cpp b/Source/IO/ReadFromWRFInput.cpp index 50724fb46..f1a6aa63a 100644 --- a/Source/IO/ReadFromWRFInput.cpp +++ b/Source/IO/ReadFromWRFInput.cpp @@ -20,8 +20,11 @@ read_from_wrfinput (int lev, FArrayBox& NC_QVAPOR_fab, FArrayBox& NC_QCLOUD_fab, FArrayBox& NC_QRAIN_fab, - FArrayBox& NC_PH_fab , FArrayBox& NC_PHB_fab, - FArrayBox& NC_ALB_fab , FArrayBox& NC_PB_fab, + FArrayBox& NC_PH_fab, + FArrayBox& NC_P_fab, + FArrayBox& NC_PHB_fab, + FArrayBox& NC_ALB_fab, + FArrayBox& NC_PB_fab, MoistureType moisture_type) { amrex::Print() << "Loading initial data from NetCDF file at level " << lev << std::endl; @@ -39,20 +42,21 @@ read_from_wrfinput (int lev, NC_fabs.push_back(&NC_PH_fab); NC_names.push_back("PH"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 6 NC_fabs.push_back(&NC_PHB_fab); NC_names.push_back("PHB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 7 NC_fabs.push_back(&NC_PB_fab); NC_names.push_back("PB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 9 - NC_fabs.push_back(&NC_ALB_fab); NC_names.push_back("ALB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 0 - NC_fabs.push_back(&NC_MUB_fab); NC_names.push_back("MUB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 10 - NC_fabs.push_back(&NC_MSFU_fab); NC_names.push_back("MAPFAC_UY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 11 - NC_fabs.push_back(&NC_MSFV_fab); NC_names.push_back("MAPFAC_VY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 12 - NC_fabs.push_back(&NC_MSFM_fab); NC_names.push_back("MAPFAC_MY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 13 - NC_fabs.push_back(&NC_SST_fab); NC_names.push_back("SST"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 14 - NC_fabs.push_back(&NC_C1H_fab); NC_names.push_back("C1H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 15 - NC_fabs.push_back(&NC_C2H_fab); NC_names.push_back("C2H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 16 - NC_fabs.push_back(&NC_RDNW_fab); NC_names.push_back("RDNW"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 17 + NC_fabs.push_back(&NC_P_fab); NC_names.push_back("P"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 10 + NC_fabs.push_back(&NC_ALB_fab); NC_names.push_back("ALB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 11 + NC_fabs.push_back(&NC_MUB_fab); NC_names.push_back("MUB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 12 + NC_fabs.push_back(&NC_MSFU_fab); NC_names.push_back("MAPFAC_UY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 13 + NC_fabs.push_back(&NC_MSFV_fab); NC_names.push_back("MAPFAC_VY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 14 + NC_fabs.push_back(&NC_MSFM_fab); NC_names.push_back("MAPFAC_MY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 15 + NC_fabs.push_back(&NC_SST_fab); NC_names.push_back("SST"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 16 + NC_fabs.push_back(&NC_C1H_fab); NC_names.push_back("C1H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 17 + NC_fabs.push_back(&NC_C2H_fab); NC_names.push_back("C2H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 18 + NC_fabs.push_back(&NC_RDNW_fab); NC_names.push_back("RDNW"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 19 if (moisture_type != MoistureType::None) { - NC_fabs.push_back(&NC_QVAPOR_fab); NC_names.push_back("QVAPOR"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 18 - NC_fabs.push_back(&NC_QCLOUD_fab); NC_names.push_back("QCLOUD"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 19 - NC_fabs.push_back(&NC_QRAIN_fab); NC_names.push_back("QRAIN"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 20 + NC_fabs.push_back(&NC_QVAPOR_fab); NC_names.push_back("QVAPOR"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 20 + NC_fabs.push_back(&NC_QCLOUD_fab); NC_names.push_back("QCLOUD"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 21 + NC_fabs.push_back(&NC_QRAIN_fab); NC_names.push_back("QRAIN"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 22 } // Read the netcdf file and fill these FABs diff --git a/Source/IndexDefines.H b/Source/IndexDefines.H index 10c894d59..b16d8a695 100644 --- a/Source/IndexDefines.H +++ b/Source/IndexDefines.H @@ -53,6 +53,17 @@ namespace BCVars { }; } +namespace RealBdyVars { + enum { + U = 0, + V = 1, + R = 2, + T = 3, + QV, + NumTypes + }; +} + namespace WRFBdyVars { enum { U = 0, diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index c40a3dee0..8134cbc19 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -70,10 +70,14 @@ ERF::init_from_metgrid (int lev) Vector NC_dx; NC_dx.resize( ntimes); Vector NC_dy; NC_dy.resize( ntimes); - 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; + // Define the arena to be used for data allocation + Arena* Arena_Used = The_Arena(); +#ifdef AMREX_USE_GPU + // Make sure this lives on CPU and GPU + Arena_Used = The_Pinned_Arena(); #endif + + for (int it = 0; it < ntimes; it++) { read_from_metgrid(lev, boxes_at_level[lev][0], nc_init_file[lev][it], NC_dateTime[it], NC_epochTime[it], flag_psfc[it], flag_msfu[it], flag_msfv[it], flag_msfm[it], @@ -84,20 +88,6 @@ ERF::init_from_metgrid (int lev) NC_ght_fab[it], NC_hgt_fab[it], NC_psfc_fab[it], 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; - // NC_dateTime is only on the IOProcessor. - amrex::Print() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << ": NC_dateTime \t" << NC_dateTime[it] << std::endl; -#endif } // it // Verify that files in nc_init_file[lev] are ordered from earliest to latest. @@ -114,9 +104,6 @@ ERF::init_from_metgrid (int lev) // Verify that met_em files have even spacing in time. 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; -#endif if (NC_dt != bdy_time_interval) amrex::Error("Time interval between consecutive met_em files must be consistent."); } @@ -127,13 +114,8 @@ ERF::init_from_metgrid (int lev) Vector theta_fab; theta_fab.resize(ntimes); for (int it = 0; it < ntimes; it++) { Box NC_box_unstag = NC_rhum_fab[it].box(); -#ifdef AMREX_USE_GPU - mxrat_fab[it].resize(NC_box_unstag, 1, The_Pinned_Arena()); - theta_fab[it].resize(NC_box_unstag, 1, The_Pinned_Arena()); -#else - mxrat_fab[it].resize(NC_box_unstag, 1); - theta_fab[it].resize(NC_box_unstag, 1); -#endif + mxrat_fab[it].resize(NC_box_unstag, 1, Arena_Used); + theta_fab[it].resize(NC_box_unstag, 1, Arena_Used); } auto& lev_new = vars_new[lev]; @@ -254,11 +236,7 @@ ERF::init_from_metgrid (int lev) auto ng = lev_new[Vars::cons].nGrowVect(); gdomain = amrex::grow(ldomain,ng); } -#ifdef AMREX_USE_GPU - fabs_for_bcs[it][nvar].resize(gdomain, 1, The_Pinned_Arena()); -#else - fabs_for_bcs[it][nvar].resize(gdomain, 1); -#endif + fabs_for_bcs[it][nvar].resize(gdomain, 1, Arena_Used); fabs_for_bcs[it][nvar].template setVal(0.0); } } @@ -366,9 +344,9 @@ ERF::init_from_metgrid (int lev) // NOTE: We must guarantee one halo cell in the bdy file. // Otherwise, we make the total width match the set width. - if (metgrid_bdy_width-1 <= metgrid_bdy_set_width) metgrid_bdy_width = metgrid_bdy_set_width; - amrex::Print() << "Running with specification width: " << metgrid_bdy_set_width - << " and relaxation width: " << metgrid_bdy_width - metgrid_bdy_set_width << std::endl; + if (real_width-1 <= real_set_width) real_width = real_set_width; + amrex::Print() << "Running with specification width: " << real_set_width + << " and relaxation width: " << real_width - real_set_width << std::endl; // Set up boxes for lateral boundary arrays. bdy_data_xlo.resize(ntimes); @@ -382,13 +360,13 @@ ERF::init_from_metgrid (int lev) amrex::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]; + phi[0] = lo[0]+real_width-1; phi[1] = hi[1]; phi[2] = hi[2]; const Box pbx_xlo(plo, phi); Box xlo_plane_no_stag(pbx_xlo); Box xlo_plane_x_stag = pbx_xlo; xlo_plane_x_stag.shiftHalf(0,-1); Box xlo_plane_y_stag = convert(pbx_xlo, {0, 1, 0}); - plo[0] = hi[0]-metgrid_bdy_width+1; plo[1] = lo[1]; plo[2] = lo[2]; + plo[0] = hi[0]-real_width+1; plo[1] = lo[1]; plo[2] = lo[2]; phi[0] = hi[0]; phi[1] = hi[1]; phi[2] = hi[2]; const Box pbx_xhi(plo, phi); Box xhi_plane_no_stag(pbx_xhi); @@ -396,42 +374,19 @@ ERF::init_from_metgrid (int lev) Box xhi_plane_y_stag = convert(pbx_xhi, {0, 1, 0}); plo[1] = lo[1]; plo[0] = lo[0]; plo[2] = lo[2]; - phi[1] = lo[1]+metgrid_bdy_width-1; phi[0] = hi[0]; phi[2] = hi[2]; + phi[1] = lo[1]+real_width-1; phi[0] = hi[0]; phi[2] = hi[2]; const Box pbx_ylo(plo, phi); Box ylo_plane_no_stag(pbx_ylo); Box ylo_plane_x_stag = convert(pbx_ylo, {1, 0, 0}); Box ylo_plane_y_stag = pbx_ylo; ylo_plane_y_stag.shiftHalf(1,-1); - plo[1] = hi[1]-metgrid_bdy_width+1; plo[0] = lo[0]; plo[2] = lo[2]; + plo[1] = hi[1]-real_width+1; plo[0] = lo[0]; plo[2] = lo[2]; phi[1] = hi[1]; phi[0] = hi[0]; phi[2] = hi[2]; const Box pbx_yhi(plo, phi); Box yhi_plane_no_stag(pbx_yhi); Box yhi_plane_x_stag = convert(pbx_yhi, {1, 0, 0}); 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; -#endif - for (int ivar(MetGridBdyVars::U); ivar < MetGridBdyEnd; ivar++) { for (int it(0); it < ntimes; it++) { if (ivar == MetGridBdyVars::U) { @@ -550,11 +505,6 @@ init_terrain_from_metgrid (FArrayBox& z_phys_nd_fab, int ntimes = 1; // Use terrain from the first met_em file. 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; -#endif - // This copies from NC_zphys on z-faces to z_phys_nd on nodes const Array4& z_arr = z_phys_nd_fab.array(); const Array4& nc_hgt_arr = NC_hgt_fab[it].const_array(); @@ -574,11 +524,6 @@ init_terrain_from_metgrid (FArrayBox& z_phys_nd_fab, Box bxu = bx; bxu.growHi(0,1); 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; -#endif - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { int ii = std::max(std::min(i,ihi-1),ilo+1); @@ -863,6 +808,13 @@ init_base_state_from_metgrid (const bool use_moisture, Gpu::copy(Gpu::hostToDevice, flag_psfc.begin(), flag_psfc.end(), flag_psfc_d.begin()); int* flag_psfc_vec = flag_psfc_d.data(); + // Define the arena to be used for data allocation + Arena* Arena_Used = The_Arena(); +#ifdef AMREX_USE_GPU + // Make sure this lives on CPU and GPU + Arena_Used = The_Pinned_Arena(); +#endif + { // set pressure and density at initialization. const Array4& r_hse_arr = r_hse_fab.array(); const Array4& p_hse_arr = p_hse_fab.array(); @@ -965,11 +917,7 @@ init_base_state_from_metgrid (const bool use_moisture, for (int it=0; it& NC_QVAPOR_fab, const Vector& NC_QCLOUD_fab, @@ -80,11 +86,19 @@ init_terrain_from_wrfinput (int lev, const Real& z_top, const Vector& NC_PHB_fab); void -init_base_state_from_wrfinput (int lev, const Box& bx, Real l_rdOcp, - FArrayBox& p_hse, FArrayBox& pi_hse, +init_base_state_from_wrfinput (int lev, + const Box& gtbx, + const Box& domain, + Real l_rdOcp, + MoistureType moisture_type, + const int& n_qstate, + FArrayBox& cons_fab, + FArrayBox& p_hse, + FArrayBox& pi_hse, FArrayBox& r_hse, const Vector& NC_ALB_fab, - const Vector& NC_PB_fab); + const Vector& NC_PB_fab, + const Vector& NC_P_fab); /** * ERF function that initializes data from a WRF dataset @@ -110,6 +124,7 @@ ERF::init_from_wrfinput (int lev) Vector NC_C2H_fab ; NC_C2H_fab.resize(num_boxes_at_level[lev]); Vector NC_RDNW_fab ; NC_RDNW_fab.resize(num_boxes_at_level[lev]); Vector NC_PH_fab ; NC_PH_fab.resize(num_boxes_at_level[lev]); + Vector NC_P_fab ; NC_P_fab.resize(num_boxes_at_level[lev]); Vector NC_PHB_fab ; NC_PHB_fab.resize(num_boxes_at_level[lev]); Vector NC_ALB_fab ; NC_ALB_fab.resize(num_boxes_at_level[lev]); Vector NC_PB_fab ; NC_PB_fab.resize(num_boxes_at_level[lev]); @@ -124,12 +139,13 @@ ERF::init_from_wrfinput (int lev) for (int idx = 0; idx < num_boxes_at_level[lev]; idx++) { read_from_wrfinput(lev, boxes_at_level[lev][idx], nc_init_file[lev][idx], - NC_xvel_fab[idx], NC_yvel_fab[idx], NC_zvel_fab[idx], NC_rho_fab[idx], - NC_rhop_fab[idx], NC_rhoth_fab[idx], NC_MUB_fab[idx], - NC_MSFU_fab[idx], NC_MSFV_fab[idx], NC_MSFM_fab[idx], - NC_SST_fab[idx], NC_C1H_fab[idx], NC_C2H_fab[idx], NC_RDNW_fab[idx], + NC_xvel_fab[idx] , NC_yvel_fab[idx] , NC_zvel_fab[idx] , NC_rho_fab[idx], + NC_rhop_fab[idx] , NC_rhoth_fab[idx] , NC_MUB_fab[idx] , + NC_MSFU_fab[idx] , NC_MSFV_fab[idx] , NC_MSFM_fab[idx] , + NC_SST_fab[idx] , NC_C1H_fab[idx] , NC_C2H_fab[idx] , NC_RDNW_fab[idx], NC_QVAPOR_fab[idx], NC_QCLOUD_fab[idx], NC_QRAIN_fab[idx], - NC_PH_fab[idx],NC_PHB_fab[idx],NC_ALB_fab[idx],NC_PB_fab[idx], + NC_PH_fab[idx] , NC_P_fab[idx] , NC_PHB_fab[idx] , + NC_ALB_fab[idx] , NC_PB_fab[idx] , solverChoice.moisture_type); } @@ -171,7 +187,7 @@ ERF::init_from_wrfinput (int lev) } // mf const Box& domain = geom[lev].Domain(); - const Real& z_top = geom[lev].ProbHi(2); + const Real& z_top = geom[lev].ProbHi(2); if (solverChoice.use_terrain) { std::unique_ptr& z_phys = z_phys_nd[lev]; @@ -190,20 +206,27 @@ ERF::init_from_wrfinput (int lev) MultiFab p_hse (base_state[lev], make_alias, 1, 1); // p_0 is second component MultiFab pi_hse(base_state[lev], make_alias, 2, 1); // pi_0 is third component + IntVect ng = p_hse.nGrowVect(); const Real l_rdOcp = solverChoice.rdOcp; if (init_type == "real") { for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + FArrayBox& cons_fab = lev_new[Vars::cons][mfi]; FArrayBox& p_hse_fab = p_hse[mfi]; FArrayBox& pi_hse_fab = pi_hse[mfi]; FArrayBox& r_hse_fab = r_hse[mfi]; - const Box valid_bx = mfi.validbox(); - init_base_state_from_wrfinput(lev, valid_bx, l_rdOcp, - p_hse_fab, pi_hse_fab, r_hse_fab, - NC_ALB_fab, NC_PB_fab); + const Box gtbx = mfi.tilebox(IntVect(0), ng); + init_base_state_from_wrfinput(lev, gtbx, domain, l_rdOcp, solverChoice.moisture_type, n_qstate, + cons_fab, p_hse_fab, pi_hse_fab, r_hse_fab, + NC_ALB_fab, NC_PB_fab, NC_P_fab); } + + // FillBoundary to populate the internal halo cells + r_hse.FillBoundary(geom[lev].periodicity()); + p_hse.FillBoundary(geom[lev].periodicity()); + pi_hse.FillBoundary(geom[lev].periodicity()); } if (init_type == "real" && (lev == 0)) { @@ -211,28 +234,28 @@ ERF::init_from_wrfinput (int lev) amrex::Error("NetCDF boundary file name must be provided via input"); bdy_time_interval = read_from_wrfbdy(nc_bdy_file,geom[0].Domain(), bdy_data_xlo,bdy_data_xhi,bdy_data_ylo,bdy_data_yhi, - wrfbdy_width, start_bdy_time); + real_width, start_bdy_time); - amrex::Print() << "Read in boundary data with width " << wrfbdy_width << std::endl; - amrex::Print() << "Running with specification width: " << wrfbdy_set_width - << " and relaxation width: " << wrfbdy_width - wrfbdy_set_width << std::endl; + amrex::Print() << "Read in boundary data with width " << real_width << std::endl; + amrex::Print() << "Running with specification width: " << real_set_width + << " and relaxation width: " << real_width - real_set_width << std::endl; convert_wrfbdy_data(0,domain,bdy_data_xlo, - NC_MUB_fab[0], NC_PH_fab[0] , NC_PHB_fab[0], - NC_C1H_fab[0], NC_C2H_fab[0], NC_RDNW_fab[0], - NC_xvel_fab[0],NC_yvel_fab[0],NC_rho_fab[0],NC_rhoth_fab[0]); + NC_MUB_fab[0] , NC_PH_fab[0] , NC_PHB_fab[0] , + NC_C1H_fab[0] , NC_C2H_fab[0] , NC_RDNW_fab[0], + NC_xvel_fab[0], NC_yvel_fab[0], NC_rho_fab[0] , NC_rhoth_fab[0], NC_QVAPOR_fab[0]); convert_wrfbdy_data(1,domain,bdy_data_xhi, - NC_MUB_fab[0], NC_PH_fab[0] , NC_PHB_fab[0], - NC_C1H_fab[0], NC_C2H_fab[0], NC_RDNW_fab[0], - NC_xvel_fab[0],NC_yvel_fab[0],NC_rho_fab[0],NC_rhoth_fab[0]); + NC_MUB_fab[0] , NC_PH_fab[0] , NC_PHB_fab[0] , + NC_C1H_fab[0] , NC_C2H_fab[0] , NC_RDNW_fab[0], + NC_xvel_fab[0], NC_yvel_fab[0], NC_rho_fab[0] , NC_rhoth_fab[0], NC_QVAPOR_fab[0]); convert_wrfbdy_data(2,domain,bdy_data_ylo, - NC_MUB_fab[0], NC_PH_fab[0] , NC_PHB_fab[0], - NC_C1H_fab[0], NC_C2H_fab[0], NC_RDNW_fab[0], - NC_xvel_fab[0],NC_yvel_fab[0],NC_rho_fab[0],NC_rhoth_fab[0]); + NC_MUB_fab[0] , NC_PH_fab[0] , NC_PHB_fab[0] , + NC_C1H_fab[0] , NC_C2H_fab[0] , NC_RDNW_fab[0], + NC_xvel_fab[0], NC_yvel_fab[0], NC_rho_fab[0] , NC_rhoth_fab[0], NC_QVAPOR_fab[0]); convert_wrfbdy_data(3,domain,bdy_data_yhi, - NC_MUB_fab[0], NC_PH_fab[0] , NC_PHB_fab[0], - NC_C1H_fab[0], NC_C2H_fab[0], NC_RDNW_fab[0], - NC_xvel_fab[0],NC_yvel_fab[0],NC_rho_fab[0],NC_rhoth_fab[0]); + NC_MUB_fab[0] , NC_PH_fab[0] , NC_PHB_fab[0] , + NC_C1H_fab[0] , NC_C2H_fab[0] , NC_RDNW_fab[0], + NC_xvel_fab[0], NC_yvel_fab[0], NC_rho_fab[0] , NC_rhoth_fab[0], NC_QVAPOR_fab[0]); } // Start at the earliest time (read_from_wrfbdy) @@ -360,11 +383,23 @@ init_msfs_from_wrfinput (int lev, FArrayBox& msfu_fab, * @param NC_PB_fab Vector of FArrayBox objects containing WRF data specifying pressure */ void -init_base_state_from_wrfinput (int lev, const Box& valid_bx, const Real l_rdOcp, - FArrayBox& p_hse, FArrayBox& pi_hse, FArrayBox& r_hse, +init_base_state_from_wrfinput (int lev, + const Box& gtbx, + const Box& domain, + const Real l_rdOcp, + MoistureType moisture_type, + const int& n_qstate, + FArrayBox& cons_fab, + FArrayBox& p_hse, + FArrayBox& pi_hse, + FArrayBox& r_hse, const Vector& NC_ALB_fab, - const Vector& NC_PB_fab) + const Vector& NC_PB_fab, + const Vector& NC_P_fab) { + const auto& dom_lo = lbound(domain); + const auto& dom_hi = ubound(domain); + int nboxes = NC_ALB_fab.size(); for (int idx = 0; idx < nboxes; idx++) { @@ -372,16 +407,42 @@ init_base_state_from_wrfinput (int lev, const Box& valid_bx, const Real l_rdOcp, // FArrayBox to FArrayBox copy does "copy on intersection" // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks // - const Array4& p_hse_arr = p_hse.array(); - const Array4& pi_hse_arr = pi_hse.array(); const Array4& r_hse_arr = r_hse.array(); - const Array4& alpha_arr = NC_ALB_fab[idx].const_array(); - const Array4& nc_pb_arr = NC_PB_fab[idx].const_array(); - - ParallelFor(valid_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - p_hse_arr(i,j,k) = nc_pb_arr(i,j,k); + const Array4& cons_arr = cons_fab.array(); + const Array4& p_hse_arr = p_hse.array(); + const Array4& pi_hse_arr = pi_hse.array(); + const Array4& r_hse_arr = r_hse.array(); + const Array4& alpha_arr = NC_ALB_fab[idx].const_array(); + const Array4& nc_pb_arr = NC_PB_fab[idx].const_array(); + const Array4& nc_p_arr = NC_P_fab[idx].const_array(); + + ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Base state needs ghost cells filled, protect FAB access + int ii = std::max(i , dom_lo.x); + ii = std::min(ii, dom_hi.x); + int jj = std::max(j , dom_lo.y); + jj = std::min(jj, dom_hi.y); + int kk = std::max(k , dom_lo.z); + kk = std::min(kk, dom_hi.z); + + // Base plus perturbational pressure + Real Ptot = nc_pb_arr(ii,jj,kk) + nc_p_arr(ii,jj,kk); + + // Compute pressure from EOS + Real Qv = (moisture_type != MoistureType::None) ? + cons_arr(ii,jj,kk,RhoQ1_comp) / cons_arr(ii,jj,kk,Rho_comp) : 0.0; + Real RT = cons_arr(ii,jj,kk,RhoTheta_comp); + Real P_eos = getPgivenRTh(RT, Qv); + Real DelP = std::fabs(Ptot - P_eos); + AMREX_ASSERT_WITH_MESSAGE((DelP < 1.0), "Initial state is inconsistent with EOS!"); + + // Compute rhse + Real Rhse_Sum = cons_arr(ii,jj,kk,Rho_comp); + for (int q_offset(0); q_offset khi) { - Real z_khi = 0.25 * ( nc_ph_arr (ii,jj ,khi ) + nc_ph_arr (ii-1,jj ,khi ) + - nc_ph_arr (ii,jj-1,khi ) + nc_ph_arr (ii-1,jj-1,khi) + - nc_phb_arr(ii,jj ,khi ) + nc_phb_arr(ii-1,jj ,khi ) + - nc_phb_arr(ii,jj-1,khi ) + nc_phb_arr(ii-1,jj-1,khi) ) / CONST_GRAV; Real z_khim1 = 0.25 * ( nc_ph_arr (ii,jj ,khi-1) + nc_ph_arr (ii-1,jj ,khi-1) + 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) + diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index af2df7443..7c9887e46 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -431,6 +431,11 @@ void erf_slow_rhs_post (int level, int finest_level, AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0); amrex::Real oma = 1.0 - alpha; + /* + // UNIT TEST DEBUG + oma = 1.0; alpha = 0.0; + */ + // Boundary data at fixed time intervals const auto& bdatxlo_n = bdy_data_xlo[n_time ][WRFBdyVars::QV].const_array(); const auto& bdatxlo_np1 = bdy_data_xlo[n_time+1][WRFBdyVars::QV].const_array(); @@ -540,6 +545,41 @@ void erf_slow_rhs_post (int level, int finest_level, arr_xlo, arr_xhi, arr_ylo, arr_yhi, new_cons, cell_rhs); } + + /* + // UNIT TEST DEBUG + compute_interior_ghost_bxs_xy(tbx, domain, width+1, 0, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi); + ParallelFor(bx_xlo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (arr_xlo(i,j,k) != new_cons(i,j,k,RhoQ1_comp)) { + amrex::Print() << "ERROR XLO: " << RhoQ1_comp << ' ' << IntVect(i,j,k) << "\n"; + exit(0); + } + }); + ParallelFor(bx_xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (arr_xhi(i,j,k) != new_cons(i,j,k,RhoQ1_comp)) { + amrex::Print() << "ERROR XHI: " << RhoQ1_comp<< ' ' << IntVect(i,j,k) << "\n"; + exit(0); + } + }); + ParallelFor(bx_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (arr_ylo(i,j,k) != new_cons(i,j,k,RhoQ1_comp)) { + amrex::Print() << "ERROR YLO: " << RhoQ1_comp << ' ' << IntVect(i,j,k) << "\n"; + exit(0); + } + }); + ParallelFor(bx_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (arr_yhi(i,j,k) != new_cons(i,j,k,RhoQ1_comp)) { + amrex::Print() << "ERROR YHI: " << RhoQ1_comp << ' ' << IntVect(i,j,k) << "\n"; + exit(0); + } + }); + */ } // moist_set_rhs #endif diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index 5edcc03c4..c5ae536a0 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -204,24 +204,14 @@ #ifdef ERF_USE_NETCDF // Populate RHS for relaxation zones if using real bcs - if (use_real_bcs) { - if (((init_type == "real") || (init_type == "metgrid")) && level == 0) { - int width, set_width; - if (init_type == "real") { - width = wrfbdy_width; - set_width = wrfbdy_set_width; - } else if (init_type == "metgrid") { - width = metgrid_bdy_width; - set_width = metgrid_bdy_set_width; - } - - if (width>0) - wrfbdy_compute_interior_ghost_rhs(init_type, - bdy_time_interval, start_bdy_time, new_stage_time, slow_dt, - width, set_width, fine_geom, - S_rhs, S_old, S_data, - bdy_data_xlo, bdy_data_xhi, - bdy_data_ylo, bdy_data_yhi); + if (use_real_bcs && (level == 0)) { + if (real_width>0) { + realbdy_compute_interior_ghost_rhs(init_type, + bdy_time_interval, start_bdy_time, new_stage_time, slow_dt, + real_width, real_set_width, fine_geom, + S_rhs, S_old, S_data, + bdy_data_xlo, bdy_data_xhi, + bdy_data_ylo, bdy_data_yhi); } } #endif @@ -283,20 +273,13 @@ Real slow_dt = new_stage_time - old_step_time; #if defined(ERF_USE_NETCDF) - int width, set_width; bool moist_set_rhs = false; - if ( use_real_bcs && (solverChoice.moisture_type != MoistureType::None) && (level==0) ) + if ( use_real_bcs && + (level==0) && + (real_set_width > 0) && + (solverChoice.moisture_type != MoistureType::None) ) { - if (init_type=="real") { - width = wrfbdy_width; - set_width = wrfbdy_set_width; - if (wrfbdy_set_width > 0) moist_set_rhs = true; - } - if (init_type=="metgrid") { - width = metgrid_bdy_width; - set_width = metgrid_bdy_set_width; - if (metgrid_bdy_set_width > 0) moist_set_rhs = true; - } + moist_set_rhs = true; } #endif @@ -327,7 +310,7 @@ mapfac_m[level], mapfac_u[level], mapfac_v[level], #if defined(ERF_USE_NETCDF) moist_set_rhs, bdy_time_interval, start_bdy_time, new_stage_time, - width, set_width, + real_width, real_set_width, bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, #endif fr_as_crse, fr_as_fine @@ -343,7 +326,7 @@ mapfac_m[level], mapfac_u[level], mapfac_v[level], #if defined(ERF_USE_NETCDF) moist_set_rhs, bdy_time_interval, start_bdy_time, new_stage_time, - width, set_width, + real_width, real_set_width, bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, #endif fr_as_crse, fr_as_fine @@ -385,24 +368,14 @@ #ifdef ERF_USE_NETCDF // Populate RHS for relaxation zones if using real bcs - if (use_real_bcs) { - if (((init_type == "real") || (init_type == "metgrid")) && level == 0) { - int width, set_width; - if (init_type == "real") { - width = wrfbdy_width; - set_width = wrfbdy_set_width; - } else if (init_type == "metgrid") { - width = metgrid_bdy_width; - set_width = metgrid_bdy_set_width; - } - - if (width>0) - wrfbdy_compute_interior_ghost_rhs(init_type, - bdy_time_interval, start_bdy_time, new_stage_time, slow_dt, - width, set_width, fine_geom, - S_rhs, S_old, S_data, - bdy_data_xlo, bdy_data_xhi, - bdy_data_ylo, bdy_data_yhi); + if (use_real_bcs && (level == 0)) { + if (real_width>0) { + realbdy_compute_interior_ghost_rhs(init_type, + bdy_time_interval, start_bdy_time, new_stage_time, slow_dt, + real_width, real_set_width, fine_geom, + S_rhs, S_old, S_data, + bdy_data_xlo, bdy_data_xhi, + bdy_data_ylo, bdy_data_yhi); } } #endif diff --git a/Source/Utils/InteriorGhostCells.cpp b/Source/Utils/InteriorGhostCells.cpp index 8c3a5bdd8..835c3a25a 100644 --- a/Source/Utils/InteriorGhostCells.cpp +++ b/Source/Utils/InteriorGhostCells.cpp @@ -105,21 +105,21 @@ compute_interior_ghost_bxs_xy (const Box& bx, * @param[in] start_bdy_time time of the first boundary data read in */ void -wrfbdy_compute_interior_ghost_rhs (const std::string& init_type, - const Real& bdy_time_interval, - const Real& start_bdy_time, - const Real& time, - const Real& delta_t, - int width, - int set_width, - const Geometry& geom, - Vector& S_rhs, - Vector& S_old_data, - Vector& S_cur_data, - Vector>& bdy_data_xlo, - Vector>& bdy_data_xhi, - Vector>& bdy_data_ylo, - Vector>& bdy_data_yhi) +realbdy_compute_interior_ghost_rhs (const std::string& init_type, + const Real& bdy_time_interval, + const Real& start_bdy_time, + const Real& time, + const Real& delta_t, + int width, + int set_width, + const Geometry& geom, + Vector& S_rhs, + Vector& S_old_data, + Vector& S_cur_data, + Vector>& bdy_data_xlo, + Vector>& bdy_data_xhi, + Vector>& bdy_data_ylo, + Vector>& bdy_data_yhi) { BL_PROFILE_REGION("wrfbdy_compute_interior_ghost_RHS()"); @@ -143,6 +143,11 @@ wrfbdy_compute_interior_ghost_rhs (const std::string& init_type, AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0); amrex::Real oma = 1.0 - alpha; + /* + // UNIT TEST DEBUG + oma = 1.0; alpha = 0.0; + */ + // Temporary FABs for storage (owned/filled on all ranks) FArrayBox U_xlo, U_xhi, U_ylo, U_yhi; FArrayBox V_xlo, V_xhi, V_ylo, V_yhi; @@ -156,20 +161,12 @@ wrfbdy_compute_interior_ghost_rhs (const std::string& init_type, // Variable icomp map Vector comp_map = {0, 0, Rho_comp, RhoTheta_comp}; - int BdyEnd, ivarU, ivarV, ivarR, ivarT; - if (init_type == "real") { - BdyEnd = WRFBdyVars::NumTypes-3; - ivarU = WRFBdyVars::U; - ivarV = WRFBdyVars::V; - ivarR = WRFBdyVars::R; - ivarT = WRFBdyVars::T; - } else if (init_type == "metgrid") { - BdyEnd = MetGridBdyVars::NumTypes-1; - ivarU = MetGridBdyVars::U; - ivarV = MetGridBdyVars::V; - ivarR = MetGridBdyVars::R; - ivarT = MetGridBdyVars::T; - } + // Indices + int BdyEnd = RealBdyVars::NumTypes-1; + int ivarU = RealBdyVars::U; + int ivarV = RealBdyVars::V; + int ivarR = RealBdyVars::R; + int ivarT = RealBdyVars::T; // Size the FABs @@ -501,6 +498,41 @@ wrfbdy_compute_interior_ghost_rhs (const std::string& init_type, tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi, arr_xlo, arr_xhi, arr_ylo, arr_yhi, data_arr, rhs_arr); + + /* + // UNIT TEST DEBUG + compute_interior_ghost_bxs_xy(tbx, domain, width+1, 0, + tbx_xlo, tbx_xhi, + tbx_ylo, tbx_yhi); + ParallelFor(tbx_xlo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (arr_xlo(i,j,k) != data_arr(i,j,k,icomp)) { + amrex::Print() << "ERROR XLO: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n"; + exit(0); + } + }); + ParallelFor(tbx_xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (arr_xhi(i,j,k) != data_arr(i,j,k,icomp)) { + amrex::Print() << "ERROR XHI: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n"; + exit(0); + } + }); + ParallelFor(tbx_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (arr_ylo(i,j,k) != data_arr(i,j,k,icomp)) { + amrex::Print() << "ERROR YLO: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n"; + exit(0); + } + }); + ParallelFor(tbx_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (arr_yhi(i,j,k) != data_arr(i,j,k,icomp)) { + amrex::Print() << "ERROR YHI: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n"; + exit(0); + } + }); + */ } // mfi } // ivar } // width diff --git a/Source/Utils/Utils.H b/Source/Utils/Utils.H index 68f41a690..8b9549cf6 100644 --- a/Source/Utils/Utils.H +++ b/Source/Utils/Utils.H @@ -67,21 +67,21 @@ void compute_interior_ghost_bxs_xy (const amrex::Box& bx, /* * Compute relaxation region RHS with wrfbdy */ -void wrfbdy_compute_interior_ghost_rhs (const std::string& init_type, - const amrex::Real& bdy_time_interval, - const amrex::Real& start_bdy_time, - const amrex::Real& time, - const amrex::Real& delta_t, - int width, - int set_width, - const amrex::Geometry& geom, - amrex::Vector& S_rhs, - amrex::Vector& S_old_data, - amrex::Vector& S_cur_data, - amrex::Vector>& bdy_data_xlo, - amrex::Vector>& bdy_data_xhi, - amrex::Vector>& bdy_data_ylo, - amrex::Vector>& bdy_data_yhi); +void realbdy_compute_interior_ghost_rhs (const std::string& init_type, + const amrex::Real& bdy_time_interval, + const amrex::Real& start_bdy_time, + const amrex::Real& time, + const amrex::Real& delta_t, + int width, + int set_width, + const amrex::Geometry& geom, + amrex::Vector& S_rhs, + amrex::Vector& S_old_data, + amrex::Vector& S_cur_data, + amrex::Vector>& bdy_data_xlo, + amrex::Vector>& bdy_data_xhi, + amrex::Vector>& bdy_data_ylo, + amrex::Vector>& bdy_data_yhi); /* * Compute relaxation region RHS at fine-crse interface