Skip to content

Commit

Permalink
Output SGS fluxes and also correct their index type. (erf-model#1713)
Browse files Browse the repository at this point in the history 
* Output SGS fluxes and also correct their index type.

* comment unused var.

* Fixes for not using moisture cases.

* Another fix for no moisture.

* Add separate outputs for Q1 (qv) and Q2 (qc)

* Output SGS Q1,Q2 fluxes if staggered

* SGS heat, qv, and qc fluxes should be staggered quantities

* Update doc

* Fix sign and output kinematic (not dynamic) fluxes

---------

Co-authored-by: Eliot Quon <[email protected]>
  • Loading branch information
@AMLattanzi @ewquon
AMLattanzi and ewquon authored Jul 30, 2024
1 parent 7987a91 commit 71431cc
Show file tree
Hide file tree
Showing 15 changed files with 206 additions and 58 deletions.
6 changes: 5 additions & 1 deletion Docs/sphinx_doc/Inputs.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -692,7 +692,11 @@ The requested output files have the following columns:

#. SGS stress tensor component, :math:`\tau_{33}` (m2/s2)

#. SGS heat flux, :math:`\tau_{\theta w}` (K m/s)
#. *SGS heat flux*, :math:`\tau_{\theta w}` (K m/s) -- *staggered*

#. *SGS water vapor flux*, :math:`\tau_{q_v w}` (K m/s) -- *staggered*

#. *SGS cloud water flux*, :math:`\tau_{q_c w}` (K m/s) -- *staggered*

#. SGS turbulence dissipation, :math:`\epsilon` (m2/s3)

Expand Down
2 changes: 1 addition & 1 deletion Source/BoundaryConditions/ABLMost.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -317,7 +317,7 @@ ABLMost::compute_most_bcs (const int& lev,
lsm_flux_arr(i,j,klo) = Tflux;
}
else if ((k == klo-1) && vbx.contains(i,j,k) && exp_most) {
hfx_arr(i,j,klo-1) = Tflux;
hfx_arr(i,j,klo) = Tflux;
}
});

Expand Down
4 changes: 4 additions & 0 deletions Source/Diffusion/Diffusion.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,8 @@ void DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain,
const amrex::Array4<const amrex::Real>& mf_u,
const amrex::Array4<const amrex::Real>& mf_v ,
amrex::Array4< amrex::Real>& hfx_z,
amrex::Array4< amrex::Real>& qfx1_z,
amrex::Array4< amrex::Real>& qfx2_z,
amrex::Array4< amrex::Real>& diss,
const amrex::Array4<const amrex::Real>& mu_turb,
const DiffChoice& diffChoice,
Expand Down Expand Up @@ -88,6 +90,8 @@ void DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain,
const amrex::Array4<const amrex::Real>& mf_u,
const amrex::Array4<const amrex::Real>& mf_v ,
amrex::Array4< amrex::Real>& hfx_z,
amrex::Array4< amrex::Real>& qfx1_z,
amrex::Array4< amrex::Real>& qfx2_z,
amrex::Array4< amrex::Real>& diss,
const amrex::Array4<const amrex::Real>& mu_turb,
const DiffChoice& diffChoice,
Expand Down
74 changes: 67 additions & 7 deletions Source/Diffusion/DiffusionSrcForState_N.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,9 @@ using namespace amrex;
* @param[in] mf_m map factor at cell center
* @param[in] mf_u map factor at x-face
* @param[in] mf_v map factor at y-face
* @param[in] hfx_z heat flux in z-dir
* @param[inout] hfx_z heat flux in z-dir
* @param[inout] qfx1_z heat flux in z-dir
* @param[out] qfx2_z heat flux in z-dir
* @param[in] diss dissipation of TKE
* @param[in] mu_turb turbulent viscosity
* @param[in] diffChoice container of diffusion parameters
Expand Down Expand Up @@ -52,6 +54,8 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
const Array4<const Real>& mf_u,
const Array4<const Real>& mf_v,
Array4< Real>& hfx_z,
Array4< Real>& qfx1_z,
Array4< Real>& qfx2_z,
Array4< Real>& diss,
const Array4<const Real>& mu_turb,
const DiffChoice &diffChoice,
Expand Down Expand Up @@ -236,10 +240,24 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
- 3. * cell_prim(i, j, k , prim_index)
+ (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
} else if (most_on_zlo && (qty_index == RhoTheta_comp)) {
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,-1);
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,0);
} else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
zflux(i,j,k,qty_index) = -rhoFace * qfx1_z(i,j,0);
} else {
zflux(i,j,k,qty_index) = rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
}

