Default Galerkin with Implicit Time Integration for Known Types

#include "pressio/type_traits.hpp"
#include "pressio/ops.hpp"
#include "pressio/rom_galerkin.hpp"

struct TrivialFomVelocityAndJacobianEigen
{
  using scalar_type       = double;
  using state_type        = Eigen::VectorXd;
  using velocity_type     = state_type;
  int N_ = {};

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

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

  template<class OperandType>
  OperandType createApplyJacobianResult(const OperandType & B) const
  {
    OperandType A(B.rows(), B.cols());
    return A;
  }

  // computes: A = Jac B
  template<class OperandType>
  void applyJacobian(const state_type & state,
                     const OperandType & B,
                     const scalar_type & time,
                     OperandType & A) const
  {
    A = B;
    A.array() += time;
  }

  void velocity(const state_type & u, const scalar_type time, velocity_type & f) const
  {
    for (auto i=0; i<f.rows(); ++i){
     f(i) = u(i) + time;
    }
  }
};

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

  namespace pls    = pressio::linearsolvers;
  namespace pnonls = pressio::nonlinearsolvers;
  namespace prom   = pressio::rom;

  using fom_t = TrivialFomVelocityAndJacobianEigen;
  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.setZero();

  using phi_t = Eigen::MatrixXd;
  phi_t phi(N, 3);
  phi.col(0).setConstant(0.);
  phi.col(1).setConstant(1.);
  phi.col(2).setConstant(2.);
  auto decoder = prom::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::BDF1;
  auto problem = pressio::rom::galerkin::create_default_implicit_problem
    (odescheme, fomSystem, decoder, romState, fomReferenceState);
  using problem_t = decltype(problem);

  using galerkin_jacobian_t = typename pressio::Traits<problem_t>::galerkin_jacobian_type;
  using lin_solver_t = pls::Solver<pls::iterative::LSCG, galerkin_jacobian_t>;
  lin_solver_t linearSolverObj;
  auto nonLinSolver = pnonls::create_newton_raphson(problem, romState, linearSolverObj);
  nonLinSolver.setMaxIterations(1);

  scalar_t dt = 2.;
  pressio::ode::advance_n_steps(problem, romState, 0.0, dt, 1, nonLinSolver);
  std::cout << romState << std::endl;

  pressio::log::finalize();
}