Skip to content

Commit

Permalink
[DAEmode] improve event handling and discrete loops
Browse files Browse the repository at this point in the history 
Belonging to [master]:
  - OpenModelica/OMCompiler#2291
  - OpenModelica/OpenModelica-testsuite#884
  • Loading branch information
@OpenModelica-Hudson
Willi Braun authored and OpenModelica-Hudson committed Mar 19, 2018
1 parent 90be434 commit c2feb3f
Show file tree
Hide file tree
Showing 4 changed files with 113 additions and 42 deletions.
82 changes: 74 additions & 8 deletions Compiler/BackEnd/DAEMode.mo
View file
Expand Up @@ -46,6 +46,7 @@ import BackendDAE;
protected

import Absyn;
import Array;
import BackendDAEOptimize;
import BackendDAEUtil;
import BackendDAEFunc;
Expand All @@ -67,6 +68,7 @@ import Flags;
import Global;
import Initialization;
import List;
import Matching;
import StackOverflow;
import Util;

Expand Down Expand Up @@ -196,6 +198,7 @@ uniontype TraverseEqnAryFold
BackendDAE.Variables systemVars;
DAE.FunctionTree functionTree;
Boolean recursiveStrongComponentRun;
BackendDAE.Shared shared;
end TRAVERSER_CREATE_DAE;
end TraverseEqnAryFold;

Expand All @@ -221,7 +224,7 @@ algorithm
systemSize := BackendDAEUtil.systemSize(syst);
newDAEVars := BackendVariable.emptyVars();
newDAEEquations := BackendEquation.emptyEqnsSized(systemSize);
travArgs := TRAVERSER_CREATE_DAE(globalDAEData, newDAEVars, newDAEEquations, syst.orderedVars, shared.functionTree, false);
travArgs := TRAVERSER_CREATE_DAE(globalDAEData, newDAEVars, newDAEEquations, syst.orderedVars, shared.functionTree, false, shared);
if debug then BackendDump.printEqSystem(syst); end if;
// for every equation create corresponding residiual variable(s)
travArgs := BackendDAEUtil.traverseEqSystemStrongComponents(syst, traverserStrongComponents, travArgs);
Expand Down Expand Up @@ -269,8 +272,11 @@ algorithm
(traverserArgs) :=
matchcontinue(inEqns, traverserArgs.recursiveStrongComponentRun, isStateVarInvoled)
local
BackendDAE.EqSystem syst;
BackendDAE.IncidenceMatrix adjMatrix;

list<BackendDAE.Equation> newResEqns;
list<BackendDAE.Var> newResVars, newAuxVars;
list<BackendDAE.Var> newResVars, newAuxVars, discVars, contVars;
BackendDAE.Var var;
BackendDAE.Variables systemVars;
BackendDAE.Var dummyVar;
Expand All @@ -282,6 +288,8 @@ algorithm
DAE.FunctionTree funcsTree;
list<DAE.ComponentRef> crlst;
Boolean b1, b2;
list<BackendDAE.Equation> discEqns;
list<BackendDAE.Equation> contEqns;

DAE.Algorithm alg;
DAE.ElementSource source;
Expand Down Expand Up @@ -482,14 +490,26 @@ algorithm

case(_, false, _)
algorithm
vars := inVars;
for e in inEqns loop

(discVars, contVars) := List.splitOnTrue(inVars, BackendVariable.isVarDiscrete);
(discEqns, contEqns) := getDiscAndContEqns(inVars, inEqns, discVars, contVars, traverserArgs.shared.functionTree);

// create discrete
for e in discEqns loop
size := BackendEquation.equationSize(e);
newAuxVars := List.firstN(discVars, size);
traverserArgs := traverserStrongComponents({e}, newAuxVars, {}, {}, traverserArgs);
discVars := List.stripN(discVars, size);
end for;

// create continuous
for e in contEqns loop
size := BackendEquation.equationSize(e);
newAuxVars := List.firstN(vars, size);
newAuxVars := List.firstN(contVars, size);
traverserArgs.recursiveStrongComponentRun := true;
traverserArgs := traverserStrongComponents({e}, newAuxVars, {}, {}, traverserArgs);
traverserArgs.recursiveStrongComponentRun := false;
vars := List.stripN(vars, size);
contVars := List.stripN(contVars, size);
end for;
then
(traverserArgs);
Expand All @@ -504,6 +524,54 @@ algorithm
end matchcontinue;
end traverserStrongComponents;

