diff --git a/Source/hydro/Castro_mol.cpp b/Source/hydro/Castro_mol.cpp index eb6314769c..fe2004a60c 100644 --- a/Source/hydro/Castro_mol.cpp +++ b/Source/hydro/Castro_mol.cpp @@ -8,6 +8,10 @@ #include #endif +#ifdef HYBRID_MOMENTUM +#include +#endif + #include #include @@ -321,3 +325,143 @@ Castro::mol_diffusive_flux(const Box& bx, }); } + + +void +Castro::compute_flux_from_q(const Box& bx, + Array4 const& qint, + Array4 const& F, + const int idir) { + + // given a primitive state, compute the flux in direction idir + // + + int iu, iv1, iv2; + int im1, im2, im3; + + auto coord = geom.Coord(); + auto mom_check = mom_flux_has_p(idir, idir, coord); + + if (idir == 0) { + iu = QU; + iv1 = QV; + iv2 = QW; + im1 = UMX; + im2 = UMY; + im3 = UMZ; + + } else if (idir == 1) { + iu = QV; + iv1 = QU; + iv2 = QW; + im1 = UMY; + im2 = UMX; + im3 = UMZ; + + } else { + iu = QW; + iv1 = QU; + iv2 = QV; + im1 = UMZ; + im2 = UMX; + im3 = UMY; + } + +#ifdef HYBRID_MOMENTUM + GeometryData geomdata = geom.data(); +#endif + + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + + Real u_adv = qint(i,j,k,iu); + Real rhoeint = qint(i,j,k,QREINT); + + // Compute fluxes, order as conserved state (not q) + F(i,j,k,URHO) = qint(i,j,k,QRHO)*u_adv; + + F(i,j,k,im1) = F(i,j,k,URHO)*qint(i,j,k,iu); + if (mom_check) { + F(i,j,k,im1) += qint(i,j,k,QPRES); + } + F(i,j,k,im2) = F(i,j,k,URHO)*qint(i,j,k,iv1); + F(i,j,k,im3) = F(i,j,k,URHO)*qint(i,j,k,iv2); + + Real rhoetot = rhoeint + 0.5_rt * qint(i,j,k,QRHO)* + (qint(i,j,k,iu)*qint(i,j,k,iu) + + qint(i,j,k,iv1)*qint(i,j,k,iv1) + + qint(i,j,k,iv2)*qint(i,j,k,iv2)); + + F(i,j,k,UEDEN) = u_adv*(rhoetot + qint(i,j,k,QPRES)); + F(i,j,k,UEINT) = u_adv*rhoeint; + + F(i,j,k,UTEMP) = 0.0; +#ifdef SHOCK_VAR + F(i,j,k,USHK) = 0.0; +#endif + +#ifdef NSE_NET + F(i,j,k,UMUP) = 0.0; + F(i,j,k,UMUN) = 0.0; +#endif + // passively advected quantities + for (int ipassive = 0; ipassive < npassive; ipassive++) { + int n = upassmap(ipassive); + int nqp = qpassmap(ipassive); + + F(i,j,k,n) = F(i,j,k,URHO)*qint(i,j,k,nqp); + } + +#ifdef HYBRID_MOMENTUM + // the hybrid routine uses the Godunov indices, not the full NQ state + GpuArray qgdnv_zone; + qgdnv_zone[GDRHO] = qint(i,j,k,QRHO); + qgdnv_zone[GDU] = qint(i,j,k,QU); + qgdnv_zone[GDV] = qint(i,j,k,QV); + qgdnv_zone[GDW] = qint(i,j,k,QW); + qgdnv_zone[GDPRES] = qint(i,j,k,QPRES); + GpuArray F_zone; + for (int n = 0; n < NUM_STATE; n++) { + F_zone[n] = F(i,j,k,n); + } + compute_hybrid_flux(qgdnv_zone, geomdata, idir, i, j, k, F_zone); + for (int n = 0; n < NUM_STATE; n++) { + F(i,j,k,n) = F_zone[n]; + } +#endif + }); +} + +void +Castro::store_godunov_state(const Box& bx, + Array4 const& qint, +#ifdef RADIATION + Array4 const& lambda, +#endif + Array4 const& qgdnv) { + + // this copies the full interface state (NQ -- one for each primitive + // variable) over to a smaller subset of size NGDNV for use later in the + // hydro advancement. + + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + + +#ifdef HYBRID_MOMENTUM + qgdnv(i,j,k,GDRHO) = qint(i,j,k,QRHO); +#endif + qgdnv(i,j,k,GDU) = qint(i,j,k,QU); + qgdnv(i,j,k,GDV) = qint(i,j,k,QV); + qgdnv(i,j,k,GDW) = qint(i,j,k,QW); + qgdnv(i,j,k,GDPRES) = qint(i,j,k,QPRES); +#ifdef RADIATION + for (int g = 0; g < NGROUPS; g++) { + qgdnv(i,j,k,GDLAMS+g) = lambda(i,j,k,g); + qgdnv(i,j,k,GDERADS+g) = qint(i,j,k,QRAD+g); + } +#endif + }); +} diff --git a/Source/hydro/Make.package b/Source/hydro/Make.package index b4d67b48fa..6b6ba6c7e8 100644 --- a/Source/hydro/Make.package +++ b/Source/hydro/Make.package @@ -30,7 +30,6 @@ CEXE_sources += riemann.cpp CEXE_headers += HLL_solvers.H CEXE_headers += riemann_solvers.H CEXE_headers += riemann_2shock_solvers.H -CEXE_sources += riemann_util.cpp CEXE_headers += riemann_type.H CEXE_headers += slope.H CEXE_headers += reconstruction.H diff --git a/Source/hydro/riemann_util.cpp b/Source/hydro/riemann_util.cpp deleted file mode 100644 index 2801dad8c2..0000000000 --- a/Source/hydro/riemann_util.cpp +++ /dev/null @@ -1,154 +0,0 @@ -#include -#include - -#ifdef RADIATION -#include -#include -#endif - -#ifdef HYBRID_MOMENTUM -#include -#endif - -#include - -using namespace amrex; - -void -Castro::compute_flux_from_q(const Box& bx, - Array4 const& qint, - Array4 const& F, - const int idir) { - - // given a primitive state, compute the flux in direction idir - // - - int iu, iv1, iv2; - int im1, im2, im3; - - auto coord = geom.Coord(); - auto mom_check = mom_flux_has_p(idir, idir, coord); - - if (idir == 0) { - iu = QU; - iv1 = QV; - iv2 = QW; - im1 = UMX; - im2 = UMY; - im3 = UMZ; - - } else if (idir == 1) { - iu = QV; - iv1 = QU; - iv2 = QW; - im1 = UMY; - im2 = UMX; - im3 = UMZ; - - } else { - iu = QW; - iv1 = QU; - iv2 = QV; - im1 = UMZ; - im2 = UMX; - im3 = UMY; - } - -#ifdef HYBRID_MOMENTUM - GeometryData geomdata = geom.data(); -#endif - - amrex::ParallelFor(bx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - - Real u_adv = qint(i,j,k,iu); - Real rhoeint = qint(i,j,k,QREINT); - - // Compute fluxes, order as conserved state (not q) - F(i,j,k,URHO) = qint(i,j,k,QRHO)*u_adv; - - F(i,j,k,im1) = F(i,j,k,URHO)*qint(i,j,k,iu); - if (mom_check) { - F(i,j,k,im1) += qint(i,j,k,QPRES); - } - F(i,j,k,im2) = F(i,j,k,URHO)*qint(i,j,k,iv1); - F(i,j,k,im3) = F(i,j,k,URHO)*qint(i,j,k,iv2); - - Real rhoetot = rhoeint + 0.5_rt * qint(i,j,k,QRHO)* - (qint(i,j,k,iu)*qint(i,j,k,iu) + - qint(i,j,k,iv1)*qint(i,j,k,iv1) + - qint(i,j,k,iv2)*qint(i,j,k,iv2)); - - F(i,j,k,UEDEN) = u_adv*(rhoetot + qint(i,j,k,QPRES)); - F(i,j,k,UEINT) = u_adv*rhoeint; - - F(i,j,k,UTEMP) = 0.0; -#ifdef SHOCK_VAR - F(i,j,k,USHK) = 0.0; -#endif - -#ifdef NSE_NET - F(i,j,k,UMUP) = 0.0; - F(i,j,k,UMUN) = 0.0; -#endif - // passively advected quantities - for (int ipassive = 0; ipassive < npassive; ipassive++) { - int n = upassmap(ipassive); - int nqp = qpassmap(ipassive); - - F(i,j,k,n) = F(i,j,k,URHO)*qint(i,j,k,nqp); - } - -#ifdef HYBRID_MOMENTUM - // the hybrid routine uses the Godunov indices, not the full NQ state - GpuArray qgdnv_zone; - qgdnv_zone[GDRHO] = qint(i,j,k,QRHO); - qgdnv_zone[GDU] = qint(i,j,k,QU); - qgdnv_zone[GDV] = qint(i,j,k,QV); - qgdnv_zone[GDW] = qint(i,j,k,QW); - qgdnv_zone[GDPRES] = qint(i,j,k,QPRES); - GpuArray F_zone; - for (int n = 0; n < NUM_STATE; n++) { - F_zone[n] = F(i,j,k,n); - } - compute_hybrid_flux(qgdnv_zone, geomdata, idir, i, j, k, F_zone); - for (int n = 0; n < NUM_STATE; n++) { - F(i,j,k,n) = F_zone[n]; - } -#endif - }); -} - -void -Castro::store_godunov_state(const Box& bx, - Array4 const& qint, -#ifdef RADIATION - Array4 const& lambda, -#endif - Array4 const& qgdnv) { - - // this copies the full interface state (NQ -- one for each primitive - // variable) over to a smaller subset of size NGDNV for use later in the - // hydro advancement. - - amrex::ParallelFor(bx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - - -#ifdef HYBRID_MOMENTUM - qgdnv(i,j,k,GDRHO) = qint(i,j,k,QRHO); -#endif - qgdnv(i,j,k,GDU) = qint(i,j,k,QU); - qgdnv(i,j,k,GDV) = qint(i,j,k,QV); - qgdnv(i,j,k,GDW) = qint(i,j,k,QW); - qgdnv(i,j,k,GDPRES) = qint(i,j,k,QPRES); -#ifdef RADIATION - for (int g = 0; g < NGROUPS; g++) { - qgdnv(i,j,k,GDLAMS+g) = lambda(i,j,k,g); - qgdnv(i,j,k,GDERADS+g) = qint(i,j,k,QRAD+g); - } -#endif - }); -}