ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
SHOCInterface Class Reference

#include <ERF_ShocInterface.H>

Collaboration diagram for SHOCInterface:

Classes

struct  Buffer
 
struct  SHOCPostprocess
 
struct  SHOCPreprocess
 

Public Member Functions

 SHOCInterface (const int &lev, SolverChoice &sc)
 
void set_grids (int &level, const amrex::BoxArray &ba, amrex::Geometry &geom, amrex::MultiFab *cons, amrex::MultiFab *xvel, amrex::MultiFab *yvel, amrex::MultiFab *zvel, amrex::Real *w_subsid, amrex::MultiFab *tau13, amrex::MultiFab *tau23, amrex::MultiFab *hfx3, amrex::MultiFab *qfx3, amrex::MultiFab *eddyDiffs, amrex::MultiFab *z_phys)
 
void initialize_impl ()
 
void run_impl (const Real dt)
 
void finalize_impl ()
 
void alloc_buffers ()
 
void dealloc_buffers ()
 
void mf_to_kokkos_buffers ()
 
void kokkos_buffers_to_mf ()
 
std::string name () const
 

Protected Member Functions

void check_flux_state_consistency (const double dt)
 
void apply_turbulent_mountain_stress ()
 
void set_computed_group_impl ()
 
size_t requested_buffer_size_in_bytes () const
 
void init_buffers ()
 

Protected Attributes

amrex::Vector< int > m_col_offsets
 
Int m_num_cols = 0
 
Int m_num_layers = 0
 
Int m_npbl
 
Int m_nadv
 
Int m_num_tracers = 3
 
Int m_num_vel_comp = 2
 
Int hdtime
 
int m_lev
 
int m_step
 
amrex::Geometry m_geom
 
amrex::BoxArray m_ba
 
amrex::MultiFab * m_cons = nullptr
 
amrex::MultiFab * m_xvel = nullptr
 
amrex::MultiFab * m_yvel = nullptr
 
amrex::MultiFab * m_zvel = nullptr
 
amrex::Realm_w_subsid = nullptr
 
amrex::MultiFab * m_tau13 = nullptr
 
amrex::MultiFab * m_tau23 = nullptr
 
amrex::MultiFab * m_hfx3 = nullptr
 
amrex::MultiFab * m_qfx3 = nullptr
 
amrex::MultiFab * m_mu = nullptr
 
amrex::MultiFab * m_z_phys = nullptr
 
bool m_first_step = true
 
Buffer m_buffer
 
SHF::SHOCInput input
 
SHF::SHOCInputOutput input_output
 
SHF::SHOCOutput output
 
SHF::SHOCHistoryOutput history_output
 
SHF::SHOCRuntime runtime_options
 
SHOCPreprocess shoc_preprocess
 
SHOCPostprocess shoc_postprocess
 
ekat::WorkspaceManager< Spack, KT::Device > workspace_mgr
 
view_1d tot_buff_view
 
bool apply_tms = true
 
bool check_flux_state = false
 
bool extra_shoc_diags = false
 
bool column_conservation_check = false
 
view_2d omega
 
view_1d surf_sens_flux
 
sview_2d surf_mom_flux
 
view_1d surf_evap
 
view_2d T_mid
 
view_2d qv
 
view_1d surf_drag_coeff_tms
 
view_2d p_mid
 
view_2d p_int
 
view_2d pseudo_dens
 
view_1d phis
 
view_3d horiz_wind
 
view_2d sgs_buoy_flux
 
view_2d tk
 
view_2d cldfrac_liq
 
view_2d tke
 
view_2d qc
 
view_1d pblh
 
view_2d inv_qc_relvar
 
view_2d tkh
 
view_2d w_sec
 
view_2d cldfrac_liq_prev
 
view_1d ustar
 
view_1d obklen
 
view_2d brunt
 
view_2d shoc_mix
 
view_2d isotropy
 
view_2d shoc_cond
 
view_2d shoc_evap
 
view_2d wthl_sec
 
view_2d thl_sec
 
view_2d wqw_sec
 
view_2d qw_sec
 
view_2d uw_sec
 
view_2d vw_sec
 
view_2d w3
 
view_3d_strided qtracers
 
view_1d vapor_flux
 
view_1d water_flux
 
view_1d ice_flux
 
view_1d heat_flux
 

Private Types

using SHF = scream::shoc::Functions< Real, KokkosDefaultDevice >
 
using PF = scream::PhysicsFunctions< KokkosDefaultDevice >
 
using C = scream::physics::Constants< Real >
 
using KT = ekat::KokkosTypes< KokkosDefaultDevice >
 
using SC = scream::shoc::Constants< Real >
 
using Spack = typename SHF::Spack
 
using IntSmallPack = typename SHF::IntSmallPack
 
using Smask = typename SHF::Smask
 
using view_1d_int = typename KT::template view_1d< Int >
 
using view_1d = typename SHF::view_1d< Real >
 
using view_1d_const = typename SHF::view_1d< const Real >
 
using view_2d = typename SHF::view_2d< SHF::Spack >
 
using view_2d_const = typename SHF::view_2d< const Spack >
 
using sview_2d = typename KT::template view_2d< Real >
 
using sview_2d_const = typename KT::template view_2d< const Real >
 
using view_3d = typename SHF::view_3d< Spack >
 
using view_3d_const = typename SHF::view_3d< const Spack >
 
using view_3d_strided = typename SHF::view_3d_strided< Spack >
 
using WSM = ekat::WorkspaceManager< Spack, KT::Device >
 
template<typename ScalarT >
using uview_1d = ekat::Unmanaged< typename KT::template view_1d< ScalarT > >
 
template<typename ScalarT >
using uview_2d = ekat::Unmanaged< typename KT::template view_2d< ScalarT > >
 

Member Typedef Documentation

◆ C

using SHOCInterface::C = scream::physics::Constants<Real>
private

◆ IntSmallPack

using SHOCInterface::IntSmallPack = typename SHF::IntSmallPack
private

◆ KT

using SHOCInterface::KT = ekat::KokkosTypes<KokkosDefaultDevice>
private

◆ PF

using SHOCInterface::PF = scream::PhysicsFunctions<KokkosDefaultDevice>
private

◆ SC

using SHOCInterface::SC = scream::shoc::Constants<Real>
private

◆ SHF

using SHOCInterface::SHF = scream::shoc::Functions<Real, KokkosDefaultDevice>
private

◆ Smask

using SHOCInterface::Smask = typename SHF::Smask
private

◆ Spack

using SHOCInterface::Spack = typename SHF::Spack
private

◆ sview_2d

using SHOCInterface::sview_2d = typename KT::template view_2d<Real>
private

◆ sview_2d_const

using SHOCInterface::sview_2d_const = typename KT::template view_2d<const Real>
private

◆ uview_1d

template<typename ScalarT >
using SHOCInterface::uview_1d = ekat::Unmanaged<typename KT::template view_1d<ScalarT> >
private

◆ uview_2d

template<typename ScalarT >
using SHOCInterface::uview_2d = ekat::Unmanaged<typename KT::template view_2d<ScalarT> >
private

◆ view_1d

using SHOCInterface::view_1d = typename SHF::view_1d<Real>
private

◆ view_1d_const

using SHOCInterface::view_1d_const = typename SHF::view_1d<const Real>
private

◆ view_1d_int

using SHOCInterface::view_1d_int = typename KT::template view_1d<Int>
private

◆ view_2d

using SHOCInterface::view_2d = typename SHF::view_2d<SHF::Spack>
private

◆ view_2d_const

using SHOCInterface::view_2d_const = typename SHF::view_2d<const Spack>
private

◆ view_3d

using SHOCInterface::view_3d = typename SHF::view_3d<Spack>
private

◆ view_3d_const

using SHOCInterface::view_3d_const = typename SHF::view_3d<const Spack>
private

◆ view_3d_strided

using SHOCInterface::view_3d_strided = typename SHF::view_3d_strided<Spack>
private

◆ WSM

using SHOCInterface::WSM = ekat::WorkspaceManager<Spack, KT::Device>
private

Constructor & Destructor Documentation

◆ SHOCInterface()

