MNASolverDirect.h

/* Copyright 2017-2021 Institute for Automation of Complex Power Systems,
 *                     EONERC, RWTH Aachen University
 *
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at https://mozilla.org/MPL/2.0/.
 *********************************************************************************/
 
#pragma once
 
#include <bitset>
#include <iostream>
#include <list>
#include <memory>
#include <unordered_map>
#include <vector>
 
#include <dpsim/Config.h>
#include <dpsim/DataLogger.h>
#include <dpsim/DenseLUAdapter.h>
#include <dpsim/DirectLinearSolver.h>
#include <dpsim/DirectLinearSolverConfiguration.h>
#include <dpsim/Solver.h>
#ifdef WITH_SPARSE
#include <dpsim/SparseLUAdapter.h>
#endif
#ifdef WITH_KLU
#include <dpsim/KLUAdapter.h>
#endif
#ifdef WITH_CUDA
#include <dpsim/GpuDenseAdapter.h>
#ifdef WITH_CUDA_SPARSE
#include <dpsim/GpuSparseAdapter.h>
#endif
#ifdef WITH_MAGMA
#include <dpsim/GpuMagmaAdapter.h>
#endif
#endif
#include <dpsim-models/AttributeList.h>
#include <dpsim-models/SimPowerComp.h>
#include <dpsim-models/SimSignalComp.h>
#include <dpsim-models/Solver/MNASwitchInterface.h>
#include <dpsim-models/Solver/MNAVariableCompInterface.h>
#include <dpsim/MNASolver.h>
 
