7 #ifndef WARPX_MusclHancock_H_
8 #define WARPX_MusclHancock_H_
19 amrex::Real
F_r (amrex::Real r, amrex::Real u_r, amrex::Real u_theta, amrex::Real u_z, amrex::Real
dt)
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;
28 amrex::Real
F_theta (amrex::Real r, amrex::Real u_r, amrex::Real u_theta, amrex::Real u_z, amrex::Real
dt)
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;
37 using namespace amrex::literals;
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;
44 amrex::Real
minmod (amrex::Real a, amrex::Real b)
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);
57 amrex::Real
min3 (amrex::Real a, amrex::Real b, amrex::Real
c)
59 return std::min(a, std::min(
b,
c) );
63 amrex::Real
max3 (amrex::Real a, amrex::Real b, amrex::Real
c)
65 return std::max(a, std::max(
b,
c) );
69 amrex::Real
minmod3 (amrex::Real a, amrex::Real b , amrex::Real
c)
71 using namespace amrex::literals;
72 if (a > 0.0_rt &&
b > 0.0_rt &&
c > 0.0_rt) {
74 }
else if (a < 0.0_rt &&
b < 0.0_rt &&
c < 0.0_rt) {
82 amrex::Real
maxmod (amrex::Real a, amrex::Real b)
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);
96 int i,
int j,
int k, amrex::Real Vm, amrex::Real Vp)
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));
105 int i,
int j,
int k, amrex::Real Vm, amrex::Real Vp)
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));
115 int i,
int j,
int k, amrex::Real Vm, amrex::Real Vp)
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));
125 int i,
int j,
int k, amrex::Real Vm, amrex::Real Vp)
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));
136 using namespace amrex::literals;
137 constexpr
auto sigma =
static_cast<amrex::Real
>(2.0*0.732050807568877);
139 return minmod3( (a+
b)/2.0_rt, sigma*a, sigma*
b );
146 amrex::Real
ave (amrex::Real a, amrex::Real b)
148 using namespace amrex::literals;
150 return minmod3( (a+
b)/2.0_rt, 2.0_rt*a, 2.0_rt*
b );
159 using namespace amrex::literals;
168 amrex::Real
ave_stage2 (amrex::Real dQ, amrex::Real a, amrex::Real b, amrex::Real
c)
170 using namespace amrex::literals;
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) );
180 using namespace amrex::literals;
182 #if defined(WARPX_DIM_3D)
184 ip =
i - 1; jp = j; kp = k;
185 }
else if (comp == 1){
186 ip =
i; jp = j - 1; kp = k;
187 }
else if (comp == 2){
188 ip =
i; jp = j; kp = k - 1;
190 #elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
192 ip =
i - 1; jp = j; kp = k;
193 }
else if (comp == 2){
194 ip =
i; jp = j - 1; kp = k;
198 ip =
i - 1; jp = j; kp = k;
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)
208 using namespace amrex::literals;
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;
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;
232 using namespace amrex::literals;
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;
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;
256 int i,
int j,
int k,
amrex::Box box, amrex::Real Ux, amrex::Real Uy, amrex::Real Uz,
260 using namespace amrex::literals;
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;
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;
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;
300 using namespace amrex::literals;
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);
313 using namespace amrex::literals;
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);
326 using namespace amrex::literals;
328 #if defined(WARPX_DIM_3D)
329 return N(
i,j,k) - N(
i,j-1,k);
339 using namespace amrex::literals;
341 #if defined(WARPX_DIM_3D)
342 return N(
i,j+1,k) - N(
i,j,k);
352 using namespace amrex::literals;
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);
359 return N(
i,j,k) - N(
i-1,j,k);
366 using namespace amrex::literals;
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);
373 return N(
i+1,j,k) - N(
i,j,k);
383 using namespace amrex::literals;
385 #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
388 if (N(
i-1,j,k) > 0) { U_m = NU(
i-1,j,k)/N(
i-1,j,k); }
400 using namespace amrex::literals;
402 #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
405 if (N(
i+1,j,k) > 0) { U_p = NU(
i+1,j,k)/N(
i+1,j,k); }
418 using namespace amrex::literals;
420 #if defined(WARPX_DIM_3D)
423 if (N(
i,j-1,k) > 0) { U_m = NU(
i,j-1,k)/N(
i,j-1,k); }
435 using namespace amrex::literals;
437 #if defined(WARPX_DIM_3D)
440 if (N(
i,j+1,k) > 0) { U_p = NU(
i,j+1,k)/N(
i,j+1,k); }
453 using namespace amrex::literals;
455 amrex::Real U_m = 0_rt;
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); }
463 if (N(
i-1,j,k) > 0) { U_m = NU(
i-1,j,k)/N(
i-1,j,k); }
474 using namespace amrex::literals;
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); }
484 if (N(
i+1,j,k) > 0) { U_p = NU(
i+1,j,k)/N(
i+1,j,k); }
497 using namespace amrex::literals;
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);
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);
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);
#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