SHOCInterface::SHOCInterface ( const int &  lev,
SolverChoice sc 
)
8 {
9  //
10  // Defaults set from E3SM/components/eamxx/cime_config/namelist_defaults_eamxx.xml
11  //
12  // Turn off SGS variability in SHOC, effectively reducing it to a 1.5 TKE closure?
13  bool def_shoc_1p5tke = false;
14  // Minimum value of stability correction
15  Real def_lambda_low = 0.001;
16  // Maximum value of stability correction
17  Real def_lambda_high = 0.08;
18  // Slope of change from lambda_low to lambda_high
19  Real def_lambda_slope = 0.08;
20  // stability threshold for which to apply more stability correction
21  Real def_lambda_thresh = 0.02;
22  // Temperature variance tuning factor
23  Real def_thl2tune = 1.0;
24  // Moisture variance tuning factor
25  Real def_qw2tune = 1.0;
26  // Temperature moisture covariance
27  Real def_qwthl2tune = 1.0;
28  // Vertical velocity variance
29  Real def_w2tune = 1.0;
30  // Length scale factor
31  Real def_length_fac = 0.5;
32  // Third moment vertical velocity damping factor
33  Real def_c_diag_3rd_mom = 7.0;
34  // Eddy diffusivity coefficient for heat
35  Real def_coeff_kh = 0.1;
36  // Eddy diffusivity coefficient for momentum
37  Real def_coeff_km = 0.1;
38 
39  runtime_options.lambda_low = def_lambda_low;
40  runtime_options.lambda_high = def_lambda_high;
41  runtime_options.lambda_slope = def_lambda_slope;
42  runtime_options.lambda_thresh = def_lambda_thresh;
43 
44  runtime_options.thl2tune = def_thl2tune;
45  runtime_options.qwthl2tune = def_qwthl2tune;
46  runtime_options.qw2tune = def_qw2tune;
47  runtime_options.w2tune = def_w2tune;
48 
49  runtime_options.length_fac = def_length_fac;
50  runtime_options.c_diag_3rd_mom = def_c_diag_3rd_mom;
51  runtime_options.Ckh = def_coeff_kh;
52  runtime_options.Ckm = def_coeff_km;
53  runtime_options.shoc_1p5tke = def_shoc_1p5tke;
54  runtime_options.extra_diags = extra_shoc_diags;
55 
56  // Construct parser object for following reads
57  ParmParse pp("erf.shoc");
58 
59  // Parse runtime inputs at start up
60  pp.query("lambda_low" , runtime_options.lambda_low );
61  pp.query("lambda_high" , runtime_options.lambda_high );
62  pp.query("lambda_slope" , runtime_options.lambda_slope );
63  pp.query("lambda_thresh" , runtime_options.lambda_thresh );
64  pp.query("thl2tune" , runtime_options.thl2tune );
65  pp.query("qw2tune" , runtime_options.qw2tune );
66  pp.query("qwthl2tune" , runtime_options.qwthl2tune );
67  pp.query("w2tune" , runtime_options.w2tune );
68  pp.query("length_fac" , runtime_options.length_fac );
69  pp.query("c_diag_3rd_mom" , runtime_options.c_diag_3rd_mom);
70  pp.query("coeff_kh" , runtime_options.Ckh );
71  pp.query("coeff_km" , runtime_options.Ckm );
72  pp.query("shoc_1p5tke" , runtime_options.shoc_1p5tke );
73  pp.query("extra_shoc_diags", runtime_options.extra_diags );
74 
75  // Set to default but allow us to change it through the inputs file
76  pp.query("apply_tms", apply_tms);
77  pp.query("check_flux_state", check_flux_state);
78  pp.query("extra_shoc_diags", extra_shoc_diags);
79  pp.query("column_conservation_check", column_conservation_check);
80 }
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real pp(amrex::Real y)
Definition: ERF_MicrophysicsUtils.H:233
amrex::Real Real
Definition: ERF_ShocInterface.H:19
bool check_flux_state
Definition: ERF_ShocInterface.H:625
bool extra_shoc_diags
Definition: ERF_ShocInterface.H:626
bool apply_tms
Definition: ERF_ShocInterface.H:624
SHF::SHOCRuntime runtime_options
Definition: ERF_ShocInterface.H:606
bool column_conservation_check
Definition: ERF_ShocInterface.H:627
Here is the call graph for this function:

Member Function Documentation

◆ alloc_buffers()

void SHOCInterface::alloc_buffers ( )
187 {
188  // Interface data structures
189  //=======================================================
190  omega = view_2d("Omega" , m_num_cols, m_num_layers );
191  surf_sens_flux = view_1d("Sfc sens flux" , m_num_cols);
192  surf_mom_flux = sview_2d("Sfc mom flux" , m_num_cols, m_num_vel_comp);
193  surf_evap = view_1d("Sfc evap" , m_num_cols);
194  T_mid = view_2d("T_mid" , m_num_cols, m_num_layers );
195  qv = view_2d("Qv" , m_num_cols, m_num_layers );
196  if (apply_tms) { surf_drag_coeff_tms = view_1d("surf_drag_coeff" , m_num_cols); }
197 
198  // Input data structures
199  //=======================================================
200  p_mid = view_2d("P_mid" , m_num_cols, m_num_layers );
201  p_int = view_2d("P_int" , m_num_cols, m_num_layers+1);
202  pseudo_dens = view_2d("Pseudo density" , m_num_cols, m_num_layers );
203  phis = view_1d("Phis" , m_num_cols);
204 
205  // Input/Output data structures
206  //=======================================================
208  cldfrac_liq = view_2d("Cld_frac_liq" , m_num_cols, m_num_layers );
209  tke = view_2d("Tke" , m_num_cols, m_num_layers );
210  qc = view_2d("Qc" , m_num_cols, m_num_layers );
211 
212  // Output data structures
213  //=======================================================
214  pblh = view_1d("pbl_height" , m_num_cols);
215  inv_qc_relvar = view_2d("inv_qc_relvar" , m_num_cols, m_num_layers);
216  tkh = view_2d("eddy_diff_heat" , m_num_cols, m_num_layers);
217  w_sec = view_2d("w_sec" , m_num_cols, m_num_layers);
218  cldfrac_liq_prev = view_2d("cld_frac_liq_prev" , m_num_cols, m_num_layers );
219  ustar = view_1d("ustar" , m_num_cols);
220  obklen = view_1d("obklen" , m_num_cols);
221 
222  // Extra diagnostic data structures
223  //=======================================================
224  if (extra_shoc_diags) {
225  brunt = view_2d("brunt" , m_num_cols, m_num_layers);
226  shoc_mix = view_2d("shoc_mix" , m_num_cols, m_num_layers);
227  isotropy = view_2d("isotropy" , m_num_cols, m_num_layers);
228  shoc_cond = view_2d("shoc_cond", m_num_cols, m_num_layers);
229  shoc_evap = view_2d("shoc_evap", m_num_cols, m_num_layers);
230 
231  wthl_sec = view_2d("wthl_sec" , m_num_cols, m_num_layers+1);
232  thl_sec = view_2d("thl_sec" , m_num_cols, m_num_layers+1);
233  wqw_sec = view_2d("wqw_sec" , m_num_cols, m_num_layers+1);
234  qw_sec = view_2d("qw_sec" , m_num_cols, m_num_layers+1);
235  uw_sec = view_2d("uw_sec" , m_num_cols, m_num_layers+1);
236  vw_sec = view_2d("vw_sec" , m_num_cols, m_num_layers+1);
237  w3 = view_2d("w3" , m_num_cols, m_num_layers+1);
238  }
239 
240  // Tracer data structures
241  //=======================================================
242  // NOTE: Use layoutright format
243  Kokkos::LayoutStride layout(m_num_cols , m_num_cols*m_num_layers, // stride for dim0
244  m_num_layers , m_num_layers, // stride for dim1
245  m_num_tracers, 1 // stride for dim2
246  );
247  qtracers = view_3d_strided("Qtracers" , layout);
248 
249  // Boundary flux data structures
250  //=======================================================
252  vapor_flux = view_1d("vapor_flux", m_num_cols);
253  water_flux = view_1d("water_flux", m_num_cols);
254  ice_flux = view_1d("ice_flux" , m_num_cols);
255  heat_flux = view_1d("heat_flux" , m_num_cols);
256  }
257 }
view_2d p_int
Definition: ERF_ShocInterface.H:642
view_2d wqw_sec
Definition: ERF_ShocInterface.H:675
view_1d surf_drag_coeff_tms
Definition: ERF_ShocInterface.H:637
view_2d thl_sec
Definition: ERF_ShocInterface.H:674
view_2d tkh
Definition: ERF_ShocInterface.H:659
view_2d cldfrac_liq
Definition: ERF_ShocInterface.H:651
typename KT::template view_2d< Real > sview_2d
Definition: ERF_ShocInterface.H:51
Int m_num_layers
Definition: ERF_ShocInterface.H:553
view_2d vw_sec
Definition: ERF_ShocInterface.H:678
view_1d ice_flux
Definition: ERF_ShocInterface.H:689
view_2d shoc_cond
Definition: ERF_ShocInterface.H:670
view_2d brunt
Definition: ERF_ShocInterface.H:667
view_1d heat_flux
Definition: ERF_ShocInterface.H:690
Int m_num_vel_comp
Definition: ERF_ShocInterface.H:557
view_1d surf_sens_flux
Definition: ERF_ShocInterface.H:632
view_1d pblh
Definition: ERF_ShocInterface.H:657
Int m_num_cols
Definition: ERF_ShocInterface.H:552
view_2d p_mid
Definition: ERF_ShocInterface.H:641
view_2d uw_sec
Definition: ERF_ShocInterface.H:677
view_2d qc
Definition: ERF_ShocInterface.H:653
view_2d shoc_mix
Definition: ERF_ShocInterface.H:668
view_2d qw_sec
Definition: ERF_ShocInterface.H:676
view_2d cldfrac_liq_prev
Definition: ERF_ShocInterface.H:661
view_2d omega
Definition: ERF_ShocInterface.H:631
typename SHF::view_3d_strided< Spack > view_3d_strided
Definition: ERF_ShocInterface.H:55
view_1d vapor_flux
Definition: ERF_ShocInterface.H:687
view_1d obklen
Definition: ERF_ShocInterface.H:663
view_2d pseudo_dens
Definition: ERF_ShocInterface.H:643
view_2d inv_qc_relvar
Definition: ERF_ShocInterface.H:658
view_1d surf_evap
Definition: ERF_ShocInterface.H:634
view_2d w_sec
Definition: ERF_ShocInterface.H:660
view_1d phis
Definition: ERF_ShocInterface.H:644
view_2d shoc_evap
Definition: ERF_ShocInterface.H:671
view_3d horiz_wind
Definition: ERF_ShocInterface.H:648
typename SHF::view_1d< Real > view_1d
Definition: ERF_ShocInterface.H:47
view_3d_strided qtracers
Definition: ERF_ShocInterface.H:683
view_1d ustar
Definition: ERF_ShocInterface.H:662
view_2d w3
Definition: ERF_ShocInterface.H:679
Int m_num_tracers
Definition: ERF_ShocInterface.H:556
sview_2d surf_mom_flux
Definition: ERF_ShocInterface.H:633
typename SHF::view_2d< SHF::Spack > view_2d
Definition: ERF_ShocInterface.H:49
view_2d qv
Definition: ERF_ShocInterface.H:636
view_2d wthl_sec
Definition: ERF_ShocInterface.H:673
view_2d tke
Definition: ERF_ShocInterface.H:652
view_1d water_flux
Definition: ERF_ShocInterface.H:688
view_2d T_mid
Definition: ERF_ShocInterface.H:635
typename SHF::view_3d< Spack > view_3d
Definition: ERF_ShocInterface.H:53
view_2d isotropy
Definition: ERF_ShocInterface.H:669

