diff --git a/Exec/DevTests/MetGrid/inputs_fitch_flowmas b/Exec/DevTests/MetGrid/inputs_fitch_flowmas index 0f4d9ae8b..375b57458 100644 --- a/Exec/DevTests/MetGrid/inputs_fitch_flowmas +++ b/Exec/DevTests/MetGrid/inputs_fitch_flowmas @@ -34,12 +34,12 @@ amr.max_level = 0 # maximum level number allowed # CHECKPOINT FILES erf.check_file = chk # root name of checkpoint file -erf.check_int = 1 # number of timesteps between checkpoints +erf.check_int = 25 # 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_int_1 = 25 # number of timesteps between plotfiles +erf.plot_vars_1 = density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres # SOLVER CHOICE erf.alpha_T = 1.0 @@ -50,10 +50,10 @@ erf.molec_diff_type = "None" erf.les_type = "Smagorinsky" erf.Cs = 0.1 -erf.moisture_model = "Kessler" +erf.moisture_model = "NullMoist" #"Kessler" # INITIALIZATION WITH ATM DATA -erf.metgrid_bdy_width = 5 +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" diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index cd5c53e64..2ad2e41e5 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -198,21 +198,55 @@ } #ifdef ERF_USE_NETCDF - // Populate RHS for relaxation zones if using real bcs + // Populate RHS for set zones if using real bcs if (use_real_bcs) { if (((init_type == "real") || (init_type == "metgrid")) && (level==0)) { - int width, set_width; + int 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; } - 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_data, bdy_data_xlo, bdy_data_xhi, - bdy_data_ylo, bdy_data_yhi); +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(S_rhs[IntVar::cons],amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const Array4 rhs_cons = S_rhs[IntVar::cons].array(mfi); + const Array4 rhs_xmom = S_rhs[IntVar::xmom].array(mfi); + const Array4 rhs_ymom = S_rhs[IntVar::ymom].array(mfi); + const Array4 rhs_zmom = S_rhs[IntVar::zmom].array(mfi); + { + Box tbx = mfi.tilebox(); + Box domain = geom[level].Domain(); + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi); + wrfbdy_zero_rhs_in_set_region(Rho_comp, 2, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_cons); + } + + { + Box tbx = mfi.tilebox(IntVect(1,0,0)); + Box domain = geom[level].Domain(); + domain.convert(S_rhs[IntVar::xmom].boxArray().ixType()); + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi); + wrfbdy_zero_rhs_in_set_region(0, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_xmom); + } + + { + Box tbx = mfi.tilebox(IntVect(0,1,0)); + Box domain = geom[level].Domain(); + domain.convert(S_rhs[IntVar::ymom].boxArray().ixType()); + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi); + wrfbdy_zero_rhs_in_set_region(0, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_ymom); + } + } // mfi } } #endif @@ -332,6 +366,29 @@ fr_as_crse, fr_as_fine ); } + +#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 > set_width + 1) + wrfbdy_compute_interior_ghost_rhs(init_type, + bdy_time_interval, start_bdy_time, new_stage_time, + width, set_width, fine_geom, S_new, + bdy_data_xlo, bdy_data_xhi, + bdy_data_ylo, bdy_data_yhi); + } + } +#endif + }; // end slow_rhs_fun_post #ifdef ERF_USE_POISSON_SOLVE @@ -368,20 +425,56 @@ dptr_rayleigh_thetabar); #ifdef ERF_USE_NETCDF - // Populate RHS for relaxation zones - 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; + // Populate RHS for set zones if using real bcs + if (use_real_bcs) { + if (((init_type == "real") || (init_type == "metgrid")) && (level==0)) { + int set_width; + if (init_type == "real") { + set_width = wrfbdy_set_width; + } else if (init_type == "metgrid") { + set_width = metgrid_bdy_set_width; + } +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(S_rhs[IntVar::cons],amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const Array4 rhs_cons = S_rhs[IntVar::cons].array(mfi); + const Array4 rhs_xmom = S_rhs[IntVar::xmom].array(mfi); + const Array4 rhs_ymom = S_rhs[IntVar::ymom].array(mfi); + const Array4 rhs_zmom = S_rhs[IntVar::zmom].array(mfi); + { + Box tbx = mfi.tilebox(); + Box domain = geom[level].Domain(); + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi); + wrfbdy_zero_rhs_in_set_region(Rho_comp, 2, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_cons); + } + + { + Box tbx = mfi.tilebox(IntVect(1,0,0)); + Box domain = geom[level].Domain(); + domain.convert(S_rhs[IntVar::xmom].boxArray().ixType()); + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi); + wrfbdy_zero_rhs_in_set_region(0, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_xmom); + } + + { + Box tbx = mfi.tilebox(IntVect(0,1,0)); + Box domain = geom[level].Domain(); + domain.convert(S_rhs[IntVar::ymom].boxArray().ixType()); + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi); + wrfbdy_zero_rhs_in_set_region(0, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_ymom); + } + } // mfi } - wrfbdy_compute_interior_ghost_rhs(init_type, bdy_time_interval, start_bdy_time, new_stage_time, slow_dt, - width-1, set_width, fine_geom, - S_rhs, 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 e2a8974b0..ec820b818 100644 --- a/Source/Utils/InteriorGhostCells.cpp +++ b/Source/Utils/InteriorGhostCells.cpp @@ -109,11 +109,9 @@ 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_data, Vector>& bdy_data_xlo, Vector>& bdy_data_xhi, @@ -122,374 +120,316 @@ wrfbdy_compute_interior_ghost_rhs (const std::string& init_type, { BL_PROFILE_REGION("wrfbdy_compute_interior_ghost_RHS()"); - // Zero RHS in set region - //========================================================== -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(S_rhs[IntVar::cons],amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const Array4 rhs_cons = S_rhs[IntVar::cons].array(mfi); - const Array4 rhs_xmom = S_rhs[IntVar::xmom].array(mfi); - const Array4 rhs_ymom = S_rhs[IntVar::ymom].array(mfi); - const Array4 rhs_zmom = S_rhs[IntVar::zmom].array(mfi); - { - Box tbx = mfi.tilebox(); - Box domain = geom.Domain(); - Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; - compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0, - bx_xlo, bx_xhi, - bx_ylo, bx_yhi); - wrfbdy_zero_rhs_in_set_region(Rho_comp, 2, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_cons); - } + // Last cell is a ghost cell + width -= 1; - { - Box tbx = mfi.tilebox(IntVect(1,0,0)); - Box domain = geom.Domain(); - domain.convert(S_rhs[IntVar::xmom].boxArray().ixType()); - Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; - compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0, - bx_xlo, bx_xhi, - bx_ylo, bx_yhi); - wrfbdy_zero_rhs_in_set_region(0, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_xmom); - } + // Relaxation constants + Real F1 = 1./(10.); + Real F2 = 1./(50.); + + // 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; + + // 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; + FArrayBox R_xlo, R_xhi, R_ylo, R_yhi; + FArrayBox T_xlo, T_xhi, T_ylo, T_yhi; + + // Variable index map (WRFBdyVars -> Vars) + Vector var_map = {Vars::xvel, Vars::yvel, Vars::cons, Vars::cons}; + Vector ivar_map = {IntVar::xmom, IntVar::ymom, IntVar::cons, IntVar::cons}; + + // 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; + } - { - Box tbx = mfi.tilebox(IntVect(0,1,0)); - Box domain = geom.Domain(); - domain.convert(S_rhs[IntVar::ymom].boxArray().ixType()); - Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; - compute_interior_ghost_bxs_xy(tbx, domain, set_width, 0, - bx_xlo, bx_xhi, - bx_ylo, bx_yhi); - wrfbdy_zero_rhs_in_set_region(0, 1, bx_xlo, bx_xhi, bx_ylo, bx_yhi, rhs_ymom); - } - } // mfi - - if (width>set_width+1) { - // Last cell is a ghost cell - width -= 1; - - // Relaxation constants - Real F1 = 1./(10.*delta_t); - Real F2 = 1./(50.*delta_t); - - // 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; - - // 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; - FArrayBox R_xlo, R_xhi, R_ylo, R_yhi; - FArrayBox T_xlo, T_xhi, T_ylo, T_yhi; - - // Variable index map (WRFBdyVars -> Vars) - Vector var_map = {Vars::xvel, Vars::yvel, Vars::cons, Vars::cons}; - Vector ivar_map = {IntVar::xmom, IntVar::ymom, IntVar::cons, IntVar::cons}; - - // 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; - } + // Size the FABs + //========================================================== + for (int ivar(ivarU); ivar < BdyEnd; ivar++) { + int var_idx = var_map[ivar]; + Box domain = geom.Domain(); + domain.convert(S_data[var_idx].boxArray().ixType()); + + // Grown domain to get the 4 halo boxes w/ ghost cells + // NOTE: 2 ghost cells needed for U -> rho*U + IntVect ng_vect{2,2,0}; + Box gdom(domain); gdom.grow(ng_vect); + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + compute_interior_ghost_bxs_xy(gdom, domain, width, set_width, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi, + ng_vect, true); // Size the FABs - //========================================================== - for (int ivar(ivarU); ivar < BdyEnd; ivar++) - { - int var_idx = var_map[ivar]; - Box domain = geom.Domain(); - domain.convert(S_data[var_idx].boxArray().ixType()); - - // Grown domain to get the 4 halo boxes w/ ghost cells - // NOTE: 2 ghost cells needed for U -> rho*U - IntVect ng_vect{2,2,0}; - Box gdom(domain); gdom.grow(ng_vect); - Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; - compute_interior_ghost_bxs_xy(gdom, domain, width, set_width, - bx_xlo, bx_xhi, - bx_ylo, bx_yhi, - ng_vect, true); - - // Size the FABs - if (ivar == ivarU) { - U_xlo.resize(bx_xlo,1); U_xhi.resize(bx_xhi,1); - U_ylo.resize(bx_ylo,1); U_yhi.resize(bx_yhi,1); - } else if (ivar == ivarV) { - V_xlo.resize(bx_xlo,1); V_xhi.resize(bx_xhi,1); - V_ylo.resize(bx_ylo,1); V_yhi.resize(bx_yhi,1); - } else if (ivar == ivarR) { - R_xlo.resize(bx_xlo,1); R_xhi.resize(bx_xhi,1); - R_ylo.resize(bx_ylo,1); R_yhi.resize(bx_yhi,1); - } else if (ivar == ivarT){ - T_xlo.resize(bx_xlo,1); T_xhi.resize(bx_xhi,1); - T_ylo.resize(bx_ylo,1); T_yhi.resize(bx_yhi,1); - } else { - continue; - } - } // ivar + if (ivar == ivarU) { + U_xlo.resize(bx_xlo,1); U_xhi.resize(bx_xhi,1); + U_ylo.resize(bx_ylo,1); U_yhi.resize(bx_yhi,1); + } else if (ivar == ivarV) { + V_xlo.resize(bx_xlo,1); V_xhi.resize(bx_xhi,1); + V_ylo.resize(bx_ylo,1); V_yhi.resize(bx_yhi,1); + } else if (ivar == ivarR) { + R_xlo.resize(bx_xlo,1); R_xhi.resize(bx_xhi,1); + R_ylo.resize(bx_ylo,1); R_yhi.resize(bx_yhi,1); + } else if (ivar == ivarT){ + T_xlo.resize(bx_xlo,1); T_xhi.resize(bx_xhi,1); + T_ylo.resize(bx_ylo,1); T_yhi.resize(bx_yhi,1); + } else { + continue; + } + } // ivar - // Elixir to avoid destruction before kernel end - Elixir U_xlo_eli = U_xlo.elixir(); Elixir U_xhi_eli = U_xhi.elixir(); - Elixir U_ylo_eli = U_ylo.elixir(); Elixir U_yhi_eli = U_yhi.elixir(); + // Elixir to avoid destruction before kernel end + Elixir U_xlo_eli = U_xlo.elixir(); Elixir U_xhi_eli = U_xhi.elixir(); + Elixir U_ylo_eli = U_ylo.elixir(); Elixir U_yhi_eli = U_yhi.elixir(); - Elixir V_xlo_eli = V_xlo.elixir(); Elixir V_xhi_eli = V_xhi.elixir(); - Elixir V_ylo_eli = V_ylo.elixir(); Elixir V_yhi_eli = V_yhi.elixir(); + Elixir V_xlo_eli = V_xlo.elixir(); Elixir V_xhi_eli = V_xhi.elixir(); + Elixir V_ylo_eli = V_ylo.elixir(); Elixir V_yhi_eli = V_yhi.elixir(); - Elixir R_xlo_eli = R_xlo.elixir(); Elixir R_xhi_eli = R_xhi.elixir(); - Elixir R_ylo_eli = R_ylo.elixir(); Elixir R_yhi_eli = R_yhi.elixir(); + Elixir R_xlo_eli = R_xlo.elixir(); Elixir R_xhi_eli = R_xhi.elixir(); + Elixir R_ylo_eli = R_ylo.elixir(); Elixir R_yhi_eli = R_yhi.elixir(); - Elixir T_xlo_eli = T_xlo.elixir(); Elixir T_xhi_eli = T_xhi.elixir(); - Elixir T_ylo_eli = T_ylo.elixir(); Elixir T_yhi_eli = T_yhi.elixir(); + Elixir T_xlo_eli = T_xlo.elixir(); Elixir T_xhi_eli = T_xhi.elixir(); + Elixir T_ylo_eli = T_ylo.elixir(); Elixir T_yhi_eli = T_yhi.elixir(); - // Populate FABs from boundary interpolation - //========================================================== - for (int ivar(ivarU); ivar < BdyEnd; ivar++) - { - int var_idx = var_map[ivar]; - Box domain = geom.Domain(); - domain.convert(S_data[var_idx].boxArray().ixType()); - const auto& dom_lo = lbound(domain); - const auto& dom_hi = ubound(domain); - - // Grown domain to get the 4 halo boxes w/ ghost cells - // NOTE: 2 ghost cells needed for U -> rho*U - IntVect ng_vect{2,2,0}; - Box gdom(domain); gdom.grow(ng_vect); - Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; - compute_interior_ghost_bxs_xy(gdom, domain, width, set_width, - bx_xlo, bx_xhi, - bx_ylo, bx_yhi, - ng_vect, true); + // Populate FABs from boundary interpolation + //========================================================== + for (int ivar(ivarU); ivar < BdyEnd; ivar++) { + int var_idx = var_map[ivar]; + Box domain = geom.Domain(); + domain.convert(S_data[var_idx].boxArray().ixType()); + const auto& dom_lo = lbound(domain); + const auto& dom_hi = ubound(domain); + + // Grown domain to get the 4 halo boxes w/ ghost cells + // NOTE: 2 ghost cells needed for U -> rho*U + IntVect ng_vect{2,2,0}; + Box gdom(domain); gdom.grow(ng_vect); + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + compute_interior_ghost_bxs_xy(gdom, domain, width, set_width, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi, + ng_vect, true); + + Array4 arr_xlo; Array4 arr_xhi; + Array4 arr_ylo; Array4 arr_yhi; + if (ivar == ivarU) { + arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array(); + arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array(); + } else if (ivar == ivarV) { + arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array(); + arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array(); + } else if (ivar == ivarR) { + arr_xlo = R_xlo.array(); arr_xhi = R_xhi.array(); + arr_ylo = R_ylo.array(); arr_yhi = R_yhi.array(); + } else if (ivar == ivarT){ + arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array(); + arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array(); + } else { + continue; + } - Array4 arr_xlo; Array4 arr_xhi; - Array4 arr_ylo; Array4 arr_yhi; - if (ivar == ivarU) { - arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array(); - arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array(); - } else if (ivar == ivarV) { - arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array(); - arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array(); - } else if (ivar == ivarR) { - arr_xlo = R_xlo.array(); arr_xhi = R_xhi.array(); - arr_ylo = R_ylo.array(); arr_yhi = R_yhi.array(); - } else if (ivar == ivarT){ - arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array(); - arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array(); - } else { - continue; - } + // Boundary data at fixed time intervals + 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(); + + // Populate with interpolation (protect from ghost cells) + ParallelFor(bx_xlo, bx_xhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int ii = std::max(i , dom_lo.x); + ii = std::min(ii, dom_lo.x+width); + int jj = std::max(j , dom_lo.y); + jj = std::min(jj, dom_hi.y); + arr_xlo(i,j,k) = oma * bdatxlo_n (ii,jj,k,0) + + alpha * bdatxlo_np1(ii,jj,k,0); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int ii = std::max(i , dom_hi.x-width); + ii = std::min(ii, dom_hi.x); + int jj = std::max(j , dom_lo.y); + jj = std::min(jj, dom_hi.y); + arr_xhi(i,j,k) = oma * bdatxhi_n (ii,jj,k,0) + + alpha * bdatxhi_np1(ii,jj,k,0); + }); + + ParallelFor(bx_ylo, bx_yhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + 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_lo.y+width); + arr_ylo(i,j,k) = oma * bdatylo_n (ii,jj,k,0) + + alpha * bdatylo_np1(ii,jj,k,0); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int ii = std::max(i , dom_lo.x); + ii = std::min(ii, dom_hi.x); + int jj = std::max(j , dom_hi.y-width); + jj = std::min(jj, dom_hi.y); + arr_yhi(i,j,k) = oma * bdatyhi_n (ii,jj,k,0) + + alpha * bdatyhi_np1(ii,jj,k,0); + }); + } // ivar + + // Velocity to momentum + //========================================================== + for (int ivar(ivarU); ivar <= ivarV; ivar++) { + int ivar_idx = ivar_map[ivar]; + Box domain = geom.Domain(); + domain.convert(S_data[ivar_idx].boxArray().ixType()); + + // Grown domain to get the 4 halo boxes w/ ghost cells + // NOTE: 1 ghost cell needed for Laplacian + IntVect ng_vect{1,1,0}; + Box gdom(domain); gdom.grow(ng_vect); + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + compute_interior_ghost_bxs_xy(gdom, domain, width, set_width, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi, + ng_vect, true); + + Array4 rarr_xlo = R_xlo.array(); Array4 rarr_xhi = R_xhi.array(); + Array4 rarr_ylo = R_ylo.array(); Array4 rarr_yhi = R_yhi.array(); + + Array4 arr_xlo; Array4 arr_xhi; + Array4 arr_ylo; Array4 arr_yhi; + if (ivar == ivarU) { + arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array(); + arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array(); - // Boundary data at fixed time intervals - 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(); - - // Populate with interpolation (protect from ghost cells) ParallelFor(bx_xlo, bx_xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int ii = std::max(i , dom_lo.x); - ii = std::min(ii, dom_lo.x+width); - int jj = std::max(j , dom_lo.y); - jj = std::min(jj, dom_hi.y); - arr_xlo(i,j,k) = oma * bdatxlo_n (ii,jj,k,0) - + alpha * bdatxlo_np1(ii,jj,k,0); + Real rho_interp = 0.5 * ( rarr_xlo(i-1,j,k) + rarr_xlo(i,j,k) ); + arr_xlo(i,j,k) *= rho_interp; }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int ii = std::max(i , dom_hi.x-width); - ii = std::min(ii, dom_hi.x); - int jj = std::max(j , dom_lo.y); - jj = std::min(jj, dom_hi.y); - arr_xhi(i,j,k) = oma * bdatxhi_n (ii,jj,k,0) - + alpha * bdatxhi_np1(ii,jj,k,0); + Real rho_interp = 0.5 * ( rarr_xhi(i-1,j,k) + rarr_xhi(i,j,k) ); + arr_xhi(i,j,k) *= rho_interp; }); ParallelFor(bx_ylo, bx_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - 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_lo.y+width); - arr_ylo(i,j,k) = oma * bdatylo_n (ii,jj,k,0) - + alpha * bdatylo_np1(ii,jj,k,0); + Real rho_interp = 0.5 * ( rarr_ylo(i-1,j,k) + rarr_ylo(i,j,k) ); + arr_ylo(i,j,k) *= rho_interp; }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int ii = std::max(i , dom_lo.x); - ii = std::min(ii, dom_hi.x); - int jj = std::max(j , dom_hi.y-width); - jj = std::min(jj, dom_hi.y); - arr_yhi(i,j,k) = oma * bdatyhi_n (ii,jj,k,0) - + alpha * bdatyhi_np1(ii,jj,k,0); + Real rho_interp = 0.5 * ( rarr_yhi(i-1,j,k) + rarr_yhi(i,j,k) ); + arr_yhi(i,j,k) *= rho_interp; }); - } // ivar - - // Velocity to momentum - //========================================================== - for (int ivar(ivarU); ivar <= ivarV; ivar++) - { - int ivar_idx = ivar_map[ivar]; - Box domain = geom.Domain(); - domain.convert(S_data[ivar_idx].boxArray().ixType()); - - // Grown domain to get the 4 halo boxes w/ ghost cells - // NOTE: 1 ghost cell needed for Laplacian - IntVect ng_vect{1,1,0}; - Box gdom(domain); gdom.grow(ng_vect); - Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; - compute_interior_ghost_bxs_xy(gdom, domain, width, set_width, - bx_xlo, bx_xhi, - bx_ylo, bx_yhi, - ng_vect, true); - - Array4 rarr_xlo = R_xlo.array(); Array4 rarr_xhi = R_xhi.array(); - Array4 rarr_ylo = R_ylo.array(); Array4 rarr_yhi = R_yhi.array(); - - Array4 arr_xlo; Array4 arr_xhi; - Array4 arr_ylo; Array4 arr_yhi; - if (ivar == ivarU) { - arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array(); - arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array(); - - ParallelFor(bx_xlo, bx_xhi, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real rho_interp = 0.5 * ( rarr_xlo(i-1,j,k) + rarr_xlo(i,j,k) ); - arr_xlo(i,j,k) *= rho_interp; - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real rho_interp = 0.5 * ( rarr_xhi(i-1,j,k) + rarr_xhi(i,j,k) ); - arr_xhi(i,j,k) *= rho_interp; - }); - - ParallelFor(bx_ylo, bx_yhi, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real rho_interp = 0.5 * ( rarr_ylo(i-1,j,k) + rarr_ylo(i,j,k) ); - arr_ylo(i,j,k) *= rho_interp; - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real rho_interp = 0.5 * ( rarr_yhi(i-1,j,k) + rarr_yhi(i,j,k) ); - arr_yhi(i,j,k) *= rho_interp; - }); - } else { - arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array(); - arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array(); - - ParallelFor(bx_xlo, bx_xhi, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real rho_interp = 0.5 * ( rarr_xlo(i,j-1,k) + rarr_xlo(i,j,k) ); - arr_xlo(i,j,k) *= rho_interp; - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real rho_interp = 0.5 * ( rarr_xhi(i,j-1,k) + rarr_xhi(i,j,k) ); - arr_xhi(i,j,k) *= rho_interp; - }); + } else { + arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array(); + arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array(); - ParallelFor(bx_ylo, bx_yhi, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real rho_interp = 0.5 * ( rarr_ylo(i,j-1,k) + rarr_ylo(i,j,k) ); - arr_ylo(i,j,k) *= rho_interp; - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real rho_interp = 0.5 * ( rarr_yhi(i,j-1,k) + rarr_yhi(i,j,k) ); - arr_yhi(i,j,k) *= rho_interp; - }); - } - } // ivar + ParallelFor(bx_xlo, bx_xhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real rho_interp = 0.5 * ( rarr_xlo(i,j-1,k) + rarr_xlo(i,j,k) ); + arr_xlo(i,j,k) *= rho_interp; + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real rho_interp = 0.5 * ( rarr_xhi(i,j-1,k) + rarr_xhi(i,j,k) ); + arr_xhi(i,j,k) *= rho_interp; + }); + ParallelFor(bx_ylo, bx_yhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real rho_interp = 0.5 * ( rarr_ylo(i,j-1,k) + rarr_ylo(i,j,k) ); + arr_ylo(i,j,k) *= rho_interp; + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real rho_interp = 0.5 * ( rarr_yhi(i,j-1,k) + rarr_yhi(i,j,k) ); + arr_yhi(i,j,k) *= rho_interp; + }); + } + } // ivar - // Compute RHS in relaxation region - //========================================================== - for (int ivar(ivarU); ivar < BdyEnd; ivar++) - { - int ivar_idx = ivar_map[ivar]; - int icomp = comp_map[ivar]; - Box domain = geom.Domain(); - domain.convert(S_data[ivar_idx].boxArray().ixType()); - const auto& dom_hi = ubound(domain); - const auto& dom_lo = lbound(domain); + // Compute RHS in relaxation region + //========================================================== + for (int ivar(ivarU); ivar < BdyEnd; ivar++) { + int ivar_idx = ivar_map[ivar]; + int icomp = comp_map[ivar]; - // For Laplacian stencil - S_rhs[ivar_idx].FillBoundary(geom.periodicity()); + Box domain = geom.Domain(); + domain.convert(S_data[ivar_idx].boxArray().ixType()); + const auto& dom_hi = ubound(domain); + const auto& dom_lo = lbound(domain); #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for ( MFIter mfi(S_data[ivar_idx],amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box tbx = mfi.tilebox(); - Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi; - compute_interior_ghost_bxs_xy(tbx, domain, width, 0, - tbx_xlo, tbx_xhi, - tbx_ylo, tbx_yhi); - - Array4 rhs_arr; Array4 data_arr; - Array4 arr_xlo; Array4 arr_xhi; - Array4 arr_ylo; Array4 arr_yhi; - if (ivar == ivarU) { - arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array(); - arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array(); - rhs_arr = S_rhs[IntVar::xmom].array(mfi); - data_arr = S_data[IntVar::xmom].array(mfi); - } else if (ivar == ivarV) { - arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array(); - arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array(); - rhs_arr = S_rhs[IntVar::ymom].array(mfi); - data_arr = S_data[IntVar::ymom].array(mfi); - } else if (ivar == ivarR) { - arr_xlo = R_xlo.array(); arr_xhi = R_xhi.array(); - arr_ylo = R_ylo.array(); arr_yhi = R_yhi.array(); - rhs_arr = S_rhs[IntVar::cons].array(mfi); - data_arr = S_data[IntVar::cons].array(mfi); - } else if (ivar == ivarT){ - arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array(); - arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array(); - rhs_arr = S_rhs[IntVar::cons].array(mfi); - data_arr = S_data[IntVar::cons].array(mfi); - } else { - continue; - } + for ( MFIter mfi(S_data[ivar_idx],amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + Box tbx = mfi.tilebox(); + Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi; + compute_interior_ghost_bxs_xy(tbx, domain, width, 0, + tbx_xlo, tbx_xhi, + tbx_ylo, tbx_yhi); - wrfbdy_compute_laplacian_relaxation(delta_t, icomp, 1, width, set_width, dom_lo, dom_hi, F1, F2, - tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi, - arr_xlo, arr_xhi, arr_ylo, arr_yhi, - data_arr, rhs_arr); - } // mfi - } // ivar - } // width + Array4 data_arr; + Array4 arr_xlo; Array4 arr_xhi; + Array4 arr_ylo; Array4 arr_yhi; + if (ivar == ivarU) { + arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array(); + arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array(); + data_arr = S_data[IntVar::xmom].array(mfi); + } else if (ivar == ivarV) { + arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array(); + arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array(); + data_arr = S_data[IntVar::ymom].array(mfi); + } else if (ivar == ivarR) { + arr_xlo = R_xlo.array(); arr_xhi = R_xhi.array(); + arr_ylo = R_ylo.array(); arr_yhi = R_yhi.array(); + data_arr = S_data[IntVar::cons].array(mfi); + } else if (ivar == ivarT){ + arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array(); + arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array(); + data_arr = S_data[IntVar::cons].array(mfi); + } else { + continue; + } + + wrfbdy_compute_laplacian_relaxation(icomp, 1, width, set_width, dom_lo, dom_hi, F1, F2, + tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi, + arr_xlo, arr_xhi, arr_ylo, arr_yhi, + data_arr); + } // mfi + } // ivar } /** diff --git a/Source/Utils/Utils.H b/Source/Utils/Utils.H index 516317a02..cf6ff69b9 100644 --- a/Source/Utils/Utils.H +++ b/Source/Utils/Utils.H @@ -71,11 +71,9 @@ 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_data, amrex::Vector>& bdy_data_xlo, amrex::Vector>& bdy_data_xhi, @@ -166,8 +164,7 @@ wrfbdy_zero_rhs_in_set_region (const int& icomp, AMREX_GPU_HOST AMREX_FORCE_INLINE void -wrfbdy_compute_laplacian_relaxation (const amrex::Real& delta_t, - const int& icomp, +wrfbdy_compute_laplacian_relaxation (const int& icomp, const int& num_var, const int& width, const int& set_width, @@ -183,15 +180,14 @@ wrfbdy_compute_laplacian_relaxation (const amrex::Real& delta_t, const amrex::Array4& arr_xhi, const amrex::Array4& arr_ylo, const amrex::Array4& arr_yhi, - const amrex::Array4& data_arr, - const amrex::Array4& rhs_arr) + const amrex::Array4& data_arr) { // RHS computation int Spec_z = set_width; int Relax_z = width - Spec_z + 1; amrex::Real num = amrex::Real(Spec_z + Relax_z); amrex::Real denom = amrex::Real(Relax_z - 1); - amrex::Real SpecExp = 0.5; + amrex::Real SpecExp = 1.0; amrex::ParallelFor(bx_xlo, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { // Corners with x boxes @@ -199,23 +195,17 @@ wrfbdy_compute_laplacian_relaxation (const amrex::Real& delta_t, int j_hi = std::min(dom_hi.y-j,width-1); int jj = std::min(j_lo,j_hi); int n_ind = std::min(i-dom_lo.x,jj) + 1; - if (n_ind <= Spec_z) { - rhs_arr(i,j,k,n+icomp) = 0.0; - } else { - amrex::Real Factor = (num - amrex::Real(n_ind))/denom - * std::exp(-SpecExp * amrex::Real(n_ind - Spec_z - 1)); - amrex::Real d = data_arr(i ,j ,k ,n+icomp) + delta_t*rhs_arr(i , j , k ,n+icomp); - amrex::Real d_ip1 = data_arr(i+1,j ,k ,n+icomp) + delta_t*rhs_arr(i+1, j , k ,n+icomp); - amrex::Real d_im1 = data_arr(i-1,j ,k ,n+icomp) + delta_t*rhs_arr(i-1, j , k ,n+icomp); - amrex::Real d_jp1 = data_arr(i ,j+1,k ,n+icomp) + delta_t*rhs_arr(i , j+1, k ,n+icomp); - amrex::Real d_jm1 = data_arr(i ,j-1,k ,n+icomp) + delta_t*rhs_arr(i , j-1, k ,n+icomp); - amrex::Real delta = arr_xlo(i ,j ,k,n) - d; - amrex::Real delta_xp = arr_xlo(i+1,j ,k,n) - d_ip1; - amrex::Real delta_xm = arr_xlo(i-1,j ,k,n) - d_im1; - amrex::Real delta_yp = arr_xlo(i ,j+1,k,n) - d_jp1; - amrex::Real delta_ym = arr_xlo(i ,j-1,k,n) - d_jm1; + if (n_ind > Spec_z) { + amrex::Real Factor = (num - amrex::Real(n_ind))/denom + * std::exp(-SpecExp * amrex::Real(n_ind - Spec_z - 1)); + amrex::Real delta = arr_xlo(i ,j ,k,n) - data_arr(i ,j ,k ,n+icomp); + amrex::Real delta_xp = arr_xlo(i+1,j ,k,n) - data_arr(i+1,j ,k ,n+icomp); + amrex::Real delta_xm = arr_xlo(i-1,j ,k,n) - data_arr(i-1,j ,k ,n+icomp); + amrex::Real delta_yp = arr_xlo(i ,j+1,k,n) - data_arr(i ,j+1,k ,n+icomp); + amrex::Real delta_ym = arr_xlo(i ,j-1,k,n) - data_arr(i ,j-1,k ,n+icomp); amrex::Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta; - rhs_arr(i,j,k,n+icomp) += (F1*delta - F2*Laplacian) * Factor; + amrex::Real RHS_nudge = (F1*delta - F2*Laplacian) * Factor; + data_arr(i,j,k,n+icomp) += RHS_nudge; } }, bx_xhi, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept @@ -225,23 +215,17 @@ wrfbdy_compute_laplacian_relaxation (const amrex::Real& delta_t, int j_hi = std::min(dom_hi.y-j,width-1); int jj = std::min(j_lo,j_hi); int n_ind = std::min(dom_hi.x-i,jj) + 1; - if (n_ind <= Spec_z) { - rhs_arr(i,j,k,n+icomp) = 0.0; - } else { - amrex::Real Factor = (num - amrex::Real(n_ind))/denom - * std::exp(-SpecExp * amrex::Real(n_ind - Spec_z - 1)); - amrex::Real d = data_arr(i ,j ,k ,n+icomp) + delta_t*rhs_arr(i , j , k ,n+icomp); - amrex::Real d_ip1 = data_arr(i+1,j ,k ,n+icomp) + delta_t*rhs_arr(i+1, j , k ,n+icomp); - amrex::Real d_im1 = data_arr(i-1,j ,k ,n+icomp) + delta_t*rhs_arr(i-1, j , k ,n+icomp); - amrex::Real d_jp1 = data_arr(i ,j+1,k ,n+icomp) + delta_t*rhs_arr(i , j+1, k ,n+icomp); - amrex::Real d_jm1 = data_arr(i ,j-1,k ,n+icomp) + delta_t*rhs_arr(i , j-1, k ,n+icomp); - amrex::Real delta = arr_xhi(i ,j ,k,n) - d; - amrex::Real delta_xp = arr_xhi(i+1,j ,k,n) - d_ip1; - amrex::Real delta_xm = arr_xhi(i-1,j ,k,n) - d_im1; - amrex::Real delta_yp = arr_xhi(i ,j+1,k,n) - d_jp1; - amrex::Real delta_ym = arr_xhi(i ,j-1,k,n) - d_jm1; + if (n_ind > Spec_z) { + amrex::Real Factor = (num - amrex::Real(n_ind))/denom + * std::exp(-SpecExp * amrex::Real(n_ind - Spec_z - 1)); + amrex::Real delta = arr_xhi(i ,j ,k,n) - data_arr(i ,j ,k ,n+icomp); + amrex::Real delta_xp = arr_xhi(i+1,j ,k,n) - data_arr(i+1,j ,k ,n+icomp); + amrex::Real delta_xm = arr_xhi(i-1,j ,k,n) - data_arr(i-1,j ,k ,n+icomp); + amrex::Real delta_yp = arr_xhi(i ,j+1,k,n) - data_arr(i ,j+1,k ,n+icomp); + amrex::Real delta_ym = arr_xhi(i ,j-1,k,n) - data_arr(i ,j-1,k ,n+icomp); amrex::Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta; - rhs_arr(i,j,k,n+icomp) += (F1*delta - F2*Laplacian) * Factor; + amrex::Real RHS_nudge = (F1*delta - F2*Laplacian) * Factor; + data_arr(i,j,k,n+icomp) += RHS_nudge; } }); @@ -249,46 +233,34 @@ wrfbdy_compute_laplacian_relaxation (const amrex::Real& delta_t, { // No corners for y boxes int n_ind = j - dom_lo.y + 1; - if (n_ind <= Spec_z) { - rhs_arr(i,j,k,n+icomp) = 0.0; - } else { - amrex::Real Factor = (num - amrex::Real(n_ind))/denom - * std::exp(-SpecExp * amrex::Real(n_ind - Spec_z - 1)); - amrex::Real d = data_arr(i ,j ,k ,n+icomp) + delta_t*rhs_arr(i , j , k ,n+icomp); - amrex::Real d_ip1 = data_arr(i+1,j ,k ,n+icomp) + delta_t*rhs_arr(i+1, j , k ,n+icomp); - amrex::Real d_im1 = data_arr(i-1,j ,k ,n+icomp) + delta_t*rhs_arr(i-1, j , k ,n+icomp); - amrex::Real d_jp1 = data_arr(i ,j+1,k ,n+icomp) + delta_t*rhs_arr(i , j+1, k ,n+icomp); - amrex::Real d_jm1 = data_arr(i ,j-1,k ,n+icomp) + delta_t*rhs_arr(i , j-1, k ,n+icomp); - amrex::Real delta = arr_ylo(i ,j ,k,n) - d; - amrex::Real delta_xp = arr_ylo(i+1,j ,k,n) - d_ip1; - amrex::Real delta_xm = arr_ylo(i-1,j ,k,n) - d_im1; - amrex::Real delta_yp = arr_ylo(i ,j+1,k,n) - d_jp1; - amrex::Real delta_ym = arr_ylo(i ,j-1,k,n) - d_jm1; + if (n_ind > Spec_z) { + amrex::Real Factor = (num - amrex::Real(n_ind))/denom + * std::exp(-SpecExp * amrex::Real(n_ind - Spec_z - 1)); + amrex::Real delta = arr_ylo(i ,j ,k,n) - data_arr(i ,j ,k ,n+icomp); + amrex::Real delta_xp = arr_ylo(i+1,j ,k,n) - data_arr(i+1,j ,k ,n+icomp); + amrex::Real delta_xm = arr_ylo(i-1,j ,k,n) - data_arr(i-1,j ,k ,n+icomp); + amrex::Real delta_yp = arr_ylo(i ,j+1,k,n) - data_arr(i ,j+1,k ,n+icomp); + amrex::Real delta_ym = arr_ylo(i ,j-1,k,n) - data_arr(i ,j-1,k ,n+icomp); amrex::Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta; - rhs_arr(i,j,k,n+icomp) += (F1*delta - F2*Laplacian) * Factor; + amrex::Real RHS_nudge = (F1*delta - F2*Laplacian) * Factor; + data_arr(i,j,k,n+icomp) += RHS_nudge; } }, bx_yhi, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { // No corners for y boxes int n_ind = dom_hi.y - j + 1; - if (n_ind <= Spec_z) { - rhs_arr(i,j,k,n+icomp) = 0.0; - } else { - amrex::Real Factor = (num - amrex::Real(n_ind))/denom - * std::exp(-SpecExp * amrex::Real(n_ind - Spec_z - 1)); - amrex::Real d = data_arr(i ,j ,k ,n+icomp) + delta_t*rhs_arr(i , j , k ,n+icomp); - amrex::Real d_ip1 = data_arr(i+1,j ,k ,n+icomp) + delta_t*rhs_arr(i+1, j , k ,n+icomp); - amrex::Real d_im1 = data_arr(i-1,j ,k ,n+icomp) + delta_t*rhs_arr(i-1, j , k ,n+icomp); - amrex::Real d_jp1 = data_arr(i ,j+1,k ,n+icomp) + delta_t*rhs_arr(i , j+1, k ,n+icomp); - amrex::Real d_jm1 = data_arr(i ,j-1,k ,n+icomp) + delta_t*rhs_arr(i , j-1, k ,n+icomp); - amrex::Real delta = arr_yhi(i ,j ,k,n) - d; - amrex::Real delta_xp = arr_yhi(i+1,j ,k,n) - d_ip1; - amrex::Real delta_xm = arr_yhi(i-1,j ,k,n) - d_im1; - amrex::Real delta_yp = arr_yhi(i ,j+1,k,n) - d_jp1; - amrex::Real delta_ym = arr_yhi(i ,j-1,k,n) - d_jm1; + if (n_ind > Spec_z) { + amrex::Real Factor = (num - amrex::Real(n_ind))/denom + * std::exp(-SpecExp * amrex::Real(n_ind - Spec_z - 1)); + amrex::Real delta = arr_yhi(i ,j ,k,n) - data_arr(i ,j ,k ,n+icomp); + amrex::Real delta_xp = arr_yhi(i+1,j ,k,n) - data_arr(i+1,j ,k ,n+icomp); + amrex::Real delta_xm = arr_yhi(i-1,j ,k,n) - data_arr(i-1,j ,k ,n+icomp); + amrex::Real delta_yp = arr_yhi(i ,j+1,k,n) - data_arr(i ,j+1,k ,n+icomp); + amrex::Real delta_ym = arr_yhi(i ,j-1,k,n) - data_arr(i ,j-1,k ,n+icomp); amrex::Real Laplacian = delta_xp + delta_xm + delta_yp + delta_ym - 4.0*delta; - rhs_arr(i,j,k,n+icomp) += (F1*delta - F2*Laplacian) * Factor; + amrex::Real RHS_nudge = (F1*delta - F2*Laplacian) * Factor; + data_arr(i,j,k,n+icomp) += RHS_nudge; } }); }