diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index cf83e3829..5ac55fe24 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -402,7 +402,6 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/SourceTerms/ERF_MakeGradP.cpp ${SRC_DIR}/SourceTerms/ERF_MakeMomSources.cpp ${SRC_DIR}/SourceTerms/ERF_MakeSources.cpp - ${SRC_DIR}/SourceTerms/ERF_MoistSetRhs.cpp ${SRC_DIR}/SourceTerms/ERF_NumericalDiffusion.cpp ${SRC_DIR}/SourceTerms/ERF_ForestDrag.cpp ${SRC_DIR}/SourceTerms/ERF_ApplySurfaceTreatment_BulkCoeff.cpp diff --git a/Source/SourceTerms/ERF_AddMoistNudgingTerms.cpp b/Source/SourceTerms/ERF_AddMoistNudgingTerms.cpp index 0688d65cb..d69877ed6 100644 --- a/Source/SourceTerms/ERF_AddMoistNudgingTerms.cpp +++ b/Source/SourceTerms/ERF_AddMoistNudgingTerms.cpp @@ -1,6 +1,7 @@ #if defined(ERF_USE_NETCDF) #include +#include using namespace amrex; @@ -26,13 +27,14 @@ using namespace amrex; void add_moist_nudging_terms (const MultiFab& S_data, MultiFab & source, + const int n_qstate, const Real& dt, - const Real& old_stage_time_total, + const Real& time, const Real& start_bdy_time, const Real& final_bdy_time, const Real& bdy_time_interval, - const Real& bdy_factor, - int width, + const Real& nudge_factor, + const int width, const Geometry& geom, Vector>& bdy_data_xlo, Vector>& bdy_data_xhi, @@ -44,23 +46,167 @@ void add_moist_nudging_terms (const MultiFab& S_data, const Box domain = geom.Domain(); + int bdy_comp = BCVars::RhoQ1_bc_comp; + Array4 bdatxlo, bdatxhi, bdatylo, bdatyhi; + + // Relaxation constants + Real F1 = one/(nudge_factor*dt); + + // Domain bounds + const auto& dom_hi = ubound(domain); + const auto& dom_lo = lbound(domain); + auto dx = geom.CellSizeArray(); + auto ProbHi = geom.ProbHiArray(); + auto ProbLo = geom.ProbLoArray(); + + // Time interpolation + Real dT = bdy_time_interval; + + // + // Note that time (= start_time+old_stage_time) is measured as total time + // start_bdy_time and final_bdy_time are also measured as total time + // + + int n_time = static_cast( (time-start_bdy_time) / dT); + int n_time_p1 = n_time + 1; + Real alpha = ((time-start_bdy_time) - n_time * dT) / dT; + + // Do not over run the last bdy file + if (time >= final_bdy_time) { + n_time = static_cast( (final_bdy_time - start_bdy_time)/ dT); + n_time_p1 = n_time; + alpha = zero; + } + + AMREX_ALWAYS_ASSERT( alpha >= zero && alpha <= one); + Real oma = one - alpha; + + // Limiting offset + int offset = width - 1; + for (MFIter mfi(S_data,TilingIfNotGPU()); mfi.isValid(); ++mfi) { Box tbx = mfi.tilebox(); - const Array4& new_cons_const = S_data.const_array(mfi); - const Array4< Real>& src_arr = source.array(mfi); + const Array4& cons_arr = S_data.const_array(mfi); + + // We will add to source terms for moist variables and + // to source term for (rho theta) + const Array4< Real>& src_arr = source.array(mfi); + // - // Note that old_stage_time_total = start_time+old_stage_time is total time + // Note that time = start_time+old_stage_time is total time // start_bdy_time and final_bdy_time are total time // - moist_set_rhs(geom, tbx, new_cons_const, src_arr, - old_stage_time_total, dt, - start_bdy_time, final_bdy_time, bdy_time_interval, - bdy_factor, width, domain, - bdy_data_xlo, bdy_data_xhi, - bdy_data_ylo, bdy_data_yhi, - m_r2d); - } + // moist_set_rhs(geom, tbx, cons_arr, src_arr, n_qstate, + // old_stage_time_total, dt, + // start_bdy_time, final_bdy_time, bdy_time_interval, + // bdy_factor, width, domain, + // bdy_data_xlo, bdy_data_xhi, + // bdy_data_ylo, bdy_data_yhi, + // m_r2d); + + // Get bndry data + if (m_r2d) { + Vector>& bndry_data = m_r2d->interp_in_time(time); + bdatxlo = (*bndry_data[0])[0].array(); + bdatylo = (*bndry_data[1])[0].array(); + bdatxhi = (*bndry_data[3])[0].array(); + bdatyhi = (*bndry_data[4])[0].array(); + } + + // NOTE: The sizing of the temporary BDY FABS is + // GLOBAL and occurs over the entire BDY region. + + // Size the FABs + //========================================================== + // NOTE: No ghost cells, we force mask to be idx type (0,0,0) + IntVect ng_vect(0); + Box gdom(domain); gdom.grow(ng_vect); + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + realbdy_interior_bxs_xy(gdom, domain, width, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi, + ng_vect, true); + + // Temporary FABs for storage (owned/filled on all ranks) + FArrayBox QV_xlo, QV_xhi, QV_ylo, QV_yhi; + QV_xlo.resize(bx_xlo,1,The_Async_Arena()); QV_xhi.resize(bx_xhi,1,The_Async_Arena()); + QV_ylo.resize(bx_ylo,1,The_Async_Arena()); QV_yhi.resize(bx_yhi,1,The_Async_Arena()); + + // Populate FABs from bdy interpolation (primitive vars) + //========================================================== + const auto& bdatxlo_n = bdy_data_xlo[n_time ][WRFBdyVars::QV].const_array(); + const auto& bdatxlo_np1 = bdy_data_xlo[n_time_p1][WRFBdyVars::QV].const_array(); + const auto& bdatxhi_n = bdy_data_xhi[n_time ][WRFBdyVars::QV].const_array(); + const auto& bdatxhi_np1 = bdy_data_xhi[n_time_p1][WRFBdyVars::QV].const_array(); + const auto& bdatylo_n = bdy_data_ylo[n_time ][WRFBdyVars::QV].const_array(); + const auto& bdatylo_np1 = bdy_data_ylo[n_time_p1][WRFBdyVars::QV].const_array(); + const auto& bdatyhi_n = bdy_data_yhi[n_time ][WRFBdyVars::QV].const_array(); + const auto& bdatyhi_np1 = bdy_data_yhi[n_time_p1][WRFBdyVars::QV].const_array(); + + // Get Array4 of interpolated values + Array4 arr_xlo = QV_xlo.array(); Array4 arr_xhi = QV_xhi.array(); + Array4 arr_ylo = QV_ylo.array(); Array4 arr_yhi = QV_yhi.array(); + + Box gtbx = grow(tbx,ng_vect); + Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi; + realbdy_interior_bxs_xy(gtbx, domain, width, + tbx_xlo, tbx_xhi, + tbx_ylo, tbx_yhi, + ng_vect, true); + + // Populate with interpolation (protect from ghost cells) + ParallelFor(tbx_xlo, tbx_xhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int ii = std::min(std::max(i , dom_lo.x), dom_lo.x+offset); + int jj = std::min(std::max(j , dom_lo.y), dom_hi.y ); + arr_xlo(i,j,k) = (bdatxlo) ? cons_arr(i,j,k,Rho_comp) * bdatxlo(ii,jj,k,bdy_comp) : + cons_arr(i,j,k,Rho_comp) * ( oma * bdatxlo_n (ii,jj,k) + + alpha * bdatxlo_np1(ii,jj,k) ); + } , + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int ii = std::min(std::max(i , dom_hi.x-offset), dom_hi.x); + int jj = std::min(std::max(j , dom_lo.y ), dom_hi.y); + arr_xhi(i,j,k) = (bdatxhi) ? cons_arr(i,j,k,Rho_comp) * bdatxhi(ii,jj,k,bdy_comp) : + cons_arr(i,j,k,Rho_comp) * ( oma * bdatxhi_n (ii,jj,k) + + alpha * bdatxhi_np1(ii,jj,k) ); + }); + + ParallelFor(tbx_ylo, tbx_yhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int ii = std::min(std::max(i , dom_lo.x), dom_hi.x ); + int jj = std::min(std::max(j , dom_lo.y), dom_lo.y+offset); + arr_ylo(i,j,k) = (bdatylo) ? cons_arr(i,j,k,Rho_comp) * bdatylo(ii,jj,k,bdy_comp) : + cons_arr(i,j,k,Rho_comp) * ( oma * bdatylo_n (ii,jj,k) + + alpha * bdatylo_np1(ii,jj,k) ); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int ii = std::min(std::max(i , dom_lo.x ), dom_hi.x); + int jj = std::min(std::max(j , dom_hi.y-offset), dom_hi.y); + jj = std::min(jj, dom_hi.y); + arr_yhi(i,j,k) = (bdatyhi) ? cons_arr(i,j,k,Rho_comp) * bdatyhi(ii,jj,k,bdy_comp) : + cons_arr(i,j,k,Rho_comp) * ( oma * bdatyhi_n (ii,jj,k) + + alpha * bdatyhi_np1(ii,jj,k) ); + }); + + realbdy_interior_bxs_xy(tbx, domain, width, + tbx_xlo, tbx_xhi, + tbx_ylo, tbx_yhi, + ng_vect); + + // + // Add relaxation terms for moist variables and (rho theta) to existing source terms + // + realbdy_compute_relaxation(RhoQ1_comp, n_qstate, + width, dx, ProbLo, ProbHi, F1, + tbx_xlo , tbx_xhi , tbx_ylo , tbx_yhi , + arr_xlo , arr_xhi , arr_ylo , arr_yhi , + cons_arr, src_arr); + } // mfi } #endif diff --git a/Source/SourceTerms/ERF_MoistSetRhs.cpp b/Source/SourceTerms/ERF_MoistSetRhs.cpp deleted file mode 100644 index 39633b4ab..000000000 --- a/Source/SourceTerms/ERF_MoistSetRhs.cpp +++ /dev/null @@ -1,169 +0,0 @@ -#if defined(ERF_USE_NETCDF) - -#include -#include - -using namespace amrex; - -/** - * Function for setting the slow variables in the "specified" zones at the domain boundary - * -*/ - -void -moist_set_rhs (const Geometry& geom, - const Box& tbx, - const Array4& new_cons, - const Array4& cell_rhs, - const Real& time, - const Real& dt, - const Real& start_bdy_time, - const Real& final_bdy_time, - const Real& bdy_time_interval, - const Real& nudge_factor, - int width, - const Box& domain, - Vector>& bdy_data_xlo, - Vector>& bdy_data_xhi, - Vector>& bdy_data_ylo, - Vector>& bdy_data_yhi, - std::unique_ptr& m_r2d) -{ - // Get bndry data - int bdy_comp = BCVars::RhoQ1_bc_comp; - Array4 bdatxlo, bdatxhi, bdatylo, bdatyhi; - if (m_r2d) { - Vector>& bndry_data = m_r2d->interp_in_time(time); - bdatxlo = (*bndry_data[0])[0].array(); - bdatylo = (*bndry_data[1])[0].array(); - bdatxhi = (*bndry_data[3])[0].array(); - bdatyhi = (*bndry_data[4])[0].array(); - } - - // - // Note that time (= start_time+old_stage_time) is measured as total time - // start_bdy_time and final_bdy_time are also measured as total time - // - - // Relaxation constants - Real F1 = one/(nudge_factor*dt); - - // Domain bounds - const auto& dom_hi = ubound(domain); - const auto& dom_lo = lbound(domain); - auto dx = geom.CellSizeArray(); - auto ProbHi = geom.ProbHiArray(); - auto ProbLo = geom.ProbLoArray(); - - // Time interpolation - Real dT = bdy_time_interval; - - int n_time = static_cast( (time-start_bdy_time) / dT); - int n_time_p1 = n_time + 1; - Real alpha = ((time-start_bdy_time) - n_time * dT) / dT; - - // Do not over run the last bdy file - if (time >= final_bdy_time) { - n_time = static_cast( (final_bdy_time - start_bdy_time)/ dT); - n_time_p1 = n_time; - alpha = zero; - } - - AMREX_ALWAYS_ASSERT( alpha >= zero && alpha <= one); - Real oma = one - alpha; - - // NOTE: The sizing of the temporary BDY FABS is - // GLOBAL and occurs over the entire BDY region. - - // Size the FABs - //========================================================== - // NOTE: No ghost cells, we force mask to be idx type (0,0,0) - IntVect ng_vect(0); - Box gdom(domain); gdom.grow(ng_vect); - Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; - realbdy_interior_bxs_xy(gdom, domain, width, - bx_xlo, bx_xhi, - bx_ylo, bx_yhi, - ng_vect, true); - - // Temporary FABs for storage (owned/filled on all ranks) - FArrayBox QV_xlo, QV_xhi, QV_ylo, QV_yhi; - QV_xlo.resize(bx_xlo,1,The_Async_Arena()); QV_xhi.resize(bx_xhi,1,The_Async_Arena()); - QV_ylo.resize(bx_ylo,1,The_Async_Arena()); QV_yhi.resize(bx_yhi,1,The_Async_Arena()); - - // Populate FABs from bdy interpolation (primitive vars) - //========================================================== - const auto& bdatxlo_n = bdy_data_xlo[n_time ][WRFBdyVars::QV].const_array(); - const auto& bdatxlo_np1 = bdy_data_xlo[n_time_p1][WRFBdyVars::QV].const_array(); - const auto& bdatxhi_n = bdy_data_xhi[n_time ][WRFBdyVars::QV].const_array(); - const auto& bdatxhi_np1 = bdy_data_xhi[n_time_p1][WRFBdyVars::QV].const_array(); - const auto& bdatylo_n = bdy_data_ylo[n_time ][WRFBdyVars::QV].const_array(); - const auto& bdatylo_np1 = bdy_data_ylo[n_time_p1][WRFBdyVars::QV].const_array(); - const auto& bdatyhi_n = bdy_data_yhi[n_time ][WRFBdyVars::QV].const_array(); - const auto& bdatyhi_np1 = bdy_data_yhi[n_time_p1][WRFBdyVars::QV].const_array(); - - // Get Array4 of interpolated values - Array4 arr_xlo = QV_xlo.array(); Array4 arr_xhi = QV_xhi.array(); - Array4 arr_ylo = QV_ylo.array(); Array4 arr_yhi = QV_yhi.array(); - - Box gtbx = grow(tbx,ng_vect); - Box tbx_xlo, tbx_xhi, tbx_ylo, tbx_yhi; - realbdy_interior_bxs_xy(gtbx, domain, width, - tbx_xlo, tbx_xhi, - tbx_ylo, tbx_yhi, - ng_vect, true); - - // Limiting offset - int offset = width - 1; - - // Populate with interpolation (protect from ghost cells) - ParallelFor(tbx_xlo, tbx_xhi, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - int ii = std::min(std::max(i , dom_lo.x), dom_lo.x+offset); - int jj = std::min(std::max(j , dom_lo.y), dom_hi.y ); - arr_xlo(i,j,k) = (bdatxlo) ? new_cons(i,j,k,Rho_comp) * bdatxlo(ii,jj,k,bdy_comp) : - new_cons(i,j,k,Rho_comp) * ( oma * bdatxlo_n (ii,jj,k) - + alpha * bdatxlo_np1(ii,jj,k) ); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - int ii = std::min(std::max(i , dom_hi.x-offset), dom_hi.x); - int jj = std::min(std::max(j , dom_lo.y ), dom_hi.y); - arr_xhi(i,j,k) = (bdatxhi) ? new_cons(i,j,k,Rho_comp) * bdatxhi(ii,jj,k,bdy_comp) : - new_cons(i,j,k,Rho_comp) * ( oma * bdatxhi_n (ii,jj,k) - + alpha * bdatxhi_np1(ii,jj,k) ); - }); - - ParallelFor(tbx_ylo, tbx_yhi, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - int ii = std::min(std::max(i , dom_lo.x), dom_hi.x ); - int jj = std::min(std::max(j , dom_lo.y), dom_lo.y+offset); - arr_ylo(i,j,k) = (bdatylo) ? new_cons(i,j,k,Rho_comp) * bdatylo(ii,jj,k,bdy_comp) : - new_cons(i,j,k,Rho_comp) * ( oma * bdatylo_n (ii,jj,k) - + alpha * bdatylo_np1(ii,jj,k) ); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - int ii = std::min(std::max(i , dom_lo.x ), dom_hi.x); - int jj = std::min(std::max(j , dom_hi.y-offset), dom_hi.y); - jj = std::min(jj, dom_hi.y); - arr_yhi(i,j,k) = (bdatyhi) ? new_cons(i,j,k,Rho_comp) * bdatyhi(ii,jj,k,bdy_comp) : - new_cons(i,j,k,Rho_comp) * ( oma * bdatyhi_n (ii,jj,k) - + alpha * bdatyhi_np1(ii,jj,k) ); - }); - - // Compute RHS in relaxation region - //========================================================== - realbdy_interior_bxs_xy(tbx, domain, width, - tbx_xlo, tbx_xhi, - tbx_ylo, tbx_yhi, - ng_vect); - realbdy_compute_relaxation(RhoQ1_comp, 1, - width, dx, ProbLo, ProbHi, F1, - tbx_xlo , tbx_xhi , tbx_ylo , tbx_yhi , - arr_xlo , arr_xhi , arr_ylo , arr_yhi , - new_cons, cell_rhs); -} // moist_set_rhs -#endif diff --git a/Source/SourceTerms/ERF_SrcHeaders.H b/Source/SourceTerms/ERF_SrcHeaders.H index 9f7f03f63..a248e3af9 100644 --- a/Source/SourceTerms/ERF_SrcHeaders.H +++ b/Source/SourceTerms/ERF_SrcHeaders.H @@ -98,13 +98,14 @@ void make_sources (int level, int nrk, void add_moist_nudging_terms (const amrex::MultiFab& S_data, amrex::MultiFab& source, + const int n_qstate, const amrex::Real& dt, const amrex::Real& old_stage_time_total, const amrex::Real& start_bdy_time, const amrex::Real& final_bdy_time, const amrex::Real& bdy_time_interval, const amrex::Real& bdy_factor, - int width, + const int width, const amrex::Geometry& geom, amrex::Vector>& bdy_data_xlo, amrex::Vector>& bdy_data_xhi, @@ -162,13 +163,14 @@ moist_set_rhs (const amrex::Geometry& geom, const amrex::Box& tbx, const amrex::Array4& new_cons, const amrex::Array4& cell_rhs, + const int n_qstate, const amrex::Real& time, const amrex::Real& dt, const amrex::Real& start_bdy_time, const amrex::Real& final_bdy_time, const amrex::Real& bdy_time_interval, const amrex::Real& nudge_factor, - int width, + const int width, const amrex::Box& domain, amrex::Vector>& bdy_data_xlo, amrex::Vector>& bdy_data_xhi, diff --git a/Source/SourceTerms/Make.package b/Source/SourceTerms/Make.package index 1bf01c966..197c74662 100644 --- a/Source/SourceTerms/Make.package +++ b/Source/SourceTerms/Make.package @@ -12,7 +12,6 @@ CEXE_sources += ERF_ApplySurfaceTreatment_BulkCoeff.cpp ifeq ($(USE_NETCDF),TRUE) CEXE_sources += ERF_AddMoistNudgingTerms.cpp -CEXE_sources += ERF_MoistSetRhs.cpp endif CEXE_headers += ERF_NumericalDiffusion.H diff --git a/Source/TimeIntegration/ERF_TI_slow_rhs_pre.H b/Source/TimeIntegration/ERF_TI_slow_rhs_pre.H index 272575664..398f607f9 100644 --- a/Source/TimeIntegration/ERF_TI_slow_rhs_pre.H +++ b/Source/TimeIntegration/ERF_TI_slow_rhs_pre.H @@ -120,7 +120,10 @@ (solverChoice.moisture_type != MoistureType::None) ) { Real bdy_factor = solverChoice.bdy_nudge_factor; - add_moist_nudging_terms(S_data[IntVars::cons], cc_src, slow_dt, + int num_q = micro->Get_Qstate_Moist_Size() - micro->Get_Qstate_Moist_NumConc_Size(); + AMREX_ALWAYS_ASSERT(solverChoice.moisture_type != MoistureType::Morrison_NoIce && + solverChoice.moisture_type != MoistureType::SAM_NoIce); + add_moist_nudging_terms(S_data[IntVars::cons], cc_src, num_q, slow_dt, start_time+old_stage_time, start_bdy_time, final_bdy_time, bdy_time_interval, bdy_factor, real_width, geom[level], diff --git a/Source/Utils/ERF_Utils.H b/Source/Utils/ERF_Utils.H index ca63e078e..b29db461e 100644 --- a/Source/Utils/ERF_Utils.H +++ b/Source/Utils/ERF_Utils.H @@ -4,6 +4,7 @@ #include "AMReX.H" #include "AMReX_MultiFab.H" #include "AMReX_BCRec.H" +#include "ERF_Constants.H" #include "ERF_DataStruct.H" #include "ERF_IndexDefines.H" #include "ERF_SurfaceLayer.H" @@ -207,7 +208,6 @@ Time_Avg_Vel_atCC (const amrex::Real& dt, * @param[in] dom_lo low bound of domain * @param[in] dom_hi high bound of domain * @param[in] F1 drift relaxation parameter - * @param[in] F2 Laplacian relaxation parameter * @param[in] bx_xlo box for low x relaxation * @param[in] bx_xhi box for high x relaxation * @param[in] bx_ylo box for low y relaxation @@ -243,7 +243,15 @@ realbdy_compute_relaxation (const int& icomp, amrex::IntVect iv = bx_xlo.type(); amrex::Real ioff = (iv[0]==1) ? zero : myhalf; amrex::Real joff = (iv[1]==1) ? zero : myhalf; - amrex::ParallelFor(bx_xlo, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + + // These are defined in ERF_Constants.H + amrex::Real cond_fac = lcond; // condensation + amrex::Real sub_fac = lsub; // sublimation + + int nq = num_var; + + amrex::ParallelFor( bx_xlo, bx_xhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { // Corners with x boxes amrex::Real x = ProbLo[0] + (i + ioff) * dx[0]; @@ -255,13 +263,33 @@ realbdy_compute_relaxation (const int& icomp, amrex::Real eta_lo = (y < y_end ) ? (y_end - y) / (y_end - ProbLo[1]) : zero; amrex::Real eta_hi = (y > y_strt) ? (y - y_strt) / (ProbHi[1] - y_strt) : zero; amrex::Real eta = std::max(eta_lo,eta_hi); - amrex::Real Factor = std::max(xi*xi,eta*eta); - - amrex::Real delta = arr_xlo(i ,j ,k,n) - data_arr(i ,j ,k ,n+icomp); - amrex::Real Temp = Factor*F1*delta; - rhs_arr(i,j,k,n+icomp) += Temp; + amrex::Real Factor = std::max(xi*xi,eta*eta) * F1; + + if (icomp == RhoQ1_comp) { + // qv + amrex::Real delta = arr_xlo(i,j,k,0) - data_arr(i,j,k,icomp); + rhs_arr(i,j,k,icomp) += Factor*delta; + + // qc and heat source term due to evaporation + delta = -data_arr(i,j,k,icomp+1); // This effectively nudges to 0 + rhs_arr(i,j,k,icomp+1 ) += Factor*delta; + rhs_arr(i,j,k,RhoTheta_comp) += Factor * delta * cond_fac; + + for (int n = 2; n < nq; n++) { + delta = -data_arr(i,j,k,icomp+n); // This effectively nudges to 0 + rhs_arr(i,j,k,icomp+n) += Factor*delta; + if (nq > 3 && n == 2) { // nq > 3 guarantees that component 2 is ice + rhs_arr(i,j,k,RhoTheta_comp) += Factor * delta * sub_fac; + } + } // n + } else { + for (int n = 0; n < nq; n++) { + amrex::Real delta = arr_xlo(i,j,k,n) - data_arr(i,j,k,n+icomp); + rhs_arr(i,j,k,n+icomp) += Factor*delta; + } + } }, - bx_xhi, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { // Corners with x boxes amrex::Real x = ProbLo[0] + (i + ioff) * dx[0]; @@ -273,36 +301,97 @@ realbdy_compute_relaxation (const int& icomp, amrex::Real eta_lo = (y < y_end ) ? (y_end - y) / (y_end - ProbLo[1]) : zero; amrex::Real eta_hi = (y > y_strt) ? (y - y_strt) / (ProbHi[1] - y_strt) : zero; amrex::Real eta = std::max(eta_lo,eta_hi); - amrex::Real Factor = std::max(xi*xi,eta*eta); - - amrex::Real delta = arr_xhi(i ,j ,k,n) - data_arr(i ,j ,k ,n+icomp); - amrex::Real Temp = Factor*F1*delta; - rhs_arr(i,j,k,n+icomp) += Temp; + amrex::Real Factor = std::max(xi*xi,eta*eta) * F1; + + if (icomp == RhoQ1_comp) { + // qv + amrex::Real delta = arr_xhi(i,j,k,0) - data_arr(i,j,k,icomp); + rhs_arr(i,j,k,icomp) += Factor*delta; + + // qc and heat source term due to evaporation + delta = -data_arr(i,j,k,icomp+1); // This effectively nudges to 0 + rhs_arr(i,j,k,icomp+1 ) += Factor*delta; + rhs_arr(i,j,k,RhoTheta_comp) += Factor * delta * cond_fac; + + for (int n = 2; n < nq; n++) { + delta = -data_arr(i,j,k,icomp+n); // This effectively nudges to 0 + rhs_arr(i,j,k,icomp+n) += Factor*delta; + if (nq > 3 && n == 2) { // nq > 3 guarantees that component 2 is ice + rhs_arr(i,j,k,RhoTheta_comp) += Factor * delta * sub_fac; + } + } // n + } else { + for (int n = 0; n < nq; n++) { + amrex::Real delta = arr_xhi(i,j,k,n) - data_arr(i,j,k,n+icomp); + rhs_arr(i,j,k,n+icomp) += Factor*delta; + } + } }); - amrex::ParallelFor(bx_ylo, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + amrex::ParallelFor( bx_ylo, bx_yhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { // No corners for y boxes amrex::Real y = ProbLo[1] + (j + joff) * dx[1]; amrex::Real y_end = ProbLo[1] + width * dx[1]; amrex::Real eta = (y_end - y) / (y_end - ProbLo[1]); - amrex::Real Factor = eta*eta; - - amrex::Real delta = arr_ylo(i ,j ,k,n) - data_arr(i ,j ,k ,n+icomp); - amrex::Real Temp = Factor*F1*delta; - rhs_arr(i,j,k,n+icomp) += Temp; + amrex::Real Factor = eta*eta * F1; + + if (icomp == RhoQ1_comp) { + // qv + amrex::Real delta = arr_ylo(i,j,k,0) - data_arr(i,j,k,icomp); + rhs_arr(i,j,k,icomp) += Factor*delta; + + // qc and heat source term due to evaporation + delta = -data_arr(i,j,k,icomp+1); // This effectively nudges to 0 + rhs_arr(i,j,k,icomp+1 ) += Factor*delta; + rhs_arr(i,j,k,RhoTheta_comp) += Factor * delta * cond_fac; + + for (int n = 2; n < nq; n++) { + delta = -data_arr(i,j,k,icomp+n); // This effectively nudges to 0 + rhs_arr(i,j,k,icomp+n) += Factor*delta; + if (nq > 3 && n == 2) { // nq > 3 guarantees that component 2 is ice + rhs_arr(i,j,k,RhoTheta_comp) += Factor * delta * sub_fac; + } + } // n + } else { + for (int n = 0; n < nq; n++) { + amrex::Real delta = arr_ylo(i,j,k,n) - data_arr(i,j,k,n+icomp); + rhs_arr(i,j,k,n+icomp) += Factor*delta; + } + } }, - bx_yhi, num_var, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { // No corners for y boxes amrex::Real y = ProbLo[1] + (j + joff) * dx[1]; amrex::Real y_strt = ProbHi[1] - width * dx[1]; amrex::Real eta = (y - y_strt) / (ProbHi[1] - y_strt); - amrex::Real Factor = eta*eta; - - amrex::Real delta = arr_yhi(i ,j ,k,n) - data_arr(i ,j ,k ,n+icomp); - amrex::Real Temp = Factor*F1*delta; - rhs_arr(i,j,k,n+icomp) += Temp; + amrex::Real Factor = eta*eta * F1; + + if (icomp == RhoQ1_comp) { + // qv + amrex::Real delta = arr_yhi(i,j,k,0) - data_arr(i,j,k,icomp); + rhs_arr(i,j,k,icomp) += Factor*delta; + + // qc and heat source term due to evaporation + delta = -data_arr(i,j,k,icomp+1); // This effectively nudges to 0 + rhs_arr(i,j,k,icomp+1 ) += Factor*delta; + rhs_arr(i,j,k,RhoTheta_comp) += Factor * delta * cond_fac; + + for (int n = 2; n < nq; n++) { + delta = -data_arr(i,j,k,icomp+n); // This effectively nudges to 0 + rhs_arr(i,j,k,icomp+n) += Factor*delta; + if (nq > 3 && n == 2) { // nq > 3 guarantees that component 2 is ice + rhs_arr(i,j,k,RhoTheta_comp) += Factor * delta * sub_fac; + } + } // n + } else { + for (int n = 0; n < nq; n++) { + amrex::Real delta = arr_yhi(i,j,k,n) - data_arr(i,j,k,n+icomp); + rhs_arr(i,j,k,n+icomp) += Factor*delta; + } + } }); }