Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

#680: Remove/simplify constants #685

Merged
merged 4 commits into from
Sep 26, 2024
Merged

#680: Remove/simplify constants #685

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f0076f7
#680: move constants to ode; static_cast in place when appropriate
cwschilly Sep 23, 2024
aba8bd7
#680: fix comment
cwschilly Sep 25, 2024
ae2d28e
#680: test with wrong constants
cwschilly Sep 25, 2024
b577bfd
Revert "#680: test with wrong constants"
cwschilly Sep 25, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion include/pressio/ode/impl/ode_advance_to_target_time.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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