WarpX
MusclHancockUtils.H
Go to the documentation of this file.
1 /* Copyright 2023 Grant Johnson, Remi Lehe
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 #ifndef WARPX_MusclHancock_H_
8 #define WARPX_MusclHancock_H_
9 
10 #include <AMReX.H>
11 #include <AMReX_Array4.H>
12 #include <AMReX_Gpu.H>
13 #include <AMReX_REAL.H>
14 
15 
16 // Euler push for momentum source (r-direction)
17 // Note: assumes U normalized by c
19 amrex::Real F_r (amrex::Real r, amrex::Real u_r, amrex::Real u_theta, amrex::Real u_z, amrex::Real dt)
20 {
21  using namespace amrex::literals;
22  return dt*(-u_theta*u_theta/r)/std::sqrt(1.0_rt + u_r*u_r + u_theta*u_theta + u_z*u_z) + u_r;
23 }
24 
25 // Euler push for momentum source (theta-direction)
26 // Note: assumes U normalized by c
28 amrex::Real F_theta (amrex::Real r, amrex::Real u_r, amrex::Real u_theta, amrex::Real u_z, amrex::Real dt)
29 {
30  using namespace amrex::literals;
31  return dt*(u_theta*u_r/r)/std::sqrt(1.0_rt + u_r*u_r + u_theta*u_theta + u_z*u_z) + u_theta;
32 }
33 // Velocity at the half step
35 amrex::Real V_calc (const amrex::Array4<amrex::Real>& U, int i, int j, int k, int comp, amrex::Real c)
36 {
37  using namespace amrex::literals;
38  // comp -> x, y, z -> 0, 1, 2, return Vx, Vy, or Vz:
39  const amrex::Real gamma = std::sqrt(1.0_rt + (U(i,j,k,1)*U(i,j,k,1) + U(i,j,k,2)*U(i,j,k,2) + U(i,j,k,3)*U(i,j,k,3))/(c*c));
40  return U(i,j,k,comp+1)/gamma;
41 }
42 // mindmod
44 amrex::Real minmod (amrex::Real a, amrex::Real b)
45 {
46  using namespace amrex::literals;
47  if (a > 0.0_rt && b > 0.0_rt) {
48  return std::min(a, b);
49  } else if (a < 0.0_rt && b < 0.0_rt) {
50  return std::max(a, b);
51  } else {
52  return 0.0_rt;
53  }
54 }
55 // Min of 3 inputs
57 amrex::Real min3 (amrex::Real a, amrex::Real b, amrex::Real c)
58 {
59  return std::min(a, std::min(b, c) );
60 }
61 // Max of 3 inputs
63 amrex::Real max3 (amrex::Real a, amrex::Real b, amrex::Real c)
64 {
65  return std::max(a, std::max(b, c) );
66 }
67 // mindmod
69 amrex::Real minmod3 (amrex::Real a, amrex::Real b , amrex::Real c)
70 {
71  using namespace amrex::literals;
72  if (a > 0.0_rt && b > 0.0_rt && c > 0.0_rt) {
73  return min3(a,b,c);
74  } else if (a < 0.0_rt && b < 0.0_rt && c < 0.0_rt) {
75  return max3(a,b,c);
76  } else {
77  return 0.0;
78  }
79 }
80 //maxmod
82 amrex::Real maxmod (amrex::Real a, amrex::Real b)
83 {
84  using namespace amrex::literals;
85  if (a > 0.0_rt && b > 0.0_rt) {
86  return std::max(a, b);
87  } else if (a < 0.0_rt && b < 0.0_rt) {
88  return std::min(a, b);
89  } else {
90  return 0.0_rt;
91  }
92 }
93 // Rusanov Flux (density)
96 int i, int j, int k, amrex::Real Vm, amrex::Real Vp)
97 {
98  using namespace amrex::literals;
99  const amrex::Real c = std::max( std::abs(Vm) , std::abs(Vp) );
100  return 0.5_rt*(Vm*Um(i,j,k,0) + Vp*Up(i,j,k,0)) - (0.5_rt*c)*(Up(i,j,k,0) - Um(i,j,k,0));
101 }
102 // Rusanov Flux (Momentum density x)
105 int i, int j, int k, amrex::Real Vm, amrex::Real Vp)
106 {
107  using namespace amrex::literals;
108  const amrex::Real c = std::max( std::abs(Vm) , std::abs(Vp) );
109  return 0.5_rt*(Vm*Um(i,j,k,0)*Um(i,j,k,1) + Vp*Up(i,j,k,0)*Up(i,j,k,1))
110  - (0.5_rt*c)*(Up(i,j,k,0)*Up(i,j,k,1) - Um(i,j,k,0)*Um(i,j,k,1));
111 }
112 // Rusanov Flux (Momentum density y)
115 int i, int j, int k, amrex::Real Vm, amrex::Real Vp)
116 {
117  using namespace amrex::literals;
118  const amrex::Real c = std::max( std::abs(Vm) , std::abs(Vp) );
119  return 0.5_rt*(Vm*Um(i,j,k,0)*Um(i,j,k,2) + Vp*Up(i,j,k,0)*Up(i,j,k,2))
120  - (0.5_rt*c)*(Up(i,j,k,0)*Up(i,j,k,2) - Um(i,j,k,0)*Um(i,j,k,2));
121 }
122 // Rusanov Flux (Momentum density z)
125 int i, int j, int k, amrex::Real Vm, amrex::Real Vp)
126 {
127  using namespace amrex::literals;
128  const amrex::Real c = std::max( std::abs(Vm) , std::abs(Vp) );
129  return 0.5_rt*(Vm*Um(i,j,k,0)*Um(i,j,k,3) + Vp*Up(i,j,k,0)*Up(i,j,k,3))
130  - (0.5_rt*c)*(Up(i,j,k,0)*Up(i,j,k,3) - Um(i,j,k,0)*Um(i,j,k,3));
131 }
132 // ave_minmod high diffusivity, sigma can be between [1,2]
134 amrex::Real ave_adjustable_diff (amrex::Real a, amrex::Real b)
135 {
136  using namespace amrex::literals;
137  constexpr auto sigma = static_cast<amrex::Real>(2.0*0.732050807568877);
138  if (a*b > 0.0_rt) {
139  return minmod3( (a+b)/2.0_rt, sigma*a, sigma*b );
140  } else {
141  return 0.0_rt;
142  }
143 }
144 // ave_minmod Low diffusivity
146 amrex::Real ave (amrex::Real a, amrex::Real b)
147 {
148  using namespace amrex::literals;
149  if (a*b > 0.0_rt) {
150  return minmod3( (a+b)/2.0_rt, 2.0_rt*a, 2.0_rt*b );
151  } else {
152  return 0.0_rt;
153  }
154 }
155 // ave_superbee
157 amrex::Real ave_superbee (amrex::Real a, amrex::Real b)
158 {
159  using namespace amrex::literals;
160  if (a*b > 0.0_rt) {
161  return minmod( maxmod(a,b), minmod(2.0_rt*a,2.0_rt*b));
162  } else {
163  return 0.0_rt;
164  }
165 }
166 // stage2 slope limiting
168 amrex::Real ave_stage2 (amrex::Real dQ, amrex::Real a, amrex::Real b, amrex::Real c)
169 {
170  using namespace amrex::literals;
171  // sigma = sqrt(3) -1
172  constexpr auto sigma = 0.732050807568877_rt;
173  const amrex::Real dq_min = 2.0_rt*std::min( b - min3(a,b,c), max3(a,b,c) - b);
174  return ( std::abs(dQ)/dQ ) * std::min( std::abs(dQ) , sigma*std::abs(dq_min) );
175 }
176 // Returns the offset indices for the "plus" grid
178 void plus_index_offsets (int i, int j, int k, int& ip, int& jp, int& kp, int comp)
179 {
180  using namespace amrex::literals;
181  // Find the correct offsets
182 #if defined(WARPX_DIM_3D)
183  if (comp == 0) { //x
184  ip = i - 1; jp = j; kp = k;
185  } else if (comp == 1){ //y
186  ip = i; jp = j - 1; kp = k;
187  } else if (comp == 2){ //z
188  ip = i; jp = j; kp = k - 1;
189  }
190  #elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
191  if (comp == 0) { //x
192  ip = i - 1; jp = j; kp = k;
193  } else if (comp == 2){ //z
194  ip = i; jp = j - 1; kp = k;
195  }
196 #else
197  if (comp == 2) { //z
198  ip = i - 1; jp = j; kp = k;
199  }
200 #endif
201 }
202 // Compute the zero edges
204 void compute_U_edges (const amrex::Array4<amrex::Real>& Um, const amrex::Array4<amrex::Real>& Up, int i, int j, int k, amrex::Box box,
205 amrex::Real U_tilde0, amrex::Real U_tilde1, amrex::Real U_tilde2, amrex::Real U_tilde3,
206 amrex::Real dU0x, amrex::Real dU1x, amrex::Real dU2x, amrex::Real dU3x, int comp)
207 {
208  using namespace amrex::literals;
209  // comp -> x, y, z -> 0, 1, 2
210  int ip, jp, kp;
211  plus_index_offsets(i, j, k, ip, jp, kp, comp);
212 
213  if ( box.contains(i,j,k) ) {
214  Um(i,j,k,0) = U_tilde0 + dU0x/2.0_rt;
215  Um(i,j,k,1) = U_tilde1 + dU1x/2.0_rt;
216  Um(i,j,k,2) = U_tilde2 + dU2x/2.0_rt;
217  Um(i,j,k,3) = U_tilde3 + dU3x/2.0_rt;
218  }
219 
220  if ( box.contains(ip,jp,kp) ) {
221  Up(ip,jp,kp,0) = U_tilde0 - dU0x/2.0_rt;
222  Up(ip,jp,kp,1) = U_tilde1 - dU1x/2.0_rt;
223  Up(ip,jp,kp,2) = U_tilde2 - dU2x/2.0_rt;
224  Up(ip,jp,kp,3) = U_tilde3 - dU3x/2.0_rt;
225  }
226 }
227 // Compute the zero edges
230 const amrex::Array4<amrex::Real>& Up, int i, int j, int k, amrex::Box box, int comp)
231 {
232  using namespace amrex::literals;
233  // comp -> x, y, z -> 0, 1, 2
234  int ip, jp, kp;
235  plus_index_offsets(i, j, k, ip, jp, kp, comp);
236 
237  if ( box.contains(i,j,k) ) {
238  Um(i,j,k,0) = 0.0_rt;
239  Um(i,j,k,1) = 0.0_rt;
240  Um(i,j,k,2) = 0.0_rt;
241  Um(i,j,k,3) = 0.0_rt;
242  }
243 
244  if ( box.contains(ip,jp,kp) ) {
245  Up(ip,jp,kp,0) = 0.0_rt;
246  Up(ip,jp,kp,1) = 0.0_rt;
247  Up(ip,jp,kp,2) = 0.0_rt;
248  Up(ip,jp,kp,3) = 0.0_rt;
249  }
250 }
251 // Positivity Limiter
252 // if Q_minus or Q_plus is zero for the density (i.e. component 0 of Q_minus/Q_plus), set dQ to 0 and recompute Q_minus / Q_plus
255 const amrex::Array4<amrex::Real>& U_edge_minus, const amrex::Array4<amrex::Real>& N_arr,
256 int i, int j, int k, amrex::Box box, amrex::Real Ux, amrex::Real Uy, amrex::Real Uz,
257 int comp)
258 {
259 
260  using namespace amrex::literals;
261  // comp -> x, y, z -> 0, 1, 2
262  int ip, jp, kp;
263  plus_index_offsets(i, j, k, ip, jp, kp, comp);
264 
265  // Set the edges to zero. If one edge in a cell is zero, we must self-consistently
266  // set the slope to zero (hence why we have the three cases, the first is when
267  // both points exist, and the second two are are edge cases)
268  if (( box.contains(i,j,k) ) && ( box.contains(ip,jp,kp) )) {
269  if ((U_edge_minus(i,j,k,0) < 0.0_rt) || (U_edge_plus(ip,jp,kp,0) < 0.0_rt)) {
270  U_edge_minus(i,j,k,0) = N_arr(i,j,k);
271  U_edge_minus(i,j,k,1) = Ux;
272  U_edge_minus(i,j,k,2) = Uy;
273  U_edge_minus(i,j,k,3) = Uz;
274  U_edge_plus(ip,jp,kp,0) = N_arr(i,j,k);
275  U_edge_plus(ip,jp,kp,1) = Ux;
276  U_edge_plus(ip,jp,kp,2) = Uy;
277  U_edge_plus(ip,jp,kp,3) = Uz;
278  }
279  } else if (( box.contains(i,j,k) ) && ( box.contains(ip,jp,kp) != 1)) {
280  if (U_edge_minus(i,j,k,0) < 0.0_rt) {
281  U_edge_minus(i,j,k,0) = N_arr(i,j,k);
282  U_edge_minus(i,j,k,1) = Ux;
283  U_edge_minus(i,j,k,2) = Uy;
284  U_edge_minus(i,j,k,3) = Uz;
285  }
286  } else if (( box.contains(i,j,k) != 1 ) && ( box.contains(ip,jp,kp) )) {
287  if (U_edge_plus(ip,jp,kp,0) < 0.0_rt){
288  U_edge_plus(ip,jp,kp,0) = N_arr(i,j,k);
289  U_edge_plus(ip,jp,kp,1) = Ux;
290  U_edge_plus(ip,jp,kp,2) = Uy;
291  U_edge_plus(ip,jp,kp,3) = Uz;
292  }
293  }
294 }
295 
296 // Compute the difference in N (down-x)
298 amrex::Real DownDx_N (const amrex::Array4<amrex::Real>& N, int i, int j, int k)
299 {
300  using namespace amrex::literals;
301  // Write the correct differences
302 #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
303  return N(i,j,k) - N(i-1,j,k);
304 #else
305  amrex::ignore_unused(N, i, j, k);
306  return 0.0_rt;
307 #endif
308 }
309 // Compute the difference in N (up-x)
311 amrex::Real UpDx_N (const amrex::Array4<amrex::Real>& N, int i, int j, int k)
312 {
313  using namespace amrex::literals;
314  // Write the correct differences
315 #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
316  return N(i+1,j,k) - N(i,j,k);
317 #else
318  amrex::ignore_unused(N, i, j, k);
319  return 0.0_rt;
320 #endif
321 }
322 // Compute the difference in N (down-y)
324 amrex::Real DownDy_N (const amrex::Array4<amrex::Real>& N, int i, int j, int k)
325 {
326  using namespace amrex::literals;
327  // Write the correct differences
328 #if defined(WARPX_DIM_3D)
329  return N(i,j,k) - N(i,j-1,k);
330 #else
331  amrex::ignore_unused(N, i, j, k);
332  return 0.0_rt;
333 #endif
334 }
335 // Compute the difference in N (up-y)
337 amrex::Real UpDy_N (const amrex::Array4<amrex::Real>& N, int i, int j, int k)
338 {
339  using namespace amrex::literals;
340  // Write the correct differences
341 #if defined(WARPX_DIM_3D)
342  return N(i,j+1,k) - N(i,j,k);
343 #else
344  amrex::ignore_unused(N, i, j, k);
345  return 0.0_rt;
346 #endif
347 }
348 // Compute the difference in N (down-z)
350 amrex::Real DownDz_N (const amrex::Array4<amrex::Real>& N, int i, int j, int k)
351 {
352  using namespace amrex::literals;
353  // Write the correct differences
354 #if defined(WARPX_DIM_3D)
355  return N(i,j,k) - N(i,j,k-1);
356 #elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
357  return N(i,j,k) - N(i,j-1,k);
358 #else
359  return N(i,j,k) - N(i-1,j,k);
360 #endif
361 }
362 // Compute the difference in N (up-z)
364 amrex::Real UpDz_N (const amrex::Array4<amrex::Real>& N, int i, int j, int k)
365 {
366  using namespace amrex::literals;
367  // Write the correct differences
368 #if defined(WARPX_DIM_3D)
369  return N(i,j,k+1) - N(i,j,k);
370 #elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
371  return N(i,j+1,k) - N(i,j,k);
372 #else
373  return N(i+1,j,k) - N(i,j,k);
374 #endif
375 }
376 
377 
378 // Compute the difference in U (down-x)
380 amrex::Real DownDx_U (const amrex::Array4<amrex::Real>& N,
381 const amrex::Array4<amrex::Real>& NU, amrex::Real& U, int i, int j, int k)
382 {
383  using namespace amrex::literals;
384  // Write the correct differences
385 #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
386  // U is zero if N is zero, Check positivity before dividing
387  amrex::Real U_m = 0;
388  if (N(i-1,j,k) > 0) { U_m = NU(i-1,j,k)/N(i-1,j,k); }
389  return U - U_m;
390 #else
391  amrex::ignore_unused(N, NU, U, i, j, k);
392  return 0.0_rt;
393 #endif
394 }
395 // Compute the difference in U (up-x)
397 amrex::Real UpDx_U (const amrex::Array4<amrex::Real>& N,
398 const amrex::Array4<amrex::Real>& NU, amrex::Real& U, int i, int j, int k)
399 {
400  using namespace amrex::literals;
401  // Write the correct differences
402 #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
403  // U is zero if N is zero, Check positivity before dividing
404  amrex::Real U_p = 0;
405  if (N(i+1,j,k) > 0) { U_p = NU(i+1,j,k)/N(i+1,j,k); }
406  return U_p - U;
407 #else
408  amrex::ignore_unused(N, NU, U, i, j, k);
409  return 0.0_rt;
410 #endif
411 }
412 
413 // Compute the difference in U (down-y)
415 amrex::Real DownDy_U (const amrex::Array4<amrex::Real>& N,
416 const amrex::Array4<amrex::Real>& NU, amrex::Real& U, int i, int j, int k)
417 {
418  using namespace amrex::literals;
419  // Write the correct differences
420 #if defined(WARPX_DIM_3D)
421  // U is zero if N is zero, Check positivity before dividing
422  amrex::Real U_m = 0;
423  if (N(i,j-1,k) > 0) { U_m = NU(i,j-1,k)/N(i,j-1,k); }
424  return U - U_m;
425 #else
426  amrex::ignore_unused(N, NU, U, i, j, k);
427  return 0.0_rt;
428 #endif
429 }
430 // Compute the difference in U (up-y)
432 amrex::Real UpDy_U (const amrex::Array4<amrex::Real>& N,
433 const amrex::Array4<amrex::Real>& NU, amrex::Real& U, int i, int j, int k)
434 {
435  using namespace amrex::literals;
436  // Write the correct differences
437 #if defined(WARPX_DIM_3D)
438  // U is zero if N is zero, Check positivity before dividing
439  amrex::Real U_p = 0;
440  if (N(i,j+1,k) > 0) { U_p = NU(i,j+1,k)/N(i,j+1,k); }
441  return U_p - U;
442 #else
443  amrex::ignore_unused(N, NU, U, i, j, k);
444  return 0.0_rt;
445 #endif
446 }
447 
448 // Compute the difference in U (down-z)
450 amrex::Real DownDz_U (const amrex::Array4<amrex::Real>& N,
451 const amrex::Array4<amrex::Real>& NU, amrex::Real& U, int i, int j, int k)
452 {
453  using namespace amrex::literals;
454  // Write the correct differences
455  amrex::Real U_m = 0_rt;
456 
457  // U is zero if N is zero, Check positivity before dividing
458 #if defined(WARPX_DIM_3D)
459  if (N(i,j,k-1) > 0) { U_m = NU(i,j,k-1)/N(i,j,k-1); }
460 #elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
461  if (N(i,j-1,k) > 0) { U_m = NU(i,j-1,k)/N(i,j-1,k); }
462 #else
463  if (N(i-1,j,k) > 0) { U_m = NU(i-1,j,k)/N(i-1,j,k); }
464 #endif
465 
466  // Return the difference
467  return U - U_m;
468 }
469 // Compute the difference in U (up-z)
471 amrex::Real UpDz_U (const amrex::Array4<amrex::Real>& N,
472 const amrex::Array4<amrex::Real>& NU, amrex::Real& U, int i, int j, int k)
473 {
474  using namespace amrex::literals;
475  // Write the correct differences
476  amrex::Real U_p = 0;
477 
478  // U is zero if N is zero, Check positivity before dividing
479 #if defined(WARPX_DIM_3D)
480  if (N(i,j,k+1) > 0) { U_p = NU(i,j,k+1)/N(i,j,k+1); }
481 #elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
482  if (N(i,j+1,k) > 0) { U_p = NU(i,j+1,k)/N(i,j+1,k); }
483 #else
484  if (N(i+1,j,k) > 0) { U_p = NU(i+1,j,k)/N(i+1,j,k); }
485 #endif
486 
487  // Return the difference
488  return U_p - U;
489 }
490 
491 
492 // Flux difference calculation
494 amrex::Real dF (const amrex::Array4<amrex::Real>& U_minus,
495 const amrex::Array4<amrex::Real>& U_plus,int i,int j,int k,amrex::Real clight, int comp, int dir)
496 {
497  using namespace amrex::literals;
498  // dir -> x, y, z -> 0, 1, 2
499  int ip, jp, kp;
500  plus_index_offsets(i, j, k, ip, jp, kp, dir);
501 
502  const amrex::Real V_L_minus = V_calc(U_minus,ip,jp,kp,dir,clight);
503  const amrex::Real V_I_minus = V_calc(U_minus,i,j,k,dir,clight);
504  const amrex::Real V_L_plus = V_calc(U_plus,ip,jp,kp,dir,clight);
505  const amrex::Real V_I_plus = V_calc(U_plus,i,j,k,dir,clight);
506 
507  // Flux differences depending on the component to compute
508  if (comp == 0){
509  return flux_N( U_minus, U_plus, i, j, k, V_I_minus, V_I_plus) - flux_N( U_minus, U_plus, ip, jp, kp, V_L_minus, V_L_plus);
510  } else if (comp == 1){
511  return flux_NUx( U_minus, U_plus, i, j, k, V_I_minus, V_I_plus) - flux_NUx( U_minus, U_plus, ip, jp, kp, V_L_minus, V_L_plus);
512  } else if (comp == 2){
513  return flux_NUy( U_minus, U_plus, i, j, k, V_I_minus, V_I_plus) - flux_NUy( U_minus, U_plus, ip, jp, kp, V_L_minus, V_L_plus);
514  } else { //if (comp == 3)
515  return flux_NUz( U_minus, U_plus, i, j, k, V_I_minus, V_I_plus) - flux_NUz( U_minus, U_plus, ip, jp, kp, V_L_minus, V_L_plus);
516  }
517 }
518 
519 #endif /*WARPX_MusclHancock_H_*/
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real max3(amrex::Real a, amrex::Real b, amrex::Real c)
Definition: MusclHancockUtils.H:63
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real DownDx_U(const amrex::Array4< amrex::Real > &N, const amrex::Array4< amrex::Real > &NU, amrex::Real &U, int i, int j, int k)
Definition: MusclHancockUtils.H:380
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void positivity_limiter(const amrex::Array4< amrex::Real > &U_edge_plus, const amrex::Array4< amrex::Real > &U_edge_minus, const amrex::Array4< amrex::Real > &N_arr, int i, int j, int k, amrex::Box box, amrex::Real Ux, amrex::Real Uy, amrex::Real Uz, int comp)
Definition: MusclHancockUtils.H:254
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real dF(const amrex::Array4< amrex::Real > &U_minus, const amrex::Array4< amrex::Real > &U_plus, int i, int j, int k, amrex::Real clight, int comp, int dir)
Definition: MusclHancockUtils.H:494
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void compute_U_edges(const amrex::Array4< amrex::Real > &Um, const amrex::Array4< amrex::Real > &Up, int i, int j, int k, amrex::Box box, amrex::Real U_tilde0, amrex::Real U_tilde1, amrex::Real U_tilde2, amrex::Real U_tilde3, amrex::Real dU0x, amrex::Real dU1x, amrex::Real dU2x, amrex::Real dU3x, int comp)
Definition: MusclHancockUtils.H:204
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real DownDz_N(const amrex::Array4< amrex::Real > &N, int i, int j, int k)
Definition: MusclHancockUtils.H:350
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real UpDx_U(const amrex::Array4< amrex::Real > &N, const amrex::Array4< amrex::Real > &NU, amrex::Real &U, int i, int j, int k)
Definition: MusclHancockUtils.H:397
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real UpDy_N(const amrex::Array4< amrex::Real > &N, int i, int j, int k)
Definition: MusclHancockUtils.H:337
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real V_calc(const amrex::Array4< amrex::Real > &U, int i, int j, int k, int comp, amrex::Real c)
Definition: MusclHancockUtils.H:35
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real minmod3(amrex::Real a, amrex::Real b, amrex::Real c)
Definition: MusclHancockUtils.H:69
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ave_superbee(amrex::Real a, amrex::Real b)
Definition: MusclHancockUtils.H:157
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real UpDy_U(const amrex::Array4< amrex::Real > &N, const amrex::Array4< amrex::Real > &NU, amrex::Real &U, int i, int j, int k)
Definition: MusclHancockUtils.H:432
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ave(amrex::Real a, amrex::Real b)
Definition: MusclHancockUtils.H:146
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real UpDx_N(const amrex::Array4< amrex::Real > &N, int i, int j, int k)
Definition: MusclHancockUtils.H:311
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void plus_index_offsets(int i, int j, int k, int &ip, int &jp, int &kp, int comp)
Definition: MusclHancockUtils.H:178
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real minmod(amrex::Real a, amrex::Real b)
Definition: MusclHancockUtils.H:44
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real DownDy_U(const amrex::Array4< amrex::Real > &N, const amrex::Array4< amrex::Real > &NU, amrex::Real &U, int i, int j, int k)
Definition: MusclHancockUtils.H:415
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real min3(amrex::Real a, amrex::Real b, amrex::Real c)
Definition: MusclHancockUtils.H:57
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real flux_NUz(const amrex::Array4< amrex::Real > &Um, const amrex::Array4< amrex::Real > &Up, int i, int j, int k, amrex::Real Vm, amrex::Real Vp)
Definition: MusclHancockUtils.H:124
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ave_adjustable_diff(amrex::Real a, amrex::Real b)
Definition: MusclHancockUtils.H:134
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real DownDx_N(const amrex::Array4< amrex::Real > &N, int i, int j, int k)
Definition: MusclHancockUtils.H:298
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real F_theta(amrex::Real r, amrex::Real u_r, amrex::Real u_theta, amrex::Real u_z, amrex::Real dt)
Definition: MusclHancockUtils.H:28
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real ave_stage2(amrex::Real dQ, amrex::Real a, amrex::Real b, amrex::Real c)
Definition: MusclHancockUtils.H:168
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real UpDz_N(const amrex::Array4< amrex::Real > &N, int i, int j, int k)
Definition: MusclHancockUtils.H:364
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real maxmod(amrex::Real a, amrex::Real b)
Definition: MusclHancockUtils.H:82
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real flux_NUy(const amrex::Array4< amrex::Real > &Um, const amrex::Array4< amrex::Real > &Up, int i, int j, int k, amrex::Real Vm, amrex::Real Vp)
Definition: MusclHancockUtils.H:114
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real flux_N(const amrex::Array4< amrex::Real > &Um, const amrex::Array4< amrex::Real > &Up, int i, int j, int k, amrex::Real Vm, amrex::Real Vp)
Definition: MusclHancockUtils.H:95
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real DownDy_N(const amrex::Array4< amrex::Real > &N, int i, int j, int k)
Definition: MusclHancockUtils.H:324
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void set_U_edges_to_zero(const amrex::Array4< amrex::Real > &Um, const amrex::Array4< amrex::Real > &Up, int i, int j, int k, amrex::Box box, int comp)
Definition: MusclHancockUtils.H:229
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real flux_NUx(const amrex::Array4< amrex::Real > &Um, const amrex::Array4< amrex::Real > &Up, int i, int j, int k, amrex::Real Vm, amrex::Real Vp)
Definition: MusclHancockUtils.H:104
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real UpDz_U(const amrex::Array4< amrex::Real > &N, const amrex::Array4< amrex::Real > &NU, amrex::Real &U, int i, int j, int k)
Definition: MusclHancockUtils.H:471
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real F_r(amrex::Real r, amrex::Real u_r, amrex::Real u_theta, amrex::Real u_z, amrex::Real dt)
Definition: MusclHancockUtils.H:19
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real DownDz_U(const amrex::Array4< amrex::Real > &N, const amrex::Array4< amrex::Real > &NU, amrex::Real &U, int i, int j, int k)
Definition: MusclHancockUtils.H:450
AMREX_GPU_HOST_DEVICE bool contains(const IntVect &p) const noexcept
static constexpr auto c
vacuum speed of light [m/s]
Definition: constant.H:44
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
i
Definition: check_interp_points_and_weights.py:174
float dt
Definition: stencil.py:442
int gamma
boosted frame
Definition: stencil.py:431
float clight
Definition: video_yt.py:41