Skip to content

Commit

Permalink
generalize to allow for WENO_JS
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren committed Jul 29, 2024
1 parent a5d7828 commit 56d61f4
Show file tree
Hide file tree
Showing 6 changed files with 114 additions and 39 deletions.
6 changes: 6 additions & 0 deletions Godunov/hydro_godunov_edge_state_2D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,12 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
AMREX_D_DECL(Ipx,Ipy,Ipz),
AMREX_D_DECL(umac,vmac,wmac),
q,geom,l_dt,pbc,ncomp,limiter);
} 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 {
auto limiter = PPM::nolimiter();
PPM::PredictStateOnFaces(bxg1,AMREX_D_DECL(Imx,Imy,Imz),
Expand Down
6 changes: 6 additions & 0 deletions Godunov/hydro_godunov_edge_state_3D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,12 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
AMREX_D_DECL(Ipx,Ipy,Ipz),
AMREX_D_DECL(umac,vmac,wmac),
q,geom,l_dt,pbc,ncomp,limiter);
} 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 {
auto limiter = PPM::nolimiter();
PPM::PredictStateOnFaces(bxg1,AMREX_D_DECL(Imx,Imy,Imz),
Expand Down
6 changes: 6 additions & 0 deletions Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,12 @@ Godunov::ExtrapVelToFaces ( MultiFab const& a_vel,
Imx, Imy, Ipx, Ipy,
vel, vel,
geom, l_dt, d_bcrec, limiter);
} else if (limiter_type == PPM::WENO_JS) {
auto limiter = PPM::weno_js();
PPM::PredictVelOnFaces( bxg1,
Imx, Imy, Ipx, Ipy,
vel, vel,
geom, l_dt, d_bcrec, limiter);
} else {
auto limiter = PPM::nolimiter();
PPM::PredictVelOnFaces( bxg1,
Expand Down
6 changes: 6 additions & 0 deletions Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,12 @@ Godunov::ExtrapVelToFaces ( MultiFab const& a_vel,
Imx, Imy, Imz, Ipx, Ipy, Ipz,
vel, vel,
geom, l_dt, d_bcrec, limiter);
} else if (limiter_type == PPM::WENO_JS) {
auto limiter = PPM::weno_js();
PPM::PredictVelOnFaces( bxg1,
Imx, Imy, Imz, Ipx, Ipy, Ipz,
vel, vel,
geom, l_dt, d_bcrec, limiter);
} else {
auto limiter = PPM::nolimiter();
PPM::PredictVelOnFaces( bxg1,
Expand Down
125 changes: 88 additions & 37 deletions Godunov/hydro_godunov_ppm.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,12 +14,14 @@

namespace PPM {

enum limiters {VanLeer, WENOZ, NoLimiter};
enum limiters {VanLeer, WENOZ, WENO_JS, NoLimiter};

static constexpr int default_limiter = VanLeer;

struct nolimiter {

static constexpr bool do_limiting = false;

[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real
sedge1(const amrex::Real sm2,
Expand Down Expand Up @@ -64,6 +66,8 @@ struct nolimiter {

struct vanleer {

static constexpr bool do_limiting = true;

[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real vanLeer(const amrex::Real a, const amrex::Real b, const amrex::Real c)
{
Expand Down Expand Up @@ -151,6 +155,8 @@ struct vanleer {

struct wenoz {

static constexpr bool do_limiting = false;

[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real
sedge2(const amrex::Real sm2,
Expand Down Expand Up @@ -206,6 +212,69 @@ struct wenoz {
}
};

struct weno_js {

static constexpr bool do_limiting = false;

[[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)
{
constexpr auto eps = amrex::Real(1.0e-6);

const amrex::Real beta1 =
amrex::Real(13.0) / amrex::Real(12.0) * (sm2 - amrex::Real(2.0) * sm1 + s0) *
(sm2 - amrex::Real(2.0) * sm1 + s0) +
amrex::Real(0.25) * (sm2 - amrex::Real(4.0) * sm1 + amrex::Real(3.0) * s0) *
(sm2 - amrex::Real(4.0) * sm1 + amrex::Real(3.0) * s0);
const amrex::Real beta2 =
amrex::Real(13.0) / amrex::Real(12.0) * (sm1 - amrex::Real(2.0) * s0 + sp1) *
(sm1 - amrex::Real(2.0) * s0 + sp1) +
amrex::Real(0.25) * (sm1 - sp1) * (sm1 - sp1);
const amrex::Real beta3 =
amrex::Real(13.0) / amrex::Real(12.0) * (s0 - amrex::Real(2.0) * sp1 + sp2) *
(s0 - amrex::Real(2.0) * sp1 + sp2) +
amrex::Real(0.25) * (amrex::Real(3.0) * s0 - amrex::Real(4.0) * sp1 + sp2) *
(amrex::Real(3.0) * s0 - amrex::Real(4.0) * sp1 + sp2);

const amrex::Real t5 = amrex::Math::abs(beta3 - beta1);
const amrex::Real omega1 = amrex::Real(0.1) / (eps + beta1));
const amrex::Real omega2 = amrex::Real(0.6) / (eps + beta2));
const amrex::Real omega3 = amrex::Real(0.3) / (eps + beta3));

const amrex::Real omega = omega1 + omega2 + omega3;

const amrex::Real v_1 = amrex::Real(2.0) * sm2 - amrex::Real(7.0) * sm1 + amrex::Real(11.0) * s0;
const amrex::Real v_2 = -sm1 + amrex::Real(5.0) * s0 + amrex::Real(2.0) * sp1;
const amrex::Real v_3 = amrex::Real(2.0) * s0 + amrex::Real(5.0) * sp1 - sp2;

return (omega1 * v_1 + omega2 * v_2 + omega3 * v_3) / (amrex::Real(6.0) * omega);
}

[[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 sedge2(sp2,sp1,s0,sm1,sm2); // NOLINT(readability-suspicious-call-argument)
}

[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::GpuTuple<amrex::Real,amrex::Real>
sm_sp(const amrex::Real /*s0*/,
const amrex::Real sedge1,
const amrex::Real sedge2)
{
return amrex::makeTuple(sedge1, sedge2);
}
};

template <typename Limiter>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Expand All @@ -224,8 +293,7 @@ void SetXBCs ( const int i, const int j, const int k, const int n,
{
sedge2 = -amrex::Real(0.2)*s(domlo-1,j,k,n) + amrex::Real(0.75)*s(domlo,j, k, n)
+amrex::Real(0.5)*s(domlo+1,j,k,n) - amrex::Real(0.05)*s(domlo+2,j,k,n);
if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
sedge2 = amrex::max(sedge2, amrex::min(s(domlo+1,j,k,n), s(domlo,j,k,n)));
sedge2 = amrex::min(sedge2, amrex::max(s(domlo+1,j,k,n), s(domlo,j,k,n)));
}
Expand All @@ -237,17 +305,15 @@ void SetXBCs ( const int i, const int j, const int k, const int n,

sedge1 = -amrex::Real(0.2)*s(domlo-1,j,k,n) + amrex::Real(0.75)*s(domlo ,j,k,n)
+ amrex::Real(0.5)*s(domlo+1,j,k,n) - amrex::Real(0.05)*s(domlo+2,j,k,n);
if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
sedge1 = amrex::max(sedge1, amrex::min(s(domlo+1,j,k,n), s(domlo,j,k,n)));
sedge1 = amrex::min(sedge1, amrex::max(s(domlo+1,j,k,n), s(domlo,j,k,n)));
}

sp = sedge2;
sm = sedge1;

if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
if ( (sp - s(domlo+1,j,k,n))*(s(domlo+1,j,k,n) - sm) <= amrex::Real(0.0))
{
sp = s(domlo+1,j,k,n);
Expand All @@ -267,8 +333,7 @@ void SetXBCs ( const int i, const int j, const int k, const int n,
{
sedge1 = -amrex::Real(0.2)*s(domhi+1,j,k,n) + amrex::Real(0.75)*s(domhi,j,k, n)
+amrex::Real(0.5)*s(domhi-1,j,k,n) - amrex::Real(0.05)*s(domhi-2,j,k,n);
if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
sedge1 = amrex::max(sedge1, amrex::min(s(domhi-1,j,k,n), s(domhi,j,k,n)));
sedge1 = amrex::min(sedge1, amrex::max(s(domhi-1,j,k,n), s(domhi,j,k,n)));
}
Expand All @@ -280,17 +345,15 @@ void SetXBCs ( const int i, const int j, const int k, const int n,

sedge2 = -amrex::Real(0.2)*s(domhi+1,j,k,n) + amrex::Real(0.75)*s(domhi ,j,k,n)
+amrex::Real(0.5)*s(domhi-1,j,k,n) - amrex::Real(0.05)*s(domhi-2,j,k,n);
if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
sedge2 = amrex::max(sedge2, amrex::min(s(domhi-1,j,k,n), s(domhi,j,k,n)));
sedge2 = amrex::min(sedge2, amrex::max(s(domhi-1,j,k,n), s(domhi,j,k,n)));
}

sp = sedge2;
sm = sedge1;

if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
if( (sp - s(domhi-1,j,k,n))*(s(domhi-1,j,k,n) - sm) <= amrex::Real(0.0))
{
sp = s(domhi-1,j,k,n);
Expand Down Expand Up @@ -322,8 +385,7 @@ void SetYBCs ( const int i, const int j, const int k, const int n,
{
sedge2 = -amrex::Real(0.2)*s(i,domlo-1,k,n) + amrex::Real(0.75)*s(i,domlo ,k,n)
+amrex::Real(0.5)*s(i,domlo+1,k,n) - amrex::Real(0.05)*s(i,domlo+2,k,n);
if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
sedge2 = amrex::max(sedge2, amrex::min(s(i,domlo+1,k,n), s(i,domlo,k,n)));
sedge2 = amrex::min(sedge2, amrex::max(s(i,domlo+1,k,n), s(i,domlo,k,n)));
}
Expand All @@ -335,17 +397,15 @@ void SetYBCs ( const int i, const int j, const int k, const int n,

sedge1 = -amrex::Real(0.2)*s(i,domlo-1,k,n) + amrex::Real(0.75)*s(i,domlo ,k,n)
+amrex::Real(0.5)*s(i,domlo+1,k,n) - amrex::Real(0.05)*s(i,domlo+2,k,n);
if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
sedge1 = amrex::max(sedge1, amrex::min(s(i,domlo+1,k,n), s(i,domlo,k,n)));
sedge1 = amrex::min(sedge1, amrex::max(s(i,domlo+1,k,n), s(i,domlo,k,n)));
}

sp = sedge2;
sm = sedge1;

if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
if ( (sp - s(i,domlo+1,k,n))*(s(i,domlo+1,k,n) - sm) <= amrex::Real(0.0))
{
sp = s(i,domlo+1,k,n);
Expand All @@ -365,8 +425,7 @@ void SetYBCs ( const int i, const int j, const int k, const int n,
{
sedge1 = -amrex::Real(0.2)*s(i,domhi+1,k,n) + amrex::Real(0.75)*s(i,domhi ,k,n)
+amrex::Real(0.5)*s(i,domhi-1,k,n) - amrex::Real(0.05)*s(i,domhi-2,k,n);
if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
sedge1 = amrex::max(sedge1, amrex::min(s(i,domhi-1,k,n), s(i,domhi,k,n)));
sedge1 = amrex::min(sedge1, amrex::max(s(i,domhi-1,k,n), s(i,domhi,k,n)));
}
Expand All @@ -378,17 +437,15 @@ void SetYBCs ( const int i, const int j, const int k, const int n,

sedge2 = -amrex::Real(0.2)*s(i,domhi+1,k,n) + amrex::Real(0.75)*s(i,domhi ,k,n)
+amrex::Real(0.5)*s(i,domhi-1,k,n) - amrex::Real(0.05)*s(i,domhi-2,k,n);
if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
sedge2 = amrex::max(sedge2, amrex::min(s(i,domhi-1,k,n), s(i,domhi,k,n)));
sedge2 = amrex::min(sedge2, amrex::max(s(i,domhi-1,k,n), s(i,domhi,k,n)));
}

sp = sedge2;
sm = sedge1;

if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
if( (sp - s(i,domhi-1,k,n))*(s(i,domhi-1,k,n) - sm) <= amrex::Real(0.0)){
sp = s(i,domhi-1,k,n);
sm = s(i,domhi-1,k,n);
Expand Down Expand Up @@ -420,8 +477,7 @@ void SetZBCs ( const int i, const int j, const int k, const int n,
{
sedge2 = -amrex::Real(0.2)*s(i,j,domlo-1,n) + amrex::Real(0.75)*s(i,j,domlo ,n)
+amrex::Real(0.5)*s(i,j,domlo+1,n) - amrex::Real(0.05)*s(i,j,domlo+2,n);
if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
sedge2 = amrex::max(sedge2, amrex::min(s(i,j,domlo+1,n), s(i,j,domlo,n)));
sedge2 = amrex::min(sedge2, amrex::max(s(i,j,domlo+1,n), s(i,j,domlo,n)));
}
Expand All @@ -433,17 +489,15 @@ void SetZBCs ( const int i, const int j, const int k, const int n,

sedge1 = -amrex::Real(0.2)*s(i,j,domlo-1,n) + amrex::Real(0.75)*s(i,j,domlo ,n)
+amrex::Real(0.5)*s(i,j,domlo+1,n) - amrex::Real(0.05)*s(i,j,domlo+2,n);
if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
sedge1 = amrex::max(sedge1, amrex::min(s(i,j,domlo+1,n), s(i,j,domlo,n)));
sedge1 = amrex::min(sedge1, amrex::max(s(i,j,domlo+1,n), s(i,j,domlo,n)));
}

sp = sedge2;
sm = sedge1;

if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
if ( (sp - s(i,j,domlo+1,n))*(s(i,j,domlo+1,n) - sm) <= 0. )
{
sp = s(i,j,domlo+1,n);
Expand All @@ -463,8 +517,7 @@ void SetZBCs ( const int i, const int j, const int k, const int n,
{
sedge1 = -amrex::Real(0.2)*s(i,j,domhi+1,n) + amrex::Real(0.75)*s(i,j,domhi ,n)
+amrex::Real(0.5)*s(i,j,domhi-1,n) - amrex::Real(0.05)*s(i,j,domhi-2,n);
if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
sedge1 = amrex::max(sedge1, amrex::min(s(i,j,domhi-1,n), s(i,j,domhi,n)));
sedge1 = amrex::min(sedge1, amrex::max(s(i,j,domhi-1,n), s(i,j,domhi,n)));
}
Expand All @@ -476,17 +529,15 @@ void SetZBCs ( const int i, const int j, const int k, const int n,

sedge2 = -amrex::Real(0.2)*s(i,j,domhi+1,n) + amrex::Real(0.75)*s(i,j,domhi ,n)
+amrex::Real(0.5)*s(i,j,domhi-1,n) - amrex::Real(0.05)*s(i,j,domhi-2,n);
if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
sedge2 = amrex::max(sedge2, amrex::min(s(i,j,domhi-1,n), s(i,j,domhi,n)));
sedge2 = amrex::min(sedge2, amrex::max(s(i,j,domhi-1,n), s(i,j,domhi,n)));
}

sp = sedge2;
sm = sedge1;

if constexpr ( (! std::is_same_v<nolimiter,Limiter>) &&
(! std::is_same_v<wenoz,Limiter>) ) {
if constexpr (Limiter::do_limiting) {
if ( (sp - s(i,j,domhi-1,n))*(s(i,j,domhi-1,n) - sm) <= 0. )
{
sp = s(i,j,domhi-1,n);
Expand Down
4 changes: 2 additions & 2 deletions Utils/hydro_utils.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,8 +130,8 @@ ComputeFluxesOnBoxFromState ( amrex::Box const& bx, int ncomp, amrex::MFIter& mf
bool godunov_use_ppm, bool godunov_use_forces_in_trans,
bool is_velocity, bool fluxes_are_area_weighted,
std::string const& advection_type,
int limiter_type = PPM::default_limiter,
amrex::Array4<int const> const& bc_arr = {});
int limiter_type = PPM::default_limiter,
amrex::Array4<int const> const& bc_arr = {});
#endif

#ifdef AMREX_USE_EB
Expand Down

0 comments on commit 56d61f4

Please sign in to comment.