◆ apply_turbulent_mountain_stress()

void SHOCInterface::apply_turbulent_mountain_stress ( )
protected
935 {
936  auto rrho_i = m_buffer.rrho_i;
937  auto upwp_sfc = m_buffer.upwp_sfc;
938  auto vpwp_sfc = m_buffer.vpwp_sfc;
939 
940  const int nlev_v = (m_num_layers-1)/Spack::n;
941  const int nlev_p = (m_num_layers-1)%Spack::n;
942  const int nlevi_v = m_num_layers/Spack::n;
943  const int nlevi_p = m_num_layers%Spack::n;
944 
945  Kokkos::parallel_for("apply_tms", KT::RangePolicy(0, m_num_cols), KOKKOS_LAMBDA (const int i)
946  {
947  upwp_sfc(i) -= surf_drag_coeff_tms(i)*horiz_wind(i,0,nlev_v)[nlev_p]/rrho_i(i,nlevi_v)[nlevi_p];
948  vpwp_sfc(i) -= surf_drag_coeff_tms(i)*horiz_wind(i,1,nlev_v)[nlev_p]/rrho_i(i,nlevi_v)[nlevi_p];
949  });
950 }
Buffer m_buffer
Definition: ERF_ShocInterface.H:599
uview_1d< Real > upwp_sfc
Definition: ERF_ShocInterface.H:458
uview_1d< Real > vpwp_sfc
Definition: ERF_ShocInterface.H:459
uview_2d< Spack > rrho_i
Definition: ERF_ShocInterface.H:479

◆ check_flux_state_consistency()

void SHOCInterface::check_flux_state_consistency ( const double  dt)
protected
955 {
956  using PC = scream::physics::Constants<Real>;
957  using RU = ekat::ReductionUtils<KT::ExeSpace>;
958  using TPF = ekat::TeamPolicyFactory<KT::ExeSpace>;
959 
960  const Real gravit = PC::gravit;
961  const Real qmin = 1e-12; // minimum permitted constituent concentration (kg/kg)
962 
963  const auto& pseudo_density = pseudo_dens;
964 
965  const auto nlevs = m_num_layers;
966  const auto nlev_packs = ekat::npack<Spack>(nlevs);
967  const auto last_pack_idx = (nlevs-1)/Spack::n;
968  const auto last_pack_entry = (nlevs-1)%Spack::n;
969  const auto policy = TPF::get_default_team_policy(m_num_cols, nlev_packs);
970  Kokkos::parallel_for("check_flux_state_consistency",
971  policy,
972  KOKKOS_LAMBDA (const KT::MemberType& team)
973  {
974  const auto i = team.league_rank();
975 
976  const auto& pseudo_density_i = ekat::subview(pseudo_density, i);
977  const auto& qv_i = ekat::subview(qv, i);
978 
979  // reciprocal of pseudo_density at the bottom layer
980  const auto rpdel = 1.0/pseudo_density_i(last_pack_idx)[last_pack_entry];
981 
982  // Check if the negative surface latent heat flux can exhaust
983  // the moisture in the lowest model level. If so, apply fixer.
984  const auto condition = surf_evap(i) - (qmin - qv_i(last_pack_idx)[last_pack_entry])/(dt*gravit*rpdel);
985  if (condition < 0) {
986  const auto cc = std::fabs(surf_evap(i)*dt*gravit);
987 
988  auto tracer_mass = [&](const int k)
989  {
990  return qv_i(k)*pseudo_density_i(k);
991  };
992  Real mm = RU::view_reduction(team, 0, nlevs, tracer_mass);
993 
994  EKAT_KERNEL_ASSERT_MSG(mm >= cc, "Error! Total mass of column vapor should be greater than mass of surf_evap.\n");
995 
996  Kokkos::parallel_for(Kokkos::TeamVectorRange(team, nlev_packs), [&](const int& k)
997  {
998  const auto adjust = cc*qv_i(k)*pseudo_density_i(k)/mm;
999  qv_i(k) = (qv_i(k)*pseudo_density_i(k) - adjust)/pseudo_density_i(k);
1000  });
1001 
1002  surf_evap(i) = 0;
1003  }
1004  });
1005 }
@ PC
Definition: ERF_IndexDefines.H:123

◆ dealloc_buffers()

void SHOCInterface::dealloc_buffers ( )
262 {
263  // Contiguous memory buffer view
264  //=======================================================
266 
267  // Interface data structures
268  //=======================================================
269  omega = view_2d();
272  surf_evap = view_1d();
273  T_mid = view_2d();
274  qv = view_2d();
276 
277  // Input data structures
278  //=======================================================
279  p_mid = view_2d();
280  p_int = view_2d();
281  pseudo_dens = view_2d();
282  phis = view_1d();
283 
284  // Input/Output data structures
285  //=======================================================
286  horiz_wind = view_3d();
287  cldfrac_liq = view_2d();
288  tke = view_2d();
289  qc = view_2d();
290 
291  // Output data structures
292  //=======================================================
293  pblh = view_1d();
295  tkh = view_2d();
296  w_sec = view_2d();
298  ustar = view_1d();
299  obklen = view_1d();
300 
301  // Extra diagnostic data structures
302  //=======================================================
303  if (extra_shoc_diags) {
304  brunt = view_2d();
305  shoc_mix = view_2d();
306  isotropy = view_2d();
307  shoc_cond = view_2d();
308  shoc_evap = view_2d();
309 
310  wthl_sec = view_2d();
311  thl_sec = view_2d();
312  wqw_sec = view_2d();
313  qw_sec = view_2d();
314  uw_sec = view_2d();
315  vw_sec = view_2d();
316  w3 = view_2d();
317  }
318 
319  // Tracer data structures
320  //=======================================================
322 
323  // Boundary flux data structures
324  //=======================================================
326  vapor_flux = view_1d();
327  water_flux = view_1d();
328  ice_flux = view_1d();
329  heat_flux = view_1d();
330  }
331 }
view_1d tot_buff_view
Definition: ERF_ShocInterface.H:620

◆ finalize_impl()

void SHOCInterface::finalize_impl ( )
922 {
923  // Do nothing (per SHOCMacrophysics::finalize_impl())
924 
925  // Fill the AMReX MFs from Kokkos Views
927 
928  // Deallocate the buffer arrays
929  dealloc_buffers();
930 }
void kokkos_buffers_to_mf()
Definition: ERF_ShocInterface.cpp:485
void dealloc_buffers()
Definition: ERF_ShocInterface.cpp:261

◆ init_buffers()

