ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_MakeFastCoeffs.cpp File Reference
#include <AMReX.H>
#include <ERF_TI_fast_headers.H>
Include dependency graph for ERF_MakeFastCoeffs.cpp:

Functions

void make_fast_coeffs (int, MultiFab &fast_coeffs, Vector< MultiFab > &S_stage_data, const MultiFab &S_stage_prim, const MultiFab &pi_stage, const amrex::Geometry geom, bool l_use_moisture, MeshType mesh_type, Real gravity, Real c_p, std::unique_ptr< MultiFab > &detJ_cc, const MultiFab *r0, const MultiFab *pi0, Real dtau, Real beta_s, amrex::GpuArray< ERF_BC, AMREX_SPACEDIM *2 > &phys_bc_type)
 

Function Documentation

◆ make_fast_coeffs()

void make_fast_coeffs ( int  ,
MultiFab &  fast_coeffs,
Vector< MultiFab > &  S_stage_data,
const MultiFab &  S_stage_prim,
const MultiFab &  pi_stage,
const amrex::Geometry  geom,
bool  l_use_moisture,
MeshType  mesh_type,
Real  gravity,
Real  c_p,
std::unique_ptr< MultiFab > &  detJ_cc,
const MultiFab *  r0,
const MultiFab *  pi0,
Real  dtau,
Real  beta_s,
amrex::GpuArray< ERF_BC, AMREX_SPACEDIM *2 > &  phys_bc_type 
)

Function for computing the coefficients for the tridiagonal solver used in the fast integrator (the acoustic substepping).