if (qty_index == RhoTheta_comp) {
if (!most_on_zlo) {
hfx_z(i,j,k) = -zflux(i,j,k,qty_index) / rhoFace;
}
} else if (qty_index == RhoQ1_comp) {
if (!most_on_zlo) {
qfx1_z(i,j,k) = -zflux(i,j,k,qty_index) / rhoFace;
}
} else if (qty_index == RhoQ2_comp) {
qfx2_z(i,j,k) = -zflux(i,j,k,qty_index) / rhoFace;
}
});
} else if (l_turb) {
// with MolecDiffType::Constant or None
Expand Down Expand Up @@ -270,6 +288,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
const int qty_index = start_comp + n;
const int prim_index = qty_index - qty_offset;

Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
Real rhoAlpha = d_alpha_eff[prim_index];
rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index])
+ mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) );
Expand All @@ -288,11 +307,24 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
- 3. * cell_prim(i, j, k , prim_index)
+ (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
} else if (most_on_zlo && (qty_index == RhoTheta_comp)) {
Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,-1);
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,0);
} else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
zflux(i,j,k,qty_index) = -rhoFace * qfx1_z(i,j,0);
} else {
zflux(i,j,k,qty_index) = rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
}

if (qty_index == RhoTheta_comp) {
if (!most_on_zlo) {
hfx_z(i,j,k) = -zflux(i,j,k,qty_index) / rhoFace;
}
} else if (qty_index == RhoQ1_comp) {
if (!most_on_zlo) {
qfx1_z(i,j,k) = -zflux(i,j,k,qty_index) / rhoFace;
}
} else if (qty_index == RhoQ2_comp) {
qfx2_z(i,j,k) = -zflux(i,j,k,qty_index) / rhoFace;
}
});
} else if(l_consA) {
// without an LES/PBL model
Expand Down Expand Up @@ -338,10 +370,24 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
- 3. * cell_prim(i, j, k , prim_index)
+ (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
} else if (most_on_zlo && (qty_index == RhoTheta_comp)) {
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,-1);
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,0);
} else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
zflux(i,j,k,qty_index) = -rhoFace * qfx1_z(i,j,0);
} else {
zflux(i,j,k,qty_index) = rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
}

if (qty_index == RhoTheta_comp) {
if (!most_on_zlo) {
hfx_z(i,j,k) = -zflux(i,j,k,qty_index) / rhoFace;
}
} else if (qty_index == RhoQ1_comp) {
if (!most_on_zlo) {
qfx1_z(i,j,k) = -zflux(i,j,k,qty_index) / rhoFace;
}
} else if (qty_index == RhoQ2_comp) {
qfx2_z(i,j,k) = -zflux(i,j,k,qty_index) / rhoFace;
}
});
} else {
// with MolecDiffType::Constant or None
Expand Down Expand Up @@ -369,6 +415,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
const int qty_index = start_comp + n;
const int prim_index = qty_index - qty_offset;

Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
Real rhoAlpha = d_alpha_eff[prim_index];

bool ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) && k == 0);
Expand All @@ -385,11 +432,24 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
- 3. * cell_prim(i, j, k , prim_index)
+ (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv;
} else if (most_on_zlo && (qty_index == RhoTheta_comp)) {
Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) );
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,-1);
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,0);
} else if (most_on_zlo && (qty_index == RhoQ1_comp)) {
zflux(i,j,k,qty_index) = -rhoFace * qfx1_z(i,j,0);
} else {
zflux(i,j,k,qty_index) = rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv;
}

if (qty_index == RhoTheta_comp) {
if (!most_on_zlo) {
hfx_z(i,j,k) = -zflux(i,j,k,qty_index) / rhoFace;
}
} else if (qty_index == RhoQ1_comp) {
if (!most_on_zlo) {
qfx1_z(i,j,k) = -zflux(i,j,k,qty_index) / rhoFace;
}
} else if (qty_index == RhoQ2_comp) {
qfx2_z(i,j,k) = -zflux(i,j,k,qty_index) / rhoFace;
}
});
}

