diff --git a/Source/hydro/reconstruction.H b/Source/hydro/reconstruction.H index 9a41710663..50c2a95062 100644 --- a/Source/hydro/reconstruction.H +++ b/Source/hydro/reconstruction.H @@ -87,20 +87,27 @@ void add_geometric_rho_source(amrex::Array4 const& q_arr, amrex::Array4 const& dloga, const int i, const int j, const int k, - amrex::Real* s) { + const int ncomp , amrex::Real* s) { using namespace reconstruction; + // For idir == 0, i.e. r-direction: // this takes the form: -alpha rho u / r // where alpha = 1 for cylindrical and 2 for spherical + // note dloga(idir==0) = alpha/r - // note: this is assumed to be working only in the x-direction + // For idir == 1, + // i.e. theta-direction for spherical and z-direction for cylindrical + // this takes the form: -rho v cot(theta) / r for spherical and 0 for cylindrical + // note: dloga(idir==1) = cot(theta)/r for spherical and 0 for cylindrical. + + // ncomp should be QU for idir == 0 and QV for idir == 1 - s[im2] += -dloga(i-2,j,k) * q_arr(i-2,j,k,QRHO) * q_arr(i-2,j,k,QU); - s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QRHO) * q_arr(i-1,j,k,QU); - s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QRHO) * q_arr(i,j,k,QU); - s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QRHO) * q_arr(i+1,j,k,QU); - s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QRHO) * q_arr(i+2,j,k,QU); + s[im2] += -dloga(i-2,j,k) * q_arr(i-2,j,k,QRHO) * q_arr(i-2,j,k,ncomp); + s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QRHO) * q_arr(i-1,j,k,ncomp); + s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QRHO) * q_arr(i,j,k,ncomp); + s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QRHO) * q_arr(i+1,j,k,ncomp); + s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QRHO) * q_arr(i+2,j,k,ncomp); } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -108,20 +115,27 @@ void add_geometric_rhoe_source(amrex::Array4 const& q_arr, amrex::Array4 const& dloga, const int i, const int j, const int k, - amrex::Real* s) { + const int ncomp, amrex::Real* s) { using namespace reconstruction; - // this takes the form: -alpha (rho e + p) u / r + // For idir == 0, i.e. r-direction: + // this takes the form: -alpha (rho e + p ) u / r // where alpha = 1 for cylindrical and 2 for spherical + // note dloga(idir==0) = alpha/r - // note: this is assumed to be working only in the x-direction + // For idir == 1, + // i.e. theta-direction for spherical and z-direction for cylindrical + // this takes the form: - (rho e + p) v cot(theta) / r for spherical and 0 for cylindrical + // note: dloga(idir==1) = cot(theta)/r for spherical and 0 for cylindrical. - 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); - s[im1] += -dloga(i-1,j,k) * (q_arr(i-1,j,k,QREINT) + q_arr(i-1,j,k,QPRES)) * q_arr(i-1,j,k,QU); - s[i0] += -dloga(i,j,k) * (q_arr(i,j,k,QREINT) + q_arr(i,j,k,QPRES)) * q_arr(i,j,k,QU); - s[ip1] += -dloga(i+1,j,k) * (q_arr(i+1,j,k,QREINT) + q_arr(i+1,j,k,QPRES)) * q_arr(i+1,j,k,QU); - s[ip2] += -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); + // ncomp should be QU for idir == 0 and QV for idir == 1 + + 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,ncomp); + s[im1] += -dloga(i-1,j,k) * (q_arr(i-1,j,k,QREINT) + q_arr(i-1,j,k,QPRES)) * q_arr(i-1,j,k,comp); + s[i0] += -dloga(i,j,k) * (q_arr(i,j,k,QREINT) + q_arr(i,j,k,QPRES)) * q_arr(i,j,k,ncomp); + s[ip1] += -dloga(i+1,j,k) * (q_arr(i+1,j,k,QREINT) + q_arr(i+1,j,k,QPRES)) * q_arr(i+1,j,k,ncomp); + s[ip2] += -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,ncomp); } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -130,20 +144,32 @@ 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, - amrex::Real* s) { + const int ncomp amrex::Real* s) { using namespace reconstruction; + // For idir == 0, i.e. r-direction: + // this takes the form: -alpha Gamma1 p u / r + // where alpha = 1 for cylindrical and 2 for spherical + // note: dloga(idir==0) = alpha/r + + // For idir == 1, + // i.e. theta-direction for spherical and z-direction for cylindrical + // this takes the form: - Gamma1 p v cot(theta) / r for spherical and 0 for cylindrical + // note: dloga(idir==1) = cot(theta)/r for spherical and 0 for cylindrical. + + // ncomp should be QU for idir == 0 and QV for idir == 1 + // 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,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); + 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,ncomp); + 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,ncomp); + s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QPRES) * qaux_arr(i,j,k,QGAMC) * q_arr(i,j,k,ncomp); + 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,ncomp); + 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,ncomp); } diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index 0795500825..63b160e32e 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -55,10 +55,12 @@ Castro::trace_ppm(const Box& bx, const auto dx = geom.CellSizeArray(); + const auto problo = geom.ProbLoArray(); const int coord = geom.Coord(); Real hdt = 0.5_rt * dt; - Real dtdx = dt / dx[idir]; + Real dtdL = dt / dx[idir]; + Real dL = dx[idir]; #ifndef AMREX_USE_GPU auto lo = bx.loVect3d(); @@ -84,7 +86,8 @@ Castro::trace_ppm(const Box& bx, do_source_trace[n] = 0; // geometric source terms in r-direction need tracing - if (coord > 0 && idir == 0 && (n == QRHO || n == QPRES || n == QREINT)) { + if ((coord == 2 || (idir == 0 && coord == 1)) + && (n == QRHO || n == QPRES || n == QREINT)) { do_source_trace[n] = 1; continue; } @@ -154,6 +157,12 @@ Castro::trace_ppm(const Box& bx, amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + // Want dt/(rdtheta) instead of dt/dtheta for 2d Spherical + if (coord == 2 && idir == 1) { + Real r = problo[0] + static_cast(i + 0.5_rt) * dx[0]; + dL = r * dx[1]; + dtdL = dt / dL; + } Real cc = qaux_arr(i,j,k,QC); Real un = q_arr(i,j,k,QUN); @@ -196,7 +205,7 @@ Castro::trace_ppm(const Box& bx, load_stencil(q_arr, idir, i, j, k, QRHO, s); ppm_reconstruct(s, flat, sm, sp); - ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_rho, Im_rho); + ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_rho, Im_rho); // reconstruct normal velocity @@ -207,8 +216,8 @@ Castro::trace_ppm(const Box& bx, load_stencil(q_arr, idir, i, j, k, QUN, s); ppm_reconstruct(s, flat, sm, sp); - ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_un_0, Im_un_0); - ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_un_2, Im_un_2); + ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdL, Ip_un_0, Im_un_0); + ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdL, Ip_un_2, Im_un_2); // reconstruct pressure @@ -224,12 +233,12 @@ Castro::trace_ppm(const Box& bx, load_stencil(q_arr, idir, i, j, k, QRHO, trho); load_stencil(srcQ, idir, i, j, k, QUN, src); - ppm_reconstruct_pslope(trho, s, src, flat, dx[idir], sm, sp); + ppm_reconstruct_pslope(trho, s, src, flat, dL, sm, sp); } else { ppm_reconstruct(s, flat, sm, sp); } - ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_p, Im_p); + ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_p, Im_p); // reconstruct rho e @@ -238,7 +247,7 @@ Castro::trace_ppm(const Box& bx, load_stencil(q_arr, idir, i, j, k, QREINT, s); ppm_reconstruct(s, flat, sm, sp); - ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_rhoe, Im_rhoe); + ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_rhoe, Im_rhoe); // reconstruct transverse velocities @@ -249,11 +258,11 @@ Castro::trace_ppm(const Box& bx, load_stencil(q_arr, idir, i, j, k, QUT, s); ppm_reconstruct(s, flat, sm, sp); - ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_ut_1, Im_ut_1); + ppm_int_profile_single(sm, sp, s[i0], un, dtdL, Ip_ut_1, Im_ut_1); load_stencil(q_arr, idir, i, j, k, QUTT, s); ppm_reconstruct(s, flat, sm, sp); - ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_utt_1, Im_utt_1); + ppm_int_profile_single(sm, sp, s[i0], un, dtdL, Ip_utt_1, Im_utt_1); // gamma_c @@ -264,8 +273,8 @@ Castro::trace_ppm(const Box& bx, load_stencil(qaux_arr, idir, i, j, k, QGAMC, s); ppm_reconstruct(s, flat, sm, sp); - ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_gc_0, Im_gc_0); - ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_gc_2, Im_gc_2); + ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdL, Ip_gc_0, Im_gc_0); + ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdL, Ip_gc_2, Im_gc_2); // source terms -- we only trace if the terms in the stencil are non-zero @@ -279,7 +288,7 @@ Castro::trace_ppm(const Box& bx, #ifndef AMREX_USE_GPU do_trace = do_source_trace[QRHO]; #else - if (idir == 0 && coord > 0) { + if (coord == 2 || (idir == 0 && coord == 1)) { do_trace = 1; } else { do_trace = check_trace_source(srcQ, idir, i, j, k, QRHO); @@ -289,12 +298,12 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QRHO, s); #if AMREX_SPACEDIM <= 2 - if (idir == 0 && coord > 0) { - add_geometric_rho_source(q_arr, dloga, i, j, k, s); + if (coord == 2 || (idir == 0 && coord == 1)) { + add_geometric_rho_source(q_arr, dloga, i, j, k, QUN, s); } #endif ppm_reconstruct(s, flat, sm, sp); - ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_rho, Im_src_rho); + ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_src_rho, Im_src_rho); } // normal velocity @@ -313,8 +322,8 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUN, s); ppm_reconstruct(s, flat, sm, sp); - ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_src_un_0, Im_src_un_0); - ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_src_un_2, Im_src_un_2); + ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdL, Ip_src_un_0, Im_src_un_0); + ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdL, Ip_src_un_2, Im_src_un_2); } // pressure @@ -325,7 +334,7 @@ Castro::trace_ppm(const Box& bx, #ifndef AMREX_USE_GPU do_trace = do_source_trace[QPRES]; #else - if (idir == 0 && coord > 0) { + if (coord == 2 || (idir == 0 && coord == 1)) { do_trace = 1; } else { do_trace = check_trace_source(srcQ, idir, i, j, k, QPRES); @@ -335,12 +344,12 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QPRES, s); #if AMREX_SPACEDIM <= 2 - if (idir == 0 && coord > 0) { - add_geometric_p_source(q_arr, qaux_arr, dloga, i, j, k, s); + if (coord == 2 || (idir == 0 && coord == 1)) { + add_geometric_p_source(q_arr, qaux_arr, dloga, i, j, k, QUN, s); } #endif ppm_reconstruct(s, flat, sm, sp); - ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_p, Im_src_p); + ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_src_p, Im_src_p); } // rho e @@ -351,7 +360,7 @@ Castro::trace_ppm(const Box& bx, #ifndef AMREX_USE_GPU do_trace = do_source_trace[QREINT]; #else - if (idir == 0 && coord > 0) { + if (coord == 2 || (idir == 0 && coord == 1)) { do_trace = 1; } else { do_trace = check_trace_source(srcQ, idir, i, j, k, QREINT); @@ -361,12 +370,12 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QREINT, s); #if AMREX_SPACEDIM <= 2 - if (idir == 0 && coord > 0) { - add_geometric_rhoe_source(q_arr, dloga, i, j, k, s); + if (coord == 2 || (idir == 0 && coord == 1)) { + add_geometric_rhoe_source(q_arr, dloga, i, j, k, QUN, s); } #endif ppm_reconstruct(s, flat, sm, sp); - ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_rhoe, Im_src_rhoe); + ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_src_rhoe, Im_src_rhoe); } // transverse velocities @@ -383,7 +392,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUT, s); ppm_reconstruct(s, flat, sm, sp); - ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_src_ut_1, Im_src_ut_1); + ppm_int_profile_single(sm, sp, s[i0], un, dtdL, Ip_src_ut_1, Im_src_ut_1); } Real Ip_src_utt_1 = 0.0_rt; @@ -398,7 +407,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUTT, s); ppm_reconstruct(s, flat, sm, sp); - ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_src_utt_1, Im_src_utt_1); + ppm_int_profile_single(sm, sp, s[i0], un, dtdL, Ip_src_utt_1, Im_src_utt_1); } @@ -416,7 +425,7 @@ Castro::trace_ppm(const Box& bx, load_passive_stencil(U_arr, rho_inv_arr, idir, i, j, k, nc, s); ppm_reconstruct(s, flat, sm, sp); - ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_passive, Im_passive); + ppm_int_profile_single(sm, sp, s[i0], un, dtdL, Ip_passive, Im_passive); // Plus state on face i