From 29372b3b6343d3c3d9abfae537941564c8114895 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Fri, 16 Aug 2024 17:47:32 -0700 Subject: [PATCH 01/10] add options for upwind and minmod, used by Amr-Wind --- Godunov/hydro_godunov_edge_state_2D.cpp | 30 +- Godunov/hydro_godunov_edge_state_3D.cpp | 38 ++- Godunov/hydro_godunov_ppm.H | 409 ++++++++++++++++++------ 3 files changed, 368 insertions(+), 109 deletions(-) diff --git a/Godunov/hydro_godunov_edge_state_2D.cpp b/Godunov/hydro_godunov_edge_state_2D.cpp index 3c9ba88fd..389148a4a 100644 --- a/Godunov/hydro_godunov_edge_state_2D.cpp +++ b/Godunov/hydro_godunov_edge_state_2D.cpp @@ -81,25 +81,45 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, PPM::PredictStateOnFaces(bxg1,AMREX_D_DECL(Imx,Imy,Imz), AMREX_D_DECL(Ipx,Ipy,Ipz), AMREX_D_DECL(umac,vmac,wmac), - q,geom,l_dt,pbc,ncomp,limiter); + q,geom,l_dt,pbc,ncomp,limiter, + limiter_type); } else if ( limiter_type == PPM::WENOZ) { auto limiter = PPM::wenoz(); PPM::PredictStateOnFaces(bxg1,AMREX_D_DECL(Imx,Imy,Imz), AMREX_D_DECL(Ipx,Ipy,Ipz), AMREX_D_DECL(umac,vmac,wmac), - q,geom,l_dt,pbc,ncomp,limiter); + q,geom,l_dt,pbc,ncomp,limiter, + limiter_type); } else if ( limiter_type == PPM::WENO_JS) { auto limiter = PPM::weno_js(); PPM::PredictStateOnFaces(bxg1,AMREX_D_DECL(Imx,Imy,Imz), AMREX_D_DECL(Ipx,Ipy,Ipz), AMREX_D_DECL(umac,vmac,wmac), - q,geom,l_dt,pbc,ncomp,limiter); - } else { + q,geom,l_dt,pbc,ncomp,limiter, + limiter_type); + } else if ( limiter_type == PPM::NoLimiter) { auto limiter = PPM::nolimiter(); PPM::PredictStateOnFaces(bxg1,AMREX_D_DECL(Imx,Imy,Imz), AMREX_D_DECL(Ipx,Ipy,Ipz), AMREX_D_DECL(umac,vmac,wmac), - q,geom,l_dt,pbc,ncomp,limiter); + q,geom,l_dt,pbc,ncomp,limiter, + limiter_type); + } else if ( limiter_type == PPM::UPWIND) { + auto limiter = PPM::upwind(); + PPM::PredictStateOnFaces(bxg1,AMREX_D_DECL(Imx,Imy,Imz), + AMREX_D_DECL(Ipx,Ipy,Ipz), + AMREX_D_DECL(umac,vmac,wmac), + q,geom,l_dt,pbc,ncomp,limiter, + limiter_type); + } else if ( limiter_type == PPM::MINMOD) { + auto limiter = PPM::minmod(); + PPM::PredictStateOnFaces(bxg1,AMREX_D_DECL(Imx,Imy,Imz), + AMREX_D_DECL(Ipx,Ipy,Ipz), + AMREX_D_DECL(umac,vmac,wmac), + q,geom,l_dt,pbc,ncomp,limiter, + limiter_type); + } else { + amrex::Abort("Unknown limiter_type in hydro_godunov_edge_state_2D.cpp"); } // Use PLM to generate Im and Ip */ } diff --git a/Godunov/hydro_godunov_edge_state_3D.cpp b/Godunov/hydro_godunov_edge_state_3D.cpp index faa74d3dc..2247c16e3 100644 --- a/Godunov/hydro_godunov_edge_state_3D.cpp +++ b/Godunov/hydro_godunov_edge_state_3D.cpp @@ -96,31 +96,51 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, PPM::PredictStateOnFaces(bxg1,AMREX_D_DECL(Imx,Imy,Imz), AMREX_D_DECL(Ipx,Ipy,Ipz), AMREX_D_DECL(umac,vmac,wmac), - q,geom,l_dt,pbc,ncomp,limiter); + q,geom,l_dt,pbc,ncomp,limiter, + limiter_type); } else if ( limiter_type == PPM::WENOZ) { auto limiter = PPM::wenoz(); PPM::PredictStateOnFaces(bxg1,AMREX_D_DECL(Imx,Imy,Imz), AMREX_D_DECL(Ipx,Ipy,Ipz), AMREX_D_DECL(umac,vmac,wmac), - q,geom,l_dt,pbc,ncomp,limiter); + q,geom,l_dt,pbc,ncomp,limiter, + limiter_type); } else if ( limiter_type == PPM::WENO_JS) { auto limiter = PPM::weno_js(); PPM::PredictStateOnFaces(bxg1,AMREX_D_DECL(Imx,Imy,Imz), AMREX_D_DECL(Ipx,Ipy,Ipz), AMREX_D_DECL(umac,vmac,wmac), - q,geom,l_dt,pbc,ncomp,limiter); - } else { + q,geom,l_dt,pbc,ncomp,limiter, + limiter_type); + } else if ( limiter_type == PPM::NoLimiter) { auto limiter = PPM::nolimiter(); PPM::PredictStateOnFaces(bxg1,AMREX_D_DECL(Imx,Imy,Imz), AMREX_D_DECL(Ipx,Ipy,Ipz), AMREX_D_DECL(umac,vmac,wmac), - q,geom,l_dt,pbc,ncomp,limiter); + q,geom,l_dt,pbc,ncomp,limiter, + limiter_type); + } else if ( limiter_type == PPM::UPWIND) { + auto limiter = PPM::upwind(); + PPM::PredictStateOnFaces(bxg1,AMREX_D_DECL(Imx,Imy,Imz), + AMREX_D_DECL(Ipx,Ipy,Ipz), + AMREX_D_DECL(umac,vmac,wmac), + q,geom,l_dt,pbc,ncomp,limiter, + limiter_type); + } else if ( limiter_type == PPM::MINMOD) { + auto limiter = PPM::minmod(); + PPM::PredictStateOnFaces(bxg1,AMREX_D_DECL(Imx,Imy,Imz), + AMREX_D_DECL(Ipx,Ipy,Ipz), + AMREX_D_DECL(umac,vmac,wmac), + q,geom,l_dt,pbc,ncomp,limiter, + limiter_type); + } else { + amrex::Abort("Unknown limiter_type in hydro_godunov_edge_state_3D.cpp"); } - // Use PLM to generate Im and Ip */ - } - else - { + } else { + // + // Use PLM to generate Im and Ip */ + // amrex::ParallelFor(xebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { diff --git a/Godunov/hydro_godunov_ppm.H b/Godunov/hydro_godunov_ppm.H index 739270189..88ae51b94 100644 --- a/Godunov/hydro_godunov_ppm.H +++ b/Godunov/hydro_godunov_ppm.H @@ -5,17 +5,73 @@ * @{ */ -#ifndef HYDRO_PPM_GODUNOV_H -#define HYDRO_PPM_GODUNOV_H +#ifndef HYDRO_GODUNOV_PPM_H +#define HYDRO_GODUNOV_PPM_H #include #include #include #include +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void +minmod_fn (const amrex::Real sm1, + const amrex::Real s0, + const amrex::Real sp1, + amrex::Real& dsm, + amrex::Real& dsp) +{ + // Calculate gradients on both sides + dsp = sp1 - s0; + dsm = s0 - sm1; + + if (!(dsp * dsm < 0.0)) { + // Select the smaller slope if same sign + if (std::abs(dsp) < std::abs(dsm)) { + dsm = dsp; + } else { + dsp = dsm; + } + } else { + // Set to zero if opposite sign + dsp = 0.0; + dsm = 0.0; + } +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void +Godunov_minmod_bc (const int n, + const amrex::Real sm1, + const amrex::Real s0, + const amrex::Real sp1, + amrex::Real& dsm, + amrex::Real& dsp, + const int bclo, + const int bchi, + const int domlo, + const int domhi) +{ + using namespace amrex; + + if (bclo == BCType::ext_dir || bclo == BCType::hoextrap) { + if (n == domlo) { + // Ensure that left-side slope is used unchanged + dsm = s0 - sm1; + } + } + + if (bchi == BCType::ext_dir || bchi == BCType::hoextrap) { + if (n == domhi) { + // Ensure that the right-side slope is used unchanged + dsp = sp1 - s0; + } + } +} + namespace PPM { -enum limiters {VanLeer, WENOZ, WENO_JS, NoLimiter}; +enum limiters {VanLeer, WENOZ, WENO_JS, NoLimiter, UPWIND, MINMOD}; static constexpr int default_limiter = VanLeer; @@ -37,7 +93,6 @@ struct nolimiter { amrex::Real sedge = amrex::Real(0.5)*(s0 + sm1) - amrex::Real(1./6.)*(d1 - d2); return sedge; - } [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -63,7 +118,7 @@ struct nolimiter { { return amrex::makeTuple(sedge1, sedge2); } -}; +}; // nolimiter struct vanleer { @@ -152,7 +207,7 @@ struct vanleer { // return amrex::makeTuple(sm_,sp_); } -}; +}; // vanleer struct wenoz { @@ -211,7 +266,7 @@ struct wenoz { { return amrex::makeTuple(sedge1, sedge2); } -}; +}; // wenoa struct weno_js { @@ -274,7 +329,79 @@ struct weno_js { { return amrex::makeTuple(sedge1, sedge2); } -}; +}; // weno_js + +struct upwind { + + static constexpr bool do_limiting = false; + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real + sedge1(const amrex::Real sm2, + const amrex::Real sm1, + const amrex::Real s0, + const amrex::Real sp1, + const amrex::Real /*sp2*/) + { + return 0; + } + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real + sedge2(const amrex::Real sm2, + const amrex::Real sm1, + const amrex::Real s0, + const amrex::Real sp1, + const amrex::Real /*sp2*/) + { + return 0; + } + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::GpuTuple + sm_sp(const amrex::Real /*s0*/, + const amrex::Real sedge1, + const amrex::Real sedge2) + { + return amrex::makeTuple(sedge1, sedge2); + } +}; // upwind + +struct minmod { + + static constexpr bool do_limiting = false; + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real + sedge1(const amrex::Real sm2, + const amrex::Real sm1, + const amrex::Real s0, + const amrex::Real sp1, + const amrex::Real /*sp2*/) + { + return 0; + } + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real + sedge2(const amrex::Real sm2, + const amrex::Real sm1, + const amrex::Real s0, + const amrex::Real sp1, + const amrex::Real /*sp2*/) + { + return 0; + } + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::GpuTuple + sm_sp(const amrex::Real /*s0*/, + const amrex::Real sedge1, + const amrex::Real sedge2) + { + return amrex::makeTuple(sedge1, sedge2); + } +}; // minmod template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -578,6 +705,7 @@ void PredictVelOnXFace ( const int i, const int j, const int k, const int n, const Limiter& /*limiter*/) { constexpr auto half{amrex::Real(0.5)}; + constexpr auto one{amrex::Real(1.0)}; constexpr auto two3rds{amrex::Real(2.0/3.0)}; amrex::Real sm2 = S(i-2,j,k,n); @@ -601,13 +729,13 @@ void PredictVelOnXFace ( const int i, const int j, const int k, const int n, if (v_ad > small_vel) { - Ip(i,j,k,n) = sp - (half*sigma)*((sp-sm) - (amrex::Real(1.0) - two3rds*sigma)*s6); + Ip(i,j,k,n) = sp - (half*sigma)*((sp-sm) - (one - two3rds*sigma)*s6); Im(i,j,k,n) = S(i,j,k,n); } else if (v_ad < -small_vel) { Ip(i,j,k,n) = S(i,j,k,n); - Im(i,j,k,n) = sm + (half*sigma)*((sp-sm) + (amrex::Real(1.0) - two3rds*sigma)*s6); + Im(i,j,k,n) = sm + (half*sigma)*((sp-sm) + (one - two3rds*sigma)*s6); } else { Ip(i,j,k,n) = S(i,j,k,n); @@ -626,6 +754,7 @@ void PredictVelOnYFace ( const int i, const int j, const int k, const int n, const Limiter& /*limiter*/) { constexpr auto half{amrex::Real(0.5)}; + constexpr auto one{amrex::Real(1.0)}; constexpr auto two3rds{amrex::Real(2.0/3.0)}; amrex::Real sm2 = S(i,j-2,k,n); @@ -649,13 +778,13 @@ void PredictVelOnYFace ( const int i, const int j, const int k, const int n, if (v_ad > small_vel) { - Ip(i,j,k,n) = sp - (half*sigma)*((sp-sm) - (amrex::Real(1.0) - two3rds*sigma)*s6); + Ip(i,j,k,n) = sp - (half*sigma)*((sp-sm) - (one - two3rds*sigma)*s6); Im(i,j,k,n) = S(i,j,k,n); } else if (v_ad < -small_vel) { Ip(i,j,k,n) = S(i,j,k,n); - Im(i,j,k,n) = sm + (half*sigma)*((sp-sm) + (amrex::Real(1.0) - two3rds*sigma)*s6); + Im(i,j,k,n) = sm + (half*sigma)*((sp-sm) + (one - two3rds*sigma)*s6); } else { Ip(i,j,k,n) = S(i,j,k,n); @@ -675,6 +804,7 @@ void PredictVelOnZFace ( const int i, const int j, const int k, const int n, const Limiter& /*limiter*/) { constexpr auto half{amrex::Real(0.5)}; + constexpr auto one{amrex::Real(1.0)}; constexpr auto two3rds{amrex::Real(2.0/3.0)}; amrex::Real sm2 = S(i,j,k-2,n); @@ -698,13 +828,13 @@ void PredictVelOnZFace ( const int i, const int j, const int k, const int n, if (v_ad > small_vel) { - Ip(i,j,k,n) = sp - (half*sigma)*((sp-sm) - (amrex::Real(1.0) - two3rds*sigma)*s6); + Ip(i,j,k,n) = sp - (half*sigma)*((sp-sm) - (one - two3rds*sigma)*s6); Im(i,j,k,n) = S(i,j,k,n); } else if (v_ad < -small_vel) { Ip(i,j,k,n) = S(i,j,k,n); - Im(i,j,k,n) = sm + (half*sigma)*((sp-sm) + (amrex::Real(1.0) - two3rds*sigma)*s6); + Im(i,j,k,n) = sm + (half*sigma)*((sp-sm) + (one - two3rds*sigma)*s6); } else { Ip(i,j,k,n) = S(i,j,k,n); @@ -726,41 +856,68 @@ void PredictStateOnXFace ( const int i, const int j, const int k, const int n, const amrex::Array4 &vel_edge, const amrex::BCRec bc, const int domlo, const int domhi, - const Limiter& /*limiter*/) + const Limiter& limiter, + int limiter_type) { - constexpr auto half{amrex::Real(0.5)}; - constexpr auto two3rds{amrex::Real(2.0/3.0)}; - amrex::Real sm2 = S(i-2,j,k,n); amrex::Real sm1 = S(i-1,j,k,n); amrex::Real s0 = S(i ,j,k,n); amrex::Real sp1 = S(i+1,j,k,n); amrex::Real sp2 = S(i+2,j,k,n); - amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2); - amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2); - - amrex::Real sm, sp; - amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2); - - SetXBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k), - bc.lo(0), bc.hi(0), domlo, domhi); - - amrex::Real s6 = 6.0*s0 - 3.0*(sm + sp); + constexpr auto half{amrex::Real(0.5)}; + constexpr auto one{amrex::Real(1.0)}; amrex::Real sigmap = amrex::Math::abs(vel_edge(i+1,j,k))*dt/dx; amrex::Real sigmam = amrex::Math::abs(vel_edge(i ,j,k))*dt/dx; - if (vel_edge(i+1,j,k) > small_vel) { - Ip = sp - (half*sigmap)*((sp - sm) - (amrex::Real(1.0) -two3rds*sigmap)*s6); - } else { + if (limiter_type == PPM::UPWIND) { + Im = S(i,j,k,n); Ip = S(i,j,k,n); - } - if(vel_edge(i,j,k) < -small_vel) { - Im = sm + (half*sigmam)*((sp-sm) + (amrex::Real(1.0) - two3rds*sigmam)*s6); + } else if (limiter_type == PPM::MINMOD) { + amrex::Real dsp = 0.0; + amrex::Real dsm = 0.0; + minmod_fn(sm1, s0, sp1, dsm, dsp); + Godunov_minmod_bc(i, sm1, s0, sp1, dsm, dsp, bc.lo(0), bc.hi(0), domlo, domhi); + + if (vel_edge(i + 1, j, k) > small_vel) { + Ip = s0 + half * (one - sigmap) * dsp; + } else { + Ip = S(i,j,k,n); + } + + if (vel_edge(i, j, k) < -small_vel) { + Im = s0 - half * (one - sigmam) * dsm; + } else { + Im = S(i,j,k,n); + } + } else { - Im = S(i,j,k,n); + constexpr auto two3rds{amrex::Real(2.0/3.0)}; + + amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2); + amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2); + + amrex::Real sm, sp; + amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2); + + SetXBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k), + bc.lo(0), bc.hi(0), domlo, domhi); + + amrex::Real s6 = 6.0*s0 - 3.0*(sm + sp); + + if (vel_edge(i+1,j,k) > small_vel) { + Ip = sp - (half*sigmap)*((sp - sm) - (one -two3rds*sigmap)*s6); + } else { + Ip = S(i,j,k,n); + } + + if(vel_edge(i,j,k) < -small_vel) { + Im = sm + (half*sigmam)*((sp-sm) + (one - two3rds*sigmam)*s6); + } else { + Im = S(i,j,k,n); + } } } @@ -773,42 +930,73 @@ void PredictStateOnYFace ( const int i, const int j, const int k, const int n, const amrex::Array4 &vel_edge, const amrex::BCRec bc, const int domlo, const int domhi, - const Limiter& /*limiter*/) + const Limiter& limiter, + int limiter_type) { - constexpr auto half{amrex::Real(0.5)}; - constexpr auto two3rds{amrex::Real(2.0/3.0)}; - amrex::Real sm2 = S(i,j-2,k,n); amrex::Real sm1 = S(i,j-1,k,n); amrex::Real s0 = S(i,j ,k,n); amrex::Real sp1 = S(i,j+1,k,n); amrex::Real sp2 = S(i,j+2,k,n); - amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2); - amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2); - - amrex::Real sm, sp; - amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2); - - SetYBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k), - bc.lo(1), bc.hi(1), domlo, domhi); + constexpr auto half{amrex::Real(0.5)}; + constexpr auto one{amrex::Real(1.0)}; - amrex::Real s6 = 6.0*s0- 3.0*(sm + sp); + amrex::Real s6; amrex::Real sigmap = amrex::Math::abs(vel_edge(i,j+1,k))*dt/dx; amrex::Real sigmam = amrex::Math::abs(vel_edge(i,j ,k))*dt/dx; - if (vel_edge(i,j+1,k) > small_vel) { - Ip = sp - (half*sigmap)*((sp - sm) - (amrex::Real(1.0) -two3rds*sigmap)*s6); - } else { + if (limiter_type == PPM::UPWIND) { + + Im = S(i,j,k,n); Ip = S(i,j,k,n); - } - if (vel_edge(i,j,k) < -small_vel) { - Im = sm + (half*sigmam)*((sp-sm) + (amrex::Real(1.0) - two3rds*sigmam)*s6); + } else if (limiter_type == PPM::MINMOD) { + + amrex::Real dsp = 0.0; + amrex::Real dsm = 0.0; + minmod_fn(sm1, s0, sp1, dsm, dsp); + Godunov_minmod_bc(i, sm1, s0, sp1, dsm, dsp, bc.lo(0), bc.hi(0), domlo, domhi); + + if (vel_edge(i,j+1,k) > small_vel) { + Ip = s0 + half * (one - sigmap) * dsp; + } else { + Ip = S(i,j,k,n); + } + if (vel_edge(i,j,k) < -small_vel) { + Im = s0 - half * (one - sigmam) * dsm; + } else { + Im = S(i,j,k,n); + } + } else { - Im = S(i,j,k,n); - } + + constexpr auto two3rds{amrex::Real(2.0/3.0)}; + + amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2); + amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2); + + amrex::Real sm, sp; + amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2); + + SetYBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k), + bc.lo(1), bc.hi(1), domlo, domhi); + + s6 = 6.0*s0- 3.0*(sm + sp); + + if (vel_edge(i,j+1,k) > small_vel) { + Ip = sp - (half*sigmap)*((sp - sm) - (one -two3rds*sigmap)*s6); + } else { + Ip = S(i,j,k,n); + } + + if (vel_edge(i,j,k) < -small_vel) { + Im = sm + (half*sigmam)*((sp-sm) + (one - two3rds*sigmam)*s6); + } else { + Im = S(i,j,k,n); + } + } // not upwind or minmod } @@ -822,41 +1010,71 @@ void PredictStateOnZFace ( const int i, const int j, const int k, const int n, const amrex::Array4 &vel_edge, const amrex::BCRec bc, const int domlo, const int domhi, - const Limiter& /*limiter*/) + const Limiter& limiter, + int limiter_type) { - - constexpr auto half{amrex::Real(0.5)}; - constexpr auto two3rds{amrex::Real(2.0/3.0)}; - amrex::Real sm2 = S(i,j,k-2,n); amrex::Real sm1 = S(i,j,k-1,n); amrex::Real s0 = S(i,j,k ,n); amrex::Real sp1 = S(i,j,k+1,n); amrex::Real sp2 = S(i,j,k+2,n); - amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2); - amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2); - - amrex::Real sm, sp; - amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2); - - SetZBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k), - bc.lo(2), bc.hi(2), domlo, domhi); + constexpr auto half{amrex::Real(0.5)}; + constexpr auto one{amrex::Real(1.0)}; - amrex::Real s6 = 6.0*s0- 3.0*(sm + sp); amrex::Real sigmap = amrex::Math::abs(vel_edge(i,j,k+1))*dt/dx; amrex::Real sigmam = amrex::Math::abs(vel_edge(i,j,k ))*dt/dx; - if(vel_edge(i,j,k+1) > small_vel) { - Ip = sp - (half*sigmap)*((sp-sm) - (amrex::Real(1.0) -two3rds*sigmap)*s6); - } else { + if (limiter_type == PPM::UPWIND) { + + Im = S(i,j,k,n); Ip = S(i,j,k,n); - } - if(vel_edge(i,j,k) < -small_vel) { - Im = sm + (half*sigmam)*((sp-sm) + (amrex::Real(1.0) - two3rds*sigmam)*s6); + } else if (limiter_type == PPM::MINMOD) { + + amrex::Real dsp = 0.0; + amrex::Real dsm = 0.0; + minmod_fn(sm1, s0, sp1, dsm, dsp); + Godunov_minmod_bc(k, sm1, s0, sp1, dsm, dsp, bc.lo(2), bc.hi(2), domlo, domhi); + + if (vel_edge(i, j, k + 1) > small_vel) { + Ip = s0 + half * (one - sigmap) * dsp; + } else { + Ip = S(i,j,k,n); + } + + if (vel_edge(i, j, k) < -small_vel) { + Im = s0 - half * (one - sigmam) * dsm; + } else { + Im = S(i,j,k,n); + } + } else { - Im = S(i,j,k,n); + constexpr auto half{amrex::Real(0.5)}; + constexpr auto two3rds{amrex::Real(2.0/3.0)}; + + amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2); + amrex::Real sedge2 = Limiter::sedge2(sm2,sm1,s0,sp1,sp2); + + amrex::Real sm, sp; + amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2); + + SetZBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k), + bc.lo(2), bc.hi(2), domlo, domhi); + + amrex::Real s6 = 6.0*s0- 3.0*(sm + sp); + + if(vel_edge(i,j,k+1) > small_vel) { + Ip = sp - (half*sigmap)*((sp-sm) - (amrex::Real(1.0) -two3rds*sigmap)*s6); + } else { + Ip = S(i,j,k,n); + } + + if(vel_edge(i,j,k) < -small_vel) { + Im = sm + (half*sigmam)*((sp-sm) + (amrex::Real(1.0) - two3rds*sigmam)*s6); + } else { + Im = S(i,j,k,n); + } } } #endif @@ -898,21 +1116,22 @@ void PredictVelOnFaces (amrex::Box const& bx, template void PredictStateOnFaces (amrex::Box const& bx, - AMREX_D_DECL(amrex::Array4 const& Imx, - amrex::Array4 const& Imy, - amrex::Array4 const& Imz), - AMREX_D_DECL(amrex::Array4 const& Ipx, - amrex::Array4 const& Ipy, - amrex::Array4 const& Ipz), - AMREX_D_DECL(amrex::Array4 const& umac, - amrex::Array4 const& vmac, - amrex::Array4 const& wmac), - amrex::Array4 const& q, - amrex::Geometry geom, - amrex::Real l_dt, - amrex::BCRec const* pbc, - const int ncomp, - const Limiter& limiter) + AMREX_D_DECL(amrex::Array4 const& Imx, + amrex::Array4 const& Imy, + amrex::Array4 const& Imz), + AMREX_D_DECL(amrex::Array4 const& Ipx, + amrex::Array4 const& Ipy, + amrex::Array4 const& Ipz), + AMREX_D_DECL(amrex::Array4 const& umac, + amrex::Array4 const& vmac, + amrex::Array4 const& wmac), + amrex::Array4 const& q, + amrex::Geometry geom, + amrex::Real l_dt, + amrex::BCRec const* pbc, + const int ncomp, + const Limiter& limiter, + int limiter_type) { const amrex::Box& domain = geom.Domain(); const amrex::Dim3 dlo = amrex::lbound(domain); @@ -926,12 +1145,12 @@ void PredictStateOnFaces (amrex::Box const& bx, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { PPM::PredictStateOnXFace(i, j, k, n, l_dt, dx, Imx(i,j,k,n), Ipx(i,j,k,n), - q, umac, pbc[n], dlo.x, dhi.x,limiter); + q, umac, pbc[n], dlo.x, dhi.x, limiter, limiter_type); PPM::PredictStateOnYFace(i, j, k, n, l_dt, dy, Imy(i,j,k,n), Ipy(i,j,k,n), - q, vmac, pbc[n], dlo.y, dhi.y,limiter); + q, vmac, pbc[n], dlo.y, dhi.y, limiter, limiter_type); #if (AMREX_SPACEDIM==3) PPM::PredictStateOnZFace(i, j, k, n, l_dt, dz, Imz(i,j,k,n), Ipz(i,j,k,n), - q, wmac, pbc[n], dlo.z, dhi.z,limiter); + q, wmac, pbc[n], dlo.z, dhi.z, limiter, limiter_type); #endif }); } From 8f577a4fcf4865f5ab78b69bc587735f94c0fa1c Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Fri, 16 Aug 2024 18:01:12 -0700 Subject: [PATCH 02/10] remove unused --- Godunov/hydro_godunov_ppm.H | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/Godunov/hydro_godunov_ppm.H b/Godunov/hydro_godunov_ppm.H index 88ae51b94..f0f6d0a89 100644 --- a/Godunov/hydro_godunov_ppm.H +++ b/Godunov/hydro_godunov_ppm.H @@ -337,10 +337,10 @@ struct upwind { [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real - sedge1(const amrex::Real sm2, - const amrex::Real sm1, - const amrex::Real s0, - const amrex::Real sp1, + sedge1(const amrex::Real /*sm2*/, + const amrex::Real /*sm1*/, + const amrex::Real /*s0*/, + const amrex::Real /*sp1*/, const amrex::Real /*sp2*/) { return 0; @@ -348,10 +348,10 @@ struct upwind { [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real - sedge2(const amrex::Real sm2, - const amrex::Real sm1, - const amrex::Real s0, - const amrex::Real sp1, + sedge2(const amrex::Real /*sm2*/, + const amrex::Real /*sm1*/, + const amrex::Real /*s0*/, + const amrex::Real /*sp1*/, const amrex::Real /*sp2*/) { return 0; @@ -373,10 +373,10 @@ struct minmod { [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real - sedge1(const amrex::Real sm2, - const amrex::Real sm1, - const amrex::Real s0, - const amrex::Real sp1, + sedge1(const amrex::Real /*sm2*/, + const amrex::Real /*sm1*/, + const amrex::Real /*s0*/, + const amrex::Real /*sp1*/, const amrex::Real /*sp2*/) { return 0; @@ -384,10 +384,10 @@ struct minmod { [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real - sedge2(const amrex::Real sm2, - const amrex::Real sm1, - const amrex::Real s0, - const amrex::Real sp1, + sedge2(const amrex::Real /*sm2*/, + const amrex::Real /*sm1*/, + const amrex::Real /*s0*/, + const amrex::Real /*sp1*/, const amrex::Real /*sp2*/) { return 0; From 650d8d7aaa813f48a87e497c5706c180a5bb448b Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Fri, 16 Aug 2024 18:21:55 -0700 Subject: [PATCH 03/10] remove unused --- Godunov/hydro_godunov_ppm.H | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Godunov/hydro_godunov_ppm.H b/Godunov/hydro_godunov_ppm.H index f0f6d0a89..8b08442e1 100644 --- a/Godunov/hydro_godunov_ppm.H +++ b/Godunov/hydro_godunov_ppm.H @@ -856,7 +856,7 @@ void PredictStateOnXFace ( const int i, const int j, const int k, const int n, const amrex::Array4 &vel_edge, const amrex::BCRec bc, const int domlo, const int domhi, - const Limiter& limiter, + const Limiter& /*limiter*/, int limiter_type) { amrex::Real sm2 = S(i-2,j,k,n); @@ -930,7 +930,7 @@ void PredictStateOnYFace ( const int i, const int j, const int k, const int n, const amrex::Array4 &vel_edge, const amrex::BCRec bc, const int domlo, const int domhi, - const Limiter& limiter, + const Limiter& /*limiter*/, int limiter_type) { amrex::Real sm2 = S(i,j-2,k,n); @@ -1010,7 +1010,7 @@ void PredictStateOnZFace ( const int i, const int j, const int k, const int n, const amrex::Array4 &vel_edge, const amrex::BCRec bc, const int domlo, const int domhi, - const Limiter& limiter, + const Limiter& /*limiter*/, int limiter_type) { amrex::Real sm2 = S(i,j,k-2,n); From 12543e862e5183a5331f750e5a5338709c5deeca Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 17 Aug 2024 08:42:54 -0700 Subject: [PATCH 04/10] remove duplicate --- Godunov/hydro_godunov_ppm.H | 1 - 1 file changed, 1 deletion(-) diff --git a/Godunov/hydro_godunov_ppm.H b/Godunov/hydro_godunov_ppm.H index 8b08442e1..8ce638579 100644 --- a/Godunov/hydro_godunov_ppm.H +++ b/Godunov/hydro_godunov_ppm.H @@ -1050,7 +1050,6 @@ void PredictStateOnZFace ( const int i, const int j, const int k, const int n, } } else { - constexpr auto half{amrex::Real(0.5)}; constexpr auto two3rds{amrex::Real(2.0/3.0)}; amrex::Real sedge1 = Limiter::sedge1(sm2,sm1,s0,sp1,sp2); From bf298ff5e043914dd423893385987456488d570e Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 17 Aug 2024 09:15:41 -0700 Subject: [PATCH 05/10] fix typo --- Godunov/hydro_godunov_ppm.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Godunov/hydro_godunov_ppm.H b/Godunov/hydro_godunov_ppm.H index 8ce638579..958fb8dc2 100644 --- a/Godunov/hydro_godunov_ppm.H +++ b/Godunov/hydro_godunov_ppm.H @@ -266,7 +266,7 @@ struct wenoz { { return amrex::makeTuple(sedge1, sedge2); } -}; // wenoa +}; // wenoz struct weno_js { From 06a526aad6114f2839e0397e6dea2a7d42e153d7 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 17 Aug 2024 16:20:06 -0700 Subject: [PATCH 06/10] more fixes for amr-wind --- Utils/hydro_utils.H | 6 ++-- Utils/hydro_utils.cpp | 78 +++++++++++++++++++++++++++++++------------ 2 files changed, 60 insertions(+), 24 deletions(-) diff --git a/Utils/hydro_utils.H b/Utils/hydro_utils.H index 64704b935..122baaee1 100644 --- a/Utils/hydro_utils.H +++ b/Utils/hydro_utils.H @@ -317,7 +317,8 @@ void ComputeFluxes ( amrex::Box const& bx, amrex::Array4 const& yed, amrex::Array4 const& zed), amrex::Geometry const& geom, int ncomp, - bool fluxes_are_area_weighted); + bool fluxes_are_area_weighted, + int const* iconserv); /** * \brief Compute divergence. @@ -350,7 +351,8 @@ void EB_ComputeFluxes ( amrex::Box const& bx, amrex::Array4 const& apz), amrex::Geometry const& geom, int ncomp, amrex::Array4 const& flag, - bool fluxes_are_area_weighted); + bool fluxes_are_area_weighted, + int const* iconserv); void EB_ComputeDivergence ( amrex::Box const& bx, diff --git a/Utils/hydro_utils.cpp b/Utils/hydro_utils.cpp index b6b8cdb5c..49269efef 100644 --- a/Utils/hydro_utils.cpp +++ b/Utils/hydro_utils.cpp @@ -23,7 +23,8 @@ HydroUtils::ComputeFluxes ( Box const& bx, Array4 const& yed, Array4 const& zed), Geometry const& geom, int ncomp, - bool fluxes_are_area_weighted ) + bool fluxes_are_area_weighted, + int const* iconserv) { #if (AMREX_SPACEDIM == 2) if (geom.IsRZ()) { @@ -50,8 +51,10 @@ HydroUtils::ComputeFluxes ( Box const& bx, { if (fluxes_are_area_weighted) { fx(i,j,k,n) = xed(i,j,k,n) * umac(i,j,k) * ax(i,j,k); - } else { + } else if (iconserv[n] != -1) { fx(i,j,k,n) = xed(i,j,k,n) * umac(i,j,k); + } else { + fx(i,j,k,n) = xed(i,j,k,n); } }); @@ -64,8 +67,10 @@ HydroUtils::ComputeFluxes ( Box const& bx, { if (fluxes_are_area_weighted) { fy(i,j,k,n) = yed(i,j,k,n) * vmac(i,j,k) * ay(i,j,k); - } else { + } else if (iconserv[n] != -1) { fy(i,j,k,n) = yed(i,j,k,n) * vmac(i,j,k); + } else { + fy(i,j,k,n) = yed(i,j,k,n); } }); } else @@ -95,20 +100,28 @@ HydroUtils::ComputeFluxes ( Box const& bx, // X flux // const Box& xbx = amrex::surroundingNodes(bx,0); - amrex::ParallelFor(xbx, ncomp, [fx, umac, xed, area] + amrex::ParallelFor(xbx, ncomp, [fx, umac, xed, area, iconserv] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - fx(i,j,k,n) = xed(i,j,k,n) * umac(i,j,k) * area[0]; + if (iconserv[n] != -1) { + fx(i,j,k,n) = xed(i,j,k,n) * umac(i,j,k) * area[0]; + } else { + fx(i,j,k,n) = xed(i,j,k,n); + } }); // // Y flux // const Box& ybx = amrex::surroundingNodes(bx,1); - amrex::ParallelFor(ybx, ncomp, [fy, vmac, yed, area] + amrex::ParallelFor(ybx, ncomp, [fy, vmac, yed, area, iconserv] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - fy(i,j,k,n) = yed(i,j,k,n) * vmac(i,j,k) * area[1]; + if (iconserv[n] != -1) { + fy(i,j,k,n) = yed(i,j,k,n) * vmac(i,j,k) * area[1]; + } else { + fy(i,j,k,n) = yed(i,j,k,n); + } }); #if ( AMREX_SPACEDIM == 3 ) @@ -116,10 +129,14 @@ HydroUtils::ComputeFluxes ( Box const& bx, // Z flux // const Box& zbx = amrex::surroundingNodes(bx,2); - amrex::ParallelFor(zbx, ncomp, [fz, wmac, zed, area] + amrex::ParallelFor(zbx, ncomp, [fz, wmac, zed, area, iconserv] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - fz(i,j,k,n) = zed(i,j,k,n) * wmac(i,j,k) * area[2]; + if (iconserv[n] != -1) { + fz(i,j,k,n) = zed(i,j,k,n) * wmac(i,j,k) * area[2]; + } else { + fz(i,j,k,n) = zed(i,j,k,n); + } }); #endif } @@ -225,8 +242,9 @@ HydroUtils::ComputeConvectiveTerm(Box const& bx, int num_comp, amrex::ParallelFor(bx, num_comp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - if (!iconserv[n]) + if (iconserv[n] == 0) { convTerm(i,j,k,n) += q(i,j,k,n)*divu(i,j,k); + } }); } else if (advection_type == "Godunov" || advection_type == "BDS") @@ -242,7 +260,7 @@ HydroUtils::ComputeConvectiveTerm(Box const& bx, int num_comp, amrex::ParallelFor(bx, num_comp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - if (!iconserv[n]) + if (iconserv[n] == 0) { Real qavg = q_on_face_x(i,j,k,n) + q_on_face_x(i+1,j,k,n) + q_on_face_y(i,j,k,n) + q_on_face_y(i,j+1,k,n); @@ -268,7 +286,7 @@ HydroUtils::ComputeConvectiveTerm(Box const& bx, int num_comp, amrex::ParallelFor(bx, num_comp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - if (!iconserv[n] && vfrac_arr(i,j,k) > 0.) + if (iconserv[n] == 0 && vfrac_arr(i,j,k) > 0.) { Real qavg = apx_arr(i,j,k)*q_on_face_x(i,j,k,n) + apx_arr(i+1,j,k)*q_on_face_x(i+1,j,k,n); qavg += apy_arr(i,j,k)*q_on_face_y(i,j,k,n) + apy_arr(i,j+1,k)*q_on_face_y(i,j+1,k,n); @@ -414,7 +432,8 @@ HydroUtils::EB_ComputeFluxes ( Box const& bx, Array4 const& apz), Geometry const& geom, int ncomp, Array4 const& flag, - bool fluxes_are_area_weighted ) + bool fluxes_are_area_weighted, + int const* iconserv) { const auto dx = geom.CellSizeArray(); @@ -448,10 +467,15 @@ HydroUtils::EB_ComputeFluxes ( Box const& bx, amrex::ParallelFor(xbx, ncomp, [fx, umac, xed, area, apx, flag] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - if (flag(i,j,k).isConnected(-1,0,0)) - fx(i,j,k,n) = xed(i,j,k,n) * umac(i,j,k) * apx(i,j,k) * area[0]; - else + if (flag(i,j,k).isConnected(-1,0,0)) { + if (iconserv[n] == -1) { + fx(i,j,k,n) = xed(i,j,k,n); + } else { + fx(i,j,k,n) = xed(i,j,k,n) * umac(i,j,k) * apx(i,j,k) * area[0]; + } + } else { fx(i,j,k,n) = 0.; + } }); // @@ -462,10 +486,15 @@ HydroUtils::EB_ComputeFluxes ( Box const& bx, amrex::ParallelFor(ybx, ncomp, [fy, vmac, yed, area, apy, flag] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - if (flag(i,j,k).isConnected(0,-1,0)) - fy(i,j,k,n) = yed(i,j,k,n) * vmac(i,j,k) * apy(i,j,k) * area[1]; - else + if (flag(i,j,k).isConnected(0,-1,0)) { + if (iconserv[n] == -1) { + fy(i,j,k,n) = yed(i,j,k,n); + } else { + fy(i,j,k,n) = yed(i,j,k,n) * vmac(i,j,k) * apy(i,j,k) * area[1]; + } + } else { fy(i,j,k,n) = 0.; + } }); #if (AMREX_SPACEDIM==3) @@ -477,10 +506,15 @@ HydroUtils::EB_ComputeFluxes ( Box const& bx, amrex::ParallelFor(zbx, ncomp, [fz, wmac, zed, area, apz, flag] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - if (flag(i,j,k).isConnected(0,0,-1)) - fz(i,j,k,n) = zed(i,j,k,n) * wmac(i,j,k) * apz(i,j,k) * area[2]; - else + if (flag(i,j,k).isConnected(0,0,-1)) { + if (iconserv[n] == -1) { + fz(i,j,k,n) = zed(i,j,k,n); + } else { + fz(i,j,k,n) = zed(i,j,k,n) * wmac(i,j,k) * apz(i,j,k) * area[2]; + } + } else { fz(i,j,k,n) = 0.; + } }); #endif From 4f1fdc81c60d583cecb2de27dc665d58b60de261 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 17 Aug 2024 16:23:34 -0700 Subject: [PATCH 07/10] fix oops --- Utils/hydro_utils.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Utils/hydro_utils.cpp b/Utils/hydro_utils.cpp index 49269efef..968e04248 100644 --- a/Utils/hydro_utils.cpp +++ b/Utils/hydro_utils.cpp @@ -464,7 +464,7 @@ HydroUtils::EB_ComputeFluxes ( Box const& bx, // const Box& xbx = amrex::surroundingNodes(bx,0); - amrex::ParallelFor(xbx, ncomp, [fx, umac, xed, area, apx, flag] + amrex::ParallelFor(xbx, ncomp, [fx, umac, xed, area, apx, flag, iconserv] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { if (flag(i,j,k).isConnected(-1,0,0)) { @@ -483,7 +483,7 @@ HydroUtils::EB_ComputeFluxes ( Box const& bx, // const Box& ybx = amrex::surroundingNodes(bx,1); - amrex::ParallelFor(ybx, ncomp, [fy, vmac, yed, area, apy, flag] + amrex::ParallelFor(ybx, ncomp, [fy, vmac, yed, area, apy, flag, iconserv] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { if (flag(i,j,k).isConnected(0,-1,0)) { @@ -503,7 +503,7 @@ HydroUtils::EB_ComputeFluxes ( Box const& bx, // const Box& zbx = amrex::surroundingNodes(bx,2); - amrex::ParallelFor(zbx, ncomp, [fz, wmac, zed, area, apz, flag] + amrex::ParallelFor(zbx, ncomp, [fz, wmac, zed, area, apz, flag, iconserv] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { if (flag(i,j,k).isConnected(0,0,-1)) { From c735be90037b8b26c56b8516def5e76d08cc6d91 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 17 Aug 2024 16:33:32 -0700 Subject: [PATCH 08/10] fix another oops --- Utils/hydro_utils.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Utils/hydro_utils.cpp b/Utils/hydro_utils.cpp index 968e04248..afd362337 100644 --- a/Utils/hydro_utils.cpp +++ b/Utils/hydro_utils.cpp @@ -46,7 +46,7 @@ HydroUtils::ComputeFluxes ( Box const& bx, // X flux // const Box& xbx = amrex::surroundingNodes(bx,0); - amrex::ParallelFor(xbx, ncomp, [fx, umac, xed, ax, fluxes_are_area_weighted] + amrex::ParallelFor(xbx, ncomp, [fx, umac, xed, ax, fluxes_are_area_weightediconserv] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { if (fluxes_are_area_weighted) { @@ -62,7 +62,7 @@ HydroUtils::ComputeFluxes ( Box const& bx, // Y flux // const Box& ybx = amrex::surroundingNodes(bx,1); - amrex::ParallelFor(ybx, ncomp, [fy, vmac, yed, ay, fluxes_are_area_weighted] + amrex::ParallelFor(ybx, ncomp, [fy, vmac, yed, ay, fluxes_are_area_weighted, iconserv] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { if (fluxes_are_area_weighted) { From 09674a1ecfe06c7a789b6a798ee31e7b2f66142f Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 17 Aug 2024 16:37:49 -0700 Subject: [PATCH 09/10] fix typo --- Utils/hydro_utils.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Utils/hydro_utils.cpp b/Utils/hydro_utils.cpp index afd362337..ba52086d6 100644 --- a/Utils/hydro_utils.cpp +++ b/Utils/hydro_utils.cpp @@ -46,7 +46,7 @@ HydroUtils::ComputeFluxes ( Box const& bx, // X flux // const Box& xbx = amrex::surroundingNodes(bx,0); - amrex::ParallelFor(xbx, ncomp, [fx, umac, xed, ax, fluxes_are_area_weightediconserv] + amrex::ParallelFor(xbx, ncomp, [fx, umac, xed, ax, fluxes_are_area_weighted, iconserv] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { if (fluxes_are_area_weighted) { From fe0dbf10f654924a016fc596b34e14d554b38b75 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 17 Aug 2024 16:40:49 -0700 Subject: [PATCH 10/10] need to pass iconserv --- Utils/hydro_compute_edgestate_and_flux.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Utils/hydro_compute_edgestate_and_flux.cpp b/Utils/hydro_compute_edgestate_and_flux.cpp index 5d020d946..8992c4204 100644 --- a/Utils/hydro_compute_edgestate_and_flux.cpp +++ b/Utils/hydro_compute_edgestate_and_flux.cpp @@ -452,7 +452,8 @@ HydroUtils::ComputeFluxesOnBoxFromState (Box const& bx, int ncomp, MFIter& mfi, AMREX_D_DECL(face_x,face_y,face_z), AMREX_D_DECL(apx,apy,apz), geom, ncomp, - flag, fluxes_are_area_weighted); + flag, fluxes_are_area_weighted, + iconserv); } else #endif { @@ -460,7 +461,8 @@ HydroUtils::ComputeFluxesOnBoxFromState (Box const& bx, int ncomp, MFIter& mfi, AMREX_D_DECL(flux_x,flux_y,flux_z), AMREX_D_DECL(u_flux,v_flux,w_flux), AMREX_D_DECL(face_x,face_y,face_z), - geom, ncomp, fluxes_are_area_weighted ); + geom, ncomp, fluxes_are_area_weighted, + iconserv); } }