Skip to content

Commit

Permalink
fix for mixed system iteration in cpp runtime
Browse files Browse the repository at this point in the history 
fix in CVode for calculation of event state


git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@11211 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
@niklwors
niklwors committed Feb 27, 2012
1 parent 83b4203 commit 4c24901
Show file tree
Hide file tree
Showing 9 changed files with 183 additions and 35 deletions.
28 changes: 21 additions & 7 deletions Compiler/susan_codegen/SimCode/CodegenCpp.tpl
View file
Expand Up @@ -223,6 +223,7 @@ case SIMCODE(modelInfo = MODELINFO(__)) then
<%giveZeroFunc1(zeroCrossings,simCode)%>
<%giveConditions(simCode)%>
<%setConditions(simCode)%>
<%saveConditions(simCode)%>
<%isODE(simCode)%>
<%DimZeroFunc(simCode)%>
<%DimInitEquations(simCode)%>
Expand Down Expand Up @@ -1222,7 +1223,7 @@ case SIMCODE(modelInfo = MODELINFO(__)) then
<%initALgloopSolvers%>
//initialize equations
<%initBoundParameters%>
saveConditions();
}
>>
end init;
Expand Down Expand Up @@ -2086,6 +2087,7 @@ case SIMCODE(modelInfo = MODELINFO(__)) then
virtual void setConditions(bool* c);
//Called to check conditions for event-handling
virtual void checkConditions(unsigned int, bool all);
virtual void saveConditions();
//Called to handle all events occured at same time
virtual void handleSystemEvents(const bool* events,update_events_type update_event);
//Called to handle an event
Expand Down Expand Up @@ -3522,7 +3524,7 @@ case SES_MIXED(__) then
bool new_disc_vars<%num%>[<%numDiscVarsStr%>];
bool restart<%num%> = true;
int iter<%num%>=0;
int max_iter<%num%> = <%valuesLenStr%> / <%numDiscVarsStr%>;
int max_iter<%num%> = (<%valuesLenStr%> / <%numDiscVarsStr%>)+1;
while(restart<%num%> && !(iter<%num%> > max_iter<%num%>))
{
<%discVars |> SIMVAR(__) hasindex i0 => 'pre_disc_vars<%num%>[<%i0%>] = <%cref(name)%>;' ;separator="\n"%>
Expand Down Expand Up @@ -5749,19 +5751,19 @@ template handleSystemEvents(list<ZeroCrossing> zeroCrossings,list<SimWhenClause>
<%varDecls%>
bool restart=true;
int iter=0;
while(restart && !(iter++ > 10))
<%zeroCrossingsCode%>
while(restart && !(iter++ > _dimZeroFunc))
{
saveAll();
<%zeroCrossingsCode%>
double h[<%helpvarlength(simCode)%>];
<%helpvarvector(whenClauses,simCode)%>
_event_handling.setHelpVars(h);
//iterate and handle all events inside the eventqueue
restart=_event_handling.IterateEventQueue(events,update_event);
saveAll();
}
saveConditions();
saveConditions();
resetTimeEvents();
if(iter>10){
if(iter>_dimZeroFunc && restart ){
throw std::runtime_error("Number of event iteration steps exceeded. ");}
}
>>
Expand Down Expand Up @@ -5843,6 +5845,18 @@ template setConditions(SimCode simCode)
>>
end setConditions;

template saveConditions(SimCode simCode)
::=
match simCode
case SIMCODE(modelInfo = MODELINFO(__)) then
<<
void <%lastIdentOfPath(modelInfo.name)%>::saveConditions()
{
SystemDefaultImplementation::saveConditions();
}
>>
end saveConditions;

template giveZeroFunc2(list<ZeroCrossing> zeroCrossings, Text &varDecls /*BUFP*/,SimCode simCode)
::=

Expand Down
78 changes: 71 additions & 7 deletions SimulationRuntime/cpp/Source/Solver/CVode/Implementation/CVode.cpp
View file
Expand Up @@ -32,8 +32,11 @@ Cvode::Cvode(IDAESystem* system, ISolverSettings* settings)
, _doubleZero (false)
, _f0 (NULL)
, _f1 (NULL)
, _zeroSign (NULL)
, _zeroFound (false)
,_updateCalls (0)
, _Cond (NULL)
, _zeroInit (NULL)
{
_data = ((void*)this);
}
Expand All @@ -56,6 +59,12 @@ Cvode::~Cvode()
delete [] _f0;
if(_f1)
delete [] _f1;
if(_zeroSign)
delete [] _zeroSign;
if(_zeroInit)
delete [] _zeroInit;
if(_Cond)
delete [] _Cond;

N_VDestroy_Serial(_CV_y0);
N_VDestroy_Serial(_CV_y);
Expand All @@ -69,6 +78,7 @@ void Cvode::init()
{
ISystemProperties* properties = dynamic_cast<ISystemProperties*>(_system);
IContinous* continous_system = dynamic_cast<IContinous*>(_system);
IEvent* event_system = dynamic_cast<IEvent*>(_system);
// Kennzeichnung, dass init()() (vor der Integration) aufgerufen wurde
_idid = 5000;

Expand All @@ -79,6 +89,8 @@ void Cvode::init()
SolverDefaultImplementation::init();

_dimSys = continous_system->getDimVars(IContinous::ALL_STATES);
_dimZeroFunc = event_system->getDimZeroFunc();

int dimAEq = continous_system->getDimVars(IContinous::DIFF_INDEX3);

if(_dimSys <= 0 || dimAEq > 0)
Expand All @@ -94,6 +106,9 @@ void Cvode::init()
if(_zLastSucess) delete [] _zLastSucess;
if(_zLargeStep) delete [] _zLargeStep;
if(_zWrite) delete [] _zWrite;
if(_zeroSign) delete [] _zeroSign;
if(_Cond) delete [] _Cond;
if(_zeroInit) delete [] _zeroInit;

_z = new double[_dimSys];
_zInit = new double[_dimSys];
Expand All @@ -102,6 +117,9 @@ void Cvode::init()
_zWrite = new double[_dimSys];
_f0 = new double[_dimSys];
_f1 = new double[_dimSys];
_zeroSign = new int[_dimZeroFunc];
_Cond = new bool[_dimZeroFunc];
_zeroInit = new bool[_dimZeroFunc];

memset(_z,0,_dimSys*sizeof(double));
memset(_zInit,0,_dimSys*sizeof(double));
Expand Down Expand Up @@ -338,23 +356,67 @@ void Cvode::CVodeCore()
{
IContinous* continous_system = dynamic_cast<IContinous*>(_system);
IEvent* event_system = dynamic_cast<IEvent*>(_system);

// Nullstellen identifizieren, bei denen zu Beginn der Simulation der Zustand nicht entschieden werden kann
// Nullstellen identifizieren, bei denen zu Beginn der Simulation der Zustand nicht entschieden werden kann
event_system->giveZeroFunc(_zeroValLastSuccess,dynamic_cast<ISolverSettings*>(_cvodesettings)->getZeroTol());
event_system->giveConditions(_Cond);
for(int i=0;i<_dimZeroFunc;i++)
{
_Cond[i] = !_Cond[i];
}
event_system->setConditions(_Cond);
event_system->giveZeroFunc(_zeroVal,dynamic_cast<ISolverSettings*>(_cvodesettings)->getZeroTol());
for(int i=0;i<_dimZeroFunc;i++)
{
if(_zeroValLastSuccess[i] == _zeroVal[i])
{
_zeroInit[i] = true;
}else
{
_zeroInit[i] = false;
}
}
for(int i=0;i<_dimZeroFunc;i++)
{
_Cond[i] = !_Cond[i];
}
event_system->setConditions(_Cond);

//


//
while(_solverStatus & IDAESolver::CONTINUE)
{
_zeroFound = false;

//CVode für einen Schritt rufen
_idid = CVode(_cvodeMem, _tEnd, _CV_y, &_tHelp, CV_ONE_STEP);

// Check, ob Schritt erfolgreich
if(check_flag(&_idid, "CVode", 1))
{
_solverStatus = IDAESolver::SOLVERERROR;
break;
}

_accStps++;
_locStps++;

//Beachtung von kritischen (mit Null initialisierten) Nullstellen
if(_tCurrent == 0.0)
{
for(int i=0;i<_dimZeroFunc;i++)
{
if(_zeroInit[i])
{
event_system->checkConditions(i,false);
}
}
event_system->saveConditions();
}


// A root is found
if(_idid == 2)
Expand All @@ -377,7 +439,7 @@ void Cvode::CVodeCore()
_tCurrent = _tEnd;
_zeroFound = false;
}

// Eventiteration starten
if(_zeroFound)
{
Expand All @@ -396,9 +458,11 @@ void Cvode::CVodeCore()
SolverDefaultImplementation::writeToFile(_locStps, _tCurrent, _h);
}

//_idid = CVodeGetRootInfo(_cvodeMem, _zeroSign);
event_system->giveZeroFunc(_zeroVal,dynamic_cast<ISolverSettings*>(_cvodesettings)->getZeroTol());

_idid = CVodeGetRootInfo(_cvodeMem, _zeroSign);
for(int i=0;i<_dimZeroFunc;i++)
_events[i] = bool(_zeroSign[i]);
//event_system->giveZeroFunc(_zeroVal,dynamic_cast<ISolverSettings*>(_cvodesettings)->getZeroTol());

//Event Iteration starten
update_events_type update_event = boost::bind(&SolverDefaultImplementation::updateEventState, this);
event_system->handleSystemEvents(_events,boost::ref(update_event));
Expand Down
8 changes: 6 additions & 2 deletions SimulationRuntime/cpp/Source/Solver/CVode/Implementation/CVode.h
View file
Expand Up @@ -155,7 +155,9 @@ class Cvode

int
_outStps, ///< Output - Total number of output-steps
_updateCalls; ///< Output - Total number of update calls
_updateCalls, ///< Output - Total number of update calls
*_zeroSign;


double
*_z, ///< Output - (Current) State vector
Expand Down Expand Up @@ -189,7 +191,9 @@ class Cvode
bool
_doubleZero, ///< Temp - Flag to denote two zeros in intervall
_zeroFound, ///< Temp - Flag to denote a root in he last step
_bWritten; ///< Temp - Is output already written
_bWritten, ///< Temp - Is output already written
*_Cond,
*_zeroInit;

N_Vector
_CV_y0, ///< Temp - Initial values in the Cvode Format
Expand Down
63 changes: 57 additions & 6 deletions SimulationRuntime/cpp/Source/Solver/Euler/Implementation/Euler.cpp
View file
Expand Up @@ -36,6 +36,8 @@ Euler::Euler(IDAESystem* system, ISolverSettings* settings)
, _f0 (NULL)
, _f1 (NULL)
, _zeroSignIter (NULL)
, _Cond (NULL)
,_zeroInit (NULL)
{
}

Expand All @@ -53,6 +55,8 @@ Euler::~Euler()
delete [] _zLastSucess;
if(_zLargeStep)
delete [] _zLargeStep;
if(_zWrite)
delete [] _zWrite;
if(_denseOutPolynominal)
delete [] _denseOutPolynominal;
if(_zeroSignIter)
Expand All @@ -61,6 +65,10 @@ Euler::~Euler()
delete [] _f0;
if(_f1)
delete [] _f1;
if(_Cond)
delete [] _Cond;
if(_zeroInit)
delete [] _zeroInit;

}

Expand Down Expand Up @@ -90,7 +98,9 @@ void Euler::init()
if(_zInit) delete [] _zInit;
if(_zLastSucess) delete [] _zLastSucess;
if(_zLargeStep) delete [] _zLargeStep;
//if(_zWrite) delete [] _zWrite;
if(_Cond) delete [] _Cond;
if(_zWrite) delete [] _zWrite;
if(_zeroInit) delete [] _zeroInit;

_z = new double[_dimSys];
_zInit = new double[_dimSys];
Expand All @@ -100,6 +110,8 @@ void Euler::init()
_f0 = new double[_dimSys];
_f1 = new double[_dimSys];
_zeroSignIter = new int[_dimZeroFunc];
_Cond = new bool[_dimZeroFunc];
_zeroInit = new bool[_dimZeroFunc];

memset(_z,0,_dimSys*sizeof(double));
memset(_f0,0,_dimSys*sizeof(double));
Expand Down Expand Up @@ -160,6 +172,7 @@ void Euler::solve(const SOLVERCALL command)

if (_eulerSettings && _system)
{
IEvent* event_system = dynamic_cast<IEvent*>(_system);

// Prepare solver and system for integration
if (command & IDAESolver::FIRST_CALL)
Expand Down Expand Up @@ -205,6 +218,33 @@ void Euler::solve(const SOLVERCALL command)
// Get initial values from system, write out initial state vector
solverOutput(_accStps,_tCurrent,_z,_h);

// Nullstellen identifizieren, bei denen zu Beginn der Simulation der Zustand nicht entschieden werden kann
event_system->giveZeroFunc(_zeroValLastSuccess,dynamic_cast<ISolverSettings*>(_eulerSettings)->getZeroTol());
event_system->giveConditions(_Cond);
for(int i=0;i<_dimZeroFunc;i++)
{
_Cond[i] = !_Cond[i];
}
event_system->setConditions(_Cond);
event_system->giveZeroFunc(_zeroVal,dynamic_cast<ISolverSettings*>(_eulerSettings)->getZeroTol());
for(int i=0;i<_dimZeroFunc;i++)
{
if(_zeroValLastSuccess[i] == _zeroVal[i])
{
_zeroInit[i] = true;
}else
{
_zeroInit[i] = false;
}
}
for(int i=0;i<_dimZeroFunc;i++)
{
_Cond[i] = !_Cond[i];
}
event_system->setConditions(_Cond);

//

// Choose integration method
if (_eulerSettings->getEulerMethod() == EulerSettings::EULERFORWARD)
doEulerForward();
Expand Down Expand Up @@ -356,7 +396,7 @@ void Euler::doEulerForward()
void Euler::doEulerBackward()
{
IEvent* event_system = dynamic_cast<IEvent*>(_system);

int numberOfIterations = 0;
double tHelp;
double nu,
Expand Down Expand Up @@ -482,17 +522,28 @@ void Euler::doEulerBackward()
for(int i = 0; i < _dimSys; ++i)
_z[i] += Z[i];

calcFunction(_tCurrent,_z,_f1);
calcFunction(tHelp,_z,_f1);
memcpy(_z1,_z,_dimSys*sizeof(double));

//Beachtung von kritischen (mit Null initialisierten) Nullstellen
if(_tCurrent == 0.0)
{
for(int i=0;i<_dimZeroFunc;i++)
{
if(_zeroInit[i])
{
event_system->checkConditions(i,false);
}
}
event_system->saveConditions();
}

if (_idid != 0)
throw invalid_argument("Euler::doEulerBackward() error" );

++_totStps;
++_accStps;

memcpy(_z1,_z,_dimSys*sizeof(double));

solverOutput(_accStps,tHelp,_z,_h);
doMyZeroSearch();

Expand Down Expand Up @@ -1113,7 +1164,7 @@ void Euler::solverOutput(const int& stp, const double& t, double* z, const doubl


if (_zeroVal)
{
{
// read values of zero functions
event_system->giveZeroFunc(_zeroVal,dynamic_cast<ISolverSettings*>(_eulerSettings)->getZeroTol());

Expand Down

0 comments on commit 4c24901

Please sign in to comment.