Skip to content

Commit

Permalink
Merge pull request #685 from Pressio/680-removesimplify-constants
Browse files Browse the repository at this point in the history 
#680: Remove/simplify constants
  • Loading branch information
@fnrizzi
fnrizzi authored Sep 26, 2024
2 parents fd938ba + b577bfd commit 825d617
Show file tree
Hide file tree
Showing 61 changed files with 281 additions and 395 deletions.
2 changes: 1 addition & 1 deletion include/pressio/ode/impl/ode_advance_to_target_time.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,7 +135,7 @@ void to_target_time_with_step_size_policy(StepperType & stepper,
}

if (enableTimeStepRecovery){
if (dtScalingFactor.get() <= ::pressio::utils::Constants<IndVarType>::one()){
if (dtScalingFactor.get() <= static_cast<IndVarType>(1)){
// need to change this to use some notion of identity
throw std::runtime_error("The time step size reduction factor must be > 1.");
}
Expand Down
46 changes: 18 additions & 28 deletions include/pressio/ode/impl/ode_explicit_stepper_with_mass_matrix.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -230,7 +230,7 @@ class ExplicitStepperWithMassMatrixImpl
// need to do: y_n+1 = y_n + stepSize*x
// odeState already contains y_n
using scalar_type = typename ::pressio::Traits<StateType>::scalar_type;
constexpr auto one = ::pressio::utils::Constants<scalar_type>::one();
constexpr auto one = static_cast<scalar_type>(1);
::pressio::ops::update(odeState, one, x, stepSize);
}

Expand All @@ -249,7 +249,7 @@ class ExplicitStepperWithMassMatrixImpl
*/

using scalar_type = typename ::pressio::Traits<StateType>::scalar_type;
constexpr auto one = ::pressio::utils::Constants<scalar_type>::one();
constexpr auto one = static_cast<scalar_type>(1);

auto & fn = rhsInstance_;

Expand All @@ -273,8 +273,10 @@ class ExplicitStepperWithMassMatrixImpl
rhsObserver(stepNumber, ::pressio::ode::IntermediateStepCount(0), stepStartVal, fn);
solver.solve(massMatrix_, xn, fn);

const auto cfn = ::pressio::utils::Constants<scalar_type>::threeOvTwo()*stepSize;
const auto cfnm1 = ::pressio::utils::Constants<scalar_type>::negOneHalf()*stepSize;
using cnst = ::pressio::ode::constants::Constants<scalar_type>;
const auto cfn = cnst::threeOvTwo()*stepSize;
const auto cfnm1 = cnst::negOneHalf()*stepSize;

::pressio::ops::update(odeState, one, xn, cfn, xnm1, cfnm1);
}

Expand All @@ -292,23 +294,15 @@ class ExplicitStepperWithMassMatrixImpl
PRESSIOLOG_DEBUG("ssprk3 stepper: do step");

using scalar_type = typename ::pressio::Traits<StateType>::scalar_type;
constexpr auto zero = ::pressio::utils::Constants<scalar_type>::zero();
constexpr auto one = ::pressio::utils::Constants<scalar_type>::one();
constexpr auto two = ::pressio::utils::Constants<scalar_type>::two();
constexpr auto three = ::pressio::utils::Constants<scalar_type>::three();
constexpr auto four = ::pressio::utils::Constants<scalar_type>::four();
constexpr auto oneOvThree = one/three;
constexpr auto twoOvThree = two/three;
constexpr auto threeOvFour = three/four;
constexpr auto fourInv = one/four;
using cnst = ::pressio::ode::constants::Constants<scalar_type>;

auto & rhs = rhsInstance_;
auto & x = xInstances_[0];
auto & auxState = xInstances_[1];

// see e.g. https://gkeyll.readthedocs.io/en/latest/dev/ssp-rk.html

const scalar_type stepSize_half{stepSize/two};
const scalar_type stepSize_half{stepSize/cnst::two()};
const independent_variable_type t_phalf{stepStartTime + stepSize_half};
const independent_variable_type t_next{stepStartTime + stepSize};

Expand All @@ -317,21 +311,21 @@ class ExplicitStepperWithMassMatrixImpl
rhsObserver(stepNumber, ::pressio::ode::IntermediateStepCount(0), stepStartTime, rhs);
solver.solve(massMatrix_, x, rhs);
// u_1 = u_n + stepSize * x
::pressio::ops::update(auxState, zero, odeState, one, x, stepSize);
::pressio::ops::update(auxState, cnst::zero(), odeState, cnst::one(), x, stepSize);

// rhs(u_1, t_n+stepSize)
systemObj_.get().massMatrixAndRhs(auxState, t_next, massMatrix_, rhs);
rhsObserver(stepNumber, ::pressio::ode::IntermediateStepCount(1), t_next, rhs);
solver.solve(massMatrix_, x, rhs);
// u_2 = 3/4*u_n + 1/4*u_1 + 1/4*stepSize*x
::pressio::ops::update(auxState, fourInv, odeState, threeOvFour, x, fourInv*stepSize);
::pressio::ops::update(auxState, cnst::fourInv(), odeState, cnst::threeOvFour(), x, cnst::fourInv()*stepSize);

// rhs(u_2, t_n + 0.5*stepSize)
systemObj_.get().massMatrixAndRhs(auxState, t_phalf, massMatrix_, rhs);
rhsObserver(stepNumber, ::pressio::ode::IntermediateStepCount(2), t_phalf, rhs);
solver.solve(massMatrix_, x, rhs);
// u_n+1 = 1/3*u_n + 2/3*u_2 + 2/3*stepSize*rhs(u_2, t_n+0.5*stepSize)
::pressio::ops::update(odeState, oneOvThree, auxState, twoOvThree, x, twoOvThree*stepSize);
::pressio::ops::update(odeState, cnst::oneOvThree(), auxState, cnst::twoOvThree(), x, cnst::twoOvThree()*stepSize);
}

template<class LinearSolver, class RhsObserverType>
Expand All @@ -354,16 +348,13 @@ class ExplicitStepperWithMassMatrixImpl
auto & auxState = xInstances_[4];

using scalar_type = typename ::pressio::Traits<StateType>::scalar_type;
constexpr auto one = ::pressio::utils::Constants<scalar_type>::one();
constexpr auto two = ::pressio::utils::Constants<scalar_type>::two();
constexpr auto three = ::pressio::utils::Constants<scalar_type>::three();
constexpr auto six = two * three;
using cnst = ::pressio::ode::constants::Constants<scalar_type>;

const scalar_type stepSize_half{stepSize / two};
const scalar_type stepSize_half{stepSize / cnst::two()};
const independent_variable_type t_phalf{stepStartTime + stepSize_half};
const independent_variable_type t_next{stepStartTime + stepSize};
const scalar_type stepSize6{stepSize / six};
const scalar_type stepSize3{stepSize / three};
const scalar_type stepSize6{stepSize / cnst::six()};
const scalar_type stepSize3{stepSize / cnst::three()};

// stage 1:
// rhs1 = rhs(y_n, t_n)
Expand Down Expand Up @@ -395,7 +386,7 @@ class ExplicitStepperWithMassMatrixImpl
rhsObserver(stepNumber, ::pressio::ode::IntermediateStepCount(3), t_next, rhs);
solver.solve(massMatrix_, x4, rhs);

::pressio::ops::update(odeState, one,
::pressio::ops::update(odeState, cnst::one(),
x1, stepSize6, x2, stepSize3,
x3, stepSize3, x4, stepSize6);
}
Expand All @@ -407,9 +398,8 @@ class ExplicitStepperWithMassMatrixImpl
const FactorType & rhsFactor)
{
using scalar_type = typename ::pressio::Traits<StateType>::scalar_type;
constexpr auto zero = ::pressio::utils::Constants<scalar_type>::zero();
constexpr auto one = ::pressio::utils::Constants<scalar_type>::one();
::pressio::ops::update(yIn, zero, stateIn, one, xIn, rhsFactor);
using cnst = ::pressio::ode::constants::Constants<scalar_type>;
::pressio::ops::update(yIn, cnst::zero(), stateIn, cnst::one(), xIn, rhsFactor);
}

};
Expand Down
65 changes: 25 additions & 40 deletions include/pressio/ode/impl/ode_explicit_stepper_without_mass_matrix.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -181,7 +181,7 @@ class ExplicitStepperNoMassMatrixImpl{

// y = y + stepSize * rhs
using scalar_type = typename ::pressio::Traits<StateType>::scalar_type;
constexpr auto one = ::pressio::utils::Constants<scalar_type>::one();
constexpr auto one = static_cast<scalar_type>(1);
::pressio::ops::update(odeState, one, rhs, stepSize);
}

