From 8828c4bf1950ea37d5da36d4073af63eb2537bc0 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 27 Jul 2024 10:46:45 -0400 Subject: [PATCH] fix pressure source --- Source/hydro/reconstruction.H | 52 +++++++++++++++++++++-------------- 1 file changed, 32 insertions(+), 20 deletions(-) diff --git a/Source/hydro/reconstruction.H b/Source/hydro/reconstruction.H index 65d6fc31dc..9a41710663 100644 --- a/Source/hydro/reconstruction.H +++ b/Source/hydro/reconstruction.H @@ -1,6 +1,8 @@ #ifndef CASTRO_RECONSTRUCTION_H #define CASTRO_RECONSTRUCTION_H +#include + namespace reconstruction { enum slope_indices { im2 = 0, @@ -14,9 +16,9 @@ namespace reconstruction { AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -load_stencil(Array4 const& q_arr, const int idir, +load_stencil(amrex::Array4 const& q_arr, const int idir, const int i, const int j, const int k, const int ncomp, - Real* s) { + amrex::Real* s) { using namespace reconstruction; @@ -47,9 +49,11 @@ load_stencil(Array4 const& q_arr, const int idir, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -load_passive_stencil(Array4 const& U_arr, Array4 const& rho_inv_arr, const int idir, +load_passive_stencil(amrex::Array4 const& U_arr, + amrex::Array4 const& rho_inv_arr, + const int idir, const int i, const int j, const int k, const int ncomp, - Real* s) { + amrex::Real* s) { using namespace reconstruction; @@ -80,13 +84,16 @@ load_passive_stencil(Array4 const& U_arr, Array4 const& AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -add_geometric_rho_source(Array4 const& q_arr, - Array4 const& dloga, +add_geometric_rho_source(amrex::Array4 const& q_arr, + amrex::Array4 const& dloga, const int i, const int j, const int k, - Real* s) { + amrex::Real* s) { using namespace reconstruction; + // this takes the form: -alpha rho u / r + // where alpha = 1 for cylindrical and 2 for spherical + // note: this is assumed to be working only in the x-direction s[im2] += -dloga(i-2,j,k) * q_arr(i-2,j,k,QRHO) * q_arr(i-2,j,k,QU); @@ -98,13 +105,16 @@ add_geometric_rho_source(Array4 const& q_arr, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -add_geometric_rhoe_source(Array4 const& q_arr, - Array4 const& dloga, +add_geometric_rhoe_source(amrex::Array4 const& q_arr, + amrex::Array4 const& dloga, const int i, const int j, const int k, - Real* s) { + amrex::Real* s) { using namespace reconstruction; + // this takes the form: -alpha (rho e + p) u / r + // where alpha = 1 for cylindrical and 2 for spherical + // note: this is assumed to be working only in the x-direction s[im2] += -dloga(i-2,j,k) * (q_arr(i-2,j,k,QREINT) + q_arr(i-2,j,k,QPRES)) * q_arr(i-2,j,k,QU); @@ -116,23 +126,25 @@ add_geometric_rhoe_source(Array4 const& q_arr, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -add_geometric_p_source(Array4 const& q_arr, - Array4 const& qaux_arr, - Array4 const& dloga, +add_geometric_p_source(amrex::Array4 const& q_arr, + amrex::Array4 const& qaux_arr, + amrex::Array4 const& dloga, const int i, const int j, const int k, - Real* s) { + amrex::Real* s) { using namespace reconstruction; + // this takes the form: -alpha Gamma1 p u / r + // where alpha = 1 for cylindrical and 2 for spherical + // note: this is assumed to be working only in the x-direction - s[im2] += -dloga(i-2,j,k) * q_arr(i-2,j,k,QRHO) * qaux_arr(i-2,j,k,QC) * q_arr(i-2,j,k,QU); - s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QRHO) * qaux_arr(i-1,j,k,QC) * q_arr(i-1,j,k,QU); - s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QRHO) * qaux_arr(i,j,k,QC) * q_arr(i,j,k,QU); - s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QRHO) * qaux_arr(i+1,j,k,QC) * q_arr(i+1,j,k,QU); - s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QRHO) * qaux_arr(i+2,j,k,QC) * q_arr(i+2,j,k,QU); + s[im2] += -dloga(i-2,j,k) * q_arr(i-2,j,k,QPRES) * qaux_arr(i-2,j,k,QGAMC) * q_arr(i-2,j,k,QU); + s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QPRES) * qaux_arr(i-1,j,k,QGAMC) * q_arr(i-1,j,k,QU); + s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QPRES) * qaux_arr(i,j,k,QGAMC) * q_arr(i,j,k,QU); + s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QPRES) * qaux_arr(i+1,j,k,QGAMC) * q_arr(i+1,j,k,QU); + s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QPRES) * qaux_arr(i+2,j,k,QGAMC) * q_arr(i+2,j,k,QU); } #endif -