namespace DPsim {
 
enum DirectLinearSolverImpl {
  Undef = 0,
  KLU,
  SparseLU,
  DenseLU,
  CUDADense,
  CUDASparse,
  CUDAMagma,
  Plugin
};

template <typename VarType> class MnaSolverDirect : public MnaSolver<VarType> {
 
protected:
  // #### Data structures for precomputed switch matrices (optionally with parallel frequencies) ####
  std::unordered_map<std::bitset<SWITCH_NUM>, std::vector<SparseMatrix>>
      mSwitchedMatrices;
  std::unordered_map<std::bitset<SWITCH_NUM>,
                     std::vector<std::shared_ptr<DirectLinearSolver>>>
      mDirectLinearSolvers;
 
  // #### Data structures for system recomputation over time ####
  SparseMatrix mBaseSystemMatrix;
  SparseMatrix mVariableSystemMatrix;
  std::shared_ptr<DirectLinearSolver> mDirectLinearSolverVariableSystemMatrix;
  DirectLinearSolverImpl mImplementationInUse;
  DirectLinearSolverConfiguration mConfigurationInUse;

  Bool mVariableComponentChanged = false;
  Bool mSystemMatrixChanged = false;
 
  using MnaSolver<VarType>::mSwitches;
  using MnaSolver<VarType>::mMNAIntfSwitches;
  using MnaSolver<VarType>::mMNAComponents;
  using MnaSolver<VarType>::mVariableComps;
  using MnaSolver<VarType>::mMNAIntfVariableComps;
  using MnaSolver<VarType>::mRightSideVector;
  using MnaSolver<VarType>::mLeftSideVector;
  using MnaSolver<VarType>::mCurrentSwitchStatus;
  using MnaSolver<VarType>::mRightVectorStamps;
  using MnaSolver<VarType>::mNumNetNodes;
  using MnaSolver<VarType>::mNodes;
  using MnaSolver<VarType>::mIsInInitialization;
  using MnaSolver<VarType>::mRightSideVectorHarm;
  using MnaSolver<VarType>::mLeftSideVectorHarm;
  using MnaSolver<VarType>::mFrequencyParallel;
  using MnaSolver<VarType>::mSLog;
  using MnaSolver<VarType>::mSystemMatrixRecomputation;
  using MnaSolver<VarType>::hasVariableComponentChanged;
  using MnaSolver<VarType>::mNumRecomputations;
  using MnaSolver<VarType>::mSyncGen;
  using MnaSolver<VarType>::mFactorizeTimes;
  using MnaSolver<VarType>::mSolveTimes;
  using MnaSolver<VarType>::mRecomputationTimes;
  using MnaSolver<VarType>::mListVariableSystemMatrixEntries;
  using MnaSolver<VarType>::mStateSpaceExtraction;
  using MnaSolver<VarType>::mStateSpaceExtractor;
 
  // #### General
  void createEmptySystemMatrix() override;
 
  // #### Methods for precomputed switch matrices (optionally with parallel frequencies) ####
  void switchedMatrixEmpty(std::size_t index) override;
  void switchedMatrixEmpty(std::size_t swIdx, Int freqIdx) override;
  void switchedMatrixStamp(
      std::size_t index,
      std::vector<std::shared_ptr<CPS::MNAInterface>> &comp) override;
 
  // #### Methods for system recomputation over time ####
  void stampVariableSystemMatrix() override;
  void solveWithSystemMatrixRecomputation(Real time,
                                          Int timeStepCount) override;
  std::shared_ptr<CPS::Task> createSolveTaskRecomp() override;
  virtual void recomputeSystemMatrix(Real time);
  void extractStateSpace(Real time);
 
  // #### Scheduler Task Methods ####
  std::shared_ptr<CPS::Task> createSolveTask() override;
  std::shared_ptr<CPS::Task> createLogTask() override;
  std::shared_ptr<CPS::Task> createSolveTaskHarm(UInt freqIdx) override;
  std::shared_ptr<CPS::Task> createStateSpaceExtractionTask() override;
  void logSystemMatrices() override;
  void solve(Real time, Int timeStepCount) override;
  void solveWithHarmonics(Real time, Int timeStepCount, Int freqIdx) override;

  void logSolveTime();
  void logFactorizationTime();
  void logRecomputationTime();

  std::shared_ptr<DirectLinearSolver>
  createDirectSolverImplementation(CPS::Logger::Log mSLog);
 
public:
  MnaSolverDirect(String name, CPS::Domain domain = CPS::Domain::DP,
                  CPS::Logger::Level logLevel = CPS::Logger::Level::info);

  virtual ~MnaSolverDirect() = default;

  void
  setDirectLinearSolverImplementation(DirectLinearSolverImpl implementation);

  void setDirectLinearSolverConfiguration(
      DirectLinearSolverConfiguration &configuration) override;

  void logLUTimes() override;

  int mIter = 0;
 
  // #### MNA Solver Tasks ####
  class SolveTask : public CPS::Task {
  public:
    SolveTask(MnaSolverDirect<VarType> &solver)
        : Task(solver.mName + ".Solve"), mSolver(solver) {
 
      for (auto it : solver.mMNAComponents) {
        if (it->getRightVector()->get().size() != 0)
          mAttributeDependencies.push_back(it->getRightVector());
      }
      for (auto node : solver.mNodes) {
        mModifiedAttributes.push_back(node->mVoltage);
      }
      mModifiedAttributes.push_back(solver.mLeftSideVector);
    }
 
    void execute(Real time, Int timeStepCount) {
      mSolver.solve(time, timeStepCount);
    }
 
  private:
    MnaSolverDirect<VarType> &mSolver;
  };

  class SolveTaskHarm : public CPS::Task {
  public:
    SolveTaskHarm(MnaSolverDirect<VarType> &solver, UInt freqIdx)
        : Task(solver.mName + ".Solve"), mSolver(solver), mFreqIdx(freqIdx) {
 
      for (auto it : solver.mMNAComponents) {
        if (it->getRightVector()->get().size() != 0)
          mAttributeDependencies.push_back(it->getRightVector());
      }
      for (auto node : solver.mNodes) {
        mModifiedAttributes.push_back(node->mVoltage);
      }
      for (auto leftVec : solver.mLeftSideVectorHarm) {
        mModifiedAttributes.push_back(leftVec);
      }
    }
 
    void execute(Real time, Int timeStepCount) {
      mSolver.solveWithHarmonics(time, timeStepCount, mFreqIdx);
    }
 
  private:
    MnaSolverDirect<VarType> &mSolver;
    UInt mFreqIdx;
  };

  class SolveTaskRecomp : public CPS::Task {
  public:
    SolveTaskRecomp(MnaSolverDirect<VarType> &solver)
        : Task(solver.mName + ".Solve"), mSolver(solver) {
 
      for (auto it : solver.mMNAComponents) {
        if (it->getRightVector()->get().size() != 0)
          mAttributeDependencies.push_back(it->getRightVector());
      }
      for (auto it : solver.mMNAIntfVariableComps) {
        if (it->getRightVector()->get().size() != 0)
          mAttributeDependencies.push_back(it->getRightVector());
      }
      for (auto node : solver.mNodes) {
        mModifiedAttributes.push_back(node->mVoltage);
      }
      mModifiedAttributes.push_back(solver.mLeftSideVector);
    }
 
    void execute(Real time, Int timeStepCount) {
      mSolver.solveWithSystemMatrixRecomputation(time, timeStepCount);
      mSolver.log(time, timeStepCount);
    }
 
  private:
    MnaSolverDirect<VarType> &mSolver;
  };

  class StateSpaceExtractionTask : public CPS::Task {
  public:
    StateSpaceExtractionTask(MnaSolverDirect<VarType> &solver)
        : Task(solver.mName + ".StateSpaceExtraction"), mSolver(solver) {
      mAttributeDependencies.push_back(solver.mLeftSideVector);
      mModifiedAttributes.push_back(Scheduler::external);
    }
 
    void execute(Real time, Int timeStepCount) {
      mSolver.extractStateSpace(time);
    }
 
  private:
    MnaSolverDirect<VarType> &mSolver;
  };

  class LogTask : public CPS::Task {
  public:
    LogTask(MnaSolverDirect<VarType> &solver)
        : Task(solver.mName + ".Log"), mSolver(solver) {
      mAttributeDependencies.push_back(solver.mLeftSideVector);
      mModifiedAttributes.push_back(Scheduler::external);
    }
 
    void execute(Real time, Int timeStepCount) {
      mSolver.log(time, timeStepCount);
    }
 
  private:
    MnaSolverDirect<VarType> &mSolver;
  };
};
} // namespace DPsim