Expand Down
14 changes: 9 additions & 5 deletions Source/Diffusion/DiffusionSrcForState_T.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,9 @@ using namespace amrex;
* @param[in] mf_m map factor at cell center
* @param[in] mf_u map factor at x-face
* @param[in] mf_v map factor at y-face
* @param[in] hfx_z heat flux in z-dir
* @param[inout] hfx_z heat flux in z-dir
* @param[inout] qfx1_z heat flux in z-dir
* @param[out] qfx2_z heat flux in z-dir
* @param[in] diss dissipation of TKE
* @param[in] mu_turb turbulent viscosity
* @param[in] diffChoice container of diffusion parameters
Expand Down Expand Up @@ -60,6 +62,8 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
const Array4<const Real>& mf_u,
const Array4<const Real>& mf_v,
Array4< Real>& hfx_z,
Array4< Real>& /*qfx1_z*/,
Array4< Real>& /*qfx2_z*/,
Array4< Real>& diss,
const Array4<const Real>& mu_turb,
const DiffChoice &diffChoice,
Expand Down Expand Up @@ -267,7 +271,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,

if (most_on_zlo && (qty_index == RhoTheta_comp)) {
// set the exact value from MOST, don't need finite diff
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,-1);
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,0);
} else {
zflux(i,j,k,qty_index) = rhoAlpha * GradCz / met_h_zeta;
}
Expand Down Expand Up @@ -343,7 +347,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
if (most_on_zlo && (qty_index == RhoTheta_comp)) {
// set the exact value from MOST, don't need finite diff
Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) );
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,-1);
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,0);
} else {
zflux(i,j,k,qty_index) = rhoAlpha * GradCz / met_h_zeta;
}
Expand Down Expand Up @@ -414,7 +418,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,

if (most_on_zlo && (qty_index == RhoTheta_comp)) {
// set the exact value from MOST, don't need finite diff
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,-1);
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,0);
} else {
zflux(i,j,k,qty_index) = rhoAlpha * GradCz / met_h_zeta;
}
Expand Down Expand Up @@ -486,7 +490,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
if (most_on_zlo && (qty_index == RhoTheta_comp)) {
// set the exact value from MOST, don't need finite diff
Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) );
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,-1);
zflux(i,j,k,qty_index) = -rhoFace * hfx_z(i,j,0);
} else {
zflux(i,j,k,qty_index) = rhoAlpha * GradCz / met_h_zeta;
}
Expand Down
8 changes: 6 additions & 2 deletions Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -212,11 +212,13 @@ public:
void derive_stress_profiles (amrex::Gpu::HostVector<amrex::Real>& h_avg_tau11, amrex::Gpu::HostVector<amrex::Real>& h_avg_tau12,
amrex::Gpu::HostVector<amrex::Real>& h_avg_tau13, amrex::Gpu::HostVector<amrex::Real>& h_avg_tau22,
amrex::Gpu::HostVector<amrex::Real>& h_avg_tau23, amrex::Gpu::HostVector<amrex::Real>& h_avg_tau33,
amrex::Gpu::HostVector<amrex::Real>& h_avg_hfx3, amrex::Gpu::HostVector<amrex::Real>& h_avg_diss);
amrex::Gpu::HostVector<amrex::Real>& h_avg_hfx3, amrex::Gpu::HostVector<amrex::Real>& h_avg_q1fx3,
amrex::Gpu::HostVector<amrex::Real>& h_avg_q2fx3, amrex::Gpu::HostVector<amrex::Real>& h_avg_diss);
void derive_stress_profiles_stag (amrex::Gpu::HostVector<amrex::Real>& h_avg_tau11, amrex::Gpu::HostVector<amrex::Real>& h_avg_tau12,
amrex::Gpu::HostVector<amrex::Real>& h_avg_tau13, amrex::Gpu::HostVector<amrex::Real>& h_avg_tau22,
amrex::Gpu::HostVector<amrex::Real>& h_avg_tau23, amrex::Gpu::HostVector<amrex::Real>& h_avg_tau33,
amrex::Gpu::HostVector<amrex::Real>& h_avg_hfx3, amrex::Gpu::HostVector<amrex::Real>& h_avg_diss);
amrex::Gpu::HostVector<amrex::Real>& h_avg_hfx3, amrex::Gpu::HostVector<amrex::Real>& h_avg_q1fx3,
amrex::Gpu::HostVector<amrex::Real>& h_avg_q2fx3, amrex::Gpu::HostVector<amrex::Real>& h_avg_diss);

