ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_TI_slow_rhs_post.H
Go to the documentation of this file.
1 #include "ERF_SrcHeaders.H"
2 
3  auto slow_rhs_fun_post = [&](Vector<MultiFab>& S_rhs,
4  Vector<MultiFab>& S_old,
5  Vector<MultiFab>& S_new,
6  Vector<MultiFab>& S_data,
7  Vector<MultiFab>& S_scratch,
8  const Real old_step_time,
9  const Real old_stage_time,
10  const Real new_stage_time,
11  const int nrk)
12  {
13  // Note that the "old" and "new" metric terms correspond to
14  // t^n and the RK stage (either t^*, t^** or t^{n+1} that this source
15  // will be used to advance to
16  Real slow_dt = new_stage_time - old_step_time;
17 
18  if (verbose) amrex::Print() << "Time integration of scalars at level " << level
19  << std::setprecision(timeprecision)
20  << " from " << old_step_time << " to " << new_stage_time
21  << " with dt = " << slow_dt
22  << " using RHS created at " << old_stage_time << std::endl;
23 
24 #if defined(ERF_USE_NETCDF)
25  bool moist_set_rhs = false;
26  if ( solverChoice.use_real_bcs && (level==0) && (real_set_width > 0) &&
27  (solverChoice.moisture_type != MoistureType::None) )
28  {
29  moist_set_rhs = true;
30  }
31 #endif
32 
33  // *************************************************************************
34  // Set up flux registers if using two_way coupling
35  // *************************************************************************
36  YAFluxRegister* fr_as_crse = nullptr;
37  YAFluxRegister* fr_as_fine = nullptr;
38  if (solverChoice.coupling_type == CouplingType::TwoWay && finest_level > 0)
39  {
40  if (level < finest_level) {
41  fr_as_crse = getAdvFluxReg(level+1);
42  }
43  if (level > 0) {
44  fr_as_fine = getAdvFluxReg(level);
45  }
46  }
47 
48  MultiFab* new_detJ =
49  (solverChoice.terrain_type == TerrainType::MovingFittedMesh) ? (detJ_cc_new[level].get()) : (detJ_cc[level].get());
50 
51  erf_slow_rhs_post(level, finest_level, nrk, slow_dt, micro->Get_Qstate_Moist_Size(),
52  S_rhs, S_old, S_new, S_data, S_prim, S_scratch,
53  xvel_new, yvel_new, zvel_new, cc_src, SmnSmn, eddyDiffs,
54  Hfx1, Hfx2, Hfx3, Q1fx1, Q1fx2, Q1fx3, Q2fx3, Diss,
55  fine_geom, solverChoice, m_SurfaceLayer, domain_bcs_type_d, domain_bcs_type,
56  z_phys_nd[level], z_phys_cc[level], ax[level], ay[level], az[level],
57  detJ_cc[level], new_detJ, stretched_dz_d[level], mapfac[level], EBFactory(level),
58 #if defined(ERF_USE_NETCDF)
59  moist_set_rhs, bdy_time_interval, new_stage_time, stop_time-start_time,
60  real_width, real_set_width,
61  bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi,
62 #endif
63  fr_as_crse, fr_as_fine);
64 
65  // Apply state redistribution for cons states
66 
67  if (solverChoice.terrain_type == TerrainType::EB)
68  {
69  Vector<int> is_valid_slow_var; is_valid_slow_var.resize(RhoQ1_comp+1,0);
70  if (solverChoice.turbChoice[level].use_tke) {is_valid_slow_var[ RhoKE_comp] = 1;}
71  is_valid_slow_var[RhoScalar_comp] = 1;
72  if (solverChoice.moisture_type != MoistureType::None) {
73  is_valid_slow_var[RhoQ1_comp] = 1;
74  }
75  const int num_comp_total = S_rhs[IntVars::cons].nComp();
76  const int num_grow = S_rhs[IntVars::cons].nGrow();
77  const int nvars = S_data[IntVars::cons].nComp();
78 
79  MultiFab dUdt_tmp(ba, dm, num_comp_total, num_grow, MFInfo(), EBFactory(level));
80  dUdt_tmp.setVal(0.0, 0, num_comp_total, num_grow);
81 
82  int start_comp;
83  int num_comp;
84 
85  for (int ivar(RhoKE_comp); ivar<= RhoQ1_comp; ++ivar)
86  {
87  if (is_valid_slow_var[ivar])
88  {
89  start_comp = ivar;
90  num_comp = 1;
91  if (ivar == RhoQ1_comp) {
92  num_comp = nvars - RhoQ1_comp;
93  } else if (ivar == RhoScalar_comp) {
94  num_comp = NSCALARS;
95  }
96  MultiFab::Copy(dUdt_tmp, S_rhs[IntVars::cons], start_comp, start_comp, num_comp, 0);
97  }
98  }
99  dUdt_tmp.FillBoundary(fine_geom.periodicity());
100  dUdt_tmp.setDomainBndry(1.234e10, 0, num_comp_total, fine_geom);
101 
102  const BCRec* bc_ptr_d = domain_bcs_type_d.data();
103 
104  // Update S_rhs by Redistribution.
105  // To-do: Currently, redistributing all the scalar variables.
106  // This needs to be redistributed only for num_comp variables starting from ivar, for efficiency.
107  redistribute_term ( num_comp_total, fine_geom, S_rhs[IntVars::cons], dUdt_tmp,
108  S_old[IntVars::cons], EBFactory(level), bc_ptr_d, slow_dt);
109 
110  // Update state using the updated S_rhs. (NOTE: redistribute_term returns RHS not state variables.)
111  for ( MFIter mfi(S_new[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
112  {
113  Box tbx = mfi.tilebox();
114  const Array4<Real>& snew = S_new[IntVars::cons].array(mfi);
115  const Array4<Real>& sold = S_old[IntVars::cons].array(mfi);
116  const Array4<Real>& srhs = S_rhs[IntVars::cons].array(mfi);
117  Array4<const Real> detJ_arr = EBFactory(level).getVolFrac().const_array(mfi);
118 
119  for (int ivar(RhoKE_comp); ivar<= RhoQ1_comp; ++ivar)
120  {
121  if (is_valid_slow_var[ivar])
122  {
123  start_comp = ivar;
124  num_comp = 1;
125  if (ivar == RhoQ1_comp) {
126  num_comp = nvars - RhoQ1_comp;
127  } else if (ivar == RhoScalar_comp) {
128  num_comp = NSCALARS;
129  }
130  ParallelFor(tbx, num_comp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn)
131  {
132  if (detJ_arr(i,j,k) > 0.0) {
133  const int n = start_comp + nn;
134  snew(i,j,k,n) = sold(i,j,k,n) + slow_dt * srhs(i,j,k,n);
135  }
136  });
137  }
138  }
139  }
140  } // EB
141 
142  // Apply boundary conditions on all the state variables that have been updated
143  // in both the fast and slow integrators
144  apply_bcs(S_new, new_stage_time, S_new[IntVars::cons].nGrow(), S_new[IntVars::xmom].nGrow(),
145  fast_only=false, vel_and_mom_synced=false);
146 
147  if (solverChoice.moisture_tight_coupling) {
148  // TODO: need iteration var for lagrangian microphysics
149  // call signature in ERF::Advance() is
150  //advance_microphysics(lev, S_new, dt_lev, iteration, time);
151  advance_microphysics(level, S_new[0], slow_dt, 123456789, old_step_time);
152  }
153  }; // end slow_rhs_fun_post
@ nvars
Definition: ERF_DataStruct.H:91
void redistribute_term(int ncomp, const Geometry &geom, MultiFab &result, MultiFab &result_tmp, MultiFab const &state, EBFArrayBoxFactory const &ebfact, BCRec const *bc, Real const local_dt)
Definition: ERF_EBRedistribute.cpp:13
#define RhoScalar_comp
Definition: ERF_IndexDefines.H:40
#define NSCALARS
Definition: ERF_IndexDefines.H:16
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
#define RhoKE_comp
Definition: ERF_IndexDefines.H:38
amrex::Real Real
Definition: ERF_ShocInterface.H:19
void erf_slow_rhs_post(int level, int finest_level, int nrk, Real dt, int n_qstate, Vector< MultiFab > &S_rhs, Vector< MultiFab > &S_old, Vector< MultiFab > &S_new, Vector< MultiFab > &S_data, const MultiFab &S_prim, Vector< MultiFab > &S_scratch, const MultiFab &xvel, const MultiFab &yvel, const MultiFab &, const MultiFab &source, const MultiFab *SmnSmn, const 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, std::unique_ptr< MultiFab > &z_phys_nd, std::unique_ptr< MultiFab > &z_phys_cc, std::unique_ptr< MultiFab > &ax, std::unique_ptr< MultiFab > &ay, std::unique_ptr< MultiFab > &az, std::unique_ptr< MultiFab > &detJ, MultiFab *detJ_new, Gpu::DeviceVector< Real > &stretched_dz_d, Vector< std::unique_ptr< MultiFab >> &mapfac, amrex::EBFArrayBoxFactory const &ebfact, YAFluxRegister *fr_as_crse, YAFluxRegister *fr_as_fine)
Definition: ERF_SlowRhsPost.cpp:45
auto slow_rhs_fun_post
Definition: ERF_TI_slow_rhs_post.H:3
auto apply_bcs
Definition: ERF_TI_utils.H:71
@ cons
Definition: ERF_IndexDefines.H:158
@ xmom
Definition: ERF_IndexDefines.H:159