8 #ifndef HYDRO_GODUNOV_PPM_H
9 #define HYDRO_GODUNOV_PPM_H
20 const amrex::Real sp1,
28 if (!(dsp * dsm < 0.0)) {
30 if (std::abs(dsp) < std::abs(dsm)) {
45 const amrex::Real sm1,
47 const amrex::Real vel_edge,
54 if ( bclo == amrex::BCType::ext_dir || bclo == amrex::BCType::hoextrap ||
55 (bclo == amrex::BCType::direction_dependent && vel_edge >= 0.0) )
66 const amrex::Real sp1,
67 const amrex::Real vel_edge,
74 if ( bchi == amrex::BCType::ext_dir || bchi == amrex::BCType::hoextrap ||
75 (bchi == amrex::BCType::direction_dependent && vel_edge <= 0.0) )
95 const amrex::Real sm1,
97 const amrex::Real sp1,
101 amrex::Real d1 = amrex::Real(0.5) * (sp1 - sm1);
102 amrex::Real d2 = amrex::Real(0.5) * (s0 - sm2);
104 amrex::Real sedge = amrex::Real(0.5)*(s0 + sm1) - amrex::Real(1./6.)*(d1 - d2);
111 const amrex::Real sm1,
112 const amrex::Real s0,
113 const amrex::Real sp1,
114 const amrex::Real sp2)
116 amrex::Real d1 = amrex::Real(0.5) * (sp2 - s0);
117 amrex::Real d2 = amrex::Real(0.5) * (sp1 - sm1);
119 amrex::Real sedge = amrex::Real(0.5)*(sp1 + s0) - amrex::Real(1./6.)*(d1 - d2);
138 static amrex::Real
vanLeer(
const amrex::Real a,
const amrex::Real b,
const amrex::Real c)
140 constexpr
auto small_qty_sq = amrex::Real(1.e-10) * amrex::Real(1.e-10);
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);
153 const amrex::Real sm1,
154 const amrex::Real s0,
155 const amrex::Real sp1,
159 amrex::Real d1 =
vanLeer(s0,sp1,sm1);
160 amrex::Real d2 =
vanLeer(sm1,s0,sm2);
162 amrex::Real sedge = amrex::Real(0.5)*(s0 + sm1) - amrex::Real(1./6.)*(d1 - d2);
170 const amrex::Real sm1,
171 const amrex::Real s0,
172 const amrex::Real sp1,
173 const amrex::Real sp2)
176 amrex::Real d1 =
vanLeer(sp1,sp2,s0);
177 amrex::Real d2 =
vanLeer(s0,sp1,sm1);
179 amrex::Real sedge = amrex::Real(0.5)*(sp1 + s0) - amrex::Real(1./6.)*(d1 - d2);
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;
229 const amrex::Real sm1,
230 const amrex::Real s0,
231 const amrex::Real sp1,
232 const amrex::Real sp2)
234 constexpr
auto eps = amrex::Real(1.0e-6);
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);
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));
251 const amrex::Real omega = omega1 + omega2 + omega3;
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;
257 return (omega1 * v_1 + omega2 * v_2 + omega3 * v_3) / (amrex::Real(6.0) * omega);
263 const amrex::Real sm1,
264 const amrex::Real s0,
265 const amrex::Real sp1,
266 const amrex::Real sp2)
268 return sedge2(sp2,sp1,s0,sm1,sm2);
288 const amrex::Real sm1,
289 const amrex::Real s0,
290 const amrex::Real sp1,
291 const amrex::Real sp2)
293 constexpr
auto eps = amrex::Real(1.0e-6);
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);
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);
314 const amrex::Real omega = omega1 + omega2 + omega3;
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;
320 return (omega1 * v_1 + omega2 * v_2 + omega3 * v_3) / (amrex::Real(6.0) * omega);
326 const amrex::Real sm1,
327 const amrex::Real s0,
328 const amrex::Real sp1,
329 const amrex::Real sp2)
331 return sedge2(sp2,sp1,s0,sm1,sm2);
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)
426 using namespace amrex;
428 if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) ||
429 (bclo == BCType::direction_dependent && velp >= 0.0) )
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) {
440 sm = s(domlo-1,j,k,n);
443 }
else if (i == domlo+1) {
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) {
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))
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;
469 if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) ||
470 (bchi == BCType::direction_dependent && velm <= 0.0) )
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) {
481 sp = s(domhi+1,j,k,n);
484 }
else if (i == domhi-1) {
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) {
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))
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;
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)
521 using namespace amrex;
523 if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) ||
524 (bclo == BCType::direction_dependent && velm >= 0.0) )
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) {
535 sm = s(i,domlo-1,k,n);
538 }
else if (j == domlo+1) {
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) {
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))
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;
564 if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) ||
565 (bchi == BCType::direction_dependent && velp <= 0.0) )
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) {
576 sp = s(i,domhi+1,k, n);
579 }
else if (j == domhi-1) {
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) {
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;
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)
616 using namespace amrex;
618 if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) ||
619 (bclo == BCType::direction_dependent && velm >= 0.0) )
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) {
630 sm = s(i,j,domlo-1,n);
633 }
else if (k == domlo+1) {
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) {
645 if constexpr (Limiter::do_limiting) {
646 if ( (sp - s(i,j,domlo+1,n))*(s(i,j,domlo+1,n) - sm) <= 0. )
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;
659 if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) ||
660 (bchi == BCType::direction_dependent && velp <= 0.0) )
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) {
671 sp = s(i,j,domhi+1,n);
674 }
else if (k == domhi-1) {
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) {
686 if constexpr (Limiter::do_limiting) {
687 if ( (sp - s(i,j,domhi-1,n))*(s(i,j,domhi-1,n) - sm) <= 0. )
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;
707 template <
typename Limiter>
710 const amrex::Real dtdx,
const amrex::Real v_ad,
714 const amrex::BCRec bc,
const int domlo,
const int domhi,
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)};
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);
727 amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
728 amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
731 amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
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);
736 amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp);
738 amrex::Real sigma = amrex::Math::abs(v_ad)*dtdx;
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);
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);
751 Ip(i,j,k,n) = S(i,j,k,n);
752 Im(i,j,k,n) = S(i,j,k,n);
756 template <
typename Limiter>
759 const amrex::Real dtdy,
const amrex::Real v_ad,
763 const amrex::BCRec bc,
const int domlo,
const int domhi,
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)};
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);
776 amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
777 amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
780 amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
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);
785 amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp);
787 amrex::Real sigma = amrex::Math::abs(v_ad)*dtdy;
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);
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);
800 Ip(i,j,k,n) = S(i,j,k,n);
801 Im(i,j,k,n) = S(i,j,k,n);
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,
813 const amrex::BCRec bc,
const int domlo,
const int domhi,
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)};
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);
826 amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
827 amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
830 amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
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);
835 amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp);
837 amrex::Real sigma = amrex::Math::abs(v_ad)*dtdz;
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);
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);
850 Ip(i,j,k,n) = S(i,j,k,n);
851 Im(i,j,k,n) = S(i,j,k,n);
860 template <
typename Limiter>
863 const amrex::Real dt,
const amrex::Real dx,
864 amrex::Real& Im, amrex::Real& Ip,
868 const int domlo,
const int domhi,
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);
876 constexpr
auto half{amrex::Real(0.5)};
877 constexpr
auto one{amrex::Real(1.0)};
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;
887 amrex::Real dsp = 0.0;
888 amrex::Real dsm = 0.0;
894 Ip = s0 +
half * (
one - sigmap) * dsp;
900 Im = s0 -
half * (
one - sigmam) * dsm;
906 constexpr
auto two3rds{amrex::Real(2.0/3.0)};
908 amrex::Real sm2 = S(i-2,j,k,n);
909 amrex::Real sp2 = S(i+2,j,k,n);
911 amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
912 amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
915 amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
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);
920 amrex::Real s6 = 6.0*s0 - 3.0*(sm + sp);
923 Ip = sp - (
half*sigmap)*((sp - sm) - (
one -two3rds*sigmap)*s6);
929 Im = sm + (
half*sigmam)*((sp-sm) + (
one - two3rds*sigmam)*s6);
936 template <
typename Limiter>
939 const amrex::Real dt,
const amrex::Real dx,
940 amrex::Real& Im, amrex::Real& Ip,
944 const int domlo,
const int domhi,
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);
952 constexpr
auto half{amrex::Real(0.5)};
953 constexpr
auto one{amrex::Real(1.0)};
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;
967 amrex::Real dsp = 0.0;
968 amrex::Real dsm = 0.0;
974 Ip = s0 +
half * (
one - sigmap) * dsp;
979 Im = s0 -
half * (
one - sigmam) * dsm;
986 constexpr
auto two3rds{amrex::Real(2.0/3.0)};
988 amrex::Real sm2 = S(i,j-2,k,n);
989 amrex::Real sp2 = S(i,j+2,k,n);
991 amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
992 amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
995 amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
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);
1000 s6 = 6.0*s0- 3.0*(sm + sp);
1003 Ip = sp - (
half*sigmap)*((sp - sm) - (
one -two3rds*sigmap)*s6);
1009 Im = sm + (
half*sigmam)*((sp-sm) + (
one - two3rds*sigmam)*s6);
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,
1026 const int domlo,
const int domhi,
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);
1034 constexpr
auto half{amrex::Real(0.5)};
1035 constexpr
auto one{amrex::Real(1.0)};
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;
1047 amrex::Real dsp = 0.0;
1048 amrex::Real dsm = 0.0;
1053 if (vel_edge(i, j, k + 1) >
small_vel) {
1054 Ip = s0 +
half * (
one - sigmap) * dsp;
1060 Im = s0 -
half * (
one - sigmam) * dsm;
1066 constexpr
auto two3rds{amrex::Real(2.0/3.0)};
1068 amrex::Real sm2 = S(i,j,k-2,n);
1069 amrex::Real sp2 = S(i,j,k+2,n);
1071 amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2);
1072 amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2);
1075 amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
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);
1080 amrex::Real s6 = 6.0*s0- 3.0*(sm + sp);
1083 Ip = sp - (
half*sigmap)*((sp-sm) - (amrex::Real(1.0) -two3rds*sigmap)*s6);
1089 Im = sm + (
half*sigmam)*((sp-sm) + (amrex::Real(1.0) - two3rds*sigmam)*s6);
1097 template <
typename Limiter>
1110 const Limiter& limiter)
1118 amrex::Real l_dtdy = dt / dx[1];,
1119 amrex::Real l_dtdz = dt / dx[2];);
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);
1132 template <
typename Limiter>
1148 const Limiter& limiter,
1157 const auto dz = geom.
CellSize(2););
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);
#define AMREX_FORCE_INLINE
#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