void SHOCInterface::init_buffers ( )
protected
618 {
619 
620  // Buffer of contiguous memory
621  auto buffer_size = requested_buffer_size_in_bytes();
622  tot_buff_view = view_1d("contiguous shoc_buffer",buffer_size);
623  Real* mem = reinterpret_cast<Real*>(tot_buff_view.data());
624 
625  // 1d scalar views
626  using scalar_view_t = decltype(m_buffer.wpthlp_sfc);
627  scalar_view_t* _1d_scalar_view_ptrs[Buffer::num_1d_scalar_ncol] =
629 #ifdef SCREAM_SHOC_SMALL_KERNELS
630  , &m_buffer.se_b, &m_buffer.ke_b, &m_buffer.wv_b, &m_buffer.wl_b
631  , &m_buffer.se_a, &m_buffer.ke_a, &m_buffer.wv_a, &m_buffer.wl_a
632  , &m_buffer.kbfs, &m_buffer.ustar2, &m_buffer.wstar
633 #endif
634  };
635  for (int i = 0; i < Buffer::num_1d_scalar_ncol; ++i) {
636  *_1d_scalar_view_ptrs[i] = scalar_view_t(mem, m_num_cols);
637  mem += _1d_scalar_view_ptrs[i]->size();
638  }
639 
640  Spack* s_mem = reinterpret_cast<Spack*>(mem);
641 
642  // 2d packed views
643  const int nlev_packs = ekat::npack<Spack>(m_num_layers);
644  const int nlevi_packs = ekat::npack<Spack>(m_num_layers+1);
645  const int num_tracer_packs = ekat::npack<Spack>(m_num_tracers);
646 
647  m_buffer.pref_mid = decltype(m_buffer.pref_mid)(s_mem, nlev_packs);
648  s_mem += m_buffer.pref_mid.size();
649 
650  using spack_2d_view_t = decltype(m_buffer.z_mid);
651  spack_2d_view_t* _2d_spack_mid_view_ptrs[Buffer::num_2d_vector_mid] =
655 #ifdef SCREAM_SHOC_SMALL_KERNELS
656  , &m_buffer.rho_zt, &m_buffer.shoc_qv, &m_buffer.tabs, &m_buffer.dz_zt
657 #endif
658  };
659 
660  spack_2d_view_t* _2d_spack_int_view_ptrs[Buffer::num_2d_vector_int] =
664 #ifdef SCREAM_SHOC_SMALL_KERNELS
665  , &m_buffer.dz_zi
666 #endif
667  };
668 
669  for (int i = 0; i < Buffer::num_2d_vector_mid; ++i) {
670  *_2d_spack_mid_view_ptrs[i] = spack_2d_view_t(s_mem, m_num_cols, nlev_packs);
671  s_mem += _2d_spack_mid_view_ptrs[i]->size();
672  }
673 
674  for (int i = 0; i < Buffer::num_2d_vector_int; ++i) {
675  *_2d_spack_int_view_ptrs[i] = spack_2d_view_t(s_mem, m_num_cols, nlevi_packs);
676  s_mem += _2d_spack_int_view_ptrs[i]->size();
677  }
678  m_buffer.wtracer_sfc = decltype(m_buffer.wtracer_sfc)(s_mem, m_num_cols, num_tracer_packs);
679  s_mem += m_buffer.wtracer_sfc.size();
680 
681  // WSM data
682  m_buffer.wsm_data = s_mem;
683 }
size_t requested_buffer_size_in_bytes() const
Definition: ERF_ShocInterface.cpp:592
typename SHF::Spack Spack
Definition: ERF_ShocInterface.H:43
uview_2d< Spack > w3
Definition: ERF_ShocInterface.H:504
uview_2d< Spack > qc_copy
Definition: ERF_ShocInterface.H:491
uview_2d< Spack > dse
Definition: ERF_ShocInterface.H:489
uview_2d< Spack > vw_sec
Definition: ERF_ShocInterface.H:503
uview_1d< Real > wpthlp_sfc
Definition: ERF_ShocInterface.H:456
uview_2d< Spack > thv
Definition: ERF_ShocInterface.H:480
uview_2d< Spack > wqls_sec
Definition: ERF_ShocInterface.H:505
uview_2d< Spack > wthl_sec
Definition: ERF_ShocInterface.H:499
uview_2d< Spack > thlm
Definition: ERF_ShocInterface.H:487
uview_2d< Spack > uw_sec
Definition: ERF_ShocInterface.H:502
uview_2d< Spack > z_mid
Definition: ERF_ShocInterface.H:476
uview_2d< Spack > rrho
Definition: ERF_ShocInterface.H:478
uview_2d< Spack > wtke_sec
Definition: ERF_ShocInterface.H:501
uview_2d< Spack > unused
Definition: ERF_ShocInterface.H:475
Spack * wsm_data
Definition: ERF_ShocInterface.H:516
static constexpr int num_2d_vector_int
Definition: ERF_ShocInterface.H:449
uview_2d< Spack > wtracer_sfc
Definition: ERF_ShocInterface.H:484
uview_2d< Spack > inv_exner
Definition: ERF_ShocInterface.H:486
uview_2d< Spack > qwthl_sec
Definition: ERF_ShocInterface.H:498
uview_2d< Spack > qw
Definition: ERF_ShocInterface.H:488
uview_2d< Spack > brunt
Definition: ERF_ShocInterface.H:506
uview_2d< Spack > shoc_mix
Definition: ERF_ShocInterface.H:493
uview_2d< Spack > tke_copy
Definition: ERF_ShocInterface.H:490
static constexpr int num_1d_scalar_ncol
Definition: ERF_ShocInterface.H:442
uview_2d< Spack > isotropy
Definition: ERF_ShocInterface.H:494
uview_2d< Spack > zt_grid
Definition: ERF_ShocInterface.H:482
uview_2d< Spack > w_sec
Definition: ERF_ShocInterface.H:495
uview_2d< Spack > thl_sec
Definition: ERF_ShocInterface.H:496
uview_2d< Spack > wm_zt
Definition: ERF_ShocInterface.H:485
uview_2d< Spack > qw_sec
Definition: ERF_ShocInterface.H:497
uview_2d< Spack > wqw_sec
Definition: ERF_ShocInterface.H:500
uview_2d< Spack > z_int
Definition: ERF_ShocInterface.H:477
uview_2d< Spack > shoc_ql2
Definition: ERF_ShocInterface.H:492
uview_1d< Real > wprtp_sfc
Definition: ERF_ShocInterface.H:457
uview_2d< Spack > zi_grid
Definition: ERF_ShocInterface.H:483
static constexpr int num_2d_vector_mid
Definition: ERF_ShocInterface.H:448
uview_1d< Spack > pref_mid
Definition: ERF_ShocInterface.H:473
uview_2d< Spack > dz
Definition: ERF_ShocInterface.H:481

◆ initialize_impl()

