AMReX-Hydro
AMReX-based hydro routines for low Mach number flows
 
Loading...
Searching...
No Matches
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
17void
18minmod_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
43void
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
63void
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
82namespace PPM {
83
85
86static constexpr int default_limiter = VanLeer;
87
88struct 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 {
130 }
131}; // nolimiter
132
133struct 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
222struct 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 {
278 }
279}; // wenoz
280
281struct 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 {
341 }
342}; // weno_js
343
344struct 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 {
377 }
378}; // upwind
379
380struct 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 {
413 }
414}; // minmod
415
416template <typename Limiter>
418void 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
511template <typename Limiter>
513void 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)
606template <typename Limiter>
608void 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
707template <typename Limiter>
709void PredictVelOnXFace ( const int i, const int j, const int k, const int n,
710 const amrex::Real dtdx, const amrex::Real v_ad,
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
755template <typename Limiter>
757void PredictVelOnYFace ( const int i, const int j, const int k, const int n,
758 const amrex::Real dtdy, const amrex::Real v_ad,
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)
804template <typename Limiter>
806void PredictVelOnZFace ( const int i, const int j, const int k, const int n,
807 const amrex::Real dtdz, const amrex::Real v_ad,
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
857template <typename Limiter>
859void 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
907template <>
909void 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
924template <>
926void 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
965template <typename Limiter>
967void 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
1016template <>
1018void 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
1033template <>
1035void 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)
1075template <typename Limiter>
1077void 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
1124template <>
1126void 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
1140template <>
1142void 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
1182template <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
1217template <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),
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)
__host__ __device__ const int * hi() const &noexcept
__host__ __device__ const int * lo() const &noexcept
const Real * CellSize() const noexcept
GpuArray< Real, 3 > 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 PredictStateOnXFace< PPM::minmod >(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 PPM::minmod &, int)
Definition hydro_godunov_ppm.H:926
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 PredictStateOnXFace< PPM::upwind >(const int i, const int j, const int k, const int n, const amrex::Real, const amrex::Real, amrex::Real &Im, amrex::Real &Ip, const amrex::Array4< const amrex::Real > &S, const amrex::Array4< const amrex::Real > &, const amrex::BCRec, const int, const int, const PPM::upwind &, int)
Definition hydro_godunov_ppm.H:909
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 PredictStateOnYFace< PPM::minmod >(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 PPM::minmod &, int)
Definition hydro_godunov_ppm.H:1035
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< PPM::upwind >(const int i, const int j, const int k, const int n, const amrex::Real, const amrex::Real, amrex::Real &Im, amrex::Real &Ip, const amrex::Array4< const amrex::Real > &S, const amrex::Array4< const amrex::Real > &, const amrex::BCRec, const int, const int, const PPM::upwind &, int)
Definition hydro_godunov_ppm.H:1018
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
__host__ __device__ Dim3 ubound(Array4< T > const &a) noexcept
__host__ __device__ constexpr GpuTuple< detail::tuple_decay_t< Ts >... > makeTuple(Ts &&... args)
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)
__host__ __device__ constexpr const T & min(const T &a, const T &b) noexcept
__host__ __device__ constexpr const T & max(const T &a, const T &b) noexcept
__host__ __device__ 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
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
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::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::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::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 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