function getDiscAndContEqns
input list<BackendDAE.Var> inAllVars;
input list<BackendDAE.Equation> inAllEqns;
input list<BackendDAE.Var> inDiscVars;
input list<BackendDAE.Var> inContVars;
input DAE.FunctionTree functionTree;
output list<BackendDAE.Equation> discEqns;
output list<BackendDAE.Equation> contEqns;
protected
BackendDAE.EqSystem syst;
BackendDAE.IncidenceMatrix adjMatrix;

list<Integer> varsIndex, eqnIndex;
array<Integer> assignVarEqn "eqn := assignVarEqn[var]";
array<Integer> assignEqnVar "var := assignEqnVar[eqn]";

array<Integer> mapEqnScalarArray;
constant Boolean debug = false;
algorithm
try
// create syst for a matching
syst := BackendDAEUtil.createEqSystem(BackendVariable.listVar1(inAllVars), BackendEquation.listEquation(inAllEqns) );
if debug then BackendDump.printEqSystem(syst); end if;
(adjMatrix, _, _, mapEqnScalarArray) := BackendDAEUtil.incidenceMatrixScalar(syst, BackendDAE.NORMAL(), SOME(functionTree));
if debug then BackendDump.dumpIncidenceMatrix(adjMatrix); end if;
(assignVarEqn, assignEqnVar, true) := Matching.RegularMatching(adjMatrix, BackendDAEUtil.systemSize(syst), BackendDAEUtil.systemSize(syst));
if debug then BackendDump.dumpMatching(assignVarEqn); end if;

// get discrete vars indexes and then the equations
varsIndex := BackendVariable.getVarIndexFromVars(inDiscVars, syst.orderedVars);
if debug then print("discVarsIndex: "); BackendDump.dumpIncidenceRow(varsIndex); end if;
eqnIndex := List.map1(varsIndex, Array.getIndexFirst, assignVarEqn);
if debug then print("discEqnIndex: "); BackendDump.dumpIncidenceRow(eqnIndex); end if;
eqnIndex := List.unique(list(mapEqnScalarArray[i] for i in eqnIndex));
discEqns := BackendEquation.getList(eqnIndex, syst.orderedEqs);
if debug then BackendDump.equationListString(discEqns, "Discrete Equations"); end if;

// get continuous equations
varsIndex := BackendVariable.getVarIndexFromVars(inContVars, syst.orderedVars);
eqnIndex := List.map1(varsIndex, Array.getIndexFirst, assignVarEqn);
eqnIndex := List.unique(list(mapEqnScalarArray[i] for i in eqnIndex));
if debug then print("contEqnIndex: "); BackendDump.dumpIncidenceRow(eqnIndex); end if;
contEqns := BackendEquation.getList(eqnIndex, syst.orderedEqs);
if debug then BackendDump.equationListString(contEqns, "Continuous Equations"); end if;
else
fail();
end try;
end getDiscAndContEqns;

function addVarsGlobalData
input output BackendDAE.BackendDAEModeData globalDAEData;
Expand All @@ -513,8 +581,6 @@ protected
algorithm
// prepare algebraic states
vars := List.filterOnTrue(inVars, BackendVariable.isNonStateVar);
// check for discrete vars
{} := List.filterOnTrue(vars, BackendVariable.isVarDiscrete);
vars := list(BackendVariable.setVarKind(v, BackendDAE.ALG_STATE()) for v in vars);
//print("Alg vars: " + BackendDump.varListString(vars, "") + "\n");
globalDAEData.algStateVars := listAppend(vars, globalDAEData.algStateVars);
Expand Down
4 changes: 2 additions & 2 deletions SimulationRuntime/c/simulation/simulation_runtime.cpp
View file
Expand Up @@ -455,7 +455,7 @@ int startNonInteractiveSimulation(int argc, char**argv, DATA* data, threadData_t
}
/* if daeMode is turned on than use also the ida solver */
if (omc_flag[FLAG_DAE_MODE] && std::string("ida") != data->simulationInfo->solverMethod) {
data->simulationInfo->solverMethod = std::string("ida").c_str();
data->simulationInfo->solverMethod = GC_strdup(std::string("ida").c_str());
infoStreamPrint(LOG_SIMULATION, 0, "overwrite solver method: %s [DAEmode works only with IDA solver]", data->simulationInfo->solverMethod);
}

