diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index cc0e652ef..6f4bae92e 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -17,10 +17,6 @@ function(build_erf_lib erf_lib_name) target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_MOISTURE) - if(ERF_ENABLE_EXPLICIT_MOST_STRESS) - target_compile_definitions(${erf_lib_name} PUBLIC ERF_EXPLICIT_MOST_STRESS) - endif() - if(ERF_ENABLE_MULTIBLOCK) target_sources(${erf_lib_name} PRIVATE ${SRC_DIR}/MultiBlock/MultiBlockContainer.cpp) diff --git a/CMakeLists.txt b/CMakeLists.txt index b5195b573..377e30b96 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,8 +19,6 @@ option(ERF_ENABLE_PARTICLES "Enable Lagrangian particles" OFF) option(ERF_ENABLE_FCOMPARE "Enable building fcompare when not testing" OFF) set(ERF_PRECISION "DOUBLE" CACHE STRING "Floating point precision SINGLE or DOUBLE") -option(ERF_ENABLE_EXPLICIT_MOST_STRESS "Explicitly include MOST stresses in momentum diffusion" OFF) - option(ERF_ENABLE_MOISTURE "Enable Full Moisture" ON) option(ERF_ENABLE_WARM_NO_PRECIP "Enable Warm Moisture" OFF) option(ERF_ENABLE_RRTMGP "Enable RTE-RRTMGP Radiation" OFF) diff --git a/Docs/sphinx_doc/MOST.rst b/Docs/sphinx_doc/MOST.rst index 13eee85a0..e644529a1 100644 --- a/Docs/sphinx_doc/MOST.rst +++ b/Docs/sphinx_doc/MOST.rst @@ -134,7 +134,14 @@ In ERF, when the MOST boundary condition is applied, velocity and temperature in .. math:: - (\rho \theta)_{i,j,-n} = \rho_{i,j,-n} \left[ \frac{(\rho\theta)_{i,j,0}}{\rho_{i,j,0}} - \left. \frac{\tau_{\theta z}}{\rho} \right|_{i,j,0} \frac{\rho_{i,j,0}}{K_{\theta,v,(i,j,0)}} n \Delta z \right] + (\rho \theta)_{i,j,-n} = \rho_{i,j,-n} \left[ \frac{(\rho\theta)_{i,j,0}}{\rho_{i,j,0}} - \left. \frac{\tau_{\theta z}}{\rho} \right|_{i,j,0} \frac{\rho_{i,j,0}}{K_{\theta,v,(i,j,0)}} n \Delta z \right]. + + The above implementation explicitly sets the ghost cells so that the local stresses in (6) are recovered. This formulation will depend upon the eddy diffusivity :math:`K_{\phi,v}` in the near-wall region. Since :math:`K_{\phi,v}` may be a function of near-wall gradients, circular dependencies may occur. An **explicit MOST** formulation has also been implemented where the stress tensors are directly populated with the values computed for :math:`\tau_{\phi z}` and the ghost cells are filled according the recommendation made in `Moeng, Journal of the Atmospheric Sciences, 1984 `_; see below. To enable the **explicit MOST** formulation, users may add the line ``erf.use_explicit_most = true``. + + .. math:: + + (\rho \theta)_{z} = \frac{(\rho \theta)_{i,j,1} - (\rho \theta)_{i,j,0}}{\Delta z} + (\rho \theta)_{i,j,-n} = (\rho \theta)_{i,j,0} - (\rho \theta)_{z} n \Delta z . Finally, it must be noted that complex terrain will modify the surface normal and tangent vectors. Consequently, the MOST implentation with terrain will require local vector rotations. While the ERF dycore accounts for terrain metrics when computing fluxes (e.g. for advection, diffusion, etc.), the impact of terrain metrics on MOST is still a work in progress. Therefore, running with terrain (``erf.use_terrain = true``) and with MOST diff --git a/Exec/ABL/GNUmakefile b/Exec/ABL/GNUmakefile index 296351cfa..c85179014 100644 --- a/Exec/ABL/GNUmakefile +++ b/Exec/ABL/GNUmakefile @@ -25,7 +25,6 @@ TEST = TRUE USE_ASSERTION = TRUE #USE_POISSON_SOLVE = TRUE -USE_EXPLICIT_MOST_STRESS = FALSE # GNU Make Bpack := ./Make.package diff --git a/Exec/Make.ERF b/Exec/Make.ERF index 328613d09..82311228b 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -170,10 +170,6 @@ ERF_LSM_MM5_DIR = $(ERF_SOURCE_DIR)/LandSurfaceModel/MM5 include $(ERF_LSM_MM5_DIR)/Make.package VPATH_LOCATIONS += $(ERF_LSM_MM5_DIR) INCLUDE_LOCATIONS += $(ERF_LSM_MM5_DIR) - -ifeq ($(USE_EXPLICIT_MOST_STRESS), TRUE) - DEFINES += -DERF_EXPLICIT_MOST_STRESS -endif ifeq ($(COMPUTE_ERROR), TRUE) DEFINES += -DERF_COMPUTE_ERROR diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H index 18147d453..1e57be63f 100644 --- a/Source/BoundaryConditions/ABLMost.H +++ b/Source/BoundaryConditions/ABLMost.H @@ -30,6 +30,7 @@ public: // Constructor explicit ABLMost (const amrex::Vector& geom, + bool& use_exp_most, amrex::Vector>& vars_old, amrex::Vector>& Theta_prim, amrex::Vector>& Qv_prim, @@ -40,7 +41,8 @@ public: amrex::Vector> lsm_flux, amrex::Real start_bdy_time = 0.0, amrex::Real bdy_time_interval = 0.0) - : m_start_bdy_time(start_bdy_time), + : m_exp_most(use_exp_most), + m_start_bdy_time(start_bdy_time), m_bdy_time_interval(bdy_time_interval), m_geom(geom), m_ma(geom,vars_old,Theta_prim,Qv_prim,z_phys_nd) @@ -220,26 +222,20 @@ public: void impose_most_bcs (const int& lev, const amrex::Vector& mfs, -#ifdef ERF_EXPLICIT_MOST_STRESS amrex::MultiFab* xzmom_flux, amrex::MultiFab* zxmom_flux, amrex::MultiFab* yzmom_flux, amrex::MultiFab* zymom_flux, amrex::MultiFab* heat_flux, -#else amrex::MultiFab* eddyDiffs, -#endif amrex::MultiFab* z_phys); template void compute_most_bcs (const int& lev, const amrex::Vector& mfs, -#ifdef ERF_EXPLICIT_MOST_STRESS - amrex::MultiFab* xzmom_flux, amrex::MultiFab* zxmom_flux, - amrex::MultiFab* yzmom_flux, amrex::MultiFab* zymom_flux, + amrex::MultiFab* xzmom_flux, amrex::MultiFab* zxmom_flux, + amrex::MultiFab* yzmom_flux, amrex::MultiFab* zymom_flux, amrex::MultiFab* heat_flux, -#else amrex::MultiFab* eddyDiffs, -#endif amrex::MultiFab* z_phys, const amrex::Real& dz_no_terrain, const FluxCalc& flux_comp); @@ -319,6 +315,7 @@ public: private: bool use_moisture; + bool m_exp_most = false; amrex::Real z0_const; amrex::Real surf_temp; amrex::Real surf_heating_rate{0}; diff --git a/Source/BoundaryConditions/ABLMost.cpp b/Source/BoundaryConditions/ABLMost.cpp index 1bf77fa71..90d352515 100644 --- a/Source/BoundaryConditions/ABLMost.cpp +++ b/Source/BoundaryConditions/ABLMost.cpp @@ -121,48 +121,36 @@ ABLMost::compute_fluxes (const int& lev, void ABLMost::impose_most_bcs (const int& lev, const Vector& mfs, -#ifdef ERF_EXPLICIT_MOST_STRESS MultiFab* xzmom_flux, MultiFab* zxmom_flux, MultiFab* yzmom_flux, MultiFab* zymom_flux, MultiFab* heat_flux, -#else MultiFab* eddyDiffs, -#endif MultiFab* z_phys) { const int klo = 0; if (flux_type == FluxCalcType::MOENG) { moeng_flux flux_comp(klo); compute_most_bcs(lev, mfs, -#ifdef ERF_EXPLICIT_MOST_STRESS xzmom_flux, xzmom_flux, yzmom_flux, yzmom_flux, heat_flux, -#else eddyDiffs, -#endif z_phys, m_geom[lev].CellSize(2), flux_comp); } else if (flux_type == FluxCalcType::DONELAN) { donelan_flux flux_comp(klo); compute_most_bcs(lev, mfs, -#ifdef ERF_EXPLICIT_MOST_STRESS xzmom_flux, xzmom_flux, yzmom_flux, yzmom_flux, heat_flux, -#else eddyDiffs, -#endif z_phys, m_geom[lev].CellSize(2), flux_comp); } else { custom_flux flux_comp(klo); compute_most_bcs(lev, mfs, -#ifdef ERF_EXPLICIT_MOST_STRESS xzmom_flux, xzmom_flux, yzmom_flux, yzmom_flux, heat_flux, -#else eddyDiffs, -#endif z_phys, m_geom[lev].CellSize(2), flux_comp); } } @@ -180,13 +168,10 @@ template void ABLMost::compute_most_bcs (const int& lev, const Vector& mfs, -#ifdef ERF_EXPLICIT_MOST_STRESS MultiFab* xzmom_flux, MultiFab* zxmom_flux, MultiFab* yzmom_flux, MultiFab* zymom_flux, MultiFab* heat_flux, -#else MultiFab* eddyDiffs, -#endif MultiFab* z_phys, const Real& dz_no_terrain, const FluxCalc& flux_comp) @@ -199,36 +184,38 @@ ABLMost::compute_most_bcs (const int& lev, // Valid CC box Box vbx = mfi.validbox(); vbx.makeSlab(2,klo-1); -#ifdef ERF_EXPLICIT_MOST_STRESS Box vbxx = surroundingNodes(vbx,0); Box vbxy = surroundingNodes(vbx,1); -#endif + + // Expose for GPU + bool exp_most = m_exp_most; // Get field arrays const auto cons_arr = mfs[Vars::cons]->array(mfi); const auto velx_arr = mfs[Vars::xvel]->array(mfi); const auto vely_arr = mfs[Vars::yvel]->array(mfi); -#ifdef ERF_EXPLICIT_MOST_STRESS - auto t13_arr = xzmom_flux->array(mfi); - auto t23_arr = yzmom_flux->array(mfi); - auto t31_arr = (zxmom_flux) ? zxmom_flux->array(mfi) : Array4{}; - auto t32_arr = (zymom_flux) ? zymom_flux->array(mfi) : Array4{}; - auto hfx_arr = heat_flux->array(mfi); -#else + + auto t13_arr = (m_exp_most) ? xzmom_flux->array(mfi) : Array4{}; + auto t23_arr = (m_exp_most) ? yzmom_flux->array(mfi) : Array4{}; + auto t31_arr = (zxmom_flux && m_exp_most) ? zxmom_flux->array(mfi) : Array4{}; + auto t32_arr = (zymom_flux && m_exp_most) ? zymom_flux->array(mfi) : Array4{}; + auto hfx_arr = (m_exp_most) ? heat_flux->array(mfi) : Array4{}; + const auto eta_arr = eddyDiffs->array(mfi); -#endif + const auto zphys_arr = (z_phys) ? z_phys->const_array(mfi) : Array4{}; // Get average arrays const auto *const u_mean = m_ma.get_average(lev,0); const auto *const v_mean = m_ma.get_average(lev,1); const auto *const t_mean = m_ma.get_average(lev,2); - // const auto *const q_mean = m_ma.get_average(lev,3); + const auto *const q_mean = m_ma.get_average(lev,3); const auto *const u_mag_mean = m_ma.get_average(lev,4); const auto um_arr = u_mean->array(mfi); const auto vm_arr = v_mean->array(mfi); const auto tm_arr = t_mean->array(mfi); + const auto qm_arr = q_mean->array(mfi); const auto umm_arr = u_mag_mean->array(mfi); // Get derived arrays @@ -254,21 +241,12 @@ ABLMost::compute_most_bcs (const int& lev, int n = RhoTheta_comp; ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo) - zphys_arr(i,j,klo-1) ) : dz_no_terrain; - -#ifdef ERF_EXPLICIT_MOST_STRESS - Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo) ) : dz_no_terrain; - Real Tflux = flux_comp.compute_t_flux(i, j, k, n, icomp, dz, dz1, - cons_arr, velx_arr, vely_arr, - umm_arr, tm_arr, u_star_arr, t_star_arr, t_surf_arr, - dest_arr); -#else - Real Tflux = flux_comp.compute_t_flux(i, j, k, n, icomp, dz, + Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo ) - zphys_arr(i,j,klo-1) ) : dz_no_terrain; + Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo ) ) : dz_no_terrain; + Real Tflux = flux_comp.compute_t_flux(i, j, k, n, icomp, dz, dz1, exp_most, eta_arr, cons_arr, velx_arr, vely_arr, - eta_arr, umm_arr, tm_arr, u_star_arr, t_star_arr, t_surf_arr, dest_arr); -#endif // TODO: make sure not to double-count surface heat flux if using a LSM @@ -276,11 +254,9 @@ ABLMost::compute_most_bcs (const int& lev, if (is_land && lsm_flux_arr && vbx.contains(i,j,k)) { lsm_flux_arr(i,j,klo) = Tflux; } -#ifdef ERF_EXPLICIT_MOST_STRESS - else if ((k == klo-1) && vbx.contains(i,j,k)) { + else if ((k == klo-1) && vbx.contains(i,j,k) && exp_most) { hfx_arr(i,j,klo-1) = Tflux; } -#endif }); // TODO: Generalize MOST q flux with MOENG & DONELAN flux types @@ -288,20 +264,12 @@ ABLMost::compute_most_bcs (const int& lev, n = RhoQ1_comp; ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo) - zphys_arr(i,j,klo-1) ) : dz_no_terrain; -#ifdef ERF_EXPLICIT_MOST_STRESS - Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo) ) : dz_no_terrain; - Real Qflux = flux_comp.compute_q_flux(i, j, k, n, icomp, dz, dz1, + Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo ) - zphys_arr(i,j,klo-1) ) : dz_no_terrain; + Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo ) ) : dz_no_terrain; + Real Qflux = flux_comp.compute_q_flux(i, j, k, n, icomp, dz, dz1, exp_most, eta_arr, cons_arr, velx_arr, vely_arr, - umm_arr, tm_arr, u_star_arr, q_star_arr, t_surf_arr, + umm_arr, qm_arr, u_star_arr, q_star_arr, t_surf_arr, dest_arr); -#else - Real Qflux = flux_comp.compute_q_flux(i, j, k, n, icomp, dz, - cons_arr, velx_arr, vely_arr, - eta_arr, - umm_arr, tm_arr, u_star_arr, q_star_arr, t_surf_arr, - dest_arr); -#endif amrex::ignore_unused(Qflux); }); } @@ -313,24 +281,16 @@ ABLMost::compute_most_bcs (const int& lev, ParallelFor(xb2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo) - zphys_arr(i,j,klo-1) ) : dz_no_terrain; -#ifdef ERF_EXPLICIT_MOST_STRESS - Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo) ) : dz_no_terrain; - Real stressx = flux_comp.compute_u_flux(i, j, k, icomp, dz, dz1, + Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo ) - zphys_arr(i,j,klo-1) ) : dz_no_terrain; + Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo ) ) : dz_no_terrain; + Real stressx = flux_comp.compute_u_flux(i, j, k, icomp, dz, dz1, exp_most, eta_arr, cons_arr, velx_arr, vely_arr, umm_arr, um_arr, u_star_arr, dest_arr); - if ((k == klo-1) && vbxx.contains(i,j,k)) { + if ((k == klo-1) && vbxx.contains(i,j,k) && exp_most) { t13_arr(i,j,klo) = -stressx; if (t31_arr) t31_arr(i,j,klo) = -stressx; } -#else - Real stressx = flux_comp.compute_u_flux(i, j, k, icomp, dz, - cons_arr, velx_arr, vely_arr, - eta_arr, - umm_arr, um_arr, u_star_arr, - dest_arr); -#endif amrex::ignore_unused(stressx); }); @@ -341,24 +301,16 @@ ABLMost::compute_most_bcs (const int& lev, ParallelFor(yb2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo) - zphys_arr(i,j,klo-1) ) : dz_no_terrain; -#ifdef ERF_EXPLICIT_MOST_STRESS - Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo) ) : dz_no_terrain; - Real stressy = flux_comp.compute_v_flux(i, j, k, icomp, dz, dz1, + Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo ) - zphys_arr(i,j,klo-1) ) : dz_no_terrain; + Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo ) ) : dz_no_terrain; + Real stressy = flux_comp.compute_v_flux(i, j, k, icomp, dz, dz1, exp_most, eta_arr, cons_arr, velx_arr, vely_arr, umm_arr, vm_arr, u_star_arr, dest_arr); - if ((k == klo-1) && vbxy.contains(i,j,k)) { + if ((k == klo-1) && vbxy.contains(i,j,k) && exp_most) { t23_arr(i,j,klo) = -stressy; if (t32_arr) t32_arr(i,j,klo) = -stressy; } -#else - Real stressy = flux_comp.compute_v_flux(i, j, k, icomp, dz, - cons_arr, velx_arr, vely_arr, - eta_arr, - umm_arr, vm_arr, u_star_arr, - dest_arr); -#endif amrex::ignore_unused(stressy); }); } diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index d0b786c31..a63253297 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -383,13 +383,10 @@ ERF::FillIntermediatePatch (int lev, Real time, // MOST boundary conditions if (!(cons_only && ncomp_cons == 1) && m_most && allow_most_bcs) { m_most->impose_most_bcs(lev,mfs_vel, -#ifdef ERF_EXPLICIT_MOST_STRESS Tau13_lev[lev].get(), Tau31_lev[lev].get(), Tau23_lev[lev].get(), Tau32_lev[lev].get(), SFS_hfx3_lev[lev].get(), -#else eddyDiffs_lev[lev].get(), -#endif z_phys_nd[lev].get()); } diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H index 874f9ffb8..db2c40773 100644 --- a/Source/BoundaryConditions/MOSTStress.H +++ b/Source/BoundaryConditions/MOSTStress.H @@ -600,17 +600,14 @@ struct moeng_flux const int& n, const int& icomp, const amrex::Real& dz, -#ifdef ERF_EXPLICIT_MOST_STRESS const amrex::Real& dz1, -#endif + const bool& exp_most, + const amrex::Array4& eta_arr, const amrex::Array4& cons_arr, const amrex::Array4& /*velx_arr*/, const amrex::Array4& /*vely_arr*/, -#ifndef ERF_EXPLICIT_MOST_STRESS - const amrex::Array4& eta_arr, -#endif const amrex::Array4& /*umm_arr*/, - const amrex::Array4& /*tm_arr*/, + const amrex::Array4& /*qm_arr*/, const amrex::Array4& /*u_star_arr*/, const amrex::Array4& /*q_star_arr*/, const amrex::Array4& /*t_surf_arr*/, @@ -631,9 +628,8 @@ struct moeng_flux amrex::Real moflux = 0.0; amrex::Real deltaz = dz * (zlo - k); -#ifdef ERF_EXPLICIT_MOST_STRESS - amrex::Abort("Explicit MOST stress not implemented for Moeng moisture flux"); -#else + if (exp_most) amrex::Abort("Explicit MOST stress not implemented for Moeng moisture flux"); + int ie, je; ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; @@ -642,7 +638,6 @@ struct moeng_flux eta = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s] eta = amrex::max(eta,eta_eps); dest_arr(i,j,k,icomp+n) = rho*(qv - moflux*rho/eta*deltaz); -#endif return moflux; } @@ -656,15 +651,12 @@ struct moeng_flux const int& n, const int& icomp, const amrex::Real& dz, -#ifdef ERF_EXPLICIT_MOST_STRESS const amrex::Real& dz1, -#endif + const bool& exp_most, + const amrex::Array4& eta_arr, const amrex::Array4& cons_arr, const amrex::Array4& velx_arr, const amrex::Array4& vely_arr, -#ifndef ERF_EXPLICIT_MOST_STRESS - const amrex::Array4& eta_arr, -#endif const amrex::Array4& umm_arr, const amrex::Array4& tm_arr, const amrex::Array4& u_star_arr, @@ -708,23 +700,23 @@ struct moeng_flux -tstar*ustar*(num1+num2)/((theta_mean-theta_surf)*wsp_mean) : 0.0; amrex::Real deltaz = dz * (zlo - k); -#ifdef ERF_EXPLICIT_MOST_STRESS - // surface gradient equal to gradient at first zface - amrex::Real thetagrad = (cons_arr(ic,jc,zlo+1,RhoTheta_comp) - cons_arr(ic,jc,zlo,RhoTheta_comp)) / (0.5*(dz+dz1)); - dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - thetagrad * deltaz; -#else - int ie, je; - ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; - je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; - ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; - je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; - eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s] - eta = amrex::max(eta,eta_eps); - // Note: Kh = eta/rho - // hfx = -Kh dT/dz ==> +ve hfx corresponds to heating from the surface - // Extrapolate from klo to ghost cell a distance of -deltaz; negative signs cancel - dest_arr(i,j,k,icomp+n) = rho*(theta + moflux*rho/eta*deltaz); -#endif + if (exp_most) { + // surface gradient equal to gradient at first zface + amrex::Real thetagrad = (cons_arr(ic,jc,zlo+1,RhoTheta_comp) - cons_arr(ic,jc,zlo,RhoTheta_comp)) / (0.5*(dz+dz1)); + dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - thetagrad * deltaz; + } else { + int ie, je; + ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; + je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; + ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; + je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; + eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s] + eta = amrex::max(eta,eta_eps); + // Note: Kh = eta/rho + // hfx = -Kh dT/dz ==> +ve hfx corresponds to heating from the surface + // Extrapolate from klo to ghost cell a distance of -deltaz; negative signs cancel + dest_arr(i,j,k,icomp+n) = rho*(theta + moflux*rho/eta*deltaz); + } return moflux; } @@ -737,15 +729,12 @@ struct moeng_flux const int& k, const int& icomp, const amrex::Real& dz, -#ifdef ERF_EXPLICIT_MOST_STRESS const amrex::Real& dz1, -#endif + const bool& exp_most, + const amrex::Array4& eta_arr, const amrex::Array4& cons_arr, const amrex::Array4& velx_arr, const amrex::Array4& vely_arr, -#ifndef ERF_EXPLICIT_MOST_STRESS - const amrex::Array4& eta_arr, -#endif const amrex::Array4& umm_arr, const amrex::Array4& um_arr, const amrex::Array4& u_star_arr, @@ -785,21 +774,21 @@ struct moeng_flux amrex::Real stressx = rho*ustar*ustar * (num1+num2)/(wsp_mean*wsp_mean); amrex::Real deltaz = dz * (zlo - k); -#ifdef ERF_EXPLICIT_MOST_STRESS - // surface gradient equal to gradient at first zface - amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx) / (0.5*(dz+dz1)); - dest_arr(i,j,k,icomp) = velx - ugrad * deltaz; -#else - int ie, je; - ie = i < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i; - je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; - ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; - je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; - eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v) - + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) ); - eta = amrex::max(eta,eta_eps); - dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressx/eta*deltaz; -#endif + if (exp_most) { + // surface gradient equal to gradient at first zface + amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx) / (0.5*(dz+dz1)); + dest_arr(i,j,k,icomp) = velx - ugrad * deltaz; + } else { + int ie, je; + ie = i < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i; + je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; + ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; + je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; + eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v) + + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) ); + eta = amrex::max(eta,eta_eps); + dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressx/eta*deltaz; + } return stressx; } @@ -812,15 +801,12 @@ struct moeng_flux const int& k, const int& icomp, const amrex::Real& dz, -#ifdef ERF_EXPLICIT_MOST_STRESS const amrex::Real& dz1, -#endif + const bool& exp_most, + const amrex::Array4& eta_arr, const amrex::Array4& cons_arr, const amrex::Array4& velx_arr, const amrex::Array4& vely_arr, -#ifndef ERF_EXPLICIT_MOST_STRESS - const amrex::Array4& eta_arr, -#endif const amrex::Array4& umm_arr, const amrex::Array4& vm_arr, const amrex::Array4& u_star_arr, @@ -860,21 +846,21 @@ struct moeng_flux amrex::Real stressy = rho*ustar*ustar * (num1+num2)/(wsp_mean*wsp_mean); amrex::Real deltaz = dz * (zlo - k); -#ifdef ERF_EXPLICIT_MOST_STRESS - // surface gradient equal to gradient at first zface - amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely) / (0.5*(dz+dz1)); - dest_arr(i,j,k,icomp) = vely - vgrad * deltaz; -#else - int ie, je; - ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; - je = j < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j; - ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; - je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; - eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v) - + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) ); - eta = amrex::max(eta,eta_eps); - dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressy/eta*deltaz; -#endif + if (exp_most) { + // surface gradient equal to gradient at first zface + amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely) / (0.5*(dz+dz1)); + dest_arr(i,j,k,icomp) = vely - vgrad * deltaz; + } else { + int ie, je; + ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; + je = j < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j; + ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; + je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; + eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v) + + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) ); + eta = amrex::max(eta,eta_eps); + dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressy/eta*deltaz; + } return stressy; } @@ -904,17 +890,14 @@ struct donelan_flux const int& n, const int& icomp, const amrex::Real& dz, -#ifdef ERF_EXPLICIT_MOST_STRESS const amrex::Real& dz1, -#endif + const bool& exp_most, + const amrex::Array4& eta_arr, const amrex::Array4& cons_arr, const amrex::Array4& /*velx_arr*/, const amrex::Array4& /*vely_arr*/, -#ifndef ERF_EXPLICIT_MOST_STRESS - const amrex::Array4& eta_arr, -#endif const amrex::Array4& /*umm_arr*/, - const amrex::Array4& /*tm_arr*/, + const amrex::Array4& /*qm_arr*/, const amrex::Array4& /*u_star_arr*/, const amrex::Array4& /*q_star_arr*/, const amrex::Array4& /*t_surf_arr*/, @@ -935,9 +918,8 @@ struct donelan_flux amrex::Real moflux = 0.0; amrex::Real deltaz = dz * (zlo - k); -#ifdef ERF_EXPLICIT_MOST_STRESS - amrex::Abort("Explicit MOST stress not implemented for Donelan moisture flux"); -#else + if (exp_most) amrex::Abort("Explicit MOST stress not implemented for Donelan moisture flux"); + int ie, je; ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; @@ -946,7 +928,6 @@ struct donelan_flux eta = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s] eta = amrex::max(eta,eta_eps); dest_arr(i,j,k,icomp+n) = rho*(qv - moflux*rho/eta*deltaz); -#endif return moflux; } @@ -960,15 +941,12 @@ struct donelan_flux const int& n, const int& icomp, const amrex::Real& dz, -#ifdef ERF_EXPLICIT_MOST_STRESS const amrex::Real& dz1, -#endif + const bool& exp_most, + const amrex::Array4& eta_arr, const amrex::Array4& cons_arr, const amrex::Array4& /*velx_arr*/, const amrex::Array4& /*vely_arr*/, -#ifndef ERF_EXPLICIT_MOST_STRESS - const amrex::Array4& eta_arr, -#endif const amrex::Array4& umm_arr, const amrex::Array4& tm_arr, const amrex::Array4& /*u_star_arr*/, @@ -994,23 +972,23 @@ struct donelan_flux amrex::Real moflux = Cd * wsp_mean * (theta_surf - theta_mean); amrex::Real deltaz = dz * (zlo - k); -#ifdef ERF_EXPLICIT_MOST_STRESS - // surface gradient equal to gradient at first zface - amrex::Real thetagrad = (cons_arr(ic,jc,zlo+1,RhoTheta_comp) - cons_arr(ic,jc,zlo,RhoTheta_comp)) / (0.5*(dz+dz1)); - dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - thetagrad * deltaz; -#else - int ie, je; - ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; - je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; - ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; - je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; - eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s] - eta = amrex::max(eta,eta_eps); - // Note: Kh = eta/rho - // hfx = -Kh dT/dz ==> +ve hfx corresponds to heating from the surface - // Extrapolate from klo to ghost cell a distance of -deltaz; negative signs cancel - dest_arr(i,j,k,icomp+n) = rho*(theta + moflux*rho/eta*deltaz); -#endif + if (exp_most) { + // surface gradient equal to gradient at first zface + amrex::Real thetagrad = (cons_arr(ic,jc,zlo+1,RhoTheta_comp) - cons_arr(ic,jc,zlo,RhoTheta_comp)) / (0.5*(dz+dz1)); + dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - thetagrad * deltaz; + } else { + int ie, je; + ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; + je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; + ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; + je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; + eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s] + eta = amrex::max(eta,eta_eps); + // Note: Kh = eta/rho + // hfx = -Kh dT/dz ==> +ve hfx corresponds to heating from the surface + // Extrapolate from klo to ghost cell a distance of -deltaz; negative signs cancel + dest_arr(i,j,k,icomp+n) = rho*(theta + moflux*rho/eta*deltaz); + } return moflux; } @@ -1023,15 +1001,12 @@ struct donelan_flux const int& k, const int& icomp, const amrex::Real& dz, -#ifdef ERF_EXPLICIT_MOST_STRESS const amrex::Real& dz1, -#endif + const bool& exp_most, + const amrex::Array4& eta_arr, const amrex::Array4& cons_arr, const amrex::Array4& velx_arr, const amrex::Array4& vely_arr, -#ifndef ERF_EXPLICIT_MOST_STRESS - const amrex::Array4& eta_arr, -#endif const amrex::Array4& umm_arr, const amrex::Array4& /*um_arr*/, const amrex::Array4& /*u_star_arr*/, @@ -1071,21 +1046,21 @@ struct donelan_flux amrex::Real stressx = rho * Cd * velx * wsp; amrex::Real deltaz = dz * (zlo - k); -#ifdef ERF_EXPLICIT_MOST_STRESS - // surface gradient equal to gradient at first zface - amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx) / (0.5*(dz+dz1)); - dest_arr(i,j,k,icomp) = velx - ugrad * deltaz; -#else - int ie, je; - ie = i < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i; - je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; - ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; - je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; - eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v) - + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) ); - eta = amrex::max(eta,eta_eps); - dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressx/eta*deltaz; -#endif + if (exp_most) { + // surface gradient equal to gradient at first zface + amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx) / (0.5*(dz+dz1)); + dest_arr(i,j,k,icomp) = velx - ugrad * deltaz; + } else { + int ie, je; + ie = i < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i; + je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; + ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; + je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; + eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v) + + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) ); + eta = amrex::max(eta,eta_eps); + dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressx/eta*deltaz; + } return stressx; } @@ -1098,15 +1073,12 @@ struct donelan_flux const int& k, const int& icomp, const amrex::Real& dz, -#ifdef ERF_EXPLICIT_MOST_STRESS const amrex::Real& dz1, -#endif + const bool& exp_most, + const amrex::Array4& eta_arr, const amrex::Array4& cons_arr, const amrex::Array4& velx_arr, const amrex::Array4& vely_arr, -#ifndef ERF_EXPLICIT_MOST_STRESS - const amrex::Array4& eta_arr, -#endif const amrex::Array4& umm_arr, const amrex::Array4& /*vm_arr*/, const amrex::Array4& /*u_star_arr*/, @@ -1146,21 +1118,21 @@ struct donelan_flux amrex::Real stressy = rho * Cd * vely * wsp; amrex::Real deltaz = dz * (zlo - k); -#ifdef ERF_EXPLICIT_MOST_STRESS - // surface gradient equal to gradient at first zface - amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely) / (0.5*(dz+dz1)); - dest_arr(i,j,k,icomp) = vely - vgrad * deltaz; -#else - int ie, je; - ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; - je = j < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j; - ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; - je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; - eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v) - + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) ); - eta = amrex::max(eta,eta_eps); - dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressy/eta*deltaz; -#endif + if (exp_most) { + // surface gradient equal to gradient at first zface + amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely) / (0.5*(dz+dz1)); + dest_arr(i,j,k,icomp) = vely - vgrad * deltaz; + } else { + int ie, je; + ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; + je = j < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j; + ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; + je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; + eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v) + + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) ); + eta = amrex::max(eta,eta_eps); + dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressy/eta*deltaz; + } return stressy; } @@ -1189,17 +1161,14 @@ struct custom_flux const int& n, const int& icomp, const amrex::Real& dz, -#ifdef ERF_EXPLICIT_MOST_STRESS const amrex::Real& dz1, -#endif + const bool& exp_most, + const amrex::Array4& eta_arr, const amrex::Array4& cons_arr, const amrex::Array4& /*velx_arr*/, const amrex::Array4& /*vely_arr*/, -#ifndef ERF_EXPLICIT_MOST_STRESS - const amrex::Array4& eta_arr, -#endif const amrex::Array4& /*umm_arr*/, - const amrex::Array4& /*tm_arr*/, + const amrex::Array4& /*qm_arr*/, const amrex::Array4& /*u_star_arr*/, const amrex::Array4& q_star_arr, const amrex::Array4& /*t_surf_arr*/, @@ -1220,9 +1189,8 @@ struct custom_flux amrex::Real moflux = (std::abs(qstar) > eps) ? -rho * qstar : 0.0; amrex::Real deltaz = dz * (zlo - k); -#ifdef ERF_EXPLICIT_MOST_STRESS - amrex::Abort("Explicit MOST stress not implemented for custom moisture flux"); -#else + if (exp_most) amrex::Abort("Explicit MOST stress not implemented for custom moisture flux"); + int ie, je; ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; @@ -1231,7 +1199,6 @@ struct custom_flux eta = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s] eta = amrex::max(eta,eta_eps); dest_arr(i,j,k,icomp+n) = rho*(qv - moflux*rho/eta*deltaz); -#endif return moflux; } @@ -1245,15 +1212,12 @@ struct custom_flux const int& n, const int& icomp, const amrex::Real& dz, -#ifdef ERF_EXPLICIT_MOST_STRESS const amrex::Real& dz1, -#endif + const bool& exp_most, + const amrex::Array4& eta_arr, const amrex::Array4& cons_arr, const amrex::Array4& /*velx_arr*/, const amrex::Array4& /*vely_arr*/, -#ifndef ERF_EXPLICIT_MOST_STRESS - const amrex::Array4& eta_arr, -#endif const amrex::Array4& /*umm_arr*/, const amrex::Array4& /*tm_arr*/, const amrex::Array4& u_star_arr, @@ -1277,24 +1241,24 @@ struct custom_flux amrex::Real moflux = (std::abs(tstar) > eps) ? -tstar*ustar : 0.0; amrex::Real deltaz = dz * (zlo - k); -#ifdef ERF_EXPLICIT_MOST_STRESS - // surface gradient equal to gradient at first zface - amrex::Real thetagrad = ( cons_arr(ic,jc,zlo+1,RhoTheta_comp)/cons_arr(ic,jc,zlo+1,Rho_comp) - - cons_arr(ic,jc,zlo ,RhoTheta_comp)/cons_arr(ic,jc,zlo ,Rho_comp)) / (0.5*(dz+dz1)); - dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - rho*thetagrad * deltaz; -#else - int ie, je; - ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; - je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; - ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; - je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; - eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s] - eta = amrex::max(eta,eta_eps); - // Note: Kh = eta/rho - // hfx = -Kh dT/dz ==> +ve hfx corresponds to heating from the surface - // Extrapolate from klo to ghost cell a distance of -deltaz; negative signs cancel - dest_arr(i,j,k,icomp+n) = rho*(theta + moflux*rho/eta*deltaz); -#endif + if (exp_most) { + // surface gradient equal to gradient at first zface + amrex::Real thetagrad = ( cons_arr(ic,jc,zlo+1,RhoTheta_comp)/cons_arr(ic,jc,zlo+1,Rho_comp) + - cons_arr(ic,jc,zlo ,RhoTheta_comp)/cons_arr(ic,jc,zlo ,Rho_comp)) / (0.5*(dz+dz1)); + dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - rho*thetagrad * deltaz; + } else { + int ie, je; + ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; + je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; + ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; + je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; + eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s] + eta = amrex::max(eta,eta_eps); + // Note: Kh = eta/rho + // hfx = -Kh dT/dz ==> +ve hfx corresponds to heating from the surface + // Extrapolate from klo to ghost cell a distance of -deltaz; negative signs cancel + dest_arr(i,j,k,icomp+n) = rho*(theta + moflux*rho/eta*deltaz); + } return moflux; } @@ -1307,15 +1271,12 @@ struct custom_flux const int& k, const int& icomp, const amrex::Real& dz, -#ifdef ERF_EXPLICIT_MOST_STRESS const amrex::Real& dz1, -#endif + const bool& exp_most, + const amrex::Array4& eta_arr, const amrex::Array4& cons_arr, const amrex::Array4& velx_arr, const amrex::Array4& vely_arr, -#ifndef ERF_EXPLICIT_MOST_STRESS - const amrex::Array4& eta_arr, -#endif const amrex::Array4& /*umm_arr*/, const amrex::Array4& /*um_arr*/, const amrex::Array4& u_star_arr, @@ -1346,21 +1307,21 @@ struct custom_flux amrex::Real stressx = rho * ustar * ustar * velx / wsp; amrex::Real deltaz = dz * (zlo - k); -#ifdef ERF_EXPLICIT_MOST_STRESS - // surface gradient equal to gradient at first zface - amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx) / (0.5*(dz+dz1)); - dest_arr(i,j,k,icomp) = velx - ugrad * deltaz; -#else - int ie, je; - ie = i < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i; - je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; - ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; - je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; - eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v) - + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) ); - eta = amrex::max(eta,eta_eps); - dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressx/eta*deltaz; -#endif + if (exp_most) { + // surface gradient equal to gradient at first zface + amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx) / (0.5*(dz+dz1)); + dest_arr(i,j,k,icomp) = velx - ugrad * deltaz; + } else { + int ie, je; + ie = i < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i; + je = j < lbound(eta_arr).y ? lbound(eta_arr).y : j; + ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; + je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; + eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v) + + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) ); + eta = amrex::max(eta,eta_eps); + dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressx/eta*deltaz; + } return stressx; } @@ -1373,15 +1334,12 @@ struct custom_flux const int& k, const int& icomp, const amrex::Real& dz, -#ifdef ERF_EXPLICIT_MOST_STRESS const amrex::Real& dz1, -#endif + const bool& exp_most, + const amrex::Array4& eta_arr, const amrex::Array4& cons_arr, const amrex::Array4& velx_arr, const amrex::Array4& vely_arr, -#ifndef ERF_EXPLICIT_MOST_STRESS - const amrex::Array4& eta_arr, -#endif const amrex::Array4& /*umm_arr*/, const amrex::Array4& /*vm_arr*/, const amrex::Array4& u_star_arr, @@ -1412,21 +1370,21 @@ struct custom_flux amrex::Real stressy = rho * ustar * ustar * vely / wsp; amrex::Real deltaz = dz * (zlo - k); -#ifdef ERF_EXPLICIT_MOST_STRESS - // surface gradient equal to gradient at first zface - amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely) / (0.5*(dz+dz1)); - dest_arr(i,j,k,icomp) = vely - vgrad * deltaz; -#else - int ie, je; - ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; - je = j < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j; - ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; - je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; - eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v) - + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) ); - eta = amrex::max(eta,eta_eps); - dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressy/eta*deltaz; -#endif + if (exp_most) { + // surface gradient equal to gradient at first zface + amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely) / (0.5*(dz+dz1)); + dest_arr(i,j,k,icomp) = vely - vgrad * deltaz; + } else { + int ie, je; + ie = i < lbound(eta_arr).x ? lbound(eta_arr).x : i; + je = j < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j; + ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie; + je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je; + eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v) + + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) ); + eta = amrex::max(eta,eta_eps); + dest_arr(i,j,k,icomp) = dest_arr(i,j,zlo,icomp) - stressy/eta*deltaz; + } return stressy; } diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index d9108830a..5a092b237 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -190,6 +190,9 @@ struct SolverChoice { pp.query("rayleigh_damp_W", rayleigh_damp_W); pp.query("rayleigh_damp_T", rayleigh_damp_T); + // Flag to do explicit MOST formulation + pp.query("use_explicit_most",use_explicit_most); + // Which external forcings? static std::string abl_driver_type_string = "None"; pp.query("abl_driver_type",abl_driver_type_string); @@ -434,6 +437,9 @@ struct SolverChoice { bool custom_geostrophic_profile = false; bool custom_forcing_prim_vars = false; + // User specified MOST BC type + bool use_explicit_most = false; + // User wishes to output time averaged velocity fields bool time_avg_vel = false; diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp index 403fc97ae..44c90a505 100644 --- a/Source/Diffusion/ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp @@ -48,7 +48,8 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, const Geometry& geom, const MultiFab& mapfac_u, const MultiFab& mapfac_v, const std::unique_ptr& z_phys_nd, - const TurbChoice& turbChoice, const Real const_grav, std::unique_ptr& most) + const TurbChoice& turbChoice, const Real const_grav, + std::unique_ptr& most, const bool& exp_most) { const GpuArray cellSizeInv = geom.InvCellSizeArray(); const Box& domain = geom.Domain(); @@ -88,7 +89,7 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real SmnSmn = ComputeSmnSmn(i,j,k,tau11,tau22,tau33,tau12,tau13,tau23,klo,use_most); + Real SmnSmn = ComputeSmnSmn(i,j,k,tau11,tau22,tau33,tau12,tau13,tau23,klo,use_most,exp_most); Real dxInv = cellSizeInv[0]; Real dyInv = cellSizeInv[1]; Real dzInv = cellSizeInv[2]; @@ -152,17 +153,17 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, Real eps = std::numeric_limits::epsilon(); Real dtheta_dz; if (use_most && k==klo) { -#ifdef ERF_EXPLICIT_MOST_STRESS - dtheta_dz = ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) - - cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) )*dzInv; -#else - dtheta_dz = 0.5 * (-3 * cell_data(i,j,k ,RhoTheta_comp) - / cell_data(i,j,k ,Rho_comp) - + 4 * cell_data(i,j,k+1,RhoTheta_comp) - / cell_data(i,j,k+1,Rho_comp) - - cell_data(i,j,k+2,RhoTheta_comp) - / cell_data(i,j,k+2,Rho_comp) ) * dzInv; -#endif + if (exp_most) { + dtheta_dz = ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) + - cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) )*dzInv; + } else { + dtheta_dz = 0.5 * (-3 * cell_data(i,j,k ,RhoTheta_comp) + / cell_data(i,j,k ,Rho_comp) + + 4 * cell_data(i,j,k+1,RhoTheta_comp) + / cell_data(i,j,k+1,Rho_comp) + - cell_data(i,j,k+2,RhoTheta_comp) + / cell_data(i,j,k+2,Rho_comp) ) * dzInv; + } } else { dtheta_dz = 0.5 * ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) - cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dzInv; @@ -413,6 +414,7 @@ void ComputeTurbulentViscosity (const MultiFab& xvel , const MultiFab& yvel , const std::unique_ptr& z_phys_nd, const TurbChoice& turbChoice, const Real const_grav, std::unique_ptr& most, + const bool& exp_most, int level, const BCRec* bc_ptr, bool vert_only) @@ -445,7 +447,8 @@ void ComputeTurbulentViscosity (const MultiFab& xvel , const MultiFab& yvel , cons_in, eddyViscosity, Hfx1, Hfx2, Hfx3, Diss, geom, mapfac_u, mapfac_v, - z_phys_nd, turbChoice, const_grav, most); + z_phys_nd, turbChoice, const_grav, + most, exp_most); } if (turbChoice.pbl_type != PBLType::None) { diff --git a/Source/Diffusion/Diffusion.H b/Source/Diffusion/Diffusion.H index ce4aefc81..e62d86e6b 100644 --- a/Source/Diffusion/Diffusion.H +++ b/Source/Diffusion/Diffusion.H @@ -42,6 +42,7 @@ void DiffusionSrcForMom_T (const amrex::Box& bxx, const amrex::Box& bxy, const a void DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain, int start_comp, int num_comp, + const bool& exp_most, const amrex::Array4& u, const amrex::Array4& v, const amrex::Array4& cell_data, @@ -67,6 +68,7 @@ void DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain, void DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain, int start_comp, int num_comp, + const bool& exp_most, const amrex::Array4& u, const amrex::Array4& v, const amrex::Array4& cell_data, diff --git a/Source/Diffusion/DiffusionSrcForState_N.cpp b/Source/Diffusion/DiffusionSrcForState_N.cpp index 5e85aa285..d1d949680 100644 --- a/Source/Diffusion/DiffusionSrcForState_N.cpp +++ b/Source/Diffusion/DiffusionSrcForState_N.cpp @@ -37,6 +37,7 @@ using namespace amrex; void DiffusionSrcForState_N (const Box& bx, const Box& domain, int start_comp, int num_comp, + const bool& exp_most, const Array4& u, const Array4& v, const Array4& cell_data, @@ -222,12 +223,9 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, bool ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) && k == 0); bool ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir) && k == dom_hi.z+1); -#ifdef ERF_EXPLICIT_MOST_STRESS - bool most_on_zlo = ( use_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); -#else - bool most_on_zlo = false; -#endif + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + if (ext_dir_on_zlo) { zflux(i,j,k,qty_index) = rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) @@ -277,12 +275,9 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, bool ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) && k == 0); bool ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir) && k == dom_hi.z+1); -#ifdef ERF_EXPLICIT_MOST_STRESS - bool most_on_zlo = ( use_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); -#else - bool most_on_zlo = false; -#endif + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + if (ext_dir_on_zlo) { zflux(i,j,k,qty_index) = rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) @@ -330,12 +325,9 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, bool ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) && k == 0); bool ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir) && k == dom_hi.z+1); -#ifdef ERF_EXPLICIT_MOST_STRESS - bool most_on_zlo = ( use_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); -#else - bool most_on_zlo = false; -#endif + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + if (ext_dir_on_zlo) { zflux(i,j,k,qty_index) = rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) @@ -380,12 +372,9 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, bool ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) && k == 0); bool ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir) && k == dom_hi.z+1); -#ifdef ERF_EXPLICIT_MOST_STRESS - bool most_on_zlo = ( use_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); -#else - bool most_on_zlo = false; -#endif + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + if (ext_dir_on_zlo) { zflux(i,j,k,qty_index) = rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) diff --git a/Source/Diffusion/DiffusionSrcForState_T.cpp b/Source/Diffusion/DiffusionSrcForState_T.cpp index ebd653050..43d41f506 100644 --- a/Source/Diffusion/DiffusionSrcForState_T.cpp +++ b/Source/Diffusion/DiffusionSrcForState_T.cpp @@ -40,6 +40,7 @@ using namespace amrex; void DiffusionSrcForState_T (const Box& bx, const Box& domain, int start_comp, int num_comp, + const bool& exp_most, const Array4& u, const Array4& v, const Array4& cell_data, @@ -249,12 +250,9 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, Real GradCz; bool ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) && k == 0); bool ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir) && k == dom_hi.z+1); -#ifdef ERF_EXPLICIT_MOST_STRESS - bool most_on_zlo = ( use_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); -#else - bool most_on_zlo = false; -#endif + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + if (ext_dir_on_zlo) { GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) @@ -327,12 +325,9 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, Real GradCz; bool ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) && k == 0); bool ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir) && k == dom_hi.z+1); -#ifdef ERF_EXPLICIT_MOST_STRESS - bool most_on_zlo = ( use_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); -#else - bool most_on_zlo = false; -#endif + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + if (ext_dir_on_zlo) { GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) @@ -402,12 +397,9 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, Real GradCz; bool ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) && k == 0); bool ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir) && k == dom_hi.z+1); -#ifdef ERF_EXPLICIT_MOST_STRESS - bool most_on_zlo = ( use_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); -#else - bool most_on_zlo = false; -#endif + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + if (ext_dir_on_zlo) { GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) @@ -476,12 +468,9 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, Real GradCz; bool ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) && k == 0); bool ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir) && k == dom_hi.z+1); -#ifdef ERF_EXPLICIT_MOST_STRESS - bool most_on_zlo = ( use_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); -#else - bool most_on_zlo = false; -#endif + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + if (ext_dir_on_zlo) { GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) diff --git a/Source/Diffusion/EddyViscosity.H b/Source/Diffusion/EddyViscosity.H index 725059157..94607eeee 100644 --- a/Source/Diffusion/EddyViscosity.H +++ b/Source/Diffusion/EddyViscosity.H @@ -19,6 +19,7 @@ ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::MultiFab& const std::unique_ptr& z_phys_nd, const TurbChoice& turbChoice, const amrex::Real const_grav, std::unique_ptr& most, + const bool& exp_most, int level, const amrex::BCRec* bc_ptr, bool vert_only = false); @@ -34,7 +35,8 @@ ComputeSmnSmn (int& i, int& j, int& k, const amrex::Array4& tau13, const amrex::Array4& tau23, const int& klo, - const bool& use_most) + const bool& use_most, + const bool& exp_most) { amrex::Real s11bar = tau11(i,j,k); amrex::Real s22bar = tau22(i,j,k); @@ -54,24 +56,20 @@ ComputeSmnSmn (int& i, int& j, int& k, // - If ERF_EXPLICIT_MOST_STRESS _is_ used, then strains may be correctly // evaluated at the surface and no assumptions are needed. amrex::Real s13bar; -#ifndef ERF_EXPLICIT_MOST_STRESS - if (use_most && k==klo) { + if (use_most && k==klo && !exp_most) { s13bar = 0.5 * ( tau13(i , j , k+1) + tau13(i+1, j , k+1) ); } else -#endif { s13bar = 0.25 * ( tau13(i , j , k ) + tau13(i , j , k+1) + tau13(i+1, j , k ) + tau13(i+1, j , k+1) ); } amrex::Real s23bar; -#ifndef ERF_EXPLICIT_MOST_STRESS - if (use_most && k==klo) { + if (use_most && k==klo && !exp_most) { s23bar = 0.5 * ( tau23(i , j , k+1) + tau23(i , j+1, k+1) ); } else -#endif { s23bar = 0.25 * ( tau23(i , j , k ) + tau23(i , j , k+1) + tau23(i , j+1, k ) + tau23(i , j+1, k+1) ); @@ -82,79 +80,4 @@ ComputeSmnSmn (int& i, int& j, int& k, return SmnSmn; } - -#if 0 -/* This alternative formulation first calculates the product of the strain rate - * tensor components, then averages afterwards. - */ -AMREX_GPU_DEVICE -AMREX_FORCE_INLINE -amrex::Real -ComputeSmnSmn (int& i, int& j, int& k, - const amrex::Array4& tau11, - const amrex::Array4& tau22, - const amrex::Array4& tau33, - const amrex::Array4& tau12, - const amrex::Array4& tau13, - const amrex::Array4& tau23, - const int& klo, - const bool& use_most) -{ - amrex::Real avg_s12_sq = 0.25 * ( tau12(i , j , k )*tau12(i , j , k ) - + tau12(i , j+1, k )*tau12(i , j+1, k ) - + tau12(i+1, j , k )*tau12(i+1, j , k ) - + tau12(i+1, j+1, k )*tau12(i+1, j+1, k ) ); - - // NOTES: - // - If ERF_EXPLICIT_MOST_STRESS is not used, then we do not use the - // strains lying on the bottom boundary with MOST. These values are - // corrupted with the reflect odd BC used to fill the ghost cells. - // - Neglecting the strains on the bottom boundary and using a one- - // sided average implies that the strains are equal at klo and klo+1. - // du/dz and dv/dz can be assumed to be equal through the first cell. - // However, whereas dw/dx = dw/dy = 0 on the surface, this is not true at - // klo+1. Therefore, the strains are not equal in general. - // - If ERF_EXPLICIT_MOST_STRESS _is_ used, then strains may be correctly - // evaluated at the surface and no assumptions are needed. - amrex::Real avg_s13_sq; -#ifndef ERF_EXPLICIT_MOST_STRESS - if (use_most && k==klo) { - avg_s13_sq = 0.5 * ( tau13(i , j , k+1)*tau13(i , j , k+1) - + tau13(i+1, j , k+1)*tau13(i+1, j , k+1) ); - } - else -#endif - { - avg_s13_sq = 0.25 * ( tau13(i , j , k )*tau13(i , j , k ) - + tau13(i , j , k+1)*tau13(i , j , k+1) - + tau13(i+1, j , k )*tau13(i+1, j , k ) - + tau13(i+1, j , k+1)*tau13(i+1, j , k+1) ); - } - - amrex::Real avg_s23_sq; -#ifndef ERF_EXPLICIT_MOST_STRESS - if (use_most && k==klo) { - avg_s23_sq = 0.5 * ( tau23(i , j , k+1)*tau23(i , j , k+1) - + tau23(i , j+1, k+1)*tau23(i , j+1, k+1) ); - } - else -#endif - { - avg_s23_sq = 0.25 * ( tau23(i , j , k )*tau23(i , j , k ) - + tau23(i , j , k+1)*tau23(i , j , k+1) - + tau23(i , j+1, k )*tau23(i , j+1, k ) - + tau23(i , j+1, k+1)*tau23(i , j+1, k+1) ); - } - - amrex::Real SmnSmn = tau11(i,j,k)*tau11(i,j,k) - + tau22(i,j,k)*tau22(i,j,k) - + tau33(i,j,k)*tau33(i,j,k) - + 2.0*avg_s12_sq - + 2.0*avg_s13_sq - + 2.0*avg_s23_sq; - - return SmnSmn; -} -#endif - #endif diff --git a/Source/ERF.H b/Source/ERF.H index bddfc3765..aa85b3117 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -625,14 +625,12 @@ private: std::unique_ptr micro; amrex::Vector> qmoist; // (lev,ncomp) This has up to 8 components: qt, qv, qc, qi, qp, qr, qs, qg -//#if defined(ERF_USE_WINDFARM) amrex::Vector Nturb; amrex::Vector vars_fitch; // Vabs, Vabsdt, dudt, dvdt, dTKEdt amrex::Vector vars_ewp; // dudt, dvdt, dTKEdt amrex::Real hub_height, rotor_dia, thrust_coeff_standing, nominal_power; amrex::Vector wind_speed, thrust_coeff, power; -//#endif LandSurface lsm; amrex::Vector> lsm_data; // (lev,ncomp) Components: theta, q1, q2 diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 321046877..e1991f4e9 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -892,11 +892,12 @@ ERF::InitData () // FillPatch does not call MOST, FillIntermediatePatch does. if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST) { -#ifdef ERF_EXPLICIT_MOST_STRESS - Print() << "Using MOST with explicitly included surface stresses" << std::endl; -#endif + bool use_exp_most = solverChoice.use_explicit_most; + if (use_exp_most) { + Print() << "Using MOST with explicitly included surface stresses" << std::endl; + } - m_most = std::make_unique(geom, vars_old, Theta_prim, Qv_prim, z_phys_nd, + m_most = std::make_unique(geom, use_exp_most, vars_old, Theta_prim, Qv_prim, z_phys_nd, sst_lev, lmask_lev, lsm_data, lsm_flux #ifdef ERF_USE_NETCDF ,start_bdy_time, bdy_time_interval diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 35670f1bb..13102fc22 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -328,7 +328,7 @@ ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMap } if (l_use_kturb) { - eddyDiffs_lev[lev] = std::make_unique( ba, dm, EddyDiff::NumDiffs, 1 ); + eddyDiffs_lev[lev] = std::make_unique( ba, dm, EddyDiff::NumDiffs, 1 ); eddyDiffs_lev[lev]->setVal(0.0); if(l_use_ddorf) { SmnSmn_lev[lev] = std::make_unique( ba, dm, 1, 0 ); diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index 5f434b02d..ae74c83d9 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -92,6 +92,7 @@ void ERF::advance_dycore(int level, (tc.pbl_type != PBLType::None) ); const bool use_most = (m_most != nullptr); + const bool exp_most = (solverChoice.use_explicit_most); amrex::ignore_unused(use_most); const BoxArray& ba = state_old[IntVars::cons].boxArray(); @@ -199,7 +200,7 @@ void ERF::advance_dycore(int level, *eddyDiffs, *Hfx1, *Hfx2, *Hfx3, *Diss, // to be updated fine_geom, *mapfac_u[level], *mapfac_v[level], z_phys_nd[level], tc, solverChoice.gravity, - m_most, level, bc_ptr_d); + m_most, exp_most, level, bc_ptr_d); } // *********************************************************************************************** diff --git a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp index a36c2ed9f..c6df98725 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp @@ -127,6 +127,7 @@ void erf_slow_rhs_inc (int level, int nrk, tc.pbl_type == PBLType::YSU ); const bool use_most = (most != nullptr); + const bool exp_most = (solverChoice.use_explicit_most); const Box& domain = geom.Domain(); const int domhi_z = domain.bigEnd(2); @@ -305,17 +306,15 @@ void erf_slow_rhs_inc (int level, int nrk, }); } -#ifdef ERF_EXPLICIT_MOST_STRESS // We've updated the strains at all locations including the // surface. This is required to get the correct strain-rate // magnitude. Now, update the stress everywhere but the surface // to retain the values set by MOST. - if (use_most) { + if (use_most && exp_most) { // Don't overwrite modeled total stress value at boundary tbxxz.setSmall(2,1); tbxyz.setSmall(2,1); } -#endif // ***************************************************************************** // Stress tensor compute terrain @@ -415,6 +414,16 @@ void erf_slow_rhs_inc (int level, int nrk, }); } + // We've updated the strains at all locations including the + // surface. This is required to get the correct strain-rate + // magnitude. Now, update the stress everywhere but the surface + // to retain the values set by MOST. + if (use_most && exp_most) { + // Don't overwrite modeled total stress value at boundary + tbxxz.setSmall(2,1); + tbxyz.setSmall(2,1); + } + // ***************************************************************************** // Stress tensor compute no terrain // ***************************************************************************** diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index c67ef2df8..22e4067c2 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -120,6 +120,7 @@ void erf_slow_rhs_post (int level, int finest_level, tc.les_type == LESType::Deardorff || tc.pbl_type == PBLType::MYNN25 || tc.pbl_type == PBLType::YSU ); + const bool exp_most = (solverChoice.use_explicit_most); const Box& domain = geom.Domain(); @@ -352,7 +353,7 @@ void erf_slow_rhs_post (int level, int finest_level, const Array4 tm_arr = t_mean_mf ? t_mean_mf->const_array(mfi) : Array4{}; if (l_use_terrain) { - DiffusionSrcForState_T(tbx, domain, start_comp, num_comp, u, v, + DiffusionSrcForState_T(tbx, domain, start_comp, num_comp, exp_most, u, v, cur_cons, cur_prim, cell_rhs, diffflux_x, diffflux_y, diffflux_z, z_nd, ax_arr, ay_arr, az_arr, detJ_arr, @@ -360,7 +361,7 @@ void erf_slow_rhs_post (int level, int finest_level, hfx_z, diss, mu_turb, dc, tc, tm_arr, grav_gpu, bc_ptr_d, use_most); } else { - DiffusionSrcForState_N(tbx, domain, start_comp, num_comp, u, v, + DiffusionSrcForState_N(tbx, domain, start_comp, num_comp, exp_most, u, v, cur_cons, cur_prim, cell_rhs, diffflux_x, diffflux_y, diffflux_z, dxInv, SmnSmn_a, mf_m, mf_u, mf_v, diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index e1d57d88e..b613b79cb 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -143,6 +143,7 @@ void erf_slow_rhs_pre (int level, int finest_level, const bool use_moisture = (solverChoice.moisture_type != MoistureType::None); const bool use_most = (most != nullptr); + const bool exp_most = (solverChoice.use_explicit_most); const Box& domain = geom.Domain(); const int domhi_z = domain.bigEnd(2); @@ -318,21 +319,19 @@ void erf_slow_rhs_pre (int level, int finest_level, SmnSmn_a = SmnSmn->array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23,domlo_z,use_most); + SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23,domlo_z,use_most,exp_most); }); } -#ifdef ERF_EXPLICIT_MOST_STRESS // We've updated the strains at all locations including the // surface. This is required to get the correct strain-rate // magnitude. Now, update the stress everywhere but the surface // to retain the values set by MOST. - if (use_most) { + if (use_most && exp_most) { // Don't overwrite modeled total stress value at boundary tbxxz.setSmall(2,1); tbxyz.setSmall(2,1); } -#endif // ***************************************************************************** // Stress tensor compute terrain @@ -428,21 +427,19 @@ void erf_slow_rhs_pre (int level, int finest_level, SmnSmn_a = SmnSmn->array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23,domlo_z,use_most); + SmnSmn_a(i,j,k) = ComputeSmnSmn(i,j,k,s11,s22,s33,s12,s13,s23,domlo_z,use_most,exp_most); }); } -#ifdef ERF_EXPLICIT_MOST_STRESS // We've updated the strains at all locations including the // surface. This is required to get the correct strain-rate // magnitude. Now, update the stress everywhere but the surface // to retain the values set by MOST. - if (use_most) { + if (use_most && exp_most) { // Don't overwrite modeled total stress value at boundary tbxxz.setSmall(2,1); tbxyz.setSmall(2,1); } -#endif // ***************************************************************************** // Stress tensor compute no terrain @@ -732,7 +729,7 @@ void erf_slow_rhs_pre (int level, int finest_level, int n_comp = end_comp - n_start + 1; if (l_use_terrain) { - DiffusionSrcForState_T(bx, domain, n_start, n_comp, u, v, + DiffusionSrcForState_T(bx, domain, n_start, n_comp, exp_most, u, v, cell_data, cell_prim, cell_rhs, diffflux_x, diffflux_y, diffflux_z, z_nd, ax_arr, ay_arr, az_arr, detJ_arr, @@ -740,7 +737,7 @@ void erf_slow_rhs_pre (int level, int finest_level, hfx_z, diss, mu_turb, dc, tc, tm_arr, grav_gpu, bc_ptr_d, use_most); } else { - DiffusionSrcForState_N(bx, domain, n_start, n_comp, u, v, + DiffusionSrcForState_N(bx, domain, n_start, n_comp, exp_most, u, v, cell_data, cell_prim, cell_rhs, diffflux_x, diffflux_y, diffflux_z, dxInv, SmnSmn_a, mf_m, mf_u, mf_v,