Commit a4c4690f authored by Markus Mirz's avatar Markus Mirz
Browse files

Simulation: matrix recompute solver and refactor

parent 5731db15
......@@ -62,7 +62,6 @@ int main(int argc, char *argv[]) {
sim.setTimeStep(0.0001);
sim.setFinalTime(0.1);
sim.doSteadyStateInit(true);
sim.doHarmonicParallelization(false);
sim.addLogger(logger);
sim.run();
......
......@@ -65,7 +65,6 @@ int main(int argc, char *argv[]) {
sim.setTimeStep(0.0001);
sim.setFinalTime(2);
sim.doSteadyStateInit(true);
sim.doHarmonicParallelization(false);
sim.addLogger(logger);
sim.run();
......
......@@ -89,7 +89,7 @@ int main(int argc, char* argv[]) {
sim.setSystem(sys);
sim.setTimeStep(timeStep);
sim.setFinalTime(finalTime);
sim.doHarmonicParallelization(false);
sim.doFrequencyParallelization(false);
// Logging
auto logger = DataLogger::make(simName);
......
......@@ -97,7 +97,7 @@ int main(int argc, char* argv[]) {
sim.setSystem(sys);
sim.setTimeStep(timeStep);
sim.setFinalTime(finalTime);
sim.doHarmonicParallelization(true);
sim.doFrequencyParallelization(true);
if (threads > 0) {
auto scheduler = std::make_shared<ThreadLevelScheduler>(threads);
sim.setScheduler(scheduler);
......
......@@ -95,7 +95,7 @@ int main(int argc, char* argv[]) {
sim.setSystem(sys);
sim.setTimeStep(timeStep);
sim.setFinalTime(finalTime);
sim.doHarmonicParallelization(true);
sim.doFrequencyParallelization(true);
sim.run();
sim.logStepTimes(simName + "_step_times");
......
......@@ -40,8 +40,7 @@ namespace DPsim {
/// \brief The Simulation holds a SystemTopology and a Solver.
///
/// Every time step, the Simulation calls the step function of the Solver.
class Simulation :
public CPS::AttributeList {
class Simulation : public CPS::AttributeList {
public:
typedef std::shared_ptr<Simulation> Ptr;
......@@ -61,7 +60,6 @@ namespace DPsim {
/// System list
CPS::SystemTopology mSystem;
// #### Logging ####
/// Simulation log level
CPS::Logger::Level mLogLevel;
......@@ -77,6 +75,8 @@ namespace DPsim {
Solver::List mSolvers;
///
Bool mPowerFlowInit = false;
/// Enable recomputation of system matrix during simulation
Bool mSystemMatrixRecomputation = false;
/// Determines if the network should be split
/// into subnetworks at decoupling lines.
......@@ -91,11 +91,11 @@ namespace DPsim {
/// This can only be done if the network is composed
/// of linear components that do no create cross
/// frequency coupling.
Bool mHarmParallel = false;
Bool mFreqParallel = false;
///
Bool mInitialized = false;
// #### steady state initialization ####
// #### Initialization ####
/// steady state initialization time limit
Real mSteadStIniTimeLimit = 10;
/// steady state initialization accuracy limit
......@@ -141,10 +141,15 @@ namespace DPsim {
CPS::Domain domain = CPS::Domain::DP,
CPS::Logger::Level logLevel = CPS::Logger::Level::info);
/// Create solvers depending on simulation settings
template <typename VarType>
void createSolvers(CPS::SystemTopology& system, CPS::IdentifiedObject::List& tearComponents);
void createSolvers();
/// Subroutine for MNA only because there are many MNA options
template <typename VarType>
void createMNASolver();
/// Prepare schedule for simulation
void prepSchedule();
public:
/// Simulation logger
CPS::Logger::Log mLog;
......@@ -189,8 +194,12 @@ namespace DPsim {
void setScheduler(std::shared_ptr<Scheduler> scheduler) {
mScheduler = scheduler;
}
/// Compute phasors of different frequencies in parallel
void doFrequencyParallelization(Bool value) { mFreqParallel = value; }
///
void doSystemMatrixRecomputation(Bool value) { mSystemMatrixRecomputation = value; }
// #### steady state initialization ####
// #### Initialization ####
/// activate steady state initialization
void doSteadyStateInit(Bool f) { mSteadyStateInit = f; }
/// set steady state initialization time limit
......@@ -198,9 +207,6 @@ namespace DPsim {
/// set steady state initialization accuracy limit
void setSteadStIniAccLimit(Real v) { mSteadStIniAccLimit = v; }
///
void doHarmonicParallelization(Bool parallel) { mHarmParallel = parallel; }
// #### Simulation Control ####
/// Create solver instances etc.
void initialize();
......
......@@ -24,6 +24,7 @@
#include <dpsim/Utils.h>
#include <cps/Utils.h>
#include <dpsim/MNASolver.h>
#include <dpsim/MNASolverSysRecomp.h>
#include <dpsim/PFSolverPowerPolar.h>
#include <dpsim/DiakopticsSolver.h>
......@@ -94,14 +95,13 @@ void Simulation::initialize() {
mSolvers.clear();
switch (mDomain) {
case Domain::SP:
// Treat SP as DP
case Domain::DP:
createSolvers<Complex>(mSystem, mTearComponents);
createSolvers<Complex>();
break;
case Domain::EMT:
createSolvers<Real>(mSystem, mTearComponents);
break;
case Domain::SP:
createSolvers<Complex>(mSystem, mTearComponents);
createSolvers<Real>();
break;
}
......@@ -114,77 +114,96 @@ void Simulation::initialize() {
}
template <typename VarType>
void Simulation::createSolvers(
CPS::SystemTopology& system,
IdentifiedObject::List& tearComponents) {
void Simulation::createSolvers() {
Solver::Ptr solver;
switch (mSolverType) {
case Solver::Type::MNA:
createMNASolver<VarType>();
break;
#ifdef WITH_SUNDIALS
case Solver::Type::DAE:
solver = std::make_shared<DAESolver>(mName, mSystem, mTimeStep, 0.0);
mSolvers.push_back(solver);
break;
#endif /* WITH_SUNDIALS */
case Solver::Type::NRP:
solver = std::make_shared<PFSolverPowerPolar>(mName, mSystem, mTimeStep, mLogLevel);
solver->doPowerFlowInit(mPowerFlowInit);
solver->initialize();
mSolvers.push_back(solver);
break;
default:
throw UnsupportedSolverException();
}
// Some components require a dedicated ODE solver.
// This solver is independent of the system solver.
#ifdef WITH_SUNDIALS
for (auto comp : mSystem.mComponents) {
auto odeComp = std::dynamic_pointer_cast<ODEInterface>(comp);
if (odeComp) {
// TODO explicit / implicit integration
auto odeSolver = std::make_shared<ODESolver>(
odeComp->attribute<String>("name")->get() + "_ODE", odeComp, false, mTimeStep);
mSolvers.push_back(odeSolver);
}
}
#endif /* WITH_SUNDIALS */
}
template <typename VarType>
void Simulation::createMNASolver() {
Solver::Ptr solver;
std::vector<SystemTopology> subnets;
// The Diakoptics solver splits the system at a later point.
// That is why the system is not split here if tear components exist.
if (mSplitSubnets && tearComponents.size() == 0)
system.splitSubnets<VarType>(subnets);
if (mSplitSubnets && mTearComponents.size() == 0)
mSystem.splitSubnets<VarType>(subnets);
else
subnets.push_back(system);
subnets.push_back(mSystem);
for (UInt net = 0; net < subnets.size(); net++) {
String copySuffix;
if (subnets.size() > 1)
copySuffix = "_" + std::to_string(net);
// TODO in the future, here we could possibly even use different
// TODO: In the future, here we could possibly even use different
// solvers for different subnets if deemed useful
Solver::Ptr solver;
switch (mSolverType) {
case Solver::Type::MNA:
if (tearComponents.size() > 0) {
solver = std::make_shared<DiakopticsSolver<VarType>>(mName,
subnets[net], tearComponents, mTimeStep, mLogLevel);
} else {
if (mTearComponents.size() > 0) {
// Tear components available, use diakoptics
solver = std::make_shared<DiakopticsSolver<VarType>>(mName,
subnets[net], mTearComponents, mTimeStep, mLogLevel);
}
else if (mSystemMatrixRecomputation) {
// Recompute system matrix if switches or other components change
solver = std::make_shared<MnaSolverSysRecomp<VarType>>(
mName + copySuffix, mDomain, mLogLevel);
solver->setTimeStep(mTimeStep);
solver->doSteadyStateInit(mSteadyStateInit);
solver->setSteadStIniTimeLimit(mSteadStIniTimeLimit);
solver->setSteadStIniAccLimit(mSteadStIniAccLimit);
solver->setSystem(subnets[net]);
solver->initialize();
}
else {
// Default case with precomputed system matrices for different configurations
#ifdef WITH_CUDA
solver = std::make_shared<MnaSolverGpu<VarType>>(
mName + copySuffix, mDomain, mLogLevel);
solver = std::make_shared<MnaSolverGpu<VarType>>(
mName + copySuffix, mDomain, mLogLevel);
#else
solver = std::make_shared<MnaSolver<VarType>>(
mName + copySuffix, mDomain, mLogLevel);
solver = std::make_shared<MnaSolver<VarType>>(
mName + copySuffix, mDomain, mLogLevel);
#endif /* WITH_CUDA */
solver->setTimeStep(mTimeStep);
solver->doSteadyStateInit(mSteadyStateInit);
solver->doFrequencyParallelization(mHarmParallel);
solver->setSteadStIniTimeLimit(mSteadStIniTimeLimit);
solver->setSteadStIniAccLimit(mSteadStIniAccLimit);
solver->setSystem(subnets[net]);
solver->initialize();
}
break;
#ifdef WITH_SUNDIALS
case Solver::Type::DAE:
solver = std::make_shared<DAESolver>(mName + copySuffix, subnets[net], mTimeStep, 0.0);
break;
#endif /* WITH_SUNDIALS */
case Solver::Type::NRP:
solver = std::make_shared<PFSolverPowerPolar>(mName, mSystem, mTimeStep, mLogLevel);
solver->doPowerFlowInit(mPowerFlowInit);
solver->initialize();
break;
default:
throw UnsupportedSolverException();
solver->setTimeStep(mTimeStep);
solver->doSteadyStateInit(mSteadyStateInit);
solver->doFrequencyParallelization(mFreqParallel);
solver->setSteadStIniTimeLimit(mSteadStIniTimeLimit);
solver->setSteadStIniAccLimit(mSteadStIniAccLimit);
solver->setSystem(subnets[net]);
solver->initialize();
}
mSolvers.push_back(solver);
}
// Some components require a dedicated ODE solver.
// This solver is independet of the system solver.
#ifdef WITH_SUNDIALS
for (auto comp : system.mComponents) {
auto odeComp = std::dynamic_pointer_cast<ODEInterface>(comp);
if (odeComp) {
// TODO explicit / implicit integration
auto odeSolver = std::make_shared<ODESolver>(
odeComp->attribute<String>("name")->get() + "_ODE", odeComp, false, mTimeStep);
mSolvers.push_back(odeSolver);
}
}
#endif /* WITH_SUNDIALS */
}
void Simulation::sync() {
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment