AMReX-Hydro
AMReX-based hydro routines for low Mach number flows
hydro_godunov_ppm.H
Go to the documentation of this file.
1 /*
2  * \file hydro_godunov_ppm.H
3  *
4  * \addtogroup Godunov
5  * @{
6  */
7 
8 #ifndef HYDRO_GODUNOV_PPM_H
9 #define HYDRO_GODUNOV_PPM_H
10 
11 #include <AMReX_MultiFab.H>
12 #include <AMReX_BCRec.H>
13 #include <AMReX_BC_TYPES.H>
14 #include <hydro_constants.H>
15 
17 void
18 minmod_fn (const amrex::Real sm1,
19  const amrex::Real s0,
20  const amrex::Real sp1,
21  amrex::Real& dsm,
22  amrex::Real& dsp)
23 {
24  // Calculate gradients on both sides
25  dsp = sp1 - s0;
26  dsm = s0 - sm1;
27 
28  if (!(dsp * dsm < 0.0)) {
29  // Select the smaller slope if same sign
30  if (std::abs(dsp) < std::abs(dsm)) {
31  dsm = dsp;
32  } else {
33  dsp = dsm;
34  }
35  } else {
36  // Set to zero if opposite sign
37  dsp = 0.0;
38  dsm = 0.0;
39  }
40 }
41 
43 void
44 Godunov_minmod_bc_lo (const int n,
45  const amrex::Real sm1,
46  const amrex::Real s0,
47  const amrex::Real vel_edge,
48  amrex::Real& dsm,
49  const int bclo,
50  const int domlo)
51 {
52  // Ensure that left-side slope is used unchanged
53  if (n == domlo) {
54  if ( bclo == amrex::BCType::ext_dir || bclo == amrex::BCType::hoextrap ||
55  (bclo == amrex::BCType::direction_dependent && vel_edge >= 0.0) )
56  {
57  dsm = s0 - sm1;
58  }
59  }
60 }
61 
63 void
64 Godunov_minmod_bc_hi (const int n,
65  const amrex::Real s0,
66  const amrex::Real sp1,
67  const amrex::Real vel_edge,
68  amrex::Real& dsp,
69  const int bchi,
70  const int domhi)
71 {
72  // Ensure that the right-side slope is used unchanged
73  if (n == domhi) {
74  if ( bchi == amrex::BCType::ext_dir || bchi == amrex::BCType::hoextrap ||
75  (bchi == amrex::BCType::direction_dependent && vel_edge <= 0.0) )
76  {
77  dsp = sp1 - s0;
78  }
79  }
80 }
81 
82 namespace PPM {
83 
85 
86 static constexpr int default_limiter = VanLeer;
87 
88 struct nolimiter {
89 
90  static constexpr bool do_limiting = false;
91 
93  static amrex::Real
94  sedge1(const amrex::Real sm2,
95  const amrex::Real sm1,
96  const amrex::Real s0,
97  const amrex::Real sp1,
98  const amrex::Real /*sp2*/)
99  {
100 
101  amrex::Real d1 = amrex::Real(0.5) * (sp1 - sm1);
102  amrex::Real d2 = amrex::Real(0.5) * (s0 - sm2);
103 
104  amrex::Real sedge = amrex::Real(0.5)*(s0 + sm1) - amrex::Real(1./6.)*(d1 - d2);
105  return sedge;
106  }
107 
109  static amrex::Real
110  sedge2(const amrex::Real /*sm2*/,
111  const amrex::Real sm1,
112  const amrex::Real s0,
113  const amrex::Real sp1,
114  const amrex::Real sp2)
115  {
116  amrex::Real d1 = amrex::Real(0.5) * (sp2 - s0);
117  amrex::Real d2 = amrex::Real(0.5) * (sp1 - sm1);
118 
119  amrex::Real sedge = amrex::Real(0.5)*(sp1 + s0) - amrex::Real(1./6.)*(d1 - d2);
120  return sedge;
121  }
122 
125  sm_sp(const amrex::Real /*s0*/,
126  const amrex::Real sedge1,
127  const amrex::Real sedge2)
128  {
129  return amrex::makeTuple(sedge1, sedge2);
130  }
131 }; // nolimiter
132 
133 struct vanleer {
134 
135  static constexpr bool do_limiting = true;
136 
138  static amrex::Real vanLeer(const amrex::Real a, const amrex::Real b, const amrex::Real c)
139  {
140  constexpr auto small_qty_sq = amrex::Real(1.e-10) * amrex::Real(1.e-10);
141 
142  amrex::Real dsc = amrex::Real(0.5)*(b - c);
143  amrex::Real dsl = amrex::Real(2.0)*(a - c);
144  amrex::Real dsr = amrex::Real(2.0)*(b - a);
145  return (dsl*dsr > small_qty_sq) ?
146  std::copysign(amrex::Real(1.0), dsc)*amrex::min(amrex::Math::abs(dsc),amrex::min(amrex::Math::abs(dsl),
147  amrex::Math::abs(dsr))) : amrex::Real(0.0);
148  }
149 
151  static amrex::Real
152  sedge1(const amrex::Real sm2,
153  const amrex::Real sm1,
154  const amrex::Real s0,
155  const amrex::Real sp1,
156  const amrex::Real /*sp2*/)
157  {
158 
159  amrex::Real d1 = vanLeer(s0,sp1,sm1);
160  amrex::Real d2 = vanLeer(sm1,s0,sm2);
161 
162  amrex::Real sedge = amrex::Real(0.5)*(s0 + sm1) - amrex::Real(1./6.)*(d1 - d2);
163  return amrex::min(amrex::max(sedge, amrex::min(s0, sm1)),amrex::max(s0,sm1));
164 
165  }
166 
168  static amrex::Real
169  sedge2(const amrex::Real /*sm2*/,
170  const amrex::Real sm1,
171  const amrex::Real s0,
172  const amrex::Real sp1,
173  const amrex::Real sp2)
174  {
175 
176  amrex::Real d1 = vanLeer(sp1,sp2,s0);
177  amrex::Real d2 = vanLeer(s0,sp1,sm1);
178 
179  amrex::Real sedge = amrex::Real(0.5)*(sp1 + s0) - amrex::Real(1./6.)*(d1 - d2);
180  return amrex::min(amrex::max(sedge, amrex::min(s0, sp1)),amrex::max(s0,sp1));
181 
182 
183  }
184 
187  sm_sp(const amrex::Real s0,
188  const amrex::Real sedge1,
189  const amrex::Real sedge2)
190  {
191  amrex::Real sm_ = sedge1;
192  amrex::Real sp_ = sedge2;
193 
194  if ((sedge2-s0)*(s0-sedge1) < 0.0) {
195  sm_ = s0;
196  sp_ = s0;
197  } else if (amrex::Math::abs(sedge2-s0) >= amrex::Real(2.0)*amrex::Math::abs(sedge1-s0)) {
198  sp_ = amrex::Real(3.0)*s0 - amrex::Real(2.0)*sedge1;
199  } else if (amrex::Math::abs(sedge1-s0) >= amrex::Real(2.0)*amrex::Math::abs(sedge2-s0)) {
200  sm_ = amrex::Real(3.0)*s0 - amrex::Real(2.0)*sedge2;
201  }
202  return amrex::makeTuple(sm_,sp_);
203 
204  // this might be the correct way to do this in case the else if are both == 0.0
205  // could also change the >= to just > for each
206 // if ((sedge2-s0)*(s0-sedge1) < amrex::Real(0.0)) {
207 // return amrex::makeTuple(s0,s0);
208 // }
209 // amrex::Real sm_ = sedge1;
210 // amrex::Real sp_ = sedge2;
211 // if (amrex::Math::abs(sedge2-s0) >= amrex::Real(2.0)*amrex::Math::abs(sedge1-s0)) {
212 // sp_ = amrex::Real(3.0)*s0 - amrex::Real(2.0)*sedge1;
213 // }
214 // if (amrex::Math::abs(sedge1-s0) >= amrex::Real(2.0)*amrex::Math::abs(sedge2-s0)) {
215 // sm_ = amrex::Real(3.0)*s0 - amrex::Real(2.0)*sedge2;
216 // }
217 // return amrex::makeTuple(sm_,sp_);
218 
219  }
220 }; // vanleer
221 
222 struct wenoz {
223 
224  static constexpr bool do_limiting = false;
225 
227  static amrex::Real
228  sedge2(const amrex::Real sm2,
229  const amrex::Real sm1,
230  const amrex::Real s0,
231  const amrex::Real sp1,
232  const amrex::Real sp2)
233  {
234  constexpr auto eps = amrex::Real(1.0e-6);
235 
236  const amrex::Real beta1 =
237  amrex::Real(13.0) / amrex::Real(12.0) * (sm2 - amrex::Real(2.0) * sm1 + s0) * (sm2 - amrex::Real(2.0) * sm1 + s0) +
238  amrex::Real(0.25) * (sm2 - amrex::Real(4.0) * sm1 + amrex::Real(3.0) * s0) * (sm2 - amrex::Real(4.0) * sm1 + amrex::Real(3.0) * s0);
239  const amrex::Real beta2 =
240  amrex::Real(13.0) / amrex::Real(12.0) * (sm1 - amrex::Real(2.0) * s0 + sp1) * (sm1 - amrex::Real(2.0) * s0 + sp1) +
241  amrex::Real(0.25) * (sm1 - sp1) * (sm1 - sp1);
242  const amrex::Real beta3 =
243  amrex::Real(13.0) / amrex::Real(12.0) * (s0 - amrex::Real(2.0) * sp1 + sp2) * (s0 - amrex::Real(2.0) * sp1 + sp2) +
244  amrex::Real(0.25) * (amrex::Real(3.0) * s0 - amrex::Real(4.0) * sp1 + sp2) * (amrex::Real(3.0) * s0 - amrex::Real(4.0) * sp1 + sp2);
245 
246  const amrex::Real t5 = amrex::Math::abs(beta3 - beta1);
247  const amrex::Real omega1 = amrex::Real(0.1) * (amrex::Real(1.0) + t5 / (eps + beta1));
248  const amrex::Real omega2 = amrex::Real(0.6) * (amrex::Real(1.0) + t5 / (eps + beta2));
249  const amrex::Real omega3 = amrex::Real(0.3) * (amrex::Real(1.0) + t5 / (eps + beta3));
250 
251  const amrex::Real omega = omega1 + omega2 + omega3;
252 
253  const amrex::Real v_1 = amrex::Real(2.0) * sm2 - amrex::Real(7.0) * sm1 + amrex::Real(11.0) * s0;
254  const amrex::Real v_2 = -sm1 + amrex::Real(5.0) * s0 + amrex::Real(2.0) * sp1;
255  const amrex::Real v_3 = amrex::Real(2.0) * s0 + amrex::Real(5.0) * sp1 - sp2;
256 
257  return (omega1 * v_1 + omega2 * v_2 + omega3 * v_3) / (amrex::Real(6.0) * omega);
258  }
259 
261  static amrex::Real
262  sedge1(const amrex::Real sm2,
263  const amrex::Real sm1,
264  const amrex::Real s0,
265  const amrex::Real sp1,
266  const amrex::Real sp2)
267  {
268  return sedge2(sp2,sp1,s0,sm1,sm2); // NOLINT(readability-suspicious-call-argument)
269  }
270 
273  sm_sp(const amrex::Real /*s0*/,
274  const amrex::Real sedge1,
275  const amrex::Real sedge2)
276  {
277  return amrex::makeTuple(sedge1, sedge2);
278  }
279 }; // wenoz
280 
281 struct weno_js {
282 
283  static constexpr bool do_limiting = false;
284 
286  static amrex::Real
287  sedge2(const amrex::Real sm2,
288  const amrex::Real sm1,
289  const amrex::Real s0,
290  const amrex::Real sp1,
291  const amrex::Real sp2)
292  {
293  constexpr auto eps = amrex::Real(1.0e-6);
294 
295  const amrex::Real beta1 =
296  amrex::Real(13.0) / amrex::Real(12.0) * (sm2 - amrex::Real(2.0) * sm1 + s0) *
297  (sm2 - amrex::Real(2.0) * sm1 + s0) +
298  amrex::Real(0.25) * (sm2 - amrex::Real(4.0) * sm1 + amrex::Real(3.0) * s0) *
299  (sm2 - amrex::Real(4.0) * sm1 + amrex::Real(3.0) * s0);
300  const amrex::Real beta2 =
301  amrex::Real(13.0) / amrex::Real(12.0) * (sm1 - amrex::Real(2.0) * s0 + sp1) *
302  (sm1 - amrex::Real(2.0) * s0 + sp1) +
303  amrex::Real(0.25) * (sm1 - sp1) * (sm1 - sp1);
304  const amrex::Real beta3 =
305  amrex::Real(13.0) / amrex::Real(12.0) * (s0 - amrex::Real(2.0) * sp1 + sp2) *
306  (s0 - amrex::Real(2.0) * sp1 + sp2) +
307  amrex::Real(0.25) * (amrex::Real(3.0) * s0 - amrex::Real(4.0) * sp1 + sp2) *
308  (amrex::Real(3.0) * s0 - amrex::Real(4.0) * sp1 + sp2);
309 
310  const amrex::Real omega1 = amrex::Real(0.1) / (eps + beta1);
311  const amrex::Real omega2 = amrex::Real(0.6) / (eps + beta2);
312  const amrex::Real omega3 = amrex::Real(0.3) / (eps + beta3);
313 
314  const amrex::Real omega = omega1 + omega2 + omega3;
315 
316  const amrex::Real v_1 = amrex::Real(2.0) * sm2 - amrex::Real(7.0) * sm1 + amrex::Real(11.0) * s0;
317  const amrex::Real v_2 = -sm1 + amrex::Real(5.0) * s0 + amrex::Real(2.0) * sp1;
318  const amrex::Real v_3 = amrex::Real(2.0) * s0 + amrex::Real(5.0) * sp1 - sp2;
319 
320  return (omega1 * v_1 + omega2 * v_2 + omega3 * v_3) / (amrex::Real(6.0) * omega);
321  }
322 
324  static amrex::Real
325  sedge1(const amrex::Real sm2,
326  const amrex::Real sm1,
327  const amrex::Real s0,
328  const amrex::Real sp1,
329  const amrex::Real sp2)
330  {
331  return sedge2(sp2,sp1,s0,sm1,sm2); // NOLINT(readability-suspicious-call-argument)
332  }
333 
336  sm_sp(const amrex::Real /*s0*/,
337  const amrex::Real sedge1,
338  const amrex::Real sedge2)
339  {
340  return amrex::makeTuple(sedge1, sedge2);
341  }
342 }; // weno_js
343 
344 struct upwind {
345 
346  static constexpr bool do_limiting = false;
347 
349  static amrex::Real
350  sedge1(const amrex::Real /*sm2*/,
351  const amrex::Real /*sm1*/,
352  const amrex::Real /*s0*/,
353  const amrex::Real /*sp1*/,
354  const amrex::Real /*sp2*/)
355  {
356  return 0;
357  }
358 
360  static amrex::Real
361  sedge2(const amrex::Real /*sm2*/,
362  const amrex::Real /*sm1*/,
363  const amrex::Real /*s0*/,
364  const amrex::Real /*sp1*/,
365  const amrex::Real /*sp2*/)
366  {
367  return 0;
368  }
369 
372  sm_sp(const amrex::Real /*s0*/,
373  const amrex::Real sedge1,
374  const amrex::Real sedge2)
375  {
376  return amrex::makeTuple(sedge1, sedge2);
377  }
378 }; // upwind
379 
380 struct minmod {
381 
382  static constexpr bool do_limiting = false;
383 
385  static amrex::Real
386  sedge1(const amrex::Real /*sm2*/,
387  const amrex::Real /*sm1*/,
388  const amrex::Real /*s0*/,
389  const amrex::Real /*sp1*/,
390  const amrex::Real /*sp2*/)
391  {
392  return 0;
393  }
394 
396  static amrex::Real
397  sedge2(const amrex::Real /*sm2*/,
398  const amrex::Real /*sm1*/,
399  const amrex::Real /*s0*/,
400  const amrex::Real /*sp1*/,
401  const amrex::Real /*sp2*/)
402  {
403  return 0;
404  }
405 
408  sm_sp(const amrex::Real /*s0*/,
409  const amrex::Real sedge1,
410  const amrex::Real sedge2)
411  {
412  return amrex::makeTuple(sedge1, sedge2);
413  }
414 }; // minmod
415 
416 template <typename Limiter>
418 void SetXBCs ( const int i, const int j, const int k, const int n,
419  amrex::Real &sm, amrex::Real &sp,
420  amrex::Real &sedge1, amrex::Real &sedge2,
422  const amrex::Real velm, const amrex::Real velp,
423  const int bclo, const int bchi,
424  const int domlo, const int domhi)
425 {
426  using namespace amrex;
427 
428  if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) ||
429  (bclo == BCType::direction_dependent && velp >= 0.0) )
430  {
431  if (i == domlo)
432  {
433  sedge2 = -amrex::Real(0.2)*s(domlo-1,j,k,n) + amrex::Real(0.75)*s(domlo,j, k, n)
434  +amrex::Real(0.5)*s(domlo+1,j,k,n) - amrex::Real(0.05)*s(domlo+2,j,k,n);
435  if constexpr (Limiter::do_limiting) {
436  sedge2 = amrex::max(sedge2, amrex::min(s(domlo+1,j,k,n), s(domlo,j,k,n)));
437  sedge2 = amrex::min(sedge2, amrex::max(s(domlo+1,j,k,n), s(domlo,j,k,n)));
438  }
439 
440  sm = s(domlo-1,j,k,n);
441  sp = sedge2;
442 
443  } else if (i == domlo+1) {
444 
445  sedge1 = -amrex::Real(0.2)*s(domlo-1,j,k,n) + amrex::Real(0.75)*s(domlo ,j,k,n)
446  + amrex::Real(0.5)*s(domlo+1,j,k,n) - amrex::Real(0.05)*s(domlo+2,j,k,n);
447  if constexpr (Limiter::do_limiting) {
448  sedge1 = amrex::max(sedge1, amrex::min(s(domlo+1,j,k,n), s(domlo,j,k,n)));
449  sedge1 = amrex::min(sedge1, amrex::max(s(domlo+1,j,k,n), s(domlo,j,k,n)));
450  }
451 
452  sp = sedge2;
453  sm = sedge1;
454 
455  if constexpr (Limiter::do_limiting) {
456  if ( (sp - s(domlo+1,j,k,n))*(s(domlo+1,j,k,n) - sm) <= amrex::Real(0.0))
457  {
458  sp = s(domlo+1,j,k,n);
459  sm = s(domlo+1,j,k,n);
460  } else if(amrex::Math::abs(sp - s(domlo+1,j,k,n)) >= amrex::Real(2.0)*amrex::Math::abs(sm - s(domlo+1,j,k,n))) {
461  sp = amrex::Real(3.0)*s(domlo+1,j,k,n) - amrex::Real(2.0)*sm;
462  } else if(amrex::Math::abs(sm - s(domlo+1,j,k,n)) >= amrex::Real(2.0)*amrex::Math::abs(sp - s(domlo+1,j,k,n))) {
463  sm = amrex::Real(3.0)*s(domlo+1,j,k,n) - amrex::Real(2.0)*sp;
464  }
465  }
466  }
467  }
468 
469  if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) ||
470  (bchi == BCType::direction_dependent && velm <= 0.0) )
471  {
472  if (i == domhi)
473  {
474  sedge1 = -amrex::Real(0.2)*s(domhi+1,j,k,n) + amrex::Real(0.75)*s(domhi,j,k, n)
475  +amrex::Real(0.5)*s(domhi-1,j,k,n) - amrex::Real(0.05)*s(domhi-2,j,k,n);
476  if constexpr (Limiter::do_limiting) {
477  sedge1 = amrex::max(sedge1, amrex::min(s(domhi-1,j,k,n), s(domhi,j,k,n)));
478  sedge1 = amrex::min(sedge1, amrex::max(s(domhi-1,j,k,n), s(domhi,j,k,n)));
479  }
480 
481  sp = s(domhi+1,j,k,n);
482  sm = sedge1;
483 
484  } else if (i == domhi-1) {
485 
486  sedge2 = -amrex::Real(0.2)*s(domhi+1,j,k,n) + amrex::Real(0.75)*s(domhi ,j,k,n)
487  +amrex::Real(0.5)*s(domhi-1,j,k,n) - amrex::Real(0.05)*s(domhi-2,j,k,n);
488  if constexpr (Limiter::do_limiting) {
489  sedge2 = amrex::max(sedge2, amrex::min(s(domhi-1,j,k,n), s(domhi,j,k,n)));
490  sedge2 = amrex::min(sedge2, amrex::max(s(domhi-1,j,k,n), s(domhi,j,k,n)));
491  }
492 
493  sp = sedge2;
494  sm = sedge1;
495 
496  if constexpr (Limiter::do_limiting) {
497  if( (sp - s(domhi-1,j,k,n))*(s(domhi-1,j,k,n) - sm) <= amrex::Real(0.0))
498  {
499  sp = s(domhi-1,j,k,n);
500  sm = s(domhi-1,j,k,n);
501  } else if(amrex::Math::abs(sp - s(domhi-1,j,k,n)) >= 2.*amrex::Math::abs(sm - s(domhi-1,j,k,n))) {
502  sp = amrex::Real(3.0)*s(domhi-1,j,k,n) - amrex::Real(2.0)*sm;
503  } else if(amrex::Math::abs(sm - s(domhi-1,j,k,n)) >= 2.*amrex::Math::abs(sp - s(domhi-1,j,k,n))) {
504  sm = amrex::Real(3.0)*s(domhi-1,j,k,n) - amrex::Real(2.0)*sp;
505  }
506  }
507  }
508  }
509 }
510 
511 template <typename Limiter>
513 void SetYBCs ( const int i, const int j, const int k, const int n,
514  amrex::Real &sm, amrex::Real &sp,
515  amrex::Real &sedge1, amrex::Real &sedge2,
517  const amrex::Real velm, const amrex::Real velp,
518  const int bclo, const int bchi,
519  const int domlo, const int domhi)
520 {
521  using namespace amrex;
522 
523  if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) ||
524  (bclo == BCType::direction_dependent && velm >= 0.0) )
525  {
526  if (j == domlo)
527  {
528  sedge2 = -amrex::Real(0.2)*s(i,domlo-1,k,n) + amrex::Real(0.75)*s(i,domlo ,k,n)
529  +amrex::Real(0.5)*s(i,domlo+1,k,n) - amrex::Real(0.05)*s(i,domlo+2,k,n);
530  if constexpr (Limiter::do_limiting) {
531  sedge2 = amrex::max(sedge2, amrex::min(s(i,domlo+1,k,n), s(i,domlo,k,n)));
532  sedge2 = amrex::min(sedge2, amrex::max(s(i,domlo+1,k,n), s(i,domlo,k,n)));
533  }
534 
535  sm = s(i,domlo-1,k,n);
536  sp = sedge2;
537 
538  } else if (j == domlo+1) {
539 
540  sedge1 = -amrex::Real(0.2)*s(i,domlo-1,k,n) + amrex::Real(0.75)*s(i,domlo ,k,n)
541  +amrex::Real(0.5)*s(i,domlo+1,k,n) - amrex::Real(0.05)*s(i,domlo+2,k,n);
542  if constexpr (Limiter::do_limiting) {
543  sedge1 = amrex::max(sedge1, amrex::min(s(i,domlo+1,k,n), s(i,domlo,k,n)));
544  sedge1 = amrex::min(sedge1, amrex::max(s(i,domlo+1,k,n), s(i,domlo,k,n)));
545  }
546 
547  sp = sedge2;
548  sm = sedge1;
549 
550  if constexpr (Limiter::do_limiting) {
551  if ( (sp - s(i,domlo+1,k,n))*(s(i,domlo+1,k,n) - sm) <= amrex::Real(0.0))
552  {
553  sp = s(i,domlo+1,k,n);
554  sm = s(i,domlo+1,k,n);
555  } else if(amrex::Math::abs(sp - s(i,domlo+1,k,n)) >= amrex::Real(2.0)*amrex::Math::abs(sm - s(i,domlo+1,k,n))) {
556  sp = amrex::Real(3.0)*s(i,domlo+1,k,n) - amrex::Real(2.0)*sm;
557  } else if(amrex::Math::abs(sm - s(i,domlo+1,k,n)) >= amrex::Real(2.0)*amrex::Math::abs(sp - s(i,domlo+1,k,n))) {
558  sm = amrex::Real(3.0)*s(i,domlo+1,k,n) - amrex::Real(2.0)*sp;
559  }
560  }
561  }
562  }
563 
564  if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) ||
565  (bchi == BCType::direction_dependent && velp <= 0.0) )
566  {
567  if (j == domhi)
568  {
569  sedge1 = -amrex::Real(0.2)*s(i,domhi+1,k,n) + amrex::Real(0.75)*s(i,domhi ,k,n)
570  +amrex::Real(0.5)*s(i,domhi-1,k,n) - amrex::Real(0.05)*s(i,domhi-2,k,n);
571  if constexpr (Limiter::do_limiting) {
572  sedge1 = amrex::max(sedge1, amrex::min(s(i,domhi-1,k,n), s(i,domhi,k,n)));
573  sedge1 = amrex::min(sedge1, amrex::max(s(i,domhi-1,k,n), s(i,domhi,k,n)));
574  }
575 
576  sp = s(i,domhi+1,k, n);
577  sm = sedge1;
578 
579  } else if (j == domhi-1) {
580 
581  sedge2 = -amrex::Real(0.2)*s(i,domhi+1,k,n) + amrex::Real(0.75)*s(i,domhi ,k,n)
582  +amrex::Real(0.5)*s(i,domhi-1,k,n) - amrex::Real(0.05)*s(i,domhi-2,k,n);
583  if constexpr (Limiter::do_limiting) {
584  sedge2 = amrex::max(sedge2, amrex::min(s(i,domhi-1,k,n), s(i,domhi,k,n)));
585  sedge2 = amrex::min(sedge2, amrex::max(s(i,domhi-1,k,n), s(i,domhi,k,n)));
586  }
587 
588  sp = sedge2;
589  sm = sedge1;
590 
591  if constexpr (Limiter::do_limiting) {
592  if( (sp - s(i,domhi-1,k,n))*(s(i,domhi-1,k,n) - sm) <= amrex::Real(0.0)){
593  sp = s(i,domhi-1,k,n);
594  sm = s(i,domhi-1,k,n);
595  } else if(amrex::Math::abs(sp - s(i,domhi-1,k,n)) >= 2.*amrex::Math::abs(sm - s(i,domhi-1,k,n))) {
596  sp = amrex::Real(3.0)*s(i,domhi-1,k,n) - amrex::Real(2.0)*sm;
597  } else if(amrex::Math::abs(sm - s(i,domhi-1,k,n)) >= 2.*amrex::Math::abs(sp - s(i,domhi-1,k,n))) {
598  sm = amrex::Real(3.0)*s(i,domhi-1,k,n) - amrex::Real(2.0)*sp;
599  }
600  }
601  }
602  }
603 }
604 
605 #if (AMREX_SPACEDIM==3)
606 template <typename Limiter>
608 void SetZBCs ( const int i, const int j, const int k, const int n,
609  amrex::Real &sm, amrex::Real &sp,
610  amrex::Real &sedge1, amrex::Real &sedge2,
612  const amrex::Real velm, const amrex::Real velp,
613  const int bclo, const int bchi,
614  const int domlo, const int domhi)
615 {
616  using namespace amrex;
617 
618  if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) ||
619  (bclo == BCType::direction_dependent && velm >= 0.0) )
620  {
621  if (k == domlo)
622  {
623  sedge2 = -amrex::Real(0.2)*s(i,j,domlo-1,n) + amrex::Real(0.75)*s(i,j,domlo ,n)
624  +amrex::Real(0.5)*s(i,j,domlo+1,n) - amrex::Real(0.05)*s(i,j,domlo+2,n);
625  if constexpr (Limiter::do_limiting) {
626  sedge2 = amrex::max(sedge2, amrex::min(s(i,j,domlo+1,n), s(i,j,domlo,n)));
627  sedge2 = amrex::min(sedge2, amrex::max(s(i,j,domlo+1,n), s(i,j,domlo,n)));
628  }
629 
630  sm = s(i,j,domlo-1,n);
631  sp = sedge2;
632 
633  } else if (k == domlo+1) {
634 
635  sedge1 = -amrex::Real(0.2)*s(i,j,domlo-1,n) + amrex::Real(0.75)*s(i,j,domlo ,n)
636  +amrex::Real(0.5)*s(i,j,domlo+1,n) - amrex::Real(0.05)*s(i,j,domlo+2,n);
637  if constexpr (Limiter::do_limiting) {
638  sedge1 = amrex::max(sedge1, amrex::min(s(i,j,domlo+1,n), s(i,j,domlo,n)));
639  sedge1 = amrex::min(sedge1, amrex::max(s(i,j,domlo+1,n), s(i,j,domlo,n)));
640  }
641 
642  sp = sedge2;
643  sm = sedge1;
644 
645  if constexpr (Limiter::do_limiting) {
646  if ( (sp - s(i,j,domlo+1,n))*(s(i,j,domlo+1,n) - sm) <= 0. )
647  {
648  sp = s(i,j,domlo+1,n);
649  sm = s(i,j,domlo+1,n);
650  } else if(amrex::Math::abs(sp - s(i,j,domlo+1,n)) >= 2.*amrex::Math::abs(sm - s(i,j,domlo+1,n))) {
651  sp = amrex::Real(3.0)*s(i,j,domlo+1,n) - amrex::Real(2.0)*sm;
652  } else if(amrex::Math::abs(sm - s(i,j,domlo+1,n)) >= 2.*amrex::Math::abs(sp - s(i,j,domlo+1,n))) {
653  sm = amrex::Real(3.0)*s(i,j,domlo+1,n) - amrex::Real(2.0)*sp;
654  }
655  }
656  }
657  }
658 
659  if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) ||
660  (bchi == BCType::direction_dependent && velp <= 0.0) )
661  {
662  if (k == domhi)
663  {
664  sedge1 = -amrex::Real(0.2)*s(i,j,domhi+1,n) + amrex::Real(0.75)*s(i,j,domhi ,n)
665  +amrex::Real(0.5)*s(i,j,domhi-1,n) - amrex::Real(0.05)*s(i,j,domhi-2,n);
666  if constexpr (Limiter::do_limiting) {
667  sedge1 = amrex::max(sedge1, amrex::min(s(i,j,domhi-1,n), s(i,j,domhi,n)));
668  sedge1 = amrex::min(sedge1, amrex::max(s(i,j,domhi-1,n), s(i,j,domhi,n)));
669  }
670 
671  sp = s(i,j,domhi+1,n);
672  sm = sedge1;
673 
674  } else if (k == domhi-1) {
675 
676  sedge2 = -amrex::Real(0.2)*s(i,j,domhi+1,n) + amrex::Real(0.75)*s(i,j,domhi ,n)
677  +amrex::Real(0.5)*s(i,j,domhi-1,n) - amrex::Real(0.05)*s(i,j,domhi-2,n);
678  if constexpr (Limiter::do_limiting) {
679  sedge2 = amrex::max(sedge2, amrex::min(s(i,j,domhi-1,n), s(i,j,domhi,n)));
680  sedge2 = amrex::min(sedge2, amrex::max(s(i,j,domhi-1,n), s(i,j,domhi,n)));
681  }
682 
683  sp = sedge2;
684  sm = sedge1;
685 
686  if constexpr (Limiter::do_limiting) {
687  if ( (sp - s(i,j,domhi-1,n))*(s(i,j,domhi-1,n) - sm) <= 0. )
688  {
689  sp = s(i,j,domhi-1,n);
690  sm = s(i,j,domhi-1,n);
691  } else if(amrex::Math::abs(sp - s(i,j,domhi-1,n)) >= 2.*amrex::Math::abs(sm - s(i,j,domhi-1,n))) {
692  sp = amrex::Real(3.0)*s(i,j,domhi-1,n) - amrex::Real(2.0)*sm;
693  } else if(amrex::Math::abs(sm - s(i,j,domhi-1,n)) >= 2.*amrex::Math::abs(sp - s(i,j,domhi-1,n))) {
694  sm = amrex::Real(3.0)*s(i,j,domhi-1,n) - amrex::Real(2.0)*sp;
695  }
696  }
697 
698  }
699  }
700 }
701 #endif
702 
703 
704 // Right now only ppm type 1 is supported on GPU
705 // This version is called before the MAC projection, when we use the cell-centered velocity
706 // for upwinding
707 template <typename Limiter>
709 void PredictVelOnXFace ( const int i, const int j, const int k, const int n,
710  const amrex::Real dtdx, const amrex::Real v_ad,
712  const amrex::Array4<amrex::Real> &Im,
713  const amrex::Array4<amrex::Real> &Ip,
714  const amrex::BCRec bc, const int domlo, const int domhi,
715  const Limiter& /*limiter*/)
716 {
717  constexpr auto half{amrex::Real(0.5)};
718  constexpr auto one{amrex::Real(1.0)};
719  constexpr auto two3rds{amrex::Real(2.0/3.0)};
720 
721  amrex::Real sm2 = S(i-2,j,k,n);
722  amrex::Real sm1 = S(i-1,j,k,n);
723  amrex::Real s0 = S(i ,j,k,n);
724  amrex::Real sp1 = S(i+1,j,k,n);
725  amrex::Real sp2 = S(i+2,j,k,n);
726 
727  amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
728  amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
729 
730  amrex::Real sm, sp;
731  amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
732 
733  SetXBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, v_ad, v_ad,
734  bc.lo(0), bc.hi(0), domlo, domhi);
735 
736  amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp);
737 
738  amrex::Real sigma = amrex::Math::abs(v_ad)*dtdx;
739 
740  if (v_ad > small_vel)
741  {
742  Ip(i,j,k,n) = sp - (half*sigma)*((sp-sm) - (one - two3rds*sigma)*s6);
743  Im(i,j,k,n) = S(i,j,k,n);
744  }
745  else if (v_ad < -small_vel)
746  {
747  Ip(i,j,k,n) = S(i,j,k,n);
748  Im(i,j,k,n) = sm + (half*sigma)*((sp-sm) + (one - two3rds*sigma)*s6);
749  } else
750  {
751  Ip(i,j,k,n) = S(i,j,k,n);
752  Im(i,j,k,n) = S(i,j,k,n);
753  }
754 }
755 
756 template <typename Limiter>
758 void PredictVelOnYFace ( const int i, const int j, const int k, const int n,
759  const amrex::Real dtdy, const amrex::Real v_ad,
761  const amrex::Array4<amrex::Real> &Im,
762  const amrex::Array4<amrex::Real> &Ip,
763  const amrex::BCRec bc, const int domlo, const int domhi,
764  const Limiter& /*limiter*/)
765 {
766  constexpr auto half{amrex::Real(0.5)};
767  constexpr auto one{amrex::Real(1.0)};
768  constexpr auto two3rds{amrex::Real(2.0/3.0)};
769 
770  amrex::Real sm2 = S(i,j-2,k,n);
771  amrex::Real sm1 = S(i,j-1,k,n);
772  amrex::Real s0 = S(i,j ,k,n);
773  amrex::Real sp1 = S(i,j+1,k,n);
774  amrex::Real sp2 = S(i,j+2,k,n);
775 
776  amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
777  amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
778 
779  amrex::Real sm, sp;
780  amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
781 
782  SetYBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, v_ad, v_ad,
783  bc.lo(1), bc.hi(1), domlo, domhi);
784 
785  amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp);
786 
787  amrex::Real sigma = amrex::Math::abs(v_ad)*dtdy;
788 
789  if (v_ad > small_vel)
790  {
791  Ip(i,j,k,n) = sp - (half*sigma)*((sp-sm) - (one - two3rds*sigma)*s6);
792  Im(i,j,k,n) = S(i,j,k,n);
793  }
794  else if (v_ad < -small_vel)
795  {
796  Ip(i,j,k,n) = S(i,j,k,n);
797  Im(i,j,k,n) = sm + (half*sigma)*((sp-sm) + (one - two3rds*sigma)*s6);
798  } else
799  {
800  Ip(i,j,k,n) = S(i,j,k,n);
801  Im(i,j,k,n) = S(i,j,k,n);
802  }
803 }
804 
805 #if (AMREX_SPACEDIM==3)
806 template <typename Limiter>
808 void PredictVelOnZFace ( const int i, const int j, const int k, const int n,
809  const amrex::Real dtdz, const amrex::Real v_ad,
811  const amrex::Array4<amrex::Real> &Im,
812  const amrex::Array4<amrex::Real> &Ip,
813  const amrex::BCRec bc, const int domlo, const int domhi,
814  const Limiter& /*limiter*/)
815 {
816  constexpr auto half{amrex::Real(0.5)};
817  constexpr auto one{amrex::Real(1.0)};
818  constexpr auto two3rds{amrex::Real(2.0/3.0)};
819 
820  amrex::Real sm2 = S(i,j,k-2,n);
821  amrex::Real sm1 = S(i,j,k-1,n);
822  amrex::Real s0 = S(i,j,k ,n);
823  amrex::Real sp1 = S(i,j,k+1,n);
824  amrex::Real sp2 = S(i,j,k+2,n);
825 
826  amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
827  amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
828 
829  amrex::Real sm, sp;
830  amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
831 
832  SetZBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, v_ad, v_ad,
833  bc.lo(2), bc.hi(2), domlo, domhi);
834 
835  amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp);
836 
837  amrex::Real sigma = amrex::Math::abs(v_ad)*dtdz;
838 
839  if (v_ad > small_vel)
840  {
841  Ip(i,j,k,n) = sp - (half*sigma)*((sp-sm) - (one - two3rds*sigma)*s6);
842  Im(i,j,k,n) = S(i,j,k,n);
843  }
844  else if (v_ad < -small_vel)
845  {
846  Ip(i,j,k,n) = S(i,j,k,n);
847  Im(i,j,k,n) = sm + (half*sigma)*((sp-sm) + (one - two3rds*sigma)*s6);
848  } else
849  {
850  Ip(i,j,k,n) = S(i,j,k,n);
851  Im(i,j,k,n) = S(i,j,k,n);
852  }
853 }
854 
855 #endif
856 
857 // Right now only ppm type 1 is supported on GPU
858 // This version is called after the MAC projection, when we use the MAC-projected velocity
859 // for upwinding
860 template <typename Limiter>
862 void PredictStateOnXFace ( const int i, const int j, const int k, const int n,
863  const amrex::Real dt, const amrex::Real dx,
864  amrex::Real& Im, amrex::Real& Ip,
866  const amrex::Array4<const amrex::Real> &vel_edge,
867  const amrex::BCRec bc,
868  const int domlo, const int domhi,
869  const Limiter& /*limiter*/,
870  int limiter_type)
871 {
872  amrex::Real sm1 = S(i-1,j,k,n);
873  amrex::Real s0 = S(i ,j,k,n);
874  amrex::Real sp1 = S(i+1,j,k,n);
875 
876  constexpr auto half{amrex::Real(0.5)};
877  constexpr auto one{amrex::Real(1.0)};
878 
879  amrex::Real sigmap = amrex::Math::abs(vel_edge(i+1,j,k))*dt/dx;
880  amrex::Real sigmam = amrex::Math::abs(vel_edge(i ,j,k))*dt/dx;
881 
882  if (limiter_type == PPM::UPWIND) {
883  Im = S(i,j,k,n);
884  Ip = S(i,j,k,n);
885 
886  } else if (limiter_type == PPM::MINMOD) {
887  amrex::Real dsp = 0.0;
888  amrex::Real dsm = 0.0;
889  minmod_fn(sm1, s0, sp1, dsm, dsp);
890  Godunov_minmod_bc_lo(i, sm1, s0, vel_edge(i ,j,k), dsm, bc.lo(0), domlo);
891  Godunov_minmod_bc_hi(i, s0, sp1, vel_edge(i+1,j,k), dsp, bc.hi(0), domhi);
892 
893  if (vel_edge(i + 1, j, k) > small_vel) {
894  Ip = s0 + half * (one - sigmap) * dsp;
895  } else {
896  Ip = S(i,j,k,n);
897  }
898 
899  if (vel_edge(i, j, k) < -small_vel) {
900  Im = s0 - half * (one - sigmam) * dsm;
901  } else {
902  Im = S(i,j,k,n);
903  }
904 
905  } else {
906  constexpr auto two3rds{amrex::Real(2.0/3.0)};
907 
908  amrex::Real sm2 = S(i-2,j,k,n);
909  amrex::Real sp2 = S(i+2,j,k,n);
910 
911  amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
912  amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
913 
914  amrex::Real sm, sp;
915  amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
916 
917  SetXBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k),
918  bc.lo(0), bc.hi(0), domlo, domhi);
919 
920  amrex::Real s6 = 6.0*s0 - 3.0*(sm + sp);
921 
922  if (vel_edge(i+1,j,k) > small_vel) {
923  Ip = sp - (half*sigmap)*((sp - sm) - (one -two3rds*sigmap)*s6);
924  } else {
925  Ip = S(i,j,k,n);
926  }
927 
928  if(vel_edge(i,j,k) < -small_vel) {
929  Im = sm + (half*sigmam)*((sp-sm) + (one - two3rds*sigmam)*s6);
930  } else {
931  Im = S(i,j,k,n);
932  }
933  }
934 }
935 
936 template <typename Limiter>
938 void PredictStateOnYFace ( const int i, const int j, const int k, const int n,
939  const amrex::Real dt, const amrex::Real dx,
940  amrex::Real& Im, amrex::Real& Ip,
942  const amrex::Array4<const amrex::Real> &vel_edge,
943  const amrex::BCRec bc,
944  const int domlo, const int domhi,
945  const Limiter& /*limiter*/,
946  int limiter_type)
947 {
948  amrex::Real sm1 = S(i,j-1,k,n);
949  amrex::Real s0 = S(i,j ,k,n);
950  amrex::Real sp1 = S(i,j+1,k,n);
951 
952  constexpr auto half{amrex::Real(0.5)};
953  constexpr auto one{amrex::Real(1.0)};
954 
955  amrex::Real s6;
956 
957  amrex::Real sigmap = amrex::Math::abs(vel_edge(i,j+1,k))*dt/dx;
958  amrex::Real sigmam = amrex::Math::abs(vel_edge(i,j ,k))*dt/dx;
959 
960  if (limiter_type == PPM::UPWIND) {
961 
962  Im = S(i,j,k,n);
963  Ip = S(i,j,k,n);
964 
965  } else if (limiter_type == PPM::MINMOD) {
966 
967  amrex::Real dsp = 0.0;
968  amrex::Real dsm = 0.0;
969  minmod_fn(sm1, s0, sp1, dsm, dsp);
970  Godunov_minmod_bc_lo(j, sm1, s0, vel_edge(i,j ,k), dsm, bc.lo(1), domlo);
971  Godunov_minmod_bc_hi(j, s0, sp1, vel_edge(i,j+1,k), dsp, bc.hi(1), domhi);
972 
973  if (vel_edge(i,j+1,k) > small_vel) {
974  Ip = s0 + half * (one - sigmap) * dsp;
975  } else {
976  Ip = S(i,j,k,n);
977  }
978  if (vel_edge(i,j,k) < -small_vel) {
979  Im = s0 - half * (one - sigmam) * dsm;
980  } else {
981  Im = S(i,j,k,n);
982  }
983 
984  } else {
985 
986  constexpr auto two3rds{amrex::Real(2.0/3.0)};
987 
988  amrex::Real sm2 = S(i,j-2,k,n);
989  amrex::Real sp2 = S(i,j+2,k,n);
990 
991  amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
992  amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
993 
994  amrex::Real sm, sp;
995  amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
996 
997  SetYBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k),
998  bc.lo(1), bc.hi(1), domlo, domhi);
999 
1000  s6 = 6.0*s0- 3.0*(sm + sp);
1001 
1002  if (vel_edge(i,j+1,k) > small_vel) {
1003  Ip = sp - (half*sigmap)*((sp - sm) - (one -two3rds*sigmap)*s6);
1004  } else {
1005  Ip = S(i,j,k,n);
1006  }
1007 
1008  if (vel_edge(i,j,k) < -small_vel) {
1009  Im = sm + (half*sigmam)*((sp-sm) + (one - two3rds*sigmam)*s6);
1010  } else {
1011  Im = S(i,j,k,n);
1012  }
1013  } // not upwind or minmod
1014 }
1015 
1016 
1017 #if (AMREX_SPACEDIM==3)
1018 template <typename Limiter>
1020 void PredictStateOnZFace ( const int i, const int j, const int k, const int n,
1021  const amrex::Real dt, const amrex::Real dx,
1022  amrex::Real& Im, amrex::Real& Ip,
1024  const amrex::Array4<const amrex::Real> &vel_edge,
1025  const amrex::BCRec bc,
1026  const int domlo, const int domhi,
1027  const Limiter& /*limiter*/,
1028  int limiter_type)
1029 {
1030  amrex::Real sm1 = S(i,j,k-1,n);
1031  amrex::Real s0 = S(i,j,k ,n);
1032  amrex::Real sp1 = S(i,j,k+1,n);
1033 
1034  constexpr auto half{amrex::Real(0.5)};
1035  constexpr auto one{amrex::Real(1.0)};
1036 
1037  amrex::Real sigmap = amrex::Math::abs(vel_edge(i,j,k+1))*dt/dx;
1038  amrex::Real sigmam = amrex::Math::abs(vel_edge(i,j,k ))*dt/dx;
1039 
1040  if (limiter_type == PPM::UPWIND) {
1041 
1042  Im = S(i,j,k,n);
1043  Ip = S(i,j,k,n);
1044 
1045  } else if (limiter_type == PPM::MINMOD) {
1046 
1047  amrex::Real dsp = 0.0;
1048  amrex::Real dsm = 0.0;
1049  minmod_fn(sm1, s0, sp1, dsm, dsp);
1050  Godunov_minmod_bc_lo(k, sm1, s0, vel_edge(i,j,k ), dsm, bc.lo(2), domlo);
1051  Godunov_minmod_bc_hi(k, s0, sp1, vel_edge(i,j,k+1), dsp, bc.hi(2), domhi);
1052 
1053  if (vel_edge(i, j, k + 1) > small_vel) {
1054  Ip = s0 + half * (one - sigmap) * dsp;
1055  } else {
1056  Ip = S(i,j,k,n);
1057  }
1058 
1059  if (vel_edge(i, j, k) < -small_vel) {
1060  Im = s0 - half * (one - sigmam) * dsm;
1061  } else {
1062  Im = S(i,j,k,n);
1063  }
1064 
1065  } else {
1066  constexpr auto two3rds{amrex::Real(2.0/3.0)};
1067 
1068  amrex::Real sm2 = S(i,j,k-2,n);
1069  amrex::Real sp2 = S(i,j,k+2,n);
1070 
1071  amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
1072  amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
1073 
1074  amrex::Real sm, sp;
1075  amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
1076 
1077  SetZBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k),
1078  bc.lo(2), bc.hi(2), domlo, domhi);
1079 
1080  amrex::Real s6 = 6.0*s0- 3.0*(sm + sp);
1081 
1082  if(vel_edge(i,j,k+1) > small_vel) {
1083  Ip = sp - (half*sigmap)*((sp-sm) - (amrex::Real(1.0) -two3rds*sigmap)*s6);
1084  } else {
1085  Ip = S(i,j,k,n);
1086  }
1087 
1088  if(vel_edge(i,j,k) < -small_vel) {
1089  Im = sm + (half*sigmam)*((sp-sm) + (amrex::Real(1.0) - two3rds*sigmam)*s6);
1090  } else {
1091  Im = S(i,j,k,n);
1092  }
1093  }
1094 }
1095 #endif
1096 
1097 template <typename Limiter>
1100  amrex::Array4<amrex::Real> const& Imy,
1101  amrex::Array4<amrex::Real> const& Imz),
1103  amrex::Array4<amrex::Real> const& Ipy,
1104  amrex::Array4<amrex::Real> const& Ipz),
1107  amrex::Geometry geom,
1108  amrex::Real dt,
1109  amrex::BCRec const* pbc,
1110  const Limiter& limiter)
1111 {
1112  const amrex::Box& domain = geom.Domain();
1113  const amrex::Dim3 dlo = amrex::lbound(domain);
1114  const amrex::Dim3 dhi = amrex::ubound(domain);
1115 
1116  const auto dx = geom.CellSizeArray();
1117  AMREX_D_TERM( amrex::Real l_dtdx = dt / dx[0];,
1118  amrex::Real l_dtdy = dt / dx[1];,
1119  amrex::Real l_dtdz = dt / dx[2];);
1120 
1121  amrex::ParallelFor(bx, AMREX_SPACEDIM,
1122  [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
1123  {
1124  PredictVelOnXFace(i,j,k,n,l_dtdx,vel(i,j,k,0),q,Imx,Ipx,pbc[n],dlo.x,dhi.x,limiter);
1125  PredictVelOnYFace(i,j,k,n,l_dtdy,vel(i,j,k,1),q,Imy,Ipy,pbc[n],dlo.y,dhi.y,limiter);
1126 #if (AMREX_SPACEDIM==3)
1127  PredictVelOnZFace(i,j,k,n,l_dtdz,vel(i,j,k,2),q,Imz,Ipz,pbc[n],dlo.z,dhi.z,limiter);
1128 #endif
1129  });
1130 }
1131 
1132 template <typename Limiter>
1135  amrex::Array4<amrex::Real> const& Imy,
1136  amrex::Array4<amrex::Real> const& Imz),
1138  amrex::Array4<amrex::Real> const& Ipy,
1139  amrex::Array4<amrex::Real> const& Ipz),
1142  amrex::Array4<amrex::Real const> const& wmac),
1144  amrex::Geometry geom,
1145  amrex::Real l_dt,
1146  amrex::BCRec const* pbc,
1147  const int ncomp,
1148  const Limiter& limiter,
1149  int limiter_type)
1150 {
1151  const amrex::Box& domain = geom.Domain();
1152  const amrex::Dim3 dlo = amrex::lbound(domain);
1153  const amrex::Dim3 dhi = amrex::ubound(domain);
1154 
1155  AMREX_D_TERM (const auto dx = geom.CellSize(0);,
1156  const auto dy = geom.CellSize(1);,
1157  const auto dz = geom.CellSize(2););
1158 
1159  amrex::ParallelFor(bx, ncomp,
1160  [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
1161  {
1162  PPM::PredictStateOnXFace(i, j, k, n, l_dt, dx, Imx(i,j,k,n), Ipx(i,j,k,n),
1163  q, umac, pbc[n], dlo.x, dhi.x, limiter, limiter_type);
1164  PPM::PredictStateOnYFace(i, j, k, n, l_dt, dy, Imy(i,j,k,n), Ipy(i,j,k,n),
1165  q, vmac, pbc[n], dlo.y, dhi.y, limiter, limiter_type);
1166 #if (AMREX_SPACEDIM==3)
1167  PPM::PredictStateOnZFace(i, j, k, n, l_dt, dz, Imz(i,j,k,n), Ipz(i,j,k,n),
1168  q, wmac, pbc[n], dlo.z, dhi.z, limiter, limiter_type);
1169 #endif
1170  });
1171 }
1172 } // namespace
1173 
1174 #endif
1175 /** @} */
#define AMREX_FORCE_INLINE
#define AMREX_GPU_DEVICE
#define AMREX_GPU_HOST_DEVICE
#define AMREX_D_TERM(a, b, c)
#define AMREX_D_DECL(a, b, c)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE const int * hi() const &noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE const int * lo() const &noexcept
GpuArray< Real, AMREX_SPACEDIM > CellSizeArray() const noexcept
const Real * CellSize() const noexcept
const Box & Domain() const noexcept
static constexpr amrex::Real small_vel
Definition: hydro_constants.H:37
constexpr amrex::Real eps
Definition: hydro_enforce_inout_solvability.cpp:9
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_minmod_bc_lo(const int n, const amrex::Real sm1, const amrex::Real s0, const amrex::Real vel_edge, amrex::Real &dsm, const int bclo, const int domlo)
Definition: hydro_godunov_ppm.H:44
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void minmod_fn(const amrex::Real sm1, const amrex::Real s0, const amrex::Real sp1, amrex::Real &dsm, amrex::Real &dsp)
Definition: hydro_godunov_ppm.H:18
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_minmod_bc_hi(const int n, const amrex::Real s0, const amrex::Real sp1, const amrex::Real vel_edge, amrex::Real &dsp, const int bchi, const int domhi)
Definition: hydro_godunov_ppm.H:64
Definition: hydro_godunov_ppm.H:82
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void PredictVelOnYFace(const int i, const int j, const int k, const int n, const amrex::Real dtdy, const amrex::Real v_ad, const amrex::Array4< const amrex::Real > &S, const amrex::Array4< amrex::Real > &Im, const amrex::Array4< amrex::Real > &Ip, const amrex::BCRec bc, const int domlo, const int domhi, const Limiter &)
Definition: hydro_godunov_ppm.H:758
static constexpr int default_limiter
Definition: hydro_godunov_ppm.H:86
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void PredictStateOnXFace(const int i, const int j, const int k, const int n, const amrex::Real dt, const amrex::Real dx, amrex::Real &Im, amrex::Real &Ip, const amrex::Array4< const amrex::Real > &S, const amrex::Array4< const amrex::Real > &vel_edge, const amrex::BCRec bc, const int domlo, const int domhi, const Limiter &, int limiter_type)
Definition: hydro_godunov_ppm.H:862
void PredictStateOnFaces(amrex::Box const &bx, AMREX_D_DECL(amrex::Array4< amrex::Real > const &Imx, amrex::Array4< amrex::Real > const &Imy, amrex::Array4< amrex::Real > const &Imz), AMREX_D_DECL(amrex::Array4< amrex::Real > const &Ipx, amrex::Array4< amrex::Real > const &Ipy, amrex::Array4< amrex::Real > const &Ipz), AMREX_D_DECL(amrex::Array4< amrex::Real const > const &umac, amrex::Array4< amrex::Real const > const &vmac, amrex::Array4< amrex::Real const > const &wmac), amrex::Array4< amrex::Real const > const &q, amrex::Geometry geom, amrex::Real l_dt, amrex::BCRec const *pbc, const int ncomp, const Limiter &limiter, int limiter_type)
Definition: hydro_godunov_ppm.H:1133
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void SetXBCs(const int i, const int j, const int k, const int n, amrex::Real &sm, amrex::Real &sp, amrex::Real &sedge1, amrex::Real &sedge2, const amrex::Array4< const amrex::Real > &s, const amrex::Real velm, const amrex::Real velp, const int bclo, const int bchi, const int domlo, const int domhi)
Definition: hydro_godunov_ppm.H:418
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void SetYBCs(const int i, const int j, const int k, const int n, amrex::Real &sm, amrex::Real &sp, amrex::Real &sedge1, amrex::Real &sedge2, const amrex::Array4< const amrex::Real > &s, const amrex::Real velm, const amrex::Real velp, const int bclo, const int bchi, const int domlo, const int domhi)
Definition: hydro_godunov_ppm.H:513
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void PredictVelOnXFace(const int i, const int j, const int k, const int n, const amrex::Real dtdx, const amrex::Real v_ad, const amrex::Array4< const amrex::Real > &S, const amrex::Array4< amrex::Real > &Im, const amrex::Array4< amrex::Real > &Ip, const amrex::BCRec bc, const int domlo, const int domhi, const Limiter &)
Definition: hydro_godunov_ppm.H:709
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void PredictStateOnYFace(const int i, const int j, const int k, const int n, const amrex::Real dt, const amrex::Real dx, amrex::Real &Im, amrex::Real &Ip, const amrex::Array4< const amrex::Real > &S, const amrex::Array4< const amrex::Real > &vel_edge, const amrex::BCRec bc, const int domlo, const int domhi, const Limiter &, int limiter_type)
Definition: hydro_godunov_ppm.H:938
void PredictVelOnFaces(amrex::Box const &bx, AMREX_D_DECL(amrex::Array4< amrex::Real > const &Imx, amrex::Array4< amrex::Real > const &Imy, amrex::Array4< amrex::Real > const &Imz), AMREX_D_DECL(amrex::Array4< amrex::Real > const &Ipx, amrex::Array4< amrex::Real > const &Ipy, amrex::Array4< amrex::Real > const &Ipz), amrex::Array4< amrex::Real const > const &q, amrex::Array4< amrex::Real const > const &vel, amrex::Geometry geom, amrex::Real dt, amrex::BCRec const *pbc, const Limiter &limiter)
Definition: hydro_godunov_ppm.H:1098
limiters
Definition: hydro_godunov_ppm.H:84
@ MINMOD
Definition: hydro_godunov_ppm.H:84
@ VanLeer
Definition: hydro_godunov_ppm.H:84
@ WENOZ
Definition: hydro_godunov_ppm.H:84
@ WENO_JS
Definition: hydro_godunov_ppm.H:84
@ NoLimiter
Definition: hydro_godunov_ppm.H:84
@ UPWIND
Definition: hydro_godunov_ppm.H:84
real(kind=amrex_real), parameter half
real(kind=amrex_real), parameter one
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... > ctos, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
constexpr AMREX_GPU_HOST_DEVICE GpuTuple< Args &... > Tie(Args &... args) noexcept
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & min(const T &a, const T &b) noexcept
constexpr AMREX_GPU_HOST_DEVICE GpuTuple< detail::tuple_decay_t< Ts >... > makeTuple(Ts &&... args)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 ubound(Array4< T > const &a) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
Definition: hydro_godunov_ppm.H:380
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::GpuTuple< amrex::Real, amrex::Real > sm_sp(const amrex::Real, const amrex::Real sedge1, const amrex::Real sedge2)
Definition: hydro_godunov_ppm.H:408
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real sedge1(const amrex::Real, const amrex::Real, const amrex::Real, const amrex::Real, const amrex::Real)
Definition: hydro_godunov_ppm.H:386
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real sedge2(const amrex::Real, const amrex::Real, const amrex::Real, const amrex::Real, const amrex::Real)
Definition: hydro_godunov_ppm.H:397
static constexpr bool do_limiting
Definition: hydro_godunov_ppm.H:382
Definition: hydro_godunov_ppm.H:88
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::GpuTuple< amrex::Real, amrex::Real > sm_sp(const amrex::Real, const amrex::Real sedge1, const amrex::Real sedge2)
Definition: hydro_godunov_ppm.H:125
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real sedge2(const amrex::Real, const amrex::Real sm1, const amrex::Real s0, const amrex::Real sp1, const amrex::Real sp2)
Definition: hydro_godunov_ppm.H:110
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real sedge1(const amrex::Real sm2, const amrex::Real sm1, const amrex::Real s0, const amrex::Real sp1, const amrex::Real)
Definition: hydro_godunov_ppm.H:94
static constexpr bool do_limiting
Definition: hydro_godunov_ppm.H:90
Definition: hydro_godunov_ppm.H:344
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real sedge1(const amrex::Real, const amrex::Real, const amrex::Real, const amrex::Real, const amrex::Real)
Definition: hydro_godunov_ppm.H:350
static constexpr bool do_limiting
Definition: hydro_godunov_ppm.H:346
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::GpuTuple< amrex::Real, amrex::Real > sm_sp(const amrex::Real, const amrex::Real sedge1, const amrex::Real sedge2)
Definition: hydro_godunov_ppm.H:372
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real sedge2(const amrex::Real, const amrex::Real, const amrex::Real, const amrex::Real, const amrex::Real)
Definition: hydro_godunov_ppm.H:361
Definition: hydro_godunov_ppm.H:133
static constexpr bool do_limiting
Definition: hydro_godunov_ppm.H:135
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real sedge2(const amrex::Real, const amrex::Real sm1, const amrex::Real s0, const amrex::Real sp1, const amrex::Real sp2)
Definition: hydro_godunov_ppm.H:169
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real vanLeer(const amrex::Real a, const amrex::Real b, const amrex::Real c)
Definition: hydro_godunov_ppm.H:138
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::GpuTuple< amrex::Real, amrex::Real > sm_sp(const amrex::Real s0, const amrex::Real sedge1, const amrex::Real sedge2)
Definition: hydro_godunov_ppm.H:187
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real sedge1(const amrex::Real sm2, const amrex::Real sm1, const amrex::Real s0, const amrex::Real sp1, const amrex::Real)
Definition: hydro_godunov_ppm.H:152
Definition: hydro_godunov_ppm.H:281
static constexpr bool do_limiting
Definition: hydro_godunov_ppm.H:283
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real sedge1(const amrex::Real sm2, const amrex::Real sm1, const amrex::Real s0, const amrex::Real sp1, const amrex::Real sp2)
Definition: hydro_godunov_ppm.H:325
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::GpuTuple< amrex::Real, amrex::Real > sm_sp(const amrex::Real, const amrex::Real sedge1, const amrex::Real sedge2)
Definition: hydro_godunov_ppm.H:336
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real sedge2(const amrex::Real sm2, const amrex::Real sm1, const amrex::Real s0, const amrex::Real sp1, const amrex::Real sp2)
Definition: hydro_godunov_ppm.H:287
Definition: hydro_godunov_ppm.H:222
static constexpr bool do_limiting
Definition: hydro_godunov_ppm.H:224
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::GpuTuple< amrex::Real, amrex::Real > sm_sp(const amrex::Real, const amrex::Real sedge1, const amrex::Real sedge2)
Definition: hydro_godunov_ppm.H:273
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real sedge2(const amrex::Real sm2, const amrex::Real sm1, const amrex::Real s0, const amrex::Real sp1, const amrex::Real sp2)
Definition: hydro_godunov_ppm.H:228
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real sedge1(const amrex::Real sm2, const amrex::Real sm1, const amrex::Real s0, const amrex::Real sp1, const amrex::Real sp2)
Definition: hydro_godunov_ppm.H:262