void SHOCInterface::initialize_impl ( )
688 {
689  // Alias local variables from temporary buffer
690  auto z_mid = m_buffer.z_mid;
691  auto z_int = m_buffer.z_int;
692  auto wpthlp_sfc = m_buffer.wpthlp_sfc;
693  auto wprtp_sfc = m_buffer.wprtp_sfc;
694  auto upwp_sfc = m_buffer.upwp_sfc;
695  auto vpwp_sfc = m_buffer.vpwp_sfc;
696  auto rrho = m_buffer.rrho;
697  auto rrho_i = m_buffer.rrho_i;
698  auto thv = m_buffer.thv;
699  auto dz = m_buffer.dz;
700  auto zt_grid = m_buffer.zt_grid;
701  auto zi_grid = m_buffer.zi_grid;
702  auto wtracer_sfc = m_buffer.wtracer_sfc;
703  auto wm_zt = m_buffer.wm_zt;
704  auto inv_exner = m_buffer.inv_exner;
705  auto thlm = m_buffer.thlm;
706  auto qw = m_buffer.qw;
707  auto dse = m_buffer.dse;
708  auto tke_copy = m_buffer.tke_copy;
709  auto qc_copy = m_buffer.qc_copy;
710  auto shoc_ql2 = m_buffer.shoc_ql2;
711 
712  // For now, set z_int(i,nlevs) = z_surf = 0
713  const Real z_surf = 0.0;
714 
715  // Set preprocess variables
718  surf_mom_flux, qtracers, qv, qc, qc_copy, tke, tke_copy, z_mid, z_int,
719  dse, rrho, rrho_i, thv, dz, zt_grid, zi_grid, wpthlp_sfc, wprtp_sfc, upwp_sfc,
720  vpwp_sfc, wtracer_sfc, wm_zt, inv_exner, thlm, qw, cldfrac_liq, cldfrac_liq_prev);
721 
722  // Input Variables:
723  input.zt_grid = shoc_preprocess.zt_grid;
724  input.zi_grid = shoc_preprocess.zi_grid;
725  input.pres = p_mid;
726  input.presi = p_int;
727  input.pdel = pseudo_dens;
728  input.thv = shoc_preprocess.thv;
729  input.w_field = shoc_preprocess.wm_zt;
730  input.wthl_sfc = shoc_preprocess.wpthlp_sfc;
731  input.wqw_sfc = shoc_preprocess.wprtp_sfc;
732  input.uw_sfc = shoc_preprocess.upwp_sfc;
733  input.vw_sfc = shoc_preprocess.vpwp_sfc;
734  input.wtracer_sfc = shoc_preprocess.wtracer_sfc;
735  input.inv_exner = shoc_preprocess.inv_exner;
736  input.phis = phis;
737 
738  // Input/Output Variables
743  input_output.horiz_wind = horiz_wind;
744  input_output.wthv_sec = sgs_buoy_flux;
746  input_output.tk = tk;
747  input_output.shoc_cldfrac = cldfrac_liq;
748  input_output.shoc_ql = qc_copy;
749 
750  // Output Variables
751  output.pblh = pblh;
752  output.shoc_ql2 = shoc_ql2;
753  output.tkh = tkh;
754  output.ustar = ustar;
755  output.obklen = obklen;
756 
757  // Output (diagnostic)
758  history_output.shoc_mix = m_buffer.shoc_mix;
759  history_output.isotropy = m_buffer.isotropy;
760  if (extra_shoc_diags) {
761  history_output.shoc_cond = shoc_cond;
762  history_output.shoc_evap = shoc_evap;
763  } else {
764  history_output.shoc_cond = m_buffer.unused;
765  history_output.shoc_evap = m_buffer.unused;
766  }
767  history_output.w_sec = w_sec;
768  history_output.thl_sec = m_buffer.thl_sec;
769  history_output.qw_sec = m_buffer.qw_sec;
770  history_output.qwthl_sec = m_buffer.qwthl_sec;
771  history_output.wthl_sec = m_buffer.wthl_sec;
772  history_output.wqw_sec = m_buffer.wqw_sec;
773  history_output.wtke_sec = m_buffer.wtke_sec;
774  history_output.uw_sec = m_buffer.uw_sec;
775  history_output.vw_sec = m_buffer.vw_sec;
777  history_output.wqls_sec = m_buffer.wqls_sec;
778  history_output.brunt = m_buffer.brunt;
779 
780 #ifdef SCREAM_SHOC_SMALL_KERNELS
781  temporaries.se_b = m_buffer.se_b;
782  temporaries.ke_b = m_buffer.ke_b;
783  temporaries.wv_b = m_buffer.wv_b;
784  temporaries.wl_b = m_buffer.wl_b;
785  temporaries.se_a = m_buffer.se_a;
786  temporaries.ke_a = m_buffer.ke_a;
787  temporaries.wv_a = m_buffer.wv_a;
788  temporaries.wl_a = m_buffer.wl_a;
789  temporaries.kbfs = m_buffer.kbfs;
790  temporaries.ustar2 = m_buffer.ustar2;
791  temporaries.wstar = m_buffer.wstar;
792 
793  temporaries.rho_zt = m_buffer.rho_zt;
794  temporaries.shoc_qv = m_buffer.shoc_qv;
795  temporaries.tabs = m_buffer.tabs;
796  temporaries.dz_zt = m_buffer.dz_zt;
797  temporaries.dz_zi = m_buffer.dz_zi;
798 #endif
799 
800  // Set postprocess variables
802  rrho, qv, qw, qc, qc_copy, tke,
803  tke_copy, qtracers, shoc_ql2,
805  T_mid, dse, z_mid, phis);
806 
807  // Setup WSM for internal local variables
808  using TPF = ekat::TeamPolicyFactory<KT::ExeSpace>;
809  const auto nlev_packs = ekat::npack<Spack>(m_num_layers);
810  const auto nlevi_packs = ekat::npack<Spack>(m_num_layers+1);
811  const int n_wind_slots = ekat::npack<Spack>(m_num_vel_comp)*Spack::n;
812  const int n_trac_slots = ekat::npack<Spack>(m_num_tracers)*Spack::n;
813  const auto default_policy = TPF::get_default_team_policy(m_num_cols, nlev_packs);
814  workspace_mgr.setup(m_buffer.wsm_data, nlevi_packs, 14+(n_wind_slots+n_trac_slots), default_policy);
815 
816  // NOTE: FIX ME!
817  // Maximum number of levels in pbl from surface
818  const int ntop_shoc = 0;
819  const int nbot_shoc = m_num_layers;
820  view_1d pref_mid("pref_mid", m_num_layers);
821  Spack* s_mem = reinterpret_cast<Spack*>(pref_mid.data());
822  SHF::view_1d<Spack> pref_mid_um(s_mem, m_num_layers);
823  const auto policy = TPF::get_default_team_policy(m_num_cols, nlev_packs);
824  Kokkos::parallel_for("pref_mid",
825  policy,
826  KOKKOS_LAMBDA (const KT::MemberType& team)
827  {
828  const auto i = team.league_rank();
829  if (i==0) {
830  const auto& pmid_i = ekat::subview(p_mid, i);
831  Kokkos::parallel_for(Kokkos::TeamVectorRange(team, nlev_packs), [&](const int& k)
832  {
833  pref_mid_um(k) = pmid_i(k);
834  });
835  }
836  });
837  Kokkos::fence();
838  m_npbl = SHF::shoc_init(nbot_shoc,ntop_shoc,pref_mid_um);
839 
840  // Cell length for input dx and dy
841  view_1d cell_length_x("cell_length_x", m_num_cols);
842  view_1d cell_length_y("cell_length_y", m_num_cols);
843  Kokkos::deep_copy(cell_length_x, m_geom.CellSize(0));
844  Kokkos::deep_copy(cell_length_y, m_geom.CellSize(1));
845  input.dx = cell_length_x;
846  input.dy = cell_length_y;
847 }
SHOCPreprocess shoc_preprocess
Definition: ERF_ShocInterface.H:612
amrex::Geometry m_geom
Definition: ERF_ShocInterface.H:567
view_2d sgs_buoy_flux
Definition: ERF_ShocInterface.H:649
SHF::SHOCOutput output
Definition: ERF_ShocInterface.H:604
ekat::WorkspaceManager< Spack, KT::Device > workspace_mgr
Definition: ERF_ShocInterface.H:616
Int m_npbl
Definition: ERF_ShocInterface.H:554
SHF::SHOCInput input
Definition: ERF_ShocInterface.H:602
SHOCPostprocess shoc_postprocess
Definition: ERF_ShocInterface.H:613
SHF::SHOCInputOutput input_output
Definition: ERF_ShocInterface.H:603
view_2d tk
Definition: ERF_ShocInterface.H:650
SHF::SHOCHistoryOutput history_output
Definition: ERF_ShocInterface.H:605
void set_variables(const int ncol_, const int nlev_, const view_2d_const &rrho_, const view_2d &qv_, const view_2d_const &qw_, const view_2d &qc_, const view_2d_const &qc_copy_, const view_2d &tke_, const view_2d_const &tke_copy_, const view_3d_strided &qtracers_, const view_2d_const &qc2_, const view_2d &cldfrac_liq_, const view_2d &inv_qc_relvar_, const view_2d &T_mid_, const view_2d_const &dse_, const view_2d_const &z_mid_, const view_1d_const phis_)
Definition: ERF_ShocInterface.H:398
view_3d_strided qtracers
Definition: ERF_ShocInterface.H:234
view_2d thlm
Definition: ERF_ShocInterface.H:256
view_2d zt_grid
Definition: ERF_ShocInterface.H:247
view_2d inv_exner
Definition: ERF_ShocInterface.H:255
view_2d wm_zt
Definition: ERF_ShocInterface.H:254
view_2d wtracer_sfc
Definition: ERF_ShocInterface.H:253
view_1d wprtp_sfc
Definition: ERF_ShocInterface.H:250
view_2d qw
Definition: ERF_ShocInterface.H:257
view_1d upwp_sfc
Definition: ERF_ShocInterface.H:251
view_1d wpthlp_sfc
Definition: ERF_ShocInterface.H:249
view_1d vpwp_sfc
Definition: ERF_ShocInterface.H:252
view_2d shoc_s
Definition: ERF_ShocInterface.H:240
view_2d tke_copy
Definition: ERF_ShocInterface.H:242
view_2d zi_grid
Definition: ERF_ShocInterface.H:248
view_2d thv
Definition: ERF_ShocInterface.H:245
void set_variables(const int ncol_, const int nlev_, const Real z_surf_, const view_2d_const &T_mid_, const view_2d_const &p_mid_, const view_2d_const &p_int_, const view_2d_const &pseudo_density_, const view_2d_const &omega_, const view_1d_const &phis_, const view_1d_const &surf_sens_flux_, const view_1d_const &surf_evap_, const sview_2d_const &surf_mom_flux_, const view_3d_strided &qtracers_, const view_2d &qv_, const view_2d_const &qc_, const view_2d &qc_copy_, const view_2d &tke_, const view_2d &tke_copy_, const view_2d &z_mid_, const view_2d &z_int_, const view_2d &dse_, const view_2d &rrho_, const view_2d &rrho_i_, const view_2d &thv_, const view_2d &dz_, const view_2d &zt_grid_, const view_2d &zi_grid_, const view_1d &wpthlp_sfc_, const view_1d &wprtp_sfc_, const view_1d &upwp_sfc_, const view_1d &vpwp_sfc_, const view_2d &wtracer_sfc_, const view_2d &wm_zt_, const view_2d &inv_exner_, const view_2d &thlm_, const view_2d &qw_, const view_2d &cldfrac_liq_, const view_2d &cldfrac_liq_prev_)
Definition: ERF_ShocInterface.H:263

◆ kokkos_buffers_to_mf()

