Skip to content

Commit

Permalink
Explicit MOST runtime (#1605)
Browse files Browse the repository at this point in the history 
* Make explicit most a runtime option.

* Fix function signature mismatch.

* Remove commented flag for windfarm.

* Add docs update.

* Fix SmnSmn since it was a ifndef.
  • Loading branch information
@AMLattanzi
AMLattanzi authored May 7, 2024
1 parent 1a5597f commit 7ed2e8f
Show file tree
Hide file tree
Showing 22 changed files with 304 additions and 485 deletions.
4 changes: 0 additions & 4 deletions CMake/BuildERFExe.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 0 additions & 2 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
9 changes: 8 additions & 1 deletion Docs/sphinx_doc/MOST.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <https://journals.ametsoc.org/view/journals/atsc/41/13/1520-0469_1984_041_2052_alesmf_2_0_co_2.xml>`_; 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
Expand Down
1 change: 0 additions & 1 deletion Exec/ABL/GNUmakefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,6 @@ TEST = TRUE
USE_ASSERTION = TRUE

#USE_POISSON_SOLVE = TRUE
USE_EXPLICIT_MOST_STRESS = FALSE

# GNU Make
Bpack := ./Make.package
Expand Down
4 changes: 0 additions & 4 deletions Exec/Make.ERF
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
15 changes: 6 additions & 9 deletions Source/BoundaryConditions/ABLMost.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ public:

// Constructor
explicit ABLMost (const amrex::Vector<amrex::Geometry>& geom,
bool& use_exp_most,
amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars_old,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Theta_prim,
amrex::Vector<std::unique_ptr<amrex::MultiFab>>& Qv_prim,
Expand All @@ -40,7 +41,8 @@ public:
amrex::Vector<amrex::Vector<amrex::MultiFab*>> 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)
Expand Down Expand Up @@ -220,26 +222,20 @@ public:
void
impose_most_bcs (const int& lev,
const amrex::Vector<amrex::MultiFab*>& 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<typename FluxCalc>
void
compute_most_bcs (const int& lev,
const amrex::Vector<amrex::MultiFab*>& 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);
Expand Down Expand Up @@ -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};
Expand Down
106 changes: 29 additions & 77 deletions Source/BoundaryConditions/ABLMost.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,48 +121,36 @@ ABLMost::compute_fluxes (const int& lev,
void
ABLMost::impose_most_bcs (const int& lev,
const Vector<MultiFab*>& 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);
}
}
Expand All @@ -180,13 +168,10 @@ template<typename FluxCalc>
void
ABLMost::compute_most_bcs (const int& lev,
const Vector<MultiFab*>& 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)
Expand All @@ -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<Real>{};
auto t32_arr = (zymom_flux) ? zymom_flux->array(mfi) : Array4<Real>{};
auto hfx_arr = heat_flux->array(mfi);
#else

auto t13_arr = (m_exp_most) ? xzmom_flux->array(mfi) : Array4<Real>{};
auto t23_arr = (m_exp_most) ? yzmom_flux->array(mfi) : Array4<Real>{};
auto t31_arr = (zxmom_flux && m_exp_most) ? zxmom_flux->array(mfi) : Array4<Real>{};
auto t32_arr = (zymom_flux && m_exp_most) ? zymom_flux->array(mfi) : Array4<Real>{};
auto hfx_arr = (m_exp_most) ? heat_flux->array(mfi) : Array4<Real>{};

const auto eta_arr = eddyDiffs->array(mfi);
#endif

const auto zphys_arr = (z_phys) ? z_phys->const_array(mfi) : Array4<const Real>{};

// 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
Expand All @@ -254,54 +241,35 @@ 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
int is_land = (lmask_arr) ? lmask_arr(i,j,klo) : 1;
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
if ((flux_type == FluxCalcType::CUSTOM) && use_moisture) {
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);
});
}
Expand All @@ -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);
});

Expand All @@ -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);
});
}
Expand Down
3 changes: 0 additions & 3 deletions Source/BoundaryConditions/ERF_FillPatch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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());
}

Expand Down
Loading

0 comments on commit 7ed2e8f

Please sign in to comment.