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 
46  const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
47  Real dzi = dxInv[2];
48 
49  const Box &domain = geom.Domain();
50 
51  MultiFab coeff_A_mf(fast_coeffs, amrex::make_alias, 0, 1);
52  MultiFab coeff_B_mf(fast_coeffs, amrex::make_alias, 1, 1);
53  MultiFab coeff_C_mf(fast_coeffs, amrex::make_alias, 2, 1);
54  MultiFab coeff_P_mf(fast_coeffs, amrex::make_alias, 3, 1);
55  MultiFab coeff_Q_mf(fast_coeffs, amrex::make_alias, 4, 1);
56 
57 
58  // *************************************************************************
59  // Set gravity as a vector
60  const Array<Real,AMREX_SPACEDIM> grav{0.0, 0.0, -gravity};
61  const GpuArray<Real,AMREX_SPACEDIM> grav_gpu{grav[0], grav[1], grav[2]};
62 
63  // *************************************************************************
64  // Define updates in the current RK stage
65  // *************************************************************************
66 #ifdef _OPENMP
67 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
68 #endif
69  {
70 
71  for ( MFIter mfi(S_stage_data[IntVars::cons],TileNoZ()); mfi.isValid(); ++mfi)
72  {
73  Box bx = mfi.tilebox();
74  Box tbz = surroundingNodes(bx,2);
75 
76  const Array4<const Real> & stage_cons = S_stage_data[IntVars::cons].const_array(mfi);
77  const Array4<const Real> & prim = S_stage_prim.const_array(mfi);
78 
79  const Array4<const Real>& detJ = (mesh_type != MeshType::ConstantDz) ? detJ_cc->const_array(mfi) : Array4<const Real>{};
80 
81  const Array4<const Real>& r0_ca = r0->const_array(mfi);
82  const Array4<const Real>& pi0_ca = pi0->const_array(mfi);
83  const Array4<const Real>& pi_stage_ca = pi_stage.const_array(mfi);
84 
85  FArrayBox gam_fab; gam_fab.resize(surroundingNodes(bx,2),1,The_Async_Arena());
86 
87  auto const& coeffA_a = coeff_A_mf.array(mfi);
88  auto const& coeffB_a = coeff_B_mf.array(mfi);
89  auto const& coeffC_a = coeff_C_mf.array(mfi);
90  auto const& coeffP_a = coeff_P_mf.array(mfi);
91  auto const& coeffQ_a = coeff_Q_mf.array(mfi);
92  auto const& gam_a = gam_fab.array();
93 
94  // *********************************************************************
95  // *********************************************************************
96  // *********************************************************************
97 
98  Box bx_shrunk_in_k = bx;
99  int klo = tbz.smallEnd(2);
100  int khi = tbz.bigEnd(2);
101  bx_shrunk_in_k.setSmall(2,klo+1);
102  bx_shrunk_in_k.setBig(2,khi-1);
103 
104  // Note that the notes use "g" to mean the magnitude of gravity, so it is positive
105  // We set grav_gpu[2] to be the vector component which is negative
106  // We define halfg to match the notes (which is why we take the absolute value)
107  Real halfg = std::abs(0.5 * grav_gpu[2]);
108 
109  //Note we don't act on the bottom or top boundaries of the domain
110  if (mesh_type != MeshType::ConstantDz)
111  {
112  ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (int i, int j, int k)
113  {
114  Real rhobar_lo, rhobar_hi, pibar_lo, pibar_hi;
115  rhobar_lo = r0_ca(i,j,k-1);
116  rhobar_hi = r0_ca(i,j,k );
117  pibar_lo = pi0_ca(i,j,k-1);
118  pibar_hi = pi0_ca(i,j,k );
119 
120  Real pi_c = 0.5 * (pi_stage_ca(i,j,k-1) + pi_stage_ca(i,j,k));
121 
122  Real detJ_on_kface = 0.5 * (detJ(i,j,k) + detJ(i,j,k-1));
123  Real inv_detJ_on_kface = 1. / detJ_on_kface;
124 
125  Real coeff_P = -Gamma * R_d * dzi * inv_detJ_on_kface
126  + halfg * R_d * rhobar_hi /
127  ( c_v * pibar_hi * stage_cons(i,j,k,RhoTheta_comp) );
128  coeff_P *= pi_c;
129 
130  Real coeff_Q = Gamma * R_d * dzi * inv_detJ_on_kface
131  + halfg * R_d * rhobar_lo /
132  ( c_v * pibar_lo * stage_cons(i,j,k-1,RhoTheta_comp) );
133  coeff_Q *= pi_c;
134 
135  coeffP_a(i,j,k) = coeff_P;
136  coeffQ_a(i,j,k) = coeff_Q;
137 
138  if (l_use_moisture) {
139  Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp)
140  +prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) );
141  coeff_P /= (1.0 + q);
142  coeff_Q /= (1.0 + q);
143  }
144 
145  Real theta_t_lo = 0.5 * ( prim(i,j,k-2,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp) );
146  Real theta_t_mid = 0.5 * ( prim(i,j,k-1,PrimTheta_comp) + prim(i,j,k ,PrimTheta_comp) );
147  Real theta_t_hi = 0.5 * ( prim(i,j,k ,PrimTheta_comp) + prim(i,j,k+1,PrimTheta_comp) );
148 
149  // LHS for tri-diagonal system
150  Real D = dtau * dtau * beta_2 * beta_2 * dzi;
151  coeffA_a(i,j,k) = D * ( halfg - coeff_Q * theta_t_lo );
152  coeffC_a(i,j,k) = D * (-halfg + coeff_P * theta_t_hi );
153 
154  coeffB_a(i,j,k) = detJ_on_kface + D * (coeff_Q - coeff_P) * theta_t_mid;
155  });
156 
157  } else {
158 
159  ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (int i, int j, int k)
160  {
161  Real rhobar_lo, rhobar_hi, pibar_lo, pibar_hi;
162  rhobar_lo = r0_ca(i,j,k-1);
163  rhobar_hi = r0_ca(i,j,k );
164  pibar_lo = pi0_ca(i,j,k-1);
165  pibar_hi = pi0_ca(i,j,k );
166 
167  Real pi_c = 0.5 * (pi_stage_ca(i,j,k-1) + pi_stage_ca(i,j,k));
168 
169  Real coeff_P = -Gamma * R_d * dzi
170  + halfg * R_d * rhobar_hi /
171  ( c_v * pibar_hi * stage_cons(i,j,k,RhoTheta_comp) );
172  coeff_P *= pi_c;
173 
174  Real coeff_Q = Gamma * R_d * dzi
175  + halfg * R_d * rhobar_lo /
176  ( c_v * pibar_lo * stage_cons(i,j,k-1,RhoTheta_comp) );
177  coeff_Q *= pi_c;
178 
179  coeffP_a(i,j,k) = coeff_P;
180  coeffQ_a(i,j,k) = coeff_Q;
181 
182  if (l_use_moisture) {
183  Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp)
184  +prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) );
185  coeff_P /= (1.0 + q);
186  coeff_Q /= (1.0 + q);
187  }
188 
189  Real theta_t_lo = 0.5 * ( prim(i,j,k-2,PrimTheta_comp) + prim(i,j,k-1,PrimTheta_comp) );
190  Real theta_t_mid = 0.5 * ( prim(i,j,k-1,PrimTheta_comp) + prim(i,j,k ,PrimTheta_comp) );
191  Real theta_t_hi = 0.5 * ( prim(i,j,k ,PrimTheta_comp) + prim(i,j,k+1,PrimTheta_comp) );
192 
193  // LHS for tri-diagonal system
194  Real D = dtau * dtau * beta_2 * beta_2 * dzi;
195  coeffA_a(i,j,k) = D * ( halfg - coeff_Q * theta_t_lo );
196  coeffC_a(i,j,k) = D * (-halfg + coeff_P * theta_t_hi );
197 
198  coeffB_a(i,j,k) = 1.0 + D * (coeff_Q - coeff_P) * theta_t_mid;
199  });
200  }
201 
202  amrex::Box b2d = tbz; // Copy constructor
203  b2d.setRange(2,0);
204 
205  auto const lo = amrex::lbound(bx);
206  auto const hi = amrex::ubound(bx);
207 
208  auto const domhi = amrex::ubound(domain);
209 
210  {
211  BL_PROFILE("make_coeffs_b2d_loop");
212 #ifdef AMREX_USE_GPU
213  ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int) {
214 
215  // If at the bottom of the grid, we will set w to a specified Dirichlet value
216  coeffA_a(i,j,lo.z) = 0.0;
217  coeffB_a(i,j,lo.z) = 1.0;
218  coeffC_a(i,j,lo.z) = 0.0;
219 
220  // If at the top of the grid, we will set w to a specified Dirichlet value
221  coeffA_a(i,j,hi.z+1) = 0.0;
222  coeffB_a(i,j,hi.z+1) = 1.0;
223  coeffC_a(i,j,hi.z+1) = 0.0;
224 
225  // UNLESS if at the top of the domain and the boundary is outflow,
226  // we will use a homogeneous Neumann condition
227  if ( (hi.z == domhi.z) &&
228  (phys_bc_type[5] == ERF_BC::outflow or phys_bc_type[5] == ERF_BC::ho_outflow) )
229  {
230  coeffA_a(i,j,hi.z+1) = -1.0;
231  }
232 
233  // w = specified Dirichlet value at k = lo.z
234  Real bet = coeffB_a(i,j,lo.z);
235 
236  for (int k = lo.z+1; k <= hi.z+1; k++) {
237  gam_a(i,j,k) = coeffC_a(i,j,k-1) / bet;
238  bet = coeffB_a(i,j,k) - coeffA_a(i,j,k)*gam_a(i,j,k);
239  coeffB_a(i,j,k) = bet;
240  }
241  });
242 #else
243  // If at the bottom of the grid, we will set w to a specified Dirichlet value
244  for (int j = lo.y; j <= hi.y; ++j) {
245  AMREX_PRAGMA_SIMD
246  for (int i = lo.x; i <= hi.x; ++i) {
247  coeffA_a(i,j,lo.z) = 0.0;
248  coeffB_a(i,j,lo.z) = 1.0;
249  coeffC_a(i,j,lo.z) = 0.0;
250  }
251  }
252  for (int j = lo.y; j <= hi.y; ++j) {
253  AMREX_PRAGMA_SIMD
254  for (int i = lo.x; i <= hi.x; ++i) {
255 
256  // If at the top of the grid, we will set w to a specified Dirichlet value
257  coeffA_a(i,j,hi.z+1) = 0.0;
258  coeffB_a(i,j,hi.z+1) = 1.0;
259  coeffC_a(i,j,hi.z+1) = 0.0;
260 
261  // UNLESS if at the top of the domain and the boundary is outflow,
262  // we will use a homogeneous Neumann condition
263  if ( (hi.z == domhi.z) &&
264  (phys_bc_type[5] == ERF_BC::outflow or phys_bc_type[5] == ERF_BC::ho_outflow) )
265  {
266  coeffA_a(i,j,hi.z+1) = -1.0;
267  }
268  }
269  }
270  for (int k = lo.z+1; k <= hi.z+1; ++k) {
271  for (int j = lo.y; j <= hi.y; ++j) {
272  AMREX_PRAGMA_SIMD
273  for (int i = lo.x; i <= hi.x; ++i) {
274  gam_a(i,j,k) = coeffC_a(i,j,k-1) / coeffB_a(i,j,k-1);
275  Real bet = coeffB_a(i,j,k) - coeffA_a(i,j,k)*gam_a(i,j,k);
276  coeffB_a(i,j,k) = bet;
277  }
278  }
279  }
280 #endif
281  } // end profile
282 
283  // In the end we save the inverse of the diagonal (B) coefficient
284  {
285  BL_PROFILE("make_coeffs_invert");
286  ParallelFor(bx_shrunk_in_k, [=] AMREX_GPU_DEVICE (int i, int j, int k)
287  {
288  coeffB_a(i,j,k) = 1.0 / coeffB_a(i,j,k);
289  });
290  } // end profile
291  } // mfi
292  } // omp
293 }
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
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:139
Here is the call graph for this function: