diff --git a/GridKit/Model/PowerElectronics/CMakeLists.txt b/GridKit/Model/PowerElectronics/CMakeLists.txt index c985d69d6..66c058cfa 100644 --- a/GridKit/Model/PowerElectronics/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/CMakeLists.txt @@ -23,10 +23,9 @@ add_subdirectory(MicrogridLine) add_subdirectory(MicrogridBusDQ) install( - FILES CircuitComponent.hpp - CircuitNode.hpp - CircuitGraph.hpp - SystemModelPowerElectronics.hpp + FILES Component.hpp + Node.hpp + SystemModel.hpp NodeBase.hpp ExternalConnection.hpp DESTINATION include/GridKit/Model/PowerElectronics) diff --git a/GridKit/Model/PowerElectronics/Capacitor/Capacitor.cpp b/GridKit/Model/PowerElectronics/Capacitor/Capacitor.cpp index 7d53e309f..1fe6e1907 100644 --- a/GridKit/Model/PowerElectronics/Capacitor/Capacitor.cpp +++ b/GridKit/Model/PowerElectronics/Capacitor/Capacitor.cpp @@ -1,5 +1,4 @@ - #include "Capacitor.hpp" #include @@ -8,146 +7,147 @@ namespace GridKit { - - /*! - * @brief Constructor for Capacitor - * - * @todo this needs to be tested on some circuit - * - * Calls default ModelEvaluatorImpl constructor. - */ - - template - Capacitor::Capacitor(IdxT id, RealT C) - : C_(C) - { - size_ = 3; - n_intern_ = 1; - n_extern_ = 2; - extern_indices_ = {0, 1}; - idc_ = id; - nnz_ = 5; - } - - template - Capacitor::~Capacitor() - { - } - - /** - * Initialization of the grid model - */ - template - int Capacitor::initialize() - { - return 0; - } - - /* - * \brief Identify differential variables - */ - template - int Capacitor::tagDifferentiable() - { - // All variables are differentials - std::fill(tag_.begin(), tag_.end(), true); - return 0; - } - - /** - * @brief Compute the absolute tolerance for each variable in the model - * - * @param rel_tol The relative tolerance which can be used to pick the - * absolute tolerance. - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return int 0 if successful, non-zero otherwise. - * - * This represents a "noise" level close to zero for which pure relative - * error cannot be used. - */ - template - int Capacitor::setAbsoluteTolerance(RealT rel_tol) - { - abs_tol_.setToConst(static_cast(rel_tol)); - return 0; - } - - /** - * @brief Evaluate the resisdual of the Capcitor - * - */ - template - int Capacitor::evaluateInternalResidual() - { - f_int_[0] = -C_ * yp_int_[0] + *y_ext_[0] - *y_ext_[1] - y_int_[0]; - return 0; - } - - template - int Capacitor::evaluateExternalResidual() - { - // input - *f_ext_[0] += C_ * yp_int_[0]; - // output - *f_ext_[1] += -C_ * yp_int_[0]; - return 0; - } - - /** - * @brief Compute the Jacobian dF/dy - a dF/dy' - * - * @tparam ScalarT - * @tparam IdxT - * @return int - */ - template - int Capacitor::evaluateJacobian() + namespace PowerElectronics { - this->zeroJacMatrix(); - // Create dF/dy - std::vector rcord{0, 1, 2, 2, 2}; - std::vector ccord{2, 2, 0, 1, 2}; - std::vector vals{C_ * alpha_, -C_ * alpha_, 1.0, -1.0, -1.0 - C_ * alpha_}; - this->setJacValues(rcord, ccord, vals); - - return 0; - } - - template - int Capacitor::evaluateIntegrand() - { - return 0; - } - - template - int Capacitor::initializeAdjoint() - { - return 0; - } - - template - int Capacitor::evaluateAdjointResidual() - { - return 0; - } - - template - int Capacitor::evaluateAdjointIntegrand() - { - return 0; - } - - template - CircuitComponent* Capacitor::clone() const - { - return new Capacitor(*this); - } - - // Available template instantiations - template class Capacitor; - template class Capacitor; - template class Capacitor; - template class Capacitor; - + /*! + * @brief Constructor for Capacitor + * + * @todo this needs to be tested on some circuit + * + * Calls default ModelEvaluatorImpl constructor. + */ + template + Capacitor::Capacitor(IdxT id, RealT C) + : C_(C) + { + size_ = 3; + n_intern_ = 1; + n_extern_ = 2; + extern_indices_ = {0, 1}; + idc_ = id; + nnz_ = 5; + } + + template + Capacitor::~Capacitor() + { + } + + /** + * Initialization of the grid model + */ + template + int Capacitor::initialize() + { + return 0; + } + + /* + * \brief Identify differential variables + */ + template + int Capacitor::tagDifferentiable() + { + // All variables are differentials + std::fill(tag_.begin(), tag_.end(), true); + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Capacitor::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Evaluate the resisdual of the Capcitor + * + */ + template + int Capacitor::evaluateInternalResidual() + { + f_int_[0] = -C_ * yp_int_[0] + *y_ext_[0] - *y_ext_[1] - y_int_[0]; + return 0; + } + + template + int Capacitor::evaluateExternalResidual() + { + // input + *f_ext_[0] += C_ * yp_int_[0]; + // output + *f_ext_[1] += -C_ * yp_int_[0]; + return 0; + } + + /** + * @brief Compute the Jacobian dF/dy - a dF/dy' + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ + template + int Capacitor::evaluateJacobian() + { + this->zeroJacMatrix(); + // Create dF/dy + std::vector rcord{0, 1, 2, 2, 2}; + std::vector ccord{2, 2, 0, 1, 2}; + std::vector vals{C_ * alpha_, -C_ * alpha_, 1.0, -1.0, -1.0 - C_ * alpha_}; + this->setJacValues(rcord, ccord, vals); + + return 0; + } + + template + int Capacitor::evaluateIntegrand() + { + return 0; + } + + template + int Capacitor::initializeAdjoint() + { + return 0; + } + + template + int Capacitor::evaluateAdjointResidual() + { + return 0; + } + + template + int Capacitor::evaluateAdjointIntegrand() + { + return 0; + } + + template + Component* Capacitor::clone() const + { + return new Capacitor(*this); + } + + // Available template instantiations + template class Capacitor; + template class Capacitor; + template class Capacitor; + template class Capacitor; + + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/Capacitor/Capacitor.hpp b/GridKit/Model/PowerElectronics/Capacitor/Capacitor.hpp index 201597fa0..79f58e5b7 100644 --- a/GridKit/Model/PowerElectronics/Capacitor/Capacitor.hpp +++ b/GridKit/Model/PowerElectronics/Capacitor/Capacitor.hpp @@ -1,70 +1,66 @@ - #pragma once -#include - -namespace GridKit -{ - template - class BaseBus; -} +#include namespace GridKit { - /*! - * @brief Declaration of a Capacitor class. - * - */ - template - class Capacitor : public CircuitComponent + namespace PowerElectronics { - using RealT = typename CircuitComponent::RealT; + /*! + * @brief Declaration of a Capacitor class. + * + */ + template + class Capacitor : public Component + { + using RealT = typename Component::RealT; - using CircuitComponent::size_; - using CircuitComponent::nnz_; - using CircuitComponent::time_; - using CircuitComponent::alpha_; - using CircuitComponent::y_ext_; - using CircuitComponent::y_int_; - using CircuitComponent::yp_ext_; - using CircuitComponent::yp_int_; - using CircuitComponent::abs_tol_; - using CircuitComponent::tag_; - using CircuitComponent::f_ext_; - using CircuitComponent::f_int_; - using CircuitComponent::g_; - using CircuitComponent::yB_; - using CircuitComponent::ypB_; - using CircuitComponent::fB_; - using CircuitComponent::gB_; - using CircuitComponent::param_; - using CircuitComponent::idc_; + using Component::size_; + using Component::nnz_; + using Component::time_; + using Component::alpha_; + using Component::y_ext_; + using Component::y_int_; + using Component::yp_ext_; + using Component::yp_int_; + using Component::abs_tol_; + using Component::tag_; + using Component::f_ext_; + using Component::f_int_; + using Component::g_; + using Component::yB_; + using Component::ypB_; + using Component::fB_; + using Component::gB_; + using Component::param_; + using Component::idc_; - using CircuitComponent::extern_indices_; - using CircuitComponent::n_extern_; - using CircuitComponent::n_intern_; + using Component::extern_indices_; + using Component::n_extern_; + using Component::n_intern_; - public: - Capacitor(IdxT id, RealT C); - virtual ~Capacitor(); + public: + Capacitor(IdxT id, RealT C); + virtual ~Capacitor(); - int initialize(); - int tagDifferentiable(); - int setAbsoluteTolerance(RealT); - int evaluateInternalResidual() final; - int evaluateExternalResidual() final; - int evaluateJacobian(); - int evaluateIntegrand(); + int initialize(); + int tagDifferentiable(); + int setAbsoluteTolerance(RealT); + int evaluateInternalResidual() final; + int evaluateExternalResidual() final; + int evaluateJacobian(); + int evaluateIntegrand(); - int initializeAdjoint(); - int evaluateAdjointResidual(); - // int evaluateAdjointJacobian(); - int evaluateAdjointIntegrand(); + int initializeAdjoint(); + int evaluateAdjointResidual(); + // int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); - CircuitComponent* clone() const; + Component* clone() const; - private: - RealT C_; - }; + private: + RealT C_; + }; + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/CircuitComponent.hpp b/GridKit/Model/PowerElectronics/CircuitComponent.hpp deleted file mode 100644 index c3e33792a..000000000 --- a/GridKit/Model/PowerElectronics/CircuitComponent.hpp +++ /dev/null @@ -1,787 +0,0 @@ - - -#pragma once - -#include -#include -#include -#include -#include - -#include -#include -#include - -namespace GridKit -{ - /*! - * @brief Declaration of a CircuitComponent class. - * - */ - template - class CircuitComponent : public Model::Evaluator - { - public: - using RealT = typename Model::Evaluator::RealT; - using CsrMatrixT = typename Model::Evaluator::CsrMatrixT; - using VectorT = typename Model::Evaluator::VectorT; - - CircuitComponent() = default; - - protected: - /** - * @brief Constructs a copy of a circuit component. - * - * Copies the component metadata, local vector data, connection-node mapping, - * and COO Jacobian storage. Dynamically allocated component-owned data is - * deep-copied so that the new component does not share ownership of this - * storage with @p other. - * - * Pointers to state, state-derivative, and residual storage supplied by a - * parent system are copied as-is. Consequently, the copied component initially - * references the same parent-system storage as @p other. The pointer arrays - * used for external variables are independently allocated, but their entries - * point to the same external state, state-derivative, and residual storage as - * the original component. - * - * Local vectors, including the state, state derivative, residual, tolerances, - * quadrature data, adjoint data, and parameter vectors, retain the values of - * the original component. - * - * @param other Component to copy. - * - * @note If the copied component is subsequently attached to a different parent - * system, its state, state-derivative, and residual pointers must be - * reassigned to the storage provided by that system before evaluation. - */ - CircuitComponent(const CircuitComponent& other) - : n_extern_(other.n_extern_), - n_intern_(other.n_intern_), - extern_indices_(other.extern_indices_), - size_(other.size_), - nnz_(other.nnz_), - size_quad_(other.size_quad_), - size_opt_(other.size_opt_), - current_jac_size_(other.current_jac_size_), - y_int_(other.y_int_), - yp_int_(other.yp_int_), - f_int_(other.f_int_), - tag_(other.tag_), - time_(other.time_), - alpha_(other.alpha_), - max_steps_(other.max_steps_), - idc_(other.idc_), - allocated_(other.allocated_) - { - - auto copyVector = [](VectorT& destination, const VectorT& source) - { - const IdxT source_size = source.getSize(); - - if (source_size == 0) - { - return; - } - - destination.resize(source_size); - destination.copyFromExternal(source); - }; - - /* - * Deep-copy the local-to-global connection mapping. - */ - if (other.connection_nodes_) - { - connection_nodes_ = std::make_unique(static_cast(size_)); - - for (size_t i = 0; i < static_cast(size_); ++i) - { - connection_nodes_[i] = other.connection_nodes_[i]; - } - } - - /* - * Deep-copy the COO Jacobian row indices. - */ - if (other.jacobian_coo_rows_) - { - jacobian_coo_rows_ = std::make_unique(static_cast(nnz_)); - - for (size_t i = 0; i < static_cast(nnz_); ++i) - { - jacobian_coo_rows_[i] = other.jacobian_coo_rows_[i]; - } - } - - /* - * Deep-copy the COO Jacobian column indices. - */ - if (other.jacobian_coo_cols_) - { - jacobian_coo_cols_ = std::make_unique(static_cast(nnz_)); - - for (size_t i = 0; i < static_cast(nnz_); ++i) - { - jacobian_coo_cols_[i] = other.jacobian_coo_cols_[i]; - } - } - - /* - * Deep-copy the COO Jacobian values. - */ - if (other.jacobian_coo_values_) - { - jacobian_coo_values_ = std::make_unique(static_cast(nnz_)); - - for (size_t i = 0; i < static_cast(nnz_); ++i) - { - jacobian_coo_values_[i] = other.jacobian_coo_values_[i]; - } - } - - if (size_ > 0) - { - y_ext_ = std::make_unique(static_cast(size_)); - yp_ext_ = std::make_unique(static_cast(size_)); - f_ext_ = std::make_unique(static_cast(size_)); - - for (size_t i = 0; i < static_cast(size_); ++i) - { - y_ext_[i] = other.y_ext_[i]; - yp_ext_[i] = other.yp_ext_[i]; - f_ext_[i] = other.f_ext_[i]; - } - } - - // State, state derivative, residual, and absolute tolerance. - copyVector(y_, other.y_); - copyVector(yp_, other.yp_); - copyVector(f_, other.f_); - copyVector(abs_tol_, other.abs_tol_); - copyVector(g_, other.g_); - copyVector(yB_, other.yB_); - copyVector(ypB_, other.ypB_); - copyVector(fB_, other.fB_); - copyVector(gB_, other.gB_); - copyVector(param_, other.param_); - copyVector(param_up_, other.param_up_); - copyVector(param_lo_, other.param_lo_); - } - - public: - /** - * @brief Create an independent copy of this component. - * - * The clone preserves the component's model configuration, parameters, - * topology, and structural data, but does not preserve bindings to - * system-owned state or residual storage. - * - * @note By default, the cloned component's state, state-derivative, and - * residual pointers are not set. The user is responsible for setting these - * pointers to the appropriate storage before evaluating the residual. - */ - virtual CircuitComponent* clone() const - { - throw std::runtime_error("clone() is not supported for this component."); - } - - /** - * @brief Indicates whether this component supports cloning. - * - * Derived components that implement clone() should override this method - * and return true. - * - * @return true if the component can be cloned, false otherwise. - */ - virtual bool isCloneable() const - { - return false; - } - - /** - * @note Cannot be marked final, since it is overriden to recurse in the system model. - */ - void updateTime(RealT t, RealT a) override - { - this->time_ = t; - this->alpha_ = a; - } - - bool hasJacobian() override - { - return true; - } - - size_t getExternSize() - { - return n_extern_; - } - - size_t getInternalSize() - { - return this->n_intern_; - } - - std::set getExternIndices() - { - return this->extern_indices_; - } - - /** - * @brief Create the mappings from local to global indices for an internal variable. - * Used for constructing system Jacobians \see connection_nodes_. - * - * @param local_index The index of the local variable - * @param global_index The index of the corresponding system variable. - * - * @pre `local_index` *must* be the index of an internal variable. Using this method for - * an external variable will not properly setup the data pointers for that variable. - */ - int setInternalConnectionNodes(size_t local_index, IdxT global_index) - { - assert(!extern_indices_.contains(static_cast(local_index))); - setConnectionNodes(local_index, global_index); - return 0; - } - - /** - * @brief Create the mappings from local to global indices for an external variable. - * External variables need extra information than internal variables - their data - * pointers \ref y_ext_, \ref yp_ext_, and \ref f_ext_. - * - * @param local_index The index of the local variable - * @param connection The necessary connection information for the variable - * - * @pre `local_index` *must* be the index of an external variable. As of now, using this method - * to set information for a local variable will silently discard the unnecessary information, but - * this may change in the future. - */ - int setExternalConnectionNodes(size_t local_index, ExternalConnection connection) - { - assert(extern_indices_.contains(local_index)); - y_ext_[local_index] = connection.y_; - yp_ext_[local_index] = connection.yp_; - f_ext_[local_index] = connection.f_; - setConnectionNodes(local_index, connection.idx_); - return 0; - } - - /** - * @brief Update the connection index for a variable. - * - * Sets only the connection index without modifying the variable's - * internal/external classification or its associated data pointers. - * - * @param local_index Index of the local variable. - * @param connection_index New connection index for the variable. - * - * @return int 0 if successful. - */ - int setConnectionNodes(size_t local_index, IdxT connection_index) - { - connection_nodes_[local_index] = connection_index; - return 0; - } - - /** - * @brief Given the location of value in the local vector map to global index - * - * f(local_index) = global_index - * - * @param local_index index of local value in vector - * @return size_t Index of the same value in the global vector - */ - IdxT getNodeConnection(size_t local_index) const - { - return connection_nodes_[local_index]; - } - - int initialize() override - { - y_.setDataUpdated(); - yp_.setDataUpdated(); - - return 0; - } - - /** - * @brief Allocates all of the internal buffers for the component. - * If a components needs a more specialized allocation (such as by having additional internal buffers), - * it should override this function and then call it in the body to ensure it stays up-to-date with - * new implementations. - * - * @pre \ref nnz_ and \ref size_ must be set. Typically these are set by the child object in its constructor. - * - * @return An error code, or 0 if success - */ - int allocate() override - { - jacobian_coo_rows_ = std::make_unique(static_cast(nnz_)); - jacobian_coo_cols_ = std::make_unique(static_cast(nnz_)); - jacobian_coo_values_ = std::make_unique(static_cast(nnz_)); - - y_ext_ = std::make_unique(static_cast(size_)); - yp_ext_ = std::make_unique(static_cast(size_)); - f_ext_ = std::make_unique(static_cast(size_)); - - connection_nodes_ = std::make_unique(static_cast(size_)); - - tag_.resize(static_cast(size_)); - - if (!allocated_) - { - allocateVectors(size_); - } - - allocated_ = true; - return 0; - } - - IdxT* jacobianCooRows() - { - return jacobian_coo_rows_.get(); - } - - const IdxT* jacobianCooRows() const - { - return jacobian_coo_rows_.get(); - } - - IdxT* jacobianCooCols() - { - return jacobian_coo_cols_.get(); - } - - const IdxT* jacobianCooCols() const - { - return jacobian_coo_cols_.get(); - } - - RealT* jacobianCooValues() - { - return jacobian_coo_values_.get(); - } - - const RealT* jacobianCooValues() const - { - return jacobian_coo_values_.get(); - } - - /** - * @brief Evaluating the residual of a CircuitComponent should be done by evaluating the - * internal residuals and external residuals. CircuitComponents should overload those - * functions for their residuals (and the system will call those function instead of this one), - * so there is no reason to overload this functionality. - * - * @return An error code, or 0 is successful. - */ - int evaluateResidual() final - { - if (int err_code = evaluateInternalResidual()) - return err_code; - - f_.setDataUpdated(); - - return evaluateExternalResidual(); - } - - /** - * @brief Evaluate all residuals for the component's internal variables, - * writing them through `f_int_`. - * - * @return An error code, or 0 if successful. - */ - virtual int evaluateInternalResidual() = 0; - - /** - * @brief Evaluate all residual contributions for the component's external variables, - * accumulating them through `f_ext_`. - * - * @return An error code, or 0 if successful. - */ - virtual int evaluateExternalResidual() = 0; - - void setInternalPointer(const ScalarT* internals) - { - y_int_ = internals; - } - - void setInternalDerivativePointer(const ScalarT* internals_p) - { - yp_int_ = internals_p; - } - - void setInternalResidualPointer(ScalarT* internal_res) - { - f_int_ = internal_res; - } - - protected: - /** - * @brief Reset the Jacobian so it can be constructed. Helper method for \ref setJacValues(). - * Sets \ref current_jac_size_ to 0 so that future calls to `setJacValues()` will override previous values. - * - */ - void zeroJacMatrix() - { - current_jac_size_ = 0; - } - - /** - * @brief Helper method for adding values to the Jacobian. Copies the rows, cols, and vals buffers and appends - * them to the end of the corresponding Jacobian buffers. Uses \ref current_jac_size_ to tell where the end of - * the Jacobian currently is. - * - * @pre `rows`, `cols`, `vals` must all be the same size - * @pre \ref allocate() must be called first. - * @pre Must call \ref zeroJacMatrix() before starting construction of a new Jacobian - * @pre The must be enough room for the values in the allocated buffers, i.e. `current_jac_size_ + rows.size() <= nnz_` - */ - void setJacValues(const std::vector& rows, const std::vector& cols, const std::vector& vals) - { - assert(rows.size() == cols.size()); - assert(rows.size() == vals.size()); - assert(current_jac_size_ + rows.size() <= static_cast(nnz_)); - - for (size_t i = 0; i < rows.size(); i++) - { - jacobian_coo_rows_[current_jac_size_] = rows[i]; - jacobian_coo_cols_[current_jac_size_] = cols[i]; - jacobian_coo_values_[current_jac_size_] = vals[i]; - - current_jac_size_++; - } - } - - public: - IdxT size() final - { - return size_; - } - - IdxT size() const - { - return size_; - } - - IdxT nnz() final - { - return nnz_; - } - - IdxT nnz() const - { - return nnz_; - } - - IdxT sizeQuadrature() final - { - return size_quad_; - } - - IdxT sizeQuadrature() const - { - return size_quad_; - } - - IdxT sizeParams() final - { - return size_opt_; - } - - IdxT sizeParams() const - { - return size_opt_; - } - - VectorT& y() final - { - return y_; - } - - const VectorT& y() const final - { - return y_; - } - - VectorT& yp() final - { - return yp_; - } - - const VectorT& yp() const final - { - return yp_; - } - - std::vector& tag() final - { - return tag_; - } - - const std::vector& tag() const final - { - return tag_; - } - - VectorT& absoluteTolerance() final - { - return abs_tol_; - } - - const VectorT& absoluteTolerance() const final - { - return abs_tol_; - } - - VectorT& yB() final - { - return yB_; - } - - const VectorT& yB() const final - { - return yB_; - } - - VectorT& ypB() final - { - return ypB_; - } - - const VectorT& ypB() const final - { - return ypB_; - } - - VectorT& param() final - { - return param_; - } - - const VectorT& param() const final - { - return param_; - } - - VectorT& param_up() final - { - return param_up_; - } - - const VectorT& param_up() const final - { - return param_up_; - } - - VectorT& param_lo() final - { - return param_lo_; - } - - const VectorT& param_lo() const final - { - return param_lo_; - } - - VectorT& getResidual() final - { - return f_; - } - - const VectorT& getResidual() const final - { - return f_; - } - - VectorT& getIntegrand() final - { - return g_; - } - - const VectorT& getIntegrand() const final - { - return g_; - } - - VectorT& getAdjointResidual() final - { - return fB_; - } - - const VectorT& getAdjointResidual() const final - { - return fB_; - } - - VectorT& getAdjointIntegrand() final - { - return gB_; - } - - const VectorT& getAdjointIntegrand() const final - { - return gB_; - } - - //@todo Fix ID naming - IdxT getIDcomponent() const - { - return idc_; - } - - /** - * @brief Check whether the component has already been allocated. - * - * @return true if allocate() has previously completed, false otherwise. - */ - bool isAllocated() const - { - return allocated_; - } - - protected: - /** - * @brief Allocate state and residual storage owned by this component. - * - * Most components do not need state and residual storages. The most notable exception - * is currently the system, so a separate flag is provided for the system. - * Systems still can't directly access \ref y_, \ref yp_, and \ref f_, so they need - * their corresponding \ref y_int_, \ref yp_int_, and \ref f_int_ set, since there isn't - * another system above them to set it. - * - * @todo This is a weird exception specifically for systems - and in a hierarchical setting - * will only be needed by the *topmost* system - subsystems shouldn't allocate and should have their - * internal pointers set by the system above them. Ideally we can remove this exception by having - * the integrator allocate these buffers instead of the system and set the internal pointers for the - * topmost system. - */ - void allocateVectors(IdxT n, bool system = false) - { - abs_tol_.resize(n); - - if (system) - { - y_.resize(n); - yp_.resize(n); - f_.resize(n); - - y_int_ = y_.getData(); - yp_int_ = yp_.getData(); - f_int_ = f_.getData(); - } - } - - /// Number of external variables in this component - ones which are referenced but not owned by this component. - size_t n_extern_; - /// Number of internal variables in this component - ones which are only referenced by this component. - size_t n_intern_; - /** - * @brief A set of variable indices which correspond to the external variables. Variables indices not in this set are internal. - * - * @invariant Must have a size of \ref n_extern_. Each element must be in the range [0, \ref size_ - 1]. Not currently verified anywhere. - */ - std::set extern_indices_; - /** - * @brief A map from local variable indices to system (global) variable indices. Used for Jacobian construction in - * \ref PowerElectronicsModel::evaluateJacobian(). - * @note If a variable does not map to a corresponding variable in the system (such as with reference nodes), a special - * sentinel value of \ref INVALID_INDEX is used. During Jacobian construction, such rows and columns will be pruned. - */ - std::unique_ptr connection_nodes_; - - /// The number of variables in this component. Should be equal to \ref n_extern_ plus \ref n_intern_. \see size() - IdxT size_{0}; - /// The number of nonzero elements in this component's Jacobian. \see nnz() - IdxT nnz_{0}; - IdxT size_quad_{0}; - IdxT size_opt_{0}; - - // COO Jacobian buffers - std::unique_ptr jacobian_coo_rows_; - std::unique_ptr jacobian_coo_cols_; - std::unique_ptr jacobian_coo_values_; - - /// The number of non-zero elements currently inserted into the Jacobian. See \ref setJacValues() - size_t current_jac_size_{0}; - - /// @brief A pointer to the internal variables of this component. - const ScalarT* y_int_; - /// @brief A pointer to the internal derivatives of this component. - const ScalarT* yp_int_; - /// @brief A pointer to the internal residuals of this component - ScalarT* f_int_; - - /** - * An array of (input) pointers to state values for external variables. - * \note The size of this array is equal to \ref size_, allowing you to index it with the index - * of the variable in question (i.e. consisten with \ref extern_indices_). Therefore, accessing - * and dereferencing the pointer in an internal variable index is undefined behavior. - * \see setExternalConnectionNodes() - */ - std::unique_ptr y_ext_; - /** - * An array of (input) pointers to derivative values for external variables. - * \note The size of this array is equal to \ref size_, allowing you to index it with the index - * of the variable in question (i.e. consisten with \ref extern_indices_). Therefore, accessing - * and dereferencing the pointer in an internal variable index is undefined behavior. - * \see setExternalConnectionNodes() - */ - std::unique_ptr yp_ext_; - /** - * An array of (output) pointers to residuals for external variables. - * \note The size of this array is equal to \ref size_, allowing you to index it with the index - * of the variable in question (i.e. consisten with \ref extern_indices_). Therefore, accessing - * and dereferencing the pointer in an internal variable index is undefined behavior. - * \see setExternalConnectionNodes() - */ - std::unique_ptr f_ext_; - - std::vector tag_; - VectorT abs_tol_; - - VectorT g_; - - VectorT yB_; - VectorT ypB_; - VectorT fB_; - VectorT gB_; - - VectorT param_; - VectorT param_up_; - VectorT param_lo_; - - RealT time_; - RealT alpha_; - - IdxT max_steps_; - - IdxT idc_; - - bool allocated_{false}; - - private: - /** - * The internal buffer for state for the component. For most components, it will be empty and shouldn't be accessed. - * Instead use \ref y_int_ for an internal variable or \ref y_ext_ for an external variable, respectively. - * For components which want an internal buffer (such as a system), make sure that \ref y_int_ points here. - * \see allocateVectors() - */ - VectorT y_; - /** - * The internal buffer for derivatives for the component. For most components, it will be empty and shouldn't be accessed. - * Instead use \ref yp_int_ for an internal variable or \ref yp_ext_ for an external variable, respectively. - * For components which want an internal buffer (such as a system), make sure that \ref yp_int_ points here. - * \see allocateVectors() - */ - VectorT yp_; - /** - * The internal buffer for state for the component. For most components, it will be empty and shouldn't be accessed. - * Instead use \ref f_int_ for an internal variable or \ref f_ext_ for an external variable, respectively. - * For components which want an internal buffer (such as a system), make sure that \ref f_int_ points here. - * \see allocateVectors() - */ - VectorT f_; - }; - -} // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/CircuitGraph.hpp b/GridKit/Model/PowerElectronics/CircuitGraph.hpp deleted file mode 100644 index 5a57104fb..000000000 --- a/GridKit/Model/PowerElectronics/CircuitGraph.hpp +++ /dev/null @@ -1,123 +0,0 @@ - - -#include -#include -#include -#include -#include - -/** - * @brief A very basic hypergraph setup for circuit representation. - * This forms the hypergraph as a bipartite graph. Doesn't allow - * removing. Can only grab sets of connections to nodes - * - * @todo should replace with something better and more efficent. - * Should replace with a libraries setup instead. This would allow - * fast and easy partitioning of circuits - * - * @todo This is to replace inserting vector size for allocating PowerElectronicsModel - * - * @todo should replace N and E with Node and Component classes respectively. - * - * @note Tested but currently not used in the rest of the code. - * - * @tparam IdxT - * @tparam Label - */ -template -class CircuitGraph -{ -private: - std::set hypernodes; - std::set hyperedges; - std::map> edgestonodes; - -public: - CircuitGraph(); - ~CircuitGraph(); - bool addHyperEdge(E he); - bool addHyperNode(N hn); - bool addConnection(N hn, E he); - std::set getHyperEdgeConnections(E he); - size_t amountHyperNodes(); - size_t amountHyperEdges(); - void printBiPartiteGraph(bool verbose = false); -}; - -template -CircuitGraph::CircuitGraph() -{ -} - -template -CircuitGraph::~CircuitGraph() -{ -} - -template -bool CircuitGraph::addHyperNode(N hn) -{ - return this->hypernodes.insert(hn).second; -} - -template -bool CircuitGraph::addHyperEdge(E he) -{ - return this->hyperedges.insert(he).second; -} - -template -bool CircuitGraph::addConnection(N hn, E he) -{ - if (this->hyperedges.count(he) == 0 || this->hypernodes.count(hn) == 0) - { - return false; - } - return this->edgestonodes[he].insert(hn).second; -} - -template -std::set CircuitGraph::getHyperEdgeConnections(E he) -{ - return this->edgestonodes[he]; -} - -template -size_t CircuitGraph::amountHyperNodes() -{ - return this->hypernodes.size(); -} - -template -size_t CircuitGraph::amountHyperEdges() -{ - return this->hyperedges.size(); -} - -/** - * @brief Print the bipartite graph - * - * @todo need to add verbose printing for connections display - * - * @tparam IdxT - * @param[in] verbose if true will print connections, - * otherwise just the number of nodes and edges - */ - -template -void CircuitGraph::printBiPartiteGraph([[maybe_unused]] bool verbose) -{ - - std::cout << "Amount of HyperNodes: " << this->amountHyperNodes() << std::endl; - std::cout << "Amount of HyperEdges: " << this->amountHyperEdges() << std::endl; - std::cout << "Connections per Edge:" << std::endl; - for (auto i : this->edgestonodes) - { - std::cout << i.first << " : {"; - for (auto j : i.second) - { - std::cout << j << ", "; - } - std::cout << "}\n"; - } -} diff --git a/GridKit/Model/PowerElectronics/CircuitNode.hpp b/GridKit/Model/PowerElectronics/CircuitNode.hpp deleted file mode 100644 index 7b1e52bf8..000000000 --- a/GridKit/Model/PowerElectronics/CircuitNode.hpp +++ /dev/null @@ -1,382 +0,0 @@ -#pragma once - -#include -#include -#include -#include - -#include - -namespace GridKit -{ - /** - * @brief Circuit node representing a connection point. - */ - template - class CircuitNode : public Model::Evaluator - { - using RealT = typename Model::Evaluator::RealT; - using VectorT = typename Model::Evaluator::VectorT; - - public: - CircuitNode() - { - size_ = 1; - } - - CircuitNode(ScalarT v0) - : V0_(v0) - { - size_ = 1; - } - - ~CircuitNode() = default; - - int setNodeID(IdxT id) - { - id_ = id; - return 0; - } - - IdxT nodeID() const - { - return id_; - } - - // Voltage accessor - ScalarT& V() - { - return y_.getData()[0]; - } - - const ScalarT& V() const - { - return y_.getData()[0]; - } - - // KCL residual accessor - ScalarT& I() - { - return f_.getData()[0]; - } - - const ScalarT& I() const - { - return f_.getData()[0]; - } - - // Allocate storage for a single-node voltage and KCL residual - int allocate() - { - size_t size = static_cast(size_); - - if (!allocated_) - { - allocateVectors(size_); - } - - tag_.resize(size); - - variable_indices_[0] = 0; - residual_indices_[0] = 0; - - allocated_ = true; - return 0; - } - - /** - * @brief Initialize node variables - */ - int initialize() - { - auto* y = y_.getData(); - auto* yp = yp_.getData(); - - y[0] = V0_; - yp[0] = 0.0; - - y_.setDataUpdated(); - yp_.setDataUpdated(); - - return 0; - } - - /** - * @brief Node variables are algebraic. - */ - int tagDifferentiable() - { - tag_[0] = false; - - return 0; - } - - /** - * @brief Compute the absolute tolerance for each variable in the model - * - * @param rel_tol The relative tolerance which can be used to pick the - * absolute tolerance. - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return int 0 if successful, non-zero otherwise. - * - * This represents a "noise" level close to zero for which pure relative - * error cannot be used. - */ - int setAbsoluteTolerance(RealT rel_tol) - { - abs_tol_.setToConst(static_cast(rel_tol)); - return 0; - } - - /** - * @brief Node does not compute residuals, so here we just reset residual values. - * - * @warning This implementation assumes node residuals are always evaluated - * _before_ component model residuals. - * - */ - int evaluateResidual() - { - auto* f = f_.getData(); - - f[0] = 0.0; - - f_.setDataUpdated(); - - return 0; - } - - bool hasJacobian() final - { - return false; - } - - /** - * @brief There is no Jacobian for node variables - */ - int evaluateJacobian() - { - return 0; - } - - int evaluateIntegrand() - { - return 0; - } - - int initializeAdjoint() - { - return 0; - } - - int evaluateAdjointResidual() - { - return 0; - } - - int evaluateAdjointIntegrand() - { - return 0; - } - - private: - IdxT id_{static_cast(-1)}; - IdxT size_{0}; - IdxT nnz_{0}; - IdxT size_quad_{0}; - IdxT size_opt_{0}; - ScalarT V0_{0.0}; - - std::map variable_indices_; - std::map residual_indices_; - - VectorT y_; - VectorT yp_; - std::vector tag_; - VectorT abs_tol_; - VectorT f_; - - VectorT g_{}; - VectorT param_{}; - VectorT param_up_{}; - VectorT param_lo_{}; - - VectorT yB_{}; - VectorT ypB_{}; - VectorT fB_{}; - VectorT gB_{}; - - RealT time_{0}; - RealT alpha_{0}; - - IdxT max_steps_{0}; - - bool allocated_{false}; - - public: - IdxT size() final - { - return size_; - } - - IdxT nnz() final - { - return nnz_; - } - - IdxT sizeQuadrature() final - { - return size_quad_; - } - - IdxT sizeParams() final - { - return size_opt_; - } - - void updateTime(RealT /* t */, RealT /* a */) final - { - // No time to update in node models - } - - VectorT& y() final - { - return y_; - } - - const VectorT& y() const final - { - return y_; - } - - VectorT& yp() final - { - return yp_; - } - - const VectorT& yp() const final - { - return yp_; - } - - std::vector& tag() final - { - return tag_; - } - - const std::vector& tag() const final - { - return tag_; - } - - VectorT& absoluteTolerance() final - { - return abs_tol_; - } - - const VectorT& absoluteTolerance() const final - { - return abs_tol_; - } - - VectorT& yB() final - { - return yB_; - } - - const VectorT& yB() const final - { - return yB_; - } - - VectorT& ypB() final - { - return ypB_; - } - - const VectorT& ypB() const final - { - return ypB_; - } - - VectorT& param() final - { - return param_; - } - - const VectorT& param() const final - { - return param_; - } - - VectorT& param_up() final - { - return param_up_; - } - - const VectorT& param_up() const final - { - return param_up_; - } - - VectorT& param_lo() final - { - return param_lo_; - } - - const VectorT& param_lo() const final - { - return param_lo_; - } - - VectorT& getResidual() final - { - return f_; - } - - const VectorT& getResidual() const final - { - return f_; - } - - VectorT& getIntegrand() final - { - return g_; - } - - const VectorT& getIntegrand() const final - { - return g_; - } - - VectorT& getAdjointResidual() final - { - return fB_; - } - - const VectorT& getAdjointResidual() const final - { - return fB_; - } - - VectorT& getAdjointIntegrand() final - { - return gB_; - } - - const VectorT& getAdjointIntegrand() const final - { - return gB_; - } - - private: - void allocateVectors(IdxT n) - { - y_.resize(n); - yp_.resize(n); - f_.resize(n); - abs_tol_.resize(n); - } - }; -} // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/Component.hpp b/GridKit/Model/PowerElectronics/Component.hpp new file mode 100644 index 000000000..aa7c5c715 --- /dev/null +++ b/GridKit/Model/PowerElectronics/Component.hpp @@ -0,0 +1,788 @@ + +#pragma once + +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace GridKit +{ + namespace PowerElectronics + { + /*! + * @brief Declaration of a Component class. + * + */ + template + class Component : public Model::Evaluator + { + public: + using RealT = typename Model::Evaluator::RealT; + using CsrMatrixT = typename Model::Evaluator::CsrMatrixT; + using VectorT = typename Model::Evaluator::VectorT; + + Component() = default; + + protected: + /** + * @brief Constructs a copy of a circuit component. + * + * Copies the component metadata, local vector data, connection-node mapping, + * and COO Jacobian storage. Dynamically allocated component-owned data is + * deep-copied so that the new component does not share ownership of this + * storage with @p other. + * + * Pointers to state, state-derivative, and residual storage supplied by a + * parent system are copied as-is. Consequently, the copied component initially + * references the same parent-system storage as @p other. The pointer arrays + * used for external variables are independently allocated, but their entries + * point to the same external state, state-derivative, and residual storage as + * the original component. + * + * Local vectors, including the state, state derivative, residual, tolerances, + * quadrature data, adjoint data, and parameter vectors, retain the values of + * the original component. + * + * @param other Component to copy. + * + * @note If the copied component is subsequently attached to a different parent + * system, its state, state-derivative, and residual pointers must be + * reassigned to the storage provided by that system before evaluation. + */ + Component(const Component& other) + : n_extern_(other.n_extern_), + n_intern_(other.n_intern_), + extern_indices_(other.extern_indices_), + size_(other.size_), + nnz_(other.nnz_), + size_quad_(other.size_quad_), + size_opt_(other.size_opt_), + current_jac_size_(other.current_jac_size_), + y_int_(other.y_int_), + yp_int_(other.yp_int_), + f_int_(other.f_int_), + tag_(other.tag_), + time_(other.time_), + alpha_(other.alpha_), + max_steps_(other.max_steps_), + idc_(other.idc_), + allocated_(other.allocated_) + { + + auto copyVector = [](VectorT& destination, const VectorT& source) + { + const IdxT source_size = source.getSize(); + + if (source_size == 0) + { + return; + } + + destination.resize(source_size); + destination.copyFromExternal(source); + }; + + /* + * Deep-copy the local-to-global connection mapping. + */ + if (other.connection_nodes_) + { + connection_nodes_ = std::make_unique(static_cast(size_)); + + for (size_t i = 0; i < static_cast(size_); ++i) + { + connection_nodes_[i] = other.connection_nodes_[i]; + } + } + + /* + * Deep-copy the COO Jacobian row indices. + */ + if (other.jacobian_coo_rows_) + { + jacobian_coo_rows_ = std::make_unique(static_cast(nnz_)); + + for (size_t i = 0; i < static_cast(nnz_); ++i) + { + jacobian_coo_rows_[i] = other.jacobian_coo_rows_[i]; + } + } + + /* + * Deep-copy the COO Jacobian column indices. + */ + if (other.jacobian_coo_cols_) + { + jacobian_coo_cols_ = std::make_unique(static_cast(nnz_)); + + for (size_t i = 0; i < static_cast(nnz_); ++i) + { + jacobian_coo_cols_[i] = other.jacobian_coo_cols_[i]; + } + } + + /* + * Deep-copy the COO Jacobian values. + */ + if (other.jacobian_coo_values_) + { + jacobian_coo_values_ = std::make_unique(static_cast(nnz_)); + + for (size_t i = 0; i < static_cast(nnz_); ++i) + { + jacobian_coo_values_[i] = other.jacobian_coo_values_[i]; + } + } + + if (size_ > 0) + { + y_ext_ = std::make_unique(static_cast(size_)); + yp_ext_ = std::make_unique(static_cast(size_)); + f_ext_ = std::make_unique(static_cast(size_)); + + for (size_t i = 0; i < static_cast(size_); ++i) + { + y_ext_[i] = other.y_ext_[i]; + yp_ext_[i] = other.yp_ext_[i]; + f_ext_[i] = other.f_ext_[i]; + } + } + + // State, state derivative, residual, and absolute tolerance. + copyVector(y_, other.y_); + copyVector(yp_, other.yp_); + copyVector(f_, other.f_); + copyVector(abs_tol_, other.abs_tol_); + copyVector(g_, other.g_); + copyVector(yB_, other.yB_); + copyVector(ypB_, other.ypB_); + copyVector(fB_, other.fB_); + copyVector(gB_, other.gB_); + copyVector(param_, other.param_); + copyVector(param_up_, other.param_up_); + copyVector(param_lo_, other.param_lo_); + } + + public: + /** + * @brief Create an independent copy of this component. + * + * The clone preserves the component's model configuration, parameters, + * topology, and structural data, but does not preserve bindings to + * system-owned state or residual storage. + * + * @note By default, the cloned component's state, state-derivative, and + * residual pointers are not set. The user is responsible for setting these + * pointers to the appropriate storage before evaluating the residual. + */ + virtual Component* clone() const + { + throw std::runtime_error("clone() is not supported for this component."); + } + + /** + * @brief Indicates whether this component supports cloning. + * + * Derived components that implement clone() should override this method + * and return true. + * + * @return true if the component can be cloned, false otherwise. + */ + virtual bool isCloneable() const + { + return false; + } + + /** + * @note Cannot be marked final, since it is overriden to recurse in the system model. + */ + void updateTime(RealT t, RealT a) override + { + this->time_ = t; + this->alpha_ = a; + } + + bool hasJacobian() override + { + return true; + } + + size_t getExternSize() + { + return n_extern_; + } + + size_t getInternalSize() + { + return this->n_intern_; + } + + std::set getExternIndices() + { + return this->extern_indices_; + } + + /** + * @brief Create the mappings from local to global indices for an internal variable. + * Used for constructing system Jacobians \see connection_nodes_. + * + * @param local_index The index of the local variable + * @param global_index The index of the corresponding system variable. + * + * @pre `local_index` *must* be the index of an internal variable. Using this method for + * an external variable will not properly setup the data pointers for that variable. + */ + int setInternalConnectionNodes(size_t local_index, IdxT global_index) + { + assert(!extern_indices_.contains(static_cast(local_index))); + setConnectionNodes(local_index, global_index); + return 0; + } + + /** + * @brief Create the mappings from local to global indices for an external variable. + * External variables need extra information than internal variables - their data + * pointers \ref y_ext_, \ref yp_ext_, and \ref f_ext_. + * + * @param local_index The index of the local variable + * @param connection The necessary connection information for the variable + * + * @pre `local_index` *must* be the index of an external variable. As of now, using this method + * to set information for a local variable will silently discard the unnecessary information, but + * this may change in the future. + */ + int setExternalConnectionNodes(size_t local_index, ExternalConnection connection) + { + assert(extern_indices_.contains(local_index)); + y_ext_[local_index] = connection.y_; + yp_ext_[local_index] = connection.yp_; + f_ext_[local_index] = connection.f_; + setConnectionNodes(local_index, connection.idx_); + return 0; + } + + /** + * @brief Update the connection index for a variable. + * + * Sets only the connection index without modifying the variable's + * internal/external classification or its associated data pointers. + * + * @param local_index Index of the local variable. + * @param connection_index New connection index for the variable. + * + * @return int 0 if successful. + */ + int setConnectionNodes(size_t local_index, IdxT connection_index) + { + connection_nodes_[local_index] = connection_index; + return 0; + } + + /** + * @brief Given the location of value in the local vector map to global index + * + * f(local_index) = global_index + * + * @param local_index index of local value in vector + * @return size_t Index of the same value in the global vector + */ + IdxT getNodeConnection(size_t local_index) const + { + return connection_nodes_[local_index]; + } + + int initialize() override + { + y_.setDataUpdated(); + yp_.setDataUpdated(); + + return 0; + } + + /** + * @brief Allocates all of the internal buffers for the component. + * If a components needs a more specialized allocation (such as by having additional internal buffers), + * it should override this function and then call it in the body to ensure it stays up-to-date with + * new implementations. + * + * @pre \ref nnz_ and \ref size_ must be set. Typically these are set by the child object in its constructor. + * + * @return An error code, or 0 if success + */ + int allocate() override + { + jacobian_coo_rows_ = std::make_unique(static_cast(nnz_)); + jacobian_coo_cols_ = std::make_unique(static_cast(nnz_)); + jacobian_coo_values_ = std::make_unique(static_cast(nnz_)); + + y_ext_ = std::make_unique(static_cast(size_)); + yp_ext_ = std::make_unique(static_cast(size_)); + f_ext_ = std::make_unique(static_cast(size_)); + + connection_nodes_ = std::make_unique(static_cast(size_)); + + tag_.resize(static_cast(size_)); + + if (!allocated_) + { + allocateVectors(size_); + } + + allocated_ = true; + return 0; + } + + IdxT* jacobianCooRows() + { + return jacobian_coo_rows_.get(); + } + + const IdxT* jacobianCooRows() const + { + return jacobian_coo_rows_.get(); + } + + IdxT* jacobianCooCols() + { + return jacobian_coo_cols_.get(); + } + + const IdxT* jacobianCooCols() const + { + return jacobian_coo_cols_.get(); + } + + RealT* jacobianCooValues() + { + return jacobian_coo_values_.get(); + } + + const RealT* jacobianCooValues() const + { + return jacobian_coo_values_.get(); + } + + /** + * @brief Evaluating the residual of a Component should be done by evaluating the + * internal residuals and external residuals. Components should overload those + * functions for their residuals (and the system will call those function instead of this one), + * so there is no reason to overload this functionality. + * + * @return An error code, or 0 is successful. + */ + int evaluateResidual() final + { + if (int err_code = evaluateInternalResidual()) + return err_code; + + f_.setDataUpdated(); + + return evaluateExternalResidual(); + } + + /** + * @brief Evaluate all residuals for the component's internal variables, + * writing them through `f_int_`. + * + * @return An error code, or 0 if successful. + */ + virtual int evaluateInternalResidual() = 0; + + /** + * @brief Evaluate all residual contributions for the component's external variables, + * accumulating them through `f_ext_`. + * + * @return An error code, or 0 if successful. + */ + virtual int evaluateExternalResidual() = 0; + + void setInternalPointer(const ScalarT* internals) + { + y_int_ = internals; + } + + void setInternalDerivativePointer(const ScalarT* internals_p) + { + yp_int_ = internals_p; + } + + void setInternalResidualPointer(ScalarT* internal_res) + { + f_int_ = internal_res; + } + + protected: + /** + * @brief Reset the Jacobian so it can be constructed. Helper method for \ref setJacValues(). + * Sets \ref current_jac_size_ to 0 so that future calls to `setJacValues()` will override previous values. + * + */ + void zeroJacMatrix() + { + current_jac_size_ = 0; + } + + /** + * @brief Helper method for adding values to the Jacobian. Copies the rows, cols, and vals buffers and appends + * them to the end of the corresponding Jacobian buffers. Uses \ref current_jac_size_ to tell where the end of + * the Jacobian currently is. + * + * @pre `rows`, `cols`, `vals` must all be the same size + * @pre \ref allocate() must be called first. + * @pre Must call \ref zeroJacMatrix() before starting construction of a new Jacobian + * @pre The must be enough room for the values in the allocated buffers, i.e. `current_jac_size_ + rows.size() <= nnz_` + */ + void setJacValues(const std::vector& rows, const std::vector& cols, const std::vector& vals) + { + assert(rows.size() == cols.size()); + assert(rows.size() == vals.size()); + assert(current_jac_size_ + rows.size() <= static_cast(nnz_)); + + for (size_t i = 0; i < rows.size(); i++) + { + jacobian_coo_rows_[current_jac_size_] = rows[i]; + jacobian_coo_cols_[current_jac_size_] = cols[i]; + jacobian_coo_values_[current_jac_size_] = vals[i]; + + current_jac_size_++; + } + } + + public: + IdxT size() final + { + return size_; + } + + IdxT size() const + { + return size_; + } + + IdxT nnz() final + { + return nnz_; + } + + IdxT nnz() const + { + return nnz_; + } + + IdxT sizeQuadrature() final + { + return size_quad_; + } + + IdxT sizeQuadrature() const + { + return size_quad_; + } + + IdxT sizeParams() final + { + return size_opt_; + } + + IdxT sizeParams() const + { + return size_opt_; + } + + VectorT& y() final + { + return y_; + } + + const VectorT& y() const final + { + return y_; + } + + VectorT& yp() final + { + return yp_; + } + + const VectorT& yp() const final + { + return yp_; + } + + std::vector& tag() final + { + return tag_; + } + + const std::vector& tag() const final + { + return tag_; + } + + VectorT& absoluteTolerance() final + { + return abs_tol_; + } + + const VectorT& absoluteTolerance() const final + { + return abs_tol_; + } + + VectorT& yB() final + { + return yB_; + } + + const VectorT& yB() const final + { + return yB_; + } + + VectorT& ypB() final + { + return ypB_; + } + + const VectorT& ypB() const final + { + return ypB_; + } + + VectorT& param() final + { + return param_; + } + + const VectorT& param() const final + { + return param_; + } + + VectorT& param_up() final + { + return param_up_; + } + + const VectorT& param_up() const final + { + return param_up_; + } + + VectorT& param_lo() final + { + return param_lo_; + } + + const VectorT& param_lo() const final + { + return param_lo_; + } + + VectorT& getResidual() final + { + return f_; + } + + const VectorT& getResidual() const final + { + return f_; + } + + VectorT& getIntegrand() final + { + return g_; + } + + const VectorT& getIntegrand() const final + { + return g_; + } + + VectorT& getAdjointResidual() final + { + return fB_; + } + + const VectorT& getAdjointResidual() const final + { + return fB_; + } + + VectorT& getAdjointIntegrand() final + { + return gB_; + } + + const VectorT& getAdjointIntegrand() const final + { + return gB_; + } + + //@todo Fix ID naming + IdxT getIDcomponent() const + { + return idc_; + } + + /** + * @brief Check whether the component has already been allocated. + * + * @return true if allocate() has previously completed, false otherwise. + */ + bool isAllocated() const + { + return allocated_; + } + + protected: + /** + * @brief Allocate state and residual storage owned by this component. + * + * Most components do not need state and residual storages. The most notable exception + * is currently the system, so a separate flag is provided for the system. + * Systems still can't directly access \ref y_, \ref yp_, and \ref f_, so they need + * their corresponding \ref y_int_, \ref yp_int_, and \ref f_int_ set, since there isn't + * another system above them to set it. + * + * @todo This is a weird exception specifically for systems - and in a hierarchical setting + * will only be needed by the *topmost* system - subsystems shouldn't allocate and should have their + * internal pointers set by the system above them. Ideally we can remove this exception by having + * the integrator allocate these buffers instead of the system and set the internal pointers for the + * topmost system. + */ + void allocateVectors(IdxT n, bool system = false) + { + abs_tol_.resize(n); + + if (system) + { + y_.resize(n); + yp_.resize(n); + f_.resize(n); + + y_int_ = y_.getData(); + yp_int_ = yp_.getData(); + f_int_ = f_.getData(); + } + } + + /// Number of external variables in this component - ones which are referenced but not owned by this component. + size_t n_extern_; + /// Number of internal variables in this component - ones which are only referenced by this component. + size_t n_intern_; + /** + * @brief A set of variable indices which correspond to the external variables. Variables indices not in this set are internal. + * + * @invariant Must have a size of \ref n_extern_. Each element must be in the range [0, \ref size_ - 1]. Not currently verified anywhere. + */ + std::set extern_indices_; + /** + * @brief A map from local variable indices to system (global) variable indices. Used for Jacobian construction in + * \ref PowerElectronics::SystemModel::evaluateJacobian(). + * @note If a variable does not map to a corresponding variable in the system (such as with reference nodes), a special + * sentinel value of \ref INVALID_INDEX is used. During Jacobian construction, such rows and columns will be pruned. + */ + std::unique_ptr connection_nodes_; + + /// The number of variables in this component. Should be equal to \ref n_extern_ plus \ref n_intern_. \see size() + IdxT size_{0}; + /// The number of nonzero elements in this component's Jacobian. \see nnz() + IdxT nnz_{0}; + IdxT size_quad_{0}; + IdxT size_opt_{0}; + + // COO Jacobian buffers + std::unique_ptr jacobian_coo_rows_; + std::unique_ptr jacobian_coo_cols_; + std::unique_ptr jacobian_coo_values_; + + /// The number of non-zero elements currently inserted into the Jacobian. See \ref setJacValues() + size_t current_jac_size_{0}; + + /// @brief A pointer to the internal variables of this component. + const ScalarT* y_int_; + /// @brief A pointer to the internal derivatives of this component. + const ScalarT* yp_int_; + /// @brief A pointer to the internal residuals of this component + ScalarT* f_int_; + + /** + * An array of (input) pointers to state values for external variables. + * \note The size of this array is equal to \ref size_, allowing you to index it with the index + * of the variable in question (i.e. consisten with \ref extern_indices_). Therefore, accessing + * and dereferencing the pointer in an internal variable index is undefined behavior. + * \see setExternalConnectionNodes() + */ + std::unique_ptr y_ext_; + /** + * An array of (input) pointers to derivative values for external variables. + * \note The size of this array is equal to \ref size_, allowing you to index it with the index + * of the variable in question (i.e. consisten with \ref extern_indices_). Therefore, accessing + * and dereferencing the pointer in an internal variable index is undefined behavior. + * \see setExternalConnectionNodes() + */ + std::unique_ptr yp_ext_; + /** + * An array of (output) pointers to residuals for external variables. + * \note The size of this array is equal to \ref size_, allowing you to index it with the index + * of the variable in question (i.e. consisten with \ref extern_indices_). Therefore, accessing + * and dereferencing the pointer in an internal variable index is undefined behavior. + * \see setExternalConnectionNodes() + */ + std::unique_ptr f_ext_; + + std::vector tag_; + VectorT abs_tol_; + + VectorT g_; + + VectorT yB_; + VectorT ypB_; + VectorT fB_; + VectorT gB_; + + VectorT param_; + VectorT param_up_; + VectorT param_lo_; + + RealT time_; + RealT alpha_; + + IdxT max_steps_; + + IdxT idc_; + + bool allocated_{false}; + + private: + /** + * The internal buffer for state for the component. For most components, it will be empty and shouldn't be accessed. + * Instead use \ref y_int_ for an internal variable or \ref y_ext_ for an external variable, respectively. + * For components which want an internal buffer (such as a system), make sure that \ref y_int_ points here. + * \see allocateVectors() + */ + VectorT y_; + /** + * The internal buffer for derivatives for the component. For most components, it will be empty and shouldn't be accessed. + * Instead use \ref yp_int_ for an internal variable or \ref yp_ext_ for an external variable, respectively. + * For components which want an internal buffer (such as a system), make sure that \ref yp_int_ points here. + * \see allocateVectors() + */ + VectorT yp_; + /** + * The internal buffer for state for the component. For most components, it will be empty and shouldn't be accessed. + * Instead use \ref f_int_ for an internal variable or \ref f_ext_ for an external variable, respectively. + * For components which want an internal buffer (such as a system), make sure that \ref f_int_ points here. + * \see allocateVectors() + */ + VectorT f_; + }; + } // namespace PowerElectronics +} // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp index db4a280a5..44f7914da 100644 --- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp +++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp @@ -1,5 +1,4 @@ - #include "DistributedGenerator.hpp" #include @@ -7,431 +6,433 @@ namespace GridKit { - - /*! - * @brief Constructor for a Distributed Generator - * @todo Maybe have parameters be templated in. Variables cannot be changed - * and are unlikely to. Allows for compile time optimizations - * - * Calls default ModelEvaluatorImpl constructor. - */ - template - DistributedGenerator::DistributedGenerator(IdxT id, - DistributedGeneratorParameters parm, - bool reference_frame, - NodeT* node_ref, - NodeT* node_bus) - : wb_(parm.wb_), - wc_(parm.wc_), - mp_(parm.mp_), - Vn_(parm.Vn_), - nq_(parm.nq_), - F_(parm.F_), - Kiv_(parm.Kiv_), - Kpv_(parm.Kpv_), - Kic_(parm.Kic_), - Kpc_(parm.Kpc_), - Cf_(parm.Cf_), - rLf_(parm.rLf_), - Lf_(parm.Lf_), - rLc_(parm.rLc_), - Lc_(parm.Lc_), - refframe_(reference_frame), - node_ref_(node_ref), - node_bus_(node_bus) - { - assert(refframe_ || node_ref_->size() == 1); - assert(node_bus_->size() == 2); - // internals [Pi, Qi, phi_di, phi_qi, gamma_di, gamma_qi, il_di, il_qi, vo_di, vo_qi, io_di, io_qi, \delta_i] - // externals [\omega_ref, vba_out, vbb_out] - size_ = refframe_ ? 15 : 16; - n_intern_ = refframe_ ? 12 : 13; - n_extern_ = 3; - idc_ = id; - nnz_ = refframe_ ? 73 : 78; - - extern_indices_ = {0, 1, 2}; - } - - template - DistributedGenerator::~DistributedGenerator() + namespace PowerElectronics { - } + /*! + * @brief Constructor for a Distributed Generator + * @todo Maybe have parameters be templated in. Variables cannot be changed + * and are unlikely to. Allows for compile time optimizations + * + * Calls default ModelEvaluatorImpl constructor. + */ + template + DistributedGenerator::DistributedGenerator(IdxT id, + DistributedGeneratorParameters parm, + bool reference_frame, + NodeT* node_ref, + NodeT* node_bus) + : wb_(parm.wb_), + wc_(parm.wc_), + mp_(parm.mp_), + Vn_(parm.Vn_), + nq_(parm.nq_), + F_(parm.F_), + Kiv_(parm.Kiv_), + Kpv_(parm.Kpv_), + Kic_(parm.Kic_), + Kpc_(parm.Kpc_), + Cf_(parm.Cf_), + rLf_(parm.rLf_), + Lf_(parm.Lf_), + rLc_(parm.rLc_), + Lc_(parm.Lc_), + refframe_(reference_frame), + node_ref_(node_ref), + node_bus_(node_bus) + { + assert(refframe_ || node_ref_->size() == 1); + assert(node_bus_->size() == 2); + // internals [Pi, Qi, phi_di, phi_qi, gamma_di, gamma_qi, il_di, il_qi, vo_di, vo_qi, io_di, io_qi, \delta_i] + // externals [\omega_ref, vba_out, vbb_out] + size_ = refframe_ ? 15 : 16; + n_intern_ = refframe_ ? 12 : 13; + n_extern_ = 3; + idc_ = id; + nnz_ = refframe_ ? 73 : 78; + + extern_indices_ = {0, 1, 2}; + } - /** - * Initialization of the grid model - */ - template - int DistributedGenerator::initialize() - { - return 0; - } - - /* - * \brief Identify differential variables - */ - template - int DistributedGenerator::tagDifferentiable() - { - // All variables are differentials - std::fill(tag_.begin(), tag_.end(), true); - return 0; - } - - /** - * @brief Compute the absolute tolerance for each variable in the model - * - * @param rel_tol The relative tolerance which can be used to pick the - * absolute tolerance. - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return int 0 if successful, non-zero otherwise. - * - * This represents a "noise" level close to zero for which pure relative - * error cannot be used. - */ - template - int DistributedGenerator::setAbsoluteTolerance(RealT rel_tol) - { - abs_tol_.setToConst(static_cast(rel_tol)); - return 0; - } - - /** - * @brief Contributes to the resisdual of the Distributed Generator - * - */ - template - int DistributedGenerator::evaluateInternalResidual() - { - ScalarT omega = wb_ - mp_ * y_int_[0]; - ScalarT delta = refframe_ ? ScalarT(0.0) : y_int_[12]; + template + DistributedGenerator::~DistributedGenerator() + { + } - // Take incoming voltages to current rotator reference frame - ScalarT vbd_in = std::cos(delta) * *y_ext_[1] + std::sin(delta) * *y_ext_[2]; - ScalarT vbq_in = -std::sin(delta) * *y_ext_[1] + std::cos(delta) * *y_ext_[2]; + /** + * Initialization of the grid model + */ + template + int DistributedGenerator::initialize() + { + return 0; + } - // ### Internal Componenets ## - // P and Q equations - f_int_[0] = -yp_int_[0] + wc_ * (y_int_[8] * y_int_[10] + y_int_[9] * y_int_[11] - y_int_[0]); - f_int_[1] = -yp_int_[1] + wc_ * (-y_int_[8] * y_int_[11] + y_int_[9] * y_int_[10] - y_int_[1]); + /* + * \brief Identify differential variables + */ + template + int DistributedGenerator::tagDifferentiable() + { + // All variables are differentials + std::fill(tag_.begin(), tag_.end(), true); + return 0; + } - // Voltage control - ScalarT vod_star = Vn_ - nq_ * y_int_[1]; - ScalarT voq_star = static_cast(0.0); + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int DistributedGenerator::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } - f_int_[2] = -yp_int_[2] + vod_star - y_int_[8]; - f_int_[3] = -yp_int_[3] + voq_star - y_int_[9]; + /** + * @brief Contributes to the resisdual of the Distributed Generator + * + */ + template + int DistributedGenerator::evaluateInternalResidual() + { + ScalarT omega = wb_ - mp_ * y_int_[0]; + ScalarT delta = refframe_ ? ScalarT(0.0) : y_int_[12]; - ScalarT ild_star = F_ * y_int_[10] - wb_ * Cf_ * y_int_[9] + Kpv_ * (vod_star - y_int_[8]) + Kiv_ * y_int_[2]; - ScalarT ilq_star = F_ * y_int_[11] + wb_ * Cf_ * y_int_[8] + Kpv_ * (voq_star - y_int_[9]) + Kiv_ * y_int_[3]; + // Take incoming voltages to current rotator reference frame + ScalarT vbd_in = std::cos(delta) * *y_ext_[1] + std::sin(delta) * *y_ext_[2]; + ScalarT vbq_in = -std::sin(delta) * *y_ext_[1] + std::cos(delta) * *y_ext_[2]; - // Current control - f_int_[4] = -yp_int_[4] + ild_star - y_int_[6]; - f_int_[5] = -yp_int_[5] + ilq_star - y_int_[7]; + // ### Internal Componenets ## + // P and Q equations + f_int_[0] = -yp_int_[0] + wc_ * (y_int_[8] * y_int_[10] + y_int_[9] * y_int_[11] - y_int_[0]); + f_int_[1] = -yp_int_[1] + wc_ * (-y_int_[8] * y_int_[11] + y_int_[9] * y_int_[10] - y_int_[1]); - ScalarT vid_star = -wb_ * Lf_ * y_int_[7] + Kpc_ * (ild_star - y_int_[6]) + Kic_ * y_int_[4]; - ScalarT viq_star = wb_ * Lf_ * y_int_[6] + Kpc_ * (ilq_star - y_int_[7]) + Kic_ * y_int_[5]; + // Voltage control + ScalarT vod_star = Vn_ - nq_ * y_int_[1]; + ScalarT voq_star = static_cast(0.0); - // Output LC Filter - f_int_[6] = -yp_int_[6] - (rLf_ / Lf_) * y_int_[6] + omega * y_int_[7] + (vid_star - y_int_[8]) / Lf_; - f_int_[7] = -yp_int_[7] - (rLf_ / Lf_) * y_int_[7] - omega * y_int_[6] + (viq_star - y_int_[9]) / Lf_; + f_int_[2] = -yp_int_[2] + vod_star - y_int_[8]; + f_int_[3] = -yp_int_[3] + voq_star - y_int_[9]; - f_int_[8] = -yp_int_[8] + omega * y_int_[9] + (y_int_[6] - y_int_[10]) / Cf_; - f_int_[9] = -yp_int_[9] - omega * y_int_[8] + (y_int_[7] - y_int_[11]) / Cf_; + ScalarT ild_star = F_ * y_int_[10] - wb_ * Cf_ * y_int_[9] + Kpv_ * (vod_star - y_int_[8]) + Kiv_ * y_int_[2]; + ScalarT ilq_star = F_ * y_int_[11] + wb_ * Cf_ * y_int_[8] + Kpv_ * (voq_star - y_int_[9]) + Kiv_ * y_int_[3]; - // Output Connector - f_int_[10] = -yp_int_[10] - (rLc_ / Lc_) * y_int_[10] + omega * y_int_[11] + (y_int_[8] - vbd_in) / Lc_; - f_int_[11] = -yp_int_[11] - (rLc_ / Lc_) * y_int_[11] - omega * y_int_[10] + (y_int_[9] - vbq_in) / Lc_; + // Current control + f_int_[4] = -yp_int_[4] + ild_star - y_int_[6]; + f_int_[5] = -yp_int_[5] + ilq_star - y_int_[7]; - // Rotor difference angle - if (!refframe_) - { - f_int_[12] = -yp_int_[12] + omega - *y_ext_[0]; - } + ScalarT vid_star = -wb_ * Lf_ * y_int_[7] + Kpc_ * (ild_star - y_int_[6]) + Kic_ * y_int_[4]; + ScalarT viq_star = wb_ * Lf_ * y_int_[6] + Kpc_ * (ilq_star - y_int_[7]) + Kic_ * y_int_[5]; - return 0; - } + // Output LC Filter + f_int_[6] = -yp_int_[6] - (rLf_ / Lf_) * y_int_[6] + omega * y_int_[7] + (vid_star - y_int_[8]) / Lf_; + f_int_[7] = -yp_int_[7] - (rLf_ / Lf_) * y_int_[7] - omega * y_int_[6] + (viq_star - y_int_[9]) / Lf_; - template - int DistributedGenerator::evaluateExternalResidual() - { - ScalarT omega = wb_ - mp_ * y_int_[0]; - ScalarT delta = refframe_ ? ScalarT(0.0) : y_int_[12]; + f_int_[8] = -yp_int_[8] + omega * y_int_[9] + (y_int_[6] - y_int_[10]) / Cf_; + f_int_[9] = -yp_int_[9] - omega * y_int_[8] + (y_int_[7] - y_int_[11]) / Cf_; + + // Output Connector + f_int_[10] = -yp_int_[10] - (rLc_ / Lc_) * y_int_[10] + omega * y_int_[11] + (y_int_[8] - vbd_in) / Lc_; + f_int_[11] = -yp_int_[11] - (rLc_ / Lc_) * y_int_[11] - omega * y_int_[10] + (y_int_[9] - vbq_in) / Lc_; - // ref common ref motor angle - if (refframe_) + // Rotor difference angle + if (!refframe_) + { + f_int_[12] = -yp_int_[12] + omega - *y_ext_[0]; + } + + return 0; + } + + template + int DistributedGenerator::evaluateExternalResidual() { - *f_ext_[0] += omega - *y_ext_[0]; + ScalarT omega = wb_ - mp_ * y_int_[0]; + ScalarT delta = refframe_ ? ScalarT(0.0) : y_int_[12]; + + // ref common ref motor angle + if (refframe_) + { + *f_ext_[0] += omega - *y_ext_[0]; + } + + // output + // current transformed to common frame + *f_ext_[1] += std::cos(delta) * y_int_[10] - std::sin(delta) * y_int_[11]; + *f_ext_[2] += std::sin(delta) * y_int_[10] + std::cos(delta) * y_int_[11]; + return 0; } - // output - // current transformed to common frame - *f_ext_[1] += std::cos(delta) * y_int_[10] - std::sin(delta) * y_int_[11]; - *f_ext_[2] += std::sin(delta) * y_int_[10] + std::cos(delta) * y_int_[11]; - return 0; - } - - /** - * @brief Compute the jacobian of the DistributedGenerator for iteration. dF/dy - \alpha dF/dy' - * - * The matrix dF/dy should be - * - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] - [-1, 0, 0, 0, -mp, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] - [ 0, 0, 0, 0, -wc, 0, 0, 0, 0, 0, 0, 0, wc*x15, wc*x16, wc*x13, wc*x14] - [ 0, 0, 0, 0, 0, -wc, 0, 0, 0, 0, 0, 0, -wc*x16, wc*x15, wc*x14, -wc*x13] - [ 0, 0, 0, 0, 0, -nq, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0] - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0] - [ 0, 0, 0, 0, 0, -Kpv*nq, Kiv, 0, 0, 0, -1, 0, -Kpv, -Cf*wb, F, 0] - [ 0, 0, 0, 0, 0, 0, 0, Kiv, 0, 0, 0, -1, Cf*wb, -Kpv, 0, F] - [ 0, 0, 0, 0, -mp*x12, -(Kpc*Kpv*nq)/Lf, (Kiv*Kpc)/Lf, 0, Kic/Lf, 0, - Kpc/Lf - rLf/Lf, -mp*x5, -(Kpc*Kpv + 1)/Lf, -(Cf*Kpc*wb)/Lf, (F*Kpc)/Lf, 0] - [ 0, 0, 0, 0, mp*x11, 0, 0, (Kiv*Kpc)/Lf, 0, Kic/Lf, mp*x5, - Kpc/Lf - rLf/Lf, (Cf*Kpc*wb)/Lf, -(Kpc*Kpv + 1)/Lf, 0, (F*Kpc)/Lf] - [ 0, 0, 0, 0, -mp*x14, 0, 0, 0, 0, 0, 1/Cf, 0, 0, wb - mp*x5, -1/Cf, 0] - [ 0, 0, 0, 0, mp*x13, 0, 0, 0, 0, 0, 0, 1/Cf, mp*x5 - wb, 0, 0, -1/Cf] - [ 0, -cos(x4)/Lc, -sin(x4)/Lc, -(x3*cos(x4) - x2*sin(x4))/Lc, -mp*x16, 0, 0, 0, 0, 0, 0, 0, 1/Lc, 0, -rLc/Lc, wb - mp*x5] - [ 0, sin(x4)/Lc, -cos(x4)/Lc, (x2*cos(x4) + x3*sin(x4))/Lc, mp*x15, 0, 0, 0, 0, 0, 0, 0, 0, 1/Lc, mp*x5 - wb, -rLc/Lc] - * 'Generated from MATLAB symbolic' - * - * @tparam ScalarT - * @tparam IdxT - * @return int - */ - template - int DistributedGenerator::evaluateJacobian() - { - this->zeroJacMatrix(); + /** + * @brief Compute the jacobian of the DistributedGenerator for iteration. dF/dy - \alpha dF/dy' + * + * The matrix dF/dy should be + * + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + [-1, 0, 0, 0, -mp, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + [ 0, 0, 0, 0, -wc, 0, 0, 0, 0, 0, 0, 0, wc*x15, wc*x16, wc*x13, wc*x14] + [ 0, 0, 0, 0, 0, -wc, 0, 0, 0, 0, 0, 0, -wc*x16, wc*x15, wc*x14, -wc*x13] + [ 0, 0, 0, 0, 0, -nq, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0] + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0] + [ 0, 0, 0, 0, 0, -Kpv*nq, Kiv, 0, 0, 0, -1, 0, -Kpv, -Cf*wb, F, 0] + [ 0, 0, 0, 0, 0, 0, 0, Kiv, 0, 0, 0, -1, Cf*wb, -Kpv, 0, F] + [ 0, 0, 0, 0, -mp*x12, -(Kpc*Kpv*nq)/Lf, (Kiv*Kpc)/Lf, 0, Kic/Lf, 0, - Kpc/Lf - rLf/Lf, -mp*x5, -(Kpc*Kpv + 1)/Lf, -(Cf*Kpc*wb)/Lf, (F*Kpc)/Lf, 0] + [ 0, 0, 0, 0, mp*x11, 0, 0, (Kiv*Kpc)/Lf, 0, Kic/Lf, mp*x5, - Kpc/Lf - rLf/Lf, (Cf*Kpc*wb)/Lf, -(Kpc*Kpv + 1)/Lf, 0, (F*Kpc)/Lf] + [ 0, 0, 0, 0, -mp*x14, 0, 0, 0, 0, 0, 1/Cf, 0, 0, wb - mp*x5, -1/Cf, 0] + [ 0, 0, 0, 0, mp*x13, 0, 0, 0, 0, 0, 0, 1/Cf, mp*x5 - wb, 0, 0, -1/Cf] + [ 0, -cos(x4)/Lc, -sin(x4)/Lc, -(x3*cos(x4) - x2*sin(x4))/Lc, -mp*x16, 0, 0, 0, 0, 0, 0, 0, 1/Lc, 0, -rLc/Lc, wb - mp*x5] + [ 0, sin(x4)/Lc, -cos(x4)/Lc, (x2*cos(x4) + x3*sin(x4))/Lc, mp*x15, 0, 0, 0, 0, 0, 0, 0, 0, 1/Lc, mp*x5 - wb, -rLc/Lc] + * 'Generated from MATLAB symbolic' + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ + template + int DistributedGenerator::evaluateJacobian() + { + this->zeroJacMatrix(); + + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + + RealT delta = refframe_ ? RealT(0.0) : static_cast(y_int_[12]); + + // r = 0 + if (refframe_) + { + rtemp.assign(2, 0); + ctemp = {0, 3}; + valtemp = {-1.0, -mp_}; + this->setJacValues(rtemp, ctemp, valtemp); + } + + // r = 1 + rtemp.assign(2, 1); + ctemp = {13, 14}; + valtemp = {std::cos(delta), -std::sin(delta)}; + if (!refframe_) + { + rtemp.push_back(1); + ctemp.push_back(15); + valtemp.push_back(-std::sin(delta) * static_cast(y_int_[10]) + - std::cos(delta) * static_cast(y_int_[11])); + } + this->setJacValues(rtemp, ctemp, valtemp); - std::vector ctemp{}; - std::vector rtemp{}; - std::vector valtemp{}; + // r = 2 + rtemp.assign(2, 2); + ctemp = {13, 14}; + valtemp = {std::sin(delta), std::cos(delta)}; + if (!refframe_) + { + rtemp.push_back(2); + ctemp.push_back(15); + valtemp.push_back(std::cos(delta) * static_cast(y_int_[10]) + - std::sin(delta) * static_cast(y_int_[11])); + } + this->setJacValues(rtemp, ctemp, valtemp); - RealT delta = refframe_ ? RealT(0.0) : static_cast(y_int_[12]); + // r = 3 + rtemp.assign(5, 3); + ctemp = {3, 11, 12, 13, 14}; + valtemp = {-wc_ - alpha_, + wc_ * static_cast(y_int_[10]), + wc_ * static_cast(y_int_[11]), + wc_ * static_cast(y_int_[8]), + wc_ * static_cast(y_int_[9])}; + this->setJacValues(rtemp, ctemp, valtemp); - // r = 0 - if (refframe_) - { - rtemp.assign(2, 0); - ctemp = {0, 3}; - valtemp = {-1.0, -mp_}; + // r = 4 + rtemp.assign(5, 4); + ctemp = {4, 11, 12, 13, 14}; + valtemp = {-wc_ - alpha_, + -wc_ * static_cast(y_int_[11]), + wc_ * static_cast(y_int_[10]), + wc_ * static_cast(y_int_[9]), + -wc_ * static_cast(y_int_[8])}; + this->setJacValues(rtemp, ctemp, valtemp); + + // r = 5 + rtemp.assign(3, 5); + ctemp = {4, 5, 11}; + valtemp = {-nq_, -alpha_, -1.0}; + this->setJacValues(rtemp, ctemp, valtemp); + + // r = 6 + rtemp.assign(2, 6); + ctemp = {6, 12}; + valtemp = {-alpha_, -1.0}; + this->setJacValues(rtemp, ctemp, valtemp); + + // r = 7 + rtemp.assign(7, 7); + ctemp = {4, 5, 7, 9, 11, 12, 13}; + valtemp = {-Kpv_ * nq_, Kiv_, -alpha_, -1.0, -Kpv_, -Cf_ * wb_, F_}; this->setJacValues(rtemp, ctemp, valtemp); + + // r = 8 + rtemp.assign(6, 8); + ctemp = {6, 8, 10, 11, 12, 14}; + valtemp = {Kiv_, -alpha_, -1.0, Cf_ * wb_, -Kpv_, F_}; + this->setJacValues(rtemp, ctemp, valtemp); + + // r = 9 + rtemp.assign(9, 9); + ctemp = {3, 4, 5, 7, 9, 10, 11, 12, 13}; + valtemp = {-mp_ * static_cast(y_int_[7]), + -(Kpc_ * Kpv_ * nq_) / Lf_, + (Kpc_ * Kiv_) / Lf_, + Kic_ / Lf_, + -(Kpc_ + rLf_) / Lf_ - alpha_, + -mp_ * static_cast(y_int_[0]), + -(Kpc_ * Kpv_ + 1.0) / Lf_, + -(Cf_ * Kpc_ * wb_) / Lf_, + (F_ * Kpc_) / Lf_}; + this->setJacValues(rtemp, ctemp, valtemp); + + // r = 10 + rtemp.assign(8, 10); + ctemp = {3, 6, 8, 9, 10, 11, 12, 14}; + valtemp = {mp_ * static_cast(y_int_[6]), + (Kiv_ * Kpc_) / Lf_, + Kic_ / Lf_, + mp_ * static_cast(y_int_[0]), + -(Kpc_ + rLf_) / Lf_ - alpha_, + (Cf_ * Kpc_ * wb_) / Lf_, + -(Kpc_ * Kpv_ + 1.0) / Lf_, + (F_ * Kpc_) / Lf_}; + this->setJacValues(rtemp, ctemp, valtemp); + + // r = 11 + rtemp.assign(5, 11); + ctemp = {3, 9, 11, 12, 13}; + valtemp = {-mp_ * static_cast(y_int_[9]), + 1.0 / Cf_, + -alpha_, + wb_ - mp_ * static_cast(y_int_[0]), + -1.0 / Cf_}; + this->setJacValues(rtemp, ctemp, valtemp); + + // r = 12 + rtemp.assign(5, 12); + ctemp = {3, 10, 11, 12, 14}; + valtemp = {mp_ * static_cast(y_int_[8]), + 1.0 / Cf_, + -wb_ + mp_ * static_cast(y_int_[0]), + -alpha_, + -1.0 / Cf_}; + this->setJacValues(rtemp, ctemp, valtemp); + + // r = 13 + rtemp.assign(6, 13); + ctemp = {1, 2, 3, 11, 13, 14}; + valtemp = { + (1.0 / Lc_) * -std::cos(delta), + (1.0 / Lc_) * -std::sin(delta), + -mp_ * static_cast(y_int_[11]), + 1.0 / Lc_, + -rLc_ / Lc_ - alpha_, + wb_ - mp_ * static_cast(y_int_[0]), + }; + if (!refframe_) + { + rtemp.push_back(13); + ctemp.push_back(15); + valtemp.push_back((1.0 / Lc_) * (std::sin(delta) * static_cast(*y_ext_[1]) - std::cos(delta) * static_cast(*y_ext_[2]))); + } + this->setJacValues(rtemp, ctemp, valtemp); + + // r = 14 + rtemp.assign(6, 14); + ctemp = {1, 2, 3, 12, 13, 14}; + valtemp = { + (1.0 / Lc_) * std::sin(delta), + (1.0 / Lc_) * -std::cos(delta), + mp_ * static_cast(y_int_[10]), + 1.0 / Lc_, + -wb_ + mp_ * static_cast(y_int_[0]), + -rLc_ / Lc_ - alpha_, + }; + if (!refframe_) + { + rtemp.push_back(14); + ctemp.push_back(15); + valtemp.push_back((1.0 / Lc_) * (std::cos(delta) * static_cast(*y_ext_[1]) + std::sin(delta) * static_cast(*y_ext_[2]))); + } + this->setJacValues(rtemp, ctemp, valtemp); + + // r = 15 + if (!refframe_) + { + rtemp.assign(3, 15); + ctemp = {0, 3, 15}; + valtemp = {-1.0, -mp_, -alpha_}; + this->setJacValues(rtemp, ctemp, valtemp); + } + + return 0; } - // r = 1 - rtemp.assign(2, 1); - ctemp = {13, 14}; - valtemp = {std::cos(delta), -std::sin(delta)}; - if (!refframe_) + template + int DistributedGenerator::allocate() { - rtemp.push_back(1); - ctemp.push_back(15); - valtemp.push_back(-std::sin(delta) * static_cast(y_int_[10]) - - std::cos(delta) * static_cast(y_int_[11])); + Component::allocate(); + + this->setExternalConnectionNodes(0, node_ref_->getNodeConnection(0)); + this->setExternalConnectionNodes(1, node_bus_->getNodeConnection(0)); + this->setExternalConnectionNodes(2, node_bus_->getNodeConnection(1)); + + return 0; } - this->setJacValues(rtemp, ctemp, valtemp); - // r = 2 - rtemp.assign(2, 2); - ctemp = {13, 14}; - valtemp = {std::sin(delta), std::cos(delta)}; - if (!refframe_) + template + int DistributedGenerator::evaluateIntegrand() { - rtemp.push_back(2); - ctemp.push_back(15); - valtemp.push_back(std::cos(delta) * static_cast(y_int_[10]) - - std::sin(delta) * static_cast(y_int_[11])); + return 0; } - this->setJacValues(rtemp, ctemp, valtemp); - - // r = 3 - rtemp.assign(5, 3); - ctemp = {3, 11, 12, 13, 14}; - valtemp = {-wc_ - alpha_, - wc_ * static_cast(y_int_[10]), - wc_ * static_cast(y_int_[11]), - wc_ * static_cast(y_int_[8]), - wc_ * static_cast(y_int_[9])}; - this->setJacValues(rtemp, ctemp, valtemp); - - // r = 4 - rtemp.assign(5, 4); - ctemp = {4, 11, 12, 13, 14}; - valtemp = {-wc_ - alpha_, - -wc_ * static_cast(y_int_[11]), - wc_ * static_cast(y_int_[10]), - wc_ * static_cast(y_int_[9]), - -wc_ * static_cast(y_int_[8])}; - this->setJacValues(rtemp, ctemp, valtemp); - - // r = 5 - rtemp.assign(3, 5); - ctemp = {4, 5, 11}; - valtemp = {-nq_, -alpha_, -1.0}; - this->setJacValues(rtemp, ctemp, valtemp); - - // r = 6 - rtemp.assign(2, 6); - ctemp = {6, 12}; - valtemp = {-alpha_, -1.0}; - this->setJacValues(rtemp, ctemp, valtemp); - - // r = 7 - rtemp.assign(7, 7); - ctemp = {4, 5, 7, 9, 11, 12, 13}; - valtemp = {-Kpv_ * nq_, Kiv_, -alpha_, -1.0, -Kpv_, -Cf_ * wb_, F_}; - this->setJacValues(rtemp, ctemp, valtemp); - - // r = 8 - rtemp.assign(6, 8); - ctemp = {6, 8, 10, 11, 12, 14}; - valtemp = {Kiv_, -alpha_, -1.0, Cf_ * wb_, -Kpv_, F_}; - this->setJacValues(rtemp, ctemp, valtemp); - - // r = 9 - rtemp.assign(9, 9); - ctemp = {3, 4, 5, 7, 9, 10, 11, 12, 13}; - valtemp = {-mp_ * static_cast(y_int_[7]), - -(Kpc_ * Kpv_ * nq_) / Lf_, - (Kpc_ * Kiv_) / Lf_, - Kic_ / Lf_, - -(Kpc_ + rLf_) / Lf_ - alpha_, - -mp_ * static_cast(y_int_[0]), - -(Kpc_ * Kpv_ + 1.0) / Lf_, - -(Cf_ * Kpc_ * wb_) / Lf_, - (F_ * Kpc_) / Lf_}; - this->setJacValues(rtemp, ctemp, valtemp); - - // r = 10 - rtemp.assign(8, 10); - ctemp = {3, 6, 8, 9, 10, 11, 12, 14}; - valtemp = {mp_ * static_cast(y_int_[6]), - (Kiv_ * Kpc_) / Lf_, - Kic_ / Lf_, - mp_ * static_cast(y_int_[0]), - -(Kpc_ + rLf_) / Lf_ - alpha_, - (Cf_ * Kpc_ * wb_) / Lf_, - -(Kpc_ * Kpv_ + 1.0) / Lf_, - (F_ * Kpc_) / Lf_}; - this->setJacValues(rtemp, ctemp, valtemp); - - // r = 11 - rtemp.assign(5, 11); - ctemp = {3, 9, 11, 12, 13}; - valtemp = {-mp_ * static_cast(y_int_[9]), - 1.0 / Cf_, - -alpha_, - wb_ - mp_ * static_cast(y_int_[0]), - -1.0 / Cf_}; - this->setJacValues(rtemp, ctemp, valtemp); - - // r = 12 - rtemp.assign(5, 12); - ctemp = {3, 10, 11, 12, 14}; - valtemp = {mp_ * static_cast(y_int_[8]), - 1.0 / Cf_, - -wb_ + mp_ * static_cast(y_int_[0]), - -alpha_, - -1.0 / Cf_}; - this->setJacValues(rtemp, ctemp, valtemp); - - // r = 13 - rtemp.assign(6, 13); - ctemp = {1, 2, 3, 11, 13, 14}; - valtemp = { - (1.0 / Lc_) * -std::cos(delta), - (1.0 / Lc_) * -std::sin(delta), - -mp_ * static_cast(y_int_[11]), - 1.0 / Lc_, - -rLc_ / Lc_ - alpha_, - wb_ - mp_ * static_cast(y_int_[0]), - }; - if (!refframe_) + + template + int DistributedGenerator::initializeAdjoint() { - rtemp.push_back(13); - ctemp.push_back(15); - valtemp.push_back((1.0 / Lc_) * (std::sin(delta) * static_cast(*y_ext_[1]) - std::cos(delta) * static_cast(*y_ext_[2]))); + return 0; } - this->setJacValues(rtemp, ctemp, valtemp); - - // r = 14 - rtemp.assign(6, 14); - ctemp = {1, 2, 3, 12, 13, 14}; - valtemp = { - (1.0 / Lc_) * std::sin(delta), - (1.0 / Lc_) * -std::cos(delta), - mp_ * static_cast(y_int_[10]), - 1.0 / Lc_, - -wb_ + mp_ * static_cast(y_int_[0]), - -rLc_ / Lc_ - alpha_, - }; - if (!refframe_) + + template + int DistributedGenerator::evaluateAdjointResidual() { - rtemp.push_back(14); - ctemp.push_back(15); - valtemp.push_back((1.0 / Lc_) * (std::cos(delta) * static_cast(*y_ext_[1]) + std::sin(delta) * static_cast(*y_ext_[2]))); + return 0; } - this->setJacValues(rtemp, ctemp, valtemp); - // r = 15 - if (!refframe_) + template + int DistributedGenerator::evaluateAdjointIntegrand() { - rtemp.assign(3, 15); - ctemp = {0, 3, 15}; - valtemp = {-1.0, -mp_, -alpha_}; - this->setJacValues(rtemp, ctemp, valtemp); + return 0; } - return 0; - } - - template - int DistributedGenerator::allocate() - { - CircuitComponent::allocate(); - - this->setExternalConnectionNodes(0, node_ref_->getNodeConnection(0)); - this->setExternalConnectionNodes(1, node_bus_->getNodeConnection(0)); - this->setExternalConnectionNodes(2, node_bus_->getNodeConnection(1)); - - return 0; - } - - template - int DistributedGenerator::evaluateIntegrand() - { - return 0; - } - - template - int DistributedGenerator::initializeAdjoint() - { - return 0; - } - - template - int DistributedGenerator::evaluateAdjointResidual() - { - return 0; - } - - template - int DistributedGenerator::evaluateAdjointIntegrand() - { - return 0; - } - - template - CircuitComponent* DistributedGenerator::clone() const - { - return new DistributedGenerator(*this); - } + template + Component* DistributedGenerator::clone() const + { + return new DistributedGenerator(*this); + } - // Available template instantiations - template class DistributedGenerator; - template class DistributedGenerator; - template class DistributedGenerator; - template class DistributedGenerator; + // Available template instantiations + template class DistributedGenerator; + template class DistributedGenerator; + template class DistributedGenerator; + template class DistributedGenerator; + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp index 11d6f9e97..eb667d67d 100644 --- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp +++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp @@ -1,114 +1,114 @@ - #pragma once -#include +#include #include namespace GridKit { - template - class BaseBus; - - template - struct DistributedGeneratorParameters + namespace PowerElectronics { - RealT wb_; - RealT wc_; - RealT mp_; - RealT Vn_; - RealT nq_; - RealT F_; - RealT Kiv_; - RealT Kpv_; - RealT Kic_; - RealT Kpc_; - RealT Cf_; - RealT rLf_; - RealT Lf_; - RealT rLc_; - RealT Lc_; - }; -} // namespace GridKit + /*! + * @brief Declaration of a DistributedGenerator parameter struct. + * + */ + template + struct DistributedGeneratorParameters + { + RealT wb_; + RealT wc_; + RealT mp_; + RealT Vn_; + RealT nq_; + RealT F_; + RealT Kiv_; + RealT Kpv_; + RealT Kic_; + RealT Kpc_; + RealT Cf_; + RealT rLf_; + RealT Lf_; + RealT rLc_; + RealT Lc_; + }; -namespace GridKit -{ - /*! - * @brief Declaration of a DistributedGenerator class. - * - */ - template - class DistributedGenerator : public CircuitComponent - { - using RealT = typename CircuitComponent::RealT; - using NodeT = typename PowerElectronics::NodeBase; + /*! + * @brief Declaration of a DistributedGenerator class. + * + */ + template + class DistributedGenerator : public Component + { + using RealT = typename Component::RealT; + using NodeT = typename PowerElectronics::NodeBase; - using CircuitComponent::size_; - using CircuitComponent::nnz_; - using CircuitComponent::time_; - using CircuitComponent::alpha_; - using CircuitComponent::y_ext_; - using CircuitComponent::y_int_; - using CircuitComponent::yp_ext_; - using CircuitComponent::yp_int_; - using CircuitComponent::abs_tol_; - using CircuitComponent::tag_; - using CircuitComponent::f_ext_; - using CircuitComponent::f_int_; - using CircuitComponent::g_; - using CircuitComponent::yB_; - using CircuitComponent::ypB_; - using CircuitComponent::fB_; - using CircuitComponent::gB_; - using CircuitComponent::param_; - using CircuitComponent::idc_; + using Component::size_; + using Component::nnz_; + using Component::time_; + using Component::alpha_; + using Component::y_ext_; + using Component::y_int_; + using Component::yp_ext_; + using Component::yp_int_; + using Component::abs_tol_; + using Component::tag_; + using Component::f_ext_; + using Component::f_int_; + using Component::g_; + using Component::yB_; + using Component::ypB_; + using Component::fB_; + using Component::gB_; + using Component::param_; + using Component::idc_; - using CircuitComponent::extern_indices_; - using CircuitComponent::n_extern_; - using CircuitComponent::n_intern_; + using Component::extern_indices_; + using Component::n_extern_; + using Component::n_intern_; - public: - DistributedGenerator(IdxT id, - DistributedGeneratorParameters parm, - bool reference_frame, - NodeT* node_ref, - NodeT* node_bus); - virtual ~DistributedGenerator(); + public: + DistributedGenerator(IdxT id, + DistributedGeneratorParameters parm, + bool reference_frame, + NodeT* node_ref, + NodeT* node_bus); + virtual ~DistributedGenerator(); - int initialize(); - int allocate() final; - int tagDifferentiable(); - int setAbsoluteTolerance(RealT); - int evaluateInternalResidual() final; - int evaluateExternalResidual() final; - int evaluateJacobian(); - int evaluateIntegrand(); - int initializeAdjoint(); - int evaluateAdjointResidual(); - // int evaluateAdjointJacobian(); - int evaluateAdjointIntegrand(); + int initialize(); + int allocate() final; + int tagDifferentiable(); + int setAbsoluteTolerance(RealT); + int evaluateInternalResidual() final; + int evaluateExternalResidual() final; + int evaluateJacobian(); + int evaluateIntegrand(); + int initializeAdjoint(); + int evaluateAdjointResidual(); + // int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); - CircuitComponent* clone() const; + Component* clone() const; - private: - RealT wb_; - RealT wc_; - RealT mp_; - RealT Vn_; - RealT nq_; - RealT F_; - RealT Kiv_; - RealT Kpv_; - RealT Kic_; - RealT Kpc_; - RealT Cf_; - RealT rLf_; - RealT Lf_; - RealT rLc_; - RealT Lc_; - bool refframe_; + private: + RealT wb_; + RealT wc_; + RealT mp_; + RealT Vn_; + RealT nq_; + RealT F_; + RealT Kiv_; + RealT Kpv_; + RealT Kic_; + RealT Kpc_; + RealT Cf_; + RealT rLf_; + RealT Lf_; + RealT rLc_; + RealT Lc_; + bool refframe_; - NodeT* node_ref_; - NodeT* node_bus_; - }; + NodeT* node_ref_; + NodeT* node_bus_; + }; + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/ExternalConnection.hpp b/GridKit/Model/PowerElectronics/ExternalConnection.hpp index 5edf9a7e8..f0a3322cd 100644 --- a/GridKit/Model/PowerElectronics/ExternalConnection.hpp +++ b/GridKit/Model/PowerElectronics/ExternalConnection.hpp @@ -2,22 +2,25 @@ namespace GridKit { - /** - * @brief The connection of a component's external variable to its system. - * Allows the component to access data about that external variable and allows - * the system to construct Jacobian information about that variable. - * @see CircuitComponent::setExternalConnectionNodes() - */ - template - struct ExternalConnection + namespace PowerElectronics { - /// A pointer to the state value of the variable. \see CircuitComponent::y_ext_ - const ScalarT* y_; - /// A pointer to the derivative value of the variable. \see CircuitComponent::yp_ext_ - const ScalarT* yp_; - /// A pointer to the residual buffer of the variable. \see CircuitComponent::f_ext_ - ScalarT* f_; - /// The corresponding system variable index. \see CircuitComponent::connection_nodes_ - IdxT idx_; - }; + /** + * @brief The connection of a component's external variable to its system. + * Allows the component to access data about that external variable and allows + * the system to construct Jacobian information about that variable. + * @see Component::setExternalConnectionNodes() + */ + template + struct ExternalConnection + { + /// A pointer to the state value of the variable. \see Component::y_ext_ + const ScalarT* y_; + /// A pointer to the derivative value of the variable. \see Component::yp_ext_ + const ScalarT* yp_; + /// A pointer to the residual buffer of the variable. \see Component::f_ext_ + ScalarT* f_; + /// The corresponding system variable index. \see Component::connection_nodes_ + IdxT idx_; + }; + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.cpp b/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.cpp index 0c7c87761..41c5c7b40 100644 --- a/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.cpp +++ b/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.cpp @@ -1,5 +1,4 @@ - #include "InductionMotor.hpp" #include @@ -8,158 +7,159 @@ namespace GridKit { - - /*! - * @brief Constructor for a constant InductionMotor model - * - * Calls default ModelEvaluatorImpl constructor. - * @todo create a test case utilizing the component. - * @todo create a unit test to check correctness of component - * - * @tparam ScalarT - data type for scalar variables in the model - * @tparam IdxT - integer index type for the model - * - * @param[in] id - unique identifier for the component - * @param[in] Lls - stator leakage inductance - */ - - template - InductionMotor::InductionMotor(IdxT id, RealT Lls, RealT Rs, RealT Llr, RealT Rr, RealT Lms, RealT RJ, RealT P) - : Lls_(Lls), - Rs_(Rs), - Llr_(Llr), - Rr_(Rr), - Lms_(Lms), - RJ_(RJ), - P_(P) - { - size_ = 10; - n_intern_ = 5; - n_extern_ = 5; - extern_indices_ = {0, 1, 2, 3, 4}; - idc_ = id; - } - - template - InductionMotor::~InductionMotor() - { - } - - /** - * Initialization of the grid model - */ - template - int InductionMotor::initialize() - { - return 0; - } - - /* - * \brief Identify differential variables - */ - template - int InductionMotor::tagDifferentiable() + namespace PowerElectronics { - // All variables are differentials - std::fill(tag_.begin(), tag_.end(), true); - return 0; - } - - /** - * @brief Compute the absolute tolerance for each variable in the model - * - * @param rel_tol The relative tolerance which can be used to pick the - * absolute tolerance. - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return int 0 if successful, non-zero otherwise. - * - * This represents a "noise" level close to zero for which pure relative - * error cannot be used. - */ - template - int InductionMotor::setAbsoluteTolerance(RealT rel_tol) - { - abs_tol_.setToConst(static_cast(rel_tol)); - return 0; - } - - /** - * @brief Contributes to the resisdual - * - */ - template - int InductionMotor::evaluateInternalResidual() - { - f_int_[0] = (1.0 / 3.0) * (2.0 * *y_ext_[0] - *y_ext_[1] - *y_ext_[2]) - Rs_ * y_int_[0] - (Lls_ + Lms_) * yp_int_[0] - Lms_ * yp_int_[1]; - f_int_[1] = (1.0 / std::sqrt(3.0)) * (-*y_ext_[1] + *y_ext_[2]) - Rs_ * y_int_[1] - (Lls_ + Lms_) * yp_int_[1] - Lms_ * yp_int_[0]; - f_int_[2] = (*y_ext_[0] + *y_ext_[1] + *y_ext_[2]) / 3.0 - Rs_ * y_int_[2] - Lls_ * yp_int_[7]; - f_int_[3] = Rr_ * y_int_[3] + (Llr_ + Lms_) * yp_int_[3] + Lms_ * yp_int_[0] - (P_ / 2.0) * *y_ext_[3] * ((Llr_ + Lms_) * y_int_[4] + Lms_ * y_int_[1]); - f_int_[4] = Rr_ * y_int_[4] + (Llr_ + Lms_) * yp_int_[4] + Lms_ * yp_int_[1] + (P_ / 2.0) * *y_ext_[3] * ((Llr_ + Lms_) * y_int_[3] + Lms_ * y_int_[0]); - return 0; - } - - template - int InductionMotor::evaluateExternalResidual() - { - *f_ext_[0] += y_int_[0] + y_int_[2]; - *f_ext_[1] += (-1.0 / 2.0) * y_int_[0] - (std::sqrt(3.0) / 2.0) * y_int_[1] + y_int_[2]; - *f_ext_[2] += (-1.0 / 2.0) * y_int_[0] + (std::sqrt(3.0) / 2.0) * y_int_[1] + y_int_[2]; - *f_ext_[3] += RJ_ * *yp_ext_[3] - (3.0 / 4.0) * P_ * Lms_ * (y_int_[0] * y_int_[4] - y_int_[1] * y_int_[3]); - *f_ext_[4] += *yp_ext_[4] - *y_ext_[3]; - return 0; - } - - /** - * @brief Compute component Jacobian - * - * @todo need to implement - * - * @tparam ScalarT - * @tparam IdxT - * @return int - */ - template - int InductionMotor::evaluateJacobian() - { - - return 0; - } - - template - int InductionMotor::evaluateIntegrand() - { - return 0; - } - - template - int InductionMotor::initializeAdjoint() - { - return 0; - } - - template - int InductionMotor::evaluateAdjointResidual() - { - return 0; - } - - template - int InductionMotor::evaluateAdjointIntegrand() - { - return 0; - } - - template - CircuitComponent* InductionMotor::clone() const - { - return new InductionMotor(*this); - } - - // Available template instantiations - template class InductionMotor; - template class InductionMotor; - template class InductionMotor; - template class InductionMotor; - + /*! + * @brief Constructor for a constant InductionMotor model + * + * Calls default ModelEvaluatorImpl constructor. + * @todo create a test case utilizing the component. + * @todo create a unit test to check correctness of component + * + * @tparam ScalarT - data type for scalar variables in the model + * @tparam IdxT - integer index type for the model + * + * @param[in] id - unique identifier for the component + * @param[in] Lls - stator leakage inductance + */ + template + InductionMotor::InductionMotor(IdxT id, RealT Lls, RealT Rs, RealT Llr, RealT Rr, RealT Lms, RealT RJ, RealT P) + : Lls_(Lls), + Rs_(Rs), + Llr_(Llr), + Rr_(Rr), + Lms_(Lms), + RJ_(RJ), + P_(P) + { + size_ = 10; + n_intern_ = 5; + n_extern_ = 5; + extern_indices_ = {0, 1, 2, 3, 4}; + idc_ = id; + } + + template + InductionMotor::~InductionMotor() + { + } + + /** + * Initialization of the grid model + */ + template + int InductionMotor::initialize() + { + return 0; + } + + /* + * \brief Identify differential variables + */ + template + int InductionMotor::tagDifferentiable() + { + // All variables are differentials + std::fill(tag_.begin(), tag_.end(), true); + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int InductionMotor::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Contributes to the resisdual + * + */ + template + int InductionMotor::evaluateInternalResidual() + { + f_int_[0] = (1.0 / 3.0) * (2.0 * *y_ext_[0] - *y_ext_[1] - *y_ext_[2]) - Rs_ * y_int_[0] - (Lls_ + Lms_) * yp_int_[0] - Lms_ * yp_int_[1]; + f_int_[1] = (1.0 / std::sqrt(3.0)) * (-*y_ext_[1] + *y_ext_[2]) - Rs_ * y_int_[1] - (Lls_ + Lms_) * yp_int_[1] - Lms_ * yp_int_[0]; + f_int_[2] = (*y_ext_[0] + *y_ext_[1] + *y_ext_[2]) / 3.0 - Rs_ * y_int_[2] - Lls_ * yp_int_[7]; + f_int_[3] = Rr_ * y_int_[3] + (Llr_ + Lms_) * yp_int_[3] + Lms_ * yp_int_[0] - (P_ / 2.0) * *y_ext_[3] * ((Llr_ + Lms_) * y_int_[4] + Lms_ * y_int_[1]); + f_int_[4] = Rr_ * y_int_[4] + (Llr_ + Lms_) * yp_int_[4] + Lms_ * yp_int_[1] + (P_ / 2.0) * *y_ext_[3] * ((Llr_ + Lms_) * y_int_[3] + Lms_ * y_int_[0]); + return 0; + } + + template + int InductionMotor::evaluateExternalResidual() + { + *f_ext_[0] += y_int_[0] + y_int_[2]; + *f_ext_[1] += (-1.0 / 2.0) * y_int_[0] - (std::sqrt(3.0) / 2.0) * y_int_[1] + y_int_[2]; + *f_ext_[2] += (-1.0 / 2.0) * y_int_[0] + (std::sqrt(3.0) / 2.0) * y_int_[1] + y_int_[2]; + *f_ext_[3] += RJ_ * *yp_ext_[3] - (3.0 / 4.0) * P_ * Lms_ * (y_int_[0] * y_int_[4] - y_int_[1] * y_int_[3]); + *f_ext_[4] += *yp_ext_[4] - *y_ext_[3]; + return 0; + } + + /** + * @brief Compute component Jacobian + * + * @todo need to implement + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ + template + int InductionMotor::evaluateJacobian() + { + + return 0; + } + + template + int InductionMotor::evaluateIntegrand() + { + return 0; + } + + template + int InductionMotor::initializeAdjoint() + { + return 0; + } + + template + int InductionMotor::evaluateAdjointResidual() + { + return 0; + } + + template + int InductionMotor::evaluateAdjointIntegrand() + { + return 0; + } + + template + Component* InductionMotor::clone() const + { + return new InductionMotor(*this); + } + + // Available template instantiations + template class InductionMotor; + template class InductionMotor; + template class InductionMotor; + template class InductionMotor; + + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.hpp b/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.hpp index 13ed162ec..a79e1ffa8 100644 --- a/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.hpp +++ b/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.hpp @@ -1,76 +1,72 @@ - #pragma once -#include - -namespace GridKit -{ - template - class BaseBus; -} +#include namespace GridKit { - /*! - * @brief Declaration of a InductionMotor class. - * - */ - template - class InductionMotor : public CircuitComponent + namespace PowerElectronics { - using RealT = typename CircuitComponent::RealT; + /*! + * @brief Declaration of a InductionMotor class. + * + */ + template + class InductionMotor : public Component + { + using RealT = typename Component::RealT; - using CircuitComponent::size_; - using CircuitComponent::nnz_; - using CircuitComponent::time_; - using CircuitComponent::alpha_; - using CircuitComponent::y_ext_; - using CircuitComponent::y_int_; - using CircuitComponent::yp_ext_; - using CircuitComponent::yp_int_; - using CircuitComponent::abs_tol_; - using CircuitComponent::tag_; - using CircuitComponent::f_ext_; - using CircuitComponent::f_int_; - using CircuitComponent::g_; - using CircuitComponent::yB_; - using CircuitComponent::ypB_; - using CircuitComponent::fB_; - using CircuitComponent::gB_; - using CircuitComponent::param_; - using CircuitComponent::idc_; + using Component::size_; + using Component::nnz_; + using Component::time_; + using Component::alpha_; + using Component::y_ext_; + using Component::y_int_; + using Component::yp_ext_; + using Component::yp_int_; + using Component::abs_tol_; + using Component::tag_; + using Component::f_ext_; + using Component::f_int_; + using Component::g_; + using Component::yB_; + using Component::ypB_; + using Component::fB_; + using Component::gB_; + using Component::param_; + using Component::idc_; - using CircuitComponent::extern_indices_; - using CircuitComponent::n_extern_; - using CircuitComponent::n_intern_; + using Component::extern_indices_; + using Component::n_extern_; + using Component::n_intern_; - public: - InductionMotor(IdxT id, RealT Lls, RealT Rs, RealT Llr, RealT Rr, RealT Lms, RealT RJ, RealT P); - virtual ~InductionMotor(); + public: + InductionMotor(IdxT id, RealT Lls, RealT Rs, RealT Llr, RealT Rr, RealT Lms, RealT RJ, RealT P); + virtual ~InductionMotor(); - int initialize(); - int tagDifferentiable(); - int setAbsoluteTolerance(RealT); - int evaluateInternalResidual() final; - int evaluateExternalResidual() final; - int evaluateJacobian(); - int evaluateIntegrand(); + int initialize(); + int tagDifferentiable(); + int setAbsoluteTolerance(RealT); + int evaluateInternalResidual() final; + int evaluateExternalResidual() final; + int evaluateJacobian(); + int evaluateIntegrand(); - int initializeAdjoint(); - int evaluateAdjointResidual(); - // int evaluateAdjointJacobian(); - int evaluateAdjointIntegrand(); + int initializeAdjoint(); + int evaluateAdjointResidual(); + // int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); - CircuitComponent* clone() const; + Component* clone() const; - private: - RealT Lls_; - RealT Rs_; - RealT Llr_; - RealT Rr_; - RealT Lms_; - RealT RJ_; - RealT P_; - }; + private: + RealT Lls_; + RealT Rs_; + RealT Llr_; + RealT Rr_; + RealT Lms_; + RealT RJ_; + RealT P_; + }; + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp b/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp index 42c15c6a5..2d367b608 100644 --- a/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp +++ b/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp @@ -1,5 +1,4 @@ - #include "Inductor.hpp" #include @@ -7,158 +6,159 @@ namespace GridKit { - - /*! - * @brief Constructor for a inductor - * - * Calls default ModelEvaluatorImpl constructor. - */ - - template - Inductor::Inductor(IdxT id, RealT L, NodeT* node1, NodeT* node2) - : L_(L), node1_(node1), node2_(node2) - { - assert(node1_->size() == 1); - assert(node2_->size() == 1); - size_ = 3; - n_intern_ = 1; - n_extern_ = 2; - extern_indices_ = {0, 1}; - idc_ = id; - nnz_ = 5; - } - - template - Inductor::~Inductor() - { - } - - /** - * Initialization of the grid model - */ - template - int Inductor::initialize() - { - return 0; - } - - /* - * \brief Identify differential variables - */ - template - int Inductor::tagDifferentiable() - { - // All variables are differentials - std::fill(tag_.begin(), tag_.end(), true); - return 0; - } - - /** - * @brief Compute the absolute tolerance for each variable in the model - * - * @param rel_tol The relative tolerance which can be used to pick the - * absolute tolerance. - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return int 0 if successful, non-zero otherwise. - * - * This represents a "noise" level close to zero for which pure relative - * error cannot be used. - */ - template - int Inductor::setAbsoluteTolerance(RealT rel_tol) - { - abs_tol_.setToConst(static_cast(rel_tol)); - return 0; - } - - /** - * @brief Compute the resisdual of the component - * - */ - template - int Inductor::evaluateInternalResidual() - { - f_int_[0] = -L_ * yp_int_[0] + *y_ext_[1] - *y_ext_[0]; - return 0; - } - - template - int Inductor::evaluateExternalResidual() - { - // input - *f_ext_[0] += -y_int_[0]; - // output - *f_ext_[1] += y_int_[0]; - return 0; - } - - /** - * @brief Evaluate the jacobian of the component - * - * @tparam ScalarT - * @tparam IdxT - * @return int - */ - template - int Inductor::evaluateJacobian() - { - this->zeroJacMatrix(); - - // Create dF/dy - std::vector rcord{0, 1, 2, 2, 2}; - std::vector ccord{2, 2, 0, 1, 2}; - std::vector vals{-1.0, 1.0, -1.0, 1.0, -L_ * alpha_}; - this->setJacValues(rcord, ccord, vals); - - return 0; - } - - template - int Inductor::allocate() - { - CircuitComponent::allocate(); - - this->setExternalConnectionNodes(0, node1_->getNodeConnection(0)); - this->setExternalConnectionNodes(1, node2_->getNodeConnection(0)); - - return 0; - } - - template - int Inductor::evaluateIntegrand() - { - return 0; - } - - template - int Inductor::initializeAdjoint() + namespace PowerElectronics { - return 0; - } - - template - int Inductor::evaluateAdjointResidual() - { - return 0; - } - - template - int Inductor::evaluateAdjointIntegrand() - { - return 0; - } - - template - CircuitComponent* Inductor::clone() const - { - return new Inductor(*this); - } - - // Available template instantiations - template class Inductor; - template class Inductor; - template class Inductor; - template class Inductor; - + /*! + * @brief Constructor for a inductor + * + * Calls default ModelEvaluatorImpl constructor. + */ + template + Inductor::Inductor(IdxT id, RealT L, NodeT* node1, NodeT* node2) + : L_(L), node1_(node1), node2_(node2) + { + assert(node1_->size() == 1); + assert(node2_->size() == 1); + size_ = 3; + n_intern_ = 1; + n_extern_ = 2; + extern_indices_ = {0, 1}; + idc_ = id; + nnz_ = 5; + } + + template + Inductor::~Inductor() + { + } + + /** + * Initialization of the grid model + */ + template + int Inductor::initialize() + { + return 0; + } + + /* + * \brief Identify differential variables + */ + template + int Inductor::tagDifferentiable() + { + // All variables are differentials + std::fill(tag_.begin(), tag_.end(), true); + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Inductor::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Compute the resisdual of the component + * + */ + template + int Inductor::evaluateInternalResidual() + { + f_int_[0] = -L_ * yp_int_[0] + *y_ext_[1] - *y_ext_[0]; + return 0; + } + + template + int Inductor::evaluateExternalResidual() + { + // input + *f_ext_[0] += -y_int_[0]; + // output + *f_ext_[1] += y_int_[0]; + return 0; + } + + /** + * @brief Evaluate the jacobian of the component + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ + template + int Inductor::evaluateJacobian() + { + this->zeroJacMatrix(); + + // Create dF/dy + std::vector rcord{0, 1, 2, 2, 2}; + std::vector ccord{2, 2, 0, 1, 2}; + std::vector vals{-1.0, 1.0, -1.0, 1.0, -L_ * alpha_}; + this->setJacValues(rcord, ccord, vals); + + return 0; + } + + template + int Inductor::allocate() + { + Component::allocate(); + + this->setExternalConnectionNodes(0, node1_->getNodeConnection(0)); + this->setExternalConnectionNodes(1, node2_->getNodeConnection(0)); + + return 0; + } + + template + int Inductor::evaluateIntegrand() + { + return 0; + } + + template + int Inductor::initializeAdjoint() + { + return 0; + } + + template + int Inductor::evaluateAdjointResidual() + { + return 0; + } + + template + int Inductor::evaluateAdjointIntegrand() + { + return 0; + } + + template + Component* Inductor::clone() const + { + return new Inductor(*this); + } + + // Available template instantiations + template class Inductor; + template class Inductor; + template class Inductor; + template class Inductor; + + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp b/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp index 7440037c9..1b4a5f07f 100644 --- a/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp +++ b/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp @@ -1,75 +1,71 @@ - #pragma once -#include +#include #include namespace GridKit { - template - class BaseBus; -} - -namespace GridKit -{ - /*! - * @brief Declaration of a Inductor class. - * - */ - template - class Inductor : public CircuitComponent + namespace PowerElectronics { - using RealT = typename CircuitComponent::RealT; - using NodeT = typename PowerElectronics::NodeBase; + /*! + * @brief Declaration of a Inductor class. + * + */ + template + class Inductor : public Component + { + using RealT = typename Component::RealT; + using NodeT = typename PowerElectronics::NodeBase; - using CircuitComponent::size_; - using CircuitComponent::nnz_; - using CircuitComponent::time_; - using CircuitComponent::alpha_; - using CircuitComponent::y_ext_; - using CircuitComponent::y_int_; - using CircuitComponent::yp_ext_; - using CircuitComponent::yp_int_; - using CircuitComponent::tag_; - using CircuitComponent::abs_tol_; - using CircuitComponent::f_ext_; - using CircuitComponent::f_int_; - using CircuitComponent::g_; - using CircuitComponent::yB_; - using CircuitComponent::ypB_; - using CircuitComponent::fB_; - using CircuitComponent::gB_; - using CircuitComponent::param_; - using CircuitComponent::idc_; + using Component::size_; + using Component::nnz_; + using Component::time_; + using Component::alpha_; + using Component::y_ext_; + using Component::y_int_; + using Component::yp_ext_; + using Component::yp_int_; + using Component::tag_; + using Component::abs_tol_; + using Component::f_ext_; + using Component::f_int_; + using Component::g_; + using Component::yB_; + using Component::ypB_; + using Component::fB_; + using Component::gB_; + using Component::param_; + using Component::idc_; - using CircuitComponent::extern_indices_; - using CircuitComponent::n_extern_; - using CircuitComponent::n_intern_; + using Component::extern_indices_; + using Component::n_extern_; + using Component::n_intern_; - public: - Inductor(IdxT id, RealT L, NodeT* node1, NodeT* node2); - virtual ~Inductor(); + public: + Inductor(IdxT id, RealT L, NodeT* node1, NodeT* node2); + virtual ~Inductor(); - int initialize(); - int allocate() final; - int tagDifferentiable(); - int setAbsoluteTolerance(RealT); - int evaluateInternalResidual() final; - int evaluateExternalResidual() final; - int evaluateJacobian(); - int evaluateIntegrand(); + int initialize(); + int allocate() final; + int tagDifferentiable(); + int setAbsoluteTolerance(RealT); + int evaluateInternalResidual() final; + int evaluateExternalResidual() final; + int evaluateJacobian(); + int evaluateIntegrand(); - int initializeAdjoint(); - int evaluateAdjointResidual(); - // int evaluateAdjointJacobian(); - int evaluateAdjointIntegrand(); + int initializeAdjoint(); + int evaluateAdjointResidual(); + // int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); - CircuitComponent* clone() const; + Component* clone() const; - private: - RealT L_; - NodeT* node1_; - NodeT* node2_; - }; + private: + RealT L_; + NodeT* node1_; + NodeT* node2_; + }; + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp b/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp index b3b909ebb..9e33b79c2 100644 --- a/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp +++ b/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp @@ -1,5 +1,4 @@ - #include "LinearTransformer.hpp" #include @@ -8,143 +7,144 @@ namespace GridKit { - - /*! - * @brief Constructor for a LinearTransformer model - * - * Calls default ModelEvaluatorImpl constructor. - * @todo Not tested in any model yet. Should be - * @todo Has not been tested for correctness - * - * @tparam ScalarT - floating point type for the model - * @tparam IdxT - integer index type for the model - * - * @param[in] id - unique identifier for the component - * @param[in] L0 - inductance 0 - * @param[in] L1 - inductance 1 - * @param[in] R0 - resistance 0 - * @param[in] R1 - resistance 1 - * @param[in] M - mutual inductance - */ - - template - LinearTransformer::LinearTransformer(IdxT id, RealT L0, RealT L1, RealT R0, RealT R1, RealT M) - : L0_(L0), - L1_(L1), - R0_(R0), - R1_(R1), - M_(M) - { - size_ = 4; - n_intern_ = 2; - n_extern_ = 2; - extern_indices_ = {0, 1}; - idc_ = id; - } - - template - LinearTransformer::~LinearTransformer() - { - } - - /** - * Initialization of the grid model - */ - template - int LinearTransformer::initialize() - { - return 0; - } - - /* - * \brief Identify differential variables - */ - template - int LinearTransformer::tagDifferentiable() - { - // All variables are differentials - std::fill(tag_.begin(), tag_.end(), true); - return 0; - } - - /** - * @brief Compute the absolute tolerance for each variable in the model - * - * @param rel_tol The relative tolerance which can be used to pick the - * absolute tolerance. - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return int 0 if successful, non-zero otherwise. - * - * This represents a "noise" level close to zero for which pure relative - * error cannot be used. - */ - template - int LinearTransformer::setAbsoluteTolerance(RealT rel_tol) - { - abs_tol_.setToConst(static_cast(rel_tol)); - return 0; - } - - /** - * @brief Computes the component resisdual - */ - template - int LinearTransformer::evaluateInternalResidual() - { - f_int_[0] = *y_ext_[0] - R0_ * y_int_[0] - L0_ * yp_int_[0] - M_ * yp_int_[1]; - f_int_[1] = *y_ext_[1] - R1_ * y_int_[1] - M_ * yp_int_[0] - L1_ * yp_int_[1]; - return 0; - } - - template - int LinearTransformer::evaluateExternalResidual() - { - *f_ext_[0] += y_int_[0]; - *f_ext_[1] += y_int_[1]; - return 0; - } - - template - int LinearTransformer::evaluateJacobian() + namespace PowerElectronics { - return 0; - } - - template - int LinearTransformer::evaluateIntegrand() - { - return 0; - } - - template - int LinearTransformer::initializeAdjoint() - { - return 0; - } - - template - int LinearTransformer::evaluateAdjointResidual() - { - return 0; - } - - template - int LinearTransformer::evaluateAdjointIntegrand() - { - return 0; - } - - template - CircuitComponent* LinearTransformer::clone() const - { - return new LinearTransformer(*this); - } - - // Available template instantiations - template class LinearTransformer; - template class LinearTransformer; - template class LinearTransformer; - template class LinearTransformer; - + /*! + * @brief Constructor for a LinearTransformer model + * + * Calls default ModelEvaluatorImpl constructor. + * @todo Not tested in any model yet. Should be + * @todo Has not been tested for correctness + * + * @tparam ScalarT - floating point type for the model + * @tparam IdxT - integer index type for the model + * + * @param[in] id - unique identifier for the component + * @param[in] L0 - inductance 0 + * @param[in] L1 - inductance 1 + * @param[in] R0 - resistance 0 + * @param[in] R1 - resistance 1 + * @param[in] M - mutual inductance + */ + template + LinearTransformer::LinearTransformer(IdxT id, RealT L0, RealT L1, RealT R0, RealT R1, RealT M) + : L0_(L0), + L1_(L1), + R0_(R0), + R1_(R1), + M_(M) + { + size_ = 4; + n_intern_ = 2; + n_extern_ = 2; + extern_indices_ = {0, 1}; + idc_ = id; + } + + template + LinearTransformer::~LinearTransformer() + { + } + + /** + * Initialization of the grid model + */ + template + int LinearTransformer::initialize() + { + return 0; + } + + /* + * \brief Identify differential variables + */ + template + int LinearTransformer::tagDifferentiable() + { + // All variables are differentials + std::fill(tag_.begin(), tag_.end(), true); + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int LinearTransformer::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Computes the component resisdual + */ + template + int LinearTransformer::evaluateInternalResidual() + { + f_int_[0] = *y_ext_[0] - R0_ * y_int_[0] - L0_ * yp_int_[0] - M_ * yp_int_[1]; + f_int_[1] = *y_ext_[1] - R1_ * y_int_[1] - M_ * yp_int_[0] - L1_ * yp_int_[1]; + return 0; + } + + template + int LinearTransformer::evaluateExternalResidual() + { + *f_ext_[0] += y_int_[0]; + *f_ext_[1] += y_int_[1]; + return 0; + } + + template + int LinearTransformer::evaluateJacobian() + { + return 0; + } + + template + int LinearTransformer::evaluateIntegrand() + { + return 0; + } + + template + int LinearTransformer::initializeAdjoint() + { + return 0; + } + + template + int LinearTransformer::evaluateAdjointResidual() + { + return 0; + } + + template + int LinearTransformer::evaluateAdjointIntegrand() + { + return 0; + } + + template + Component* LinearTransformer::clone() const + { + return new LinearTransformer(*this); + } + + // Available template instantiations + template class LinearTransformer; + template class LinearTransformer; + template class LinearTransformer; + template class LinearTransformer; + + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp b/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp index 0785dcb32..11ed78dd4 100644 --- a/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp +++ b/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp @@ -1,74 +1,70 @@ - #pragma once -#include - -namespace GridKit -{ - template - class BaseBus; -} +#include namespace GridKit { - /*! - * @brief Declaration of a LinearTransformer class. - * - */ - template - class LinearTransformer : public CircuitComponent + namespace PowerElectronics { - using RealT = typename CircuitComponent::RealT; + /*! + * @brief Declaration of a LinearTransformer class. + * + */ + template + class LinearTransformer : public Component + { + using RealT = typename Component::RealT; - using CircuitComponent::size_; - using CircuitComponent::nnz_; - using CircuitComponent::time_; - using CircuitComponent::alpha_; - using CircuitComponent::y_ext_; - using CircuitComponent::y_int_; - using CircuitComponent::yp_ext_; - using CircuitComponent::yp_int_; - using CircuitComponent::abs_tol_; - using CircuitComponent::tag_; - using CircuitComponent::f_ext_; - using CircuitComponent::f_int_; - using CircuitComponent::g_; - using CircuitComponent::yB_; - using CircuitComponent::ypB_; - using CircuitComponent::fB_; - using CircuitComponent::gB_; - using CircuitComponent::param_; - using CircuitComponent::idc_; + using Component::size_; + using Component::nnz_; + using Component::time_; + using Component::alpha_; + using Component::y_ext_; + using Component::y_int_; + using Component::yp_ext_; + using Component::yp_int_; + using Component::abs_tol_; + using Component::tag_; + using Component::f_ext_; + using Component::f_int_; + using Component::g_; + using Component::yB_; + using Component::ypB_; + using Component::fB_; + using Component::gB_; + using Component::param_; + using Component::idc_; - using CircuitComponent::extern_indices_; - using CircuitComponent::n_extern_; - using CircuitComponent::n_intern_; + using Component::extern_indices_; + using Component::n_extern_; + using Component::n_intern_; - public: - LinearTransformer(IdxT id, RealT L0, RealT L1, RealT R0, RealT R1, RealT M); - virtual ~LinearTransformer(); + public: + LinearTransformer(IdxT id, RealT L0, RealT L1, RealT R0, RealT R1, RealT M); + virtual ~LinearTransformer(); - int initialize(); - int tagDifferentiable(); - int setAbsoluteTolerance(RealT); - int evaluateInternalResidual() final; - int evaluateExternalResidual() final; - int evaluateJacobian(); - int evaluateIntegrand(); + int initialize(); + int tagDifferentiable(); + int setAbsoluteTolerance(RealT); + int evaluateInternalResidual() final; + int evaluateExternalResidual() final; + int evaluateJacobian(); + int evaluateIntegrand(); - int initializeAdjoint(); - int evaluateAdjointResidual(); - // int evaluateAdjointJacobian(); - int evaluateAdjointIntegrand(); + int initializeAdjoint(); + int evaluateAdjointResidual(); + // int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); - CircuitComponent* clone() const; + Component* clone() const; - private: - RealT L0_; - RealT L1_; - RealT R0_; - RealT R1_; - RealT M_; - }; + private: + RealT L0_; + RealT L1_; + RealT R0_; + RealT R1_; + RealT M_; + }; + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp index a918eaf9d..fec34aa9c 100644 --- a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp @@ -7,161 +7,163 @@ namespace GridKit { - - /*! - * @brief Constructor for a constant MicrogridBusDQ model - * - * Calls default ModelEvaluatorImpl constructor. - * - * In DQ space - * Each microgrid line has a virtual resistance RN - * Model is from paper: "Modeling, Analysis and Testing of Autonomous Operation - * of an Inverter-Based Microgrid", Nagaraju Pogaku, Milan Prodanovic, and - * Timothy C. Green, Section E - */ - template - MicrogridBusDQ::MicrogridBusDQ(IdxT id, RealT RN, NodeT* node1) - : RN_(RN), node1_(node1) - { - assert(node1_->size() == 2); - // externals [vbus_d, vbus_q] - size_ = 2; - n_intern_ = 0; - n_extern_ = 2; - extern_indices_ = {0, 1}; - idc_ = id; - nnz_ = 2; - } - - template - MicrogridBusDQ::~MicrogridBusDQ() - { - } - - /** - * Initialization of the grid model - */ - template - int MicrogridBusDQ::initialize() - { - return 0; - } - - /// There are no internal variables in this component, so \ref tag_ can be set arbitrarily. - template - int MicrogridBusDQ::tagDifferentiable() - { - return 0; - } - - /** - * @brief Compute the absolute tolerance for each variable in the model - * - * @param rel_tol The relative tolerance which can be used to pick the - * absolute tolerance. - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return int 0 if successful, non-zero otherwise. - * - * This represents a "noise" level close to zero for which pure relative - * error cannot be used. - */ - template - int MicrogridBusDQ::setAbsoluteTolerance(RealT) - { - return 0; - } - - template - int MicrogridBusDQ::evaluateInternalResidual() - { - return 0; - } - - /** - * @brief Evaluate residual - * This model has "Virtual resistors". The voltage of the bus divided by its virtual resistance. - * The components are external to allow for outside components to add inductances to the terms. - * - * refernce to equations in class header - * - */ - template - int MicrogridBusDQ::evaluateExternalResidual() - { - // bus voltage - *f_ext_[0] += -*y_ext_[0] / RN_; - *f_ext_[1] += -*y_ext_[1] / RN_; - - return 0; - } - - /** - * @brief Generate Jacobian - * - * @tparam ScalarT - * @tparam IdxT - * @return int - */ - template - int MicrogridBusDQ::evaluateJacobian() - { - this->zeroJacMatrix(); - - // Create dF/dy - std::vector rtemp{0, 1}; - std::vector ctemp{0, 1}; - std::vector vals{-1.0 / RN_, -1.0 / RN_}; - this->setJacValues(rtemp, ctemp, vals); - - return 0; - } - - template - int MicrogridBusDQ::allocate() + namespace PowerElectronics { - CircuitComponent::allocate(); - - this->setExternalConnectionNodes(0, node1_->getNodeConnection(0)); - this->setExternalConnectionNodes(1, node1_->getNodeConnection(1)); - - return 0; - } - - template - int MicrogridBusDQ::evaluateIntegrand() - { - return 0; - } - - template - int MicrogridBusDQ::initializeAdjoint() - { - return 0; - } - - template - int MicrogridBusDQ::evaluateAdjointResidual() - { - return 0; - } - - template - int MicrogridBusDQ::evaluateAdjointIntegrand() - { - return 0; - } - - template - CircuitComponent* MicrogridBusDQ::clone() const - { - return new MicrogridBusDQ(*this); - } - - // Available template instantiations - template class MicrogridBusDQ; - template class MicrogridBusDQ; - template class MicrogridBusDQ; - template class MicrogridBusDQ; - + /*! + * @brief Constructor for a constant MicrogridBusDQ model + * + * Calls default ModelEvaluatorImpl constructor. + * + * In DQ space + * Each microgrid line has a virtual resistance RN + * Model is from paper: "Modeling, Analysis and Testing of Autonomous Operation + * of an Inverter-Based Microgrid", Nagaraju Pogaku, Milan Prodanovic, and + * Timothy C. Green, Section E + */ + template + MicrogridBusDQ::MicrogridBusDQ(IdxT id, RealT RN, NodeT* node1) + : RN_(RN), node1_(node1) + { + assert(node1_->size() == 2); + // externals [vbus_d, vbus_q] + size_ = 2; + n_intern_ = 0; + n_extern_ = 2; + extern_indices_ = {0, 1}; + idc_ = id; + nnz_ = 2; + } + + template + MicrogridBusDQ::~MicrogridBusDQ() + { + } + + /** + * Initialization of the grid model + */ + template + int MicrogridBusDQ::initialize() + { + return 0; + } + + /// There are no internal variables in this component, so \ref tag_ can be set arbitrarily. + template + int MicrogridBusDQ::tagDifferentiable() + { + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int MicrogridBusDQ::setAbsoluteTolerance(RealT) + { + return 0; + } + + template + int MicrogridBusDQ::evaluateInternalResidual() + { + return 0; + } + + /** + * @brief Evaluate residual + * This model has "Virtual resistors". The voltage of the bus divided by its virtual resistance. + * The components are external to allow for outside components to add inductances to the terms. + * + * refernce to equations in class header + * + */ + template + int MicrogridBusDQ::evaluateExternalResidual() + { + // bus voltage + *f_ext_[0] += -*y_ext_[0] / RN_; + *f_ext_[1] += -*y_ext_[1] / RN_; + + return 0; + } + + /** + * @brief Generate Jacobian + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ + template + int MicrogridBusDQ::evaluateJacobian() + { + this->zeroJacMatrix(); + + // Create dF/dy + std::vector rtemp{0, 1}; + std::vector ctemp{0, 1}; + std::vector vals{-1.0 / RN_, -1.0 / RN_}; + this->setJacValues(rtemp, ctemp, vals); + + return 0; + } + + template + int MicrogridBusDQ::allocate() + { + Component::allocate(); + + this->setExternalConnectionNodes(0, node1_->getNodeConnection(0)); + this->setExternalConnectionNodes(1, node1_->getNodeConnection(1)); + + return 0; + } + + template + int MicrogridBusDQ::evaluateIntegrand() + { + return 0; + } + + template + int MicrogridBusDQ::initializeAdjoint() + { + return 0; + } + + template + int MicrogridBusDQ::evaluateAdjointResidual() + { + return 0; + } + + template + int MicrogridBusDQ::evaluateAdjointIntegrand() + { + return 0; + } + + template + Component* MicrogridBusDQ::clone() const + { + return new MicrogridBusDQ(*this); + } + + // Available template instantiations + template class MicrogridBusDQ; + template class MicrogridBusDQ; + template class MicrogridBusDQ; + template class MicrogridBusDQ; + + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp index 806dd2751..8ba4cf6de 100644 --- a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp @@ -1,74 +1,70 @@ - #pragma once -#include +#include #include namespace GridKit { - template - class BaseBus; -} - -namespace GridKit -{ - /*! - * @brief Declaration of a MicrogridBusDQ class. - * - */ - template - class MicrogridBusDQ : public CircuitComponent + namespace PowerElectronics { - using RealT = typename CircuitComponent::RealT; - using NodeT = typename PowerElectronics::NodeBase; + /*! + * @brief Declaration of a MicrogridBusDQ class. + * + */ + template + class MicrogridBusDQ : public Component + { + using RealT = typename Component::RealT; + using NodeT = typename PowerElectronics::NodeBase; - using CircuitComponent::size_; - using CircuitComponent::nnz_; - using CircuitComponent::time_; - using CircuitComponent::alpha_; - using CircuitComponent::y_ext_; - using CircuitComponent::y_int_; - using CircuitComponent::yp_ext_; - using CircuitComponent::yp_int_; - using CircuitComponent::tag_; - using CircuitComponent::abs_tol_; - using CircuitComponent::f_ext_; - using CircuitComponent::f_int_; - using CircuitComponent::g_; - using CircuitComponent::yB_; - using CircuitComponent::ypB_; - using CircuitComponent::fB_; - using CircuitComponent::gB_; - using CircuitComponent::param_; - using CircuitComponent::idc_; + using Component::size_; + using Component::nnz_; + using Component::time_; + using Component::alpha_; + using Component::y_ext_; + using Component::y_int_; + using Component::yp_ext_; + using Component::yp_int_; + using Component::tag_; + using Component::abs_tol_; + using Component::f_ext_; + using Component::f_int_; + using Component::g_; + using Component::yB_; + using Component::ypB_; + using Component::fB_; + using Component::gB_; + using Component::param_; + using Component::idc_; - using CircuitComponent::extern_indices_; - using CircuitComponent::n_extern_; - using CircuitComponent::n_intern_; + using Component::extern_indices_; + using Component::n_extern_; + using Component::n_intern_; - public: - MicrogridBusDQ(IdxT id, RealT RN, NodeT* node1); - virtual ~MicrogridBusDQ(); + public: + MicrogridBusDQ(IdxT id, RealT RN, NodeT* node1); + virtual ~MicrogridBusDQ(); - int initialize(); - int allocate() final; - int tagDifferentiable(); - int setAbsoluteTolerance(RealT); - int evaluateInternalResidual() final; - int evaluateExternalResidual() final; - int evaluateJacobian(); - int evaluateIntegrand(); + int initialize(); + int allocate() final; + int tagDifferentiable(); + int setAbsoluteTolerance(RealT); + int evaluateInternalResidual() final; + int evaluateExternalResidual() final; + int evaluateJacobian(); + int evaluateIntegrand(); - int initializeAdjoint(); - int evaluateAdjointResidual(); - // int evaluateAdjointJacobian(); - int evaluateAdjointIntegrand(); + int initializeAdjoint(); + int evaluateAdjointResidual(); + // int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); - CircuitComponent* clone() const; + Component* clone() const; - private: - RealT RN_; - NodeT* node1_; - }; + private: + RealT RN_; + NodeT* node1_; + }; + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp index 7dc2f6b80..4ac28f489 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp @@ -7,193 +7,194 @@ namespace GridKit { - - /*! - * @brief Constructor for a constant MicrogridLine model - * - * Calls default ModelEvaluatorImpl constructor. - * - * - * Model is from paper: "Modeling, Analysis and Testing of Autonomous Operation - * of an Inverter-Based Microgrid", Nagaraju Pogaku, Milan Prodanovic, and - * Timothy C. Green, Section C - * - * @todo Consider having \omegaref as a global constant, not a node variable. - */ - - template - MicrogridLine::MicrogridLine(IdxT id, RealT R, RealT L, NodeT* node_ref, NodeT* bus1, NodeT* bus2) - : R_(R), - L_(L), - node_ref_(node_ref), - bus1_(bus1), - bus2_(bus2) - { - assert(node_ref_->size() == 1); - assert(bus1_->size() == 2); - assert(bus2_->size() == 2); - // internals [id, iq] - // externals [\omegaref, vbd_in, vbq_in, vbd_out, vbq_out] - size_ = 7; - n_intern_ = 2; - n_extern_ = 5; - extern_indices_ = {0, 1, 2, 3, 4}; - idc_ = id; - nnz_ = 14; - } - - template - MicrogridLine::~MicrogridLine() - { - } - - /** - * Initialization of the grid model - */ - template - int MicrogridLine::initialize() - { - return 0; - } - - /* - * \brief Identify differential variables - */ - template - int MicrogridLine::tagDifferentiable() - { - // All variables are differentials - std::fill(tag_.begin(), tag_.end(), true); - return 0; - } - - /** - * @brief Compute the absolute tolerance for each variable in the model - * - * @param rel_tol The relative tolerance which can be used to pick the - * absolute tolerance. - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return int 0 if successful, non-zero otherwise. - * - * This represents a "noise" level close to zero for which pure relative - * error cannot be used. - */ - template - int MicrogridLine::setAbsoluteTolerance(RealT rel_tol) - { - abs_tol_.setToConst(static_cast(rel_tol)); - return 0; - } - - /** - * @brief Evaluate residual of microgrid line - * - */ - template - int MicrogridLine::evaluateInternalResidual() - { - f_int_[0] = -yp_int_[0] - (R_ / L_) * y_int_[0] + *y_ext_[0] * y_int_[1] + (*y_ext_[1] - *y_ext_[3]) / L_; - f_int_[1] = -yp_int_[1] - (R_ / L_) * y_int_[1] - *y_ext_[0] * y_int_[0] + (*y_ext_[2] - *y_ext_[4]) / L_; - - return 0; - } - - template - int MicrogridLine::evaluateExternalResidual() - { - // Port 1 - *f_ext_[1] += -y_int_[0]; - *f_ext_[2] += -y_int_[1]; - - // Port 2 - *f_ext_[3] += y_int_[0]; - *f_ext_[4] += y_int_[1]; - - return 0; - } - - /** - * @brief Generate Jacobian for Microgrid Line - * - * @tparam ScalarT - * @tparam IdxT - * @return int - */ - template - int MicrogridLine::evaluateJacobian() - { - this->zeroJacMatrix(); - - // Create dF/dy - std::vector rtemp{1, 2, 3, 4}; - std::vector ctemp{5, 6, 5, 6}; - std::vector valtemp{-1.0, -1.0, 1.0, 1.0}; - this->setJacValues(rtemp, ctemp, valtemp); - - std::vector ccord{0, 1, 3, 5, 6}; - - std::vector rcord(ccord.size(), 5); - std::vector vals{}; - vals = {static_cast(y_int_[1]), (1.0 / L_), -(1.0 / L_), -(R_ / L_) - alpha_, static_cast(*y_ext_[0])}; - this->setJacValues(rcord, ccord, vals); - - std::vector ccor2{0, 2, 4, 5, 6}; - std::fill(rcord.begin(), rcord.end(), 6); - vals = {-static_cast(y_int_[0]), (1.0 / L_), -(1.0 / L_), -static_cast(*y_ext_[0]), -(R_ / L_) - alpha_}; - this->setJacValues(rcord, ccor2, vals); - - return 0; - } - - template - int MicrogridLine::allocate() - { - CircuitComponent::allocate(); - - this->setExternalConnectionNodes(0, node_ref_->getNodeConnection(0)); - this->setExternalConnectionNodes(1, bus1_->getNodeConnection(0)); - this->setExternalConnectionNodes(2, bus1_->getNodeConnection(1)); - this->setExternalConnectionNodes(3, bus2_->getNodeConnection(0)); - this->setExternalConnectionNodes(4, bus2_->getNodeConnection(1)); - - return 0; - } - - template - int MicrogridLine::evaluateIntegrand() - { - return 0; - } - - template - int MicrogridLine::initializeAdjoint() + namespace PowerElectronics { - return 0; - } - - template - int MicrogridLine::evaluateAdjointResidual() - { - return 0; - } - - template - int MicrogridLine::evaluateAdjointIntegrand() - { - return 0; - } - - template - CircuitComponent* MicrogridLine::clone() const - { - return new MicrogridLine(*this); - } - - // Available template instantiations - template class MicrogridLine; - template class MicrogridLine; - template class MicrogridLine; - template class MicrogridLine; - + /*! + * @brief Constructor for a constant MicrogridLine model + * + * Calls default ModelEvaluatorImpl constructor. + * + * + * Model is from paper: "Modeling, Analysis and Testing of Autonomous Operation + * of an Inverter-Based Microgrid", Nagaraju Pogaku, Milan Prodanovic, and + * Timothy C. Green, Section C + * + * @todo Consider having \omegaref as a global constant, not a node variable. + */ + template + MicrogridLine::MicrogridLine(IdxT id, RealT R, RealT L, NodeT* node_ref, NodeT* bus1, NodeT* bus2) + : R_(R), + L_(L), + node_ref_(node_ref), + bus1_(bus1), + bus2_(bus2) + { + assert(node_ref_->size() == 1); + assert(bus1_->size() == 2); + assert(bus2_->size() == 2); + // internals [id, iq] + // externals [\omegaref, vbd_in, vbq_in, vbd_out, vbq_out] + size_ = 7; + n_intern_ = 2; + n_extern_ = 5; + extern_indices_ = {0, 1, 2, 3, 4}; + idc_ = id; + nnz_ = 14; + } + + template + MicrogridLine::~MicrogridLine() + { + } + + /** + * Initialization of the grid model + */ + template + int MicrogridLine::initialize() + { + return 0; + } + + /* + * \brief Identify differential variables + */ + template + int MicrogridLine::tagDifferentiable() + { + // All variables are differentials + std::fill(tag_.begin(), tag_.end(), true); + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int MicrogridLine::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Evaluate residual of microgrid line + * + */ + template + int MicrogridLine::evaluateInternalResidual() + { + f_int_[0] = -yp_int_[0] - (R_ / L_) * y_int_[0] + *y_ext_[0] * y_int_[1] + (*y_ext_[1] - *y_ext_[3]) / L_; + f_int_[1] = -yp_int_[1] - (R_ / L_) * y_int_[1] - *y_ext_[0] * y_int_[0] + (*y_ext_[2] - *y_ext_[4]) / L_; + + return 0; + } + + template + int MicrogridLine::evaluateExternalResidual() + { + // Port 1 + *f_ext_[1] += -y_int_[0]; + *f_ext_[2] += -y_int_[1]; + + // Port 2 + *f_ext_[3] += y_int_[0]; + *f_ext_[4] += y_int_[1]; + + return 0; + } + + /** + * @brief Generate Jacobian for Microgrid Line + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ + template + int MicrogridLine::evaluateJacobian() + { + this->zeroJacMatrix(); + + // Create dF/dy + std::vector rtemp{1, 2, 3, 4}; + std::vector ctemp{5, 6, 5, 6}; + std::vector valtemp{-1.0, -1.0, 1.0, 1.0}; + this->setJacValues(rtemp, ctemp, valtemp); + + std::vector ccord{0, 1, 3, 5, 6}; + + std::vector rcord(ccord.size(), 5); + std::vector vals{}; + vals = {static_cast(y_int_[1]), (1.0 / L_), -(1.0 / L_), -(R_ / L_) - alpha_, static_cast(*y_ext_[0])}; + this->setJacValues(rcord, ccord, vals); + + std::vector ccor2{0, 2, 4, 5, 6}; + std::fill(rcord.begin(), rcord.end(), 6); + vals = {-static_cast(y_int_[0]), (1.0 / L_), -(1.0 / L_), -static_cast(*y_ext_[0]), -(R_ / L_) - alpha_}; + this->setJacValues(rcord, ccor2, vals); + + return 0; + } + + template + int MicrogridLine::allocate() + { + Component::allocate(); + + this->setExternalConnectionNodes(0, node_ref_->getNodeConnection(0)); + this->setExternalConnectionNodes(1, bus1_->getNodeConnection(0)); + this->setExternalConnectionNodes(2, bus1_->getNodeConnection(1)); + this->setExternalConnectionNodes(3, bus2_->getNodeConnection(0)); + this->setExternalConnectionNodes(4, bus2_->getNodeConnection(1)); + + return 0; + } + + template + int MicrogridLine::evaluateIntegrand() + { + return 0; + } + + template + int MicrogridLine::initializeAdjoint() + { + return 0; + } + + template + int MicrogridLine::evaluateAdjointResidual() + { + return 0; + } + + template + int MicrogridLine::evaluateAdjointIntegrand() + { + return 0; + } + + template + Component* MicrogridLine::clone() const + { + return new MicrogridLine(*this); + } + + // Available template instantiations + template class MicrogridLine; + template class MicrogridLine; + template class MicrogridLine; + template class MicrogridLine; + + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp index 0749e5ceb..dd520cb45 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp @@ -1,77 +1,73 @@ - #pragma once -#include +#include #include namespace GridKit { - template - class BaseBus; -} - -namespace GridKit -{ - /*! - * @brief Declaration of a MicrogridLine class. - * - */ - template - class MicrogridLine : public CircuitComponent + namespace PowerElectronics { - using RealT = typename CircuitComponent::RealT; - using NodeT = typename PowerElectronics::NodeBase; + /*! + * @brief Declaration of a MicrogridLine class. + * + */ + template + class MicrogridLine : public Component + { + using RealT = typename Component::RealT; + using NodeT = typename PowerElectronics::NodeBase; - using CircuitComponent::size_; - using CircuitComponent::nnz_; - using CircuitComponent::time_; - using CircuitComponent::alpha_; - using CircuitComponent::y_ext_; - using CircuitComponent::y_int_; - using CircuitComponent::yp_ext_; - using CircuitComponent::yp_int_; - using CircuitComponent::tag_; - using CircuitComponent::abs_tol_; - using CircuitComponent::f_ext_; - using CircuitComponent::f_int_; - using CircuitComponent::g_; - using CircuitComponent::yB_; - using CircuitComponent::ypB_; - using CircuitComponent::fB_; - using CircuitComponent::gB_; - using CircuitComponent::param_; - using CircuitComponent::idc_; + using Component::size_; + using Component::nnz_; + using Component::time_; + using Component::alpha_; + using Component::y_ext_; + using Component::y_int_; + using Component::yp_ext_; + using Component::yp_int_; + using Component::tag_; + using Component::abs_tol_; + using Component::f_ext_; + using Component::f_int_; + using Component::g_; + using Component::yB_; + using Component::ypB_; + using Component::fB_; + using Component::gB_; + using Component::param_; + using Component::idc_; - using CircuitComponent::extern_indices_; - using CircuitComponent::n_extern_; - using CircuitComponent::n_intern_; + using Component::extern_indices_; + using Component::n_extern_; + using Component::n_intern_; - public: - MicrogridLine(IdxT id, RealT R, RealT L, NodeT* node_ref, NodeT* bus1, NodeT* bus2); - virtual ~MicrogridLine(); + public: + MicrogridLine(IdxT id, RealT R, RealT L, NodeT* node_ref, NodeT* bus1, NodeT* bus2); + virtual ~MicrogridLine(); - int initialize(); - int allocate() final; - int tagDifferentiable(); - int setAbsoluteTolerance(RealT); - int evaluateInternalResidual() final; - int evaluateExternalResidual() final; - int evaluateJacobian(); - int evaluateIntegrand(); + int initialize(); + int allocate() final; + int tagDifferentiable(); + int setAbsoluteTolerance(RealT); + int evaluateInternalResidual() final; + int evaluateExternalResidual() final; + int evaluateJacobian(); + int evaluateIntegrand(); - int initializeAdjoint(); - int evaluateAdjointResidual(); - // int evaluateAdjointJacobian(); - int evaluateAdjointIntegrand(); + int initializeAdjoint(); + int evaluateAdjointResidual(); + // int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); - CircuitComponent* clone() const; + Component* clone() const; - private: - RealT R_; - RealT L_; - NodeT* node_ref_; - NodeT* bus1_; - NodeT* bus2_; - }; + private: + RealT R_; + RealT L_; + NodeT* node_ref_; + NodeT* bus1_; + NodeT* bus2_; + }; + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp index 34de349d1..1d0e07f63 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp @@ -7,184 +7,185 @@ namespace GridKit { - - /*! - * @brief Constructor for a constant MicrogridLoad model - * - * Calls default ModelEvaluatorImpl constructor. - * - * - * Model is from paper: " - "Modeling, Analysis and Testing of Autonomous Operation of an Inverter-Based Microgrid" Nagaraju Pogaku, Milan Prodanovic, and Timothy C. Green" - * Section D - */ - - template - MicrogridLoad::MicrogridLoad(IdxT id, RealT R, RealT L, NodeT* node_ref, NodeT* node_bus) - : R_(R), - L_(L), - node_ref_(node_ref), - node_bus_(node_bus) - { - assert(node_ref_->size() == 1); - assert(node_bus_->size() == 2); - // internals [id, iq] - // externals [\omegaref, vbd_out, vbq_out] - size_ = 5; - n_intern_ = 2; - n_extern_ = 3; - extern_indices_ = {0, 1, 2}; - idc_ = id; - nnz_ = 10; - } - - template - MicrogridLoad::~MicrogridLoad() - { - } - - /** - * Initialization of the grid model - */ - template - int MicrogridLoad::initialize() - { - return 0; - } - - /* - * \brief Identify differential variables - */ - template - int MicrogridLoad::tagDifferentiable() - { - // All variables are differentials - std::fill(tag_.begin(), tag_.end(), true); - return 0; - } - - /** - * @brief Compute the absolute tolerance for each variable in the model - * - * @param rel_tol The relative tolerance which can be used to pick the - * absolute tolerance. - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return int 0 if successful, non-zero otherwise. - * - * This represents a "noise" level close to zero for which pure relative - * error cannot be used. - */ - template - int MicrogridLoad::setAbsoluteTolerance(RealT rel_tol) - { - abs_tol_.setToConst(static_cast(rel_tol)); - return 0; - } - - /** - * @brief Eval Micro Load - */ - template - int MicrogridLoad::evaluateInternalResidual() - { - f_int_[0] = -yp_int_[0] - (R_ / L_) * y_int_[0] + *y_ext_[0] * y_int_[1] + *y_ext_[1] / L_; - f_int_[1] = -yp_int_[1] - (R_ / L_) * y_int_[1] - *y_ext_[0] * y_int_[0] + *y_ext_[2] / L_; - - return 0; - } - - template - int MicrogridLoad::evaluateExternalResidual() - { - // only input for loads - - // input - *f_ext_[1] += -y_int_[0]; - *f_ext_[2] += -y_int_[1]; - - return 0; - } - - /** - * @brief Generate Jacobian for Micro Load - * - * @tparam ScalarT - * @tparam IdxT - * @return int - */ - template - int MicrogridLoad::evaluateJacobian() - { - this->zeroJacMatrix(); - - // Create dF/dy - std::vector rtemp{1, 2}; - std::vector ctemp{3, 4}; - std::vector valtemp{-1.0, -1.0}; - this->setJacValues(rtemp, ctemp, valtemp); - - std::vector ccord{0, 1, 3, 4}; - - std::vector rcord(ccord.size(), 3); - std::vector vals{}; - vals = {static_cast(y_int_[1]), (1.0 / L_), -(R_ / L_) - alpha_, static_cast(*y_ext_[0])}; - this->setJacValues(rcord, ccord, vals); - - std::vector ccor2{0, 2, 3, 4}; - std::fill(rcord.begin(), rcord.end(), 4); - vals = {-static_cast(y_int_[0]), (1.0 / L_), -static_cast(*y_ext_[0]), -(R_ / L_) - alpha_}; - this->setJacValues(rcord, ccor2, vals); - - return 0; - } - - template - int MicrogridLoad::allocate() - { - CircuitComponent::allocate(); - - this->setExternalConnectionNodes(0, node_ref_->getNodeConnection(0)); - this->setExternalConnectionNodes(1, node_bus_->getNodeConnection(0)); - this->setExternalConnectionNodes(2, node_bus_->getNodeConnection(1)); - - return 0; - } - - template - int MicrogridLoad::evaluateIntegrand() - { - return 0; - } - - template - int MicrogridLoad::initializeAdjoint() + namespace PowerElectronics { - return 0; - } - - template - int MicrogridLoad::evaluateAdjointResidual() - { - return 0; - } - - template - int MicrogridLoad::evaluateAdjointIntegrand() - { - return 0; - } - - template - CircuitComponent* MicrogridLoad::clone() const - { - return new MicrogridLoad(*this); - } - - // Available template instantiations - template class MicrogridLoad; - template class MicrogridLoad; - template class MicrogridLoad; - template class MicrogridLoad; - + /*! + * @brief Constructor for a constant MicrogridLoad model + * + * Calls default ModelEvaluatorImpl constructor. + * + * + * Model is from paper: " + "Modeling, Analysis and Testing of Autonomous Operation of an Inverter-Based Microgrid" Nagaraju Pogaku, Milan Prodanovic, and Timothy C. Green" + * Section D + */ + template + MicrogridLoad::MicrogridLoad(IdxT id, RealT R, RealT L, NodeT* node_ref, NodeT* node_bus) + : R_(R), + L_(L), + node_ref_(node_ref), + node_bus_(node_bus) + { + assert(node_ref_->size() == 1); + assert(node_bus_->size() == 2); + // internals [id, iq] + // externals [\omegaref, vbd_out, vbq_out] + size_ = 5; + n_intern_ = 2; + n_extern_ = 3; + extern_indices_ = {0, 1, 2}; + idc_ = id; + nnz_ = 10; + } + + template + MicrogridLoad::~MicrogridLoad() + { + } + + /** + * Initialization of the grid model + */ + template + int MicrogridLoad::initialize() + { + return 0; + } + + /* + * \brief Identify differential variables + */ + template + int MicrogridLoad::tagDifferentiable() + { + // All variables are differentials + std::fill(tag_.begin(), tag_.end(), true); + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int MicrogridLoad::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Eval Micro Load + */ + template + int MicrogridLoad::evaluateInternalResidual() + { + f_int_[0] = -yp_int_[0] - (R_ / L_) * y_int_[0] + *y_ext_[0] * y_int_[1] + *y_ext_[1] / L_; + f_int_[1] = -yp_int_[1] - (R_ / L_) * y_int_[1] - *y_ext_[0] * y_int_[0] + *y_ext_[2] / L_; + + return 0; + } + + template + int MicrogridLoad::evaluateExternalResidual() + { + // only input for loads + + // input + *f_ext_[1] += -y_int_[0]; + *f_ext_[2] += -y_int_[1]; + + return 0; + } + + /** + * @brief Generate Jacobian for Micro Load + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ + template + int MicrogridLoad::evaluateJacobian() + { + this->zeroJacMatrix(); + + // Create dF/dy + std::vector rtemp{1, 2}; + std::vector ctemp{3, 4}; + std::vector valtemp{-1.0, -1.0}; + this->setJacValues(rtemp, ctemp, valtemp); + + std::vector ccord{0, 1, 3, 4}; + + std::vector rcord(ccord.size(), 3); + std::vector vals{}; + vals = {static_cast(y_int_[1]), (1.0 / L_), -(R_ / L_) - alpha_, static_cast(*y_ext_[0])}; + this->setJacValues(rcord, ccord, vals); + + std::vector ccor2{0, 2, 3, 4}; + std::fill(rcord.begin(), rcord.end(), 4); + vals = {-static_cast(y_int_[0]), (1.0 / L_), -static_cast(*y_ext_[0]), -(R_ / L_) - alpha_}; + this->setJacValues(rcord, ccor2, vals); + + return 0; + } + + template + int MicrogridLoad::allocate() + { + Component::allocate(); + + this->setExternalConnectionNodes(0, node_ref_->getNodeConnection(0)); + this->setExternalConnectionNodes(1, node_bus_->getNodeConnection(0)); + this->setExternalConnectionNodes(2, node_bus_->getNodeConnection(1)); + + return 0; + } + + template + int MicrogridLoad::evaluateIntegrand() + { + return 0; + } + + template + int MicrogridLoad::initializeAdjoint() + { + return 0; + } + + template + int MicrogridLoad::evaluateAdjointResidual() + { + return 0; + } + + template + int MicrogridLoad::evaluateAdjointIntegrand() + { + return 0; + } + + template + Component* MicrogridLoad::clone() const + { + return new MicrogridLoad(*this); + } + + // Available template instantiations + template class MicrogridLoad; + template class MicrogridLoad; + template class MicrogridLoad; + template class MicrogridLoad; + + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp index 12935d671..b5759d609 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp @@ -1,76 +1,72 @@ - #pragma once -#include +#include #include namespace GridKit { - template - class BaseBus; -} - -namespace GridKit -{ - /*! - * @brief Declaration of a passive MicrogridLoad class. - * - */ - template - class MicrogridLoad : public CircuitComponent + namespace PowerElectronics { - using RealT = typename CircuitComponent::RealT; - using NodeT = typename PowerElectronics::NodeBase; + /*! + * @brief Declaration of a passive MicrogridLoad class. + * + */ + template + class MicrogridLoad : public Component + { + using RealT = typename Component::RealT; + using NodeT = typename PowerElectronics::NodeBase; - using CircuitComponent::size_; - using CircuitComponent::nnz_; - using CircuitComponent::time_; - using CircuitComponent::alpha_; - using CircuitComponent::y_ext_; - using CircuitComponent::y_int_; - using CircuitComponent::yp_ext_; - using CircuitComponent::yp_int_; - using CircuitComponent::tag_; - using CircuitComponent::abs_tol_; - using CircuitComponent::f_ext_; - using CircuitComponent::f_int_; - using CircuitComponent::g_; - using CircuitComponent::yB_; - using CircuitComponent::ypB_; - using CircuitComponent::fB_; - using CircuitComponent::gB_; - using CircuitComponent::param_; - using CircuitComponent::idc_; + using Component::size_; + using Component::nnz_; + using Component::time_; + using Component::alpha_; + using Component::y_ext_; + using Component::y_int_; + using Component::yp_ext_; + using Component::yp_int_; + using Component::tag_; + using Component::abs_tol_; + using Component::f_ext_; + using Component::f_int_; + using Component::g_; + using Component::yB_; + using Component::ypB_; + using Component::fB_; + using Component::gB_; + using Component::param_; + using Component::idc_; - using CircuitComponent::extern_indices_; - using CircuitComponent::n_extern_; - using CircuitComponent::n_intern_; + using Component::extern_indices_; + using Component::n_extern_; + using Component::n_intern_; - public: - MicrogridLoad(IdxT id, RealT R, RealT L, NodeT* node_ref, NodeT* node_bus); - virtual ~MicrogridLoad(); + public: + MicrogridLoad(IdxT id, RealT R, RealT L, NodeT* node_ref, NodeT* node_bus); + virtual ~MicrogridLoad(); - int initialize(); - int allocate() final; - int tagDifferentiable(); - int setAbsoluteTolerance(RealT); - int evaluateInternalResidual() final; - int evaluateExternalResidual() final; - int evaluateJacobian(); - int evaluateIntegrand(); + int initialize(); + int allocate() final; + int tagDifferentiable(); + int setAbsoluteTolerance(RealT); + int evaluateInternalResidual() final; + int evaluateExternalResidual() final; + int evaluateJacobian(); + int evaluateIntegrand(); - int initializeAdjoint(); - int evaluateAdjointResidual(); - // int evaluateAdjointJacobian(); - int evaluateAdjointIntegrand(); + int initializeAdjoint(); + int evaluateAdjointResidual(); + // int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); - CircuitComponent* clone() const; + Component* clone() const; - private: - RealT R_; - RealT L_; - NodeT* node_ref_; - NodeT* node_bus_; - }; + private: + RealT R_; + RealT L_; + NodeT* node_ref_; + NodeT* node_bus_; + }; + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/Node.hpp b/GridKit/Model/PowerElectronics/Node.hpp new file mode 100644 index 000000000..c9cd1b57b --- /dev/null +++ b/GridKit/Model/PowerElectronics/Node.hpp @@ -0,0 +1,385 @@ +#pragma once + +#include +#include +#include +#include + +#include + +namespace GridKit +{ + namespace PowerElectronics + { + /** + * @brief Circuit node representing a connection point. + */ + template + class Node : public Model::Evaluator + { + using RealT = typename Model::Evaluator::RealT; + using VectorT = typename Model::Evaluator::VectorT; + + public: + Node() + { + size_ = 1; + } + + Node(ScalarT v0) + : V0_(v0) + { + size_ = 1; + } + + ~Node() = default; + + int setNodeID(IdxT id) + { + id_ = id; + return 0; + } + + IdxT nodeID() const + { + return id_; + } + + // Voltage accessor + ScalarT& V() + { + return y_.getData()[0]; + } + + const ScalarT& V() const + { + return y_.getData()[0]; + } + + // KCL residual accessor + ScalarT& I() + { + return f_.getData()[0]; + } + + const ScalarT& I() const + { + return f_.getData()[0]; + } + + // Allocate storage for a single-node voltage and KCL residual + int allocate() + { + size_t size = static_cast(size_); + + if (!allocated_) + { + allocateVectors(size_); + } + + tag_.resize(size); + + variable_indices_[0] = 0; + residual_indices_[0] = 0; + + allocated_ = true; + return 0; + } + + /** + * @brief Initialize node variables + */ + int initialize() + { + auto* y = y_.getData(); + auto* yp = yp_.getData(); + + y[0] = V0_; + yp[0] = 0.0; + + y_.setDataUpdated(); + yp_.setDataUpdated(); + + return 0; + } + + /** + * @brief Node variables are algebraic. + */ + int tagDifferentiable() + { + tag_[0] = false; + + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + int setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Node does not compute residuals, so here we just reset residual values. + * + * @warning This implementation assumes node residuals are always evaluated + * _before_ component model residuals. + * + */ + int evaluateResidual() + { + auto* f = f_.getData(); + + f[0] = 0.0; + + f_.setDataUpdated(); + + return 0; + } + + bool hasJacobian() final + { + return false; + } + + /** + * @brief There is no Jacobian for node variables + */ + int evaluateJacobian() + { + return 0; + } + + int evaluateIntegrand() + { + return 0; + } + + int initializeAdjoint() + { + return 0; + } + + int evaluateAdjointResidual() + { + return 0; + } + + int evaluateAdjointIntegrand() + { + return 0; + } + + private: + IdxT id_{static_cast(-1)}; + IdxT size_{0}; + IdxT nnz_{0}; + IdxT size_quad_{0}; + IdxT size_opt_{0}; + ScalarT V0_{0.0}; + + std::map variable_indices_; + std::map residual_indices_; + + VectorT y_; + VectorT yp_; + std::vector tag_; + VectorT abs_tol_; + VectorT f_; + + VectorT g_{}; + VectorT param_{}; + VectorT param_up_{}; + VectorT param_lo_{}; + + VectorT yB_{}; + VectorT ypB_{}; + VectorT fB_{}; + VectorT gB_{}; + + RealT time_{0}; + RealT alpha_{0}; + + IdxT max_steps_{0}; + + bool allocated_{false}; + + public: + IdxT size() final + { + return size_; + } + + IdxT nnz() final + { + return nnz_; + } + + IdxT sizeQuadrature() final + { + return size_quad_; + } + + IdxT sizeParams() final + { + return size_opt_; + } + + void updateTime(RealT /* t */, RealT /* a */) final + { + // No time to update in node models + } + + VectorT& y() final + { + return y_; + } + + const VectorT& y() const final + { + return y_; + } + + VectorT& yp() final + { + return yp_; + } + + const VectorT& yp() const final + { + return yp_; + } + + std::vector& tag() final + { + return tag_; + } + + const std::vector& tag() const final + { + return tag_; + } + + VectorT& absoluteTolerance() final + { + return abs_tol_; + } + + const VectorT& absoluteTolerance() const final + { + return abs_tol_; + } + + VectorT& yB() final + { + return yB_; + } + + const VectorT& yB() const final + { + return yB_; + } + + VectorT& ypB() final + { + return ypB_; + } + + const VectorT& ypB() const final + { + return ypB_; + } + + VectorT& param() final + { + return param_; + } + + const VectorT& param() const final + { + return param_; + } + + VectorT& param_up() final + { + return param_up_; + } + + const VectorT& param_up() const final + { + return param_up_; + } + + VectorT& param_lo() final + { + return param_lo_; + } + + const VectorT& param_lo() const final + { + return param_lo_; + } + + VectorT& getResidual() final + { + return f_; + } + + const VectorT& getResidual() const final + { + return f_; + } + + VectorT& getIntegrand() final + { + return g_; + } + + const VectorT& getIntegrand() const final + { + return g_; + } + + VectorT& getAdjointResidual() final + { + return fB_; + } + + const VectorT& getAdjointResidual() const final + { + return fB_; + } + + VectorT& getAdjointIntegrand() final + { + return gB_; + } + + const VectorT& getAdjointIntegrand() const final + { + return gB_; + } + + private: + void allocateVectors(IdxT n) + { + y_.resize(n); + yp_.resize(n); + f_.resize(n); + abs_tol_.resize(n); + } + }; + } // namespace PowerElectronics +} // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/NodeBase.hpp b/GridKit/Model/PowerElectronics/NodeBase.hpp index b39539781..6e95322f1 100644 --- a/GridKit/Model/PowerElectronics/NodeBase.hpp +++ b/GridKit/Model/PowerElectronics/NodeBase.hpp @@ -104,7 +104,7 @@ namespace GridKit /** * @brief Create the mappings from local to global indices for a node variable (either internal or external), - * to be used from an attached component. \see CircuitComponent::setExternalConnectionNodes() + * to be used from an attached component. \see Component::setExternalConnectionNodes() * * @param local_index The index of the local variable * @param connection The necessary connection information for the variable diff --git a/GridKit/Model/PowerElectronics/README.md b/GridKit/Model/PowerElectronics/README.md index 669e82172..c4f88acec 100644 --- a/GridKit/Model/PowerElectronics/README.md +++ b/GridKit/Model/PowerElectronics/README.md @@ -1,4 +1,4 @@ -This sub-directory provides components utilized within the PowerElectronicsModel composer. All components are treated equally and only the composer can view and distribute data to components. Components have no knowledge nor require the existence of any other component. +This sub-directory provides components utilized within the PowerElectronics::SystemModel composer. All components are treated equally and only the composer can view and distribute data to components. Components have no knowledge nor require the existence of any other component. Each component evaluates a set of equations to form residuals. There are two types of variables. + Internal Variables $y$. Variables only need by the component. + External Variables $w$. Variables shared between multiple components. diff --git a/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp b/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp index 75b3ec18b..0284bafaf 100644 --- a/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp +++ b/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp @@ -1,5 +1,4 @@ - #include "Resistor.hpp" #include @@ -8,147 +7,148 @@ namespace GridKit { - - /*! - * @brief Constructor for a resistor model - * - * Calls default ModelEvaluatorImpl constructor. - */ - - template - Resistor::Resistor(IdxT id, RealT R, NodeT* node1, NodeT* node2) - : R_(R), node1_(node1), node2_(node2) - { - assert(node1_->size() == 1); - assert(node2_->size() == 1); - size_ = 2; - n_intern_ = 0; - n_extern_ = 2; - extern_indices_ = {0, 1}; - idc_ = id; - nnz_ = 4; - } - - template - Resistor::~Resistor() - { - } - - /** - * Initialization of the grid model - */ - template - int Resistor::initialize() - { - return 0; - } - - /// There are no internal variables in this component, so \ref tag_ can be set arbitrarily. - template - int Resistor::tagDifferentiable() + namespace PowerElectronics { - return 0; - } - - /** - * @brief Compute the absolute tolerance for each variable in the model - * - * @param rel_tol The relative tolerance which can be used to pick the - * absolute tolerance. - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return int 0 if successful, non-zero otherwise. - * - * This represents a "noise" level close to zero for which pure relative - * error cannot be used. - */ - template - int Resistor::setAbsoluteTolerance(RealT rel_tol) - { - abs_tol_.setToConst(static_cast(rel_tol)); - return 0; - } - - /** - * @brief Computes the resistors resisdual - * - */ - template - int Resistor::evaluateInternalResidual() - { - return 0; - } - - template - int Resistor::evaluateExternalResidual() - { - // input - *f_ext_[0] += (*y_ext_[0] - *y_ext_[1]) / R_; - // ouput - *f_ext_[1] += (*y_ext_[1] - *y_ext_[0]) / R_; - return 0; - } - - template - int Resistor::evaluateJacobian() - { - this->zeroJacMatrix(); - - // Create dF/dy - // does compiler make constant??? - std::vector rcord{0, 0, 1, 1}; - std::vector ccord{0, 1, 0, 1}; - std::vector vals{1.0 / R_, -1.0 / R_, -1.0 / R_, 1.0 / R_}; - this->setJacValues(rcord, ccord, vals); - - return 0; - } - - template - int Resistor::allocate() - { - CircuitComponent::allocate(); - - this->setExternalConnectionNodes(0, node1_->getNodeConnection(0)); - this->setExternalConnectionNodes(1, node2_->getNodeConnection(0)); - - return 0; - } - - template - int Resistor::evaluateIntegrand() - { - return 0; - } - - template - int Resistor::initializeAdjoint() - { - return 0; - } - - template - int Resistor::evaluateAdjointResidual() - { - return 0; - } - - template - int Resistor::evaluateAdjointIntegrand() - { - return 0; - } - - template - CircuitComponent* Resistor::clone() const - { - return new Resistor(*this); - } - - // Available template instantiations - template class Resistor; - template class Resistor; - template class Resistor; - template class Resistor; - + /*! + * @brief Constructor for a resistor model + * + * Calls default ModelEvaluatorImpl constructor. + */ + template + Resistor::Resistor(IdxT id, RealT R, NodeT* node1, NodeT* node2) + : R_(R), node1_(node1), node2_(node2) + { + assert(node1_->size() == 1); + assert(node2_->size() == 1); + size_ = 2; + n_intern_ = 0; + n_extern_ = 2; + extern_indices_ = {0, 1}; + idc_ = id; + nnz_ = 4; + } + + template + Resistor::~Resistor() + { + } + + /** + * Initialization of the grid model + */ + template + int Resistor::initialize() + { + return 0; + } + + /// There are no internal variables in this component, so \ref tag_ can be set arbitrarily. + template + int Resistor::tagDifferentiable() + { + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Resistor::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Computes the resistors resisdual + * + */ + template + int Resistor::evaluateInternalResidual() + { + return 0; + } + + template + int Resistor::evaluateExternalResidual() + { + // input + *f_ext_[0] += (*y_ext_[0] - *y_ext_[1]) / R_; + // ouput + *f_ext_[1] += (*y_ext_[1] - *y_ext_[0]) / R_; + return 0; + } + + template + int Resistor::evaluateJacobian() + { + this->zeroJacMatrix(); + + // Create dF/dy + // does compiler make constant??? + std::vector rcord{0, 0, 1, 1}; + std::vector ccord{0, 1, 0, 1}; + std::vector vals{1.0 / R_, -1.0 / R_, -1.0 / R_, 1.0 / R_}; + this->setJacValues(rcord, ccord, vals); + + return 0; + } + + template + int Resistor::allocate() + { + Component::allocate(); + + this->setExternalConnectionNodes(0, node1_->getNodeConnection(0)); + this->setExternalConnectionNodes(1, node2_->getNodeConnection(0)); + + return 0; + } + + template + int Resistor::evaluateIntegrand() + { + return 0; + } + + template + int Resistor::initializeAdjoint() + { + return 0; + } + + template + int Resistor::evaluateAdjointResidual() + { + return 0; + } + + template + int Resistor::evaluateAdjointIntegrand() + { + return 0; + } + + template + Component* Resistor::clone() const + { + return new Resistor(*this); + } + + // Available template instantiations + template class Resistor; + template class Resistor; + template class Resistor; + template class Resistor; + + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp b/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp index 354074bc9..0be29459e 100644 --- a/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp +++ b/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp @@ -1,75 +1,71 @@ - #pragma once -#include +#include #include namespace GridKit { - template - class BaseBus; -} - -namespace GridKit -{ - /*! - * @brief Declaration of a Resistor class. - * - */ - template - class Resistor : public CircuitComponent + namespace PowerElectronics { - using RealT = typename CircuitComponent::RealT; - using NodeT = typename PowerElectronics::NodeBase; + /*! + * @brief Declaration of a Resistor class. + * + */ + template + class Resistor : public Component + { + using RealT = typename Component::RealT; + using NodeT = typename PowerElectronics::NodeBase; - using CircuitComponent::size_; - using CircuitComponent::nnz_; - using CircuitComponent::time_; - using CircuitComponent::alpha_; - using CircuitComponent::y_ext_; - using CircuitComponent::y_int_; - using CircuitComponent::yp_ext_; - using CircuitComponent::yp_int_; - using CircuitComponent::tag_; - using CircuitComponent::abs_tol_; - using CircuitComponent::f_ext_; - using CircuitComponent::f_int_; - using CircuitComponent::g_; - using CircuitComponent::yB_; - using CircuitComponent::ypB_; - using CircuitComponent::fB_; - using CircuitComponent::gB_; - using CircuitComponent::param_; - using CircuitComponent::idc_; + using Component::size_; + using Component::nnz_; + using Component::time_; + using Component::alpha_; + using Component::y_ext_; + using Component::y_int_; + using Component::yp_ext_; + using Component::yp_int_; + using Component::tag_; + using Component::abs_tol_; + using Component::f_ext_; + using Component::f_int_; + using Component::g_; + using Component::yB_; + using Component::ypB_; + using Component::fB_; + using Component::gB_; + using Component::param_; + using Component::idc_; - using CircuitComponent::extern_indices_; - using CircuitComponent::n_extern_; - using CircuitComponent::n_intern_; + using Component::extern_indices_; + using Component::n_extern_; + using Component::n_intern_; - public: - Resistor(IdxT id, RealT R, NodeT* node1, NodeT* node2); - virtual ~Resistor(); + public: + Resistor(IdxT id, RealT R, NodeT* node1, NodeT* node2); + virtual ~Resistor(); - int initialize(); - int allocate() final; - int tagDifferentiable(); - int setAbsoluteTolerance(RealT); - int evaluateInternalResidual() final; - int evaluateExternalResidual() final; - int evaluateJacobian(); - int evaluateIntegrand(); + int initialize(); + int allocate() final; + int tagDifferentiable(); + int setAbsoluteTolerance(RealT); + int evaluateInternalResidual() final; + int evaluateExternalResidual() final; + int evaluateJacobian(); + int evaluateIntegrand(); - int initializeAdjoint(); - int evaluateAdjointResidual(); - // int evaluateAdjointJacobian(); - int evaluateAdjointIntegrand(); + int initializeAdjoint(); + int evaluateAdjointResidual(); + // int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); - CircuitComponent* clone() const; + Component* clone() const; - private: - RealT R_; - NodeT* node1_; - NodeT* node2_; - }; + private: + RealT R_; + NodeT* node1_; + NodeT* node2_; + }; + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp b/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp index 0616977e2..f1e1f6d8d 100644 --- a/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp +++ b/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp @@ -7,193 +7,194 @@ namespace GridKit { - - /*! - * @brief Constructor for a constant SynchronousMachine model - * - * Calls default ModelEvaluatorImpl constructor. - * @todo This model's equations are not finished - * @todo needs to be tested for correctness - * - * @tparam ScalarT - floating point type for the model - * @tparam IdxT - integer index type for the model - * - * @param[in] id - unique identifier for the component - * @param[in] Lls - stator leakage inductance - * @param[in] Llkq - tuple of damper leakage reactances - * @param[in] Llfd - field leakage reactance - * @param[in] Llkd - damper leakage reactance - * @param[in] Lmq - quadrature axis magnetizing reactance - * @param[in] Lmd - direct axis magnetizing reactance - * @param[in] Rs - stator resistance - * @param[in] Rkq - tuple of damper resistances - * @param[in] Rfd - field resistance - * @param[in] Rkd - damper resistance - * @param[in] RJ - rotor moment of inertia - * @param[in] P - number of poles - * @param[in] mub - rated frequency - */ - - template - SynchronousMachine::SynchronousMachine(IdxT id, RealT Lls, std::tuple Llkq, RealT Llfd, RealT Llkd, RealT Lmq, RealT Lmd, RealT Rs, std::tuple Rkq, RealT Rfd, RealT Rkd, RealT RJ, RealT P, RealT mub) - : Lls_(Lls), - Llkq_(Llkq), - Llfd_(Llfd), - Llkd_(Llkd), - Lmq_(Lmq), - Lmd_(Lmd), - Rs_(Rs), - Rkq_(Rkq), - Rfd_(Rfd), - Rkd_(Rkd), - RJ_(RJ), - P_(P), - mub_(mub) - { - size_ = 13; - n_intern_ = 6; - n_extern_ = 7; - extern_indices_ = {0, 1, 2, 3, 4}; - idc_ = id; - } - - template - SynchronousMachine::~SynchronousMachine() - { - } - - /** - * Initialization of the grid model - */ - template - int SynchronousMachine::initialize() - { - return 0; - } - - /* - * \brief Identify differential variables - */ - template - int SynchronousMachine::tagDifferentiable() - { - // All variables are differentials - std::fill(tag_.begin(), tag_.end(), true); - return 0; - } - - /** - * @brief Compute the absolute tolerance for each variable in the model - * - * @param rel_tol The relative tolerance which can be used to pick the - * absolute tolerance. - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return int 0 if successful, non-zero otherwise. - * - * This represents a "noise" level close to zero for which pure relative - * error cannot be used. - */ - template - int SynchronousMachine::setAbsoluteTolerance(RealT rel_tol) - { - abs_tol_.setToConst(static_cast(rel_tol)); - return 0; - } - - /** - * @brief Compute the resisdual of the component. - * - * @todo not finished - */ - template - int SynchronousMachine::evaluateInternalResidual() - { - ScalarT rkq1 = static_cast(std::get<0>(Rkq_)); - [[maybe_unused]] ScalarT rkq2 = static_cast(std::get<1>(Rkq_)); - ScalarT llkq1 = static_cast(std::get<0>(Llkq_)); - [[maybe_unused]] ScalarT llkq2 = static_cast(std::get<1>(Llkq_)); - - static constexpr auto pi = std::numbers::pi_v; - - ScalarT cos1 = std::cos((P_ / 2.0) * y_int_[0]); - ScalarT sin1 = std::sin((P_ / 2.0) * y_int_[0]); - ScalarT cos23m = std::cos((P_ / 2.0) * y_int_[0] - (2.0 / 3.0) * pi); - ScalarT sin23m = std::sin((P_ / 2.0) * y_int_[0] - (2.0 / 3.0) * pi); - ScalarT cos23p = std::cos((P_ / 2.0) * y_int_[0] + (2.0 / 3.0) * pi); - ScalarT sin23p = std::sin((P_ / 2.0) * y_int_[0] + (2.0 / 3.0) * pi); - - f_int_[0] = (-2.0 / 3.0) * (*y_ext_[0] * cos1 + *y_ext_[1] * cos23m + *y_ext_[2] * cos23p) + Rs_ * y_int_[1] + (Lls_ + Lmq_) * yp_int_[1] + Lmq_ * yp_int_[4] + Lmq_ * yp_int_[5] + *y_ext_[4] * (P_ / 2.0) * ((Lls_ + Lmd_) * y_int_[2] + Lmd_ * y_int_[6] + Lmd_ * y_int_[7]); - f_int_[1] = (-2.0 / 3.0) * (*y_ext_[0] * sin1 - *y_ext_[1] * sin23m - *y_ext_[2] * sin23p) + Rs_ * y_int_[2] + (Lls_ + Lmd_) * yp_int_[2] + Lmd_ * yp_int_[6] + Lmd_ * yp_int_[7] - *y_ext_[4] * (P_ / 2.0) * ((Lls_ + Lmq_) * y_int_[1] + Lmq_ * y_int_[4] + Lmq_ * y_int_[5]); - f_int_[2] = (-1.0 / 3.0) * (*y_ext_[0] + *y_ext_[1] + *y_ext_[2]) + Rs_ * y_int_[3] + Lls_ * yp_int_[3]; - f_int_[3] = rkq1 * y_int_[4] + (llkq1 + Lmq_) * yp_int_[4] + Lmq_ * yp_int_[1] + Lmq_ * yp_int_[5]; - f_int_[4] = rkq1 * y_int_[4] + (llkq1 + Lmq_) * yp_int_[4] + Lmq_ * yp_int_[1] + Lmq_ * yp_int_[5]; - return 0; - } - - template - int SynchronousMachine::evaluateExternalResidual() - { - [[maybe_unused]] ScalarT rkq2 = static_cast(std::get<1>(Rkq_)); - [[maybe_unused]] ScalarT llkq2 = static_cast(std::get<1>(Llkq_)); - - static constexpr auto pi = std::numbers::pi_v; - - ScalarT cos1 = std::cos((P_ / 2.0) * y_int_[0]); - ScalarT sin1 = std::sin((P_ / 2.0) * y_int_[0]); - ScalarT cos23m = std::cos((P_ / 2.0) * y_int_[0] - (2.0 / 3.0) * pi); - ScalarT sin23m = std::sin((P_ / 2.0) * y_int_[0] - (2.0 / 3.0) * pi); - ScalarT cos23p = std::cos((P_ / 2.0) * y_int_[0] + (2.0 / 3.0) * pi); - ScalarT sin23p = std::sin((P_ / 2.0) * y_int_[0] + (2.0 / 3.0) * pi); - - *f_ext_[0] += y_int_[1] * cos1 + y_int_[2] * sin1 + y_int_[3]; - *f_ext_[1] += y_int_[1] * cos23m + y_int_[2] * sin23m + y_int_[3]; - *f_ext_[2] += y_int_[1] * cos23p + y_int_[2] * sin23p + y_int_[3]; - *f_ext_[3] += RJ_ * *yp_ext_[4] - (3.0 / 4.0) * P_ * (Lmd_ * y_int_[1] * (y_int_[2] + y_int_[6] + y_int_[7]) - Lmq_ * y_int_[2] * (y_int_[1] + y_int_[4] + *y_ext_[0])); - *f_ext_[4] += yp_int_[0] - *y_ext_[4]; - return 0; - } - - template - int SynchronousMachine::evaluateJacobian() - { - return 0; - } - - template - int SynchronousMachine::evaluateIntegrand() + namespace PowerElectronics { - return 0; - } - - template - int SynchronousMachine::initializeAdjoint() - { - return 0; - } - - template - int SynchronousMachine::evaluateAdjointResidual() - { - return 0; - } - - template - int SynchronousMachine::evaluateAdjointIntegrand() - { - return 0; - } - - template - CircuitComponent* SynchronousMachine::clone() const - { - return new SynchronousMachine(*this); - } - - // Available template instantiations - template class SynchronousMachine; - template class SynchronousMachine; - template class SynchronousMachine; - template class SynchronousMachine; - + /*! + * @brief Constructor for a constant SynchronousMachine model + * + * Calls default ModelEvaluatorImpl constructor. + * @todo This model's equations are not finished + * @todo needs to be tested for correctness + * + * @tparam ScalarT - floating point type for the model + * @tparam IdxT - integer index type for the model + * + * @param[in] id - unique identifier for the component + * @param[in] Lls - stator leakage inductance + * @param[in] Llkq - tuple of damper leakage reactances + * @param[in] Llfd - field leakage reactance + * @param[in] Llkd - damper leakage reactance + * @param[in] Lmq - quadrature axis magnetizing reactance + * @param[in] Lmd - direct axis magnetizing reactance + * @param[in] Rs - stator resistance + * @param[in] Rkq - tuple of damper resistances + * @param[in] Rfd - field resistance + * @param[in] Rkd - damper resistance + * @param[in] RJ - rotor moment of inertia + * @param[in] P - number of poles + * @param[in] mub - rated frequency + */ + template + SynchronousMachine::SynchronousMachine(IdxT id, RealT Lls, std::tuple Llkq, RealT Llfd, RealT Llkd, RealT Lmq, RealT Lmd, RealT Rs, std::tuple Rkq, RealT Rfd, RealT Rkd, RealT RJ, RealT P, RealT mub) + : Lls_(Lls), + Llkq_(Llkq), + Llfd_(Llfd), + Llkd_(Llkd), + Lmq_(Lmq), + Lmd_(Lmd), + Rs_(Rs), + Rkq_(Rkq), + Rfd_(Rfd), + Rkd_(Rkd), + RJ_(RJ), + P_(P), + mub_(mub) + { + size_ = 13; + n_intern_ = 6; + n_extern_ = 7; + extern_indices_ = {0, 1, 2, 3, 4}; + idc_ = id; + } + + template + SynchronousMachine::~SynchronousMachine() + { + } + + /** + * Initialization of the grid model + */ + template + int SynchronousMachine::initialize() + { + return 0; + } + + /* + * \brief Identify differential variables + */ + template + int SynchronousMachine::tagDifferentiable() + { + // All variables are differentials + std::fill(tag_.begin(), tag_.end(), true); + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int SynchronousMachine::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Compute the resisdual of the component. + * + * @todo not finished + */ + template + int SynchronousMachine::evaluateInternalResidual() + { + ScalarT rkq1 = static_cast(std::get<0>(Rkq_)); + [[maybe_unused]] ScalarT rkq2 = static_cast(std::get<1>(Rkq_)); + ScalarT llkq1 = static_cast(std::get<0>(Llkq_)); + [[maybe_unused]] ScalarT llkq2 = static_cast(std::get<1>(Llkq_)); + + static constexpr auto pi = std::numbers::pi_v; + + ScalarT cos1 = std::cos((P_ / 2.0) * y_int_[0]); + ScalarT sin1 = std::sin((P_ / 2.0) * y_int_[0]); + ScalarT cos23m = std::cos((P_ / 2.0) * y_int_[0] - (2.0 / 3.0) * pi); + ScalarT sin23m = std::sin((P_ / 2.0) * y_int_[0] - (2.0 / 3.0) * pi); + ScalarT cos23p = std::cos((P_ / 2.0) * y_int_[0] + (2.0 / 3.0) * pi); + ScalarT sin23p = std::sin((P_ / 2.0) * y_int_[0] + (2.0 / 3.0) * pi); + + f_int_[0] = (-2.0 / 3.0) * (*y_ext_[0] * cos1 + *y_ext_[1] * cos23m + *y_ext_[2] * cos23p) + Rs_ * y_int_[1] + (Lls_ + Lmq_) * yp_int_[1] + Lmq_ * yp_int_[4] + Lmq_ * yp_int_[5] + *y_ext_[4] * (P_ / 2.0) * ((Lls_ + Lmd_) * y_int_[2] + Lmd_ * y_int_[6] + Lmd_ * y_int_[7]); + f_int_[1] = (-2.0 / 3.0) * (*y_ext_[0] * sin1 - *y_ext_[1] * sin23m - *y_ext_[2] * sin23p) + Rs_ * y_int_[2] + (Lls_ + Lmd_) * yp_int_[2] + Lmd_ * yp_int_[6] + Lmd_ * yp_int_[7] - *y_ext_[4] * (P_ / 2.0) * ((Lls_ + Lmq_) * y_int_[1] + Lmq_ * y_int_[4] + Lmq_ * y_int_[5]); + f_int_[2] = (-1.0 / 3.0) * (*y_ext_[0] + *y_ext_[1] + *y_ext_[2]) + Rs_ * y_int_[3] + Lls_ * yp_int_[3]; + f_int_[3] = rkq1 * y_int_[4] + (llkq1 + Lmq_) * yp_int_[4] + Lmq_ * yp_int_[1] + Lmq_ * yp_int_[5]; + f_int_[4] = rkq1 * y_int_[4] + (llkq1 + Lmq_) * yp_int_[4] + Lmq_ * yp_int_[1] + Lmq_ * yp_int_[5]; + return 0; + } + + template + int SynchronousMachine::evaluateExternalResidual() + { + [[maybe_unused]] ScalarT rkq2 = static_cast(std::get<1>(Rkq_)); + [[maybe_unused]] ScalarT llkq2 = static_cast(std::get<1>(Llkq_)); + + static constexpr auto pi = std::numbers::pi_v; + + ScalarT cos1 = std::cos((P_ / 2.0) * y_int_[0]); + ScalarT sin1 = std::sin((P_ / 2.0) * y_int_[0]); + ScalarT cos23m = std::cos((P_ / 2.0) * y_int_[0] - (2.0 / 3.0) * pi); + ScalarT sin23m = std::sin((P_ / 2.0) * y_int_[0] - (2.0 / 3.0) * pi); + ScalarT cos23p = std::cos((P_ / 2.0) * y_int_[0] + (2.0 / 3.0) * pi); + ScalarT sin23p = std::sin((P_ / 2.0) * y_int_[0] + (2.0 / 3.0) * pi); + + *f_ext_[0] += y_int_[1] * cos1 + y_int_[2] * sin1 + y_int_[3]; + *f_ext_[1] += y_int_[1] * cos23m + y_int_[2] * sin23m + y_int_[3]; + *f_ext_[2] += y_int_[1] * cos23p + y_int_[2] * sin23p + y_int_[3]; + *f_ext_[3] += RJ_ * *yp_ext_[4] - (3.0 / 4.0) * P_ * (Lmd_ * y_int_[1] * (y_int_[2] + y_int_[6] + y_int_[7]) - Lmq_ * y_int_[2] * (y_int_[1] + y_int_[4] + *y_ext_[0])); + *f_ext_[4] += yp_int_[0] - *y_ext_[4]; + return 0; + } + + template + int SynchronousMachine::evaluateJacobian() + { + return 0; + } + + template + int SynchronousMachine::evaluateIntegrand() + { + return 0; + } + + template + int SynchronousMachine::initializeAdjoint() + { + return 0; + } + + template + int SynchronousMachine::evaluateAdjointResidual() + { + return 0; + } + + template + int SynchronousMachine::evaluateAdjointIntegrand() + { + return 0; + } + + template + Component* SynchronousMachine::clone() const + { + return new SynchronousMachine(*this); + } + + // Available template instantiations + template class SynchronousMachine; + template class SynchronousMachine; + template class SynchronousMachine; + template class SynchronousMachine; + + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp b/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp index 55356f191..39c21072e 100644 --- a/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp +++ b/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp @@ -1,84 +1,80 @@ - #pragma once #include -#include - -namespace GridKit -{ - template - class BaseBus; -} +#include namespace GridKit { - /*! - * @brief Declaration of a SynchronousMachine class. - * - */ - template - class SynchronousMachine : public CircuitComponent + namespace PowerElectronics { - using RealT = typename CircuitComponent::RealT; + /*! + * @brief Declaration of a SynchronousMachine class. + * + */ + template + class SynchronousMachine : public Component + { + using RealT = typename Component::RealT; - using CircuitComponent::size_; - using CircuitComponent::nnz_; - using CircuitComponent::time_; - using CircuitComponent::alpha_; - using CircuitComponent::y_ext_; - using CircuitComponent::y_int_; - using CircuitComponent::yp_ext_; - using CircuitComponent::yp_int_; - using CircuitComponent::abs_tol_; - using CircuitComponent::tag_; - using CircuitComponent::f_ext_; - using CircuitComponent::f_int_; - using CircuitComponent::g_; - using CircuitComponent::yB_; - using CircuitComponent::ypB_; - using CircuitComponent::fB_; - using CircuitComponent::gB_; - using CircuitComponent::param_; - using CircuitComponent::idc_; + using Component::size_; + using Component::nnz_; + using Component::time_; + using Component::alpha_; + using Component::y_ext_; + using Component::y_int_; + using Component::yp_ext_; + using Component::yp_int_; + using Component::abs_tol_; + using Component::tag_; + using Component::f_ext_; + using Component::f_int_; + using Component::g_; + using Component::yB_; + using Component::ypB_; + using Component::fB_; + using Component::gB_; + using Component::param_; + using Component::idc_; - using CircuitComponent::extern_indices_; - using CircuitComponent::n_extern_; - using CircuitComponent::n_intern_; + using Component::extern_indices_; + using Component::n_extern_; + using Component::n_intern_; - public: - SynchronousMachine(IdxT id, RealT Lls, std::tuple Llkq, RealT Llfd, RealT Llkd, RealT Lmq, RealT Lmd, RealT Rs, std::tuple Rkq, RealT Rfd, RealT Rkd, RealT RJ, RealT P, RealT mub); - virtual ~SynchronousMachine(); + public: + SynchronousMachine(IdxT id, RealT Lls, std::tuple Llkq, RealT Llfd, RealT Llkd, RealT Lmq, RealT Lmd, RealT Rs, std::tuple Rkq, RealT Rfd, RealT Rkd, RealT RJ, RealT P, RealT mub); + virtual ~SynchronousMachine(); - int initialize(); - int tagDifferentiable(); - int setAbsoluteTolerance(RealT); - int evaluateInternalResidual() final; - int evaluateExternalResidual() final; - int evaluateJacobian(); - int evaluateIntegrand(); + int initialize(); + int tagDifferentiable(); + int setAbsoluteTolerance(RealT); + int evaluateInternalResidual() final; + int evaluateExternalResidual() final; + int evaluateJacobian(); + int evaluateIntegrand(); - int initializeAdjoint(); - int evaluateAdjointResidual(); - // int evaluateAdjointJacobian(); - int evaluateAdjointIntegrand(); + int initializeAdjoint(); + int evaluateAdjointResidual(); + // int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); - CircuitComponent* clone() const; + Component* clone() const; - private: - RealT Lls_; - std::tuple Llkq_; - RealT Llfd_; - RealT Llkd_; - RealT Lmq_; - RealT Lmd_; - RealT Rs_; - std::tuple Rkq_; - RealT Rfd_; - RealT Rkd_; - RealT RJ_; - RealT P_; - RealT mub_; - }; + private: + RealT Lls_; + std::tuple Llkq_; + RealT Llfd_; + RealT Llkd_; + RealT Lmq_; + RealT Lmd_; + RealT Rs_; + std::tuple Rkq_; + RealT Rfd_; + RealT Rkd_; + RealT RJ_; + RealT P_; + RealT mub_; + }; + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/SystemModel.hpp b/GridKit/Model/PowerElectronics/SystemModel.hpp new file mode 100644 index 000000000..c3a997b9c --- /dev/null +++ b/GridKit/Model/PowerElectronics/SystemModel.hpp @@ -0,0 +1,532 @@ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PowerElectronics + { + template + class SystemModel : public Component + { + using RealT = typename Component::RealT; + using CsrMatrixT = typename Component::CsrMatrixT; + using component_type = Component; + using node_type = NodeBase; + + using Component::abs_tol_; + using Component::allocated_; + using Component::allocateVectors; + using Component::alpha_; + using Component::f_ext_; + using Component::f_int_; + using Component::n_extern_; + using Component::n_intern_; + using Component::nnz_; + using Component::size_; + using Component::tag_; + using Component::time_; + using Component::y_ext_; + using Component::y_int_; + using Component::yp_ext_; + using Component::yp_int_; + + public: + /** + * @brief Default constructor for the system model + * + * @post System model parameters set as default + */ + SystemModel() + { + // By default don't use the jacobian + use_jac_ = false; + } + + /** + * @brief Constructor for the system model + * + * @param[in] use_jac Boolean to choose if to use jacobian + * + * @post System model parameters set as input + */ + SystemModel(bool use_jac = false) + { + // Can choose if to use jacobian + use_jac_ = use_jac; + } + + /** + * @brief Destructor for the system model + * + * @pre System components are allocated + * + * @post System components are deallocated + * + */ + virtual ~SystemModel() + { + for (auto comp : this->components_) + { + delete comp; + } + delete csr_jac_; + delete[] map_to_csr_; + } + + /** + * @brief Will check if each component has jacobian avalible. If one doesn't have it, return false. + * + * + * @return true if all components have jacobian + * @return false otherwise + */ + bool hasJacobian() final + { + if (!this->use_jac_) + return false; + + for (const auto& component : components_) + { + if (!component->hasJacobian()) + { + return false; + } + } + return true; + } + + /** + * @brief Allocate system vectors and construct the system CSR Jacobian + * + * @post System model vectors allocated with the computed total number of unknowns + * @post CSR Jacobian sparsity pattern is computed + * @post COO->CSR mapping is computed + * @post Every component's \ref Component::y_int_, \ref Component::yp_int_, and \ref Component::f_int_ pointers + * are set to their appropriate offsets in the system vector, allowing them to directly access their internal variables, derivatives, + * and residuals. + * + * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable + */ + int allocate() final + { + size_t component_internal_size = 0; + for (component_type* comp : components_) + { + component_internal_size += comp->getInternalSize(); + } + + size_t node_internal_size = 0; + for (node_type* node : nodes_) + { + node_internal_size += node->getInternalSize(); + } + + n_intern_ = component_internal_size + node_internal_size; + n_extern_ = 0; + size_ = n_intern_ + n_extern_; + + // Allocation always rebuilds the system Jacobian and its COO-to-CSR map. + delete csr_jac_; + csr_jac_ = nullptr; + + delete[] map_to_csr_; + map_to_csr_ = nullptr; + + if (!allocated_) + { + allocateVectors(static_cast(size_), true); + // Component and node offsets can change when topology is modified. + abs_tol_.setToZero(memory::HOST); + } + + tag_.resize(size_); + + { // Start node internal indexing after all component internals for proper KLU ordering + size_t node_internal_idx = component_internal_size; + for (node_type* node : nodes_) + { + node->allocate(); + + for (size_t i = 0; i < node->getInternalSize(); i++) + { + ExternalConnection node_connection{ + .y_ = y_int_ + node_internal_idx, + .yp_ = yp_int_ + node_internal_idx, + .f_ = f_int_ + node_internal_idx, + .idx_ = static_cast(node_internal_idx)}; + + node->setExternalConnectionNodes(i, node_connection); + node_internal_idx++; + } + } + } + + { + // The offset for each component's internal variables in the system vector. + // They start at 0, and are stacked on top of each other. + size_t component_internal_idx = 0; + for (component_type* comp : components_) + { + comp->allocate(); + + // Update component internal pointers to their correct offsets + comp->setInternalPointer(&y_int_[component_internal_idx]); + comp->setInternalDerivativePointer(&yp_int_[component_internal_idx]); + comp->setInternalResidualPointer(&f_int_[component_internal_idx]); + + const auto& external_indices = comp->getExternIndices(); + for (IdxT i = 0; i < comp->size(); i++) + { + if (!external_indices.contains(i)) + { + comp->setInternalConnectionNodes(i, component_internal_idx); + component_internal_idx++; + } + } + } + } + + // Evaluate component Jacobians to get sparsity + for (component_type* component : components_) + { + component->evaluateJacobian(); + } + + // Count the number of non-zeros + IdxT nnz_dup = 0; + for (const component_type* component : components_) + { + const IdxT* r = component->jacobianCooRows(); + const IdxT* c = component->jacobianCooCols(); + IdxT nnz = component->nnz(); + + for (IdxT i = 0; i < nnz; ++i) + { + if (component->getNodeConnection(r[i]) != neg1_ && component->getNodeConnection(c[i]) != neg1_) + { + ++nnz_dup; + } + } + } + + // Allocate COO triplet arrays (we own these until we hand off to CsrMatrix) + IdxT* rows_dup = new IdxT[nnz_dup]; + IdxT* cols_dup = new IdxT[nnz_dup]; + RealT* vals_dup = new RealT[nnz_dup]; + + IdxT counter = 0; + for (const component_type* component : components_) + { + const IdxT* r = component->jacobianCooRows(); + const IdxT* c = component->jacobianCooCols(); + const RealT* v = component->jacobianCooValues(); + IdxT nnz = component->nnz(); + + for (IdxT i = 0; i < nnz; ++i) + { + if (component->getNodeConnection(r[i]) != neg1_ && component->getNodeConnection(c[i]) != neg1_) + { + rows_dup[counter] = component->getNodeConnection(r[i]); + cols_dup[counter] = component->getNodeConnection(c[i]); + vals_dup[counter] = v[i]; + counter++; + } + } + } + + // Build the system COO Jacobian + LinearAlgebra::CooMatrix jac(size_, size_, nnz_dup, &rows_dup, &cols_dup, &vals_dup); + + // Populate CSR data with sort and deduplicate + IdxT* row_ptrs = jac.getCsrRowData(); + + // Deduplicated nnz + nnz_ = jac.getNnz(); + + // Allocate cols/vals with deduplicated nnz + IdxT* cols = new IdxT[nnz_]; + RealT* vals = new RealT[nnz_]; + + std::copy(jac.getColData(), jac.getColData() + nnz_, cols); + std::copy(jac.getValues(), jac.getValues() + nnz_, vals); + + // Create the CSR Jacobian + csr_jac_ = new CsrMatrixT(size_, size_, nnz_, &row_ptrs, &cols, &vals); + + const IdxT* map_to_sorted = jac.getMapToSorted(); + const IdxT* map_to_dedup = jac.getMapToDeduplicated(); + + // Build a mappping from original COO index to CSR index + map_to_csr_ = new IdxT[nnz_dup]; + for (IdxT i = 0; i < nnz_dup; ++i) + { + map_to_csr_[map_to_sorted[i]] = map_to_dedup[i]; + } + + allocated_ = true; + return 0; + } + + /** + * @brief Set intial y and y' of each component + * + * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable + */ + int initialize() final + { + // Initialize components + for (const auto& component : components_) + { + component->initialize(); + } + + return Component::initialize(); + } + + /** + * @brief Tags all system variables as differentiable, based on what the + * components that own those variables tag them as. + * + * Starts by asking all components to tag their differentiables. This implementation + * assumes all node variables are algebraic, and will not ask nodes to tag their differentiables. + * Sets all variables to algebraic (`false`) to start, then loops over all component internal variables. + * Re-creates the same internal variable to system variables mapping as in \ref allocate() - all + * internal variables from the same component are stored contiguously in a block, and blocks are + * stored contiguously in the same order as \ref components_, with node variables at the end. + * Each internal variable's tag in the system is set to its tag in the component. + */ + int tagDifferentiable() final + { + // Ask all component to tag their differentiables + for (size_t i = 0; i < components_.size(); i++) + { + component_type* component = components_[i]; + + // Bubble up errors if necessary + if (int err = component->tagDifferentiable()) + { + return err; + } + } + + // Fill tags with a default value (false) for node variables. Assumed to be algebraic here. + std::fill(tag_.begin(), tag_.end(), false); + + // Copy tags for internal variables from their components - going in the order as described above + size_t idx = 0; + for (component_type* comp : components_) + { + const auto& external_indices = comp->getExternIndices(); + + // Loop over all component variables - including externals + for (IdxT i = 0; i < comp->size(); i++) + { + // Discard externals + if (!external_indices.contains(i)) + { + tag_[idx] = comp->tag()[i]; + + // Ensures internal variables are contiguous, and in the same order as the component + idx++; + } + } + } + + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + int setAbsoluteTolerance(RealT rel_tol) final + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Evaluate Residuals at each component then collect them + * + * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable + */ + int evaluateInternalResidual() final + { + for (IdxT i = 0; i < size_; i++) + { + f_int_[i] = 0.0; + } + + // Update system residual vector + + // Evaluate component internal residuals - this is embarassingly parallel + for (component_type* component : components_) + { + if (int err_code = component->evaluateInternalResidual()) + return err_code; + } + + for (component_type* component : components_) + { + if (int err_code = component->evaluateExternalResidual()) + return err_code; + } + + return 0; + } + + /** + * @todo implement this for nested systems + */ + int evaluateExternalResidual() final + { + return 0; + } + + /** + * @brief Creates the system Jacobian representing \f$\alpha dF/dy' + dF/dy\f$ + * + * Updates the CSR Jacobian values using the per-component mappings + * computed during allocate(). + * + * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable + */ + int evaluateJacobian() final + { + // Zero out values + RealT* vals = csr_jac_->getValues(); + for (IdxT i = 0; i < csr_jac_->getNnz(); ++i) + { + vals[i] = 0.0; + } + + // Update CSR values from component Jacobians + IdxT counter = 0; + for (const auto& component : components_) + { + component->evaluateJacobian(); + + const IdxT* r = component->jacobianCooRows(); + const IdxT* c = component->jacobianCooCols(); + const RealT* v = component->jacobianCooValues(); + IdxT nnz = component->nnz(); + + for (IdxT i = 0; i < nnz; ++i) + { + if (component->getNodeConnection(r[i]) != neg1_ && component->getNodeConnection(c[i]) != neg1_) + { + vals[map_to_csr_[counter]] += v[i]; + ++counter; + } + } + } + + jac_call_count_++; + return 0; + } + + /** + * @brief Evaluate integrands for the system quadratures. + */ + int evaluateIntegrand() final + { + return 0; + } + + /** + * @brief Initialize system adjoint. + * + * Updates variables and optimization parameters, then initializes + * adjoints locally and copies them to the system adjoint vector. + */ + int initializeAdjoint() final + { + return 0; + } + + /** + * @brief Compute adjoint residual for the system model. + * + * + */ + int evaluateAdjointResidual() final + { + return 0; + } + + /** + * @brief Evaluate adjoint integrand for the system model. + * + * + */ + int evaluateAdjointIntegrand() final + { + return 0; + } + + /** + * @brief Distribute time and time scaling for each component + * + * @param t + * @param a + */ + void updateTime(RealT t, RealT a) final + { + for (const auto& component : components_) + { + component->updateTime(t, a); + } + time_ = t; + alpha_ = a; + } + + CsrMatrixT* getCsrJacobian() const override + { + return csr_jac_; + } + + void addComponent(component_type* component) + { + components_.push_back(component); + allocated_ = false; + } + + void addNode(node_type* node) + { + nodes_.push_back(node); + allocated_ = false; + } + + private: + static constexpr IdxT neg1_ = INVALID_INDEX; + + std::vector components_; + std::vector nodes_; + + IdxT* map_to_csr_{nullptr}; + CsrMatrixT* csr_jac_{nullptr}; + + int jac_call_count_{0}; + bool use_jac_; + + }; // class SystemModel + } // namespace PowerElectronics +} // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp deleted file mode 100644 index 51c8dcee1..000000000 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ /dev/null @@ -1,531 +0,0 @@ - - -#pragma once - -#include -#include -#include - -#include -#include -#include -#include -#include - -namespace GridKit -{ - template - class PowerElectronicsModel : public CircuitComponent - { - using RealT = typename CircuitComponent::RealT; - using CsrMatrixT = typename CircuitComponent::CsrMatrixT; - using component_type = CircuitComponent; - using node_type = PowerElectronics::NodeBase; - - using CircuitComponent::abs_tol_; - using CircuitComponent::allocated_; - using CircuitComponent::allocateVectors; - using CircuitComponent::alpha_; - using CircuitComponent::f_ext_; - using CircuitComponent::f_int_; - using CircuitComponent::n_extern_; - using CircuitComponent::n_intern_; - using CircuitComponent::nnz_; - using CircuitComponent::size_; - using CircuitComponent::tag_; - using CircuitComponent::time_; - using CircuitComponent::y_ext_; - using CircuitComponent::y_int_; - using CircuitComponent::yp_ext_; - using CircuitComponent::yp_int_; - - public: - /** - * @brief Default constructor for the system model - * - * @post System model parameters set as default - */ - PowerElectronicsModel() - { - // By default don't use the jacobian - use_jac_ = false; - } - - /** - * @brief Constructor for the system model - * - * @param[in] use_jac Boolean to choose if to use jacobian - * - * @post System model parameters set as input - */ - PowerElectronicsModel(bool use_jac = false) - { - // Can choose if to use jacobian - use_jac_ = use_jac; - } - - /** - * @brief Destructor for the system model - * - * @pre System components are allocated - * - * @post System components are deallocated - * - */ - virtual ~PowerElectronicsModel() - { - for (auto comp : this->components_) - { - delete comp; - } - delete csr_jac_; - delete[] map_to_csr_; - } - - /** - * @brief Will check if each component has jacobian avalible. If one doesn't have it, return false. - * - * - * @return true if all components have jacobian - * @return false otherwise - */ - bool hasJacobian() final - { - if (!this->use_jac_) - return false; - - for (const auto& component : components_) - { - if (!component->hasJacobian()) - { - return false; - } - } - return true; - } - - /** - * @brief Allocate system vectors and construct the system CSR Jacobian - * - * @post System model vectors allocated with the computed total number of unknowns - * @post CSR Jacobian sparsity pattern is computed - * @post COO->CSR mapping is computed - * @post Every component's \ref CircuitComponent::y_int_, \ref CircuitComponent::yp_int_, and \ref CircuitComponent::f_int_ pointers - * are set to their appropriate offsets in the system vector, allowing them to directly access their internal variables, derivatives, - * and residuals. - * - * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable - */ - int allocate() final - { - size_t component_internal_size = 0; - for (component_type* comp : components_) - { - component_internal_size += comp->getInternalSize(); - } - - size_t node_internal_size = 0; - for (node_type* node : nodes_) - { - node_internal_size += node->getInternalSize(); - } - - n_intern_ = component_internal_size + node_internal_size; - n_extern_ = 0; - size_ = n_intern_ + n_extern_; - - // Allocation always rebuilds the system Jacobian and its COO-to-CSR map. - delete csr_jac_; - csr_jac_ = nullptr; - - delete[] map_to_csr_; - map_to_csr_ = nullptr; - - if (!allocated_) - { - allocateVectors(static_cast(size_), true); - // Component and node offsets can change when topology is modified. - abs_tol_.setToZero(memory::HOST); - } - - tag_.resize(size_); - - { // Start node internal indexing after all component internals for proper KLU ordering - size_t node_internal_idx = component_internal_size; - for (node_type* node : nodes_) - { - node->allocate(); - - for (size_t i = 0; i < node->getInternalSize(); i++) - { - ExternalConnection node_connection{ - .y_ = y_int_ + node_internal_idx, - .yp_ = yp_int_ + node_internal_idx, - .f_ = f_int_ + node_internal_idx, - .idx_ = static_cast(node_internal_idx)}; - - node->setExternalConnectionNodes(i, node_connection); - node_internal_idx++; - } - } - } - - { - // The offset for each component's internal variables in the system vector. - // They start at 0, and are stacked on top of each other. - size_t component_internal_idx = 0; - for (component_type* comp : components_) - { - comp->allocate(); - - // Update component internal pointers to their correct offsets - comp->setInternalPointer(&y_int_[component_internal_idx]); - comp->setInternalDerivativePointer(&yp_int_[component_internal_idx]); - comp->setInternalResidualPointer(&f_int_[component_internal_idx]); - - const auto& external_indices = comp->getExternIndices(); - for (IdxT i = 0; i < comp->size(); i++) - { - if (!external_indices.contains(i)) - { - comp->setInternalConnectionNodes(i, component_internal_idx); - component_internal_idx++; - } - } - } - } - - // Evaluate component Jacobians to get sparsity - for (component_type* component : components_) - { - component->evaluateJacobian(); - } - - // Count the number of non-zeros - IdxT nnz_dup = 0; - for (const component_type* component : components_) - { - const IdxT* r = component->jacobianCooRows(); - const IdxT* c = component->jacobianCooCols(); - IdxT nnz = component->nnz(); - - for (IdxT i = 0; i < nnz; ++i) - { - if (component->getNodeConnection(r[i]) != neg1_ && component->getNodeConnection(c[i]) != neg1_) - { - ++nnz_dup; - } - } - } - - // Allocate COO triplet arrays (we own these until we hand off to CsrMatrix) - IdxT* rows_dup = new IdxT[nnz_dup]; - IdxT* cols_dup = new IdxT[nnz_dup]; - RealT* vals_dup = new RealT[nnz_dup]; - - IdxT counter = 0; - for (const component_type* component : components_) - { - const IdxT* r = component->jacobianCooRows(); - const IdxT* c = component->jacobianCooCols(); - const RealT* v = component->jacobianCooValues(); - IdxT nnz = component->nnz(); - - for (IdxT i = 0; i < nnz; ++i) - { - if (component->getNodeConnection(r[i]) != neg1_ && component->getNodeConnection(c[i]) != neg1_) - { - rows_dup[counter] = component->getNodeConnection(r[i]); - cols_dup[counter] = component->getNodeConnection(c[i]); - vals_dup[counter] = v[i]; - counter++; - } - } - } - - // Build the system COO Jacobian - LinearAlgebra::CooMatrix jac(size_, size_, nnz_dup, &rows_dup, &cols_dup, &vals_dup); - - // Populate CSR data with sort and deduplicate - IdxT* row_ptrs = jac.getCsrRowData(); - - // Deduplicated nnz - nnz_ = jac.getNnz(); - - // Allocate cols/vals with deduplicated nnz - IdxT* cols = new IdxT[nnz_]; - RealT* vals = new RealT[nnz_]; - - std::copy(jac.getColData(), jac.getColData() + nnz_, cols); - std::copy(jac.getValues(), jac.getValues() + nnz_, vals); - - // Create the CSR Jacobian - csr_jac_ = new CsrMatrixT(size_, size_, nnz_, &row_ptrs, &cols, &vals); - - const IdxT* map_to_sorted = jac.getMapToSorted(); - const IdxT* map_to_dedup = jac.getMapToDeduplicated(); - - // Build a mappping from original COO index to CSR index - map_to_csr_ = new IdxT[nnz_dup]; - for (IdxT i = 0; i < nnz_dup; ++i) - { - map_to_csr_[map_to_sorted[i]] = map_to_dedup[i]; - } - - allocated_ = true; - return 0; - } - - /** - * @brief Set intial y and y' of each component - * - * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable - */ - int initialize() final - { - // Initialize components - for (const auto& component : components_) - { - component->initialize(); - } - - return CircuitComponent::initialize(); - } - - /** - * @brief Tags all system variables as differentiable, based on what the - * components that own those variables tag them as. - * - * Starts by asking all components to tag their differentiables. This implementation - * assumes all node variables are algebraic, and will not ask nodes to tag their differentiables. - * Sets all variables to algebraic (`false`) to start, then loops over all component internal variables. - * Re-creates the same internal variable to system variables mapping as in \ref allocate() - all - * internal variables from the same component are stored contiguously in a block, and blocks are - * stored contiguously in the same order as \ref components_, with node variables at the end. - * Each internal variable's tag in the system is set to its tag in the component. - */ - int tagDifferentiable() final - { - // Ask all component to tag their differentiables - for (size_t i = 0; i < components_.size(); i++) - { - component_type* component = components_[i]; - - // Bubble up errors if necessary - if (int err = component->tagDifferentiable()) - { - return err; - } - } - - // Fill tags with a default value (false) for node variables. Assumed to be algebraic here. - std::fill(tag_.begin(), tag_.end(), false); - - // Copy tags for internal variables from their components - going in the order as described above - size_t idx = 0; - for (component_type* comp : components_) - { - const auto& external_indices = comp->getExternIndices(); - - // Loop over all component variables - including externals - for (IdxT i = 0; i < comp->size(); i++) - { - // Discard externals - if (!external_indices.contains(i)) - { - tag_[idx] = comp->tag()[i]; - - // Ensures internal variables are contiguous, and in the same order as the component - idx++; - } - } - } - - return 0; - } - - /** - * @brief Compute the absolute tolerance for each variable in the model - * - * @param rel_tol The relative tolerance which can be used to pick the - * absolute tolerance. - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return int 0 if successful, non-zero otherwise. - * - * This represents a "noise" level close to zero for which pure relative - * error cannot be used. - */ - int setAbsoluteTolerance(RealT rel_tol) final - { - abs_tol_.setToConst(static_cast(rel_tol)); - return 0; - } - - /** - * @brief Evaluate Residuals at each component then collect them - * - * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable - */ - int evaluateInternalResidual() final - { - for (IdxT i = 0; i < size_; i++) - { - f_int_[i] = 0.0; - } - - // Update system residual vector - - // Evaluate component internal residuals - this is embarassingly parallel - for (component_type* component : components_) - { - if (int err_code = component->evaluateInternalResidual()) - return err_code; - } - - for (component_type* component : components_) - { - if (int err_code = component->evaluateExternalResidual()) - return err_code; - } - - return 0; - } - - /** - * @todo implement this for nested systems - */ - int evaluateExternalResidual() final - { - return 0; - } - - /** - * @brief Creates the system Jacobian representing \f$\alpha dF/dy' + dF/dy\f$ - * - * Updates the CSR Jacobian values using the per-component mappings - * computed during allocate(). - * - * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable - */ - int evaluateJacobian() final - { - // Zero out values - RealT* vals = csr_jac_->getValues(); - for (IdxT i = 0; i < csr_jac_->getNnz(); ++i) - { - vals[i] = 0.0; - } - - // Update CSR values from component Jacobians - IdxT counter = 0; - for (const auto& component : components_) - { - component->evaluateJacobian(); - - const IdxT* r = component->jacobianCooRows(); - const IdxT* c = component->jacobianCooCols(); - const RealT* v = component->jacobianCooValues(); - IdxT nnz = component->nnz(); - - for (IdxT i = 0; i < nnz; ++i) - { - if (component->getNodeConnection(r[i]) != neg1_ && component->getNodeConnection(c[i]) != neg1_) - { - vals[map_to_csr_[counter]] += v[i]; - ++counter; - } - } - } - - jac_call_count_++; - return 0; - } - - /** - * @brief Evaluate integrands for the system quadratures. - */ - int evaluateIntegrand() final - { - return 0; - } - - /** - * @brief Initialize system adjoint. - * - * Updates variables and optimization parameters, then initializes - * adjoints locally and copies them to the system adjoint vector. - */ - int initializeAdjoint() final - { - return 0; - } - - /** - * @brief Compute adjoint residual for the system model. - * - * - */ - int evaluateAdjointResidual() final - { - return 0; - } - - /** - * @brief Evaluate adjoint integrand for the system model. - * - * - */ - int evaluateAdjointIntegrand() final - { - return 0; - } - - /** - * @brief Distribute time and time scaling for each component - * - * @param t - * @param a - */ - void updateTime(RealT t, RealT a) final - { - for (const auto& component : components_) - { - component->updateTime(t, a); - } - time_ = t; - alpha_ = a; - } - - CsrMatrixT* getCsrJacobian() const override - { - return csr_jac_; - } - - void addComponent(component_type* component) - { - components_.push_back(component); - allocated_ = false; - } - - void addNode(node_type* node) - { - nodes_.push_back(node); - allocated_ = false; - } - - private: - static constexpr IdxT neg1_ = INVALID_INDEX; - - std::vector components_; - std::vector nodes_; - - IdxT* map_to_csr_{nullptr}; - CsrMatrixT* csr_jac_{nullptr}; - - int jac_call_count_{0}; - bool use_jac_; - - }; // class PowerElectronicsModel - -} // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp b/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp index f3be694d4..cd120c2de 100644 --- a/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp +++ b/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp @@ -7,214 +7,215 @@ namespace GridKit { - - /*! - * @brief Constructor for a TransmissionLine model - * - * Calls default ModelEvaluatorImpl constructor. - * - * This is the Medium distance form with the use of the admittance matrix. - * Since the line is of medium length then there is no real part for shunt admittance - * @todo needs to used in a model - * @todo test for correctness - */ - - template - TransmissionLine::TransmissionLine(IdxT id, RealT R, RealT X, RealT B) - : R_(R), - X_(X), - B_(B) - { - // internals [Iret1, Iimt1, Iret2, Iimt2] - // externals [Vre11, Vim11, Vre12, Vim12, Vre21, Vim21, Vre22, Vim22] - size_ = 12; - n_intern_ = 4; - n_extern_ = 8; - extern_indices_ = {0, 1, 2, 3, 4, 5, 6, 7}; - idc_ = id; - nnz_ = 44; - - RealT magImpendence = 1.0 / (R_ * R_ + X_ * X_); - YReMat_ = magImpendence * R_; - YImMatOff_ = magImpendence * X_; - YImMatDi_ = B_ / (2.0) - YImMatOff_; - } - - template - TransmissionLine::~TransmissionLine() - { - } - - /** - * Initialization of the grid model - */ - template - int TransmissionLine::initialize() - { - return 0; - } - - /* - * \brief Identify differential variables - */ - template - int TransmissionLine::tagDifferentiable() - { - // All variables are algebraics - std::fill(tag_.begin(), tag_.end(), false); - return 0; - } - - /** - * @brief Compute the absolute tolerance for each variable in the model - * - * @param rel_tol The relative tolerance which can be used to pick the - * absolute tolerance. - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return int 0 if successful, non-zero otherwise. - * - * This represents a "noise" level close to zero for which pure relative - * error cannot be used. - */ - template - int TransmissionLine::setAbsoluteTolerance(RealT rel_tol) - { - abs_tol_.setToConst(static_cast(rel_tol)); - return 0; - } - - /** - * @brief Evaluate residual of transmission line - * - * The complex admittance matrix is: - * [[ Y/2 + 1/Z, -1/Z]; - * [ -1/Z, Y/2 + 1/Z]] = - * [[R/|Z|, -R/|Z|]; - * [-R/|Z|, R/|Z|]] + - * i [[B/2 - X/|Z|, X/|Z|]; - * [X/|Z|, B/2 - X/|Z|]] - * = Dre + i Dim - * - * Then - * Ire = Dre Vre - Dim Vim - * Iim = Dre Vim + Dim Vre - * - * To express this for Modified Nodal Analysis the Voltages of the admittance matrix are put into voltage drops - */ - template - int TransmissionLine::evaluateInternalResidual() - { - // Voltage drop accross terminals - ScalarT V1re = *y_ext_[0] - *y_ext_[4]; - ScalarT V1im = *y_ext_[1] - *y_ext_[5]; - ScalarT V2re = *y_ext_[2] - *y_ext_[6]; - ScalarT V2im = *y_ext_[3] - *y_ext_[7]; - - // Internal variables - // row 1 - f_int_[0] = YReMat_ * (V1re - V2re) - (YImMatDi_ * V1im + YImMatOff_ * V2im) - y_int_[0]; - f_int_[1] = YReMat_ * (V1im - V2im) + (YImMatDi_ * V1re + YImMatOff_ * V2re) - y_int_[1]; - - // row2 - f_int_[2] = YReMat_ * (V2re - V1re) - (YImMatOff_ * V1im + YImMatDi_ * V2im) - y_int_[2]; - f_int_[3] = YReMat_ * (V2im - V1im) + (YImMatOff_ * V1re + YImMatDi_ * V2re) - y_int_[3]; - - return 0; - } - - template - int TransmissionLine::evaluateExternalResidual() - { - // input - *f_ext_[0] += y_int_[0]; - *f_ext_[1] += y_int_[1]; - - *f_ext_[2] += y_int_[2]; - *f_ext_[3] += y_int_[3]; - // ouput - *f_ext_[4] += -y_int_[0]; - *f_ext_[5] += -y_int_[1]; - - *f_ext_[6] += -y_int_[2]; - *f_ext_[7] += -y_int_[3]; - - return 0; - } - - /** - * @brief Generate Jacobian for Transmission Line - * - * @tparam ScalarT - * @tparam IdxT - * @return int - */ - template - int TransmissionLine::evaluateJacobian() - { - this->zeroJacMatrix(); - - // Create dF/dy - std::vector rtemp{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11}; - std::vector ctemp{8, 9, 10, 11, 8, 9, 10, 11, 8, 9, 10, 11}; - std::vector vals{1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0}; - this->setJacValues(rtemp, ctemp, vals); - - std::vector ccord{0, 1, 2, 3, 4, 5, 6, 7}; - - std::vector rcord(ccord.size(), 8); - vals = {YReMat_, -YImMatDi_, -YReMat_, -YImMatOff_, -YReMat_, YImMatDi_, YReMat_, YImMatOff_}; - this->setJacValues(rtemp, ctemp, vals); - - std::fill(rcord.begin(), rcord.end(), 9); - vals = {YImMatDi_, YReMat_, YImMatOff_, -YReMat_, -YImMatDi_, -YReMat_, -YImMatOff_, YReMat_}; - this->setJacValues(rtemp, ctemp, vals); - - std::fill(rcord.begin(), rcord.end(), 10); - vals = {-YReMat_, -YImMatDi_, YReMat_, -YImMatOff_, YReMat_, YImMatDi_, -YReMat_, YImMatOff_}; - this->setJacValues(rtemp, ctemp, vals); - - std::fill(rcord.begin(), rcord.end(), 11); - vals = {YImMatDi_, -YReMat_, YImMatOff_, YReMat_, -YImMatDi_, YReMat_, -YImMatOff_, -YReMat_}; - this->setJacValues(rtemp, ctemp, vals); - - return 0; - } - - template - int TransmissionLine::evaluateIntegrand() - { - return 0; - } - - template - int TransmissionLine::initializeAdjoint() - { - return 0; - } - - template - int TransmissionLine::evaluateAdjointResidual() - { - return 0; - } - - template - int TransmissionLine::evaluateAdjointIntegrand() + namespace PowerElectronics { - return 0; - } - - template - CircuitComponent* TransmissionLine::clone() const - { - return new TransmissionLine(*this); - } - - // Available template instantiations - template class TransmissionLine; - template class TransmissionLine; - template class TransmissionLine; - template class TransmissionLine; - + /*! + * @brief Constructor for a TransmissionLine model + * + * Calls default ModelEvaluatorImpl constructor. + * + * This is the Medium distance form with the use of the admittance matrix. + * Since the line is of medium length then there is no real part for shunt admittance + * @todo needs to used in a model + * @todo test for correctness + */ + template + TransmissionLine::TransmissionLine(IdxT id, RealT R, RealT X, RealT B) + : R_(R), + X_(X), + B_(B) + { + // internals [Iret1, Iimt1, Iret2, Iimt2] + // externals [Vre11, Vim11, Vre12, Vim12, Vre21, Vim21, Vre22, Vim22] + size_ = 12; + n_intern_ = 4; + n_extern_ = 8; + extern_indices_ = {0, 1, 2, 3, 4, 5, 6, 7}; + idc_ = id; + nnz_ = 44; + + RealT magImpendence = 1.0 / (R_ * R_ + X_ * X_); + YReMat_ = magImpendence * R_; + YImMatOff_ = magImpendence * X_; + YImMatDi_ = B_ / (2.0) - YImMatOff_; + } + + template + TransmissionLine::~TransmissionLine() + { + } + + /** + * Initialization of the grid model + */ + template + int TransmissionLine::initialize() + { + return 0; + } + + /* + * \brief Identify differential variables + */ + template + int TransmissionLine::tagDifferentiable() + { + // All variables are algebraics + std::fill(tag_.begin(), tag_.end(), false); + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int TransmissionLine::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Evaluate residual of transmission line + * + * The complex admittance matrix is: + * [[ Y/2 + 1/Z, -1/Z]; + * [ -1/Z, Y/2 + 1/Z]] = + * [[R/|Z|, -R/|Z|]; + * [-R/|Z|, R/|Z|]] + + * i [[B/2 - X/|Z|, X/|Z|]; + * [X/|Z|, B/2 - X/|Z|]] + * = Dre + i Dim + * + * Then + * Ire = Dre Vre - Dim Vim + * Iim = Dre Vim + Dim Vre + * + * To express this for Modified Nodal Analysis the Voltages of the admittance matrix are put into voltage drops + */ + template + int TransmissionLine::evaluateInternalResidual() + { + // Voltage drop accross terminals + ScalarT V1re = *y_ext_[0] - *y_ext_[4]; + ScalarT V1im = *y_ext_[1] - *y_ext_[5]; + ScalarT V2re = *y_ext_[2] - *y_ext_[6]; + ScalarT V2im = *y_ext_[3] - *y_ext_[7]; + + // Internal variables + // row 1 + f_int_[0] = YReMat_ * (V1re - V2re) - (YImMatDi_ * V1im + YImMatOff_ * V2im) - y_int_[0]; + f_int_[1] = YReMat_ * (V1im - V2im) + (YImMatDi_ * V1re + YImMatOff_ * V2re) - y_int_[1]; + + // row2 + f_int_[2] = YReMat_ * (V2re - V1re) - (YImMatOff_ * V1im + YImMatDi_ * V2im) - y_int_[2]; + f_int_[3] = YReMat_ * (V2im - V1im) + (YImMatOff_ * V1re + YImMatDi_ * V2re) - y_int_[3]; + + return 0; + } + + template + int TransmissionLine::evaluateExternalResidual() + { + // input + *f_ext_[0] += y_int_[0]; + *f_ext_[1] += y_int_[1]; + + *f_ext_[2] += y_int_[2]; + *f_ext_[3] += y_int_[3]; + // ouput + *f_ext_[4] += -y_int_[0]; + *f_ext_[5] += -y_int_[1]; + + *f_ext_[6] += -y_int_[2]; + *f_ext_[7] += -y_int_[3]; + + return 0; + } + + /** + * @brief Generate Jacobian for Transmission Line + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ + template + int TransmissionLine::evaluateJacobian() + { + this->zeroJacMatrix(); + + // Create dF/dy + std::vector rtemp{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11}; + std::vector ctemp{8, 9, 10, 11, 8, 9, 10, 11, 8, 9, 10, 11}; + std::vector vals{1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0}; + this->setJacValues(rtemp, ctemp, vals); + + std::vector ccord{0, 1, 2, 3, 4, 5, 6, 7}; + + std::vector rcord(ccord.size(), 8); + vals = {YReMat_, -YImMatDi_, -YReMat_, -YImMatOff_, -YReMat_, YImMatDi_, YReMat_, YImMatOff_}; + this->setJacValues(rtemp, ctemp, vals); + + std::fill(rcord.begin(), rcord.end(), 9); + vals = {YImMatDi_, YReMat_, YImMatOff_, -YReMat_, -YImMatDi_, -YReMat_, -YImMatOff_, YReMat_}; + this->setJacValues(rtemp, ctemp, vals); + + std::fill(rcord.begin(), rcord.end(), 10); + vals = {-YReMat_, -YImMatDi_, YReMat_, -YImMatOff_, YReMat_, YImMatDi_, -YReMat_, YImMatOff_}; + this->setJacValues(rtemp, ctemp, vals); + + std::fill(rcord.begin(), rcord.end(), 11); + vals = {YImMatDi_, -YReMat_, YImMatOff_, YReMat_, -YImMatDi_, YReMat_, -YImMatOff_, -YReMat_}; + this->setJacValues(rtemp, ctemp, vals); + + return 0; + } + + template + int TransmissionLine::evaluateIntegrand() + { + return 0; + } + + template + int TransmissionLine::initializeAdjoint() + { + return 0; + } + + template + int TransmissionLine::evaluateAdjointResidual() + { + return 0; + } + + template + int TransmissionLine::evaluateAdjointIntegrand() + { + return 0; + } + + template + Component* TransmissionLine::clone() const + { + return new TransmissionLine(*this); + } + + // Available template instantiations + template class TransmissionLine; + template class TransmissionLine; + template class TransmissionLine; + template class TransmissionLine; + + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp b/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp index 7a0bfb230..7a21c4bf7 100644 --- a/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp +++ b/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp @@ -1,8 +1,7 @@ - #pragma once -#include +#include namespace GridKit { @@ -12,68 +11,71 @@ namespace GridKit namespace GridKit { - /*! - * @brief Declaration of a TransmissionLine class. - * - * Model from Adam Birchfield paper (medium distances < 2km). - * See also textbooks "Power System Analysis" by Grainger and "Power System Dynamics and Stability" by Sauer & Pai - * - * @note Not used in the Microgrid model. - */ - template - class TransmissionLine : public CircuitComponent + namespace PowerElectronics { - using RealT = typename CircuitComponent::RealT; + /*! + * @brief Declaration of a TransmissionLine class. + * + * Model from Adam Birchfield paper (medium distances < 2km). + * See also textbooks "Power System Analysis" by Grainger and "Power System Dynamics and Stability" by Sauer & Pai + * + * @note Not used in the Microgrid model. + */ + template + class TransmissionLine : public Component + { + using RealT = typename Component::RealT; - using CircuitComponent::size_; - using CircuitComponent::nnz_; - using CircuitComponent::time_; - using CircuitComponent::alpha_; - using CircuitComponent::y_ext_; - using CircuitComponent::y_int_; - using CircuitComponent::yp_ext_; - using CircuitComponent::yp_int_; - using CircuitComponent::abs_tol_; - using CircuitComponent::tag_; - using CircuitComponent::f_ext_; - using CircuitComponent::f_int_; - using CircuitComponent::g_; - using CircuitComponent::yB_; - using CircuitComponent::ypB_; - using CircuitComponent::fB_; - using CircuitComponent::gB_; - using CircuitComponent::param_; - using CircuitComponent::idc_; + using Component::size_; + using Component::nnz_; + using Component::time_; + using Component::alpha_; + using Component::y_ext_; + using Component::y_int_; + using Component::yp_ext_; + using Component::yp_int_; + using Component::abs_tol_; + using Component::tag_; + using Component::f_ext_; + using Component::f_int_; + using Component::g_; + using Component::yB_; + using Component::ypB_; + using Component::fB_; + using Component::gB_; + using Component::param_; + using Component::idc_; - using CircuitComponent::extern_indices_; - using CircuitComponent::n_extern_; - using CircuitComponent::n_intern_; + using Component::extern_indices_; + using Component::n_extern_; + using Component::n_intern_; - public: - TransmissionLine(IdxT id, RealT R, RealT X, RealT B); - virtual ~TransmissionLine(); + public: + TransmissionLine(IdxT id, RealT R, RealT X, RealT B); + virtual ~TransmissionLine(); - int initialize(); - int tagDifferentiable(); - int setAbsoluteTolerance(RealT); - int evaluateInternalResidual() final; - int evaluateExternalResidual() final; - int evaluateJacobian(); - int evaluateIntegrand(); + int initialize(); + int tagDifferentiable(); + int setAbsoluteTolerance(RealT); + int evaluateInternalResidual() final; + int evaluateExternalResidual() final; + int evaluateJacobian(); + int evaluateIntegrand(); - int initializeAdjoint(); - int evaluateAdjointResidual(); - // int evaluateAdjointJacobian(); - int evaluateAdjointIntegrand(); + int initializeAdjoint(); + int evaluateAdjointResidual(); + // int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); - CircuitComponent* clone() const; + Component* clone() const; - private: - RealT R_; - RealT X_; - RealT B_; - RealT YReMat_; - RealT YImMatDi_; - RealT YImMatOff_; - }; + private: + RealT R_; + RealT X_; + RealT B_; + RealT YReMat_; + RealT YImMatDi_; + RealT YImMatOff_; + }; + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp index 900d2fc64..e29d13de4 100644 --- a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp +++ b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp @@ -1,5 +1,4 @@ - #include "VoltageSource.hpp" #include @@ -8,150 +7,151 @@ namespace GridKit { - - /*! - * @brief Constructor for a constant VoltageSource model - * - * Calls default ModelEvaluatorImpl constructor. - */ - - template - VoltageSource::VoltageSource(IdxT id, RealT V, NodeT* node1, NodeT* node2) - : V_(V), node1_(node1), node2_(node2) - { - assert(node1_->size() == 1); - assert(node2_->size() == 1); - size_ = 3; - n_intern_ = 1; - n_extern_ = 2; - extern_indices_ = {0, 1}; - idc_ = id; - nnz_ = 4; - } - - template - VoltageSource::~VoltageSource() - { - } - - /** - * Initialization of the grid model - */ - template - int VoltageSource::initialize() - { - return 0; - } - - /* - * \brief Identify differential variables - */ - template - int VoltageSource::tagDifferentiable() - { - // All variables are algebraics - std::fill(tag_.begin(), tag_.end(), false); - return 0; - } - - /** - * @brief Compute the absolute tolerance for each variable in the model - * - * @param rel_tol The relative tolerance which can be used to pick the - * absolute tolerance. - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return int 0 if successful, non-zero otherwise. - * - * This represents a "noise" level close to zero for which pure relative - * error cannot be used. - */ - template - int VoltageSource::setAbsoluteTolerance(RealT rel_tol) - { - abs_tol_.setToConst(static_cast(rel_tol)); - return 0; - } - - /** - * @brief Evaluate resisdual of component - */ - template - int VoltageSource::evaluateInternalResidual() - { - f_int_[0] = *y_ext_[1] - *y_ext_[0] - V_; - return 0; - } - - template - int VoltageSource::evaluateExternalResidual() - { - // input - *f_ext_[0] += -y_int_[0]; - // ouput - *f_ext_[1] += y_int_[0]; - return 0; - } - - template - int VoltageSource::evaluateJacobian() - { - this->zeroJacMatrix(); - - // Create dF/dy - std::vector rcord{0, 1, 2, 2}; - std::vector ccord{2, 2, 0, 1}; - std::vector vals{-1.0, 1.0, -1.0, 1.0}; - this->setJacValues(rcord, ccord, vals); - - return 0; - } - - template - int VoltageSource::allocate() - { - CircuitComponent::allocate(); - - this->setExternalConnectionNodes(0, node1_->getNodeConnection(0)); - this->setExternalConnectionNodes(1, node2_->getNodeConnection(0)); - - return 0; - } - - template - int VoltageSource::evaluateIntegrand() - { - return 0; - } - - template - int VoltageSource::initializeAdjoint() + namespace PowerElectronics { - return 0; - } - - template - int VoltageSource::evaluateAdjointResidual() - { - return 0; - } - - template - int VoltageSource::evaluateAdjointIntegrand() - { - return 0; - } - - template - CircuitComponent* VoltageSource::clone() const - { - return new VoltageSource(*this); - } - - // Available template instantiations - template class VoltageSource; - template class VoltageSource; - template class VoltageSource; - template class VoltageSource; - + /*! + * @brief Constructor for a constant VoltageSource model + * + * Calls default ModelEvaluatorImpl constructor. + */ + template + VoltageSource::VoltageSource(IdxT id, RealT V, NodeT* node1, NodeT* node2) + : V_(V), node1_(node1), node2_(node2) + { + assert(node1_->size() == 1); + assert(node2_->size() == 1); + size_ = 3; + n_intern_ = 1; + n_extern_ = 2; + extern_indices_ = {0, 1}; + idc_ = id; + nnz_ = 4; + } + + template + VoltageSource::~VoltageSource() + { + } + + /** + * Initialization of the grid model + */ + template + int VoltageSource::initialize() + { + return 0; + } + + /* + * \brief Identify differential variables + */ + template + int VoltageSource::tagDifferentiable() + { + // All variables are algebraics + std::fill(tag_.begin(), tag_.end(), false); + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int VoltageSource::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Evaluate resisdual of component + */ + template + int VoltageSource::evaluateInternalResidual() + { + f_int_[0] = *y_ext_[1] - *y_ext_[0] - V_; + return 0; + } + + template + int VoltageSource::evaluateExternalResidual() + { + // input + *f_ext_[0] += -y_int_[0]; + // ouput + *f_ext_[1] += y_int_[0]; + return 0; + } + + template + int VoltageSource::evaluateJacobian() + { + this->zeroJacMatrix(); + + // Create dF/dy + std::vector rcord{0, 1, 2, 2}; + std::vector ccord{2, 2, 0, 1}; + std::vector vals{-1.0, 1.0, -1.0, 1.0}; + this->setJacValues(rcord, ccord, vals); + + return 0; + } + + template + int VoltageSource::allocate() + { + Component::allocate(); + + this->setExternalConnectionNodes(0, node1_->getNodeConnection(0)); + this->setExternalConnectionNodes(1, node2_->getNodeConnection(0)); + + return 0; + } + + template + int VoltageSource::evaluateIntegrand() + { + return 0; + } + + template + int VoltageSource::initializeAdjoint() + { + return 0; + } + + template + int VoltageSource::evaluateAdjointResidual() + { + return 0; + } + + template + int VoltageSource::evaluateAdjointIntegrand() + { + return 0; + } + + template + Component* VoltageSource::clone() const + { + return new VoltageSource(*this); + } + + // Available template instantiations + template class VoltageSource; + template class VoltageSource; + template class VoltageSource; + template class VoltageSource; + + } // namespace PowerElectronics } // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp index 0cc94f66c..e80ce5008 100644 --- a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp +++ b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp @@ -1,75 +1,71 @@ - #pragma once -#include +#include #include namespace GridKit { - template - class BaseBus; -} - -namespace GridKit -{ - /*! - * @brief Declaration of a VoltageSource class. - * - */ - template - class VoltageSource : public CircuitComponent + namespace PowerElectronics { - using RealT = typename CircuitComponent::RealT; - using NodeT = typename PowerElectronics::NodeBase; + /*! + * @brief Declaration of a VoltageSource class. + * + */ + template + class VoltageSource : public Component + { + using RealT = typename Component::RealT; + using NodeT = typename PowerElectronics::NodeBase; - using CircuitComponent::size_; - using CircuitComponent::nnz_; - using CircuitComponent::time_; - using CircuitComponent::alpha_; - using CircuitComponent::y_ext_; - using CircuitComponent::y_int_; - using CircuitComponent::yp_ext_; - using CircuitComponent::yp_int_; - using CircuitComponent::tag_; - using CircuitComponent::abs_tol_; - using CircuitComponent::f_ext_; - using CircuitComponent::f_int_; - using CircuitComponent::g_; - using CircuitComponent::yB_; - using CircuitComponent::ypB_; - using CircuitComponent::fB_; - using CircuitComponent::gB_; - using CircuitComponent::param_; - using CircuitComponent::idc_; + using Component::size_; + using Component::nnz_; + using Component::time_; + using Component::alpha_; + using Component::y_ext_; + using Component::y_int_; + using Component::yp_ext_; + using Component::yp_int_; + using Component::tag_; + using Component::abs_tol_; + using Component::f_ext_; + using Component::f_int_; + using Component::g_; + using Component::yB_; + using Component::ypB_; + using Component::fB_; + using Component::gB_; + using Component::param_; + using Component::idc_; - using CircuitComponent::extern_indices_; - using CircuitComponent::n_extern_; - using CircuitComponent::n_intern_; + using Component::extern_indices_; + using Component::n_extern_; + using Component::n_intern_; - public: - VoltageSource(IdxT id, RealT V, NodeT* node1, NodeT* node2); - virtual ~VoltageSource(); + public: + VoltageSource(IdxT id, RealT V, NodeT* node1, NodeT* node2); + virtual ~VoltageSource(); - int initialize(); - int allocate() final; - int tagDifferentiable(); - int setAbsoluteTolerance(RealT); - int evaluateInternalResidual() final; - int evaluateExternalResidual() final; - int evaluateJacobian() final; - int evaluateIntegrand(); + int initialize(); + int allocate() final; + int tagDifferentiable(); + int setAbsoluteTolerance(RealT); + int evaluateInternalResidual() final; + int evaluateExternalResidual() final; + int evaluateJacobian() final; + int evaluateIntegrand(); - int initializeAdjoint(); - int evaluateAdjointResidual(); - // int evaluateAdjointJacobian(); - int evaluateAdjointIntegrand(); + int initializeAdjoint(); + int evaluateAdjointResidual(); + // int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); - CircuitComponent* clone() const; + Component* clone() const; - private: - RealT V_; - NodeT* node1_; - NodeT* node2_; - }; + private: + RealT V_; + NodeT* node1_; + NodeT* node2_; + }; + } // namespace PowerElectronics } // namespace GridKit diff --git a/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp b/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp index 649146bab..7ee73f9de 100644 --- a/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp +++ b/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp @@ -9,7 +9,7 @@ #include #include -#include +#include #include /** @@ -22,7 +22,7 @@ int main(int /* argc */, char const** /* argv */) { - GridKit::DistributedGeneratorParameters parms; + GridKit::PowerElectronics::DistributedGeneratorParameters parms; // Parameters from MATLAB Microgrid code for first DG parms.wb_ = 2.0 * std::numbers::pi_v * 50.0; parms.wc_ = 31.41; @@ -46,7 +46,7 @@ int main(int /* argc */, char const** /* argv */) using Bus = GridKit::PowerElectronics::MicrogridBus; Bus bus; - GridKit::DistributedGenerator dg(0, parms, false, &dg_signal, &bus); + GridKit::PowerElectronics::DistributedGenerator dg(0, parms, false, &dg_signal, &bus); std::vector t1(16, 0.0); std::vector t2{ @@ -79,7 +79,7 @@ int main(int /* argc */, char const** /* argv */) for (size_t idx : dg.getExternIndices()) { - GridKit::ExternalConnection connection{ + GridKit::PowerElectronics::ExternalConnection connection{ .y_ = &t2[idx], .yp_ = &t1[idx], .f_ = &res[idx], diff --git a/examples/PowerElectronics/ExamplesHelper/MicrogridNetwork.hpp b/examples/PowerElectronics/ExamplesHelper/MicrogridNetwork.hpp index 4a543e7d7..8425cc87e 100644 --- a/examples/PowerElectronics/ExamplesHelper/MicrogridNetwork.hpp +++ b/examples/PowerElectronics/ExamplesHelper/MicrogridNetwork.hpp @@ -13,7 +13,7 @@ #include #include #include -#include +#include /* * Contains components and nodes that make up the scaled microgrid network. @@ -27,11 +27,11 @@ struct ScaleMicrogridNetwork using SignalNode = GridKit::PowerElectronics::SignalNode; using Bus = GridKit::PowerElectronics::MicrogridBus; - using BusDQ = GridKit::MicrogridBusDQ; - using DGGenerator = GridKit::DistributedGenerator; - using Line = GridKit::MicrogridLine; - using Load = GridKit::MicrogridLoad; - using GenParams = GridKit::DistributedGeneratorParameters; + using BusDQ = GridKit::PowerElectronics::MicrogridBusDQ; + using DGGenerator = GridKit::PowerElectronics::DistributedGenerator; + using Line = GridKit::PowerElectronics::MicrogridLine; + using Load = GridKit::PowerElectronics::MicrogridLoad; + using GenParams = GridKit::PowerElectronics::DistributedGeneratorParameters; size_t model_id_next; size_t N_size; diff --git a/examples/PowerElectronics/ExamplesHelper/SystemAssembler.hpp b/examples/PowerElectronics/ExamplesHelper/SystemAssembler.hpp index 3565ed63d..e46400a56 100644 --- a/examples/PowerElectronics/ExamplesHelper/SystemAssembler.hpp +++ b/examples/PowerElectronics/ExamplesHelper/SystemAssembler.hpp @@ -29,10 +29,11 @@ * @p network have been added to @p sys_model. * * @note This function only assembles the network into the system model. It - * does not call PowerElectronicsModel::allocate(). + * does not call GridKit::PowerElectronics::SystemModel::allocate(). */ template -void assembleSystem(ScaleMicrogridNetwork& network, GridKit::PowerElectronicsModel& sys_model) +void assembleSystem(ScaleMicrogridNetwork& network, + GridKit::PowerElectronics::SystemModel& sys_model) { size_t N_size = network.N_size; diff --git a/examples/PowerElectronics/Microgrid/Microgrid.cpp b/examples/PowerElectronics/Microgrid/Microgrid.cpp index 7dbe48096..0f568ccd4 100644 --- a/examples/PowerElectronics/Microgrid/Microgrid.cpp +++ b/examples/PowerElectronics/Microgrid/Microgrid.cpp @@ -2,7 +2,7 @@ #include #include -#include +#include #include #include @@ -19,7 +19,7 @@ int main(int /* argc */, char const** /* argv */) bool debug_output = true; // Create model - auto* sysmodel = new GridKit::PowerElectronicsModel(use_jac); + auto* sysmodel = new GridKit::PowerElectronics::SystemModel(use_jac); // Build the four-generator microgrid network. size_t N_size = 2; diff --git a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp index a802d84c9..4d3693442 100644 --- a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp +++ b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp @@ -1,5 +1,4 @@ - #include #include @@ -8,7 +7,7 @@ #include #include #include -#include +#include #include #include #include @@ -21,7 +20,7 @@ int main(int /* argc */, char const** /* argv */) // TODO:setup as named parameters // Create circuit model - GridKit::PowerElectronicsModel sysmodel(use_jac); + GridKit::PowerElectronics::SystemModel sysmodel(use_jac); size_t idoff = 0; @@ -41,17 +40,17 @@ int main(int /* argc */, char const** /* argv */) sysmodel.addNode(&bus_ir); // inductor - GridKit::Inductor* induct = new GridKit::Inductor(idoff, linit, &bus_ir, &bus_iv); + auto* induct = new GridKit::PowerElectronics::Inductor(idoff, linit, &bus_ir, &bus_iv); sysmodel.addComponent(induct); // resistor idoff++; - GridKit::Resistor* resis = new GridKit::Resistor(idoff, rinit, &bus_vr, &bus_ir); + auto* resis = new GridKit::PowerElectronics::Resistor(idoff, rinit, &bus_vr, &bus_ir); sysmodel.addComponent(resis); // voltage source idoff++; - GridKit::VoltageSource* vsource = new GridKit::VoltageSource(idoff, vinit, &bus_iv, &bus_vr); + auto* vsource = new GridKit::PowerElectronics::VoltageSource(idoff, vinit, &bus_iv, &bus_vr); sysmodel.addComponent(vsource); sysmodel.allocate(); diff --git a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp index 4c689cbc8..134358745 100644 --- a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp +++ b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp @@ -2,7 +2,7 @@ #include #include -#include +#include #include #include #include @@ -56,7 +56,7 @@ int main(int /* argc */, char const** /* argv */) */ int test(index_type Nsize, real_type error_tol, bool debug_output) { - using namespace GridKit; + using namespace GridKit::PowerElectronics; bool use_jac = true; @@ -67,7 +67,7 @@ int test(index_type Nsize, real_type error_tol, bool debug_output) real_type abs_tol = SCALE_MICROGRID_ABS_TOL; // Create circuit model - auto* sys_model = new PowerElectronicsModel(use_jac); + auto* sys_model = new SystemModel(use_jac); const std::vector* true_vec = &answer_key_N8; diff --git a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogridArbitrary.cpp b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogridArbitrary.cpp index 1edc577e2..5e6e04803 100644 --- a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogridArbitrary.cpp +++ b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogridArbitrary.cpp @@ -55,7 +55,7 @@ int main(int argc, char const* argv[]) */ int printMicrogridSystems(index_type N_size) { - using namespace GridKit; + using namespace GridKit::PowerElectronics; bool use_jac = true; @@ -65,7 +65,7 @@ int printMicrogridSystems(index_type N_size) real_type rel_tol = SCALE_MICROGRID_REL_TOL; // Create circuit model - PowerElectronicsModel sys_model(use_jac); + SystemModel sys_model(use_jac); // Ensure minimum size requirement if (N_size < 1) diff --git a/tests/UnitTests/PowerElectronics/CMakeLists.txt b/tests/UnitTests/PowerElectronics/CMakeLists.txt index 94cd4fd91..dc3c16737 100644 --- a/tests/UnitTests/PowerElectronics/CMakeLists.txt +++ b/tests/UnitTests/PowerElectronics/CMakeLists.txt @@ -1,4 +1,4 @@ -add_executable(test_power_electronics_node runCircuitNodeTests.cpp) +add_executable(test_power_electronics_node runNodeTests.cpp) target_link_libraries( test_power_electronics_node PRIVATE GridKit::power_electronics_circuit_node GridKit::testing) diff --git a/tests/UnitTests/PowerElectronics/ComponentCloneTests.hpp b/tests/UnitTests/PowerElectronics/ComponentCloneTests.hpp index 77ea9a623..d9f6d33da 100644 --- a/tests/UnitTests/PowerElectronics/ComponentCloneTests.hpp +++ b/tests/UnitTests/PowerElectronics/ComponentCloneTests.hpp @@ -13,18 +13,18 @@ namespace GridKit namespace Testing { template - class CircuitComponentCloneTests + class ComponentCloneTests { using SignalNode = PowerElectronics::SignalNode; using Bus = PowerElectronics::MicrogridBus; - using BusDQ = MicrogridBusDQ; - using Generator = DistributedGenerator; - using GeneratorParameters = DistributedGeneratorParameters; - using Line = MicrogridLine; - using Load = MicrogridLoad; + using BusDQ = PowerElectronics::MicrogridBusDQ; + using Generator = PowerElectronics::DistributedGenerator; + using GeneratorParameters = PowerElectronics::DistributedGeneratorParameters; + using Line = PowerElectronics::MicrogridLine; + using Load = PowerElectronics::MicrogridLoad; public: - CircuitComponentCloneTests() + ComponentCloneTests() { /************************************************************************** * Construct Network Nodes @@ -77,7 +77,7 @@ namespace GridKit bus_dq_->allocate(); } - ~CircuitComponentCloneTests() + ~ComponentCloneTests() { delete generator_; delete line_; @@ -125,7 +125,7 @@ namespace GridKit template bool verifyComponentClone(ComponentT& component) { - using RealT = typename CircuitComponent::RealT; + using RealT = typename PowerElectronics::Component::RealT; bool success = true; diff --git a/tests/UnitTests/PowerElectronics/CircuitNodeTests.hpp b/tests/UnitTests/PowerElectronics/NodeTests.hpp similarity index 83% rename from tests/UnitTests/PowerElectronics/CircuitNodeTests.hpp rename to tests/UnitTests/PowerElectronics/NodeTests.hpp index bc380f972..fa16ee018 100644 --- a/tests/UnitTests/PowerElectronics/CircuitNodeTests.hpp +++ b/tests/UnitTests/PowerElectronics/NodeTests.hpp @@ -2,7 +2,7 @@ #include #include -#include +#include #include #include @@ -24,10 +24,10 @@ namespace GridKit ScalarT V{1.0}; - CircuitNode* node = nullptr; + PowerElectronics::Node* node = nullptr; // Default construct - node = new CircuitNode(); + node = new PowerElectronics::Node(); node->allocate(); node->initialize(); success *= isEqual(node->V(), static_cast(0)); @@ -35,11 +35,12 @@ namespace GridKit delete node; // Construct with initial voltage - node = new CircuitNode(V); + node = new PowerElectronics::Node(V); node->allocate(); node->initialize(); success *= isEqual(node->V(), V); success *= isEqual(node->I(), static_cast(0)); + delete node; node = nullptr; @@ -55,7 +56,7 @@ namespace GridKit ScalarT V{1.0}; ScalarT I{1.0}; - CircuitNode node(V); + PowerElectronics::Node node(V); node.allocate(); node.initialize(); success *= isEqual(node.V(), V); diff --git a/tests/UnitTests/PowerElectronics/runComponentCloneTests.cpp b/tests/UnitTests/PowerElectronics/runComponentCloneTests.cpp index dfb5fa97f..83f052ec9 100644 --- a/tests/UnitTests/PowerElectronics/runComponentCloneTests.cpp +++ b/tests/UnitTests/PowerElectronics/runComponentCloneTests.cpp @@ -2,7 +2,7 @@ int main() { - GridKit::Testing::CircuitComponentCloneTests tests; + GridKit::Testing::ComponentCloneTests tests; GridKit::Testing::TestingResults result; diff --git a/tests/UnitTests/PowerElectronics/runCircuitNodeTests.cpp b/tests/UnitTests/PowerElectronics/runNodeTests.cpp similarity index 86% rename from tests/UnitTests/PowerElectronics/runCircuitNodeTests.cpp rename to tests/UnitTests/PowerElectronics/runNodeTests.cpp index 600babd50..25535b2fa 100644 --- a/tests/UnitTests/PowerElectronics/runCircuitNodeTests.cpp +++ b/tests/UnitTests/PowerElectronics/runNodeTests.cpp @@ -1,4 +1,4 @@ -#include "CircuitNodeTests.hpp" +#include "NodeTests.hpp" int main() {