diff --git a/Source/hydro/reconstruction.H b/Source/hydro/reconstruction.H index 8f791a1ad5..86241e8e73 100644 --- a/Source/hydro/reconstruction.H +++ b/Source/hydro/reconstruction.H @@ -110,6 +110,43 @@ add_geometric_rho_source(amrex::Array4 const& q_arr, 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 +void +delete_UN_geometric_pres_source(amrex::Array4 const& q_arr, + amrex::Array4 const& dloga, + const int i, const int j, const int k, + amrex::Real* s) { + + using namespace reconstruction; + + // Depending on the direction of reconstruction, delete corresponding geometric + // pressure source in the velocity equations. + + // If idir == 0, i.e. r-direction, we should delete the geometric pressure source + // that were obtained by grouping pressure into x (r) momentum flux. + + // If idir == 1, i.e. theta-direction, delete the geometric pressure source in V + // that were obtained by grouping pressure into y (theta) momentum flux. + + // For idir == 0, i.e. r-direction: + // this takes the form: -alpha p / 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: - p 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,QPRES); + s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QPRES); + s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QPRES); + s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QPRES); + s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QPRES); +} + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void add_geometric_rhoe_source(amrex::Array4 const& q_arr, diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index a824fbb56f..6d2150c43f 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -86,12 +86,19 @@ Castro::trace_ppm(const Box& bx, // geometric source terms in r-direction for nonCartesian coordinate // or theta-direction in spherical coordinate need tracing - if ((coord == 2 || (idir == 0 && coord == 1)) + if ((idir == 0 && coord == 1) && (n == QRHO || n == QPRES || n == QREINT)) { do_source_trace[n] = 1; continue; } + if (coord == 2 && + (n == QRHO || n == QPRES || n == QREINT + n == QU || n == QV)) { + do_source_trace[n] = 1; + continue; + } + for (int k = lo[2]-2*dg2; k <= hi[2]+2*dg2; k++) { for (int j = lo[1]-2*dg1; j <= hi[1]+2*dg1; j++) { for (int i = lo[0]-2; i <= hi[0]+2; i++) { @@ -337,11 +344,20 @@ Castro::trace_ppm(const Box& bx, #ifndef AMREX_USE_GPU do_trace = do_source_trace[QUN]; #else - do_trace = check_trace_source(srcQ, idir, i, j, k, QUN); + if (coord == 2 || (idir == 0 && coord = 1)) { + do_trace = 1; + } else { + do_trace = check_trace_source(srcQ, idir, i, j, k, QUN); + } #endif if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUN, s); +#if AMREX_SPACEDIM <= 2 + if (coord == 2 || (idir == 0 && coord == 1)) { + delete_UN_geometric_pres_source(q_arr, dloga, i, j, k, s); + } +#endif ppm_reconstruct(s, flat, sm, sp); 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); diff --git a/Source/sources/Castro_geom.cpp b/Source/sources/Castro_geom.cpp index 5322539870..9615846585 100644 --- a/Source/sources/Castro_geom.cpp +++ b/Source/sources/Castro_geom.cpp @@ -159,41 +159,40 @@ Castro::fill_RZ_geom_source (Real time, Real dt, MultiFab& cons_state, MultiFab& const Box& bx = mfi.tilebox(); - Array4 const U_arr = cons_state.array(mfi); + // dlogaX = 1 / r for Cylindrical + + Array4 const dlogaX = (dLogArea[0]).array(mfi); + Array4 const U_arr = cons_state.array(mfi); Array4 const src = geom_src.array(mfi); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real rhoinv = 1.0_rt / U_arr(i,j,k,QRHO); - - // radius for non-Cartesian - Real rinv = 1.0_rt / (prob_lo[0] + (static_cast(i) + 0.5_rt)*dx[0]); + Real rhoinv = 1.0_rt / U_arr(i,j,k,URHO); // radial momentum: F = rho v_phi**2 / r - src(i,j,k,UMX) = U_arr(i,j,k,UMZ) * U_arr(i,j,k,UMZ) * rhoinv * rinv; + src(i,j,k,UMX) = U_arr(i,j,k,UMZ) * U_arr(i,j,k,UMZ) * dlogaX(i,j,k) * rhoinv ; // azimuthal momentum: F = - rho v_r v_phi / r - src(i,j,k,UMZ) = - U_arr(i,j,k,UMX) * U_arr(i,j,k,UMZ) * rhoinv * rinv; + src(i,j,k,UMZ) = - U_arr(i,j,k,UMX) * U_arr(i,j,k,UMZ) * dlogaX(i,j,k) * rhoinv; // We absorb pressure gradient into div Flux term, // which causes an extra geometric pressure term. Address that here. // Find pressure first: - // eos_rep_t eos_state; - // eos_state.rho = U_arr(i,j,k,URHO); - // eos_state.e = U_arr(i,j,k,UEINT) * rhoinv; - // for (int n = 0; n < NumSpec; n++) { - // eos_state.xn[n] = U_arr(i,j,k,UFS+n) * rhoinv; - // } - // eos(eos_input_re, eos_state); - - // // radial geometric pres term: F = p/r - // src(i,j,k,UMX) = eos_state.p * rinv; + eos_rep_t eos_state; + eos_state.rho = U_arr(i,j,k,URHO); + eos_state.e = U_arr(i,j,k,UEINT) * rhoinv; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = U_arr(i,j,k,UFS+n) * rhoinv; + } + eos(eos_input_re, eos_state); + // radial geometric pres term: F = p/r + src(i,j,k,UMX) = eos_state.p * dlogaX(i,j,k); }); } } @@ -220,49 +219,51 @@ Castro::fill_RTheta_geom_source (Real time, Real dt, MultiFab& cons_state, Multi const Box& bx = mfi.tilebox(); - Array4 const U_arr = cons_state.array(mfi); + // dlogaX = 2 / r for Spherical + // dlogaY = cot(theta) / r for Spherical + + Array4 const dlogaX = (dLogArea[0]).array(mfi); + Array4 const dlogaY = (dLogArea[1]).array(mfi); + Array4 const U_arr = cons_state.array(mfi); Array4 const src = geom_src.array(mfi); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real rhoinv = 1.0_rt / U_arr(i,j,k,QRHO); - - // Cell-centered Spherical Radius and Theta - Real rinv = 1.0_rt / (prob_lo[0] + (static_cast(i) + 0.5_rt)*dx[0]); - Real cotTheta = cot(prob_lo[1] + (static_cast(j) + 0.5_rt)*dx[1]); + Real rhoinv = 1.0_rt / U_arr(i,j,k,URHO); + Real rinv = 0.5_rt * dlogaX(i,j,k); // radial momentum: F = rho (v_theta**2 + v_phi**2) / r src(i,j,k,UMX) = (U_arr(i,j,k,UMY) * U_arr(i,j,k,UMY) + U_arr(i,j,k,UMZ) * U_arr(i,j,k,UMZ)) * rhoinv * rinv; // Theta momentum F = rho v_phi**2 cot(theta) / r - rho v_r v_theta / r - src(i,j,k,UMY) = (U_arr(i,j,k,UMZ) * U_arr(i,j,k,UMZ) * cotTheta - - U_arr(i,j,k,UMX) * U_arr(i,j,k,UMY)) * rhoinv * rinv; + src(i,j,k,UMY) = (U_arr(i,j,k,UMZ) * U_arr(i,j,k,UMZ) * dlogaY(i,j,k) + + U_arr(i,j,k,UMX) * U_arr(i,j,k,UMY) * rinv) * rhoinv; // Phi momentum: F = - rho v_r v_phi / r - rho v_theta v_phi cot(theta) / r - src(i,j,k,UMZ) = (- U_arr(i,j,k,UMY) * U_arr(i,j,k,UMZ) * cotTheta - - U_arr(i,j,k,UMX) * U_arr(i,j,k,UMZ)) * rhoinv * rinv; + src(i,j,k,UMZ) = (-U_arr(i,j,k,UMY) * U_arr(i,j,k,UMZ) * dlogaY(i,j,k) - + U_arr(i,j,k,UMX) * U_arr(i,j,k,UMZ) * rinv) * rhoinv; // We absorb pressure gradient into div Flux term, // which causes an extra geometric pressure term. Address that here. // Find pressure first: - // eos_rep_t eos_state; - // eos_state.rho = U_arr(i,j,k,URHO); - // eos_state.e = U_arr(i,j,k,UEINT) * rhoinv; - // for (int n = 0; n < NumSpec; n++) { - // eos_state.xn[n] = U_arr(i,j,k,UFS+n) * rhoinv; - // } - // eos(eos_input_re, eos_state); + eos_rep_t eos_state; + eos_state.rho = U_arr(i,j,k,URHO); + eos_state.e = U_arr(i,j,k,UEINT) * rhoinv; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = U_arr(i,j,k,UFS+n) * rhoinv; + } + eos(eos_input_re, eos_state); - // // radial geometric pres term: F = 2 p / r - // src(i,j,k,UMX) = 2.0_rt * eos_state.p * rinv; + // radial geometric pres term: F = 2 p / r + src(i,j,k,UMX) += eos_state.p * dlogaX(i,j,k); - // // Theta geometric pres term: F = cot(theta) p / r - // src(i,j,k,UMY) = cotTheta * eos_state.p * rinv; + // Theta geometric pres term: F = cot(theta) p / r + src(i,j,k,UMY) += eos_state.p * dlogaY(i,j,k); }); } }