Parameters
[in]levellevel of refinement
[out]fast_coeffsthe coefficients for the tridiagonal solver computed here
[in]S_stage_datasolution at the last stage
[in]S_stage_primprimitive variables (i.e. conserved variables divided by density) at the last stage
[in]pi_stageExner function at the last stage
[in]geomContainer for geometric information
[in]mesh_typeDo we have constant dz?
[in]gravityMagnitude of gravity
[in]c_pCoefficient at constant pressure
[in]r0Reference (hydrostatically stratified) density
[in]pi0Reference (hydrostatically stratified) Exner function
[in]dtauFast time step
[in]beta_sCoefficient which determines how implicit vs explicit the solve is
39 {
40  BL_PROFILE_VAR("make_fast_coeffs()",make_fast_coeffs);
41 
42  Real beta_2 = myhalf * (one + beta_s); // multiplies implicit terms
43 
44  Real c_v = c_p - R_d;
45  Real RvOverRd = R_v / R_d;
46 
47  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
48  Real dzi = dxInv[2];
49 
50  const Box &domain = geom.Domain();
51 
52  MultiFab coeff_A_mf(fast_coeffs, amrex::make_alias, 0, 1);
53  MultiFab coeff_B_mf(fast_coeffs, amrex::make_alias, 1, 1);
54  MultiFab coeff_C_mf(fast_coeffs, amrex::make_alias, 2, 1);
55  MultiFab coeff_P_mf(fast_coeffs, amrex::make_alias, 3, 1);
56  MultiFab coeff_Q_mf(fast_coeffs, amrex::make_alias, 4, 1);
57 
58 
59  // *************************************************************************
60  // Set gravity as a vector
61  const Array<Real,AMREX_SPACEDIM> grav{zero, zero, -gravity};
62  const GpuArray<Real,AMREX_SPACEDIM> grav_gpu{grav[0], grav[1], grav[2]};
63 
64  // *************************************************************************
65  // Define updates in the current RK stage
66  // *************************************************************************
67 #ifdef _OPENMP
68 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
69 #endif
70  {
71 
72  for ( MFIter mfi(S_stage_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi)
73  {
74  Box bx = mfi.tilebox();
75  Box tbz = surroundingNodes(bx,2);
76 
77  const Array4<const Real> & stage_cons = S_stage_data[IntVars::cons].const_array(mfi);
78  const Array4<const Real> & prim = S_stage_prim.const_array(mfi);
79 
80  const Array4<const Real>& detJ = (mesh_type != MeshType::ConstantDz) ? detJ_cc->const_array(mfi) : Array4<const Real>{};
81 
82  const Array4<const Real>& r0_ca = r0->const_array(mfi);
83  const Array4<const Real>& pi0_ca = pi0->const_array(mfi);
84  const Array4<const Real>& pi_stage_ca = pi_stage.const_array(mfi);
85 
86  FArrayBox gam_fab; gam_fab.resize(surroundingNodes(bx,2),1,The_Async_Arena());
87 
88  auto const& coeffA_a = coeff_A_mf.array(mfi);
89  auto const& coeffB_a = coeff_B_mf.array(mfi);
90  auto const& coeffC_a = coeff_C_mf.array(mfi);
91  auto const& coeffP_a = coeff_P_mf.array(mfi);
92  auto const& coeffQ_a = coeff_Q_mf.array(mfi);
93  auto const& gam_a = gam_fab.array();
94 
95  // *********************************************************************
96  // *********************************************************************
97  // *********************************************************************
98 
99  Box bx_shrunk_in_k = bx;
100  int klo = tbz.smallEnd(2);
101  int khi = tbz.bigEnd(2);
102  bx_shrunk_in_k.setSmall(2,klo+1);
103  bx_shrunk_in_k.setBig(2,khi-1);
104 
105  // Note that the notes use "g" to mean the magnitude of gravity, so it is positive
106  // We set grav_gpu[2] to be the vector component which is negative
107  // We define halfg to match the notes (which is why we take the absolute value)
108  Real halfg = std::abs(myhalf * grav_gpu[2]);
109 
110  //Note we don't act on the bottom or top boundaries of the domain
111  if (mesh_type != MeshType::ConstantDz)
112  {
113  ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (int i, int j, int k)
114  {
115  Real rhobar_lo, rhobar_hi, pibar_lo, pibar_hi;
116  rhobar_lo = r0_ca(i,j,k-1);
117  rhobar_hi = r0_ca(i,j,k );
118  pibar_lo = pi0_ca(i,j,k-1);
119  pibar_hi = pi0_ca(i,j,k );
120 
121  Real pi_c = myhalf * (pi_stage_ca(i,j,k-1) + pi_stage_ca(i,j,k));
122 
123  Real detJ_on_kface = myhalf * (detJ(i,j,k) + detJ(i,j,k-1));
124  Real inv_detJ_on_kface = one / detJ_on_kface;
125 
126  Real qv_p = (l_use_moisture) ? prim(i,j,k ,PrimQ1_comp) : zero;
127  Real qv_q = (l_use_moisture) ? prim(i,j,k-1,PrimQ1_comp) : zero;
128 
129  Real coeff_P = -Gamma * R_d * dzi * inv_detJ_on_kface * pi_c * (one + RvOverRd*qv_p)
130  + halfg * R_d * rhobar_hi * pi_stage_ca(i,j,k) /
131  ( c_v * pibar_hi * stage_cons(i,j,k,RhoTheta_comp) );
132 
133  Real coeff_Q = Gamma * R_d * dzi * inv_detJ_on_kface * pi_c * (one + RvOverRd*qv_q)
134  + halfg * R_d * rhobar_lo * pi_stage_ca(i,j,k-1) /
135  ( c_v * pibar_lo * stage_cons(i,j,k-1,RhoTheta_comp) );
136 
137  coeffP_a(i,j,k) = coeff_P;
138  coeffQ_a(i,j,k) = coeff_Q;
139 
140  if (l_use_moisture) {
141  Real q = myhalf * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp)
142  + prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) );
143  coeff_P /= (one + q);
144  coeff_Q /= (one + q);
145  }
146 
147  Real theta_t_lo = myhalf * ( prim(i,j,k-2,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp) );
148  Real theta_t_mid = myhalf * ( prim(i,j,k-1,PrimTheta_comp) + prim(i,j,k ,PrimTheta_comp) );
149  Real theta_t_hi = myhalf * ( prim(i,j,k ,PrimTheta_comp) + prim(i,j,k+1,PrimTheta_comp) );
150 
151  // LHS for tri-diagonal system
152  Real D = dtau * dtau * beta_2 * beta_2 * dzi;
153  coeffA_a(i,j,k) = D * (one/detJ(i,j,k-1)) * ( halfg - coeff_Q * theta_t_lo );
154  coeffC_a(i,j,k) = D * (one/detJ(i,j,k )) * (-halfg + coeff_P * theta_t_hi );
155 
156  coeffB_a(i,j,k) = one + D * (coeff_Q/detJ(i,j,k-1) - coeff_P/detJ(i,j,k)) * theta_t_mid;
157  });
158 
159  } else {
160 
161  ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (int i, int j, int k)
162  {
163  Real rhobar_lo, rhobar_hi, pibar_lo, pibar_hi;
164  rhobar_lo = r0_ca(i,j,k-1);
165  rhobar_hi = r0_ca(i,j,k );
166  pibar_lo = pi0_ca(i,j,k-1);
167  pibar_hi = pi0_ca(i,j,k );
168 
169  Real pi_c = myhalf * (pi_stage_ca(i,j,k-1) + pi_stage_ca(i,j,k));
170 
171  Real qv_p = (l_use_moisture) ? prim(i,j,k ,PrimQ1_comp) : zero;
172  Real qv_q = (l_use_moisture) ? prim(i,j,k-1,PrimQ1_comp) : zero;
173 
174  Real coeff_P = -Gamma * R_d * dzi * pi_c * (one + RvOverRd*qv_p)
175  + halfg * R_d * rhobar_hi * pi_stage_ca(i,j,k) /
176  ( c_v * pibar_hi * stage_cons(i,j,k,RhoTheta_comp) );
177 
178  Real coeff_Q = Gamma * R_d * dzi * pi_c * (one + RvOverRd*qv_q)
179  + halfg * R_d * rhobar_lo * pi_stage_ca(i,j,k-1) /
180  ( c_v * pibar_lo * stage_cons(i,j,k-1,RhoTheta_comp) );
181 
182  coeffP_a(i,j,k) = coeff_P;
183  coeffQ_a(i,j,k) = coeff_Q;
184 
185  if (l_use_moisture) {
186  Real q = myhalf * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp)
187  + prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) );
188  coeff_P /= (one + q);
189  coeff_Q /= (one + q);
190  }
191 
192  Real theta_t_lo = myhalf * ( prim(i,j,k-2,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp) );
193  Real theta_t_mid = myhalf * ( prim(i,j,k-1,PrimTheta_comp) + prim(i,j,k ,PrimTheta_comp) );
194  Real theta_t_hi = myhalf * ( prim(i,j,k ,PrimTheta_comp) + prim(i,j,k+1,PrimTheta_comp) );
195 
196  // LHS for tri-diagonal system
197  Real D = dtau * dtau * beta_2 * beta_2 * dzi;
198  coeffA_a(i,j,k) = D * ( halfg - coeff_Q * theta_t_lo );
199  coeffC_a(i,j,k) = D * (-halfg + coeff_P * theta_t_hi );
200 
201  coeffB_a(i,j,k) = one + D * (coeff_Q - coeff_P) * theta_t_mid;
202  });
203  }
204 
205  amrex::Box b2d = tbz; // Copy constructor
206  b2d.setRange(2,0);
207 
208  auto const lo = amrex::lbound(bx);
209  auto const hi = amrex::ubound(bx);
210 
211  auto const domhi = amrex::ubound(domain);
212 
213  {
214  BL_PROFILE("make_coeffs_b2d_loop");
215 #ifdef AMREX_USE_GPU
216  ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int) {
217 
218  // If at the bottom of the grid, we will set w to a specified Dirichlet value
219  coeffA_a(i,j,lo.z) = zero;
220  coeffB_a(i,j,lo.z) = one;
221  coeffC_a(i,j,lo.z) = zero;
222 
223  // If at the top of the grid, we will set w to a specified Dirichlet value
224  coeffA_a(i,j,hi.z+1) = zero;
225  coeffB_a(i,j,hi.z+1) = one;
226  coeffC_a(i,j,hi.z+1) = zero;
227 
228  // UNLESS if at the top of the domain and the boundary is outflow,
229  // we will use a homogeneous Neumann condition
230  if ( (hi.z == domhi.z) &&
231  (phys_bc_type[5] == ERF_BC::outflow or phys_bc_type[5] == ERF_BC::ho_outflow) )
232  {
233  coeffA_a(i,j,hi.z+1) = -one;
234  }
235 
236  // w = specified Dirichlet value at k = lo.z
237  Real bet = coeffB_a(i,j,lo.z);
238 
239  for (int k = lo.z+1; k <= hi.z+1; k++) {
240  gam_a(i,j,k) = coeffC_a(i,j,k-1) / bet;
241  bet = coeffB_a(i,j,k) - coeffA_a(i,j,k)*gam_a(i,j,k);
242  coeffB_a(i,j,k) = bet;
243  }
244  });
245 #else
246  // If at the bottom of the grid, we will set w to a specified Dirichlet value
247  for (int j = lo.y; j <= hi.y; ++j) {
248  AMREX_PRAGMA_SIMD
249  for (int i = lo.x; i <= hi.x; ++i) {
250  coeffA_a(i,j,lo.z) = zero;
251  coeffB_a(i,j,lo.z) = one;
252  coeffC_a(i,j,lo.z) = zero;
253  }
254  }
255  for (int j = lo.y; j <= hi.y; ++j) {
256  AMREX_PRAGMA_SIMD
257  for (int i = lo.x; i <= hi.x; ++i) {
258 
259  // If at the top of the grid, we will set w to a specified Dirichlet value
260  coeffA_a(i,j,hi.z+1) = zero;
261  coeffB_a(i,j,hi.z+1) = one;
262  coeffC_a(i,j,hi.z+1) = zero;
263 
264  // UNLESS if at the top of the domain and the boundary is outflow,
265  // we will use a homogeneous Neumann condition
266  if ( (hi.z == domhi.z) &&
267  (phys_bc_type[5] == ERF_BC::outflow or phys_bc_type[5] == ERF_BC::ho_outflow) )
268  {
269  coeffA_a(i,j,hi.z+1) = -one;
270  }
271  }
272  }
273  for (int k = lo.z+1; k <= hi.z+1; ++k) {
274  for (int j = lo.y; j <= hi.y; ++j) {
275  AMREX_PRAGMA_SIMD
276  for (int i = lo.x; i <= hi.x; ++i) {
277  gam_a(i,j,k) = coeffC_a(i,j,k-1) / coeffB_a(i,j,k-1);
278  Real bet = coeffB_a(i,j,k) - coeffA_a(i,j,k)*gam_a(i,j,k);
279  coeffB_a(i,j,k) = bet;
280  }
281  }
282  }
283 #endif
284  } // end profile
285 
286  // In the end we save the inverse of the diagonal (B) coefficient
287  {
288  BL_PROFILE("make_coeffs_invert");
289  ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (int i, int j, int k)
290  {
291  coeffB_a(i,j,k) = one / coeffB_a(i,j,k);
292  });
293  } // end profile
294  } // mfi
295  } // omp
296 }
R_v
constexpr amrex::Real R_v
Definition: ERF_Constants.H:21
one
constexpr amrex::Real one
Definition: ERF_Constants.H:7
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:11
R_d
constexpr amrex::Real R_d
Definition: ERF_Constants.H:20
Gamma
constexpr amrex::Real Gamma
Definition: ERF_Constants.H:29
PrimQ1_comp
#define PrimQ1_comp
Definition: ERF_IndexDefines.H:53
PrimQ2_comp
#define PrimQ2_comp
Definition: ERF_IndexDefines.H:54
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
PrimTheta_comp
#define PrimTheta_comp
Definition: ERF_IndexDefines.H:50
ERF_BC::ho_outflow
@ ho_outflow
ERF_BC::outflow
@ outflow
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
khi
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
ParallelFor
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const auto prob_lo=geomdata.ProbLo();const auto dx=geomdata.CellSize();const Real x=(prob_lo[0]+(i+myhalf) *dx[0])/mf_m(i, j, 0);const Real z=z_cc(i, j, k);Real L=std::sqrt(std::pow((x - x_c)/x_r, 2)+std::pow((z - z_c)/z_r, 2));if(L<=one) { Real dT=T_pert *(std::cos(PI *L)+one)/two;Real Tbar_hse=p_hse(i, j, k)/(R_d *r_hse(i, j, k));Real theta_perturbed=(Tbar_hse+dT) *std::pow(p_0/p_hse(i, j, k), rdOcp);Real theta_0=(Tbar_hse) *std::pow(p_0/p_hse(i, j, k), rdOcp);if(const_rho) { state_pert(i, j, k, RhoTheta_comp)=r_hse(i, j, k) *(theta_perturbed - theta_0);} else { state_pert(i, j, k, Rho_comp)=getRhoThetagivenP(p_hse(i, j, k))/theta_perturbed - r_hse(i, j, k);} } })
make_fast_coeffs
void make_fast_coeffs(int, MultiFab &fast_coeffs, Vector< MultiFab > &S_stage_data, const MultiFab &S_stage_prim, const MultiFab &pi_stage, const amrex::Geometry geom, bool l_use_moisture, MeshType mesh_type, Real gravity, Real c_p, std::unique_ptr< MultiFab > &detJ_cc, const MultiFab *r0, const MultiFab *pi0, Real dtau, Real beta_s, amrex::GpuArray< ERF_BC, AMREX_SPACEDIM *2 > &phys_bc_type)
Definition: ERF_MakeFastCoeffs.cpp:26
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:176
mp_wsm6::r0
real(kind=kind_phys), parameter, private r0
Definition: ERF_module_mp_wsm6.F90:21
Here is the call graph for this function: