BaseQNAcceleration.hpp

Go to the documentation of this file.

#pragma once
 
#include <Eigen/Core>
#include <algorithm>
#include <deque>
#include <fstream>
#include <map>
#include <sstream>
#include <string>
#include <vector>
#include "acceleration/Acceleration.hpp"
#include "acceleration/impl/QRFactorization.hpp"
#include "acceleration/impl/SharedPointer.hpp"
#include "logging/Logger.hpp"
#include "time/TimeGrids.hpp"
 
/* ****************************************************************************
 *
 * A few comments concerning the present design choice.
 *
 * All the functions from the base class BAseQNAcceleration are specialized in
 * the sub classes as needed. This is done vi overwriting the base functions in
 * the specialized sub classes and calling the respective base function after
 * performing the specialized stuff (in order to perform the common, generalized
 * computations i.e. handling of V,W matrices etc.)
 * However, for the performAcceleration Method we decided (for better readability)
 * to have this method only in the base class, while introducing a function
 * performPPSecondaryData that handles all the specialized stuff concerning acceleration
 * processing for the secondary data in the sub classes.
 *
 * Another possibility would have been to introduce a bunch of functions like
 * initializeSpecialized(), removeMatrixColumnSpecialized(),
 * iterationsConvergedSpecialized(), etc
 * and call those function from the base class top down to the sub classes.
 *
 * The third possibility was to separate the approximation of the Jacobian from
 * the common stuff like handling V,W matrices in the acceleration.
 * Here, we have a class QNAcceleration that handles the V,W stuff and the basic
 * scheme of the QN update. Furthermore we have a base class (or rather interface)
 * JacobianApproximation with sub classes IQNIMVJAPX and IQNAPX that handle all the
 * specialized stuff like Jacobian approximation, handling of secondary data etc.
 * However, this approach is not feasible, as we have to call the function
 * removeMatrixColumn() down in the specialized sub classes IQNIMVJApx and IQNApx.
 * This is not possible as the function works on the V, W matrices that are
 * completely treated by QNAcceleration.
 *
 * ****************************************************************************
 */
 
// ----------------------------------------------------------- CLASS DEFINITION
 
namespace precice {
namespace io {
class TXTReader;
class TXTWriter;
} // namespace io
 
namespace acceleration {

class BaseQNAcceleration : public Acceleration {
public:
  BaseQNAcceleration(
      double                  initialRelaxation,
      bool                    forceInitialRelaxation,
      int                     maxIterationsUsed,
      int                     timeWindowsReused,
      int                     filter,
      double                  singularityLimit,
      std::vector<int>        dataIDs,
      OnBoundViolation        onBoundViolation,
      impl::PtrPreconditioner preconditioner,
      bool                    reducedTimeGrid);

  ~BaseQNAcceleration() override
  {
    // not necessary for user, only for developer, if needed, this should be configurable
    //     if (utils::IntraComm::isPrimary() || !utils::IntraComm::isParallel()) {
    //       _infostream.open("precice-accelerationInfo.log", std::ios_base::out);
    //       _infostream << std::setprecision(16);
    //       _infostream << _infostringstream.str();
    //     }
  }

  std::vector<int> getPrimaryDataIDs() const final override
  {
    return _primaryDataIDs;
  }

  void initialize(const DataMap &cplData) final override;

  void performAcceleration(DataMap &cplData, double windowStart, double windowEnd) final override;
  std::vector<DataID> checkBoundViolation(Eigen::VectorXd &data, DataMap &cplData) const;
  void onBoundViolationsClamp(Eigen::VectorXd &data, DataMap &cplData, const std::vector<DataID> &violatingIDs);
  void onBoundViolationsScale(Eigen::VectorXd &data, DataMap &cplData, const std::vector<DataID> &violatingIDs, Eigen::VectorXd &xUpdate);
  void clampToBounds(Eigen::Map<Eigen::MatrixXd> &data, const std::vector<std::optional<double>> &lowerBound, const std::vector<std::optional<double>> &upperBound);
  double computeShorteningFactor(Eigen::Map<Eigen::MatrixXd> &data, Eigen::Map<Eigen::MatrixXd> &update, const std::vector<std::optional<double>> &lowerBound, const std::vector<std::optional<double>> &upperBound);
  void iterationsConverged(const DataMap &cplData, double windowStart) final override;

