pressio/pressio

core C++ library

Overview

Pressio is an open-source computational framework aimed at advancing the field of reduced-order models (ROMs) for dynamical systems in science and engineering. We employ generic programming, and target shared and distributed-memory applications using arbitrary data-types and diverse programming models.

We believe that the key to develop such a capability from the ground up is to properly identify the building blocks and modularize them accordingly. This is the approach adopted here: the library is composed of several components that can be used independently and in a self-contained fashion, but as a whole constitute the stack foundation of the rom component. One of the main outcome of this is that regardless of your interest in ROMs, you might find useful some of the components of the library.

Click below to checkout the documentation:

Development, versioning and backward compatibility

Until we reach a stable 1.0, please be patient and do not be surprised if the API and functionalities somewhat rapidly change from one 0.x version to the next, thus affecting backward compability. Some components of pressio are more mature and stable than others, but until we can claim the same stability level for all, please keep this mind.

Questions?

Find us on Slack: https://pressioteam.slack.com and/or open an issue on github.

License and Citation

The full license is available here.

At some point we plan to publish this, for now we have an arXiv preprint at: https://arxiv.org/abs/2003.07798.

pressio's Issues

ode: API for arbitrary (user-defined) implicit stepper

The API to let user define how to do a single step.
So baiscally user is in charge of computing the time-discrete operator at each step.

class ImplicitOdeArbitraryIf{
public:
  using scalar_type = 
  using state_type = 
  using residual_type = 
  using jacobian_type = 

  template <typename step_t, typename ... Args>
  void timeDiscreteResidual(const step_t & step,
                            const scalar_type & time,
                            residual_type & R,
                            Args & ... states) const
  {
    // forward to whatever approriate impl method, e.g.
    // timeDiscreteResidualImpl<step_t>( step, time, f, std::forward<Args>(states)... );
  }

  template <typename step_t, typename ... Args>
  void timeDiscreteJacobian(const step_t & step,
                            const scalar_type & time,
                            jacobian_type & J,
                            Args & ... states) const
  {
    // forward to whatever approriate impl method, e. g.
    // timeDiscreteJacobianImpl<step_t>( step, time, f, std::forward<Args>(states)... );
  }

  template <typename step_t, typename ... Args>
  residual_type timeDiscreteResidual(const step_t & step,
                                     const scalar_type & time,
                                     Args & ... states) const
  {
    // initialize (depends on the app and data structure types used)
    residual_type R( /*whatever_it_needs_to_be_constructed*/ );
    this->timeDiscreteResidual(step, time, R, std::forward<Args>(states)...);
    return R;
  }

  template <typename step_t, typename ... Args>
  jacobian_type timeDiscreteJacobian(const step_t & step,
                                     const scalar_type & time,
                                     Args & ... states) const
  {
    // initialize (depends on the app and data structure types used)
    jacobian_type J( /*whatever_it_needs_to_be_constructed*/ );
    this->timeDiscreteJacobian(step, time, J, std::forward<Args>(states)...);
    return J;
  }

rom: enable support for linear space-time ROMs

Enable support for linear space-time ROMs inside pressio/rom.
This should reuse the preliminary work that Chi did in FY19, see #9

• create/finalize a tentative design plan
• identify which steps are the main bottlenecks
• identify main challenges of integrating the design into pressio
• draft an implementation inside an experimental namespace inside pressio/rom
• create unit tests
• create regular tests
• move functionality into main namespace

rom/tests: main files should be refactored to avoid code duplication

rom/tests: main files should be refactored to benefit reusability

rom: interface for regular LSPG using *residual API*

The proposed residual api interface for LSPG.

class RomLSPGResidualAPI{
public:
  using scalar_type = 
  using state_type = 
  using residual_type = 
  using jacobian_type = 
  using dense_matrix_type = 

public:
  template <typename step_t, typename ... Args>
  void timeDiscreteResidual(const step_t & step,
                            const scalar_type & time,
                            const scalar_type & dt,
                            residual_type & R,
                            Args && ... states) const
  {
    // forward to whatever appropriate impl method, e.g.
    // timeDiscreteResidualImpl<step_t>( step, time, dt, f, std::forward<Args>(states)... );
    // of which we have multiple overloads based on the number of states 
  }