// Perform the volume-weighted sum
amrex::Real
Expand Down Expand Up @@ -719,6 +721,8 @@ private:
// Other SFS terms
amrex::Vector<std::unique_ptr<amrex::MultiFab>> SFS_hfx1_lev, SFS_hfx2_lev, SFS_hfx3_lev;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> SFS_diss_lev;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> SFS_q1fx3_lev;
amrex::Vector<std::unique_ptr<amrex::MultiFab>> SFS_q2fx3_lev;

// Terrain / grid stretching
amrex::Vector<amrex::Real> zlevels_stag; // nominal height levels
Expand Down
4 changes: 4 additions & 0 deletions Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -226,6 +226,8 @@ ERF::ERF ()
Tau23_lev.resize(nlevs_max); Tau32_lev.resize(nlevs_max);
SFS_hfx1_lev.resize(nlevs_max); SFS_hfx2_lev.resize(nlevs_max); SFS_hfx3_lev.resize(nlevs_max);
SFS_diss_lev.resize(nlevs_max);
SFS_q1fx3_lev.resize(nlevs_max);
SFS_q2fx3_lev.resize(nlevs_max);
eddyDiffs_lev.resize(nlevs_max);
SmnSmn_lev.resize(nlevs_max);

Expand Down Expand Up @@ -1933,6 +1935,8 @@ ERF::ERF (const RealBox& rb, int max_level_in,
Tau23_lev.resize(nlevs_max); Tau32_lev.resize(nlevs_max);
SFS_hfx1_lev.resize(nlevs_max); SFS_hfx2_lev.resize(nlevs_max); SFS_hfx3_lev.resize(nlevs_max);
SFS_diss_lev.resize(nlevs_max);
SFS_q1fx3_lev.resize(nlevs_max);
SFS_q2fx3_lev.resize(nlevs_max);
eddyDiffs_lev.resize(nlevs_max);
SmnSmn_lev.resize(nlevs_max);

Expand Down
16 changes: 13 additions & 3 deletions Source/ERF_make_new_arrays.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -359,6 +359,7 @@ ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMap
bool l_use_kturb = ( (solverChoice.turbChoice[lev].les_type != LESType::None) ||
(solverChoice.turbChoice[lev].pbl_type != PBLType::None) );
bool l_use_ddorf = ( solverChoice.turbChoice[lev].les_type == LESType::Deardorff);
bool l_use_moist = ( solverChoice.moisture_type != MoistureType::None );

BoxArray ba12 = convert(ba, IntVect(1,1,0));
BoxArray ba13 = convert(ba, IntVect(1,0,1));
Expand All @@ -384,14 +385,23 @@ ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMap
Tau31_lev[lev] = nullptr;
Tau32_lev[lev] = nullptr;
}
SFS_hfx1_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
SFS_hfx2_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
SFS_hfx3_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
SFS_hfx1_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(1,0,0)), dm, 1, IntVect(1,1,1) );
SFS_hfx2_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(0,1,0)), dm, 1, IntVect(1,1,1) );
SFS_hfx3_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(0,0,1)), dm, 1, IntVect(1,1,1) );
SFS_diss_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
SFS_hfx1_lev[lev]->setVal(0.);
SFS_hfx2_lev[lev]->setVal(0.);
SFS_hfx3_lev[lev]->setVal(0.);
SFS_diss_lev[lev]->setVal(0.);
if (l_use_moist) {
SFS_q1fx3_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(0,0,1)), dm, 1, IntVect(1,1,1) );
SFS_q1fx3_lev[lev]->setVal(0.0);
SFS_q2fx3_lev[lev] = std::make_unique<MultiFab>( convert(ba,IntVect(0,0,1)), dm, 1, IntVect(1,1,1) );
SFS_q2fx3_lev[lev]->setVal(0.0);
} else {
SFS_q1fx3_lev[lev] = nullptr;
SFS_q2fx3_lev[lev] = nullptr;
}
} else {
Tau11_lev[lev] = nullptr; Tau22_lev[lev] = nullptr; Tau33_lev[lev] = nullptr;
Tau12_lev[lev] = nullptr; Tau21_lev[lev] = nullptr;
Expand Down
Loading

0 comments on commit 71431cc

Please sign in to comment.