From 12666acc6cacfb3c59ae6a7d32b70babac82c933 Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Thu, 26 Sep 2024 18:02:20 -0400 Subject: [PATCH] update primitive variable reconstruction with spherical 2d geometry. (#2964) this adds both the geometry source terms to the interface tracing and makes sure that the d/dtheta terms have the r weighting --- Source/hydro/reconstruction.H | 65 ++++++++++++++++++++----------- Source/hydro/trace_ppm.cpp | 72 ++++++++++++++++++++--------------- 2 files changed, 85 insertions(+), 52 deletions(-) diff --git a/Source/hydro/reconstruction.H b/Source/hydro/reconstruction.H index 9a41710663..8f791a1ad5 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. - 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); + // 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,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 + + // 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. - // note: this is assumed to be working only in the x-direction + // 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,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); + 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,ncomp); + 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,27 @@ 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. - // note: this is assumed to be working only in the x-direction + // 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,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..a4fa02e6e5 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -55,10 +55,10 @@ 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]; #ifndef AMREX_USE_GPU auto lo = bx.loVect3d(); @@ -83,8 +83,11 @@ Castro::trace_ppm(const Box& bx, for (int n = 0; n < NQSRC; n++) { do_source_trace[n] = 0; - // geometric source terms in r-direction need tracing - if (coord > 0 && idir == 0 && (n == QRHO || n == QPRES || n == QREINT)) { + // geometric source terms in r-direction for nonCartesian coordinate + // or theta-direction in spherical coordinate need tracing + + if ((coord == 2 || (idir == 0 && coord == 1)) + && (n == QRHO || n == QPRES || n == QREINT)) { do_source_trace[n] = 1; continue; } @@ -154,6 +157,15 @@ Castro::trace_ppm(const Box& bx, amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + Real dtdL = dt / dx[idir]; + Real dL = dx[idir]; + + // 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 +208,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 +219,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 +236,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 +250,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 +261,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 +276,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 +291,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 +301,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 +325,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 +337,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 +347,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 +363,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 +373,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 +395,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 +410,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 +428,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