Explicit Time Integration for Arbitrary Types

Custom Vector Class

Suppose that you have an application that uses a custom vector class as follows:

template<class ScalarType>
struct CustomVector
{
  using value_type = ScalarType;
  using size_type = std::size_t;

  CustomVector(std::size_t ext) : d_(ext){}

  ScalarType & operator()(int i){ return d_[i]; }
  const ScalarType & operator()(int i)const { return d_[i]; }
  ScalarType & operator[](int i){ return d_[i]; }
  const ScalarType & operator[](int i)const { return d_[i]; }

  std::size_t extent(int k)const { return (k==0) ? d_.size() : 0; }

  void fill(ScalarType value){
    std::for_each(d_.begin(), d_.end(), [](ScalarType & v){ v= 0.; });
  }

private:
  std::vector<ScalarType> d_ = {};
};

This custom vector class is very simple but the key point is that the CustomVector is unknown to Pressio. You can replace CustomVector with any type that fits your needs.

Problem Class

Here, we want to integrate in time the same system of ODEs shown in tutorial1, but with the problem implemented using the CustomVector class as follows:

struct MySystem
{
  using scalar_type   = double;
  using state_type    = CustomVector<scalar_type>;
  using velocity_type = state_type;

public:
  void velocity(const state_type & y,
                scalar_type t,
                velocity_type & rhs) const
  {
    constexpr scalar_type ten{10};
    rhs[0] = ten * y[0];
    rhs[1] = ten * y[1];
    rhs[2] = ten * y[2];
  };

  velocity_type createVelocity() const{
    return velocity_type(3);
  };
};

Specialize trait and ops

Because we are working with custom data types, we need to provide the necessary operations to do the algebra that pressio needs. This is done as follows:

namespace pressio{

template<class T> struct Traits<CustomVector<T>>{
  using scalar_type = typename CustomVector<T>::value_type;
};

namespace ops{
template<class T>
void deep_copy(CustomVector<T> & dest, const CustomVector<T> & from){
  dest = from;
}

template<class T>
CustomVector<T> clone(const CustomVector<T> & src){
  return CustomVector<T>(src.extent(0));
};

template<class T>
void update(CustomVector<T> & v,        const double a,
            const CustomVector<T> & v1, const double b)
{
  for (size_t i=0; i<v.extent(0); ++i){
    v[i] = a*v[i] + b*v1[i];
  }
}

template<class T>
void update(CustomVector<T> & v,        const double a,
            const CustomVector<T> & v1, const double b,
        const CustomVector<T> & v2, const double c)
{
  for (size_t i=0; i<v.extent(0); ++i){
    v[i] = a*v[i] + b*v1[i] + c*v2[i];
  }
}

template<class T>
void update(CustomVector<T> & v,        const double a,
            const CustomVector<T> & v1, const double b,
        const CustomVector<T> & v2, const double c,
        const CustomVector<T> & v3, const double d,
        const CustomVector<T> & v4, const double e)
{
  for (size_t i=0; i<v.extent(0); ++i){
    v[i] = a*v[i] + b*v1[i] + c*v2[i] + d*v3[i] + e*v4[i];
  }
}
}} //end namespace pressio::ops

Main

int main(int argc, char *argv[])
{
  pressio::log::initialize(pressio::logto::terminal);

  // create the system/app object
  using app_t = MySystem;
  app_t appObj;

  // create and initialize the state
  using state_t = typename app_t::state_type;
  state_t y(3);
  y[0] = 1.; y[1] = 2.; y[2] = 3.;

  // instantiate the stepper
  namespace pode = pressio::ode;
  constexpr auto scheme = pode::StepScheme::ForwardEuler;
  auto stepperObj = pode::create_explicit_stepper(scheme, y, appObj);

  // integrate in time
  const double dt = 0.1;
  const pressio::ode::step_count_type numSteps = 1;
  ::pressio::ode::advance_n_steps(stepperObj, y, 0.0, dt, numSteps);

  // check solution
  std::cout << "Computed solution: ["
            << y[0] << " " << y[1] << " " << y[2] << "] "
        << "Expected solution: [2,4,6] "
        << std::endl;

  pressio::log::finalize();
  return 0;

Full Code

The full code is available here.