Skip to content

Commit

Permalink
- multirate rk12 solver for cpp runtime
Browse files Browse the repository at this point in the history
  • Loading branch information
vwaurich committed Jan 13, 2016
1 parent f701c46 commit fd65350
Show file tree
Hide file tree
Showing 16 changed files with 1,884 additions and 4 deletions.
85 changes: 82 additions & 3 deletions Compiler/Template/CodegenCpp.tpl
Expand Up @@ -3661,11 +3661,25 @@ case SIMCODE(modelInfo = MODELINFO(__)) then
>>
end simulationCppFile;

template partitionInfoInit(Integer numPartitions, Integer numStates, list<Integer> stateActivators)
::=
let stateActs = (stateActivators |> act hasindex i0 => '_stateActivator[<%i0%>] = <%intSub(act,1)%>;' ;separator="\n")
<<
//partitioning of the system, all partitions are active at t0
_dimPartitions = <%numPartitions%>;
_partitionActivation = new bool[_dimPartitions]();
memset(_partitionActivation,true,_dimPartitions*sizeof(bool));
_stateActivator = new int[<%numStates%>]();
<%stateActs%>
>>
end partitionInfoInit;

template generateSimulationCppConstructorContent(SimCode simCode, Context context, Text& extraFuncs, Text& extraFuncsDecl, Text extraFuncsNamespace, Text stateDerVectorName /*=__zDot*/, Boolean useFlatArrayNotation)
::=
match simCode
case SIMCODE(modelInfo = MODELINFO(__)) then
case SIMCODE(modelInfo = MODELINFO(varInfo = vi as VARINFO(__)), partitionData=PARTITIONDATA(__)) then
let className = lastIdentOfPath(modelInfo.name)
let partitionInit = if Flags.isSet(Flags.MULTIRATE_PARTITION) then partitionInfoInit(partitionData.numPartitions, vi.numStateVars, partitionData.stateToActivators) else ""
<<
defineConstVals();
defineAlgVars();
Expand Down Expand Up @@ -3730,6 +3744,8 @@ match simCode
%>
//DAEs are not supported yet, Index reduction is enabled
_dimAE = 0; // algebraic equations
<%partitionInit%>

//Initialize the state vector
SystemDefaultImplementation::initialize();
//Instantiate auxiliary object for event handling functionality
Expand Down Expand Up @@ -7189,6 +7205,10 @@ match modelInfo
bool* _pointerToBoolVars;
string* _pointerToStringVars;

int _dimPartitions;
bool* _partitionActivation;
int* _stateActivator;

<%memberVariableDefinitions%>
<%memberPreVariableDefinitions%>
<%conditionvariables%>
Expand Down Expand Up @@ -7502,9 +7522,29 @@ template DefaultImplementationCode(SimCode simCode, Text& extraFuncs, Text& extr
SystemDefaultImplementation::setStringStartValue(var, val);
}

void <%lastIdentOfPath(modelInfo.name)%>::setNumPartitions(int numPartitions)
{
_dimPartitions = numPartitions;
}

int <%lastIdentOfPath(modelInfo.name)%>::getNumPartitions()
{
return _dimPartitions;
}
void <%lastIdentOfPath(modelInfo.name)%>::setPartitionActivation(bool* partitions)
{
_partitionActivation = partitions;
}

void <%lastIdentOfPath(modelInfo.name)%>::getPartitionActivation(bool* partitions)
{
partitions = _partitionActivation;
}

int <%lastIdentOfPath(modelInfo.name)%>::getActivator(int state)
{
return (int)_stateActivator[state];
}

// Provide the right hand side (according to the index)
void <%lastIdentOfPath(modelInfo.name)%>::getRHS(double* f)
Expand Down Expand Up @@ -7718,6 +7758,12 @@ case SIMCODE(modelInfo = MODELINFO(vars = vars as SIMVARS(__))) then
virtual void setIntStartValue(int& var,int val);
virtual void setStringStartValue(string& var,string val);

virtual void setNumPartitions(int numPartitions);
virtual int getNumPartitions();
virtual void setPartitionActivation(bool* partitions);
virtual void getPartitionActivation(bool* partitions);
virtual int getActivator(int state);

//Resets all time events

// Provide variables with given index to the system
Expand Down Expand Up @@ -9777,7 +9823,7 @@ end clockIntervalsInit;
template dimension1(SimCode simCode ,Text& extraFuncs,Text& extraFuncsDecl,Text extraFuncsNamespace)
::=
match simCode
case SIMCODE(modelInfo = MODELINFO(varInfo = vi as VARINFO(__)))
case SIMCODE(modelInfo = MODELINFO(varInfo = vi as VARINFO(__)), partitionData = PARTITIONDATA(__))
then
let numRealVars = numRealvars(modelInfo)
let numIntVars = numIntvars(modelInfo)
Expand All @@ -9790,6 +9836,7 @@ template dimension1(SimCode simCode ,Text& extraFuncs,Text& extraFuncsDecl,Text
_dimInteger = <%numIntVars%>;
_dimString = <%numStringVars%>;
_dimReal = <%numRealVars%>;
_dimPartitions = <%partitionData.numPartitions%>;
>>
end dimension1;

Expand Down Expand Up @@ -12847,12 +12894,18 @@ end createEvaluateConditions;

template createEvaluate(list<list<SimEqSystem>> odeEquations, SimCode simCode ,Text& extraFuncs,Text& extraFuncsDecl,Text extraFuncsNamespace, Context context, Boolean createMeasureTime)
::=
match simCode
case SIMCODE(partitionData = PARTITIONDATA(partitions = partitions, activatorsForPartitions=activatorsForPartitions)) then
//case MODELINFO(vars = vars as SIMVARS(__))
let className = lastIdentOfPathFromSimCode(simCode , &extraFuncs , &extraFuncsDecl, extraFuncsNamespace)
let &varDecls = buffer "" /*BUFD*/

let equation_ode_func_calls = (List.partition(List.flatten(odeEquations), 100) |> eqs hasindex i0 =>
let equation_ode_func_calls = if not Flags.isSet(Flags.MULTIRATE_PARTITION) then (List.partition(List.flatten(odeEquations), 100) |> eqs hasindex i0 =>
createEvaluateWithSplit(i0, context, eqs, "evaluateODE", className, simCode, &extraFuncs, &extraFuncsDecl, extraFuncsNamespace)
;separator="\n")
else ( List.intRange(partitionData.numPartitions) |> partIdx =>
createEvaluatePartitions(partIdx, context, List.flatten(odeEquations), listGet(partitions, partIdx),
listGet(activatorsForPartitions,partIdx), className,simCode, &extraFuncs, &extraFuncsDecl, extraFuncsNamespace) ;separator="\n")
<<
void <%className%>::evaluateODE(const UPDATETYPE command)
{
Expand All @@ -12865,6 +12918,31 @@ template createEvaluate(list<list<SimEqSystem>> odeEquations, SimCode simCode ,T
>>
end createEvaluate;

template createEvaluatePartitions(Integer partIdx, Context context, list<SimEqSystem> odeEquations, list<Integer> partition, list<Integer> activators, String className, SimCode simCode, Text& extraFuncs,Text& extraFuncsDecl,Text extraFuncsNamespace)
::=
let &varDecls = buffer "" /*BUFD*/
let condition = partitionCondition(activators)
let equation_func_calls = (SimCodeUtil.getSimEqSystemsByIndexLst(partition,odeEquations) |> eq =>
equation_function_call(eq, context, &varDecls /*BUFC*/, simCode, &extraFuncs, &extraFuncsDecl, extraFuncsNamespace, "evaluate")
;separator="\n")
<<
// Partition <%partIdx%>
if (<%condition%>)
{
<%varDecls%>
<%equation_func_calls%>
}
>>
end createEvaluatePartitions;

template partitionCondition(list<Integer> partitions)
::=
let bVec = (partitions |> part => "_partitionActivation[" + intSub(part,1) + "]" ;separator=" || ")
<<
<%bVec%>
>>
end partitionCondition;

template createEvaluateZeroFuncs( list<SimEqSystem> equationsForZeroCrossings, SimCode simCode ,Text& extraFuncs,Text& extraFuncsDecl,Text extraFuncsNamespace, Context context)
::=
let className = lastIdentOfPathFromSimCode(simCode , &extraFuncs , &extraFuncsDecl, extraFuncsNamespace)
Expand Down Expand Up @@ -12908,6 +12986,7 @@ template createEvaluateWithSplit(Integer sectionIndex, Context context, list<Sim
>>
end createEvaluateWithSplit;


/*
//! Evaluates only the equations whose indexs are passed to it.
bool <%className%>::evaluate_selective(const std::vector<int>& indices) {
Expand Down
7 changes: 7 additions & 0 deletions SimulationRuntime/cpp/CMakeLists.txt
Expand Up @@ -74,6 +74,7 @@ SET(PeerName ${LIBPREFIX}Peer)
SET(CppDASSLName ${LIBPREFIX}CppDASSL)
SET(RTRKName ${LIBPREFIX}RTRK)
SET(EulerName ${LIBPREFIX}Euler)
SET(RK12Name ${LIBPREFIX}RK12)
#SET(kluName ${LIBPREFIX}klu)
SET(RTEulerName ${LIBPREFIX}RTEuler)
SET(IdaName ${LIBPREFIX}Ida)
Expand Down Expand Up @@ -568,6 +569,8 @@ add_subdirectory (Core/SimController)

#add Euler solver project
add_subdirectory (Solver/Euler)
add_subdirectory (Solver/RK12)


add_subdirectory (Solver/RTEuler)
add_subdirectory (Solver/Newton)
Expand Down Expand Up @@ -650,6 +653,9 @@ GET_FILENAME_COMPONENT(libUmfPackName ${libUmfPack} NAME)
GET_TARGET_PROPERTY(libEuler ${EulerName} LOCATION)
GET_FILENAME_COMPONENT(libEulerName ${libEuler} NAME)

GET_TARGET_PROPERTY(libRK12 ${RK12Name} LOCATION)
GET_FILENAME_COMPONENT(libRK12Name ${libRK12} NAME)

#GET_TARGET_PROPERTY(libKlu ${kluName} LOCATION)
#GET_FILENAME_COMPONENT(libKluName ${libKlu} NAME)

Expand Down Expand Up @@ -691,6 +697,7 @@ GET_FILENAME_COMPONENT(libFMUName ${libFMU} NAME)

#set (KLU_LIB ${libKluName})
set (EULER_LIB ${libEulerName})
set (RK12_LIB ${libRK12Name})
set (RTEULER_LIB ${libRTEulerName})
set (SETTINGSFACTORY_LIB ${libSetFactoryName})
set (MODELICASYSTEM_LIB ${libModelicaName})
Expand Down
3 changes: 2 additions & 1 deletion SimulationRuntime/cpp/Core/SimController/SimManager.cpp
Expand Up @@ -549,7 +549,8 @@ void SimManager::runSingleProcess()
_solver->setStartTime(_tStart);
_solver->setEndTime(_tEnd);

_solver->solve(_solverTask);
//_solver->solve(_solverTask);
initialize();
_solverTask = ISolver::SOLVERCALL(_solverTask ^ ISolver::RECORDCALL);
/* Logs temporarily disabled
BOOST_LOG_SEV(simmgr_lg::get(), simmgr_normal) <<"Run single process." ; */
Expand Down
9 changes: 9 additions & 0 deletions SimulationRuntime/cpp/Include/Core/System/IContinuous.h
Expand Up @@ -109,6 +109,15 @@ class IContinuous
virtual void setBoolStartValue(bool& var,bool val) = 0;
virtual void setIntStartValue(int& var,int val) = 0;
virtual void setStringStartValue(string& var,string val) = 0;

//in case of solver-based activation of system equations
virtual void setNumPartitions(int numPartitions) = 0;
virtual int getNumPartitions() = 0;
virtual void setPartitionActivation(bool* partitions) = 0;
virtual void getPartitionActivation(bool* partitions) = 0;
virtual int getActivator(int state) = 0;


};
/** @} */ // end of coreSystem
/*
Expand Down
Expand Up @@ -80,6 +80,18 @@ struct SolverOMCFactory : public ObjectFactory<CreationPolicy>
}

}
else if(solvername.compare("rk12")==0)
{
fs::path rk12_path = ObjectFactory<CreationPolicy>::_library_path;
fs::path rk12_name(RK12_LIB);
rk12_path/=rk12_name;
std::cerr << rk12_path.string() << std::endl;
LOADERRESULT result = ObjectFactory<CreationPolicy>::_factory->LoadLibrary(rk12_path.string(),*_solver_type_map);
if (result != LOADER_SUCCESS)
{
throw ModelicaSimulationError(MODEL_FACTORY,"Failed loading RK12 solver library!");
}
}
else if(solvername.compare("peer")==0)
{
fs::path peer_path = ObjectFactory<CreationPolicy>::_library_path;
Expand Down
Expand Up @@ -61,6 +61,18 @@ struct SolverSettingsOMCFactory : public ObjectFactory<CreationPolicy>
}
solver_settings_key.assign("createEulerSettings");
}
else if(solvername.compare("rk12")==0)
{
fs::path rk12_path = ObjectFactory<CreationPolicy>::_library_path;
fs::path rk12_name(RK12_LIB);
rk12_path/=rk12_name;
LOADERRESULT result = ObjectFactory<CreationPolicy>::_factory->LoadLibrary(rk12_name.string(),*_solver_type_map);
if (result != LOADER_SUCCESS)
{
throw ModelicaSimulationError(MODEL_FACTORY,"Failed loading RK12 solver library!");
}
solver_settings_key.assign("createRK12Settings");
}
else if(solvername.compare("peer")==0)
{
fs::path peer_path = ObjectFactory<CreationPolicy>::_library_path;
Expand Down
24 changes: 24 additions & 0 deletions SimulationRuntime/cpp/Include/Solver/RK12/FactoryExport.h
@@ -0,0 +1,24 @@
#pragma once
/** @addtogroup solverEuler
*
* @{
*/
#if defined(__vxworks)
#define BOOST_EXTENSION_SOLVER_DECL
#define BOOST_EXTENSION_SOLVERSETTINGS_DECL
#elif defined(RUNTIME_STATIC_LINKING) && (defined(OMC_BUILD) || defined(SIMSTER_BUILD))
#define BOOST_EXTENSION_SOLVER_DECL
#define BOOST_EXTENSION_STATESELECT_DECL
#define BOOST_EXTENSION_SOLVERSETTINGS_DECL
#define BOOST_EXTENSION_MONITOR_DECL
#elif defined(OMC_BUILD) || defined(SIMSTER_BUILD)
#define BOOST_EXTENSION_SOLVER_DECL BOOST_EXTENSION_IMPORT_DECL
#define BOOST_EXTENSION_STATESELECT_DECL BOOST_EXTENSION_IMPORT_DECL
#define BOOST_EXTENSION_SOLVERSETTINGS_DECL BOOST_EXTENSION_IMPORT_DECL
#define BOOST_EXTENSION_MONITOR_DECL BOOST_EXTENSION_IMPORT_DECL
#else
error "operating system not supported"
#endif
/** @} */ // end of solverEuler


76 changes: 76 additions & 0 deletions SimulationRuntime/cpp/Include/Solver/RK12/IRK12Settings.h
@@ -0,0 +1,76 @@
#pragma once
/** @addtogroup solverRK12
*
* @{
*/


/*****************************************************************************/
/**
Encapsulation of settings for RK12 solver
\date October, 1st, 2008
\author
*/
/*****************************************************************************
Copyright (c) 2008, OSMC
*****************************************************************************/
class IRK12Settings
{

public:
/// Enum to choose the integration method
enum RK12METHOD
{
STEPSIZECONTROL = 0, ///< Explicit RK12
MULTIRATE = 1, ///< Implicit RK12

};

/// Enum to choose the method for zero search
enum ZEROSEARCHMETHOD
{
NO_ZERO_SEARCH = 0, ///< Ignore zero functions
BISECTION = 1, ///< Bisection method
LINEAR_INTERPOLATION = 2, ///< Linear interpolation
};
virtual ~IRK12Settings() {};

/**
Choise of solution method according to RK12METHOD ([0,1,2,3,4,5]; default: 0)
**/
virtual unsigned int getRK12Method() =0;
virtual void setRK12Method(unsigned int) =0;
/**
Choise of method for zero search according to ZEROSEARCHMETHOD ([0,1]; default: 0)
*/
virtual unsigned int getZeroSearchMethod() =0;
virtual void setZeroSearchMethod(unsigned int ) =0;

/**
Determination of number of zeros in one intervall (used only for methods [2,3]) ([true,false]; default: false)
*/
virtual bool getUseSturmSequence() =0;
virtual void setUseSturmSequence(bool) =0;
/**
For implicit methods only. Choise between fixpoint and newton-iteration kann eine Newtoniteration gewählt werden. ([false,true]; default: false = Fixpunktiteration)
*/
virtual bool getUseNewtonIteration() =0;
virtual void setUseNewtonIteration(bool) =0;
/**
Equidistant output(by interpolation polynominal) ([true,false]; default: false)
*/
virtual bool getDenseOutput() =0;
virtual void setDenseOutput(bool) =0;
/**
Tolerance for newton iteration (used when _useNewtonIteration=true) (default: 1e-8)
*/
virtual double getIterTol()=0;
virtual void setIterTol(double)=0;
virtual void load(std::string xml_file)=0;

};
/** @} */ // end of solverRK12

0 comments on commit fd65350

Please sign in to comment.