diff --git a/Docs/sphinx_doc/MOST.rst b/Docs/sphinx_doc/MOST.rst index dc87adfe6..13eee85a0 100644 --- a/Docs/sphinx_doc/MOST.rst +++ b/Docs/sphinx_doc/MOST.rst @@ -136,6 +136,10 @@ In ERF, when the MOST boundary condition is applied, velocity and temperature in (\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] + 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 + (``zlo.type = "Most"``) should be cautioned. + MOST Inputs ~~~~~~~~~~~~~~~~~~~ To evaluate the fluxes with MOST, the surface rougness parameter :math:`z_{0}` must be specified. This quantity may be considered a constant or may be parameterized through the friction velocity :math:`u_{\star}`. ERF supports three methods for parameterizing the surface roughness: ``constant``, ``charnock``, and ``modified_charnock``. The latter two methods parameterize :math:`z_{0} = f(u_{\star})` and are described in `Jimenez & Dudhia, American Meteorological Society, 2018 `_. The rougness calculation method may be specified with diff --git a/Exec/RegTests/Terrain3d_Hemisphere/inputs_most_test b/Exec/RegTests/Terrain3d_Hemisphere/inputs_most_test index 063fbf980..cca193eb2 100644 --- a/Exec/RegTests/Terrain3d_Hemisphere/inputs_most_test +++ b/Exec/RegTests/Terrain3d_Hemisphere/inputs_most_test @@ -18,7 +18,7 @@ zhi.type = "SlipWall" # MOST BOUNDARY #================================================================= zlo.type = "Most" -erf.most.average_policy = 1 # POLICY FOR AVERAGING +erf.most.average_policy = 0 # POLICY FOR AVERAGING erf.most.use_normal_vector = false # USE NORMAL VECTOR W/ TERRAIN erf.most.use_interpolation = true # USE INTERPOLATION ON DESTINATION #erf.most.time_average = true # USE TIME AVERAGING @@ -81,3 +81,7 @@ erf.alpha_T = 0.0 # [m^2/s] prob.T_0 = 300.0 prob.U_0 = 1.0 prob.rho_0 = 1.16 +prob.U_0_Pert_Mag = 0.01 +prob.V_0_Pert_Mag = 0.01 +prob.W_0_Pert_Mag = 0.0 +prob.pert_ref_height = 1.0 \ No newline at end of file diff --git a/Exec/RegTests/Terrain3d_Hemisphere/prob.H b/Exec/RegTests/Terrain3d_Hemisphere/prob.H index 64e5368e3..dcc81141e 100644 --- a/Exec/RegTests/Terrain3d_Hemisphere/prob.H +++ b/Exec/RegTests/Terrain3d_Hemisphere/prob.H @@ -12,12 +12,30 @@ struct ProbParm : ProbParmDefaults { amrex::Real T_0 = 300.0; // surface temperature == mean potential temperature amrex::Real U_0 = 0.0; amrex::Real V_0 = 0.0; + + // random initial perturbations (legacy code) + amrex::Real U_0_Pert_Mag = 0.0; + amrex::Real V_0_Pert_Mag = 0.0; + amrex::Real W_0_Pert_Mag = 0.0; + + // divergence-free initial perturbations + amrex::Real pert_deltaU = 0.0; + amrex::Real pert_deltaV = 0.0; + amrex::Real pert_periods_U = 5.0; + amrex::Real pert_periods_V = 5.0; + amrex::Real pert_ref_height = 1.0; + + // helper vars + amrex::Real aval; + amrex::Real bval; + amrex::Real ufac; + amrex::Real vfac; }; // namespace ProbParm class Problem : public ProblemBase { public: - Problem(); + Problem (const amrex::Real* problo, const amrex::Real* probhi); #include "Prob/init_density_hse_dry_terrain.H" diff --git a/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp b/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp index 0acce57b3..19a1e8818 100644 --- a/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp +++ b/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp @@ -4,25 +4,33 @@ using namespace amrex; std::unique_ptr -amrex_probinit( - const amrex_real* /*problo*/, - const amrex_real* /*probhi*/) +amrex_probinit (const amrex_real* problo, + const amrex_real* probhi) { - return std::make_unique(); + return std::make_unique(problo, probhi); } -Problem::Problem() +Problem::Problem (const amrex::Real* problo, + const amrex::Real* probhi) { // Parse params amrex::ParmParse pp("prob"); pp.query("rho_0", parms.rho_0); pp.query("U_0", parms.U_0); + pp.query("U_0_Pert_Mag", parms.U_0_Pert_Mag); + pp.query("V_0_Pert_Mag", parms.V_0_Pert_Mag); + pp.query("W_0_Pert_Mag", parms.W_0_Pert_Mag); + pp.query("pert_ref_height", parms.pert_ref_height); + parms.aval = parms.pert_periods_U * 2.0 * PI / (probhi[1] - problo[1]); + parms.bval = parms.pert_periods_V * 2.0 * PI / (probhi[0] - problo[0]); + parms.ufac = parms.pert_deltaU * std::exp(0.5) / parms.pert_ref_height; + parms.vfac = parms.pert_deltaV * std::exp(0.5) / parms.pert_ref_height; init_base_parms(parms.rho_0, parms.T_0); } void -Problem::init_custom_pert( +Problem::init_custom_pert ( const Box& bx, const Box& xbx, const Box& ybx, @@ -34,66 +42,107 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const& z_nd, - Array4 const& z_cc, + Array4 const& /*z_cc*/, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, Array4 const& /*mf_v*/, const SolverChoice& sc) { - const int khi = geomdata.Domain().bigEnd()[2]; + const int khi = geomdata.Domain().bigEnd()[2]; const bool use_moisture = (sc.moisture_type != MoistureType::None); + const bool use_terrain = sc.use_terrain; - AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); - // Geometry (note we must include these here to get the data on device) - amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - const auto prob_lo = geomdata.ProbLo(); - const auto dx = geomdata.CellSize(); - const Real x = prob_lo[0] + (i + 0.5) * dx[0]; - const Real z = z_cc(i,j,k); - - // Set scalar = 0 everywhere - state(i, j, k, RhoScalar_comp) = 0.0; + // Geometry (note we must include these here to get the data on device) + ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Set scalar = 0 everywhere + state(i, j, k, RhoScalar_comp) = 0.0; - if (use_moisture) { - state(i, j, k, RhoQ1_comp) = 0.0; - state(i, j, k, RhoQ2_comp) = 0.0; - } - }); - - // Set the x-velocity - amrex::ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - x_vel(i, j, k) = 10.0; - }); - - // Set the y-velocity - amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - y_vel(i, j, k) = 0.0; - }); - - const auto dx = geomdata.CellSize(); - amrex::GpuArray dxInv; - dxInv[0] = 1. / dx[0]; - dxInv[1] = 1. / dx[1]; - dxInv[2] = 1. / dx[2]; - - // Set the z-velocity from impenetrable condition - amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - z_vel(i, j, k) = 0.0; - }); - - amrex::Gpu::streamSynchronize(); + if (use_moisture) { + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + } + }); + // Set the x-velocity + ParallelForRNG(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept + { + const Real* prob_lo = geomdata.ProbLo(); + const Real* dx = geomdata.CellSize(); + const Real z = use_terrain ? 0.25*( z_nd(i,j ,k) + z_nd(i,j ,k+1) + + z_nd(i,j+1,k) + z_nd(i,j+1,k+1) ) + : prob_lo[2] + (k + 0.5) * dx[2]; + + x_vel(i, j, k) = 10.0; + if ((z <= parms.pert_ref_height) && (parms.U_0_Pert_Mag != 0.0)) + { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + Real x_vel_prime = (rand_double*2.0 - 1.0)*parms.U_0_Pert_Mag; + x_vel(i, j, k) += x_vel_prime; + } + }); + + // Set the y-velocity + ParallelForRNG(ybx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept + { + const Real* prob_lo = geomdata.ProbLo(); + const Real* dx = geomdata.CellSize(); + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real z = use_terrain ? 0.25*( z_nd(i ,j,k) + z_nd(i ,j,k+1) + + z_nd(i+1,j,k) + z_nd(i+1,j,k+1) ) + : prob_lo[2] + (k + 0.5) * dx[2]; + + // Set the y-velocity + y_vel(i, j, k) = 0.0; + if ((z <= parms.pert_ref_height) && (parms.V_0_Pert_Mag != 0.0)) + { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + Real y_vel_prime = (rand_double*2.0 - 1.0)*parms.V_0_Pert_Mag; + y_vel(i, j, k) += y_vel_prime; + } + if (parms.pert_deltaV != 0.0) + { + const amrex::Real xl = x - prob_lo[0]; + const amrex::Real zl = z / parms.pert_ref_height; + const amrex::Real damp = std::exp(-0.5 * zl * zl); + y_vel(i, j, k) += parms.vfac * damp * z * std::cos(parms.bval * xl); + } + }); + + const auto dx = geomdata.CellSize(); + amrex::GpuArray dxInv; + dxInv[0] = 1. / dx[0]; + dxInv[1] = 1. / dx[1]; + dxInv[2] = 1. / dx[2]; + + // Set the z-velocity from impenetrable condition + ParallelForRNG(zbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept + { + const int dom_lo_z = geomdata.Domain().smallEnd()[2]; + const int dom_hi_z = geomdata.Domain().bigEnd()[2]; + + // Set the z-velocity + if (k == dom_lo_z || k == dom_hi_z+1) + { + z_vel(i, j, k) = 0.0; + } + else if (parms.W_0_Pert_Mag != 0.0) + { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + Real z_vel_prime = (rand_double*2.0 - 1.0)*parms.W_0_Pert_Mag; + z_vel(i, j, k) = z_vel_prime; + } + }); + + amrex::Gpu::streamSynchronize(); } void -Problem::init_custom_terrain( +Problem::init_custom_terrain ( const Geometry& geom, MultiFab& z_phys_nd, const Real& /*time*/) @@ -111,7 +160,6 @@ Problem::init_custom_terrain( // User function parameters Real a = 0.5; - Real num = 8 * a * a * a; Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]); Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]); diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H index 95e93fe3f..4828890ba 100644 --- a/Source/BoundaryConditions/ABLMost.H +++ b/Source/BoundaryConditions/ABLMost.H @@ -286,9 +286,9 @@ public: MODIFIED_CHARNOCK }; - FluxCalcType flux_type; - ThetaCalcType theta_type; - RoughCalcType rough_type; + FluxCalcType flux_type{FluxCalcType::MOENG}; + ThetaCalcType theta_type{ThetaCalcType::ADIABATIC}; + RoughCalcType rough_type{RoughCalcType::CONSTANT}; private: bool use_moisture; diff --git a/Source/BoundaryConditions/ABLMost.cpp b/Source/BoundaryConditions/ABLMost.cpp index e5cc85748..a13cee3d4 100644 --- a/Source/BoundaryConditions/ABLMost.cpp +++ b/Source/BoundaryConditions/ABLMost.cpp @@ -174,7 +174,8 @@ ABLMost::compute_most_bcs (const int& lev, 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 u_mag_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); @@ -202,7 +203,6 @@ ABLMost::compute_most_bcs (const int& lev, Box b2d = bx; b2d.setBig(2,zlo-1); int n = RhoTheta_comp; - ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real dz = (zphys_arr) ? ( zphys_arr(i,j,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain; diff --git a/Source/BoundaryConditions/MOSTAverage.cpp b/Source/BoundaryConditions/MOSTAverage.cpp index 226c9f6fb..cc82f1724 100644 --- a/Source/BoundaryConditions/MOSTAverage.cpp +++ b/Source/BoundaryConditions/MOSTAverage.cpp @@ -317,13 +317,15 @@ MOSTAverage::set_k_indices_T() for (int lev(0); lev < m_maxlev; lev++) { int kmax = m_geom[lev].Domain().bigEnd(2); for (MFIter mfi(*m_k_indx[lev], TileNoZ()); mfi.isValid(); ++mfi) { - Box npbx = mfi.tilebox(); npbx.convert({1,1,0}); + Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0)); const auto z_phys_arr = m_z_phys_nd[lev]->const_array(mfi); auto k_arr = m_k_indx[lev]->array(mfi); ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { k_arr(i,j,k) = 0; - Real z_target = d_zref + z_phys_arr(i,j,k); + Real z_bot_face = 0.25 * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k) + + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) ); + Real z_target = z_bot_face + d_zref; for (int lk(0); lk<=kmax; ++lk) { Real z_lo = 0.25 * ( z_phys_arr(i,j ,lk ) + z_phys_arr(i+1,j ,lk ) + z_phys_arr(i,j+1,lk ) + z_phys_arr(i+1,j+1,lk ) ); @@ -392,7 +394,9 @@ MOSTAverage::set_norm_indices_T() j_arr(i,j,k) = j_new; // Search for k (grid is stretched in z) - Real z_target = delta_z + z_phys_arr(i,j,k); + Real z_bot_face = 0.25 * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k) + + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) ); + Real z_target = z_bot_face + delta_z; for (int lk(0); lk<=kmax; ++lk) { Real z_lo = 0.25 * ( z_phys_arr(i_new,j_new ,lk ) + z_phys_arr(i_new+1,j_new ,lk ) + z_phys_arr(i_new,j_new+1,lk ) + z_phys_arr(i_new+1,j_new+1,lk ) ); @@ -479,7 +483,7 @@ MOSTAverage::set_norm_positions_T() const auto dxInv = m_geom[lev].InvCellSizeArray(); IntVect ng = m_x_pos[lev]->nGrowVect(); ng[2]=0; for (MFIter mfi(*m_x_pos[lev], TileNoZ()); mfi.isValid(); ++mfi) { - Box npbx = mfi.tilebox(); npbx.convert({1,1,0}); + Box npbx = mfi.tilebox(IntVect(1,1,0),IntVect(1,1,0)); Box gpbx = mfi.growntilebox(ng); RealBox grb{gpbx,dx.data(),base.dataPtr()}; @@ -506,7 +510,9 @@ MOSTAverage::set_norm_positions_T() // Final position at end of vector x_pos_arr(i,j,k) = x0 + delta_x; y_pos_arr(i,j,k) = y0 + delta_y; - z_pos_arr(i,j,k) = z_phys_arr(i,j,k) + delta_z; + Real z_bot_face = 0.25 * ( z_phys_arr(i ,j ,k) + z_phys_arr(i+1,j ,k) + + z_phys_arr(i ,j+1,k) + z_phys_arr(i+1,j+1,k) ); + z_pos_arr(i,j,k) = z_bot_face + delta_z; // Destination position must be contained on the current process! Real pos[] = {x_pos_arr(i,j,k),y_pos_arr(i,j,k),0.5*dx[2]}; diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H index 2a01130d6..5e78bf855 100644 --- a/Source/BoundaryConditions/MOSTStress.H +++ b/Source/BoundaryConditions/MOSTStress.H @@ -627,7 +627,7 @@ struct moeng_flux rho = cons_arr(ic,jc,zlo,Rho_comp); qv = cons_arr(ic,jc,zlo,RhoQ1_comp) / rho; eta = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s] - eta = amrex::max(eta,eps); + eta = amrex::max(eta,eta_eps); // TODO: Integrate MOST with moisture and MOENG FLUX type amrex::Real moflux = 0.0; @@ -686,7 +686,7 @@ struct moeng_flux rho = cons_arr(ic,jc,zlo,Rho_comp); theta = cons_arr(ic,jc,zlo,RhoTheta_comp) / rho; eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s] - eta = amrex::max(eta,eps); + eta = amrex::max(eta,eta_eps); amrex::Real theta_mean = tm_arr(ic,jc,zlo); amrex::Real wsp_mean = umm_arr(ic,jc,zlo); @@ -749,7 +749,7 @@ struct moeng_flux + cons_arr(ic ,jc,zlo,Rho_comp) ); eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v) + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) ); - eta = amrex::max(eta,eps); + eta = amrex::max(eta,eta_eps); amrex::Real umean = um_arr(i,j,zlo); amrex::Real wsp_mean = 0.5 * ( umm_arr(ic-1,jc,zlo) + umm_arr(ic,jc,zlo) ); @@ -811,7 +811,7 @@ struct moeng_flux + cons_arr(ic,jc ,zlo,Rho_comp) ); eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v) + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) ); - eta = amrex::max(eta,eps); + eta = amrex::max(eta,eta_eps); amrex::Real vmean = vm_arr(i,j,zlo); amrex::Real wsp_mean = 0.5 * ( umm_arr(ic,jc-1,zlo) + umm_arr(ic,jc,zlo) ); @@ -832,7 +832,8 @@ struct moeng_flux private: int zlo; - const amrex::Real eps = 1e-16; + const amrex::Real eps = 1e-15; + const amrex::Real eta_eps = 1e-8; }; @@ -881,7 +882,7 @@ struct donelan_flux rho = cons_arr(ic,jc,zlo,Rho_comp); qv = cons_arr(ic,jc,zlo,RhoQ1_comp) / rho; eta = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s] - eta = amrex::max(eta,eps); + eta = amrex::max(eta,eta_eps); // TODO: Integrate MOST with moisture and DONELAN FLUX type amrex::Real moflux = 0.0; @@ -928,7 +929,7 @@ struct donelan_flux rho = cons_arr(ic,jc,zlo,Rho_comp); theta = cons_arr(ic,jc,zlo,RhoTheta_comp) / rho; eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s] - eta = amrex::max(eta,eps); + eta = amrex::max(eta,eta_eps); amrex::Real Cd = 0.0012; amrex::Real wsp_mean = umm_arr(ic,jc,zlo); @@ -985,7 +986,7 @@ struct donelan_flux + cons_arr(ic ,jc,zlo,Rho_comp) ); eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v) + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) ); - eta = amrex::max(eta,eps); + eta = amrex::max(eta,eta_eps); amrex::Real Cd = 0.001; const amrex::Real c = 7e-5; @@ -1047,7 +1048,7 @@ struct donelan_flux + cons_arr(ic,jc ,zlo,Rho_comp) ); eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v) + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) ); - eta = amrex::max(eta,eps); + eta = amrex::max(eta,eta_eps); amrex::Real Cd = 0.001; const amrex::Real c = 7e-5; @@ -1068,7 +1069,7 @@ struct donelan_flux private: int zlo; - const amrex::Real eps = 1e-16; + const amrex::Real eta_eps = 1e-8; }; @@ -1117,7 +1118,7 @@ struct custom_flux rho = cons_arr(ic,jc,zlo,Rho_comp); qv = cons_arr(ic,jc,zlo,RhoQ1_comp) / rho; eta = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s] - eta = amrex::max(eta,eps); + eta = amrex::max(eta,eta_eps); amrex::Real qstar = q_star_arr(ic,jc,zlo); amrex::Real moflux = (std::abs(qstar) > eps) ? -rho * qstar : 0.0; @@ -1164,7 +1165,7 @@ struct custom_flux rho = cons_arr(ic,jc,zlo,Rho_comp); theta = cons_arr(ic,jc,zlo,RhoTheta_comp) / rho; eta = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s] - eta = amrex::max(eta,eps); + eta = amrex::max(eta,eta_eps); amrex::Real tstar = t_star_arr(ic,jc,zlo); amrex::Real moflux = (std::abs(tstar) > eps) ? -rho * tstar : 0.0; @@ -1218,12 +1219,8 @@ struct custom_flux + cons_arr(ic ,jc,zlo,Rho_comp) ); eta = 0.5 *( eta_arr(ie-1,je,zlo,EddyDiff::Mom_v) + eta_arr(ie ,je,zlo,EddyDiff::Mom_v) ); - eta = amrex::max(eta,eps); + eta = amrex::max(eta,eta_eps); - // Note: The surface mean shear stress is decomposed into tau_xz by - // multiplying the modeled shear stress (rho*ustar^2) with - // a factor of umean/wsp_mean for directionality; this factor - // modifies the demoninator from what is in Moeng 1984. amrex::Real ustar = 0.5 * ( u_star_arr(ic-1,jc,zlo) + u_star_arr(ic,jc,zlo) ); amrex::Real wsp = sqrt(velx*velx+vely*vely); amrex::Real stressx = rho * ustar * ustar * velx / wsp; @@ -1275,12 +1272,8 @@ struct custom_flux + cons_arr(ic,jc ,zlo,Rho_comp) ); eta = 0.5*( eta_arr(ie,je-1,zlo,EddyDiff::Mom_v) + eta_arr(ie,je ,zlo,EddyDiff::Mom_v) ); - eta = amrex::max(eta,eps); + eta = amrex::max(eta,eta_eps); - // Note: The surface mean shear stress is decomposed into tau_yz by - // multiplying the modeled shear stress (rho*ustar^2) with - // a factor of vmean/wsp_mean for directionality; this factor - // modifies the demoninator from what is in Moeng 1984. amrex::Real ustar = 0.5 * ( u_star_arr(ic,jc-1,zlo) + u_star_arr(ic,jc,zlo) ); amrex::Real wsp = sqrt(velx*velx+vely*vely); amrex::Real stressy = rho * ustar * ustar * vely / wsp; @@ -1291,6 +1284,7 @@ struct custom_flux private: int zlo; - const amrex::Real eps = 1e-16; + const amrex::Real eps = 1e-15; + const amrex::Real eta_eps = 1e-8; }; #endif diff --git a/Source/Prob/init_density_hse_dry_terrain.H b/Source/Prob/init_density_hse_dry_terrain.H index d78598aa8..c22a2fe29 100644 --- a/Source/Prob/init_density_hse_dry_terrain.H +++ b/Source/Prob/init_density_hse_dry_terrain.H @@ -8,7 +8,7 @@ */ void erf_init_dens_hse (amrex::MultiFab& rho_hse, - std::unique_ptr& z_phys_nd, + std::unique_ptr& /*z_phys_nd*/, std::unique_ptr& z_phys_cc, amrex::Geometry const& geom) override {