  void exportState(io::TXTWriter &writer) final override;

  void importState(io::TXTReader &reader) final override;

  int getDeletedColumns() const;

  int getDroppedColumns() const;

  int getLSSystemCols() const;

  int getMaxUsedIterations() const;

  int getMaxUsedTimeWindows() const;

  void addLogEntries(io::TXTTableWriter &writer) const override;

  void writeLogEntries(io::TXTTableWriter &writer) const override;
 
protected:
  logging::Logger _log{"acceleration::BaseQNAcceleration"};

  impl::PtrPreconditioner _preconditioner;

  const double _initialRelaxation;

  const int _maxIterationsUsed;

  const int _timeWindowsReused;

  const std::vector<int> _primaryDataIDs;

  std::vector<int> _dataIDs;

  std::vector<DataID> _idsWithBounds;

  bool _firstIteration = true;
 
  /* @brief Indicates the first time window, where constant relaxation is used
   *        later, we replace the constant relaxation by a qN-update from last time window.
   */
  bool _firstTimeWindow = true;
 
  /*
   * @brief True if this process has nodes at the coupling interface
   */
  bool _hasNodesOnInterface = true;
 
  OnBoundViolation _onBoundViolation;
 
  /* @brief If true, the QN-scheme always performs a underrelaxation in the first iteration of
   *        a new time window. Otherwise, the LS system from the previous time window is used in the
   *        first iteration.
   */
  bool _forceInitialRelaxation;

  bool _resetLS = false;

  Eigen::VectorXd _oldPrimaryXTilde;

  Eigen::VectorXd _oldXTilde;

  Eigen::VectorXd _primaryResiduals;

  Eigen::VectorXd _residuals; // @todo is this member still needed? Potential refactoring.

  Eigen::MatrixXd _matrixV;

  Eigen::MatrixXd _matrixW;

  const bool _reducedTimeGrid;

  impl::QRFactorization _qrV;

  const int _filter;

  const double _singularityLimit;

  std::deque<int> _matrixCols;

  std::vector<int> _dimOffsets;

  std::vector<int> _dimOffsetsPrimary;

  std::ostringstream _infostringstream;
  std::fstream       _infostream;
 
  int getLSSystemRows() const;
  int getPrimaryLSSystemRows() const;

  virtual void specializedIterationsConverged(const DataMap &cplData) = 0;

  virtual void updateDifferenceMatrices(const DataMap &cplData);

  virtual void updateCouplingData(const DataMap &cplData, double windowStart);

  virtual void applyFilter();

  virtual void computeQNUpdate(Eigen::VectorXd &xUpdate) = 0;

  virtual void removeMatrixColumn(int columnIndex);

  void writeInfo(const std::string &s, bool allProcs = false);
 
  int its = 0, tWindows = 0;
 
private:
  void initializeVectorsAndPreconditioner(const DataMap &cplData, double windowStart);

  virtual void specializedInitializeVectorsAndPreconditioner(const DataMap &cplData) = 0;

  void concatenateCouplingData(Eigen::VectorXd &data, Eigen::VectorXd &oldData, const DataMap &cplData, std::vector<int> dataIDs, precice::time::TimeGrids timeGrids, double windowStart) const;

  std::optional<time::TimeGrids> _timeGrids;

  std::optional<time::TimeGrids> _primaryTimeGrids;

  Eigen::VectorXd _primaryValues;

  Eigen::VectorXd _oldPrimaryValues;

  Eigen::VectorXd _oldPrimaryResiduals;

  Eigen::VectorXd _values;

  Eigen::VectorXd _oldValues;

  Eigen::MatrixXd _matrixVBackup;
  Eigen::MatrixXd _matrixWBackup;
  std::deque<int> _matrixColsBackup;

  int _nbDelCols = 0;

  int _nbDropCols = 0;
};
} // namespace acceleration
} // namespace precice