From 22d6158e94282a740459ba4c0829e54d1de759da Mon Sep 17 00:00:00 2001 From: Max Katz Date: Tue, 10 Oct 2023 09:20:39 -0400 Subject: [PATCH 1/4] Convert correct_dterm to C++ (#2620) This uses inline metric calculation to support GPU offloading. Note that the vector s was unused and was previously only in the code to enable the call to sphe to fill in re for spherical geometries, which is now handled through edge_center_metric if we're in a spherical geometry. We also do a loop over dimensions to make the code a little more understandable. --- Source/radiation/CastroRad_1d.f90 | 18 ------- Source/radiation/CastroRad_2d.f90 | 30 ----------- Source/radiation/CastroRad_3d.f90 | 18 ------- Source/radiation/Make.package | 1 - Source/radiation/RAD_F.H | 11 ---- Source/radiation/RadSolve.cpp | 86 ++++++++++++++----------------- 6 files changed, 38 insertions(+), 126 deletions(-) delete mode 100644 Source/radiation/CastroRad_1d.f90 delete mode 100644 Source/radiation/CastroRad_2d.f90 delete mode 100644 Source/radiation/CastroRad_3d.f90 diff --git a/Source/radiation/CastroRad_1d.f90 b/Source/radiation/CastroRad_1d.f90 deleted file mode 100644 index 94732d9c25..0000000000 --- a/Source/radiation/CastroRad_1d.f90 +++ /dev/null @@ -1,18 +0,0 @@ -! no tiling -subroutine ca_correct_dterm(dfx, dfx_l1, dfx_h1, & - re, rc) bind(C, name="ca_correct_dterm") - - use amrex_fort_module, only : rt => amrex_real - implicit none - - integer, intent(in) :: dfx_l1, dfx_h1 - real(rt) , intent(inout) :: dfx(dfx_l1:dfx_h1) - real(rt) , intent(in) :: re(dfx_l1:dfx_h1), rc(1) - - integer :: i - - do i=dfx_l1, dfx_h1 - dfx(i) = dfx(i) / (re(i) + 1.e-50_rt) - end do - -end subroutine ca_correct_dterm diff --git a/Source/radiation/CastroRad_2d.f90 b/Source/radiation/CastroRad_2d.f90 deleted file mode 100644 index 892aadac3b..0000000000 --- a/Source/radiation/CastroRad_2d.f90 +++ /dev/null @@ -1,30 +0,0 @@ -! no tiling -subroutine ca_correct_dterm( & - dfx, dfx_l1, dfx_l2, dfx_h1, dfx_h2, & - dfy, dfy_l1, dfy_l2, dfy_h1, dfy_h2, & - re, rc) bind(C, name="ca_correct_dterm") - - use amrex_fort_module, only : rt => amrex_real - implicit none - - integer, intent(in) :: dfx_l1, dfx_l2, dfx_h1, dfx_h2 - integer, intent(in) :: dfy_l1, dfy_l2, dfy_h1, dfy_h2 - real(rt) , intent(inout) :: dfx(dfx_l1:dfx_h1,dfx_l2:dfx_h2) - real(rt) , intent(inout) :: dfy(dfy_l1:dfy_h1,dfy_l2:dfy_h2) - real(rt) , intent(in) :: re(dfx_l1:dfx_h1), rc(dfy_l1:dfy_h1) - - integer :: i, j - - do j=dfx_l2, dfx_h2 - do i=dfx_l1, dfx_h1 - dfx(i,j) = dfx(i,j) / (re(i) + 1.e-50_rt) - end do - end do - - do j=dfy_l2, dfy_h2 - do i=dfy_l1, dfy_h1 - dfy(i,j) = dfy(i,j) / rc(i) - end do - end do - -end subroutine ca_correct_dterm diff --git a/Source/radiation/CastroRad_3d.f90 b/Source/radiation/CastroRad_3d.f90 deleted file mode 100644 index c21d2945ac..0000000000 --- a/Source/radiation/CastroRad_3d.f90 +++ /dev/null @@ -1,18 +0,0 @@ -subroutine ca_correct_dterm( & - dfx, dfx_l1, dfx_l2, dfx_l3, dfx_h1, dfx_h2, dfx_h3, & - dfy, dfy_l1, dfy_l2, dfy_l3, dfy_h1, dfy_h2, dfy_h3, & - dfz, dfz_l1, dfz_l2, dfz_l3, dfz_h1, dfz_h2, dfz_h3, & - re, rc) bind(C, name="ca_correct_dterm") - - use amrex_fort_module, only : rt => amrex_real - implicit none - - integer, intent(in) :: dfx_l1, dfx_l2, dfx_l3, dfx_h1, dfx_h2, dfx_h3 - integer, intent(in) :: dfy_l1, dfy_l2, dfy_l3, dfy_h1, dfy_h2, dfy_h3 - integer, intent(in) :: dfz_l1, dfz_l2, dfz_l3, dfz_h1, dfz_h2, dfz_h3 - real(rt) , intent(inout) :: dfx(dfx_l1:dfx_h1,dfx_l2:dfx_h2,dfx_l3:dfx_h3) - real(rt) , intent(inout) :: dfy(dfy_l1:dfy_h1,dfy_l2:dfy_h2,dfy_l3:dfy_h3) - real(rt) , intent(inout) :: dfz(dfz_l1:dfz_h1,dfz_l2:dfz_h2,dfz_l3:dfz_h3) - real(rt) , intent(in) :: re(1), rc(1) - -end subroutine ca_correct_dterm diff --git a/Source/radiation/Make.package b/Source/radiation/Make.package index 977a82fee5..1def860164 100644 --- a/Source/radiation/Make.package +++ b/Source/radiation/Make.package @@ -36,7 +36,6 @@ ca_F90EXE_sources += RAD_$(DIM)D.F90 ca_F90EXE_sources += HABEC_$(DIM)D.F90 CEXE_sources += trace_ppm_rad.cpp -ca_f90EXE_sources += CastroRad_$(DIM)d.f90 ca_F90EXE_sources += rad_params.F90 ca_F90EXE_sources += Rad_nd.F90 diff --git a/Source/radiation/RAD_F.H b/Source/radiation/RAD_F.H index 079dedba62..3b6e374ba0 100644 --- a/Source/radiation/RAD_F.H +++ b/Source/radiation/RAD_F.H @@ -53,17 +53,6 @@ BL_FORT_PROC_DECL(CA_INITGROUPS3,ca_initgroups3) BL_FORT_PROC_DECL(CA_COMPUTE_KAPKAP, ca_compute_kapkap) (BL_FORT_FAB_ARG(kapkap), const BL_FORT_FAB_ARG(kap_r)); -#ifdef __cplusplus -extern "C" { -#endif - void ca_correct_dterm - (D_DECL(BL_FORT_FAB_ARG(dfx), - BL_FORT_FAB_ARG(dfy), - BL_FORT_FAB_ARG(dfz)), - const amrex::Real* re, const amrex::Real* rc); -#ifdef __cplusplus -} -#endif #ifdef __cplusplus extern "C" { diff --git a/Source/radiation/RadSolve.cpp b/Source/radiation/RadSolve.cpp index e671cf66f1..12c8fd7b5f 100644 --- a/Source/radiation/RadSolve.cpp +++ b/Source/radiation/RadSolve.cpp @@ -813,6 +813,7 @@ void RadSolve::levelDterm(int level, MultiFab& Dterm, MultiFab& Er, int igroup) const BoxArray& grids = parent->boxArray(level); const DistributionMapping& dmap = parent->DistributionMap(level); const Geometry& geom = parent->Geom(level); + const GeometryData& geomdata = geom.data(); auto dx = parent->Geom(level).CellSizeArray(); const Castro *castro = dynamic_cast(&parent->getLevel(level)); @@ -863,75 +864,64 @@ void RadSolve::levelDterm(int level, MultiFab& Dterm, MultiFab& Er, int igroup) // Correct D terms at physical and coarse-fine boundaries. hem->boundaryDterm(level, &Dterm_face[0], Er, igroup); + // Correct for metric terms (only has an effect in non-Cartesian geometries). + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { #ifdef _OPENMP #pragma omp parallel #endif - { - Vector rc, re, s; - - if (geom.IsSPHERICAL() || geom.IsRZ()) { - for (MFIter fi(Dterm_face[0]); fi.isValid(); ++fi) { // omp over boxes - int i = fi.index(); - const Box ® = grids[i]; + for (MFIter mfi(Dterm_face[dir], true); mfi.isValid(); ++mfi) { + const Box& box = mfi.tilebox(); + Array4 const d = Dterm_face[dir][mfi].array(); - parent->Geom(level).GetEdgeLoc(re, reg, 0); - parent->Geom(level).GetCellLoc(rc, reg, 0); - - if (geom.IsSPHERICAL()) { - parent->Geom(level).GetCellLoc(s, reg, 0); + amrex::ParallelFor(box, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (dir == 0 && (geomdata.Coord() == CoordSys::SPHERICAL || geomdata.Coord() == CoordSys::RZ)) { + Real r, s; + edge_center_metric(i, j, k, dir, geomdata, r, s); - const Box &dbox = Dterm_face[0][fi].box(); + d(i,j,k) = d(i,j,k) / (r + 1.0e-50_rt); + } + else if (dir == 1 && geomdata.Coord() == CoordSys::RZ) { + Real r, s; + cell_center_metric(i, j, k, geomdata, r, s); - for (int i = dbox.loVect()[0]; i <= dbox.hiVect()[0]; ++i) { - re[i] *= re[i]; - } -#if AMREX_SPACEDIM >= 2 - Real h2 = 0.5e0_rt * dx[1]; - Real d2 = 1.e0_rt / dx[1]; - for (int j = dbox.loVect()[1]; j <= dbox.hiVect()[1]; ++j) { - s[j] = d2 * (std::cos(s[j] - h2) - std::cos(s[j] + h2)); - } -#endif + d(i,j,k) = d(i,j,k) / r; } + }); + } + } - ca_correct_dterm(D_DECL(BL_TO_FORTRAN(Dterm_face[0][fi]), - BL_TO_FORTRAN(Dterm_face[1][fi]), - BL_TO_FORTRAN(Dterm_face[2][fi])), - re.dataPtr(), rc.dataPtr()); - } #ifdef _OPENMP -#pragma omp barrier +#pragma omp parallel #endif - } - - for (MFIter fi(Dterm,true); fi.isValid(); ++fi) { - const Box& bx = fi.tilebox(); + for (MFIter fi(Dterm,true); fi.isValid(); ++fi) { + const Box& bx = fi.tilebox(); - auto Dx = Dterm_face[0][fi].array(); + auto Dx = Dterm_face[0][fi].array(); #if AMREX_SPACEDIM >= 2 - auto Dy = Dterm_face[1][fi].array(); + auto Dy = Dterm_face[1][fi].array(); #endif #if AMREX_SPACEDIM == 3 - auto Dz = Dterm_face[2][fi].array(); + auto Dz = Dterm_face[2][fi].array(); #endif - auto D = Dterm[fi].array(); + auto D = Dterm[fi].array(); - amrex::ParallelFor(bx, - [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) - { + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { #if AMREX_SPACEDIM == 1 - D(i,j,k) = (Dx(i,j,k) + Dx(i+1,j,k)) * 0.5_rt; + D(i,j,k) = (Dx(i,j,k) + Dx(i+1,j,k)) * 0.5_rt; #elif AMREX_SPACEDIM == 2 - D(i,j,k) = (Dx(i,j,k) + Dx(i+1,j,k) + - Dy(i,j,k) + Dy(i,j+1,k)) * 0.25_rt; + D(i,j,k) = (Dx(i,j,k) + Dx(i+1,j,k) + + Dy(i,j,k) + Dy(i,j+1,k)) * 0.25_rt; #else - D(i,j,k) = (Dx(i,j,k) + Dx(i+1,j,k) + - Dy(i,j,k) + Dy(i,j+1,k) + - Dz(i,j,k) + Dz(i,j,k+1)) * (1.0_rt / 6.0_rt); + D(i,j,k) = (Dx(i,j,k) + Dx(i+1,j,k) + + Dy(i,j,k) + Dy(i,j+1,k) + + Dz(i,j,k) + Dz(i,j,k+1)) * (1.0_rt / 6.0_rt); #endif - }); - } + }); } } From 9b81621ad4668811be60d992bcf0db3aee699c2e Mon Sep 17 00:00:00 2001 From: Max Katz Date: Tue, 10 Oct 2023 12:31:53 -0400 Subject: [PATCH 2/4] Convert hdterm and hdterm3 to C++ (#2625) --- Source/radiation/HABEC.H | 216 +++++++++++++++++ Source/radiation/HABEC_1D.F90 | 133 ----------- Source/radiation/HABEC_2D.F90 | 206 ----------------- Source/radiation/HABEC_3D.F90 | 308 ------------------------- Source/radiation/HABEC_F.H | 37 --- Source/radiation/HypreABec.cpp | 1 - Source/radiation/HypreExtMultiABec.cpp | 38 +-- Source/radiation/HypreMultiABec.cpp | 1 - Source/radiation/Make.package | 2 - Source/radiation/RadSolve.cpp | 1 - 10 files changed, 235 insertions(+), 708 deletions(-) delete mode 100644 Source/radiation/HABEC_1D.F90 delete mode 100644 Source/radiation/HABEC_2D.F90 delete mode 100644 Source/radiation/HABEC_3D.F90 delete mode 100644 Source/radiation/HABEC_F.H diff --git a/Source/radiation/HABEC.H b/Source/radiation/HABEC.H index 2c9df6b2e3..3039e759d3 100644 --- a/Source/radiation/HABEC.H +++ b/Source/radiation/HABEC.H @@ -367,6 +367,222 @@ namespace HABEC } }); } + + AMREX_INLINE + void hdterm (Array4 const dterm, + Array4 const er, + const Box& reg, + int cdir, int bct, Real bcl, + Array4 const bcval, + Array4 const mask, + Array4 const d, + const Real* dx) + { + Real h; + + // Index shift for whichever edge we're checking. + // Negative means we're looking at the lo edge, + // positive means we're looking at the hi edge. + // The first group is for cell centers, the second + // group is for cell edges. + + int icp = 0; + int jcp = 0; + int kcp = 0; + + int iep = 0; + int jep = 0; + int kep = 0; + +#if AMREX_SPACEDIM == 1 + if (cdir == 0) { + h = dx[0]; + icp = -1; + } + else if (cdir == 1) { + h = dx[0]; + icp = 1; + iep = 1; + } +#elif AMREX_SPACEDIM == 2 + if (cdir == 0) { + h = dx[0]; + icp = -1; + } + else if (cdir == 2) { + h = dx[0]; + icp = 1; + iep = 1; + } + else if (cdir == 1) { + h = dx[1]; + jcp = -1; + } + else if (cdir == 3) { + h = dx[1]; + jcp = 1; + jep = 1; + } +#else + if (cdir == 0) { + h = dx[0]; + icp = -1; + } + else if (cdir == 3) { + h = dx[0]; + icp = 1; + iep = 1; + } + else if (cdir == 1) { + h = dx[1]; + jcp = -1; + } + else if (cdir == 4) { + h = dx[1]; + jcp = 1; + jep = 1; + } + else if (cdir == 2) { + h = dx[2]; + kcp = -1; + } + else if (cdir == 5) { + h = dx[2]; + kcp = 1; + kep = 1; + } +#endif + + amrex::LoopOnCpu(reg, [=] (int i, int j, int k) noexcept + { + if (mask.contains(i+icp,j+jcp,k+kcp)) { + if (mask(i+icp,j+jcp,k+kcp) > 0) { + if (bct == LO_DIRICHLET) { + Real d_sign = 1.0_rt; + if (iep != 0 || jep != 0 || kep != 0) { + // right edge + d_sign = -1.0_rt; + } + dterm(i+iep,j+jep,k+kep) = d(i+iep,j+jep,k+kep) * d_sign * (er(i,j,k) - bcval(i+icp,j+jcp,k+kcp)) / (0.5_rt * h + bcl); + } + } + } + }); + } + + AMREX_INLINE + void hdterm3 (Array4 const dterm, + Array4 const er, + const Box& reg, + int cdir, int bctype, + Array4 const tf, + Real bcl, + Array4 const bcval, + Array4 const mask, + Array4 const d, + const Real* const dx) + { + Real h; + + // Index shift for whichever edge we're checking. + // Negative means we're looking at the lo edge, + // positive means we're looking at the hi edge. + // The first group is for cell centers, the second + // group is for cell edges. + + int icp = 0; + int jcp = 0; + int kcp = 0; + + int iep = 0; + int jep = 0; + int kep = 0; + +#if AMREX_SPACEDIM == 1 + if (cdir == 0) { + h = dx[0]; + icp = -1; + } + else if (cdir == 1) { + h = dx[0]; + icp = 1; + iep = 1; + } +#elif AMREX_SPACEDIM == 2 + if (cdir == 0) { + h = dx[0]; + icp = -1; + } + else if (cdir == 2) { + h = dx[0]; + icp = 1; + iep = 1; + } + else if (cdir == 1) { + h = dx[1]; + jcp = -1; + } + else if (cdir == 3) { + h = dx[1]; + jcp = 1; + jep = 1; + } +#else + if (cdir == 0) { + h = dx[0]; + icp = -1; + } + else if (cdir == 3) { + h = dx[0]; + icp = 1; + iep = 1; + } + else if (cdir == 1) { + h = dx[1]; + jcp = -1; + } + else if (cdir == 4) { + h = dx[1]; + jcp = 1; + jep = 1; + } + else if (cdir == 2) { + h = dx[2]; + kcp = -1; + } + else if (cdir == 5) { + h = dx[2]; + kcp = 1; + kep = 1; + } +#endif + + amrex::LoopOnCpu(reg, [=] (int i, int j, int k) noexcept + { + if (mask.contains(i+icp,j+jcp,k+kcp)) { + if (mask(i+icp,j+jcp,k+kcp) > 0) { + int bct; + if (bctype == -1) { + bct = tf(i+icp,j+jcp,k+kcp); + } + else { + bct = bctype; + } + if (bct == LO_DIRICHLET) { + Real d_sign = 1.0_rt; + if (iep != 0 || jep != 0 || kep != 0) { + // right edge + d_sign = -1.0_rt; + } + dterm(i+iep,j+jep,k+kep) = d(i+iep,j+jep,k+kep) * d_sign * (er(i,j,k) - bcval(i+icp,j+jcp,k+kcp)) / (0.5_rt * h + bcl); + } + else if (bct == LO_NEUMANN && bcval(i+icp,j+jcp,k+kcp) == 0.0_rt) { + dterm(i+iep,j+jep,k+kep) = 0.0_rt; + } + } + } + }); + } } // namespace HABEC #endif diff --git a/Source/radiation/HABEC_1D.F90 b/Source/radiation/HABEC_1D.F90 deleted file mode 100644 index 544297195e..0000000000 --- a/Source/radiation/HABEC_1D.F90 +++ /dev/null @@ -1,133 +0,0 @@ -#include -#include - -module habec_module - - ! habec is Hypre abec, where abec is the form of the linear equation - ! we are solving: - ! - ! alpha*phi - div(beta*grad phi) + div(\vec{c}*phi) - - use amrex_fort_module, only : rt => amrex_real - implicit none - -contains - -subroutine hdterm(dterm, & - DIMS(dtbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bct, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - d, DIMS(dbox), & - dx) bind(C, name="hdterm") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(dtbox) - integer :: DIMDEC(ebox) - integer :: DIMDEC(reg) - integer :: DIMDEC(bcv) - integer :: DIMDEC(msk) - integer :: DIMDEC(dbox) - integer :: cdir, bct - real(rt) :: bcl, dx(1) - real(rt) :: dterm(DIMV(dtbox)) - real(rt) :: er(DIMV(ebox)) - real(rt) :: bcval(DIMV(bcv)) - integer :: mask(DIMV(msk)) - real(rt) :: d(DIMV(dbox)) - real(rt) :: h - integer :: i - h = dx(1) - if (bct == LO_DIRICHLET) then - if (cdir == 0) then - ! Left face of grid - i = reg_l1 - if (mask(i-1) > 0) then - dterm(i) = d(i)*(er(i) - bcval(i-1))/(0.5e0_rt*h+bcl) - endif - else if (cdir == 1) then - ! Right face of grid - i = reg_h1 - if (mask(i+1) > 0) then - dterm(i+1) = d(i+1)*(bcval(i+1)-er(i))/(0.5e0_rt*h+bcl) - endif - else - print *, "hdterm: impossible face orientation" - endif - else - print *, "hdterm: unsupported boundary type" - stop - endif -end subroutine hdterm - -subroutine hdterm3(dterm, & - DIMS(dtbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bctype, tf, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - d, DIMS(dbox), & - dx) bind(C, name="hdterm3") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(dtbox) - integer :: DIMDEC(ebox) - integer :: DIMDEC(reg) - integer :: DIMDEC(bcv) - integer :: DIMDEC(msk) - integer :: DIMDEC(dbox) - integer :: cdir, bctype, tf(DIMV(bcv)) - real(rt) :: bcl, dx(1) - real(rt) :: dterm(DIMV(dtbox)) - real(rt) :: er(DIMV(ebox)) - real(rt) :: bcval(DIMV(bcv)) - integer :: mask(DIMV(msk)) - real(rt) :: d(DIMV(dbox)) - real(rt) :: h - integer :: i, bct - h = dx(1) - if (cdir == 0) then - ! Left face of grid - i = reg_l1 - if (mask(i-1) > 0) then - if (bctype == -1) then - bct = tf(i-1) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - dterm(i) = d(i)*(er(i) - bcval(i-1))/(0.5e0_rt*h+bcl) - else if (bct == LO_NEUMANN .AND. bcval(i-1) == 0.e0_rt) then - dterm(i) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - endif - endif - else if (cdir == 1) then - ! Right face of grid - i = reg_h1 - if (mask(i+1) > 0) then - if (bctype == -1) then - bct = tf(i+1) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - dterm(i+1) = d(i+1)*(bcval(i+1)-er(i))/(0.5e0_rt*h+bcl) - else if (bct == LO_NEUMANN .AND. bcval(i+1) == 0.e0_rt) then - dterm(i+1) = 0.e0_rt - else - print *, "hbterm3: unsupported boundary type" - stop - endif - endif - else - print *, "hdterm3: impossible face orientation" - endif -end subroutine hdterm3 - -end module habec_module diff --git a/Source/radiation/HABEC_2D.F90 b/Source/radiation/HABEC_2D.F90 deleted file mode 100644 index 837e29d219..0000000000 --- a/Source/radiation/HABEC_2D.F90 +++ /dev/null @@ -1,206 +0,0 @@ -#include -#include - - -module habec_module - - ! habec is Hypre abec, where abec is the form of the linear equation - ! we are solving: - ! - ! alpha*phi - div(beta*grad phi) + div(\vec{c}*phi) - - use amrex_fort_module, only : rt => amrex_real - implicit none - -contains - -subroutine hdterm(dterm, & - DIMS(dtbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bct, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - d, DIMS(dbox), & - dx) bind(C, name="hdterm") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(dtbox) - integer :: DIMDEC(ebox) - integer :: DIMDEC(reg) - integer :: DIMDEC(bcv) - integer :: DIMDEC(msk) - integer :: DIMDEC(dbox) - integer :: cdir, bct - real(rt) :: bcl, dx(2) - real(rt) :: dterm(DIMV(dtbox)) - real(rt) :: er(DIMV(ebox)) - real(rt) :: bcval(DIMV(bcv)) - integer :: mask(DIMV(msk)) - real(rt) :: d(DIMV(dbox)) - real(rt) :: h - integer :: i, j - if (cdir == 0 .OR. cdir == 2) then - h = dx(1) - else - h = dx(2) - endif - if (bct == LO_DIRICHLET) then - if (cdir == 0) then - ! Left face of grid - i = reg_l1 - do j = reg_l2, reg_h2 - if (mask(i-1,j) > 0) then - dterm(i,j) = d(i,j)*(er(i,j)-bcval(i-1,j))/(0.5e0_rt*h+bcl) - endif - enddo - else if (cdir == 2) then - ! Right face of grid - i = reg_h1 - do j = reg_l2, reg_h2 - if (mask(i+1,j) > 0) then - dterm(i+1,j) = d(i+1,j)*(bcval(i+1,j)-er(i,j))/(0.5e0_rt*h+bcl) - endif - enddo - else if (cdir == 1) then - ! Bottom face of grid - j = reg_l2 - do i = reg_l1, reg_h1 - if (mask(i,j-1) > 0) then - dterm(i,j) = d(i,j)*(er(i,j)-bcval(i,j-1))/(0.5e0_rt*h+bcl) - endif - enddo - else if (cdir == 3) then - ! Top face of grid - j = reg_h2 - do i = reg_l1, reg_h1 - if (mask(i,j+1) > 0) then - dterm(i,j+1) = d(i,j+1)*(bcval(i,j+1)-er(i,j))/(0.5e0_rt*h+bcl) - endif - enddo - else - print *, "hdterm: impossible face orientation" - endif - else - print *, "hdterm: unsupported boundary type" - stop - endif -end subroutine hdterm - -subroutine hdterm3(dterm, & - DIMS(dtbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bctype, tf, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - d, DIMS(dbox), & - dx) bind(C, name="hdterm3") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(dtbox) - integer :: DIMDEC(ebox) - integer :: DIMDEC(reg) - integer :: DIMDEC(bcv) - integer :: DIMDEC(msk) - integer :: DIMDEC(dbox) - integer :: cdir, bctype, tf(DIMV(bcv)) - real(rt) :: bcl, dx(2) - real(rt) :: dterm(DIMV(dtbox)) - real(rt) :: er(DIMV(ebox)) - real(rt) :: bcval(DIMV(bcv)) - integer :: mask(DIMV(msk)) - real(rt) :: d(DIMV(dbox)) - real(rt) :: h - integer :: i, j, bct - if (cdir == 0 .OR. cdir == 2) then - h = dx(1) - else - h = dx(2) - endif - if (cdir == 0) then - ! Left face of grid - i = reg_l1 - do j = reg_l2, reg_h2 - if (mask(i-1,j) > 0) then - if (bctype == -1) then - bct = tf(i-1,j) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - dterm(i,j) = d(i,j)*(er(i,j)-bcval(i-1,j))/(0.5e0_rt*h+bcl) - else if (bct == LO_NEUMANN .AND. bcval(i-1,j) == 0.e0_rt) then - dterm(i,j) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - endif - endif - enddo - else if (cdir == 2) then - ! Right face of grid - i = reg_h1 - do j = reg_l2, reg_h2 - if (mask(i+1,j) > 0) then - if (bctype == -1) then - bct = tf(i+1,j) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - dterm(i+1,j) = d(i+1,j)*(bcval(i+1,j)-er(i,j))/(0.5e0_rt*h+bcl) - else if (bct == LO_NEUMANN .AND. bcval(i+1,j) == 0.e0_rt) then - dterm(i+1,j) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - endif - endif - enddo - else if (cdir == 1) then - ! Bottom face of grid - j = reg_l2 - do i = reg_l1, reg_h1 - if (mask(i,j-1) > 0) then - if (bctype == -1) then - bct = tf(i,j-1) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - dterm(i,j) = d(i,j)*(er(i,j)-bcval(i,j-1))/(0.5e0_rt*h+bcl) - else if (bct == LO_NEUMANN .AND. bcval(i,j-1) == 0.e0_rt) then - dterm(i,j) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - endif - endif - enddo - else if (cdir == 3) then - ! Top face of grid - j = reg_h2 - do i = reg_l1, reg_h1 - if (mask(i,j+1) > 0) then - if (bctype == -1) then - bct = tf(i,j+1) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - dterm(i,j+1) = d(i,j+1)*(bcval(i,j+1)-er(i,j))/(0.5e0_rt*h+bcl) - else if (bct == LO_NEUMANN .AND. bcval(i,j+1) == 0.e0_rt) then - dterm(i,j+1) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - endif - endif - enddo - else - print *, "hdterm3: impossible face orientation" - endif -end subroutine hdterm3 - -end module habec_module diff --git a/Source/radiation/HABEC_3D.F90 b/Source/radiation/HABEC_3D.F90 deleted file mode 100644 index d59edec61a..0000000000 --- a/Source/radiation/HABEC_3D.F90 +++ /dev/null @@ -1,308 +0,0 @@ -#include -#include - - -module habec_module - - ! habec is Hypre abec, where abec is the form of the linear equation - ! we are solving: - ! - ! alpha*phi - div(beta*grad phi) + div(\vec{c}*phi) - - use amrex_fort_module, only : rt => amrex_real - implicit none - -contains - -subroutine hdterm(dterm, & - DIMS(dtbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bct, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - d, DIMS(dbox), & - dx) bind(C, name="hdterm") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(dtbox) - integer :: DIMDEC(ebox) - integer :: DIMDEC(reg) - integer :: DIMDEC(bcv) - integer :: DIMDEC(msk) - integer :: DIMDEC(dbox) - integer :: cdir, bct - real(rt) :: bcl, dx(3) - real(rt) :: dterm(DIMV(dtbox)) - real(rt) :: er(DIMV(ebox)) - real(rt) :: bcval(DIMV(bcv)) - integer :: mask(DIMV(msk)) - real(rt) :: d(DIMV(dbox)) - real(rt) :: h, bfm, bfv - integer :: i, j, k - if (cdir == 0 .OR. cdir == 3) then - h = dx(1) - elseif (cdir == 1 .OR. cdir == 4) then - h = dx(2) - else - h = dx(3) - endif - if (bct == LO_DIRICHLET) then - if (cdir == 0) then - i = reg_l1 - do k = reg_l3, reg_h3 - do j = reg_l2, reg_h2 - if (mask(i-1,j,k) > 0) then - dterm(i,j,k) = d(i,j,k) * & - (er(i,j,k) - bcval(i-1,j,k)) & - / (0.5e0_rt*h + bcl) - endif - enddo - enddo - else if (cdir == 3) then - i = reg_h1 - do k = reg_l3, reg_h3 - do j = reg_l2, reg_h2 - if (mask(i+1,j,k) > 0) then - dterm(i+1,j,k) = d(i+1,j,k) * & - (bcval(i+1,j,k) - er(i,j,k)) & - / (0.5e0_rt*h + bcl) - endif - enddo - enddo - else if (cdir == 1) then - j = reg_l2 - do k = reg_l3, reg_h3 - do i = reg_l1, reg_h1 - if (mask(i,j-1,k) > 0) then - dterm(i,j,k) = d(i,j,k) * & - (er(i,j,k) - bcval(i,j-1,k)) & - / (0.5e0_rt*h + bcl) - endif - enddo - enddo - else if (cdir == 4) then - j = reg_h2 - do k = reg_l3, reg_h3 - do i = reg_l1, reg_h1 - if (mask(i,j+1,k) > 0) then - dterm(i,j+1,k) = d(i,j+1,k) * & - (bcval(i,j+1,k) - er(i,j,k)) & - / (0.5e0_rt*h + bcl) - endif - enddo - enddo - else if (cdir == 2) then - k = reg_l3 - do j = reg_l2, reg_h2 - do i = reg_l1, reg_h1 - if (mask(i,j,k-1) > 0) then - dterm(i,j,k) = d(i,j,k) * & - (er(i,j,k) - bcval(i,j,k-1)) & - / (0.5e0_rt*h + bcl) - endif - enddo - enddo - else if (cdir == 5) then - k = reg_h3 - do j = reg_l2, reg_h2 - do i = reg_l1, reg_h1 - if (mask(i,j,k+1) > 0) then - dterm(i,j,k+1) = d(i,j,k+1) * & - (bcval(i,j,k+1) - er(i,j,k)) & - / (0.5e0_rt*h + bcl) - endif - enddo - enddo - else - print *, "hdterm: impossible face orientation" - endif - else - print *, "hdterm: unsupported boundary type" - stop - endif -end subroutine hdterm - -subroutine hdterm3(dterm, & - DIMS(dtbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bctype, tf, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - d, DIMS(dbox), & - dx) bind(C, name="hdterm3") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(dtbox) - integer :: DIMDEC(ebox) - integer :: DIMDEC(reg) - integer :: DIMDEC(bcv) - integer :: DIMDEC(msk) - integer :: DIMDEC(dbox) - integer :: cdir, bctype, tf(DIMV(bcv)) - real(rt) :: bcl, dx(3) - real(rt) :: dterm(DIMV(dtbox)) - real(rt) :: er(DIMV(ebox)) - real(rt) :: bcval(DIMV(bcv)) - integer :: mask(DIMV(msk)) - real(rt) :: d(DIMV(dbox)) - real(rt) :: h, bfm, bfv - integer :: i, j, k, bct - if (cdir == 0 .OR. cdir == 3) then - h = dx(1) - elseif (cdir == 1 .OR. cdir == 4) then - h = dx(2) - else - h = dx(3) - endif - if (cdir == 0) then - i = reg_l1 - do k = reg_l3, reg_h3 - do j = reg_l2, reg_h2 - if (mask(i-1,j,k) > 0) then - if (bctype == -1) then - bct = tf(i-1,j,k) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - dterm(i,j,k) = d(i,j,k) * & - (er(i,j,k) - bcval(i-1,j,k)) & - / (0.5e0_rt*h + bcl) - else if (bct == LO_NEUMANN & - .AND. bcval(i-1,j,k) == 0.e0_rt) then - dterm(i,j,k) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - endif - endif - enddo - enddo - else if (cdir == 3) then - i = reg_h1 - do k = reg_l3, reg_h3 - do j = reg_l2, reg_h2 - if (mask(i+1,j,k) > 0) then - if (bctype == -1) then - bct = tf(i+1,j,k) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - dterm(i+1,j,k) = d(i+1,j,k) * & - (bcval(i+1,j,k) - er(i,j,k)) & - / (0.5e0_rt*h + bcl) - else if (bct == LO_NEUMANN & - .AND. bcval(i+1,j,k) == 0.e0_rt) then - dterm(i+1,j,k) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - endif - endif - enddo - enddo - else if (cdir == 1) then - j = reg_l2 - do k = reg_l3, reg_h3 - do i = reg_l1, reg_h1 - if (mask(i,j-1,k) > 0) then - if (bctype == -1) then - bct = tf(i,j-1,k) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - dterm(i,j,k) = d(i,j,k) * & - (er(i,j,k) - bcval(i,j-1,k)) & - / (0.5e0_rt*h + bcl) - else if (bct == LO_NEUMANN & - .AND. bcval(i,j-1,k) == 0.e0_rt) then - dterm(i,j,k) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - endif - endif - enddo - enddo - else if (cdir == 4) then - j = reg_h2 - do k = reg_l3, reg_h3 - do i = reg_l1, reg_h1 - if (mask(i,j+1,k) > 0) then - if (bctype == -1) then - bct = tf(i,j+1,k) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - dterm(i,j+1,k) = d(i,j+1,k) * & - (bcval(i,j+1,k) - er(i,j,k)) & - / (0.5e0_rt*h + bcl) - else if (bct == LO_NEUMANN & - .AND. bcval(i,j+1,k) == 0.e0_rt) then - dterm(i,j+1,k) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - endif - endif - enddo - enddo - else if (cdir == 2) then - k = reg_l3 - do j = reg_l2, reg_h2 - do i = reg_l1, reg_h1 - if (mask(i,j,k-1) > 0) then - if (bctype == -1) then - bct = tf(i,j,k-1) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - dterm(i,j,k) = d(i,j,k) * & - (er(i,j,k) - bcval(i,j,k-1)) & - / (0.5e0_rt*h + bcl) - else if (bct == LO_NEUMANN & - .AND. bcval(i,j,k-1) == 0.e0_rt) then - dterm(i,j,k) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - endif - endif - enddo - enddo - else if (cdir == 5) then - k = reg_h3 - do j = reg_l2, reg_h2 - do i = reg_l1, reg_h1 - if (mask(i,j,k+1) > 0) then - if (bctype == -1) then - bct = tf(i,j,k+1) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - dterm(i,j,k+1) = d(i,j,k+1) * & - (bcval(i,j,k+1) - er(i,j,k)) & - / (0.5e0_rt*h + bcl) - else if (bct == LO_NEUMANN & - .AND. bcval(i,j,k+1) == 0.e0_rt) then - dterm(i,j,k+1) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - endif - endif - enddo - enddo - else - print *, "hdterm3: impossible face orientation" - endif -end subroutine hdterm3 - -end module habec_module diff --git a/Source/radiation/HABEC_F.H b/Source/radiation/HABEC_F.H deleted file mode 100644 index 72a3f7be97..0000000000 --- a/Source/radiation/HABEC_F.H +++ /dev/null @@ -1,37 +0,0 @@ -#ifndef CASTRO_HABEC_F_H -#define CASTRO_HABEC_F_H - -#include -#include - -#ifdef __cplusplus -extern "C" { -#else -#define RadBoundCond int -#endif - - void hdterm(BL_FORT_FAB_ARG(dterm), - BL_FORT_FAB_ARG(soln), - ARLIM_P(reglo), ARLIM_P(reghi), - const int& cdir, const RadBoundCond& bct, - const amrex::Real& bcl, - const BL_FORT_FAB_ARG(bcval), - const BL_FORT_IFAB_ARG(mask), - BL_FORT_FAB_ARG(dcoefs), - const amrex::Real* dx); - - void hdterm3(BL_FORT_FAB_ARG(dterm), - BL_FORT_FAB_ARG(soln), - ARLIM_P(reglo), ARLIM_P(reghi), - const int& cdir, const int& bctype, const int* tf, - const amrex::Real& bcl, - const BL_FORT_FAB_ARG(bcval), - const BL_FORT_IFAB_ARG(mask), - BL_FORT_FAB_ARG(dcoefs), - const amrex::Real* dx); - -#ifdef __cplusplus -}; -#endif - -#endif diff --git a/Source/radiation/HypreABec.cpp b/Source/radiation/HypreABec.cpp index e3ac5a4931..638d4715bb 100644 --- a/Source/radiation/HypreABec.cpp +++ b/Source/radiation/HypreABec.cpp @@ -4,7 +4,6 @@ #include #include -#include #include #include diff --git a/Source/radiation/HypreExtMultiABec.cpp b/Source/radiation/HypreExtMultiABec.cpp index 9f958fae89..a65292751a 100644 --- a/Source/radiation/HypreExtMultiABec.cpp +++ b/Source/radiation/HypreExtMultiABec.cpp @@ -1,6 +1,6 @@ #include -#include +#include #include #include <_hypre_sstruct_mv.h> @@ -1040,31 +1040,31 @@ void HypreExtMultiABec::boundaryDterm(int level, const Mask &msk = bd[level]->bndryMasks(oitr(), i); if (reg[oitr()] == domain[oitr()]) { - const int *tfp = NULL; + Array4 tf_arr; int bctype = bct; if (bd[level]->mixedBndry(oitr())) { const BaseFab &tf = *(bd[level]->bndryTypes(oitr())[i]); - tfp = tf.dataPtr(); + tf_arr = tf.array(); bctype = -1; } - hdterm3(BL_TO_FORTRAN(Dterm[idim][mfi]), - BL_TO_FORTRAN_N(Soln[mfi], icomp), - ARLIM(reg.loVect()), ARLIM(reg.hiVect()), - cdir, bctype, tfp, bcl, - BL_TO_FORTRAN_N(bcv, bdcomp), - BL_TO_FORTRAN(msk), - BL_TO_FORTRAN((*d2coefs[level])[idim][mfi]), - geom[level].CellSize()); + HABEC::hdterm3(Dterm[idim][mfi].array(), + Soln[mfi].array(icomp), + reg, + cdir, bctype, tf_arr, bcl, + bcv.array(bdcomp), + msk.array(), + (*d2coefs[level])[idim][mfi].array(), + geom[level].CellSize()); } else { - hdterm(BL_TO_FORTRAN(Dterm[idim][mfi]), - BL_TO_FORTRAN_N(Soln[mfi], icomp), - ARLIM(reg.loVect()), ARLIM(reg.hiVect()), - cdir, bct, bcl, - BL_TO_FORTRAN_N(bcv, bdcomp), - BL_TO_FORTRAN(msk), - BL_TO_FORTRAN((*d2coefs[level])[idim][mfi]), - geom[level].CellSize()); + HABEC::hdterm(Dterm[idim][mfi].array(), + Soln[mfi].array(icomp), + reg, + cdir, bct, bcl, + bcv.array(bdcomp), + msk.array(), + (*d2coefs[level])[idim][mfi].array(), + geom[level].CellSize()); } } } diff --git a/Source/radiation/HypreMultiABec.cpp b/Source/radiation/HypreMultiABec.cpp index b548767211..311c78c0a4 100644 --- a/Source/radiation/HypreMultiABec.cpp +++ b/Source/radiation/HypreMultiABec.cpp @@ -4,7 +4,6 @@ #include #include -#include #include #include diff --git a/Source/radiation/Make.package b/Source/radiation/Make.package index 1def860164..eb9eb54a20 100644 --- a/Source/radiation/Make.package +++ b/Source/radiation/Make.package @@ -30,10 +30,8 @@ CEXE_headers += filter.H CEXE_headers += filt_prim.H FEXE_headers += RAD_F.H -FEXE_headers += HABEC_F.H ca_F90EXE_sources += RAD_$(DIM)D.F90 -ca_F90EXE_sources += HABEC_$(DIM)D.F90 CEXE_sources += trace_ppm_rad.cpp diff --git a/Source/radiation/RadSolve.cpp b/Source/radiation/RadSolve.cpp index 12c8fd7b5f..ea8e78d0df 100644 --- a/Source/radiation/RadSolve.cpp +++ b/Source/radiation/RadSolve.cpp @@ -8,7 +8,6 @@ #include #include #include -#include // only for nonsymmetric flux; may be changed? #include From d2ace11ffe41866da4e0d21f3a8975feef5df4d8 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 10 Oct 2023 20:16:55 -0400 Subject: [PATCH 3/4] document the shock burning flag (#2631) --- Docs/source/reactions.rst | 28 ++++++++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/Docs/source/reactions.rst b/Docs/source/reactions.rst index 21d6c64091..cbef46f81d 100644 --- a/Docs/source/reactions.rst +++ b/Docs/source/reactions.rst @@ -57,8 +57,32 @@ in the inputs file. Reactions are enabled by setting:: the efficiency). It is possible to set the maximum and minimum temperature and density for allowing -reactions to occur in a zone using the parameters ``castro.react_T_min``, -``castro.react_T_max``, ``castro.react_rho_min`` and ``castro.react_rho_max``. +reactions to occur in a zone using the parameters: + +* ``castro.react_T_min`` and ``castro.react_T_max`` for temperature + +* ``castro.react_rho_min`` and ``castro.react_rho_max`` for density + +.. index:: castro.disable_shock_burning, USE_SHOCK_VAR + +Burning can also be disabled inside shocks. This requires that the code be +compiled with:: + + USE_SHOCK_VAR = TRUE + +in the ``GNUmakefile``. This will allocate storage for a shock flag in the conserved +state array. This flag is computed via a multidimensional shock detection algorithm +that looks for compression (:math:`\nabla \cdot \ub < 0`) along with a pressure jump +in the direction of compression. The runtime parameter:: + + castro.disable_shock_burning = 1 + +will skip reactions in a zone where we've detected a shock. + +.. note:: + + Both the compilation with ``USE_SHOCK_VAR = TRUE`` and the runtime parameter + ``castro.disable_shock_burning = 1`` are needed to turn off burning in shocks. Reactions Flowchart =================== From f8446cd94c1d54ab003e3021622efea31c9f9b41 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 12 Oct 2023 13:40:00 -0400 Subject: [PATCH 4/4] some changes to the SDC Newton solver to make it more robust (#2586) we now reject a solve if the mass fractions are too far out of [0, 1] also do the total energy update conservatively --- .../react_converge_128_true_sdc.out | 18 ++++---- Source/sdc/Castro_sdc.cpp | 16 ++++++- Source/sdc/Castro_sdc_util.H | 11 +---- Source/sdc/sdc_newton_solve.H | 42 +++++++++++++++---- 4 files changed, 58 insertions(+), 29 deletions(-) diff --git a/Exec/reacting_tests/reacting_convergence/ci-benchmarks/react_converge_128_true_sdc.out b/Exec/reacting_tests/reacting_convergence/ci-benchmarks/react_converge_128_true_sdc.out index 2357cefc1f..35ea0d15ee 100644 --- a/Exec/reacting_tests/reacting_convergence/ci-benchmarks/react_converge_128_true_sdc.out +++ b/Exec/reacting_tests/reacting_convergence/ci-benchmarks/react_converge_128_true_sdc.out @@ -21,16 +21,16 @@ rho_e_sdc_source -6.1984065959e+22 6.3082607681e+21 Temp_sdc_source 0 0 rho_He4_sdc_source -543231.29383 78510.973919 - rho_C12_sdc_source -1694.4890575 10.073635142 + rho_C12_sdc_source -1694.4890575 10.073635144 rho_O16_sdc_source -0.0098093776861 7.1506350196e-06 rho_Fe56_sdc_source -5.4492579275e-05 7.8513247889e-06 pressure 1.4155320805e+22 4.2608130858e+22 - kineng 5.228354476e-13 1.6647276226e+18 + kineng 3.6057617076e-14 1.6647276226e+18 soundspeed 212069503.63 257404342.21 Gamma_1 1.5601126452 1.5885713572 - MachNumber 6.8192135042e-18 0.0086982771596 + MachNumber 1.7908132003e-18 0.0086982771596 entropy 348439780.9 349209883.57 - magvort 6.3108872418e-30 0.00018541051835 + magvort 0 0.00018541051835 divu -0.1163387912 0.55816306007 eint_E 4.5262357143e+16 6.9937847678e+16 eint_e 4.5262357143e+16 6.9937843728e+16 @@ -49,11 +49,11 @@ z_velocity 0 0 t_sound_t_enuc 0.00023710316663 0.0195732693 enuc 2.9131490847e+15 4.5102586513e+17 - magvel 1.446147223e-09 2067446.1363 - radvel -0.00067837194793 2067446.1363 + magvel 3.7977686648e-10 2067446.1363 + radvel -0.00067839201193 2067446.1363 circvel 0 11.820144682 - magmom 0.00072307361144 1.6422547233e+12 - angular_momentum_x -0 -0 + magmom 0.00018988843323 1.6422547233e+12 + angular_momentum_x 0 0 angular_momentum_y 0 0 - angular_momentum_z -1.2410862734e+14 1.2410862747e+14 + angular_momentum_z -1.241086281e+14 1.2410862748e+14 diff --git a/Source/sdc/Castro_sdc.cpp b/Source/sdc/Castro_sdc.cpp index 8e87268882..d82d2764e9 100644 --- a/Source/sdc/Castro_sdc.cpp +++ b/Source/sdc/Castro_sdc.cpp @@ -208,8 +208,6 @@ Castro::do_sdc_update(int m_start, int m_end, Real dt) { normalize_species_sdc(i, j, k, U_center_arr); }); - // ca_normalize_species(AMREX_INT_ANYD(bx1.loVect()), AMREX_INT_ANYD(bx1.hiVect()), - // BL_TO_FORTRAN_ANYD(U_center)); // convert the C source to cell-centers C_center.resize(bx1, NUM_STATE); @@ -234,6 +232,13 @@ Castro::do_sdc_update(int m_start, int m_end, Real dt) sdc_update_centers_o4(i, j, k, U_center_arr, U_new_center_arr, C_center_arr, dt_m, sdc_iteration); }); + // enforce that the species sum to one after the reaction solve + amrex::ParallelFor(bx1, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + normalize_species_sdc(i, j, k, U_new_center_arr); + }); + // compute R_i and in 1 ghost cell and then convert to in // place (only for the interior) R_new.resize(bx1, NUM_STATE); @@ -350,6 +355,13 @@ Castro::construct_old_react_source(MultiFab& U_state, make_cell_center(obx, U_state.array(mfi), U_center_arr, domain_lo, domain_hi); + // sometimes the Laplacian can make the species go negative near discontinuities + amrex::ParallelFor(obx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + normalize_species_sdc(i, j, k, U_center_arr); + }); + // burn, including one ghost cell R_center.resize(obx, NUM_STATE); Elixir elix_r_center = R_center.elixir(); diff --git a/Source/sdc/Castro_sdc_util.H b/Source/sdc/Castro_sdc_util.H index 3a4e71d2d4..c08b289d78 100644 --- a/Source/sdc/Castro_sdc_util.H +++ b/Source/sdc/Castro_sdc_util.H @@ -52,13 +52,6 @@ sdc_solve(const Real dt_m, int ierr; Real err_out; - // for debugging - GpuArray U_orig; - - for (int n = 0; n < NUM_STATE; ++n) { - U_orig[n] = U_old[n]; - } - if (sdc_solver == NEWTON_SOLVE) { // We are going to assume we already have a good guess // for the solve in U_new and just pass the solve onto @@ -66,7 +59,7 @@ sdc_solve(const Real dt_m, sdc_newton_subdivide(dt_m, U_old, U_new, C, sdc_iteration, err_out, ierr); // failing? - if (ierr != NEWTON_SUCCESS) { + if (ierr != newton::NEWTON_SUCCESS) { Abort("Newton subcycling failed in sdc_solve"); } } else if (sdc_solver == VODE_SOLVE) { @@ -85,7 +78,7 @@ sdc_solve(const Real dt_m, sdc_newton_subdivide(dt_m, U_old, U_new, C, sdc_iteration, err_out, ierr); // Failing? - if (ierr != NEWTON_SUCCESS) { + if (ierr != newton::NEWTON_SUCCESS) { Abort("Newton failure in sdc_solve"); } } diff --git a/Source/sdc/sdc_newton_solve.H b/Source/sdc/sdc_newton_solve.H index 313632746c..500e31092e 100644 --- a/Source/sdc/sdc_newton_solve.H +++ b/Source/sdc/sdc_newton_solve.H @@ -5,10 +5,14 @@ #include // error codes -constexpr int NEWTON_SUCCESS = 0; -constexpr int SINGULAR_MATRIX = -1; -constexpr int CONVERGENCE_FAILURE = -2; +namespace newton { + constexpr int NEWTON_SUCCESS = 0; + constexpr int SINGULAR_MATRIX = -1; + constexpr int CONVERGENCE_FAILURE = -2; + constexpr int BAD_MASS_FRACTIONS = -3; + constexpr Real species_failure_tolerance = 1.e-2_rt; +}; #ifdef REACTIONS @@ -107,7 +111,7 @@ sdc_newton_solve(const Real dt_m, const int MAX_ITER = 100; - ierr = NEWTON_SUCCESS; + ierr = newton::NEWTON_SUCCESS; // update the density and momenta for this zone -- they don't react @@ -162,7 +166,7 @@ sdc_newton_solve(const Real dt_m, dgefa(Jac, ipvt, info); #endif if (info != 0) { - ierr = SINGULAR_MATRIX; + ierr = newton::SINGULAR_MATRIX; return; } @@ -206,7 +210,7 @@ sdc_newton_solve(const Real dt_m, err_out = err; if (!converged) { - ierr = CONVERGENCE_FAILURE; + ierr = newton::CONVERGENCE_FAILURE; return; } @@ -221,7 +225,13 @@ sdc_newton_solve(const Real dt_m, } U_new[UEINT] = burn_state.y[SEINT]; - U_new[UEDEN] = U_new[UEINT] + 0.5_rt * v2 / U_new[URHO]; + + // we want to do a conservative update for (rho E), so first figure out the + // energy generation rate + + Real rho_Sdot = (U_new[UEINT] - U_old[UEINT]) / dt_m - C[UEINT]; + + U_new[UEDEN] = U_old[UEDEN] + dt_m * (C[UEDEN] + rho_Sdot); } @@ -251,13 +261,13 @@ sdc_newton_subdivide(const Real dt_m, // use the old time solution. int nsub = 1; - ierr = CONVERGENCE_FAILURE; + ierr = newton::CONVERGENCE_FAILURE; for (int n = 0; n < NUM_STATE; ++n) { U_begin[n] = U_old[n]; } - while (nsub < MAX_NSUB && ierr != NEWTON_SUCCESS) { + while (nsub < MAX_NSUB && ierr != newton::NEWTON_SUCCESS) { if (nsub > 1) { for (int n = 0; n < NUM_STATE; ++n) { U_new[n] = U_old[n]; @@ -278,6 +288,20 @@ sdc_newton_subdivide(const Real dt_m, sdc_newton_solve(dt_sub, U_begin, U_new, C, sdc_iteration, err_out, ierr); + // our solve may have resulted in mass fractions outside + // of [0, 1] -- reject if this is the case + for (int n = 0; n < NumSpec; ++n) { + if (U_new[UFS+n] < -newton::species_failure_tolerance * U_new[URHO] || + U_new[UFS+n] > (1.0_rt + newton::species_failure_tolerance) * U_new[URHO]) { + ierr = newton::BAD_MASS_FRACTIONS; + } + } + + if (ierr != newton::NEWTON_SUCCESS) { + // no point in continuing this subdivision if one stage failed + break; + } + for (int n = 0; n < NUM_STATE; ++n) { U_begin[n] = U_new[n]; }