void SHOCInterface::kokkos_buffers_to_mf ( )
486 {
487  //
488  // Expose for device capture
489  //
490 
491  // Interface data structures
492  //=======================================================
493  auto T_mid_d = T_mid;
494  auto qv_d = qv;
495 
496  // Input data structures
497  //=======================================================
498  auto p_mid_d = p_mid;
499 
500  // Input/Output data structures
501  //=======================================================
502  auto horiz_wind_d = horiz_wind;
503  //auto sgs_buoy_flux_d = sgs_buoy_flux;
504  //auto tk_d = tk;
505  //auto cldfrac_liq_d = cldfrac_liq;
506  auto tke_d = tke;
507  auto qc_d = qc;
508 
509  // Output data structures
510  //=======================================================
511  //auto pblh_d = pblh;
512  //auto inv_qc_relvar_d = inv_qc_relvar;
513  //auto tkh_d = tkh;
514  //auto w_sec_d = w_sec;
515  //auto cldfrac_liq_prev_d = cldfrac_liq_prev;
516  //auto ustar_d = ustar;
517  //auto obklen_d = obklen;
518 
519  bool moist = (m_cons->nComp() > RhoQ1_comp);
520  for (MFIter mfi(*m_cons); mfi.isValid(); ++mfi) {
521  const auto& vbx = mfi.validbox();
522  const int nx = vbx.length(0);
523  const int imin = vbx.smallEnd(0);
524  const int jmin = vbx.smallEnd(1);
525  const int offset = m_col_offsets[mfi.index()];
526 
527  const Array4<Real>& cons_arr = m_cons->array(mfi);
528 
529  const Array4<Real>& u_arr = m_xvel->array(mfi);
530  const Array4<Real>& v_arr = m_yvel->array(mfi);
531 
532  const Array4<Real>& mu_arr = m_mu->array(mfi);
533 
534  int ilo = vbx.smallEnd(0);
535  int ihi = vbx.bigEnd(0);
536  int jlo = vbx.smallEnd(1);
537  int jhi = vbx.bigEnd(1);
538 
539  ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
540  {
541  // map [i,j,k] 0-based to [icol, ilay] 0-based
542  const int icol = (j-jmin)*nx + (i-imin) + offset;
543  const int ilay = k;
544 
545  // Density at CC
546  Real r = cons_arr(i,j,k,Rho_comp);
547 
548  // Theta at CC
549  Real Th = getThgivenTandP(T_mid_d(icol,ilay)[0], p_mid_d(icol,ilay)[0], R_d/Cp_d);
550 
551  // Interface data structures
552  //=======================================================
553  if (moist) { cons_arr(i,j,k,RhoQ1_comp) = r * qv_d(icol,ilay)[0]; }
554  cons_arr(i,j,k,RhoTheta_comp) = r * Th;
555 
556  // Input/Output data structures
557  //=======================================================
558  // NOTE: Averaging for U/V to account for the CC data
559  int icolim = (j -jmin)*nx + (i-1-imin) + offset;
560  int icoljm = (j-1-jmin)*nx + (i -imin) + offset;
561  if (i==ilo) {
562  int icolip = (j -jmin)*nx + (i+1-imin) + offset;
563  u_arr(i,j,k) = 1.5 * horiz_wind_d(icol,0,ilay)[0] - 0.5 * horiz_wind_d(icolip,0,ilay)[0];
564  } else {
565  u_arr(i,j,k) = 0.5 * (horiz_wind_d(icol,0,ilay)[0] + horiz_wind_d(icolim,0,ilay)[0]);
566  if (i==ihi) {
567  u_arr(i+1,j,k) = 1.5 * horiz_wind_d(icol,0,ilay)[0] - 0.5 * horiz_wind_d(icolim,0,ilay)[0];
568  }
569  }
570  if (j==ilo) {
571  int icoljp = (j+1-jmin)*nx + (i -imin) + offset;
572  v_arr(i,j,k) = 1.5 * horiz_wind_d(icol,1,ilay)[0] - 0.5 * horiz_wind_d(icoljp,1,ilay)[0];
573  } else {
574  v_arr(i,j,k) = 0.5 * (horiz_wind_d(icol,1,ilay)[0] + horiz_wind_d(icoljm,1,ilay)[0]);
575  if (j==jhi) {
576  v_arr(i,j+1,k) = 1.5 * horiz_wind_d(icol,1,ilay)[0] - 0.5 * horiz_wind_d(icoljm,1,ilay)[0];
577  }
578  }
579  //mu_arr(i,j,k,EddyDiff::Mom_v) = tk_d(icol,ilay)[0];
580  cons_arr(i,j,k,RhoKE_comp) = std::max(r*tke_d(icol,ilay)[0],0.0);
581  if (moist) { cons_arr(i,j,k,RhoQ2_comp) = r * qc_d(icol,ilay)[0]; }
582 
583  // Output data structures
584  //=======================================================
585  //mu_arr(i,j,k,EddyDiff::Theta_v) = tkh_d(icol,ilay)[0];
586  });
587  }
588 
589 }
constexpr amrex::Real Cp_d
Definition: ERF_Constants.H:12
constexpr amrex::Real R_d
Definition: ERF_Constants.H:10
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getThgivenTandP(const amrex::Real T, const amrex::Real P, const amrex::Real rdOcp)
Definition: ERF_EOS.H:18
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:43
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
#define RhoKE_comp
Definition: ERF_IndexDefines.H:38
AMREX_FORCE_INLINE IntVect offset(const int face_dir, const int normal)
Definition: ERF_ReadBndryPlanes.cpp:28
amrex::MultiFab * m_yvel
Definition: ERF_ShocInterface.H:577
amrex::MultiFab * m_mu
Definition: ERF_ShocInterface.H:590
amrex::MultiFab * m_cons
Definition: ERF_ShocInterface.H:573
amrex::MultiFab * m_xvel
Definition: ERF_ShocInterface.H:576
amrex::Vector< int > m_col_offsets
Definition: ERF_ShocInterface.H:548
Here is the call graph for this function:

◆ mf_to_kokkos_buffers()

void SHOCInterface::mf_to_kokkos_buffers ( )
336 {
337  //
338  // Expose for device capture
339  //
340 
341  // Interface data structures
342  //=======================================================
343  auto omega_d = omega;
344  auto surf_sens_flux_d = surf_sens_flux;
345  auto surf_mom_flux_d = surf_mom_flux;
346  auto surf_evap_d = surf_evap;
347  auto T_mid_d = T_mid;
348  auto qv_d = qv;
349  auto surf_drag_coeff_tms_d = surf_drag_coeff_tms;
350 
351  // Input data structures
352  //=======================================================
353  auto p_mid_d = p_mid;
354  auto p_int_d = p_int;
355  auto pseudo_dens_d = pseudo_dens;
356  auto phis_d = phis;
357 
358  // Input/Output data structures
359  //=======================================================
360  auto horiz_wind_d = horiz_wind;
361  auto cldfrac_liq_d = cldfrac_liq;
362  auto tke_d = tke;
363  auto qc_d = qc;
364 
365  // Enforce the correct grid heights and density
366  //=======================================================
367  auto dz_d = m_buffer.dz;
368 
369  // Subsidence pointer to device vector data
370  Real* w_sub = m_w_subsid;
371 
372  int nlay = m_num_layers;
373  Real dz = m_geom.CellSize(2);
374  bool moist = (m_cons->nComp() > RhoQ1_comp);
375  auto ProbLoArr = m_geom.ProbLoArray();
376  for (MFIter mfi(*m_cons); mfi.isValid(); ++mfi) {
377  const auto& vbx = mfi.validbox();
378  const int nx = vbx.length(0);
379  const int imin = vbx.smallEnd(0);
380  const int jmin = vbx.smallEnd(1);
381  const int offset = m_col_offsets[mfi.index()];
382  const Array4<const Real>& cons_arr = m_cons->const_array(mfi);
383 
384  const Array4<const Real>& u_arr = m_xvel->const_array(mfi);
385  const Array4<const Real>& v_arr = m_yvel->const_array(mfi);
386  const Array4<const Real>& w_arr = m_zvel->const_array(mfi);
387 
388  const Array4<const Real>& t13_arr = m_tau13->const_array(mfi);
389  const Array4<const Real>& t23_arr = m_tau23->const_array(mfi);
390  const Array4<const Real>& hfx3_arr = m_hfx3->const_array(mfi);
391  const Array4<const Real>& qfx3_arr = m_qfx3->const_array(mfi);
392 
393  const Array4<const Real>& mu_arr = m_mu->const_array(mfi);
394 
395  const Array4<const Real>& z_arr = (m_z_phys) ? m_z_phys->const_array(mfi) :
396  Array4<const Real>{};
397  ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
398  {
399  // map [i,j,k] 0-based to [icol, ilay] 0-based
400  const int icol = (j-jmin)*nx + (i-imin) + offset;
401  const int ilay = k;
402 
403  // EOS input (at CC)
404  Real r = cons_arr(i,j,k,Rho_comp);
405  Real rt = cons_arr(i,j,k,RhoTheta_comp);
406  Real qv = (moist) ? cons_arr(i,j,k,RhoQ1_comp)/r : 0.0;
407  Real qc = (moist) ? cons_arr(i,j,k,RhoQ2_comp)/r : 0.0;
408 
409  // EOS avg to z-face
410  Real r_lo = cons_arr(i,j,k-1,Rho_comp);
411  Real rt_lo = cons_arr(i,j,k-1,RhoTheta_comp);
412  Real qv_lo = (moist) ? cons_arr(i,j,k-1,RhoQ1_comp)/r_lo : 0.0;
413  Real r_avg = 0.5 * (r + r_lo);
414  Real rt_avg = 0.5 * (rt + rt_lo);
415  Real qv_avg = 0.5 * (qv + qv_lo);
416 
417  // Z height
418  Real z = (z_arr) ? 0.125 * ( (z_arr(i ,j ,k+1) + z_arr(i ,j ,k))
419  + (z_arr(i+1,j ,k+1) + z_arr(i+1,j ,k))
420  + (z_arr(i ,j+1,k+1) + z_arr(i ,j+1,k))
421  + (z_arr(i+1,j+1,k+1) + z_arr(i+1,j+1,k)) ) * CONST_GRAV :
422  ProbLoArr[2] + (k + 0.5) * dz;
423 
424  // Delta z
425  Real delz = (z_arr) ? 0.25 * ( (z_arr(i ,j ,k+1) - z_arr(i ,j ,k))
426  + (z_arr(i+1,j ,k+1) - z_arr(i+1,j ,k))
427  + (z_arr(i ,j+1,k+1) - z_arr(i ,j+1,k))
428  + (z_arr(i+1,j+1,k+1) - z_arr(i+1,j+1,k)) ) : dz;
429 
430  // W at cc (cannot be 0?; inspection of shoc code...)
431  Real w_cc = 0.5 * (w_arr(i,j,k) + w_arr(i,j,k+1));
432  w_cc += (w_sub) ? w_sub[k] : 0.0;
433  Real w_limited = std::copysign(std::max(std::fabs(w_cc),1.0e-6),w_cc);
434 
435 
436  // Interface data structures
437  //=======================================================
438  // eamxx_common_physics_functions_impl.hpp: calculate_vertical_velocity
439  omega_d(icol,ilay) = -w_limited * r * CONST_GRAV;
440  if (k==0) {
441  surf_mom_flux_d(icol,0) = 0.5 * (t13_arr(i,j,k) + t13_arr(i+1,j ,k));
442  surf_mom_flux_d(icol,1) = 0.5 * (t23_arr(i,j,k) + t23_arr(i ,j+1,k));
443  // Unit conversion to W/m^2 (eamxx_shoc_process_interface.cpp L62)
444  surf_sens_flux_d(icol) = hfx3_arr(i,j,k) * r * Cp_d;
445  surf_evap(icol) = (moist) ? qfx3_arr(i,j,k) : 0.0;
446  }
447  T_mid_d(icol,ilay) = getTgivenRandRTh(r, rt, qv);
448  qv_d(icol,ilay) = qv;
449  surf_drag_coeff_tms_d(icol) = 0.0;
450 
451  // Input data structures
452  //=======================================================
453  p_mid_d(icol,ilay) = getPgivenRTh(rt, qv);
454  p_int_d(icol,ilay) = getPgivenRTh(rt_avg, qv_avg);
455  // eamxx_common_physics_functions_impl.hpp: calculate_density
456  pseudo_dens_d(icol,ilay) = r * CONST_GRAV * delz;
457  // Enforce the grid spacing
458  dz_d(icol,ilay) = delz;
459  // Geopotential
460  phis_d(icol) = CONST_GRAV * z;
461 
462  if (ilay==(nlay-1)) {
463  Real r_hi = cons_arr(i,j,k+1,Rho_comp);
464  Real rt_hi = cons_arr(i,j,k+1,RhoTheta_comp);
465  Real qv_hi = (moist) ? std::max(cons_arr(i,j,k+1,RhoQ1_comp)/r_hi,0.0) : 0.0;
466  r_avg = 0.5 * (r + r_hi);
467  rt_avg = 0.5 * (rt + rt_hi);
468  qv_avg = 0.5 * (qv + qv_hi);
469  p_int_d(icol,ilay+1) = getPgivenRTh(rt_avg, qv_avg);
470  }
471 
472  // Input/Output data structures
473  //=======================================================
474  horiz_wind_d(icol,0,ilay) = 0.5 * (u_arr(i,j,k) + u_arr(i+1,j ,k));
475  horiz_wind_d(icol,1,ilay) = 0.5 * (v_arr(i,j,k) + v_arr(i ,j+1,k));
476  cldfrac_liq_d(icol,ilay) = (qc>0.0) ? 1. : 0.;
477  tke_d(icol,ilay) = std::max(cons_arr(i,j,k,RhoKE_comp)/r, 0.0);
478  qc_d(icol,ilay) = qc;
479  });
480  }
481 }
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=0.)
Definition: ERF_EOS.H:81
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, const amrex::Real qv=0.0)
Definition: ERF_EOS.H:46
amrex::MultiFab * m_z_phys
Definition: ERF_ShocInterface.H:593
amrex::MultiFab * m_qfx3
Definition: ERF_ShocInterface.H:587
amrex::MultiFab * m_zvel
Definition: ERF_ShocInterface.H:578
amrex::MultiFab * m_tau23
Definition: ERF_ShocInterface.H:585
amrex::MultiFab * m_tau13
Definition: ERF_ShocInterface.H:584
amrex::MultiFab * m_hfx3
Definition: ERF_ShocInterface.H:586
amrex::Real * m_w_subsid
Definition: ERF_ShocInterface.H:581
Here is the call graph for this function:

