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  const amrex::Real sm2 = S(i-2,j,k,n);
722  const amrex::Real sm1 = S(i-1,j,k,n);
723  const amrex::Real s0 = S(i ,j,k,n);
724  const amrex::Real sp1 = S(i+1,j,k,n);
725  const 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  auto [sm, sp] = Limiter::sm_sp(s0, sedge1, sedge2);
731 
732  SetXBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, v_ad, v_ad,
733  bc.lo(0), bc.hi(0), domlo, domhi);
734 
735  const amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp);
736 
737  const amrex::Real sigma = amrex::Math::abs(v_ad)*dtdx;
738 
739  if (v_ad > small_vel)
740  {
741  Ip(i,j,k,n) = sp - (half*sigma)*((sp-sm) - (one - two3rds*sigma)*s6);
742  Im(i,j,k,n) = S(i,j,k,n);
743  }
744  else if (v_ad < -small_vel)
745  {
746  Ip(i,j,k,n) = S(i,j,k,n);
747  Im(i,j,k,n) = sm + (half*sigma)*((sp-sm) + (one - two3rds*sigma)*s6);
748  } else
749  {
750  Ip(i,j,k,n) = S(i,j,k,n);
751  Im(i,j,k,n) = S(i,j,k,n);
752  }
753 }
754 
755 template <typename Limiter>
757 void PredictVelOnYFace ( const int i, const int j, const int k, const int n,
758  const amrex::Real dtdy, const amrex::Real v_ad,
760  const amrex::Array4<amrex::Real> &Im,
761  const amrex::Array4<amrex::Real> &Ip,
762  const amrex::BCRec bc, const int domlo, const int domhi,
763  const Limiter& /*limiter*/)
764 {
765  constexpr auto half{amrex::Real(0.5)};
766  constexpr auto one{amrex::Real(1.0)};
767  constexpr auto two3rds{amrex::Real(2.0/3.0)};
768 
769  const amrex::Real sm2 = S(i,j-2,k,n);
770  const amrex::Real sm1 = S(i,j-1,k,n);
771  const amrex::Real s0 = S(i,j ,k,n);
772  const amrex::Real sp1 = S(i,j+1,k,n);
773  const amrex::Real sp2 = S(i,j+2,k,n);
774 
775  amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
776  amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
777 
778  auto [sm, sp] = Limiter::sm_sp(s0, sedge1, sedge2);
779 
780  SetYBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, v_ad, v_ad,
781  bc.lo(1), bc.hi(1), domlo, domhi);
782 
783  const amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp);
784 
785  const amrex::Real sigma = amrex::Math::abs(v_ad)*dtdy;
786 
787  if (v_ad > small_vel)
788  {
789  Ip(i,j,k,n) = sp - (half*sigma)*((sp-sm) - (one - two3rds*sigma)*s6);
790  Im(i,j,k,n) = S(i,j,k,n);
791  }
792  else if (v_ad < -small_vel)
793  {
794  Ip(i,j,k,n) = S(i,j,k,n);
795  Im(i,j,k,n) = sm + (half*sigma)*((sp-sm) + (one - two3rds*sigma)*s6);
796  } else
797  {
798  Ip(i,j,k,n) = S(i,j,k,n);
799  Im(i,j,k,n) = S(i,j,k,n);
800  }
801 }
802 
803 #if (AMREX_SPACEDIM==3)
804 template <typename Limiter>
806 void PredictVelOnZFace ( const int i, const int j, const int k, const int n,
807  const amrex::Real dtdz, const amrex::Real v_ad,
809  const amrex::Array4<amrex::Real> &Im,
810  const amrex::Array4<amrex::Real> &Ip,
811  const amrex::BCRec bc, const int domlo, const int domhi,
812  const Limiter& /*limiter*/)
813 {
814  constexpr auto half{amrex::Real(0.5)};
815  constexpr auto one{amrex::Real(1.0)};
816  constexpr auto two3rds{amrex::Real(2.0/3.0)};
817 
818  const amrex::Real sm2 = S(i,j,k-2,n);
819  const amrex::Real sm1 = S(i,j,k-1,n);
820  const amrex::Real s0 = S(i,j,k ,n);
821  const amrex::Real sp1 = S(i,j,k+1,n);
822  const amrex::Real sp2 = S(i,j,k+2,n);
823 
824  amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
825  amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
826 
827  auto [sm, sp] = Limiter::sm_sp(s0, sedge1, sedge2);
828 
829  SetZBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, v_ad, v_ad,
830  bc.lo(2), bc.hi(2), domlo, domhi);
831 
832  const amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp);
833 
834  const amrex::Real sigma = amrex::Math::abs(v_ad)*dtdz;
835 
836  if (v_ad > small_vel)
837  {
838  Ip(i,j,k,n) = sp - (half*sigma)*((sp-sm) - (one - two3rds*sigma)*s6);
839  Im(i,j,k,n) = S(i,j,k,n);
840  }
841  else if (v_ad < -small_vel)
842  {
843  Ip(i,j,k,n) = S(i,j,k,n);
844  Im(i,j,k,n) = sm + (half*sigma)*((sp-sm) + (one - two3rds*sigma)*s6);
845  } else
846  {
847  Ip(i,j,k,n) = S(i,j,k,n);
848  Im(i,j,k,n) = S(i,j,k,n);
849  }
850 }
851 
852 #endif
853 
854 // Right now only ppm type 1 is supported on GPU
855 // This version is called after the MAC projection, when we use the MAC-projected velocity
856 // for upwinding
857 template <typename Limiter>
859 void PredictStateOnXFace ( const int i, const int j, const int k, const int n,
860  const amrex::Real dt, const amrex::Real dx,
861  amrex::Real& Im, amrex::Real& Ip,
863  const amrex::Array4<const amrex::Real> &vel_edge,
864  const amrex::BCRec bc,
865  const int domlo, const int domhi,
866  const Limiter& /*limiter*/,
867  int /*unused*/)
868 {
869  const amrex::Real sm1 = S(i-1,j,k,n);
870  const amrex::Real s0 = S(i ,j,k,n);
871  const amrex::Real sp1 = S(i+1,j,k,n);
872 
873  constexpr auto half{amrex::Real(0.5)};
874  constexpr auto one{amrex::Real(1.0)};
875 
876  const amrex::Real sigmap = amrex::Math::abs(vel_edge(i+1,j,k))*dt/dx;
877  const amrex::Real sigmam = amrex::Math::abs(vel_edge(i ,j,k))*dt/dx;
878 
879  constexpr auto two3rds{amrex::Real(2.0/3.0)};
880 
881  const amrex::Real sm2 = S(i-2,j,k,n);
882  const amrex::Real sp2 = S(i+2,j,k,n);
883 
884  amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
885  amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
886 
887  auto [sm, sp] = Limiter::sm_sp(s0, sedge1, sedge2);
888 
889  SetXBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k),
890  bc.lo(0), bc.hi(0), domlo, domhi);
891 
892  const amrex::Real s6 = 6.0*s0 - 3.0*(sm + sp);
893 
894  if (vel_edge(i+1,j,k) > small_vel) {
895  Ip = sp - (half*sigmap)*((sp - sm) - (one -two3rds*sigmap)*s6);
896  } else {
897  Ip = S(i,j,k,n);
898  }
899 
900  if(vel_edge(i,j,k) < -small_vel) {
901  Im = sm + (half*sigmam)*((sp-sm) + (one - two3rds*sigmam)*s6);
902  } else {
903  Im = S(i,j,k,n);
904  }
905 }
906 
907 template <>
909 void PredictStateOnXFace<PPM::upwind> ( const int i, const int j, const int k, const int n,
910  const amrex::Real /*unused*/, const amrex::Real /*unused*/,
911  amrex::Real& Im, amrex::Real& Ip,
913  const amrex::Array4<const amrex::Real> & /*unused*/,
914  const amrex::BCRec /*unused*/,
915  const int /*unused*/, const int /*unused*/,
916  const PPM::upwind& /*limiter*/,
917  int /*unused*/)
918 {
919  Im = S(i,j,k,n);
920  Ip = S(i,j,k,n);
921 }
922 
923 
924 template <>
926 void PredictStateOnXFace<PPM::minmod> ( const int i, const int j, const int k, const int n,
927  const amrex::Real dt, const amrex::Real dx,
928  amrex::Real& Im, amrex::Real& Ip,
930  const amrex::Array4<const amrex::Real> &vel_edge,
931  const amrex::BCRec bc,
932  const int domlo, const int domhi,
933  const PPM::minmod& /*limiter*/,
934  int /*unused*/)
935 {
936  const amrex::Real sm1 = S(i-1,j,k,n);
937  const amrex::Real s0 = S(i ,j,k,n);
938  const amrex::Real sp1 = S(i+1,j,k,n);
939 
940  constexpr auto half{amrex::Real(0.5)};
941  constexpr auto one{amrex::Real(1.0)};
942 
943  const amrex::Real sigmap = amrex::Math::abs(vel_edge(i+1,j,k))*dt/dx;
944  const amrex::Real sigmam = amrex::Math::abs(vel_edge(i ,j,k))*dt/dx;
945 
946  amrex::Real dsp = 0.0;
947  amrex::Real dsm = 0.0;
948  minmod_fn(sm1, s0, sp1, dsm, dsp);
949  Godunov_minmod_bc_lo(i, sm1, s0, vel_edge(i ,j,k), dsm, bc.lo(0), domlo);
950  Godunov_minmod_bc_hi(i, s0, sp1, vel_edge(i+1,j,k), dsp, bc.hi(0), domhi);
951 
952  if (vel_edge(i + 1, j, k) > small_vel) {
953  Ip = s0 + half * (one - sigmap) * dsp;
954  } else {
955  Ip = S(i,j,k,n);
956  }
957 
958  if (vel_edge(i, j, k) < -small_vel) {
959  Im = s0 - half * (one - sigmam) * dsm;
960  } else {
961  Im = S(i,j,k,n);
962  }
963 }
964 
965 template <typename Limiter>
967 void PredictStateOnYFace ( const int i, const int j, const int k, const int n,
968  const amrex::Real dt, const amrex::Real dx,
969  amrex::Real& Im, amrex::Real& Ip,
971  const amrex::Array4<const amrex::Real> &vel_edge,
972  const amrex::BCRec bc,
973  const int domlo, const int domhi,
974  const Limiter& /*limiter*/,
975  int /*unused*/)
976 {
977  const amrex::Real sm1 = S(i,j-1,k,n);
978  const amrex::Real s0 = S(i,j ,k,n);
979  const amrex::Real sp1 = S(i,j+1,k,n);
980 
981  constexpr auto half{amrex::Real(0.5)};
982  constexpr auto one{amrex::Real(1.0)};
983 
984  const amrex::Real sigmap = amrex::Math::abs(vel_edge(i,j+1,k))*dt/dx;
985  const amrex::Real sigmam = amrex::Math::abs(vel_edge(i,j ,k))*dt/dx;
986 
987  constexpr auto two3rds{amrex::Real(2.0/3.0)};
988 
989  const amrex::Real sm2 = S(i,j-2,k,n);
990  const amrex::Real sp2 = S(i,j+2,k,n);
991 
992  amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
993  amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
994 
995  auto [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  const amrex::Real 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 }
1014 
1015 
1016 template <>
1018 void PredictStateOnYFace<PPM::upwind> ( const int i, const int j, const int k, const int n,
1019  const amrex::Real /*unused*/, const amrex::Real /*unused*/,
1020  amrex::Real& Im, amrex::Real& Ip,
1022  const amrex::Array4<const amrex::Real> & /*unused*/,
1023  const amrex::BCRec /*unused*/,
1024  const int /*unused*/, const int /*unused*/,
1025  const PPM::upwind& /*limiter*/,
1026  int /*unused*/)
1027 {
1028  Im = S(i,j,k,n);
1029  Ip = S(i,j,k,n);
1030 }
1031 
1032 
1033 template <>
1035 void PredictStateOnYFace<PPM::minmod> ( const int i, const int j, const int k, const int n,
1036  const amrex::Real dt, const amrex::Real dx,
1037  amrex::Real& Im, amrex::Real& Ip,
1039  const amrex::Array4<const amrex::Real> &vel_edge,
1040  const amrex::BCRec bc,
1041  const int domlo, const int domhi,
1042  const PPM::minmod& /*limiter*/,
1043  int /*unused*/)
1044 {
1045  const amrex::Real sm1 = S(i,j-1,k,n);
1046  const amrex::Real s0 = S(i,j ,k,n);
1047  const amrex::Real sp1 = S(i,j+1,k,n);
1048 
1049  constexpr auto half{amrex::Real(0.5)};
1050  constexpr auto one{amrex::Real(1.0)};
1051 
1052  const amrex::Real sigmap = amrex::Math::abs(vel_edge(i,j+1,k))*dt/dx;
1053  const amrex::Real sigmam = amrex::Math::abs(vel_edge(i,j ,k))*dt/dx;
1054 
1055  amrex::Real dsp = 0.0;
1056  amrex::Real dsm = 0.0;
1057  minmod_fn(sm1, s0, sp1, dsm, dsp);
1058  Godunov_minmod_bc_lo(j, sm1, s0, vel_edge(i,j ,k), dsm, bc.lo(1), domlo);
1059  Godunov_minmod_bc_hi(j, s0, sp1, vel_edge(i,j+1,k), dsp, bc.hi(1), domhi);
1060 
1061  if (vel_edge(i,j+1,k) > small_vel) {
1062  Ip = s0 + half * (one - sigmap) * dsp;
1063  } else {
1064  Ip = S(i,j,k,n);
1065  }
1066  if (vel_edge(i,j,k) < -small_vel) {
1067  Im = s0 - half * (one - sigmam) * dsm;
1068  } else {
1069  Im = S(i,j,k,n);
1070  }
1071 }
1072 
1073 
1074 #if (AMREX_SPACEDIM==3)
1075 template <typename Limiter>
1077 void PredictStateOnZFace ( const int i, const int j, const int k, const int n,
1078  const amrex::Real dt, const amrex::Real dx,
1079  amrex::Real& Im, amrex::Real& Ip,
1081  const amrex::Array4<const amrex::Real> &vel_edge,
1082  const amrex::BCRec bc,
1083  const int domlo, const int domhi,
1084  const Limiter& /*limiter*/,
1085  int /*unused*/)
1086 {
1087  const amrex::Real sm1 = S(i,j,k-1,n);
1088  const amrex::Real s0 = S(i,j,k ,n);
1089  const amrex::Real sp1 = S(i,j,k+1,n);
1090 
1091  constexpr auto half{amrex::Real(0.5)};
1092 
1093  const amrex::Real sigmap = amrex::Math::abs(vel_edge(i,j,k+1))*dt/dx;
1094  const amrex::Real sigmam = amrex::Math::abs(vel_edge(i,j,k ))*dt/dx;
1095 
1096  constexpr auto two3rds{amrex::Real(2.0/3.0)};
1097 
1098  const amrex::Real sm2 = S(i,j,k-2,n);
1099  const amrex::Real sp2 = S(i,j,k+2,n);
1100 
1101  amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
1102  amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
1103 
1104  auto [sm, sp] = Limiter::sm_sp(s0, sedge1, sedge2);
1105 
1106  SetZBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k),
1107  bc.lo(2), bc.hi(2), domlo, domhi);
1108 
1109  const amrex::Real s6 = 6.0*s0- 3.0*(sm + sp);
1110 
1111  if(vel_edge(i,j,k+1) > small_vel) {
1112  Ip = sp - (half*sigmap)*((sp-sm) - (amrex::Real(1.0) -two3rds*sigmap)*s6);
1113  } else {
1114  Ip = S(i,j,k,n);
1115  }
1116 
1117  if(vel_edge(i,j,k) < -small_vel) {
1118  Im = sm + (half*sigmam)*((sp-sm) + (amrex::Real(1.0) - two3rds*sigmam)*s6);
1119  } else {
1120  Im = S(i,j,k,n);
1121  }
1122 }
1123 
1124 template <>
1126 void PredictStateOnZFace<PPM::upwind> ( const int i, const int j, const int k, const int n,
1127  const amrex::Real /*unused*/, const amrex::Real /*unused*/,
1128  amrex::Real& Im, amrex::Real& Ip,
1130  const amrex::Array4<const amrex::Real> & /*unused*/,
1131  const amrex::BCRec /*unused*/,
1132  const int /*unused*/, const int /*unused*/,
1133  const PPM::upwind& /*limiter*/,
1134  int /*unused*/)
1135 {
1136  Im = S(i,j,k,n);
1137  Ip = S(i,j,k,n);
1138 }
1139 
1140 template <>
1142 void PredictStateOnZFace<PPM::minmod> ( const int i, const int j, const int k, const int n,
1143  const amrex::Real dt, const amrex::Real dx,
1144  amrex::Real& Im, amrex::Real& Ip,
1146  const amrex::Array4<const amrex::Real> &vel_edge,
1147  const amrex::BCRec bc,
1148  const int domlo, const int domhi,
1149  const PPM::minmod& /*limiter*/,
1150  int /*unused*/)
1151 {
1152  const amrex::Real sm1 = S(i,j,k-1,n);
1153  const amrex::Real s0 = S(i,j,k ,n);
1154  const amrex::Real sp1 = S(i,j,k+1,n);
1155 
1156  constexpr auto half{amrex::Real(0.5)};
1157  constexpr auto one{amrex::Real(1.0)};
1158 
1159  const amrex::Real sigmap = amrex::Math::abs(vel_edge(i,j,k+1))*dt/dx;
1160  const amrex::Real sigmam = amrex::Math::abs(vel_edge(i,j,k ))*dt/dx;
1161 
1162  amrex::Real dsp = 0.0;
1163  amrex::Real dsm = 0.0;
1164  minmod_fn(sm1, s0, sp1, dsm, dsp);
1165  Godunov_minmod_bc_lo(k, sm1, s0, vel_edge(i,j,k ), dsm, bc.lo(2), domlo);
1166  Godunov_minmod_bc_hi(k, s0, sp1, vel_edge(i,j,k+1), dsp, bc.hi(2), domhi);
1167 
1168  if (vel_edge(i, j, k + 1) > small_vel) {
1169  Ip = s0 + half * (one - sigmap) * dsp;
1170  } else {
1171  Ip = S(i,j,k,n);
1172  }
1173 
1174  if (vel_edge(i, j, k) < -small_vel) {
1175  Im = s0 - half * (one - sigmam) * dsm;
1176  } else {
1177  Im = S(i,j,k,n);
1178  }
1179 }
1180 #endif
1181 
1182 template <typename Limiter>
1185  amrex::Array4<amrex::Real> const& Imy,
1186  amrex::Array4<amrex::Real> const& Imz),
1188  amrex::Array4<amrex::Real> const& Ipy,
1189  amrex::Array4<amrex::Real> const& Ipz),
1192  amrex::Geometry geom,
1193  amrex::Real dt,
1194  amrex::BCRec const* pbc,
1195  const Limiter& limiter)
1196 {
1197  const amrex::Box& domain = geom.Domain();
1198  const amrex::Dim3 dlo = amrex::lbound(domain);
1199  const amrex::Dim3 dhi = amrex::ubound(domain);
1200 
1201  const auto dx = geom.CellSizeArray();
1202  AMREX_D_TERM( const amrex::Real l_dtdx = dt / dx[0];,
1203  const amrex::Real l_dtdy = dt / dx[1];,
1204  const amrex::Real l_dtdz = dt / dx[2];);
1205 
1206  amrex::ParallelFor(bx, AMREX_SPACEDIM,
1207  [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
1208  {
1209  PredictVelOnXFace(i,j,k,n,l_dtdx,vel(i,j,k,0),q,Imx,Ipx,pbc[n],dlo.x,dhi.x,limiter);
1210  PredictVelOnYFace(i,j,k,n,l_dtdy,vel(i,j,k,1),q,Imy,Ipy,pbc[n],dlo.y,dhi.y,limiter);
1211 #if (AMREX_SPACEDIM==3)
1212  PredictVelOnZFace(i,j,k,n,l_dtdz,vel(i,j,k,2),q,Imz,Ipz,pbc[n],dlo.z,dhi.z,limiter);
1213 #endif
1214  });
1215 }
1216 
1217 template <typename Limiter>
1220  amrex::Array4<amrex::Real> const& Imy,
1221  amrex::Array4<amrex::Real> const& Imz),
1223  amrex::Array4<amrex::Real> const& Ipy,
1224  amrex::Array4<amrex::Real> const& Ipz),
1227  amrex::Array4<amrex::Real const> const& wmac),
1229  amrex::Geometry geom,
1230  amrex::Real l_dt,
1231  amrex::BCRec const* pbc,
1232  const int ncomp,
1233  const Limiter& limiter,
1234  int limiter_type)
1235 {
1236  const amrex::Box& domain = geom.Domain();
1237  const amrex::Dim3 dlo = amrex::lbound(domain);
1238  const amrex::Dim3 dhi = amrex::ubound(domain);
1239 
1240  AMREX_D_TERM (const auto dx = geom.CellSize(0);,
1241  const auto dy = geom.CellSize(1);,
1242  const auto dz = geom.CellSize(2););
1243 
1244  amrex::ParallelFor(bx, ncomp,
1245  [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
1246  {
1247  PPM::PredictStateOnXFace(i, j, k, n, l_dt, dx, Imx(i,j,k,n), Ipx(i,j,k,n),
1248  q, umac, pbc[n], dlo.x, dhi.x, limiter, limiter_type);
1249  PPM::PredictStateOnYFace(i, j, k, n, l_dt, dy, Imy(i,j,k,n), Ipy(i,j,k,n),
1250  q, vmac, pbc[n], dlo.y, dhi.y, limiter, limiter_type);
1251 #if (AMREX_SPACEDIM==3)
1252  PPM::PredictStateOnZFace(i, j, k, n, l_dt, dz, Imz(i,j,k,n), Ipz(i,j,k,n),
1253  q, wmac, pbc[n], dlo.z, dhi.z, limiter, limiter_type);
1254 #endif
1255  });
1256 }
1257 } // namespace
1258 
1259 #endif
1260 /** @} */
#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
const Real * CellSize() const noexcept
GpuArray< Real, AMREX_SPACEDIM > CellSizeArray() const noexcept
const Box & Domain() const noexcept
constexpr amrex::Real eps
Definition: hydro_bds_edge_state_2D.cpp:13
static constexpr amrex::Real small_vel
Definition: hydro_constants.H:37
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:757
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)
Definition: hydro_godunov_ppm.H:859
static constexpr int default_limiter
Definition: hydro_godunov_ppm.H:86
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:1218
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 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)
Definition: hydro_godunov_ppm.H:967
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
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:1183
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 AMREX_FORCE_INLINE constexpr const T & min(const T &a, const T &b) noexcept
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
AMREX_GPU_HOST_DEVICE constexpr GpuTuple< detail::tuple_decay_t< Ts >... > makeTuple(Ts &&... args)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr const T & max(const T &a, const T &b) 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