Expand All @@ -198,18 +198,16 @@ class ExplicitStepperNoMassMatrixImpl{
// // y_n+1 = y_n + stepSize*[ (3/2)*f(y_n, t_n) - (1/2)*f(y_n-1, t_n-1) ]

using scalar_type = typename ::pressio::Traits<StateType>::scalar_type;
const auto cfn = ::pressio::utils::Constants<scalar_type>::threeOvTwo()*stepSize;
const auto cfnm1 = ::pressio::utils::Constants<scalar_type>::negOneHalf()*stepSize;
using cnst = ::pressio::ode::constants::Constants<scalar_type>;
const auto cfn = cnst::threeOvTwo()*stepSize;
const auto cfnm1 = cnst::negOneHalf()*stepSize;

if (stepNumber.get()==1){
// use Euler forward or we could use something else here maybe RK4
auto & rhs = rhsInstances_[0];
systemObj_.get().rhs(odeState, stepStartTime, rhs);
rhsObserver(stepNumber, ::pressio::ode::IntermediateStepCount(0), stepStartTime, rhs);

using scalar_type = typename ::pressio::Traits<StateType>::scalar_type;
constexpr auto one = ::pressio::utils::Constants<scalar_type>::one();
::pressio::ops::update(odeState, one, rhs, stepSize);
::pressio::ops::update(odeState, cnst::one(), rhs, stepSize);
}
else{
auto & fn = rhsInstances_[0];
Expand All @@ -219,10 +217,7 @@ class ExplicitStepperNoMassMatrixImpl{

systemObj_.get().rhs(odeState, stepStartTime, fn);
rhsObserver(stepNumber, ::pressio::ode::IntermediateStepCount(0), stepStartTime, fn);

using scalar_type = typename ::pressio::Traits<StateType>::scalar_type;
constexpr auto one = ::pressio::utils::Constants<scalar_type>::one();
::pressio::ops::update(odeState, one, fn, cfn, fnm1, cfnm1);
::pressio::ops::update(odeState, cnst::one(), fn, cfn, fnm1, cfnm1);
}
}

Expand All @@ -237,47 +232,41 @@ class ExplicitStepperNoMassMatrixImpl{
PRESSIOLOG_DEBUG("ssprk3 stepper: do step");

using scalar_type = typename ::pressio::Traits<StateType>::scalar_type;
constexpr auto zero = ::pressio::utils::Constants<scalar_type>::zero();
constexpr auto one = ::pressio::utils::Constants<scalar_type>::one();
constexpr auto two = ::pressio::utils::Constants<scalar_type>::two();
constexpr auto three = ::pressio::utils::Constants<scalar_type>::three();
constexpr auto four = ::pressio::utils::Constants<scalar_type>::four();
constexpr auto oneOvThree = one/three;
constexpr auto twoOvThree = two/three;
constexpr auto threeOvFour = three/four;
constexpr auto fourInv = one/four;
using cnst = ::pressio::ode::constants::Constants<scalar_type>;

auto & rhs0 = rhsInstances_[0];

// see e.g. https://gkeyll.readthedocs.io/en/latest/dev/ssp-rk.html

const scalar_type stepSize_half{stepSize/two};
const scalar_type stepSize_half{stepSize/cnst::two()};
const independent_variable_type t_phalf{stepStartTime + stepSize_half};
const independent_variable_type t_next{stepStartTime + stepSize};

// rhs(u_n, t_n)
systemObj_.get().rhs(odeState, stepStartTime, rhs0);
rhsObserver(stepNumber, ::pressio::ode::IntermediateStepCount(0), stepStartTime, rhs0);
// u_1 = u_n + stepSize * rhs(u_n, t_n)
::pressio::ops::update(auxiliaryState_, zero,
odeState, one,
::pressio::ops::update(auxiliaryState_, cnst::zero(),
odeState, cnst::one(),
rhs0, stepSize);

// rhs(u_1, t_n+stepSize)
systemObj_.get().rhs(auxiliaryState_, t_next, rhs0);
rhsObserver(stepNumber, ::pressio::ode::IntermediateStepCount(1), t_next, rhs0);
// u_2 = 3/4*u_n + 1/4*u_1 + 1/4*stepSize*rhs(u_1, t_n+stepSize)
::pressio::ops::update(auxiliaryState_, fourInv,
odeState, threeOvFour,
rhs0, fourInv*stepSize);
::pressio::ops::update(
auxiliaryState_, cnst::fourInv(), odeState,
cnst::threeOvFour(), rhs0, cnst::fourInv()*stepSize
);

// rhs(u_2, t_n + 0.5*stepSize)
systemObj_.get().rhs(auxiliaryState_, t_phalf, rhs0);
rhsObserver(stepNumber, ::pressio::ode::IntermediateStepCount(2), t_phalf, rhs0);
// u_n+1 = 1/3*u_n + 2/3*u_2 + 2/3*stepSize*rhs(u_2, t_n+0.5*stepSize)
::pressio::ops::update(odeState, oneOvThree,
auxiliaryState_, twoOvThree,
rhs0, twoOvThree*stepSize);
::pressio::ops::update(
odeState, cnst::oneOvThree(), auxiliaryState_,
cnst::twoOvThree(), rhs0, cnst::twoOvThree()*stepSize
);
}

template<class RhsObserverType>
Expand All @@ -296,16 +285,13 @@ class ExplicitStepperNoMassMatrixImpl{
auto & rhs4 = rhsInstances_[3];

using scalar_type = typename ::pressio::Traits<StateType>::scalar_type;
constexpr auto one = ::pressio::utils::Constants<scalar_type>::one();
constexpr auto two = ::pressio::utils::Constants<scalar_type>::two();
constexpr auto three = ::pressio::utils::Constants<scalar_type>::three();
constexpr auto six = two * three;
using cnst = ::pressio::ode::constants::Constants<scalar_type>;

const scalar_type stepSize_half{stepSize/two};
const scalar_type stepSize_half{stepSize/cnst::two()};
const independent_variable_type t_phalf{stepStartTime + stepSize_half};
const independent_variable_type t_next{stepStartTime + stepSize};
const scalar_type stepSize6{stepSize/six};
const scalar_type stepSize3{stepSize/three};
const scalar_type stepSize6{stepSize/cnst::six()};
const scalar_type stepSize3{stepSize/cnst::three()};

// stage 1:
// rhs1 = rhs(y_n, t_n)
Expand Down Expand Up @@ -334,7 +320,7 @@ class ExplicitStepperNoMassMatrixImpl{
rhsObserver(stepNumber, ::pressio::ode::IntermediateStepCount(3), t_next, rhs4);

// y_n += stepSize/6 * ( rhs1 + 2*rhs2 + 2*rhs3 + rhs4 )
::pressio::ops::update(odeState, one,
::pressio::ops::update(odeState, cnst::one(),
rhs1, stepSize6, rhs2, stepSize3,
rhs3, stepSize3, rhs4, stepSize6);
}
Expand All @@ -346,9 +332,8 @@ class ExplicitStepperNoMassMatrixImpl{
const FactorType & rhsFactor)
{
using scalar_type = typename ::pressio::Traits<StateType>::scalar_type;
constexpr auto zero = ::pressio::utils::Constants<scalar_type>::zero();
constexpr auto one = ::pressio::utils::Constants<scalar_type>::one();
::pressio::ops::update(yIn, zero, stateIn, one, rhsIn, rhsFactor);
using cnst = ::pressio::ode::constants::Constants<scalar_type>;
::pressio::ops::update(yIn, cnst::zero(), stateIn, cnst::one(), rhsIn, rhsFactor);
}

};
Expand Down
2 changes: 1 addition & 1 deletion include/pressio/ode/impl/ode_implicit_discrete_jacobian.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,7 @@ discrete_jacobian(::pressio::ode::BDF2,
{

using sc_t = typename ::pressio::Traits<JacobianType>::scalar_type;
// constexpr sc_t one = ::pressio::utils::Constants<sc_t>::one();
// constexpr sc_t one = static_cast<scalar_type>(1);
constexpr sc_t cnp1 = ::pressio::ode::constants::bdf2<sc_t>::c_np1_;
const sc_t cf = ::pressio::ode::constants::bdf2<sc_t>::c_f_ * dt;
::pressio::ops::update(jac, cf, M_np1, cnp1);
Expand Down
17 changes: 8 additions & 9 deletions include/pressio/ode/impl/ode_implicit_discrete_residual.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,8 +117,8 @@ discrete_residual(::pressio::ode::BDF1,
{

using sc_t = typename ::pressio::Traits<ResidualType>::scalar_type;
constexpr sc_t zero = ::pressio::utils::Constants<sc_t>::zero();
constexpr sc_t one = ::pressio::utils::Constants<sc_t>::one();
constexpr sc_t zero = static_cast<sc_t>(0);
constexpr sc_t one = static_cast<sc_t>(1);

const sc_t cf = ::pressio::ode::constants::bdf1<sc_t>::c_f_ * dt;
const auto & y_n = stencilStates(::pressio::ode::n());
Expand Down Expand Up @@ -200,8 +200,8 @@ discrete_residual(::pressio::ode::BDF2,
{

using sc_t = typename ::pressio::Traits<ResidualType>::scalar_type;
constexpr sc_t zero = ::pressio::utils::Constants<sc_t>::zero();
constexpr sc_t one = ::pressio::utils::Constants<sc_t>::one();
using cnst = ::pressio::ode::constants::Constants<sc_t>;

constexpr sc_t cnp1 = ::pressio::ode::constants::bdf2<sc_t>::c_np1_;
constexpr sc_t cn = ::pressio::ode::constants::bdf2<sc_t>::c_n_;
constexpr sc_t cnm1 = ::pressio::ode::constants::bdf2<sc_t>::c_nm1_;
Expand All @@ -211,11 +211,11 @@ discrete_residual(::pressio::ode::BDF2,
const auto & y_nm1 = stencilStates(::pressio::ode::nMinusOne());

// scratchState = y_n+1 - (4/3)*y_n + (1/3)*y_n-1
::pressio::ops::update(scratchState, zero, y_np1, cnp1, y_n, cn, y_nm1, cnm1);
::pressio::ops::update(scratchState, cnst::zero(), y_np1, cnp1, y_n, cn, y_nm1, cnm1);
// R = M_n+1 * scratchState
::pressio::ops::product(::pressio::nontranspose(), one, M_np1, scratchState, zero, R);
::pressio::ops::product(::pressio::nontranspose(), cnst::one(), M_np1, scratchState, cnst::zero(), R);
// R = R - dt * f_n+1
::pressio::ops::update(R, one, f_np1, cf);
::pressio::ops::update(R, cnst::one(), f_np1, cf);
}

/*
Expand Down Expand Up @@ -249,7 +249,6 @@ discrete_residual(::pressio::ode::CrankNicolson,

using sc_t = typename ::pressio::Traits<ResidualType>::scalar_type;

constexpr auto zero = ::pressio::utils::Constants<sc_t>::zero();
using cnst = ::pressio::ode::constants::cranknicolson<sc_t>;
constexpr sc_t cnp1 = cnst::c_np1_;
constexpr sc_t cn = cnst::c_n_;
Expand All @@ -261,7 +260,7 @@ discrete_residual(::pressio::ode::CrankNicolson,
const auto & f_n = stencilVelocities(::pressio::ode::n());
const auto & f_np1 = stencilVelocities(::pressio::ode::nPlusOne());

::pressio::ops::update(R, zero,
::pressio::ops::update(R, static_cast<sc_t>(0),
y_np1, cnp1,
y_n, cn,
f_n, cfnDt,
Expand Down
26 changes: 23 additions & 3 deletions include/pressio/ode/ode_constants.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,17 +55,37 @@ constexpr typename StepCount::value_type first_step_value = 1;

namespace constants{

template <typename scalar_t>
struct Constants
{
static constexpr scalar_t negOne(){ return static_cast<scalar_t>(-1); }
static constexpr scalar_t zero() { return static_cast<scalar_t>(0); }
static constexpr scalar_t one() { return static_cast<scalar_t>(1); }
static constexpr scalar_t two() { return static_cast<scalar_t>(2); }
static constexpr scalar_t three() { return static_cast<scalar_t>(3); }
static constexpr scalar_t four() { return static_cast<scalar_t>(4); }
static constexpr scalar_t six() { return static_cast<scalar_t>(6); }

static constexpr scalar_t negOneHalf() { return negOne()/two(); }
static constexpr scalar_t oneOvThree() { return one()/three(); }
static constexpr scalar_t twoOvThree() { return two()/three(); }
static constexpr scalar_t threeOvFour() { return three()/four(); }
static constexpr scalar_t fourOvThree() { return four()/three(); }
static constexpr scalar_t threeOvTwo() { return three()/two(); }
static constexpr scalar_t fourInv() { return one()/four(); }
};

template <typename scalar_t>
struct bdf1{
using cnst = ::pressio::utils::Constants<scalar_t>;
using cnst = Constants<scalar_t>;
static constexpr scalar_t c_np1_= cnst::one();
static constexpr scalar_t c_n_ = cnst::negOne();
static constexpr scalar_t c_f_ = cnst::negOne();
};

template <typename scalar_t>
struct bdf2{
using cnst = ::pressio::utils::Constants<scalar_t>;
using cnst = Constants<scalar_t>;
static constexpr scalar_t c_np1_ = cnst::one();
static constexpr scalar_t c_n_ = cnst::negOne()*cnst::fourOvThree();
static constexpr scalar_t c_nm1_ = cnst::oneOvThree();
Expand All @@ -74,7 +94,7 @@ struct bdf2{

template <typename scalar_t>
struct cranknicolson{
using cnst = ::pressio::utils::Constants<scalar_t>;
using cnst = Constants<scalar_t>;
static constexpr scalar_t c_np1_ = cnst::one();
static constexpr scalar_t c_n_ = cnst::negOne();
static constexpr scalar_t c_fnp1_ = cnst::negOneHalf();
Expand Down
Loading

0 comments on commit 825d617

Please sign in to comment.