◆ name()

std::string SHOCInterface::name ( ) const
inline
117 { return "shoc"; }

◆ requested_buffer_size_in_bytes()

size_t SHOCInterface::requested_buffer_size_in_bytes ( ) const
protected
593 {
594  using TPF = ekat::TeamPolicyFactory<KT::ExeSpace>;
595 
596  const int nlev_packs = ekat::npack<Spack>(m_num_layers);
597  const int nlevi_packs = ekat::npack<Spack>(m_num_layers+1);
598  const int num_tracer_packs = ekat::npack<Spack>(m_num_tracers);
599 
600  // Number of Reals needed by local views in the interface
601  const size_t interface_request = Buffer::num_1d_scalar_ncol*m_num_cols*sizeof(Real) +
602  Buffer::num_1d_scalar_nlev*nlev_packs*sizeof(Spack) +
603  Buffer::num_2d_vector_mid*m_num_cols*nlev_packs*sizeof(Spack) +
604  Buffer::num_2d_vector_int*m_num_cols*nlevi_packs*sizeof(Spack) +
605  Buffer::num_2d_vector_tr*m_num_cols*num_tracer_packs*sizeof(Spack);
606 
607  // Number of Reals needed by the WorkspaceManager passed to shoc_main
608  const auto policy = TPF::get_default_team_policy(m_num_cols, nlev_packs);
609  const int n_wind_slots = ekat::npack<Spack>(m_num_vel_comp)*Spack::n;
610  const int n_trac_slots = ekat::npack<Spack>(m_num_tracers) *Spack::n;
611  const size_t wsm_request = WSM::get_total_bytes_needed(nlevi_packs, 14+(n_wind_slots+n_trac_slots), policy);
612 
613  return ( (interface_request + wsm_request)/sizeof(Real) );
614 }
static constexpr int num_2d_vector_tr
Definition: ERF_ShocInterface.H:454
static constexpr int num_1d_scalar_nlev
Definition: ERF_ShocInterface.H:446

◆ run_impl()