  template <typename step_t, typename ... Args>
  residual_type timeDiscreteResidual(const step_t & step,
                                     const scalar_type & time,
                                     const scalar_type & dt,
                                     Args && ... states) const
  {
    // initialize (depends on the app and data structure types used)
    residual_type R( /*whatever_it_needs_to_be_constructed*/ );
    this->timeDiscreteResidual(step, time, R, std::forward<Args>(states)...);
    return R;
  }

  // computes: time discrete Jacobian, Jac, and then A = Jac B
  void applyTimeDiscreteJacobian(const state_type & state,
                                 const dense_matrix_type & B,
                                 const scalar_type & t,
                                 dense_matrix_type & A) const{
    // compute Jacobian 
    // compute time-discrete Jacobian 
    // apply time-discrete Jacobian to B and store into A 
  }

  // computes: time discrete Jacobian, Jac, and then A = Jac B
  dense_matrix_type applyTimeDiscreteJacobian(const state_type & state,
                    			      const dense_matrix_type & B,
                                              const scalar_type & t) const{
    dense_matrix_type A( /*whatever_it_needs_to_be_constructed*/ );
    this->applyTimeDiscreteJacobian(state, B, t, A);
    return A;
  }
};

ode: set to zero the StatesContainers when constructed

ode: set to zero the StatesContainers when constructed, except for the first entry which is supposed to be the initial condition which makes sense to be nonzero.

containers: increase size of static collection to something reasonable (10)

pybind11: fix specializations for when user-defined ops are void

For pybind11-related code, we need to fix specializations for when user-defined ops are void.
So that if ops are void, we use blas. Otherwise we use those ops.

ode: interface for implicit stepper (regular stepper, not user-defined)

Current API admissible to ode implicit stepper.
This is only valid when using ImplicitEnum::Euler or BDF2. NOT valid for ImplicitEnum::Arbitrary.

class ImplicitOdeIf{
public:
  using scalar_type = 
  using state_type = 
  using velocity_type = 
  using jacobian_type = 

  // compute velocity, f(x,t;...), for a given state, x
  void velocity(const state_type  & state,
                const scalar_type & time,
                velocity_type     & f) const{
    //...
  }

  // first compute the Jacobian, J, 
  // then compute A=J*B where B=jacobian of ROM decoder
  void jacobian(const state_type  & state,
                jacobian_type     & J)const{
    //...
  }

  // overload called once to construct initial object
  velocity_type velocity(const state_type  & state,
                         const scalar_type & t) const
  {
    velocity_type f( /*whatever_it_needs_to_be_constructed*/ );
    this->velocity(state, t, f);
    return f;
  }

