Skip to content

Commit

Permalink
fix in cpp template for algoop evaluation
Browse files Browse the repository at this point in the history
changed zero crossing tollernce in cpp template
change absolute tollerance in cvode 


git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@21306 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
niklwors committed Jun 27, 2014
1 parent 0d74bad commit 60b1018
Show file tree
Hide file tree
Showing 4 changed files with 47 additions and 24 deletions.
33 changes: 17 additions & 16 deletions Compiler/Template/CodegenCpp.tpl
Expand Up @@ -7383,12 +7383,14 @@ template equationLinearOrNonLinear(SimEqSystem eq, Context context,Text &varDecl
double* algloop<%index%>Vars = new double[dim<%index%>];
_algLoop<%index%>->getReal(algloop<%index%>Vars );
bool restatDiscrete<%index%>= false;
try
IContinuous::UPDATETYPE calltype = _callType;
try
{
_algLoop<%index%>->evaluate();
if( _callType == IContinuous::DISCRETE )
{
_algLoop<%index%>->evaluate();
while(restart<%index%> && !(iterations<%index%>++>500))
{
Expand Down Expand Up @@ -7420,7 +7422,6 @@ template equationLinearOrNonLinear(SimEqSystem eq, Context context,Text &varDecl
{
try
{ //workaround: try to solve algoop discrete (evaluate all zero crossing conditions) since we do not have the information which zercrossing contains a algloop var
IContinuous::UPDATETYPE calltype = _callType;
_callType = IContinuous::DISCRETE;
_algLoop<%index%>->setReal(algloop<%index%>Vars );
_algLoopSolver<%index%>->solve();
Expand Down Expand Up @@ -9886,30 +9887,30 @@ template giveZeroFunc3(Integer index1, Exp relation, Text &varDecls /*BUFP*/,Tex
case LESS(__) then
<<
if(_conditions[<%zerocrossingIndex%>])
f[<%index1%>]=(<%e1%> -EPSILON - <%e2%>);
f[<%index1%>]=(<%e1%> - 10000*EPSILON - <%e2%>);
else
f[<%index1%>]=(<%e2%> - <%e1%> - 10*EPSILON);
f[<%index1%>]=(<%e2%> - <%e1%> - 10000*EPSILON);
>>
case LESSEQ(__) then
<<
if(_conditions[<%zerocrossingIndex%>])
f[<%index1%>]=(<%e1%> - EPSILON - <%e2%>);
f[<%index1%>]=(<%e1%> - 10000*EPSILON - <%e2%>);
else
f[<%index1%>]=(<%e2%> - <%e1%> - EPSILON);
f[<%index1%>]=(<%e2%> - <%e1%> - 10000*EPSILON);
>>
case GREATER(__) then
<<
if(_conditions[<%zerocrossingIndex%>])
f[<%index1%>]=(<%e2%> - <%e1%> - EPSILON);
f[<%index1%>]=(<%e2%> - <%e1%> - 10000*EPSILON);
else
f[<%index1%>]=(<%e1%> -10*EPSILON - <%e2%>);
f[<%index1%>]=(<%e1%> - 10000*EPSILON - <%e2%>);
>>
case GREATEREQ(__) then
<<
if(_conditions[<%zerocrossingIndex%>])
f[<%index1%>]=(<%e2%> - <%e1%> - EPSILON);
f[<%index1%>]=(<%e2%> - <%e1%> - 10000*EPSILON);
else
f[<%index1%>]=(<%e1%> - EPSILON - <%e2%>);
f[<%index1%>]=(<%e1%> - 10000*EPSILON - <%e2%>);
>>
else
<<
Expand Down Expand Up @@ -9939,17 +9940,17 @@ template giveZeroFunc3(Integer index1, Exp relation, Text &varDecls /*BUFP*/,Sim
case LESSEQ(__) then
<<
if(_event_handling.pre(_condition<%zerocrossingIndex%>,"_condition<%zerocrossingIndex%>"))
f[<%index1%>]=(<%e1%>-EPSILON-<%e2%>);
f[<%index1%>]=(<%e1%>-10000*EPSILON-<%e2%>);
else
f[<%index1%>]=(<%e2%>-<%e1%>-EPSILON);
f[<%index1%>]=(<%e2%>-<%e1%>-10000*EPSILON);
>>
case GREATER(__)
case GREATEREQ(__) then
<<
if(_event_handling.pre(_condition<%zerocrossingIndex%>,"_condition<%zerocrossingIndex%>"))
f[<%index1%>]=(<%e2%>-<%e1%>-EPSILON);
f[<%index1%>]=(<%e2%>-<%e1%>-10000*EPSILON);
else
f[<%index1%>]=(<%e1%>-EPSILON-<%e2%>);
f[<%index1%>]=(<%e1%>-10000*EPSILON-<%e2%>);
>>
end match
end giveZeroFunc3;
Expand Down
21 changes: 18 additions & 3 deletions SimulationRuntime/cpp/Solver/CVode/CVode.cpp
Expand Up @@ -15,6 +15,7 @@ Cvode::Cvode(IMixedSystem* system, ISolverSettings* settings)
, _hOut (0.0)
, _tOut (0.0)
, _zeroSign (NULL)
, _absTol (NULL)
,_cvode_initialized(false)
,_tLastEvent(0.0)
,_event_n(0)
Expand All @@ -31,11 +32,14 @@ Cvode::~Cvode()
delete [] _zInit;
if(_zeroSign)
delete [] _zeroSign;
if(_absTol)
delete [] _absTol;
if(_cvode_initialized)
{
N_VDestroy_Serial(_CV_y0);
N_VDestroy_Serial(_CV_y);
N_VDestroy_Serial(_CV_yWrite);
N_VDestroy_Serial(_CV_absTol);
CVodeFree(&_cvodeMem);
}
}
Expand Down Expand Up @@ -69,12 +73,14 @@ void Cvode::initialize()
if(_zInit) delete [] _zInit;
if(_zWrite) delete [] _zWrite;
if(_zeroSign) delete [] _zeroSign;
if(_absTol) delete [] _absTol;


_z = new double[_dimSys];
_zInit = new double[_dimSys];
_zWrite = new double[_dimSys];
_zeroSign = new int[_dimZeroFunc];
_absTol = new double[_dimSys];

memset(_z,0,_dimSys*sizeof(double));
memset(_zInit,0,_dimSys*sizeof(double));
Expand Down Expand Up @@ -107,9 +113,17 @@ void Cvode::initialize()
_continuous_system->getContinuousStates(_zInit);
memcpy(_z,_zInit,_dimSys*sizeof(double));

// Get nominal values
_continuous_system->getNominalStates(_absTol);
for(int i=0;i<_dimSys;i++)
_absTol[i]*= dynamic_cast<ISolverSettings*>(_cvodesettings)->getATol();


_CV_y0 = N_VMake_Serial(_dimSys, _zInit);
_CV_y = N_VMake_Serial(_dimSys, _z);
_CV_yWrite = N_VMake_Serial(_dimSys, _zWrite);
_CV_absTol = N_VMake_Serial(_dimSys, _absTol);

if(check_flag((void*)_CV_y0, "N_VMake_Serial", 0))
{
_idid = -5;
Expand All @@ -125,7 +139,7 @@ void Cvode::initialize()
}

// Set Tolerances
_idid = CVodeSStolerances(_cvodeMem,dynamic_cast<ISolverSettings*>(_cvodesettings)->getRTol(),dynamic_cast<ISolverSettings*>(_cvodesettings)->getATol());// RTOL and ATOL
_idid = CVodeSVtolerances(_cvodeMem,dynamic_cast<ISolverSettings*>(_cvodesettings)->getRTol(),_CV_absTol);// RTOL and ATOL
if(_idid < 0)
throw std::invalid_argument("CVode::initialize()");

Expand Down Expand Up @@ -278,7 +292,7 @@ void Cvode::CVodeCore()
{
_idid = CVodeReInit(_cvodeMem, _tCurrent, _CV_y);
_idid = CVodeSetStopTime(_cvodeMem, _tEnd);
_idid = CVodeSetInitStep(_cvodeMem, 1e-6);
_idid = CVodeSetInitStep(_cvodeMem, 1e-12);
if(_idid <0)
throw std::runtime_error("CVode::ReInit");

Expand All @@ -287,7 +301,7 @@ void Cvode::CVodeCore()
_cv_rt = CVode(_cvodeMem, _tEnd, _CV_y, &_tCurrent, CV_ONE_STEP);

_idid = CVodeGetNumSteps(_cvodeMem, &_locStps);
_accStps +=_locStps;

_idid = CVodeGetLastStep(_cvodeMem,&_h);
//Ausgabe
writeCVodeOutput(_tCurrent,_h,_locStps);
Expand Down Expand Up @@ -384,6 +398,7 @@ void Cvode::CVodeCore()
_continuous_system->setContinuousStates(NV_DATA_S(_CV_y));
_continuous_system->evaluateODE(IContinuous::CONTINUOUS);
writeToFile(0, _tEnd, _h);
_accStps +=_locStps;
_solverStatus = DONE;
}
}
Expand Down
14 changes: 10 additions & 4 deletions SimulationRuntime/cpp/Solver/CVode/CVode.h
Expand Up @@ -7,6 +7,7 @@
#include <cvodes/cvodes.h>
#include <nvector/nvector_serial.h>
#include <cvodes/cvodes_dense.h>
#include <cvodes/cvodes_spgmr.h>

/*
#include <boost/log/core.hpp>
Expand Down Expand Up @@ -162,9 +163,10 @@ class Cvode


double
*_z, ///< Output - (Current) State vector
*_zInit, ///< Temp - Initial state vector
*_zWrite; ///< Temp - Zustand den das System rausschreibt
*_z, ///< Output - (Current) State vector
*_zInit, ///< Temp - Initial state vector
*_zWrite, ///< Temp - Zustand den das System rausschreibt
*_absTol; /// - Vektor für absolute Toleranzen

double
_hOut; ///< Temp - Ouput step size for dense output
Expand All @@ -189,7 +191,11 @@ double
N_Vector
_CV_y0, ///< Temp - Initial values in the Cvode Format
_CV_y, ///< Temp - State in Cvode Format
_CV_yWrite; ///< Temp - Vector for dense out
_CV_yWrite, ///< Temp - Vector for dense out
_CV_absTol;



bool _cvode_initialized;


Expand Down
3 changes: 2 additions & 1 deletion SimulationRuntime/cpp/Solver/Kinsol/Kinsol.cpp
Expand Up @@ -166,7 +166,8 @@ void Kinsol::solve()
long int dimSys = _dimSys;
long int irtrn = 0; // Retrun-flag of Fortran code _algLoop->getReal(_y);

_algLoop->getRHS(_f);
_algLoop->evaluate();
_algLoop->getRHS(_f);
_algLoop->getSystemMatrix(_jac);
dgesv_(&dimSys,&dimRHS,_jac,&dimSys,_helpArray,_f,&dimSys,&irtrn);
memcpy(_y,_f,_dimSys*sizeof(double));
Expand Down

0 comments on commit 60b1018

Please sign in to comment.