void SHOCInterface::run_impl ( const Real  dt)
852 {
853  using TPF = ekat::TeamPolicyFactory<KT::ExeSpace>;
854 
855  EKAT_REQUIRE_MSG (dt<=300,
856  "Error! SHOC is intended to run with a timestep no longer than 5 minutes.\n"
857  " Please, reduce timestep (perhaps increasing subcycling iterations).\n");
858 
859  const auto nlev_packs = ekat::npack<Spack>(m_num_layers);
860  const auto scan_policy = TPF::get_thread_range_parallel_scan_team_policy(m_num_cols, nlev_packs);
861  const auto default_policy = TPF::get_default_team_policy(m_num_cols, nlev_packs);
862 
863  // Preprocessing of SHOC inputs. Kernel contains a parallel_scan,
864  // so a special TeamPolicy is required.
865  Kokkos::parallel_for("shoc_preprocess",
866  scan_policy,
868  Kokkos::fence();
869 
870  auto wtracer_sfc = shoc_preprocess.wtracer_sfc;
871  Kokkos::deep_copy(wtracer_sfc, 0);
872 
873  if (apply_tms) {
875  }
876 
877  if (check_flux_state) {
879  }
880 
881  // For now set the host timestep to the shoc timestep. This forces
882  // number of SHOC timesteps (nadv) to be 1.
883  // TODO: input parameter?
884  hdtime = dt;
885  m_nadv = std::max(static_cast<int>(round(hdtime/dt)),1);
886 
887  // Reset internal WSM variables.
888  workspace_mgr.reset_internals();
889 
890  // Run shoc main
893 #ifdef SCREAM_SHOC_SMALL_KERNELS
894  , temporaries
895 #endif
896  );
897 
898  // Postprocessing of SHOC outputs
899  Kokkos::parallel_for("shoc_postprocess",
900  default_policy,
902  Kokkos::fence();
903 
904  // Extra SHOC output diagnostics
905  if (runtime_options.extra_diags) {
906  Kokkos::deep_copy(shoc_mix,history_output.shoc_mix);
907  Kokkos::deep_copy(brunt,history_output.brunt);
908  Kokkos::deep_copy(w3,history_output.w3);
909  Kokkos::deep_copy(isotropy,history_output.isotropy);
910  Kokkos::deep_copy(wthl_sec,history_output.wthl_sec);
911  Kokkos::deep_copy(wqw_sec,history_output.wqw_sec);
912  Kokkos::deep_copy(uw_sec,history_output.uw_sec);
913  Kokkos::deep_copy(vw_sec,history_output.vw_sec);
914  Kokkos::deep_copy(qw_sec,history_output.qw_sec);
915  Kokkos::deep_copy(thl_sec,history_output.thl_sec);
916  }
917 }
Int m_nadv
Definition: ERF_ShocInterface.H:555
void check_flux_state_consistency(const double dt)
Definition: ERF_ShocInterface.cpp:954
Int hdtime
Definition: ERF_ShocInterface.H:558
void apply_turbulent_mountain_stress()
Definition: ERF_ShocInterface.cpp:934

◆ set_computed_group_impl()

void SHOCInterface::set_computed_group_impl ( )
protected

◆ set_grids()

void SHOCInterface::set_grids ( int &  level,
const amrex::BoxArray &  ba,
amrex::Geometry &  geom,
amrex::MultiFab *  cons,
amrex::MultiFab *  xvel,
amrex::MultiFab *  yvel,
amrex::MultiFab *  zvel,
amrex::Real w_subsid,
amrex::MultiFab *  tau13,
amrex::MultiFab *  tau23,
amrex::MultiFab *  hfx3,
amrex::MultiFab *  qfx3,
amrex::MultiFab *  eddyDiffs,
amrex::MultiFab *  z_phys 
)
128 {
129  // Set data members that may change
130  m_lev = level;
131  m_geom = geom;
132  m_cons = cons;
133  m_xvel = xvel;
134  m_yvel = yvel;
135  m_zvel = zvel;
136  m_w_subsid = w_subsid;
137  m_tau13 = tau13;
138  m_tau23 = tau23;
139  m_hfx3 = hfx3;
140  m_qfx3 = qfx3;
141  m_mu = eddyDiffs;
142  m_z_phys = z_phys;
143 
144  // Ensure the boxes span klo -> khi
145  int klo = geom.Domain().smallEnd(2);
146  int khi = geom.Domain().bigEnd(2);
147 
148  // Reset vector of offsets for columnar data
149  m_num_layers = geom.Domain().length(2);
150 
151  int num_cols = 0;
152  m_col_offsets.clear();
153  m_col_offsets.resize(int(ba.size()));
154  for (amrex::MFIter mfi(*m_cons); mfi.isValid(); ++mfi) {
155  const auto& vbx = mfi.validbox();
156  AMREX_ALWAYS_ASSERT_WITH_MESSAGE((klo == vbx.smallEnd(2)) &&
157  (khi == vbx.bigEnd(2)),
158  "Vertical decomposition with shoc is not allowed.");
159  int nx = vbx.length(0);
160  int ny = vbx.length(1);
161  m_col_offsets[mfi.index()] = num_cols;
162  num_cols += nx * ny;
163  }
164 
165  // Resize variables that persist in memory
166  if (num_cols != m_num_cols) {
168  tk = view_2d();
169  sgs_buoy_flux = view_2d("sgs_buoy_flux", num_cols, m_num_layers);
170  tk = view_2d("eddy_diff_mom", num_cols, m_num_layers);
171  }
172  m_num_cols = num_cols;
173 
174  // Allocate the buffer arrays in ERF
175  alloc_buffers();
176 
177  // Allocate the m_buffer struct
178  init_buffers();
179 
180  // Fill the KOKKOS Views from AMReX MFs
182 }
@ tau23
Definition: ERF_DataStruct.H:30
@ tau13
Definition: ERF_DataStruct.H:30
void alloc_buffers()
Definition: ERF_ShocInterface.cpp:186
void init_buffers()
Definition: ERF_ShocInterface.cpp:617
int m_lev
Definition: ERF_ShocInterface.H:561
void mf_to_kokkos_buffers()
Definition: ERF_ShocInterface.cpp:335
@ xvel
Definition: ERF_IndexDefines.H:141
@ cons
Definition: ERF_IndexDefines.H:140
@ zvel
Definition: ERF_IndexDefines.H:143
@ yvel
Definition: ERF_IndexDefines.H:142

Member Data Documentation

◆ apply_tms

bool SHOCInterface::apply_tms = true
protected

◆ brunt

view_2d SHOCInterface::brunt
protected

◆ check_flux_state

bool SHOCInterface::check_flux_state = false
protected

◆ cldfrac_liq

view_2d SHOCInterface::cldfrac_liq
protected

◆ cldfrac_liq_prev

view_2d SHOCInterface::cldfrac_liq_prev
protected

◆ column_conservation_check

bool SHOCInterface::column_conservation_check = false
protected

◆ extra_shoc_diags

bool SHOCInterface::extra_shoc_diags = false
protected

◆ hdtime

Int SHOCInterface::hdtime
protected

◆ heat_flux

view_1d SHOCInterface::heat_flux
protected

◆ history_output

SHF::SHOCHistoryOutput SHOCInterface::history_output
protected

◆ horiz_wind

view_3d SHOCInterface::horiz_wind
protected

◆ ice_flux

view_1d SHOCInterface::ice_flux
protected

◆ input

SHF::SHOCInput SHOCInterface::input
protected

◆ input_output

SHF::SHOCInputOutput SHOCInterface::input_output
protected

◆ inv_qc_relvar

view_2d SHOCInterface::inv_qc_relvar
protected

◆ isotropy

view_2d SHOCInterface::isotropy
protected

◆ m_ba

amrex::BoxArray SHOCInterface::m_ba
protected

◆ m_buffer

Buffer SHOCInterface::m_buffer
protected

◆ m_col_offsets

amrex::Vector<int> SHOCInterface::m_col_offsets
protected

◆ m_cons

amrex::MultiFab* SHOCInterface::m_cons = nullptr
protected

◆ m_first_step

bool SHOCInterface::m_first_step = true
protected

◆ m_geom

amrex::Geometry SHOCInterface::m_geom
protected

◆ m_hfx3

amrex::MultiFab* SHOCInterface::m_hfx3 = nullptr
protected

◆ m_lev

int SHOCInterface::m_lev
protected

◆ m_mu

amrex::MultiFab* SHOCInterface::m_mu = nullptr
protected

◆ m_nadv

Int SHOCInterface::m_nadv
protected

◆ m_npbl

Int SHOCInterface::m_npbl
protected

◆ m_num_cols

Int SHOCInterface::m_num_cols = 0
protected

◆ m_num_layers

Int SHOCInterface::m_num_layers = 0
protected

◆ m_num_tracers

Int SHOCInterface::m_num_tracers = 3
protected

◆ m_num_vel_comp

Int SHOCInterface::m_num_vel_comp = 2
protected

◆ m_qfx3

amrex::MultiFab* SHOCInterface::m_qfx3 = nullptr
protected

◆ m_step

int SHOCInterface::m_step
protected

◆ m_tau13

amrex::MultiFab* SHOCInterface::m_tau13 = nullptr
protected

◆ m_tau23

amrex::MultiFab* SHOCInterface::m_tau23 = nullptr
protected

◆ m_w_subsid

amrex::Real* SHOCInterface::m_w_subsid = nullptr
protected

◆ m_xvel

amrex::MultiFab* SHOCInterface::m_xvel = nullptr
protected

◆ m_yvel

amrex::MultiFab* SHOCInterface::m_yvel = nullptr
protected

◆ m_z_phys

amrex::MultiFab* SHOCInterface::m_z_phys = nullptr
protected

◆ m_zvel

amrex::MultiFab* SHOCInterface::m_zvel = nullptr
protected

◆ obklen

view_1d SHOCInterface::obklen
protected

◆ omega

view_2d SHOCInterface::omega
protected

◆ output

SHF::SHOCOutput SHOCInterface::output
protected

◆ p_int

view_2d SHOCInterface::p_int
protected

◆ p_mid

view_2d SHOCInterface::p_mid
protected

◆ pblh

view_1d SHOCInterface::pblh
protected

◆ phis

view_1d SHOCInterface::phis
protected

◆ pseudo_dens

view_2d SHOCInterface::pseudo_dens
protected

◆ qc

view_2d SHOCInterface::qc
protected

◆ qtracers

view_3d_strided SHOCInterface::qtracers
protected

◆ qv

view_2d SHOCInterface::qv
protected

◆ qw_sec

view_2d SHOCInterface::qw_sec
protected

◆ runtime_options

SHF::SHOCRuntime SHOCInterface::runtime_options
protected

◆ sgs_buoy_flux

view_2d SHOCInterface::sgs_buoy_flux
protected

◆ shoc_cond

view_2d SHOCInterface::shoc_cond
protected

◆ shoc_evap

view_2d SHOCInterface::shoc_evap
protected

◆ shoc_mix

view_2d SHOCInterface::shoc_mix
protected

◆ shoc_postprocess

SHOCPostprocess SHOCInterface::shoc_postprocess
protected

◆ shoc_preprocess

SHOCPreprocess SHOCInterface::shoc_preprocess
protected

◆ surf_drag_coeff_tms

view_1d SHOCInterface::surf_drag_coeff_tms
protected

◆ surf_evap

view_1d SHOCInterface::surf_evap
protected

◆ surf_mom_flux

sview_2d SHOCInterface::surf_mom_flux
protected

◆ surf_sens_flux

view_1d SHOCInterface::surf_sens_flux
protected

◆ T_mid

view_2d SHOCInterface::T_mid
protected

◆ thl_sec

view_2d SHOCInterface::thl_sec
protected

◆ tk

view_2d SHOCInterface::tk
protected

◆ tke

view_2d SHOCInterface::tke
protected

◆ tkh

view_2d SHOCInterface::tkh
protected

◆ tot_buff_view

view_1d SHOCInterface::tot_buff_view
protected

◆ ustar

view_1d SHOCInterface::ustar
protected

◆ uw_sec

view_2d SHOCInterface::uw_sec
protected

◆ vapor_flux

view_1d SHOCInterface::vapor_flux
protected

◆ vw_sec

view_2d SHOCInterface::vw_sec
protected

◆ w3

view_2d SHOCInterface::w3
protected

◆ w_sec

view_2d SHOCInterface::w_sec
protected

◆ water_flux

view_1d SHOCInterface::water_flux
protected

◆ workspace_mgr

ekat::WorkspaceManager<Spack, KT::Device> SHOCInterface::workspace_mgr
protected

◆ wqw_sec

view_2d SHOCInterface::wqw_sec
protected

◆ wthl_sec

view_2d SHOCInterface::wthl_sec
protected

The documentation for this class was generated from the following files: