Default Galerkin with Explicit Time Integration for Custom Types

Custom Data Class

Suppose that you have an application that uses arbitrary types. For the sake of this tutorial, let's define these to be:

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_ = {};
};

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

  CustomMatrix(std::size_t nr, std::size_t nc)
    : num_rows_(nr), num_cols_(nc), d_(nr*nc){}

  std::size_t extent(int k)const { return (k==0) ? num_rows_ : num_cols_; }

  ScalarType & operator()(int i, int j){ return d_[num_cols_*i+j]; }
  const ScalarType & operator()(int i, int j) const { return d_[num_cols_*i+j]; }

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

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

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 via specialization as follows:

namespace pressio{

template<class ScalarType>
struct Traits<CustomVector<ScalarType>>{
  using scalar_type = ScalarType;
  using size_type = std::size_t;
};

template<class ScalarType>
struct Traits<CustomMatrix<ScalarType>>{
  using scalar_type = ScalarType;
  using size_type = std::size_t;
};

namespace ops{

template<class ScalarType>
std::size_t extent(CustomVector<ScalarType> & object, int i){
    return object.extent(i);
}

template<class ScalarType>
std::size_t extent(CustomMatrix<ScalarType> & object, int i){
    return object.extent(i);
}

template<class ScalarType>
void set_zero(CustomVector<ScalarType> & object){
  object.fill(0);
}

template<class ScalarType>
void set_zero(CustomMatrix<ScalarType> & object){
  object.fill(0);
}

template<class ScalarType>
void deep_copy(CustomVector<ScalarType> & dest,
               const CustomVector<ScalarType> & src){
  dest = src;
}

template<class ScalarType>
void deep_copy(CustomMatrix<ScalarType> & dest,
               const CustomMatrix<ScalarType> & src){
  dest = src;
}

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

template<class ScalarType>
CustomMatrix<ScalarType> clone(const CustomMatrix<ScalarType> & src){
  return CustomMatrix<ScalarType>(src.extent(0), src.extent(1));
}

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

// z = beta*z + alpha * A * x
// where x is indexable as x(i)
template<class x_t, class ScalarType>
void product(pressio::nontranspose,
             ScalarType alpha,
             const CustomMatrix<ScalarType> & A,
             const x_t & x,
             ScalarType beta,
             CustomVector<ScalarType> & z)
{
  // obviously not efficient, just for demonstration
  for (std::size_t i=0; i<A.extent(0); ++i)
  {
    z(i) = beta*z(i);
    for (std::size_t j=0; j<A.extent(1); ++j){
      z(i) += alpha*A(i,j)*x(j);
    }
  }
}

// z = beta*z + alpha * A^T * x
// where z is indexable as z(i)
template<class z_t, class ScalarType>
void product(pressio::transpose,
             ScalarType alpha,
             const CustomMatrix<ScalarType> & A,
             const CustomVector<ScalarType> & x,
             ScalarType beta,
             z_t & z)
{
  // obviously not efficient, just for demonstration
  for (std::size_t k=0; k<A.extent(1); ++k)
  {
    z(k) = beta*z(k);
    for (std::size_t i=0; i<A.extent(0); ++i){
      z(k) += alpha*A(i,k)*x(i);
    }
  }
}

// C = beta*C + alpha * A^T * B
// where C is indexable as C(i,j)
template<class C_t, class ScalarType>
void product(pressio::transpose,
             pressio::nontranspose,
             ScalarType alpha,
             const CustomMatrix<ScalarType> & A,
             const CustomMatrix<ScalarType> & B,
             ScalarType beta,
             C_t & C)
{
  for (std::size_t i=0; i<A.extent(1); ++i){
    for (std::size_t j=0; j<B.extent(1); ++j)
    {
      C(i,j) *= beta;
      for (std::size_t k=0; k<A.extent(0); ++k){
        C(i,j) += alpha*A(k,i)*B(k,j);
      }
    }
  }
}

}}//end namespace pressio::ops

Main

#include "pressio/rom_galerkin.hpp"

struct TrivialFomOnlyVelocityCustomTypes
{
  using scalar_type    = double;
  using state_type     = CustomVector<scalar_type>;
  using velocity_type  = state_type;
  int N_ = {};

  TrivialFomOnlyVelocityCustomTypes(int N): N_(N){}

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

  void velocity(const state_type & u, const scalar_type time, velocity_type & f) const{
    for (std::size_t i=0; i<f.extent(0); ++i){
     f(i) = u(i) + time;
    }
  }
};

struct Observer
{
  void operator()(int32_t step, double time, Eigen::VectorXd state)
  {
    std::cout << "Observer called: \n"
    << " step: " << step
    << "\n"
    << " state = ";
    for (int i=0; i<state.size(); ++i){
        std::cout << state(i) << ", ";
    }
    std::cout << "\n";
  }
};

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

  using fom_t = TrivialFomOnlyVelocityCustomTypes;
  using scalar_t    = typename fom_t::scalar_type;
  using fom_state_t = typename fom_t::state_type;

  constexpr int N = 10;
  fom_t fomSystem(N);
  fom_state_t fomReferenceState(N);
  fomReferenceState.fill(0);

  using phi_t = ::CustomMatrix<scalar_t>;
  phi_t phi(N, 3);
  for (std::size_t i=0; i<phi.extent(0); ++i){
    phi(i,0) = 0.1;
    phi(i,1) = 0.5;
    phi(i,2) = 0.8;
  }
  auto decoder = pressio::rom::create_time_invariant_linear_decoder<fom_state_t>(phi);

  Eigen::VectorXd romState(3);
  romState[0]=0.;
  romState[1]=1.;
  romState[2]=2.;

  constexpr auto odescheme = pressio::ode::StepScheme::ForwardEuler;
  auto problem = pressio::rom::galerkin::create_default_explicit_problem
    (odescheme, fomSystem, decoder, romState, fomReferenceState);

  const scalar_t dt = 1.;
  const int num_steps = 3;
  Observer obs;
  pressio::ode::advance_n_steps_and_observe(problem, romState, 0., dt, num_steps, obs);

  pressio::log::finalize();
}