Skip to content

Commit

Permalink
Provide alternate solution path for imposing MOST (#1472)
Browse files Browse the repository at this point in the history 
* Add ERF_EXPLICIT_MOST_STRESS flag

enable in CMake with ERF_ENABLE_EXPLICIT_MOST_STRESS

* Add preprocessor directives for ERF_EXPLICIT_MOST_STRESS

* Pass in second cell dz for terrain case

* Set MOST ghost cell by extrapolating from first cell

This is how the surface BC is set in SOWFA and is a simplification of Moeng
1984, in which the first-cell velocity gradient -- used to update the tke
production term -- is determined by interpolating between the resolved velocity
gradient at the first zface and the modeled velocity gradient at a test height
below the first cell center. We neglect the latter contribution for now.

* Remove dependence on eta

* Return shear stress from MOST

* Pass Tau13, Tau23, and SFS_hfx3 to the shear stress update

* Oops

* Set surface stresses based on MOST

* Fix bug and make shear stress models consistent

- custom_flux.compute_t_flux was incorrectly setting `moflux`
- off-setting sign errors on `moflux` for moeng_flux; convention is for a
  +ve heat flux to correspond to tstar and dtheta/dz < 0
- fix sign on value set by donelan_flux.compute_t_flux
- all stressx and stressy are multiplied by rho

* Note if solver was compiled with ERF_EXPLICIT_MOST_STRESS

* Set surface heat flux based on MOST

* Super minor cleanup

* Pass use_most flag to DiffusionSrcForState

* Only overwrite surface flux if compiled with ERF_EXPLICIT_MOST_STRESS

* Fix problem with valid box

* Print out custom most info

* Update SmnSmn calculation to use valid surface strains if available

* Reenable wall treatment for ComputeSmnSmn

NOTE: the current implementation assumes that dw/dx at the first zface above
the surface is 0, which isn't strictly correct. But the expectation is that
du/dz >> dw/dx in general

* Surface temperature BC based on extrapolating dtheta/dz rather than d(rhotheta)/dz

* Use first-order one-sided calc to be consistent with BC

* Cleanup

* Don't calculate strain at k=0 prior to first RK stage (unneeded)

* Set velocity on boundary, not momenta

* Fingers crossed

When calculating the diffusion source, calculate the strain rate everywhere
before resizing the box for the stress calc. This provides the correct strain
rate on the surface, without making any implicit assumption about dw/dx, dw/dy.

* Add alternative ComputeSmnSmn formulation

This differs from the original formulation, which matches what is in the arw_v4
manual. Originally, the strain rates are averaged to cell centered and then
squared. The alternative (not compiled by default) is to first calculate the
square of the strain-rate tensor components and then average them.

* Add comments about hfx

* Store surface heat flux in ghost cell

The kinematic heat flux SFS_hfx3_lev now has dual uses.

- The surface value (at k=-1) is used when calculating the diffusive z-flux at
  k=0 (zflux is staggered in z). This is the modeled flux that comes from
  MOSTStress when applying BCs.

- When using Deardorff, the interior values are the subgrid heat fluxes that go
  into the buoyancy term. These are based on resolved temperature gradients and
  modeled eddy diffusivity.

---------

Co-authored-by: Ann Almgren <[email protected]>
Co-authored-by: Aaron M. Lattanzi <[email protected]>
  • Loading branch information
@ewquon @asalmgren @AMLattanzi
3 people authored Mar 12, 2024
1 parent 78370fd commit f97cda0
Show file tree
Hide file tree
Showing 17 changed files with 703 additions and 180 deletions.
4 changes: 4 additions & 0 deletions CMake/BuildERFExe.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,10 @@ 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: 2 additions & 0 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,8 @@ 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
16 changes: 16 additions & 0 deletions Source/BoundaryConditions/ABLMost.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,10 @@ public:
if (custom_qstar != 0) {
AMREX_ASSERT_WITH_MESSAGE(use_moisture, "Specified custom MOST qv flux without moisture model!");
}
amrex::Print() << "Using specified ustar, tstar, qstar for MOST = "
<< custom_ustar << " "
<< custom_tstar << " "
<< custom_qstar << std::endl;

// Specify surface temperature or surface flux
} else {
Expand Down Expand Up @@ -216,14 +220,26 @@ 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* heat_flux,
#else
amrex::MultiFab* eddyDiffs,
#endif
amrex::MultiFab* z_phys,
const amrex::Real& dz_no_terrain,
const FluxCalc& flux_comp);
Expand Down
157 changes: 129 additions & 28 deletions Source/BoundaryConditions/ABLMost.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,19 +121,49 @@ 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 zlo = 0;
const int klo = 0;
if (flux_type == FluxCalcType::MOENG) {
moeng_flux flux_comp(zlo);
compute_most_bcs(lev, mfs, eddyDiffs, z_phys, m_geom[lev].CellSize(2), flux_comp);
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(zlo);
compute_most_bcs(lev, mfs, eddyDiffs, z_phys, m_geom[lev].CellSize(2), flux_comp);
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(zlo);
compute_most_bcs(lev, mfs, eddyDiffs, z_phys, m_geom[lev].CellSize(2), flux_comp);
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 @@ -150,24 +180,41 @@ 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)
{
const int zlo = 0;
const int klo = 0;
const int icomp = 0;
for (MFIter mfi(*mfs[0]); mfi.isValid(); ++mfi)
{
// TODO: No LSM lateral ghost cells, should this change?
// Valid CC box
Box vbx = mfi.validbox(); vbx.makeSlab(2,zlo-1);
Box vbx = mfi.validbox(); vbx.makeSlab(2,klo-1);

Box vbxx = surroundingNodes(vbx,0);
Box vbxy = surroundingNodes(vbx,1);

// 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
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
Expand Down Expand Up @@ -201,61 +248,115 @@ ABLMost::compute_most_bcs (const int& lev,

if (var_idx == Vars::cons) {
Box b2d = bx;
b2d.setBig(2,zlo-1);
b2d.setBig(2,klo-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;
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;

// This is the _kinematic_ heat flux [K m/s]
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,
cons_arr, velx_arr, vely_arr, 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


int is_land = (lmask_arr) ? lmask_arr(i,j,zlo) : 1;
// 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,zlo) = Tflux;
lsm_flux_arr(i,j,klo) = Tflux;
}
#ifdef ERF_EXPLICIT_MOST_STRESS
else if ((k == klo-1) && vbx.contains(i,j,k)) {
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,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain;
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,
cons_arr, velx_arr, vely_arr,
umm_arr, tm_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,
cons_arr, velx_arr, vely_arr,
eta_arr,
umm_arr, tm_arr, u_star_arr, q_star_arr, t_surf_arr,
dest_arr);
#endif
});
}

} else if (var_idx == Vars::xvel) {

Box xb2d = surroundingNodes(bx,0);
xb2d.setBig(2,zlo-1);
xb2d.setBig(2,klo-1);

ParallelFor(xb2d, [=] 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;
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);
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,
cons_arr, velx_arr, vely_arr,
umm_arr, um_arr, u_star_arr,
dest_arr);
if ((k == klo-1) && vbxx.contains(i,j,k)) {
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
});

} else if (var_idx == Vars::yvel) {

Box yb2d = surroundingNodes(bx,1);
yb2d.setBig(2,zlo-1);
yb2d.setBig(2,klo-1);

ParallelFor(yb2d, [=] 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;
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);
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,
cons_arr, velx_arr, vely_arr,
umm_arr, vm_arr, u_star_arr,
dest_arr);
if ((k == klo-1) && vbxy.contains(i,j,k)) {
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
});
}
} // var_idx
Expand Down
13 changes: 11 additions & 2 deletions Source/BoundaryConditions/ERF_FillPatch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -319,8 +319,17 @@ 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,eddyDiffs_lev[lev].get(),z_phys_nd[lev].get());
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());
}

// We always come in to this call with momenta so we need to leave with momenta!
// We need to make sure we convert back on all ghost cells/faces because this is
Expand Down
Loading

0 comments on commit f97cda0

Please sign in to comment.