ODESolver.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 <dpsim/Solver.h>
 
#include <dpsim-models/Solver/ODEInterface.h>
 
#include <arkode/arkode.h> // prototypes for ARKode fcts., consts. and includes sundials_types.h
#include <arkode/arkode_direct.h>      // access to ARKDls interface
#include <nvector/nvector_serial.h>    // access to serial N_Vector
#include <sunlinsol/sunlinsol_dense.h> // access to dense SUNLinearSolver
#include <sunmatrix/sunmatrix_dense.h> // access to dense SUNMatrix
 
//using namespace CPS; // led to problems
 
namespace DPsim {
class ODESolver : public Solver {
protected:
  CPS::ODEInterface::Ptr mComponent;

  Int mProbDim;
 
  // ###ARKode-specific variables ###
  void *mArkode_mem{nullptr};
  N_Vector mStates{nullptr};
 
  //for implicit solve:
  bool mImplicitIntegration;
  SUNMatrix A{nullptr};
  SUNLinearSolver LS{nullptr};

  Real mTimestep;
 
  // Same tolerance for each component regardless of system characteristics
  realtype reltol = RCONST(1.0e-6);
  realtype abstol = RCONST(1.0e-10);

  int mFlag{0};
 
  // Similar to DAE-Solver
  CPS::ODEInterface::StSpFn mStSpFunction;
  CPS::ODEInterface::JacFn mJacFunction;

  static int StateSpaceWrapper(realtype t, N_Vector y, N_Vector ydot,
                               void *user_data);
  int StateSpace(realtype t, N_Vector y, N_Vector ydot);
 
  static int JacobianWrapper(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
                             void *user_data, N_Vector tmp1, N_Vector tmp2,
                             N_Vector tmp3);
  int Jacobian(realtype t, N_Vector y, N_Vector fy, SUNMatrix J, N_Vector tmp1,
               N_Vector tmp2, N_Vector tmp3);
  int check_flag(void *flagvalue, const std::string &funcname, int opt);
 
public:
  ODESolver(String name, const CPS::ODEInterface::Ptr &comp,
            bool implicit_integration, Real timestep);
  ~ODESolver();
 
  class SolveTask : public CPS::Task {
  public:
    SolveTask(ODESolver &solver)
        : Task(solver.mName + ".Solve"), mSolver(solver) {
      mAttributeDependencies.push_back(solver.mComponent->mOdePreState);
      mModifiedAttributes.push_back(solver.mComponent->mOdePostState);
    }
 
    void execute(Real time, Int timeStepCount);
 
  private:
    ODESolver &mSolver;
  };
 
  virtual CPS::Task::List getTasks() {
    return CPS::Task::List{std::make_shared<SolveTask>(*this)};
  }

  void initialize();
  Real step(Real initial_time);
};
} // namespace DPsim