  // overload called once to construct initial object
  jacobian_type jacobian(const state_t   & state,
                         const scalar_t  & t) const
  {
    jacobian_type J( /*whatever_it_needs_to_be_constructed*/ );
    this->jacobian(state, t, J);
    return J;
  }

create a unique class for containing a set of containers

need a single class to contain a set of containers, which can be vector or something else.
This can then be inherited from to customize different functionalities, but the basic container of containers should be in a single place.
We might need two versions, one for when we know how many container objects we need at compile time and one when we don't know.

This allows us to refactor ode::odeStorage, ode::odeVelocitiesContainer and `rom::fomStatesContainers'

ode: which implicit time stepping schemes should pressio support?

@eparish1 @pjb236 @vbrunini

Currently, we have BDF1, BDF2.
Any suggestions on which other implicit ones we should have? (this is beside the functionality to let user define an arbitrary scheme), what do you guys think we should natively support?

solvers: geqrf enable support for repeated calls

solvers: geqrf enable support for repeated calls.
Right now, it is not very robust. Need to enable workspace to be adaptable if the matrix size
changes over multiple sequential calls.

rom: interface for regular LSPG using velocity API

The current interface based on the velocity api for regular LSPG is:

class LSPGVelocityAPI{
public:
  using scalar_type = 
  using state_type = 
  using velocity_type = 
  using jacobian_type = 
  using dense_matrix_type = 

public:
  void velocity(const state_type  & state, 
                const scalar_type & time, 
                velocity_type & f) const{
    //...
  }

  velocity_type velocity(const state_type & u,
			 const scalar_type & t) const{
    velocity_type f(/*whatever_it_needs_to_be_constructed*/);
    this->velocity(u, t, f);
    return f;
  }

  // computes: A = Jac B
  void applyJacobian(const state_type & state,
                     const dense_matrix_type & B,
                     const scalar_type & t,
                     dense_matrix_type & A) const{
    //...
  }

  // computes: A = Jac B
  dense_matrix_type applyJacobian(const state_type & state,
                                  const dense_matrix_type & B,
                                  const scalar_type & t) const{
    dense_matrix_type A( /*whatever_it_needs_to_be_constructed*/ );
    this->applyJacobian(state, B, t, A);
    return A;
  }

  // -------------------------------
  // for preconditioned LSPG, you need the methods below
  // -------------------------------
  // TODO: fill in
};

rom: fom states reconstruction

Currently, the FOM states when using LSPG (both with velocity and residual API) are reconstructed only by the residual policy and the Jacobian policy simply uses them (since they share the data). This implies we need to ensure the residual policy to be called BEFORE the Jacobian.
Right now this is the case but there is no check for it.
Problem: find a way to make sure residual and jacobian policiies are called in that order and if the residual policy is called before the jacobian one then either throw and error, or do the reconstruction in the Jacobian policy and then let residual one use that.

rom: need to detect user-defined ops

make eigen an optional dependency

so far we kept eigen as required dependency, but it should really be optional.

improving wiki

• general description
• page on how to do a serial build
• page on how to do a MPI build

apps: remove setup method

apps: remote setup method, it is confusing

rom: re-enable pybind11 code

re-enable pybind11 code. Commented out for now.

solvers: dumping generalized coordinates when a NaN is detected in the

This would provide a useful debugging tool, since I could reconstruct the full state that caused the NaNs with the generalized coordinates and POD basis. Either we print the generalized coordinates prior to screen or to disk. I would personally prefer printing to disk, but either way works.

Improve tests

The goal of this issue is to discuss Pressio testing and CI

merge PyGaussNewton with GaussNewton

need to merge PyGaussNewton with GaussNewton. Might need some refactoring to handle pybind11 objects but it is better than keeping two distinct solvers that do the same thing

ode: Support integrateToTargetTime

Similar to integrateNSteps but integrates until an end time is reached regardless of the number of steps it takes. This is particularly useful when using a user-provided time step or other method of adaptive time stepping.

identify best TPLs to leverage for autoencoders

identify best TPLs to leverage for autoencoders to use for manifold ROMs

rom: add static check that fom_type supports mask

solvers: fix template parameter for system

in solvers, I need to make sure the system type passed to the solve method matches the one passed to the class

solvers: exit when line search compute zero update?

Currently the line search terminates when the solution increment drops below 1e-14. If this is the case, then shouldn't we also stop the nonlinear iteration? Having a solution increment that's machine zero in magnitude means we've found a local extrema and there's no escape*. Currently the solver just keeps iterating until it reaches the maximum iteration count, which is a waste of time in my opinion.

I think we should have the line search output a boolean to tell the nonlinear solver if it found a non-zero solution increment. If it did not, then the nonlinear solver should terminate, perhaps with some kind of warning (or error?) about poor solver convergence.

*Unless we use levenberg-marquardt or something else that can escape local extrema.

rom: add static check that fom_type supports preconditioner

all: revise special member functions in classes

ode: add integrator function without dt to use for arbitrary stepper

Add integrator function without dt to use for arbitrary stepper.
This allows user to set the time step at each iteration without pressio knowing about it.

all: ``package-specific'' helper headers under ./packages/ should be checked

global includes in the packages dir should be checked. Maybe things are missing or some files do not need to be included here.

fix the min version of cmake needed

reenable some tests after moving to pure cmake

copy/mv constructors for odeStorage

Do I need the cp/mv constr and assignment operator for `odeStogage'?
With pybind11, I cannot delete them.

  OdeStorage() = delete;
  OdeStorage(const OdeStorage &) = ?;
  OdeStorage & operator=(const OdeStorage &) = ?;
  OdeStorage(OdeStorage &&) = ?;
  OdeStorage & operator=(OdeStorage &&) = ?;

solvers: reuse residual computed in line-search to avoid recomputing it

@eparish1 rightfully suggested that in the GN solvers, when the line search is on, we should reuse the residual computed in the line-search to avoid recomputing it.

all: design exception handling where missing

self-explanatory :)

move issues from old hosting here

query FOM for the time-discrete residual

This feature seems to be useful for a few cases.

• List feasibility and limitations (e.g does this seamlessly apply to ST, TC-LSRM and other ROMs?)
• create/finalize a tentative design plan
• identify which steps are the main bottlenecks
• identify main challenges of integrating the design into pressio
• draft an implementation inside an experimental namespace inside pressio/rom
• create unit tests
• create regular tests
• move functionality into main namespace

Improve documentation

This issue is mostly to discuss how documentation shall and will be improved. In particular, we should consider using the embedded Wiki.

add developers guidelines

gtest: fix mains for custom environment

gtest: fix mains for custom environment. for some reason current unique_ptr segfaults.

But the following works:

  ::testing::InitGoogleTest(&argc,argv);
  int err = 0;
  {
    // std::unique_ptr<MPIEnv> envPtr(new MPIEnv(argc, argv));
    // ::testing::AddGlobalTestEnvironment(envPtr.get());

    int rank = {};
    int mpiError = MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    //ASSERT_FALSE(mpiError);

    err = RUN_ALL_TESTS();

    mpiError = MPI_Finalize();
    //ASSERT_FALSE(mpiError);
  }

improve namespaces in rom package

improve namespaces in rom package.
Right now we have classes like

pressio::rom::LSPGUnsteadyProblem

we should make this:

pressio::rom::lspg::UnsteadyProblem

Similarly for Galerkin

gn solver for pybind: fix initialization of hessian

Fix this so that we remove the multiply call and only init the hessian using proper sizes that we can infer from jac. However, we need to find out if pybind11::array has a constructor accepting sizes

// need to figure out a way to init just using sizes
// rather than calling this multiply
hess_(pythonOps_.attr("multiply")(jac_, true, jac_, false)),

New Gauss-Newton Solver

#6
We need a new Gauss-Newton solver for an efficient WLS implementation. The solver should be the same as the current GN solver with the following differences:

1.) Instead of having routines to compute J, r, J^TJ, and J^TR, we have one single routine that computes J^TJ and J^TR; i.e., the ``system" should have a member

JTJ_and_JTR(state_t & state, JTJ_type & JTJ, JTR_type & JTR) 

The reason for this is so that we can compute these products efficiently for different systems.

2.) It is important to avoid computing anything that stores an entire ""windowed" FOM vector (i.e. the residual over the window).

The sample code would look like

namespace pressio{ namespace solvers{ namespace iterative{ namespace impl{
template <
  typename line_search_t,
  typename converged_when_tag,
  typename system_t,
  typename lin_solver_t,
  typename iteration_t,
  typename scalar_t,
  typename observer_t = utils::impl::empty
  >
void gauss_newton_neq_solve_general(const system_t & sys,
          typename system_t::state_type & y,
          typename system_t::state_type & ytrial,
          typename system_t::JTR_type & JTR, // this should be the same as state_type
          typename system_t::JTJ_type & JTJ,
          typename system_t::state_type & dy,
          lin_solver_t & linSolver,
          iteration_t maxNonLIt,
          scalar_t tolerance,
          scalar_t & normO,
          scalar_t & norm_dy,
          const observer_t * observer,
          std::string & convCondDescr)
{

  using residual_t  = typename system_t::residual_type;
  using jacobian_t  = typename system_t::jacobian_type;

  // find out which norm to use
  using norm_t = typename NormSelectorHelper<converged_when_tag>::norm_t;

  // policy for evaluating the norm of a vector
  using norm_evaluator_t = ComputeNormHelper<norm_t>;

  // policy to checking convergence
  using is_converged_t = IsConvergedHelper<converged_when_tag>;

  /* policy for computing line search factor (alpha) such that
   * the update is done with y = y + alpha dy
   * alpha = 1 default when user does not want line search
   */
  using lsearch_helper_t = LineSearchHelper<line_search_t>;

  //-------------------------------------------------------

  // alpha for taking steps
  scalar_t alpha = {};
  // storing residual norm
  scalar_t normRes = {};
  scalar_t normRes0 = {};
  // storing projected residual norm
  scalar_t normJTRes = {};
  scalar_t normJTRes0 = {};

  constexpr auto one = ::pressio::utils::constants::one<scalar_t>();
  convCondDescr = std::string(is_converged_t::description_);

  // compute the initial norm of y (the state)
  norm_evaluator_t::evaluate(y, normO);
  norm_dy = {0};

  iteration_t iStep = 0;
  while (++iStep <= maxNonLIt)
  {


    // compute everything in one go
    sys.JTJ_and_JTR(y, JTJ,JTR);

    JTR.scale(static_cast<scalar_t>(-1));
    // projected residual norm for current state
    norm_evaluator_t::evaluate(JTR, normJTRes);

    // store initial residual norm
    if (iStep==1) normJTRes0 = normJTRes;

    // solve normal equations
    linSolver.solveAllowMatOverwrite(JTJ, JTR, dy);

    // compute norm of the correction
    norm_evaluator_t::evaluate(dy, norm_dy);

    // exit with error if NaNs detected in solution update dy
    if (std::isnan(norm_dy))
    {
      throw std::runtime_error(
        "Nonlinear solver: NEQ-based Gausss Newton: NaNs detected in solution update dy");
    }

    // compute multiplicative factor if needed
    // Eric's comment: not sure about how this is done in the code currently, 
    // but we only want JTJ and JTR as inputs. The line search would then require ||JTR|| to 
    // decrease
    lsearch_helper_t::evaluate(alpha, y, ytrial, dy,JTR,JTJ, sys);

    // solution update: y = y + alpha*dy
    ::pressio::containers::ops::do_update(y, one, dy, alpha);

    // check convergence (whatever method user decided)
   const auto flag = is_converged_t::evaluate(y, dy,
                 norm_dy,normJTRes, normJTRes0,
                 iStep, maxNonLIt, tolerance);

    // store new norm into old variable
    normO = norm_dy;


  }//loop

}// end

}}}} //end namespace pressio::solvers::iterative::impl
#endif

Optimal Sample Mesh Design

Currently our sample meshes are generated randomly (with some constraints that at least 1 element lies on each sideset, etc). There is agreement that this is likely sub optimal. We should attempt to formulate more effective ways of generating the sample mesh based on the snapshot data or data from the snapshot run.

Current idea involves using the adapter class to query the product of the Jacobian and the basis matrix at different values of the state vector chosen based on the snapshots.

rom: rename the fom reconstructor

Now we have:

lspgProblem.yFomReconstructor_(xROM)

chane so that there is no ''y''.

fomStateReconstructor_(xROM)

rom: add support for WLS

Implement WLS inside pressio/rom

We have added support for WLS in Pressio. We currently support Eigen datatypes, explicit Euler, BDF1, and BDF2. We do not support other datatypes or preconditioning yet.

We are working on adding:

• Support for Kokkos
• Increased number of unit tests
• Support for hyper reduction
• Support for preconditioning
• Support for parallel in time over each window
• Support for QR solution technique
• (Maybe) Support for block Jacobi

rom: move preliminary work for space-time from old repo in gitlab

rom: move preliminary work for space-time that Chi did in FY19 from old repo in gitlab to here.
See #8

containers: add functionality to multivector for [ ] method

Add functionality such that an object of type pressio::containers::MultiVector<....>
has the method "operator[](int i)", which returns the ith column of the multivector.

ode: interface for explicit steppers

The interface for explicit steppers is:

class ExplicitOdeIf{
public:
  using scalar_type = 
  using state_type = 
  using velocity_type = 

  // compute velocity, f(x,t;...), for a given state, x
  void velocity(const state_type  & state, 
                const scalar_type & time,
                velocity_type     & f) const{
    //...
  }

  // overload called once to construct initial object
  velocity_type velocity(const state_type &  state,
                         const scalar_type & time) const
  {
    velocity_type f( /*whatever_it_needs_to_be_constructed*/ );
    this->velocity(state, time, f);
    return f;
  }

rom: redundant reconstruction of FOM prev states

rom: it seems to me we have a redundant reconstruction of FOM prev states in the residual policies. Is there a way we can make the stepper communicate to the policy when do that so that we don't repeat same things many times? The answer is yes, but I need to figure out the best way.

One way to do this would be to consider that we need to reconstruct the FOM states at previous steps only when the stepper is starting a new step. So we have the stepper tell the policy when a new step is being beginning. This would trigger the need to reconstuct previos states.