ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
ERF_TI_slow_rhs_pre.H
Go to the documentation of this file.
1 #include "ERF_SrcHeaders.H"
2 
3 /**
4  * Wrapper for calling the routine that creates the slow RHS
5  */
6  auto slow_rhs_fun_pre = [&](Vector<MultiFab>& S_rhs,
7  Vector<MultiFab>& S_old,
8  Vector<MultiFab>& S_data,
9  Vector<MultiFab>& S_scratch,
10  const Real old_step_time,
11  const Real old_stage_time,
12  const Real new_stage_time,
13  const int nrk)
14  {
15  //
16  // Define primitive variables for all later RK stages
17  // (We have already done this for the first RK step)
18  // Note that it is essential this happen before the call to make_mom_sources
19  // because some of the buoyancy routines use the primitive variables
20  //
21  if (nrk > 0) {
22  int ng_cons = S_data[IntVars::cons].nGrow();
23  cons_to_prim(S_data[IntVars::cons], ng_cons);
24  }
25 
26  BL_PROFILE("slow_rhs_fun_pre");
27  if (verbose) Print() << std::setprecision(timeprecision)
28  << "Making slow rhs at time " << old_stage_time
29  << " for fast variables advancing from " << old_step_time
30  << " to " << new_stage_time << std::endl;
31 
32  Real slow_dt = new_stage_time - old_step_time;
33 
34  // *************************************************************************
35  // Set up flux registers if using two_way coupling
36  // *************************************************************************
37  YAFluxRegister* fr_as_crse = nullptr;
38  YAFluxRegister* fr_as_fine = nullptr;
39  if (solverChoice.coupling_type == CouplingType::TwoWay && finest_level > 0) {
40  if (level < finest_level) {
41  fr_as_crse = getAdvFluxReg(level+1);
42  fr_as_crse->reset();
43  }
44  if (level > 0) {
45  fr_as_fine = getAdvFluxReg(level);
46  }
47  }
48 
49  // Canopy data for mom sources
50  MultiFab* forest_drag = (solverChoice.do_forest_drag) ?
51  m_forest_drag[level]->get_drag_field() : nullptr;
52 
53  // Immersed Forcing
54  MultiFab* terrain_blank = (solverChoice.terrain_type == TerrainType::ImmersedForcing) ?
55  terrain_blanking[level].get() : nullptr;
56 
57  // Update the total moisture variable *before* computing sources since this is used in
58  // the buoyancy calculation
59  if (solverChoice.moisture_type != MoistureType::None) {
60  make_qt(S_data[IntVars::cons], qt);
61  }
62 
63  MultiFab p0_to_use, base_to_use;
64  MultiFab *zpn_to_use, *zpc_to_use, *ax_to_use, *ay_to_use, *az_to_use, *dJ_to_use;
65 
66  // Moving terrain
67  std::unique_ptr<MultiFab> z_t_pert;
68  if ( solverChoice.terrain_type == TerrainType::MovingFittedMesh )
69  {
70  z_t_pert = std::make_unique<MultiFab>(S_data[IntVars::zmom].boxArray(), S_data[IntVars::zmom].DistributionMap(), 1, 1);
71  update_terrain_stage(level, old_step_time, old_stage_time, new_stage_time, slow_dt);
72 
73  p0_to_use = MultiFab(base_state_new[level], make_alias, BaseState::p0_comp, 1);
74  base_to_use = MultiFab(base_state_new[level], make_alias, 0, BaseState::num_comps);
75  zpn_to_use = z_phys_nd_src[level].get();
76  zpc_to_use = z_phys_cc_src[level].get();
77  ax_to_use = ax_src[level].get();
78  ay_to_use = ay_src[level].get();
79  az_to_use = az_src[level].get();
80  dJ_to_use = detJ_cc_src[level].get();
81 
82  } else {
83  p0_to_use = MultiFab(base_state[level], make_alias, BaseState::p0_comp, 1);
84  base_to_use = MultiFab(base_state[level], make_alias, 0, BaseState::num_comps);
85  zpn_to_use = z_phys_nd[level].get();
86  zpc_to_use = z_phys_cc[level].get();
87  ax_to_use = ax[level].get();
88  ay_to_use = ay[level].get();
89  az_to_use = az[level].get();
90  dJ_to_use = detJ_cc[level].get();
91  }
92 
93  // *****************************************************************************
94  // Construct the source terms for the cell-centered (conserved) variables
95  // *****************************************************************************
96  make_sources(level, nrk, slow_dt, old_stage_time,
97  S_data, S_prim, cc_src, base_state[level], zpc_to_use,
98  qheating_rates[level].get(),
99  terrain_blank, fine_geom, solverChoice,
100  mapfac[level],
101  dptr_rhotheta_src, dptr_rhoqt_src,
102  dptr_wbar_sub, d_rayleigh_ptrs_at_lev,
103  input_sounding_data, turbPert, true);
104 
105  // *****************************************************************************
106  // Define the pressure gradient
107  // *****************************************************************************
108  make_gradp_pert(level, solverChoice, fine_geom, S_data,
109  p0_to_use, *zpn_to_use, *zpc_to_use, get_eb(level), gradp[level]);
110 
111  // *****************************************************************************
112  // Define the buoyancy forcing term in the z-direction
113  // *****************************************************************************
114  make_buoyancy(S_data, S_prim, qt, buoyancy, fine_geom, solverChoice, base_to_use,
115  micro->Get_Qstate_Moist_Size(), get_eb(level), solverChoice.anelastic[level]);
116 
117  // *****************************************************************************
118  // Make remaining (not gradp or buoyancy) momentum sources
119  // *****************************************************************************
120  make_mom_sources(old_stage_time, S_data, *zpn_to_use, *zpc_to_use, stretched_dz_h[level],
121  xvel_new, yvel_new, zvel_new,
122  xmom_src, ymom_src, zmom_src,
123  base_to_use, forest_drag, terrain_blank,
124  cosPhi_m[level].get(), sinPhi_m[level].get(), fine_geom, solverChoice,
125  mapfac[level],
126  (solverChoice.have_geo_wind_profile) ? d_u_geos[level].data(): nullptr,
127  (solverChoice.have_geo_wind_profile) ? d_v_geos[level].data(): nullptr,
128  dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev,
129  input_sounding_data, true);
130 
131  // *****************************************************************************
132  // Add body sources if doing flow around a body
133  // *****************************************************************************
134  add_thin_body_sources(xmom_src, ymom_src, zmom_src,
135  xflux_imask[level], yflux_imask[level], zflux_imask[level],
136  thin_xforce[level], thin_yforce[level], thin_zforce[level]);
137 
138  // *****************************************************************************
139  // Define RHS for rho, rho_theta and momenta
140  // *****************************************************************************
141  erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_old, S_data,
142  S_prim, qt, S_scratch, xvel_new, yvel_new, zvel_new,
143  z_t_rk[level], cc_src, xmom_src, ymom_src, zmom_src, buoyancy,
144  (level > 0) ? &zmom_crse_rhs[level] : nullptr,
145  Tau[level], SmnSmn, eddyDiffs, Hfx1, Hfx2, Hfx3, Q1fx1, Q1fx2, Q1fx3, Q2fx3, Diss,
146  fine_geom, solverChoice, m_SurfaceLayer, domain_bcs_type_d, domain_bcs_type,
147  *zpn_to_use, *zpc_to_use, *ax_to_use, *ay_to_use, *az_to_use, *dJ_to_use,
148  stretched_dz_d[level], gradp[level],
149  mapfac[level], get_eb(level), fr_as_crse, fr_as_fine);
150 
151  // *****************************************************************************
152  // Update for moving terrain
153  // *****************************************************************************
154  if ( solverChoice.terrain_type == TerrainType::MovingFittedMesh )
155  {
156  MultiFab r_hse_new (base_state_new[level], make_alias, BaseState::r0_comp, 1);
157  MultiFab p_hse_new (base_state_new[level], make_alias, BaseState::p0_comp, 1);
158  MultiFab pi_hse_new (base_state_new[level], make_alias, BaseState::pi0_comp, 1);
159  MultiFab th_hse_new (base_state_new[level], make_alias, BaseState::th0_comp, 1);
160 
161  MultiFab* r0_new = &r_hse_new;
162  MultiFab* p0_new = &p_hse_new;
163  MultiFab* pi0_new = &pi_hse_new;
164  MultiFab* th0_new = &th_hse_new;
165 
166  // We define and evolve (rho theta)_0 in order to re-create p_0 in a way that is consistent
167  // with our update of (rho theta) but does NOT maintain dp_0 / dz = -rho_0 g. This is why
168  // we no longer discretize the vertical pressure gradient in perturbational form.
169  MultiFab rt0(p0->boxArray(),p0->DistributionMap(),1,1);
170  MultiFab rt0_new(p0->boxArray(),p0->DistributionMap(),1,1);
171  MultiFab r0_temp(p0->boxArray(),p0->DistributionMap(),1,1);
172 
173  // Remember this does NOT maintain dp_0 / dz = -rho_0 g, so we can no longer
174  // discretize the vertical pressure gradient in perturbational form.
175  AMREX_ALWAYS_ASSERT(solverChoice.advChoice.dycore_horiz_adv_type == AdvType::Centered_2nd);
176  AMREX_ALWAYS_ASSERT(solverChoice.advChoice.dycore_vert_adv_type == AdvType::Centered_2nd);
177 
178  Real dt_base = (new_stage_time - old_step_time);
179 
180  const GpuArray<Real, AMREX_SPACEDIM> dxInv = fine_geom.InvCellSizeArray();
181 
182  const Real l_rdOcp = solverChoice.rdOcp;
183 
184 #ifdef _OPENMP
185 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
186 #endif
187  for ( MFIter mfi(*p0,TilingIfNotGPU()); mfi.isValid(); ++mfi)
188  {
189  const Array4<Real > rt0_arr = rt0.array(mfi);
190  const Array4<Real > rt0_tmp_arr = rt0_new.array(mfi);
191 
192  const Array4<Real const> r0_arr = r0->const_array(mfi);
193  const Array4<Real > r0_new_arr = r0_new->array(mfi);
194  const Array4<Real > r0_tmp_arr = r0_temp.array(mfi);
195 
196  const Array4<Real const> p0_arr = p0->const_array(mfi);
197  const Array4<Real > p0_new_arr = p0_new->array(mfi);
198  const Array4<Real > pi0_new_arr = pi0_new->array(mfi);
199  const Array4<Real > th0_new_arr = th0_new->array(mfi);
200 
201  const Array4<Real >& z_t_arr = z_t_rk[level]->array(mfi);
202 
203  const Array4<Real const>& dJ_old_arr = detJ_cc[level]->const_array(mfi);
204  const Array4<Real const>& dJ_new_arr = detJ_cc_new[level]->const_array(mfi);
205  const Array4<Real const>& dJ_src_arr = detJ_cc_src[level]->const_array(mfi);
206 
207  Box gbx = mfi.growntilebox({1,1,1});
208  amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
209  {
210  rt0_arr(i,j,k) = getRhoThetagivenP(p0_arr(i,j,k));
211  rt0_tmp_arr(i,j,k) = getRhoThetagivenP(p0_new_arr(i,j,k));
212  r0_tmp_arr(i,j,k) = r0_new_arr(i,j,k);
213  });
214 
215  Box gbx2 = mfi.growntilebox({1,1,0});
216  amrex::ParallelFor(gbx2, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
217  {
218  Real zflux_r_lo = -z_t_arr(i,j,k ) * 0.5 * (r0_tmp_arr(i,j,k) + r0_tmp_arr(i,j,k-1));
219  Real zflux_r_hi = -z_t_arr(i,j,k+1) * 0.5 * (r0_tmp_arr(i,j,k) + r0_tmp_arr(i,j,k+1));
220 
221  Real zflux_rt_lo = zflux_r_lo * 0.5 * (rt0_tmp_arr(i,j,k)/r0_tmp_arr(i,j,k) + rt0_tmp_arr(i,j,k-1)/r0_tmp_arr(i,j,k-1));
222  Real zflux_rt_hi = zflux_r_hi * 0.5 * (rt0_tmp_arr(i,j,k)/r0_tmp_arr(i,j,k) + rt0_tmp_arr(i,j,k+1)/r0_tmp_arr(i,j,k+1));
223 
224  Real invdetJ = 1.0 / dJ_src_arr(i,j,k);
225 
226  Real src_r = - invdetJ * ( zflux_r_hi - zflux_r_lo ) * dxInv[2];
227  Real src_rt = - invdetJ * ( zflux_rt_hi - zflux_rt_lo ) * dxInv[2];
228 
229  Real rho0_new = dJ_old_arr(i,j,k) * r0_arr(i,j,k) + dt_base * dJ_src_arr(i,j,k) * src_r;
230  Real rt0_tmp_new = dJ_old_arr(i,j,k) * rt0_arr(i,j,k) + dt_base * dJ_src_arr(i,j,k) * src_rt;
231 
232  r0_new_arr(i,j,k) = rho0_new / dJ_new_arr(i,j,k);
233  rt0_tmp_new /= dJ_new_arr(i,j,k);
234 
235  p0_new_arr(i,j,k) = getPgivenRTh(rt0_tmp_new);
236  pi0_new_arr(i,j,k) = getExnergivenRTh(rt0_tmp_new, l_rdOcp);
237  th0_new_arr(i,j,k) = rt0_tmp_new / r0_new_arr(i,j,k);
238  });
239  } // MFIter
240  r0_new->FillBoundary(fine_geom.periodicity());
241  p0_new->FillBoundary(fine_geom.periodicity());
242  th0_new->FillBoundary(fine_geom.periodicity());
243  }
244 
245 #ifdef ERF_USE_NETCDF
246  // Populate RHS for relaxation zones if using real bcs
247  if (solverChoice.use_real_bcs && (level == 0)) {
248  if (real_width>0) {
249  realbdy_compute_interior_ghost_rhs(bdy_time_interval, start_bdy_time,
250  new_stage_time, slow_dt, stop_time,
251  real_width, real_set_width, fine_geom,
252  S_rhs, S_old, S_data,
253  bdy_data_xlo, bdy_data_xhi,
254  bdy_data_ylo, bdy_data_yhi);
255  }
256  }
257 #endif
258  }; // end slow_rhs_fun_pre
void add_thin_body_sources(MultiFab &xmom_src, MultiFab &ymom_src, MultiFab &zmom_src, std::unique_ptr< iMultiFab > &xflux_imask_lev, std::unique_ptr< iMultiFab > &yflux_imask_lev, std::unique_ptr< iMultiFab > &zflux_imask_lev, std::unique_ptr< MultiFab > &thin_xforce_lev, std::unique_ptr< MultiFab > &thin_yforce_lev, std::unique_ptr< MultiFab > &thin_zforce_lev)
Definition: ERF_AddThinBodySources.cpp:27
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 getExnergivenRTh(const amrex::Real rhotheta, const amrex::Real rdOcp, const amrex::Real qv=0.0)
Definition: ERF_EOS.H:156
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getRhoThetagivenP(const amrex::Real p, const amrex::Real qv=0.0)
Definition: ERF_EOS.H:172
@ Centered_2nd
void realbdy_compute_interior_ghost_rhs(const Real &bdy_time_interval, const Real &start_bdy_time, const Real &time, const Real &delta_t, const Real &stop_time, int width, int set_width, const Geometry &geom, Vector< MultiFab > &S_rhs, Vector< MultiFab > &S_old_data, Vector< MultiFab > &S_cur_data, Vector< Vector< FArrayBox >> &bdy_data_xlo, Vector< Vector< FArrayBox >> &bdy_data_xhi, Vector< Vector< FArrayBox >> &bdy_data_ylo, Vector< Vector< FArrayBox >> &bdy_data_yhi)
Definition: ERF_InteriorGhostCells.cpp:175
void make_buoyancy(const Vector< MultiFab > &S_data, const MultiFab &S_prim, const MultiFab &qt, MultiFab &buoyancy, const Geometry geom, const SolverChoice &solverChoice, const MultiFab &base_state, const int n_qstate, const eb_ &ebfact, const int anelastic)
Definition: ERF_MakeBuoyancy.cpp:32
void make_gradp_pert(int level, const SolverChoice &solverChoice, const Geometry &geom, Vector< MultiFab > &S_data, MultiFab &p0, MultiFab &z_phys_nd, MultiFab &z_phys_cc, const eb_ &ebfact, Vector< MultiFab > &gradp)
Definition: ERF_MakeGradP.cpp:28
void make_mom_sources(Real time, const Vector< MultiFab > &S_data, MultiFab &z_phys_nd, MultiFab &z_phys_cc, Vector< Real > &stretched_dz_h, const MultiFab &xvel, const MultiFab &yvel, const MultiFab &wvel, MultiFab &xmom_src, MultiFab &ymom_src, MultiFab &zmom_src, const MultiFab &base_state, MultiFab *forest_drag, MultiFab *terrain_blank, MultiFab *cosPhi_mf, MultiFab *sinPhi_mf, const Geometry geom, const SolverChoice &solverChoice, Vector< std::unique_ptr< MultiFab >> &, const Real *dptr_u_geos, const Real *dptr_v_geos, const Real *dptr_wbar_sub, const Vector< Real * > d_rayleigh_ptrs_at_lev, const Vector< Real * > d_sponge_ptrs_at_lev, InputSoundingData &input_sounding_data, bool is_slow_step)
Definition: ERF_MakeMomSources.cpp:33
void make_sources(int level, int, Real dt, Real time, const Vector< MultiFab > &S_data, const MultiFab &S_prim, MultiFab &source, const MultiFab &base_state, const MultiFab *z_phys_cc, const MultiFab *qheating_rates, MultiFab *terrain_blank, const Geometry geom, const SolverChoice &solverChoice, Vector< std::unique_ptr< MultiFab >> &mapfac, const Real *dptr_rhotheta_src, const Real *dptr_rhoqt_src, const Real *dptr_wbar_sub, const Vector< Real * > d_rayleigh_ptrs_at_lev, InputSoundingData &input_sounding_data, TurbulentPerturbation &turbPert, bool is_slow_step)
Definition: ERF_MakeSources.cpp:31
void erf_slow_rhs_pre(int level, int finest_level, int nrk, Real dt, Vector< MultiFab > &S_rhs, Vector< MultiFab > &S_old, Vector< MultiFab > &S_data, const MultiFab &S_prim, const MultiFab &qt, Vector< MultiFab > &S_scratch, const MultiFab &xvel, const MultiFab &yvel, const MultiFab &zvel, std::unique_ptr< MultiFab > &z_t_mf, const MultiFab &cc_src, const MultiFab &xmom_src, const MultiFab &ymom_src, const MultiFab &zmom_src, const MultiFab &buoyancy, const MultiFab *zmom_crse_rhs, Vector< std::unique_ptr< MultiFab >> &Tau_lev, MultiFab *SmnSmn, MultiFab *eddyDiffs, MultiFab *Hfx1, MultiFab *Hfx2, MultiFab *Hfx3, MultiFab *Q1fx1, MultiFab *Q1fx2, MultiFab *Q1fx3, MultiFab *Q2fx3, MultiFab *Diss, const Geometry geom, const SolverChoice &solverChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, const Gpu::DeviceVector< BCRec > &domain_bcs_type_d, const Vector< BCRec > &domain_bcs_type_h, const MultiFab &z_phys_nd, const MultiFab &z_phys_cc, const MultiFab &ax, const MultiFab &ay, const MultiFab &az, const MultiFab &detJ, Gpu::DeviceVector< Real > &stretched_dz_d, Vector< MultiFab > &gradp, Vector< std::unique_ptr< MultiFab >> &mapfac, const eb_ &ebfact, YAFluxRegister *fr_as_crse, YAFluxRegister *fr_as_fine)
Definition: ERF_SlowRhsPre.cpp:63
auto slow_rhs_fun_pre
Definition: ERF_TI_slow_rhs_pre.H:6
auto cons_to_prim
Definition: ERF_TI_utils.H:4
auto make_qt
Definition: ERF_TI_utils.H:50
auto update_terrain_stage
Definition: ERF_TI_utils.H:162
@ num_comps
Definition: ERF_IndexDefines.H:68
@ pi0_comp
Definition: ERF_IndexDefines.H:65
@ p0_comp
Definition: ERF_IndexDefines.H:64
@ th0_comp
Definition: ERF_IndexDefines.H:66
@ r0_comp
Definition: ERF_IndexDefines.H:63
@ cons
Definition: ERF_IndexDefines.H:158
@ zmom
Definition: ERF_IndexDefines.H:161
@ qt
Definition: ERF_Kessler.H:27
real(c_double), parameter p0
Definition: ERF_module_model_constants.F90:40