Expand Down Expand Up @@ -644,7 +644,7 @@ static int callSolver(DATA* simData, threadData_t *threadData, string init_initM
/* if no states are present, then we can
* use euler method, since it does nothing.
*/
if (simData->modelData->nStates < 1 && solverID != S_OPTIMIZATION && solverID != S_SYM_SOLVER) {
if (simData->modelData->nStates < 1 && solverID != S_OPTIMIZATION && solverID != S_SYM_SOLVER && !compiledInDAEMode) {
solverID = S_EULER;
}

Expand Down
51 changes: 26 additions & 25 deletions SimulationRuntime/c/simulation/solver/ida_solver.c
View file
Expand Up @@ -555,22 +555,21 @@ ida_solver_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo
/* configure algebraic variables as such */
if (idaData->daeMode)
{
if (!omc_flag[FLAG_NO_SUPPRESS_ALG])
if (omc_flag[FLAG_NO_SUPPRESS_ALG])
{
flag = IDASetSuppressAlg(idaData->ida_mem, TRUE);
if (checkIDAflag(flag)){
throwStreamPrint(threadData, "##IDA## Suppress algebraic variables in the local error test failed");
}
}
for(i=0; i<idaData->N; ++i)
{
tmp[i] = (i<data->modelData->nStates)? 1.0: 0.0;
}

for(i=0; i<idaData->N; ++i)
{
tmp[i] = (i<data->modelData->nStates)? 1.0: 0.0;
}

flag = IDASetId(idaData->ida_mem, N_VMake_Serial(idaData->N,tmp));
if (checkIDAflag(flag)){
throwStreamPrint(threadData, "##IDA## Mark algebraic variables as such failed!");
}
flag = IDASetId(idaData->ida_mem, N_VMake_Serial(idaData->N,tmp));
if (checkIDAflag(flag)){
throwStreamPrint(threadData, "##IDA## Mark algebraic variables as such failed!");
}
}

Expand Down Expand Up @@ -703,6 +702,7 @@ ida_event_update(DATA* data, threadData_t *threadData)
if (compiledInDAEMode == 3){
if (initializedSolver){

data->simulationInfo->needToIterate = 0;
/* get new values from data -> TODO: update all discrete */
data->simulationInfo->discreteCall = 1;

Expand Down Expand Up @@ -731,9 +731,9 @@ ida_event_update(DATA* data, threadData_t *threadData)
}

/* increase limits of the non-linear solver */
IDASetMaxNumStepsIC(idaData->ida_mem, 2*data->modelData->nStates*10);
IDASetMaxNumJacsIC(idaData->ida_mem, 2*data->modelData->nStates*10);
IDASetMaxNumItersIC(idaData->ida_mem, 2*data->modelData->nStates*10);
IDASetMaxNumStepsIC(idaData->ida_mem, 2*idaData->N*10);
IDASetMaxNumJacsIC(idaData->ida_mem, 2*idaData->N*10);
IDASetMaxNumItersIC(idaData->ida_mem, 2*idaData->N*10);
/* Calc Consistent y_algebraic and y_prime with current y */
flag = IDACalcIC(idaData->ida_mem, IDA_YA_YDP_INIT, data->localData[0]->timeValue+init_h);

Expand Down Expand Up @@ -1029,11 +1029,12 @@ ida_solver_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
/* initialize states and der(states) */
if (idaData->daeMode)
{
idaData->residualFunction(data->localData[0]->timeValue, idaData->y, idaData->yp, idaData->newdelta, idaData);
memcpy(data->localData[0]->realVars, idaData->states, sizeof(double)*data->modelData->nStates);
// and also algebraic vars
data->simulationInfo->daeModeData->setAlgebraicDAEVars(data, threadData, idaData->states + data->modelData->nStates);
memcpy(data->localData[0]->realVars + data->modelData->nStates, idaData->statesDer, sizeof(double)*data->modelData->nStates);
sData->timeValue = solverInfo->currentTime;
idaData->residualFunction(sData->timeValue, idaData->y , idaData->yp, idaData->newdelta, idaData);
}

/* sensitivity mode */
Expand Down Expand Up @@ -1149,7 +1150,6 @@ int residualFunctionIDA(double time, N_Vector yy, N_Vector yp, N_Vector res, voi
{
setContext(data, &time, CONTEXT_ODE);
}
timeBackup = data->localData[0]->timeValue;
data->localData[0]->timeValue = time;

saveJumpState = threadData->currentErrorStage;
Expand Down Expand Up @@ -1233,8 +1233,6 @@ int residualFunctionIDA(double time, N_Vector yy, N_Vector yp, N_Vector res, voi

threadData->currentErrorStage = saveJumpState;

data->localData[0]->timeValue = timeBackup;

if (data->simulationInfo->currentContext == CONTEXT_ODE){
unsetContext(data);
}
Expand All @@ -1250,8 +1248,9 @@ int rootsFunctionIDA(double time, N_Vector yy, N_Vector yp, double *gout, void*
IDA_SOLVER* idaData = (IDA_SOLVER*)userData;
DATA* data = (DATA*)(((IDA_USERDATA*)idaData->simData)->data);
threadData_t* threadData = (threadData_t*)(((IDA_USERDATA*)((IDA_SOLVER*)userData)->simData)->threadData);
double *states = N_VGetArrayPointer(yy);
double *statesDer = N_VGetArrayPointer(yp);

double timeBackup;
int saveJumpState;

infoStreamPrint(LOG_SOLVER_V, 1, "### eval rootsFunctionIDA ###");
Expand All @@ -1270,21 +1269,24 @@ int rootsFunctionIDA(double time, N_Vector yy, N_Vector yp, double *gout, void*
saveJumpState = threadData->currentErrorStage;
threadData->currentErrorStage = ERROR_EVENTSEARCH;

timeBackup = data->localData[0]->timeValue;
data->localData[0]->timeValue = time;

if (idaData->daeMode)
{
memcpy(data->localData[0]->realVars, idaData->states, sizeof(double)*data->modelData->nStates);
data->simulationInfo->daeModeData->setAlgebraicDAEVars(data, threadData, idaData->states + data->modelData->nStates);
memcpy(data->localData[0]->realVars + data->modelData->nStates, idaData->statesDer, sizeof(double)*data->modelData->nStates);
memcpy(data->localData[0]->realVars, states, sizeof(double)*data->modelData->nStates);
data->simulationInfo->daeModeData->setAlgebraicDAEVars(data, threadData, states + data->modelData->nStates);
memcpy(data->localData[0]->realVars + data->modelData->nStates, statesDer, sizeof(double)*data->modelData->nStates);
}

data->localData[0]->timeValue = time;

/* read input vars */
externalInputUpdate(data);
data->callback->input_function(data, threadData);
/* eval needed equations*/
if (idaData->daeMode && compiledInDAEMode == 1){}
else if (idaData->daeMode && compiledInDAEMode == 3)
{
idaData->residualFunction(time, yy, yp, idaData->newdelta, idaData);
}
else
{
data->callback->function_ZeroCrossingsEquations(data, threadData);
Expand All @@ -1293,7 +1295,6 @@ int rootsFunctionIDA(double time, N_Vector yy, N_Vector yp, double *gout, void*
data->callback->function_ZeroCrossings(data, threadData, gout);

threadData->currentErrorStage = saveJumpState;
data->localData[0]->timeValue = timeBackup;

/* scale data again */
if (omc_flag[FLAG_IDA_SCALING])
Expand Down
18 changes: 11 additions & 7 deletions SimulationRuntime/c/simulation/solver/solver_main.c
View file
Expand Up @@ -687,15 +687,19 @@ int solver_main(DATA* data, threadData_t *threadData, const char* init_initMetho
externalInputallocate(data);
/* set tolerance for ZeroCrossings */
setZCtol(fmin(data->simulationInfo->stepSize, data->simulationInfo->tolerance));

omc_alloc_interface.collect_a_little();
/* initialize all parts of the model */
retVal = initializeModel(data, threadData, init_initMethod, init_file, init_time);
omc_alloc_interface.collect_a_little();

if(0 == retVal) {
retVal = initializeSolverData(data, threadData, &solverInfo);
initSolverInfo = 1;
/* initialize solver data */
/* For the DAEmode we need to initialize solverData before the initialization,
* since the solver is used to obtain consistent values also via updateDiscreteSystem
*/
retVal = initializeSolverData(data, threadData, &solverInfo);
initSolverInfo = 1;

/* initialize all parts of the model */
if (0 == retVal){
retVal = initializeModel(data, threadData, init_initMethod, init_file, init_time);
omc_alloc_interface.collect_a_little();
}

#if !defined(OMC_MINIMAL_RUNTIME)
Expand Down

0 comments on commit c2feb3f

Please sign in to comment.