Skip to content

Commit

Permalink
Enhance DASSL in C++ runtime to support basic DAE mode and add test
Browse files Browse the repository at this point in the history 
DASSL is the default solver, whereas IDA fails for the simple test.
So far not supplying DAE Jacobian.
  • Loading branch information
@rfranke
rfranke committed Apr 20, 2023
1 parent f77c006 commit 5d8d697
Show file tree
Hide file tree
Showing 5 changed files with 100 additions and 28 deletions.
10 changes: 4 additions & 6 deletions OMCompiler/Compiler/Template/CodegenCppOld.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1528,9 +1528,7 @@ case SIMCODE(modelInfo=MODELINFO(__)) then
end functionStateSets;




template simulationMainRunScript(SimCode simCode ,Text& extraFuncs,Text& extraFuncsDecl,Text extraFuncsNamespace, String preRunCommandLinux, String preRunCommandWindows, String execCommandLinux)
template simulationMainRunScript(SimCode simCode, Text& extraFuncs, Text& extraFuncsDecl, Text extraFuncsNamespace, String preRunCommandLinux, String preRunCommandWindows, String execCommandLinux)
"Generates code for header file for simulation target."
::=
match simCode
Expand All @@ -1540,15 +1538,15 @@ template simulationMainRunScript(SimCode simCode ,Text& extraFuncs,Text& extraFu
let stepsize = settings.stepSize
let intervals = settings.numberOfIntervals
let tol = settings.tolerance
let solver = match simCode case SIMCODE(daeModeData=NONE()) then settings.method else 'ida' //for dae mode only ida is supported
let moLib = makefileParams.compileDir
let solver = settings.method
let moLib = makefileParams.compileDir
let home = makefileParams.omhome
let outputformat = settings.outputFormat
let execParameters = '-S <%start%> -E <%end%> -H <%stepsize%> -G <%intervals%> -P <%outputformat%> -T <%tol%> -I <%solver%> -R "<%simulationLibDir(simulationCodeTarget(),simCode , &extraFuncs , &extraFuncsDecl, extraFuncsNamespace)%>" -M "<%moLib%>" -r "<%simulationResults(Testsuite.isRunning(),simCode , &extraFuncs , &extraFuncsDecl, extraFuncsNamespace)%>"'
let outputParameter = if (stringEq(settings.outputFormat, "empty")) then "-O none" else ""
let fileNamePrefixx = fileNamePrefix

let libFolder = simulationLibDir(simulationCodeTarget(),simCode , &extraFuncs , &extraFuncsDecl, extraFuncsNamespace)
let libFolder = simulationLibDir(simulationCodeTarget(), simCode, &extraFuncs, &extraFuncsDecl, extraFuncsNamespace)
let libPaths = makefileParams.libPaths |> path => path; separator=";"

match makefileParams.platform
Expand Down
71 changes: 50 additions & 21 deletions OMCompiler/SimulationRuntime/cpp/Solver/DASSL/DASSL.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -180,7 +180,9 @@ void DASSL::initialize()
_tLastEvent = 0.0;
_event_n = 0;
SolverDefaultImplementation::initialize();
_dimSys = _continuous_system->getDimContinuousStates();
_dimStates = _continuous_system->getDimContinuousStates();
_dimAE = _continuous_system->getDimAE();
_dimSys = _dimStates + _dimAE;
_dimZeroFunc = _event_system->getDimZeroFunc();

if (_dimSys == 0)
Expand Down Expand Up @@ -243,11 +245,19 @@ void DASSL::initialize()

// Set tolerances
_info[1] = 1;
_atol[0] = _rtol[0] = 1.0; // in case of dummy state
_continuous_system->getNominalStates(_atol);
for (int i = 0; i < _dimSys; i++) {
_atol[i] = max(_atol[i] * _settings->getATol(), 1e-10);
_rtol[i] = max(_settings->getRTol(), 1e-10);
if (_dimAE == 0) {
_atol[0] = _rtol[0] = 1.0; // in case of dummy state
_continuous_system->getNominalStates(_atol);
for (int i = 0; i < _dimStates; i++) {
_atol[i] = max(_atol[i] * _settings->getATol(), 1e-10);
_rtol[i] = max(_settings->getRTol(), 1e-10);
}
}
else {
for (int i = 0; i < _dimSys; i++) {
_atol[i] = _settings->getATol();
_rtol[i] = _settings->getRTol();
}
}
LOGGER_WRITE_VECTOR("atol", _atol, _dimSys, LC_SOLVER, LL_DEBUG);
LOGGER_WRITE_VECTOR("rtol", _rtol, _dimSys, LC_SOLVER, LL_DEBUG);
Expand All @@ -256,7 +266,7 @@ void DASSL::initialize()
_info[2] = 1;

// Use supplied Jacobian function
_info[4] = 1;
_info[4] = _dimAE == 0? 1: 0;
_maxColors = _system->getAMaxColors();
if (_system->isAnalyticJacobianGenerated() && _continuous_system->getDimContinuousStates() > 0)
{
Expand Down Expand Up @@ -436,7 +446,13 @@ void DASSL::DASSLCore()
if (_idid == 3)
_time_system->setTime(_tEnd); // interpolated time point
_continuous_system->setContinuousStates(_y);
_continuous_system->evaluateAll(IContinuous::CONTINUOUS);
if (_dimAE == 0)
_continuous_system->evaluateAll(IContinuous::CONTINUOUS);
else {
_mixed_system->setAlgebraicDAEVars(_y + _dimStates);
_continuous_system->setStateDerivatives(_yPrime);
_continuous_system->evaluateDAE(IContinuous::CONTINUOUS);
}
writeToFile(_accStps, _tCurrent, _h);
}

Expand Down Expand Up @@ -503,7 +519,10 @@ void DASSL::DASSLCore()
{
try
{
_continuous_system->evaluateAll(IContinuous::CONTINUOUS);
if (_dimAE == 0)
_continuous_system->evaluateAll(IContinuous::CONTINUOUS);
else
_continuous_system->evaluateDAE(IContinuous::CONTINUOUS);
}
catch (std::exception& ex)
{
Expand All @@ -530,7 +549,10 @@ void DASSL::DASSLCore()
if (writeEventOutput)
{
// If we want to write the event-results, we should evaluate the whole system again
_continuous_system->evaluateAll(IContinuous::CONTINUOUS);
if (_dimAE == 0)
_continuous_system->evaluateAll(IContinuous::CONTINUOUS);
else
_continuous_system->evaluateDAE(IContinuous::CONTINUOUS);
writeToFile(_accStps, _tCurrent, _h);
}

Expand Down Expand Up @@ -575,16 +597,13 @@ bool DASSL::stateSelection()
int DASSL::_res(double *t, double *y, double *yp,
double *cj, double *delta, int *ires, double *rpar, int *ipar)
{
int success = ((DASSL *)ipar)->calcFunction(*t, y, delta);
int n = ((DASSL *)ipar)->_dimSys;
for (int i = 0; i < n; i++)
delta[i] -= yp[i];
int success = ((DASSL *)ipar)->calcFunction(*t, y, yp, delta);
if (!success)
*ires = -1;
return 0;
}

int DASSL::calcFunction(const double& time, const double* y, double* f)
int DASSL::calcFunction(const double& time, const double* y, const double *yp, double* f)
{
int success = 0;

Expand All @@ -601,8 +620,18 @@ int DASSL::calcFunction(const double& time, const double* y, double* f)
f[0] = 0.0; // in case of dummy state
_time_system->setTime(time);
_continuous_system->setContinuousStates(y);
_continuous_system->evaluateODE(IContinuous::CONTINUOUS);
_continuous_system->getRHS(f);
if (_dimAE == 0) {
_continuous_system->evaluateODE(IContinuous::CONTINUOUS);
_continuous_system->getRHS(f);
for (int i = 0; i < _dimStates; i++)
f[i] -= yp[i];
}
else {
_mixed_system->setAlgebraicDAEVars(y + _dimStates);
_continuous_system->setStateDerivatives(yp);
_continuous_system->evaluateDAE(IContinuous::CONTINUOUS);
_mixed_system->getResidual(f);
}
success = 1;
}
catch (std::exception & ex)
Expand Down Expand Up @@ -707,14 +736,14 @@ int DASSL::calcJacobian(double t, double *y, double *yp, double *delta,
_yJac[j] += _dyJac[j];
}

calcFunction(t, _yJac, _fJac);
calcFunction(t, _yJac, yp, _fJac);

for (int j: _system->getAColumnsOfColor(color))
{
int startOfColumn = j * _dimSys;
for (int i: _system->getADependenciesOfColumn(j))
{
pd[startOfColumn + i] = (_fJac[i] - delta[i] - yp[i]) / _dyJac[j];
pd[startOfColumn + i] = (_fJac[i] - delta[i]) / _dyJac[j];
}
_yJac[j] = y[j];
}
Expand All @@ -726,12 +755,12 @@ int DASSL::calcJacobian(double t, double *y, double *yp, double *delta,
{
_yJac[j] += _dyJac[j];

calcFunction(t, _yJac, _fJac);
calcFunction(t, _yJac, yp, _fJac);

int startOfColumn = j * _dimSys;
for (int i = 0; i < _dimSys; i++)
{
pd[startOfColumn + i] = (_fJac[i] - delta[i] - yp[i]) / _dyJac[j];
pd[startOfColumn + i] = (_fJac[i] - delta[i]) / _dyJac[j];
}

_yJac[j] = y[j];
Expand Down
6 changes: 5 additions & 1 deletion OMCompiler/SimulationRuntime/cpp/Solver/DASSL/DASSL.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ class DASSL
void DASSLCore();

/// evaluate right hand side of ODE
int calcFunction(const double& time, const double* y, double* yp);
int calcFunction(const double& time, const double* y, const double* yp, double *f);

// callback for right hand side
static int _res(double *t, double *y, double *yp, double *cj,
Expand Down Expand Up @@ -142,6 +142,10 @@ class DASSL
// Variables for Coloured Jacobians
int _maxColors;

// DAE support
int _dimStates;
int _dimAE;

ISystemProperties* _properties;
IContinuous* _continuous_system;
IEvent* _event_system;
Expand Down
1 change: 1 addition & 0 deletions testsuite/openmodelica/cppruntime/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ RefArrayDim2.mos \
solveTest.mos \
solveOneNonlinearEquationTest.mos \
testArrayEquations.mos \
testDAE.mos \
testMatrixIO.mos \
testMatrixState.mos \
testReduction.mos \
Expand Down
40 changes: 40 additions & 0 deletions testsuite/openmodelica/cppruntime/testDAE.mos
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,40 @@
// name: testDAE
// keywords: DAE mode
// status: correct
// teardown_command: rm -f *DAETest*

setCommandLineOptions("--simCodeTarget=Cpp");

loadString("
model DAETest
Real x1(start = -1);
Real x2;
Real x3;
equation
der(x1) = x3;
x2 * (1 - x2) = 0;
x1*x2 + x3*(1 - x2) = time;
annotation(experiment(StopTime = 1), __OpenModelica_commandLineOptions = \"--daeMode\");
end DAETest;
");
getErrorString();

simulate(DAETest);

val(x1, 1);
val(x2, 1);
val(x3, 1);

// Result:
// true
// true
// ""
// record SimulationResult
// resultFile = "DAETest_res.mat",
// simulationOptions = "startTime = 0.0, stopTime = 1.0, numberOfIntervals = 500, tolerance = 1e-06, method = 'dassl', fileNamePrefix = 'DAETest', options = '', outputFormat = 'mat', variableFilter = '.*', cflags = '', simflags = ''",
// messages = ""
// end SimulationResult;
// -0.5000000000000006
// 0.0
// 1.0
// endResult

0 comments on commit 5d8d697

Please sign in to comment.