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 = 0.5 * (1.0 + 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{0.0, 0.0, -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(0.5 * 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 = 0.5 * (pi_stage_ca(i,j,k-1) + pi_stage_ca(i,j,k));
122 
123  Real detJ_on_kface = 0.5 * (detJ(i,j,k) + detJ(i,j,k-1));
124  Real inv_detJ_on_kface = 1. / detJ_on_kface;
125 
126  Real qv_p = (l_use_moisture) ? prim(i,j,k ,PrimQ1_comp) : 0.0;
127  Real qv_q = (l_use_moisture) ? prim(i,j,k-1,PrimQ1_comp) : 0.0;
128 
129  Real coeff_P = -Gamma * R_d * dzi * inv_detJ_on_kface * (1.0 + RvOverRd*qv_p)
130  + halfg * R_d * rhobar_hi /
131  ( c_v * pibar_hi * stage_cons(i,j,k,RhoTheta_comp) );
132  coeff_P *= pi_c;
133 
134  Real coeff_Q = Gamma * R_d * dzi * inv_detJ_on_kface * (1.0 + RvOverRd*qv_q)
135  + halfg * R_d * rhobar_lo /
136  ( c_v * pibar_lo * stage_cons(i,j,k-1,RhoTheta_comp) );
137  coeff_Q *= pi_c;
138 
139  coeffP_a(i,j,k) = coeff_P;
140  coeffQ_a(i,j,k) = coeff_Q;
141 
142  if (l_use_moisture) {
143  Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp)
144  +prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) );
145  coeff_P /= (1.0 + q);
146  coeff_Q /= (1.0 + q);
147  }
148 
149  Real theta_t_lo = 0.5 * ( prim(i,j,k-2,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp) );
150  Real theta_t_mid = 0.5 * ( prim(i,j,k-1,PrimTheta_comp) + prim(i,j,k ,PrimTheta_comp) );
151  Real theta_t_hi = 0.5 * ( prim(i,j,k ,PrimTheta_comp) + prim(i,j,k+1,PrimTheta_comp) );
152 
153  // LHS for tri-diagonal system
154  Real D = dtau * dtau * beta_2 * beta_2 * dzi;
155  coeffA_a(i,j,k) = D * ( halfg - coeff_Q * theta_t_lo );
156  coeffC_a(i,j,k) = D * (-halfg + coeff_P * theta_t_hi );
157 
158  coeffB_a(i,j,k) = detJ_on_kface + D * (coeff_Q - coeff_P) * theta_t_mid;
159  });
160 
161  } else {
162 
163  ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (int i, int j, int k)
164  {
165  Real rhobar_lo, rhobar_hi, pibar_lo, pibar_hi;
166  rhobar_lo = r0_ca(i,j,k-1);
167  rhobar_hi = r0_ca(i,j,k );
168  pibar_lo = pi0_ca(i,j,k-1);
169  pibar_hi = pi0_ca(i,j,k );
170 
171  Real pi_c = 0.5 * (pi_stage_ca(i,j,k-1) + pi_stage_ca(i,j,k));
172 
173  Real qv_p = (l_use_moisture) ? prim(i,j,k ,PrimQ1_comp) : 0.0;
174  Real qv_q = (l_use_moisture) ? prim(i,j,k-1,PrimQ1_comp) : 0.0;
175 
176  Real coeff_P = -Gamma * R_d * dzi * (1.0 + RvOverRd*qv_p)
177  + halfg * R_d * rhobar_hi /
178  ( c_v * pibar_hi * stage_cons(i,j,k,RhoTheta_comp) );
179  coeff_P *= pi_c;
180 
181  Real coeff_Q = Gamma * R_d * dzi * (1.0 + RvOverRd*qv_q)
182  + halfg * R_d * rhobar_lo /
183  ( c_v * pibar_lo * stage_cons(i,j,k-1,RhoTheta_comp) );
184  coeff_Q *= pi_c;
185 
186  coeffP_a(i,j,k) = coeff_P;
187  coeffQ_a(i,j,k) = coeff_Q;
188 
189  if (l_use_moisture) {
190  Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp)
191  +prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) );
192  coeff_P /= (1.0 + q);
193  coeff_Q /= (1.0 + q);
194  }
195 
196  Real theta_t_lo = 0.5 * ( prim(i,j,k-2,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp) );
197  Real theta_t_mid = 0.5 * ( prim(i,j,k-1,PrimTheta_comp) + prim(i,j,k ,PrimTheta_comp) );
198  Real theta_t_hi = 0.5 * ( prim(i,j,k ,PrimTheta_comp) + prim(i,j,k+1,PrimTheta_comp) );
199 
200  // LHS for tri-diagonal system
201  Real D = dtau * dtau * beta_2 * beta_2 * dzi;
202  coeffA_a(i,j,k) = D * ( halfg - coeff_Q * theta_t_lo );
203  coeffC_a(i,j,k) = D * (-halfg + coeff_P * theta_t_hi );
204 
205  coeffB_a(i,j,k) = 1.0 + D * (coeff_Q - coeff_P) * theta_t_mid;
206  });
207  }
208 
209  amrex::Box b2d = tbz; // Copy constructor
210  b2d.setRange(2,0);
211 
212  auto const lo = amrex::lbound(bx);
213  auto const hi = amrex::ubound(bx);
214 
215  auto const domhi = amrex::ubound(domain);
216 
217  {
218  BL_PROFILE("make_coeffs_b2d_loop");
219 #ifdef AMREX_USE_GPU
220  ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int) {
221 
222  // If at the bottom of the grid, we will set w to a specified Dirichlet value
223  coeffA_a(i,j,lo.z) = 0.0;
224  coeffB_a(i,j,lo.z) = 1.0;
225  coeffC_a(i,j,lo.z) = 0.0;
226 
227  // If at the top of the grid, we will set w to a specified Dirichlet value
228  coeffA_a(i,j,hi.z+1) = 0.0;
229  coeffB_a(i,j,hi.z+1) = 1.0;
230  coeffC_a(i,j,hi.z+1) = 0.0;
231 
232  // UNLESS if at the top of the domain and the boundary is outflow,
233  // we will use a homogeneous Neumann condition
234  if ( (hi.z == domhi.z) &&
235  (phys_bc_type[5] == ERF_BC::outflow or phys_bc_type[5] == ERF_BC::ho_outflow) )
236  {
237  coeffA_a(i,j,hi.z+1) = -1.0;
238  }
239 
240  // w = specified Dirichlet value at k = lo.z
241  Real bet = coeffB_a(i,j,lo.z);
242 
243  for (int k = lo.z+1; k <= hi.z+1; k++) {
244  gam_a(i,j,k) = coeffC_a(i,j,k-1) / bet;
245  bet = coeffB_a(i,j,k) - coeffA_a(i,j,k)*gam_a(i,j,k);
246  coeffB_a(i,j,k) = bet;
247  }
248  });
249 #else
250  // If at the bottom of the grid, we will set w to a specified Dirichlet value
251  for (int j = lo.y; j <= hi.y; ++j) {
252  AMREX_PRAGMA_SIMD
253  for (int i = lo.x; i <= hi.x; ++i) {
254  coeffA_a(i,j,lo.z) = 0.0;
255  coeffB_a(i,j,lo.z) = 1.0;
256  coeffC_a(i,j,lo.z) = 0.0;
257  }
258  }
259  for (int j = lo.y; j <= hi.y; ++j) {
260  AMREX_PRAGMA_SIMD
261  for (int i = lo.x; i <= hi.x; ++i) {
262 
263  // If at the top of the grid, we will set w to a specified Dirichlet value
264  coeffA_a(i,j,hi.z+1) = 0.0;
265  coeffB_a(i,j,hi.z+1) = 1.0;
266  coeffC_a(i,j,hi.z+1) = 0.0;
267 
268  // UNLESS if at the top of the domain and the boundary is outflow,
269  // we will use a homogeneous Neumann condition
270  if ( (hi.z == domhi.z) &&
271  (phys_bc_type[5] == ERF_BC::outflow or phys_bc_type[5] == ERF_BC::ho_outflow) )
272  {
273  coeffA_a(i,j,hi.z+1) = -1.0;
274  }
275  }
276  }
277  for (int k = lo.z+1; k <= hi.z+1; ++k) {
278  for (int j = lo.y; j <= hi.y; ++j) {
279  AMREX_PRAGMA_SIMD
280  for (int i = lo.x; i <= hi.x; ++i) {
281  gam_a(i,j,k) = coeffC_a(i,j,k-1) / coeffB_a(i,j,k-1);
282  Real bet = coeffB_a(i,j,k) - coeffA_a(i,j,k)*gam_a(i,j,k);
283  coeffB_a(i,j,k) = bet;
284  }
285  }
286  }
287 #endif
288  } // end profile
289 
290  // In the end we save the inverse of the diagonal (B) coefficient
291  {
292  BL_PROFILE("make_coeffs_invert");
293  ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (int i, int j, int k)
294  {
295  coeffB_a(i,j,k) = 1.0 / coeffB_a(i,j,k);
296  });
297  } // end profile
298  } // mfi
299  } // omp
300 }
R_v
constexpr amrex::Real R_v
Definition: ERF_Constants.H:11
R_d
constexpr amrex::Real R_d
Definition: ERF_Constants.H:10
Gamma
constexpr amrex::Real Gamma
Definition: ERF_Constants.H:19
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
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:158
Here is the call graph for this function: