Commit 04644796 authored by Martin Alonso Moraga's avatar Martin Alonso Moraga

Extend DPsim components by DAE solver functions

parent 4dac5bee
......@@ -63,6 +63,8 @@ if(WITH_SUNDIALS)
DAE/DAE_EMT_RL1.cpp
DAE/DAE_EMT_RC1.cpp
DAE/DAE_EMT_RLC1.cpp
DAE/DAE_EMT_Slack_PiLine_RXLoad.cpp
DAE/MNA_EMT_Slack_PiLine_RXLoad.cpp
)
endif()
......
......@@ -58,7 +58,7 @@ int main(int argc, char* argv[])
// Define system topology
auto sys = SystemTopology(50, SystemNodeList{n1, n2, n3}, SystemComponentList{vs, r1, l1, c1});
vs->setInitialCurrent(0.99993);
// vs->setInitialCurrent(0.99993);
// Logger
auto logger = DataLogger::make(simName);
......
......@@ -45,7 +45,6 @@ namespace DPsim {
Int mNEQ;
/// Components of the Problem
CPS::DAEInterface::List mDAEComponents;
/// Nodes of the Problem
typename CPS::SimNode<VarType>::List mNodes;
......
......@@ -23,7 +23,7 @@ DAESolver<VarType>::DAESolver(String name,
mTimestep(dt) {
// Defines offset vector of the residual which is composed as follows:
// mOffset[0] = # nodal voltage equations
// mOffset[0] = # nodal voltage equations (1 for SinglePhase nodes, 3 for ThreePhase nodes)
// mOffset[1] = # of components equations (1 for inductor, cap and voltagesource)
mOffsets.push_back(0);
mOffsets.push_back(0);
......@@ -41,7 +41,12 @@ DAESolver<VarType>::DAESolver(String name,
throw CPS::Exception();
}
mNodes.push_back(node);
mNEQ +=1;
if (node->phaseType() == PhaseType::Single) {
mNEQ += 1;
}
else if (node->phaseType() == PhaseType::ABC) {
mNEQ += 3;
}
mSLog->info("Added node {:s}", node->name());;
}
}
......@@ -114,18 +119,50 @@ void DAESolver<VarType>::initialize(Real t0) {
// capturing a returned array/pointer
sval = N_VGetArrayPointer(mStateVector);
s_dtval = N_VGetArrayPointer_Serial(mDerivativeStateVector);
// Initialize nodal voltages of state vector
for (auto node : mNodes) {
// Initialize nodal voltages of state vector
Real tempVolt = std::real(node->initialSingleVoltage(node->phaseType()));
sval[counter] = tempVolt;
s_dtval[counter] = 0;
mSLog->info(
"Added node '{:s}' to state vector, init voltage = {:f}V",
node->name(), sval[counter]);
mSLog->info(
"Added derivative of the voltage node of '{:s}' to derivative state vector, initial value = {:f}V",
node->name(), s_dtval[counter]);
counter++;
if (node->phaseType() == PhaseType::Single) {
Real tempVolt = std::real(node->initialSingleVoltage(PhaseType::Single));
sval[counter] = tempVolt;
s_dtval[counter] = 0;
mSLog->info(
"Added node '{:s}' to state vector, init voltage = {:f}V"
"\nAdded derivative of the voltage node of '{:s}' to derivative state vector, initial value = {:f}",
node->name(), sval[counter], node->name(), s_dtval[counter]
);
counter++;
}
else if (node->phaseType() == PhaseType::ABC) {
Real tempVolt_phase_A = std::real(node->initialSingleVoltage(PhaseType::A));
Real tempVolt_phase_B = std::real(node->initialSingleVoltage(PhaseType::B));
Real tempVolt_phase_C = std::real(node->initialSingleVoltage(PhaseType::C));
sval[counter] = tempVolt_phase_A;
s_dtval[counter] = 0;
mSLog->info(
"Added node '{:s}'-phase_A to state vector, init voltage = {:f}V"
"\nAdded derivative of the voltage node of '{:s}'-phase_A to derivative state vector, initial value = {:f}",
node->name(), sval[counter], node->name(), s_dtval[counter]
);
counter++;
sval[counter] = tempVolt_phase_B;
s_dtval[counter] = 0;
mSLog->info(
"Added node '{:s}'-phase_B to state vector, init voltage = {:f}V"
"\nAdded derivative of the voltage node of '{:s}'-phase_B to derivative state vector, initial value = {:f}",
node->name(), sval[counter], node->name(), s_dtval[counter]
);
counter++;
sval[counter] = tempVolt_phase_C;
s_dtval[counter] = 0;
mSLog->info(
"Added node '{:s}'-phase_C to state vector, init voltage = {:f}V"
"\nAdded derivative of the voltage node of '{:s}'-phase_C to derivative state vector, initial value = {:f}",
node->name(), sval[counter], node->name(), s_dtval[counter]
);
counter++;
}
}
for (auto daeComp : mDAEComponents) {
// Initialize component voltages of state vector
......@@ -197,10 +234,19 @@ template <typename VarType>
int DAESolver<VarType>::residualFunction(realtype step_time,
N_Vector state, N_Vector dstate_dt, N_Vector resid)
{
mOffsets[0] = mNodes.size(); // Reset Offset of nodes
mOffsets[1] = 0; // Reset Offset of componentes
// Reset Offset of nodes
mOffsets[0] =0;
for (auto node : mNodes) {
if (node->phaseType() == PhaseType::Single) {
mOffsets[0] +=1;
}
else if (node->phaseType() == PhaseType::ABC) {
mOffsets[0] +=3;
}
}
mOffsets[1] = 0; // Reset Offset of componentes
//reset residual functions of nodes
//reset residual functions of nodes (nodal equations)
realtype *residual = NULL;
residual = N_VGetArrayPointer(resid);
for (int i=0; i<mOffsets[0]; i++) {
......@@ -237,10 +283,23 @@ Real DAESolver<VarType>::step(Real time) {
sval = N_VGetArrayPointer(mStateVector);
dstate_val = N_VGetArrayPointer(mDerivativeStateVector);
mOffsets[0] = 0; // Reset Offset of nodes
// Update voltage of nodes
for (auto node : mNodes) {
node->setVoltage(sval[mOffsets[0]]);
mOffsets[0] += 1;
if (node->phaseType() == PhaseType::Single) {
node->setVoltage(sval[mOffsets[0]]);
mOffsets[0] +=1;
}
else if (node->phaseType() == PhaseType::ABC) {
node->setVoltage(sval[mOffsets[0]], PhaseType::A);
mOffsets[0] +=1;
node->setVoltage(sval[mOffsets[0]], PhaseType::B);
mOffsets[0] +=1;
node->setVoltage(sval[mOffsets[0]], PhaseType::C);
mOffsets[0] +=1;
}
}
mOffsets[1] = 0; // Reset Offset of componentes
//update components
mOffsets[1] = mOffsets[0]; // Reset Offset of componentes
......
......@@ -10,6 +10,7 @@
#include <cps/SimPowerComp.h>
#include <cps/Solver/MNAInterface.h>
#include <cps/Solver/DAEInterface.h>
namespace CPS {
namespace EMT {
......@@ -25,6 +26,7 @@ namespace CPS {
class NetworkInjection :
public SimPowerComp<Real>,
public MNAInterface,
public DAEInterface,
public SharedFactory<NetworkInjection> {
private:
///
......@@ -48,6 +50,13 @@ namespace CPS {
// #### General ####
/// Initializes component from power flow data
void initializeFromPowerflow(Real frequency);
void setInitialCurrent(MatrixComp initCurrent) {
std::cout << "hola" <<std::endl;
Please register or sign in to reply
mIntfCurrent(0,0) = initCurrent(0,0).real();
mIntfCurrent(1,0) = initCurrent(1,0).real();
mIntfCurrent(2,0) = initCurrent(2,0).real();
std::cout << "hola" <<std::endl;
}
///
void setSourceValue(MatrixComp voltage);
///
......@@ -91,11 +100,16 @@ namespace CPS {
NetworkInjection& mNetworkInjection;
Attribute<Matrix>::Ptr mLeftVector;
};
//// #### DAE Section ####
///// Residual function for DAE Solver
//void daeResidual(double ttime, const double state[], const double dstate_dt[], double resid[], std::vector<int>& off);
/////Voltage Getter
//Complex daeInitialize();
// #### DAE Section ####
///
void daeInitialize(double time, double state[], double dstate_dt[], int& offset);
/// Residual function for DAE Solver
void daeResidual(double time, const double state[], const double dstate_dt[], double resid[], std::vector<int>& off);
///
void daePostStep(const double state[], const double dstate_dt[], int& offset);
///
int getNumberOfStateVariables() {return 3;}
};
}
}
......
......@@ -10,6 +10,7 @@
#include <cps/SimPowerComp.h>
#include <cps/Solver/MNAInterface.h>
#include <cps/Solver/DAEInterface.h>
#include <cps/Base/Base_Ph3_PiLine.h>
#include <cps/EMT/EMT_Ph3_Resistor.h>
#include <cps/EMT/EMT_Ph3_Inductor.h>
......@@ -25,6 +26,7 @@ namespace Ph3 {
class PiLine :
public SimPowerComp<Real>,
public MNAInterface,
public DAEInterface,
public Base::Ph3::PiLine,
public SharedFactory<PiLine> {
protected:
......@@ -116,6 +118,15 @@ namespace Ph3 {
void mnaTearApplyVoltageStamp(Matrix& voltageVector);
void mnaTearPostStep(Matrix voltage, Matrix current);*/
// #### DAE Section ####
///
void daeInitialize(double time, double state[], double dstate_dt[], int& counter);
///Residual Function for DAE Solver
void daeResidual(double time, const double state[], const double dstate_dt[], double resid[], std::vector<int>& off);
///
void daePostStep(const double state[], const double dstate_dt[], int& counter);
///
int getNumberOfStateVariables() {return 3;}
};
}
}
......
......@@ -10,6 +10,7 @@
#include <cps/SimPowerComp.h>
#include <cps/Solver/MNAInterface.h>
#include <cps/Solver/DAEInterface.h>
#include <cps/EMT/EMT_Ph3_Capacitor.h>
#include <cps/EMT/EMT_Ph3_Inductor.h>
#include <cps/EMT/EMT_Ph3_Resistor.h>
......@@ -23,6 +24,7 @@ namespace CPS {
class RXLoad :
public SimPowerComp<Real>,
public MNAInterface,
public DAEInterface,
public SharedFactory<RXLoad> {
protected:
/// Power [Watt]
......@@ -127,6 +129,18 @@ namespace CPS {
RXLoad& mLoad;
Attribute<Matrix>::Ptr mLeftVector;
};
//TODO
// #### DAE Section ####
///
void daeInitialize(double time, double state[], double dstate_dt[], int& offset);
/// Residual function for DAE Solver
void daeResidual(double time, const double state[], const double dstate_dt[], double resid[], std::vector<int>& off);
///
void daePostStep(const double state[], const double dstate_dt[], int& offset);
///
int getNumberOfStateVariables() {return 3;}
Please register or sign in to reply
};
}
}
......
......@@ -93,7 +93,8 @@ namespace CPS {
///
void setMatrixNodeIndex(UInt phase, UInt matrixNodeIndex) { mMatrixNodeIndex[phase] = matrixNodeIndex; }
///
void setVoltage(VarType newVoltage) { }
//void setVoltage(VarType newVoltage) { }
void setVoltage(VarType newVoltage, PhaseType phaseType = PhaseType::Single) { }
// #### MNA Section ####
///
......@@ -147,8 +148,8 @@ namespace CPS {
void SimNode<Complex>::mnaInitializeHarm(std::vector<Attribute<Matrix>::Ptr> leftVector);
template<>
void SimNode<Complex>::setVoltage(Complex newVoltage);
void SimNode<Complex>::setVoltage(Complex newVoltage, PhaseType phaseType);
template<>
void SimNode<Real>::setVoltage(double newVoltage);
void SimNode<Real>::setVoltage(double newVoltage, PhaseType phaseType);
}
......@@ -25,6 +25,9 @@ void EMT::Ph1::VoltageSource::setParameters(Complex voltageRef, Real srcFreq) {
attribute<Complex>("V_ref")->set(voltageRef);
attribute<Real>("f_src")->set(srcFreq);
mVoltageRef = attribute<Complex>("V_ref");
mSrcFreq = attribute<Real>("f_src");
parametersSet = true;
}
......@@ -38,8 +41,8 @@ void EMT::Ph1::VoltageSource::mnaInitialize(Real omega, Real timeStep, Attribute
MNAInterface::mnaInitialize(omega, timeStep);
updateMatrixNodeIndices();
mVoltageRef = attribute<Complex>("V_ref");
mSrcFreq = attribute<Real>("f_src");
// mVoltageRef = attribute<Complex>("V_ref");
// mSrcFreq = attribute<Real>("f_src");
mIntfVoltage(0,0) = Math::abs(mVoltageRef->get()) * cos(Math::phase(mVoltageRef->get()));
mMnaTasks.push_back(std::make_shared<MnaPreStep>(*this));
mMnaTasks.push_back(std::make_shared<MnaPostStep>(*this, leftVector));
......@@ -133,7 +136,8 @@ void EMT::Ph1::VoltageSource::daeResidual(double sim_time,
// resid[Pos1] = nodal current equation of node matrixNodeIndex(0)
// resid[Pos2] = nodal current equation of node matrixNodeIndex(1)
this->daeUpdateVoltage(sim_time);
//this->daeUpdateVoltage(sim_time);
this->updateVoltage(sim_time);
int Pos1 = matrixNodeIndex(0);
int Pos2 = matrixNodeIndex(1);
int c_offset = off[0]+off[1]; //current offset for component
......
......@@ -30,7 +30,6 @@ SimPowerComp<Real>::Ptr EMT::Ph3::NetworkInjection::clone(String name) {
return copy;
}
void EMT::Ph3::NetworkInjection::setParameters(MatrixComp voltageRef, Real srcFreq) {
attribute<MatrixComp>("V_ref")->set(voltageRef);
attribute<Real>("f_src")->set(srcFreq);
......@@ -172,32 +171,72 @@ void EMT::Ph3::NetworkInjection::mnaUpdateCurrent(const Matrix& leftVector) {
mIntfCurrent(2, 0) = Math::realFromVectorElement(leftVector, mVirtualNodes[0]->matrixNodeIndex(PhaseType::C));
}
//void EMT::Ph3::NetworkInjection::daeResidual(double ttime, const double state[], const double dstate_dt[], double resid[], std::vector<int>& off) {
// /* new state vector definintion:
// state[0]=node0_voltage
// state[1]=node1_voltage
// ....
// state[n]=noden_voltage
// state[n+1]=component0_voltage
// state[n+2]=component0_inductance (not yet implemented)
// ...
// state[m-1]=componentm_voltage
// state[m]=componentm_inductance
// */
//
// int Pos1 = matrixNodeIndex(0);
// int Pos2 = matrixNodeIndex(1);
// int c_offset = off[0] + off[1]; //current offset for component
// int n_offset_1 = c_offset + Pos1 + 1;// current offset for first nodal equation
// int n_offset_2 = c_offset + Pos2 + 1;// current offset for second nodal equation
// resid[c_offset] = (state[Pos2] - state[Pos1]) - state[c_offset]; // Voltage equation for Resistor
// //resid[++c_offset] = ; //TODO : add inductance equation
// resid[n_offset_1] += mIntfCurrent(0, 0).real();
// resid[n_offset_2] += mIntfCurrent(0, 0).real();
// off[1] += 1;
//}
//
//Complex EMT::Ph3::NetworkInjection::daeInitialize() {
// mIntfVoltage(0, 0) = mVoltageRef->get();
// return mVoltageRef->get();
//}
// #### DAE functions ####
void EMT::Ph3::NetworkInjection::daeInitialize(double time, double state[],
double dstate_dt[], int& offset) {
// offset: number of component in state, dstate_dt
// state[offset] = current through voltage source flowing into node matrixNodeIndex(1)
// dstate_dt[offset] = derivative of current through voltage source (not used yed) --> set to zero
updateMatrixNodeIndices();
state[offset] = mIntfCurrent(0,0);
dstate_dt[offset] = 0.0;
state[offset+1] = mIntfCurrent(1,0);
dstate_dt[offset+1] = 0.0;
state[offset+2] = mIntfCurrent(2,0);
dstate_dt[offset+2] = 0.0;
mSLog->info(
"\n--- daeInitialize ---"
"\nAdded current phase1 of NetworkInjection '{:s}' to state vector, initial value={:f}A"
"\nAdded current phase2 of NetworkInjection '{:s}' to state vector, initial value={:f}A"
"\nAdded current phase3 of NetworkInjection '{:s}' to state vector, initial value={:f}A"
"\nAdded derivative of current phase1 of NetworkInjection '{:s}' to derivative state vector, initial value={:f}"
"\nAdded derivative of current phase2 of NetworkInjection '{:s}' to derivative state vector, initial value={:f}"
"\nAdded derivative of current phase3 of NetworkInjection '{:s}' to derivative state vector, initial value={:f}"
"\n--- daeInitialize finished ---",
this->name(), state[offset],
this->name(), state[offset+1],
this->name(), state[offset+2],
this->name(), dstate_dt[offset],
this->name(), dstate_dt[offset+1],
this->name(), dstate_dt[offset+2]
);
mSLog->flush();
offset+=3;
}
void EMT::Ph3::NetworkInjection::daeResidual(double sim_time,
const double state[], const double dstate_dt[],
double resid[], std::vector<int>& off) {
// state[c_offset] = current through VoltageSource flowing into node matrixNodeIndex(1)
// dstate_dt[offset] = derivative of current through voltage source (not used yed)
// resid[c_offset] = v2-v1-v_s = state[Pos2] - state[Pos1] - mIntfVoltage(0,0)
// resid[Pos1] = nodal current equation of node matrixNodeIndex(0)
// resid[Pos2] = nodal current equation of node matrixNodeIndex(1)
this->updateVoltage(sim_time);
// int pos1 = matrixNodeIndex(PhaseType::A);
// int pos2 = matrixNodeIndex(PhaseType::B);
// int pos3 = matrixNodeIndex(PhaseType::C);
int c_offset = off[0]+off[1]; //current offset for component
resid[c_offset] = -mIntfVoltage(0,0);
resid[c_offset+1] = -mIntfVoltage(1,0);
resid[c_offset+2] = -mIntfVoltage(2,0);
// update nodal equations
resid[matrixNodeIndex(0, 0)] += state[c_offset];
resid[matrixNodeIndex(0, 1)] += state[c_offset+1];
resid[matrixNodeIndex(0, 2)] += state[c_offset+2];
off[1]+=3;
}
void EMT::Ph3::NetworkInjection::daePostStep(const double state[], const double dstate_dt[], int& offset) {
mIntfCurrent(0,0) = state[offset++];
mIntfCurrent(1,0) = state[offset++];
mIntfCurrent(2,0) = state[offset++];
}
\ No newline at end of file
......@@ -140,7 +140,8 @@ void EMT::Ph3::PiLine::initializeFromPowerflow(Real frequency) {
Logger::matrixToString(mIntfCurrent),
Logger::phasorToString(initialSingleVoltage(0)),
Logger::phasorToString(initialSingleVoltage(1)),
Logger::phasorToString(mVirtualNodes[0]->initialSingleVoltage()));
Logger::phasorToString(mVirtualNodes[0]->initialSingleVoltage())
);
}
void EMT::Ph3::PiLine::mnaInitialize(Real omega, Real timeStep, Attribute<Matrix>::Ptr leftVector) {
......@@ -257,3 +258,123 @@ void EMT::Ph3::PiLine::mnaUpdateCurrent(const Matrix& leftVector) {
//void EMT::Ph3::PiLine::mnaTearPostStep(Complex voltage, Complex current) {
// mSubSeriesInductor->mnaTearPostStep(voltage - current * mSeriesRes, current);
//}
// #### DAE functions ####
void EMT::Ph3::PiLine::daeInitialize(double time, double state[],
double dstate_dt[], int& offset) {
// offset: number of component in state, dstate_dt
// state[offset] = current through voltage source flowing into node matrixNodeIndex(1)
// dstate_dt[offset] = derivative of current through voltage source (not used yed) --> set to zero
updateMatrixNodeIndices();
// add seriesInductors to state vector
Matrix inductorCurrent = mSubSeriesInductor->intfCurrent();
Matrix inductorVoltage = mSubSeriesInductor->intfVoltage();
state[offset] = inductorCurrent(0,0);
dstate_dt[offset] = inductorVoltage(0,0)/mSeriesInd(0, 0);
state[offset+1] = inductorCurrent(1,0);
dstate_dt[offset+1] = inductorVoltage(1,0)/mSeriesInd(1, 1);
state[offset+2] = inductorCurrent(2,0);
dstate_dt[offset+2] = inductorVoltage(2,0)/mSeriesInd(2, 2);
// Modify initial value of derivative of node voltages
Matrix capacitor0Current = mSubParallelCapacitor0->intfCurrent();
Matrix capacitor1Current = mSubParallelCapacitor1->intfCurrent();
if (mParallelCond(0,0) > 0) {
if (terminalNotGrounded(1)) {
dstate_dt[matrixNodeIndex(1, 0)] = capacitor0Current(0,0)/mParallelCap(0,0);
dstate_dt[matrixNodeIndex(1, 1)] = capacitor0Current(1,0)/mParallelCap(1,1);
dstate_dt[matrixNodeIndex(1, 2)] = capacitor0Current(2,0)/mParallelCap(2,2);
}
if (terminalNotGrounded(0)) {
dstate_dt[matrixNodeIndex(0, 0)] = capacitor1Current(0,0)/mParallelCap(0,0);
dstate_dt[matrixNodeIndex(0, 1)] = capacitor1Current(1,0)/mParallelCap(1,1);
dstate_dt[matrixNodeIndex(0, 2)] = capacitor1Current(2,0)/mParallelCap(2,2);
}
}
mSLog->info(
"\n--- daeInitialize ---"
"\nAdded current-phase1 through the inductor of PiLine '{:s}' to state vector, initial value={:f}A"
"\nAdded current-phase2 through the inductor of PiLine '{:s}' to state vector, initial value={:f}A"
"\nAdded current-phase3 through the inductor of PiLine '{:s}' to state vector, initial value={:f}A"
"\nAdded derivative of current-phase1 through the inductor of PiLine '{:s}' to derivative state vector, initial value={:f}"
"\nAdded derivative of current-phase2 through the inductor of PiLine '{:s}' to derivative state vector, initial value={:f}"
"\nAdded derivative of current-phase3 through the inductor of PiLine '{:s}' to derivative state vector, initial value={:f}"
"\nSet initial value of derivative of node voltage of capacitor0-phaseA of PiLine '{:s}' ={:f}"
"\nSet initial value of derivative of node voltage of capacitor0-phaseB of PiLine '{:s}' ={:f}"
"\nSet initial value of derivative of node voltage of capacitor0-phaseC of PiLine '{:s}' ={:f}"
"\nSet initial value of derivative of node voltage of capacitor1-phaseA of PiLine '{:s}' ={:f}"
"\nSet initial value of derivative of node voltage of capacitor1-phaseB of PiLine '{:s}' ={:f}"
"\nSet initial value of derivative of node voltage of capacitor1-phaseC of PiLine '{:s}' ={:f}"
"\n--- daeInitialize finished ---",
this->name(), state[offset],
this->name(), state[offset+1],
this->name(), state[offset+2],
this->name(), dstate_dt[offset],
this->name(), dstate_dt[offset+1],
this->name(), dstate_dt[offset+2],
this->name(), dstate_dt[matrixNodeIndex(0, 0)],
this->name(), dstate_dt[matrixNodeIndex(0, 1)],
this->name(), dstate_dt[matrixNodeIndex(0, 2)],
this->name(), dstate_dt[matrixNodeIndex(1, 0)],
this->name(), dstate_dt[matrixNodeIndex(1, 1)],
this->name(), dstate_dt[matrixNodeIndex(1, 2)]
);
std::cout << "test\n";
mSLog->flush();
offset+=3;
}
void EMT::Ph3::PiLine::daeResidual(double sim_time,
const double state[], const double dstate_dt[],
double resid[], std::vector<int>& off) {
/*
R*i_rx + L*der(i_rx) = v1 - v2;
i1 = i_rx + G/2*v1 + C/2*der(v1);
-i2 = -i_rx + G/2*v2 + C/2*der(v2);
*/
int c_offset = off[0]+off[1]; //current offset for component
resid[c_offset] = -mSeriesInd(0, 0)*dstate_dt[c_offset] - state[c_offset]*mSeriesRes(0,0);
resid[c_offset+1] = -mSeriesInd(1, 1)*dstate_dt[c_offset+1] - state[c_offset]*mSeriesRes(1,1);
resid[c_offset+2] = -mSeriesInd(2, 2)*dstate_dt[c_offset+2] - state[c_offset]*mSeriesRes(2,2);
if (terminalNotGrounded(0)) {
resid[c_offset] -= state[matrixNodeIndex(0, 0)];
resid[c_offset+1] -= state[matrixNodeIndex(0, 1)];
resid[c_offset+2] -= state[matrixNodeIndex(0, 2)];
resid[matrixNodeIndex(0, 0)] -= state[c_offset] ;
resid[matrixNodeIndex(0, 1)] -= state[c_offset+1];
resid[matrixNodeIndex(0, 2)] -= state[c_offset+2];
resid[matrixNodeIndex(0, 0)] += -state[c_offset] + mParallelCap(0,0)*dstate_dt[matrixNodeIndex(0, 0)]/2 + state[matrixNodeIndex(0, 0)]*mParallelCond(0,0)/2;
resid[matrixNodeIndex(0, 1)] += -state[c_offset+1] + mParallelCap(1,1)*dstate_dt[matrixNodeIndex(0, 1)]/2 + state[matrixNodeIndex(0, 1)]*mParallelCond(1,1)/2;
resid[matrixNodeIndex(0, 2)] += -state[c_offset+2] + mParallelCap(2,2)*dstate_dt[matrixNodeIndex(0, 2)]/2 + state[matrixNodeIndex(0, 2)]*mParallelCond(2,2)/2;
}
if (terminalNotGrounded(1)) {
resid[c_offset] += state[matrixNodeIndex(1, 0)];
resid[c_offset+1] += state[matrixNodeIndex(1, 1)];
resid[c_offset+1] += state[matrixNodeIndex(1, 2)];
resid[matrixNodeIndex(1, 0)] += state[c_offset] + mParallelCap(0,0)*dstate_dt[matrixNodeIndex(1, 0)]/2 + state[matrixNodeIndex(1, 0)]*mParallelCond(0,0)/2;
resid[matrixNodeIndex(1, 1)] += state[c_offset+1] + mParallelCap(1,1)*dstate_dt[matrixNodeIndex(1, 1)]/2 + state[matrixNodeIndex(1, 1)]*mParallelCond(1,1)/2;
resid[matrixNodeIndex(1, 2)] += state[c_offset+2] + mParallelCap(2,2)*dstate_dt[matrixNodeIndex(1, 2)]/2 + state[matrixNodeIndex(1, 2)]*mParallelCond(2,2)/2;
}
off[1] += 3;
}
void EMT::Ph3::PiLine::daePostStep(const double state[], const double dstate_dt[], int& offset) {
mIntfVoltage(0,0) = state[matrixNodeIndex(1, 0)] - state[matrixNodeIndex(0, 0)];
mIntfVoltage(1,0) = state[matrixNodeIndex(1, 1)] - state[matrixNodeIndex(0, 1)];
mIntfVoltage(2,0) = state[matrixNodeIndex(1, 2)] - state[matrixNodeIndex(0, 2)];
mIntfCurrent(0,0) = state[offset++];
mIntfCurrent(1,0) = state[offset++];
mIntfCurrent(2,0) = state[offset++];
}
......@@ -215,3 +215,125 @@ void EMT::Ph3::RXLoad::mnaUpdateCurrent(const Matrix& leftVector) {
if (mSubCapacitor)
mIntfCurrent += mSubCapacitor->intfCurrent();
}
// #### DAE functions ####
void EMT::Ph3::RXLoad::daeInitialize(double time, double state[],
double dstate_dt[], int& offset) {
// offset: number of component in state, dstate_dt
// state[offset] = current through voltage source flowing into node matrixNodeIndex(1)
// dstate_dt[offset] = derivative of current through voltage source (not used yed) --> set to zero
updateMatrixNodeIndices();
if (mReactance(0,0) > 0) {
// state variable is the inductor current
Matrix inductorCurrent = mSubInductor->intfCurrent();
Matrix inductorVoltage = mSubInductor->intfVoltage();
state[offset] = inductorCurrent(0,0);
dstate_dt[offset] = inductorVoltage(0,0)/mInductance(0,0);
state[offset+1] = inductorCurrent(1,0);
dstate_dt[offset+1] = inductorVoltage(1,0)/mInductance(1,1);
state[offset+2] = inductorCurrent(2,0);
dstate_dt[offset+2] = inductorVoltage(2,0)/mInductance(2, 2);
mSLog->info(
"\n--- daeInitialize ---"
"\nmReactance(0,0) > 0 --> state variable are inductor currents"
"\nAdded current-phase1 through the inductor of RXLoad '{:s}' to state vector, initial value={:f}A"
"\nAdded current-phase2 through the inductor of RXLoad '{:s}' to state vector, initial value={:f}A"
"\nAdded current-phase3 through the inductor of RXLoad '{:s}' to state vector, initial value={:f}A"
"\nAdded derivative of current-phase1 through the inductor of RXLoad '{:s}' to derivative state vector, initial value={:f}"
"\nAdded derivative of current-phase2 through the inductor of RXLoad '{:s}' to derivative state vector, initial value={:f}"
"\nAdded derivative of current-phase3 through the inductor of RXLoad '{:s}' to derivative state vector, initial value={:f}"
"\n--- daeInitialize finished ---",
this->name(), state[offset],
this->name(), state[offset+1],
this->name(), state[offset+2],
this->name(), dstate_dt[offset],
this->name(), dstate_dt[offset+1],
this->name(), dstate_dt[offset+2]
);
}
else if (mReactance(0,0) < 0) {
// state variable is the voltage through capacitor
Matrix capacitorCurrent = mSubCapacitor->intfCurrent();
Matrix capacitorVoltage = mSubCapacitor->intfVoltage();
state[offset] = capacitorVoltage(0,0);
dstate_dt[offset] = capacitorCurrent(0,0)/mCapacitance(0,0);
state[offset+1] = capacitorVoltage(1,0);
dstate_dt[offset+1] = capacitorCurrent(1,0)/mCapacitance(1,1);
state[offset+2] = capacitorVoltage(2,0);
dstate_dt[offset+2] = capacitorCurrent(2,0)/mCapacitance(2, 2);
mSLog->info(
"\n--- daeInitialize ---"
"\nmReactance(0,0) < 0 --> state variable are capacitor currentvoltages"
"\nAdded voltage-phase1 through the capacitor of RXLoad '{:s}' to state vector, initial value={:f}A"
"\nAdded voltage-phase2 through the capacitor of RXLoad '{:s}' to state vector, initial value={:f}A"
"\nAdded voltage-phase3 through the capacitor of RXLoad '{:s}' to state vector, initial value={:f}A"
"\nAdded derivative of voltage-phase1 through the capacitor of RXLoad '{:s}' to derivative state vector, initial value={:f}"
"\nAdded derivative of voltage-phase2 through the capacitor of RXLoad '{:s}' to derivative state vector, initial value={:f}"
"\nAdded derivative of voltage-phase3 through the capacitor of RXLoad '{:s}' to derivative state vector, initial value={:f}"
"\n--- daeInitialize finished ---",
this->name(), state[offset],
this->name(), state[offset+1],
this->name(), state[offset+2],
this->name(), dstate_dt[offset],
this->name(), dstate_dt[offset+1],