Integrator.h

Go to the documentation of this file.

#ifndef KINETX_INTEGRATOR_H_
#define KINETX_INTEGRATOR_H_

#include "Kinetx/Network.h"

namespace Scine {
namespace Kinetx {
class Integrator {
 public:
  Integrator() = default;

  virtual void propagate(Eigen::VectorXd& concentrations, Eigen::VectorXd& yFlux, Eigen::VectorXd& rFlux,
                         Eigen::VectorXd& rForwardFlux, Eigen::VectorXd& rBackwardFlux, double& t, double& dt) const = 0;
  virtual Eigen::MatrixXd runIntegration(Eigen::VectorXd y, double tStart, double dt, Eigen::VectorXd& rFlux,
                                         Eigen::VectorXd& rForwardFlux, Eigen::VectorXd& rBackwardFlux,
                                         const unsigned int batchInterval = 1000, const unsigned int nBatches = 100000,
                                         const double convergenceConcentrationChange = 1e-10) = 0;
  virtual Eigen::MatrixXd runIntegrationByTime(Eigen::VectorXd y, double tStart, double dt, Eigen::VectorXd& rFlux,
                                               Eigen::VectorXd& rForwardFlux, Eigen::VectorXd& rBackwardFlux,
                                               const double tMax, const unsigned int batchInterval = 1000,
                                               const double convergenceConcentrationChange = 1e-10) = 0;
};

class IntegratorBase : public Integrator {
 public:
  IntegratorBase(Network& net);

  virtual void propagate(Eigen::VectorXd& concentrations, Eigen::VectorXd& yFlux, Eigen::VectorXd& rFlux,
                         Eigen::VectorXd& rForwardFlux, Eigen::VectorXd& rBackwardFlux, double& t, double& dt) const = 0;

  virtual Eigen::MatrixXd runIntegration(Eigen::VectorXd y, double tStart, double dt, Eigen::VectorXd& rFlux,
                                         Eigen::VectorXd& rForwardFlux, Eigen::VectorXd& rBackwardFlux,
                                         const unsigned int batchInterval = 1000, const unsigned int nBatches = 100000,
                                         const double convergenceConcentrationChange = 1e-10);

  virtual Eigen::MatrixXd runIntegrationByTime(Eigen::VectorXd y, double tStart, double dt, Eigen::VectorXd& rFlux,
                                               Eigen::VectorXd& rForwardFlux, Eigen::VectorXd& rBackwardFlux,
                                               const double tMax, const unsigned int batchInterval = 1000,
                                               const double convergenceConcentrationChange = 1e-10);

 protected:
  void trackVertexAndEdgeFluxes(const Eigen::VectorXd& y, const Eigen::VectorXd& yInitial, Eigen::VectorXd& yFlux,
                                Eigen::VectorXd& rFlux, Eigen::VectorXd& rForwardFlux, Eigen::VectorXd& rBackwardFlux,
                                const double& dt) const;
  Eigen::VectorXd g(const Eigen::VectorXd& concentrations) const;
  Eigen::VectorXd gFlux(const Eigen::VectorXd& concentrationsT0, const Eigen::VectorXd& concentrationsT1,
                        Eigen::VectorXd& flux, Eigen::VectorXd& forwardFlux, Eigen::VectorXd& backwardFlux) const;
  Eigen::MatrixXd fFlux(const Eigen::VectorXd& concentrationsT0, const Eigen::VectorXd& concentrationsT1) const;
  std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> fFluxDirected(const Eigen::VectorXd& concentrationsT0,
                                                                              const Eigen::VectorXd& concentrationsT1) const;
  Eigen::MatrixXd f(const Eigen::VectorXd& concentrations) const;
  std::pair<Eigen::MatrixXd, Eigen::MatrixXd> fDirected(const Eigen::VectorXd& concentrations) const;
  Eigen::SparseMatrix<double> jacobi(const Eigen::VectorXd& concentrations) const;
  Network& _net;

  void printHeader(const Eigen::VectorXd& y, const double dt, const double t);
  bool printTimeAndCheckConvergenceStep(const Eigen::VectorXd& y, const Eigen::VectorXd& yOld, Eigen::VectorXd& yMax,
                                        const unsigned int batchInterval, const unsigned int iBatch, const double dt,
                                        const double t, const double convergenceConcentrationChange);

  static double fast_pow(const double& x, const int& n) {
    if (n == 0) {
      return 1.0;
    }
    if (n == 1) {
      return x;
    }
    if (n == 2) {
      return x * x;
    }
    if (n == 3) {
      return x * x * x;
    }
    if (n == 4) {
      const double tmp = x * x;
      return tmp * tmp;
    }
    return std::pow(x, n);
  }
};

} /* namespace Kinetx */
} /* namespace Scine */

#endif // KINETX_INTEGRATOR_H_