Skip to content

Commit

Permalink
- add new nonlinear solver: kinsol (runtime-flag: -nls=kinsol)
Browse files Browse the repository at this point in the history 
git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@14383 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
@lochel
lochel committed Dec 14, 2012
1 parent 59cd639 commit 19fc2f3
Show file tree
Hide file tree
Showing 19 changed files with 984 additions and 569 deletions.
5 changes: 2 additions & 3 deletions Compiler/BackEnd/BackendDAEOptimize.mo
View file
Expand Up @@ -2773,7 +2773,7 @@ algorithm
(r,t,_,dlow_1,m_1,mT_1,v1_1,v2_1,comps_2) = tearingSystem1(dlow,dlow1,m,mT,v1,v2,comps_1);
Debug.fcall(Flags.TEARING_DUMP, BackendDump.dumpIncidenceMatrix, m_1);
Debug.fcall(Flags.TEARING_DUMP, BackendDump.dumpIncidenceMatrixT, mT_1);
Debug.fcall(Flags.TEARING_DUMP, BackendDump.dumpBackendDAE, dlow_1);
Debug.fcall(Flags.TEARING_DUMP, BackendDump.printBackendDAE, dlow_1);
Debug.fcall(Flags.TEARING_DUMP, BackendDump.dumpMatching, v1_1);
//Debug.fcall(Flags.TEARING_DUMP, BackendDump.dumpComponents, comps_2);
Debug.fcall(Flags.TEARING_DUMP, print, "==========\n");
Expand Down Expand Up @@ -2981,8 +2981,7 @@ algorithm
eqSystem = BackendDAE.EQSYSTEM(variables,eqArray,SOME(m_new),SOME(mT_new),matching);
linearDAE = BackendDAE.DAE({eqSystem},shared);
linearDAE1 = BackendDAEUtil.transformBackendDAE(linearDAE,SOME((BackendDAE.NO_INDEX_REDUCTION(), BackendDAE.EXACT())),NONE(),NONE());
print("\n---linearDAE1---\n");
BackendDump.dumpBackendDAE(linearDAE1);
BackendDump.dumpBackendDAE(linearDAE1, "linearDAE1");
BackendDAE.DAE({eqSystem},shared) = linearDAE1;
BackendDAE.EQSYSTEM(variables,eqArray,_,_,matching) = eqSystem;
BackendDAE.MATCHING(v1_new,v2_new,_) = matching;
Expand Down
16 changes: 6 additions & 10 deletions Compiler/BackEnd/BackendDAEUtil.mo
View file
Expand Up @@ -7857,8 +7857,7 @@ algorithm
matchingAlgorithm := getMatchingAlgorithm(strmatchingAlgorithm);
daeHandler := getIndexReductionMethod(strdaeHandler);

Debug.fcall(Flags.DUMP_DAE_LOW, print, "dumpdaelow:\n");
Debug.fcall(Flags.DUMP_DAE_LOW, BackendDump.dumpBackendDAE, inDAE);
Debug.fcall2(Flags.DUMP_DAE_LOW, BackendDump.dumpBackendDAE, inDAE, "dumpdaelow");
System.realtimeTick(BackendDAE.RT_CLOCK_EXECSTAT_BACKEND_MODULES);
// pre optimisation phase
// Frenkel TUD: why is this neccesarray? it only consumes time!
Expand All @@ -7880,8 +7879,7 @@ algorithm
//Debug.execStat("calculateValue",BackendDAE.RT_CLOCK_EXECSTAT_BACKEND_MODULES);
//outSODE := expandAlgorithmsbyInitStmts(sode2);
Debug.execStat("expandAlgorithmsbyInitStmts",BackendDAE.RT_CLOCK_EXECSTAT_BACKEND_MODULES);
Debug.fcall(Flags.DUMP_INDX_DAE, print, "dumpindxdae:\n");
Debug.fcall(Flags.DUMP_INDX_DAE, BackendDump.dumpBackendDAE, outSODE);
Debug.fcall2(Flags.DUMP_INDX_DAE, BackendDump.dumpBackendDAE, outSODE, "dumpindxdae");
Debug.fcall(Flags.DUMP_BACKENDDAE_INFO, BackendDump.dumpCompShort, outSODE);
Debug.fcall(Flags.DUMP_EQNINORDER, BackendDump.dumpEqnsSolved, outSODE);
checkBackendDAEWithErrorMsg(outSODE);
Expand Down Expand Up @@ -7925,7 +7923,7 @@ algorithm
dae = optModule(inDAE);
Debug.execStat("preOpt " +& moduleStr,BackendDAE.RT_CLOCK_EXECSTAT_BACKEND_MODULES);
Debug.fcall(Flags.OPT_DAE_DUMP, print, stringAppendList({"\nOptimisation Module ",moduleStr,":\n\n"}));
Debug.fcall(Flags.OPT_DAE_DUMP, BackendDump.dumpBackendDAE, dae);
Debug.fcall(Flags.OPT_DAE_DUMP, BackendDump.printBackendDAE, dae);
(dae1,status) = preoptimiseDAE(dae,rest);
then (dae1,status);
case (_,(optModule,moduleStr,b)::rest)
Expand Down Expand Up @@ -8151,7 +8149,7 @@ algorithm
dae = optModule(inDAE);
Debug.execStat("pastOpt " +& moduleStr,BackendDAE.RT_CLOCK_EXECSTAT_BACKEND_MODULES);
Debug.fcall(Flags.OPT_DAE_DUMP, print, stringAppendList({"\nOptimisation Module ",moduleStr,":\n\n"}));
Debug.fcall(Flags.OPT_DAE_DUMP, BackendDump.dumpBackendDAE, dae);
Debug.fcall(Flags.OPT_DAE_DUMP, BackendDump.printBackendDAE, dae);
dae1 = causalizeDAE(dae,NONE(),matchingAlgorithm,daeHandler,false);
(dae2,status) = pastoptimiseDAE(dae1,rest,matchingAlgorithm,daeHandler);
then
Expand Down Expand Up @@ -8260,8 +8258,7 @@ algorithm
matchingAlgorithm := getMatchingAlgorithm(strmatchingAlgorithm);
daeHandler := getIndexReductionMethod(strdaeHandler);

Debug.fcall(Flags.DUMP_DAE_LOW, print, "dumpdaelow:\n");
Debug.fcall(Flags.DUMP_DAE_LOW, BackendDump.dumpBackendDAE, inDAE);
Debug.fcall2(Flags.DUMP_DAE_LOW, BackendDump.dumpBackendDAE, inDAE, "dumpdaelow");
// pre optimisation phase
_ := traverseBackendDAEExps(inDAE,ExpressionSimplify.simplifyTraverseHelper,0) "simplify all expressions";
(optdae,Util.SUCCESS()) := preoptimiseDAE(inDAE,preOptModules);
Expand All @@ -8274,8 +8271,7 @@ algorithm
(outSODE,Util.SUCCESS()) := pastoptimiseDAE(sode,pastOptModules,matchingAlgorithm,daeHandler);
_ := traverseBackendDAEExps(outSODE,ExpressionSimplify.simplifyTraverseHelper,0) "simplify all expressions";

Debug.fcall(Flags.DUMP_INDX_DAE, print, "dumpindxdae:\n");
Debug.fcall(Flags.DUMP_INDX_DAE, BackendDump.dumpBackendDAE, outSODE);
Debug.fcall2(Flags.DUMP_INDX_DAE, BackendDump.dumpBackendDAE, outSODE, "dumpindxdae");
Debug.fcall(Flags.DUMP_BACKENDDAE_INFO, BackendDump.dumpCompShort, outSODE);
Debug.fcall(Flags.DUMP_EQNINORDER, BackendDump.dumpEqnsSolved, outSODE);
end getSolvedSystemforJacobians;
Expand Down
22 changes: 16 additions & 6 deletions Compiler/BackEnd/BackendDump.mo
View file
Expand Up @@ -77,6 +77,19 @@ protected import Util;
// These are functions, that print directly to the standard-stream.
// =============================================================================

public function printBackendDAE "function printBackendDAE
This function dumps the BackendDAE.BackendDAE representaton to stdout."
input BackendDAE.BackendDAE inBackendDAE;
protected
BackendDAE.EqSystems eqs;
BackendDAE.Shared shared;
algorithm
BackendDAE.DAE(eqs,shared) := inBackendDAE;
List.map_0(eqs,dumpEqSystem);
print("\n");
dumpShared(shared);
end printBackendDAE;

public function printVariables "function printVariables
Helper function to dump."
input BackendDAE.Variables vars;
Expand Down Expand Up @@ -237,14 +250,11 @@ end dumpEquationArray;
public function dumpBackendDAE "function dumpBackendDAE
This function dumps the BackendDAE.BackendDAE representaton to stdout."
input BackendDAE.BackendDAE inBackendDAE;
protected
BackendDAE.EqSystems eqs;
BackendDAE.Shared shared;
input String heading;
algorithm
BackendDAE.DAE(eqs,shared) := inBackendDAE;
List.map_0(eqs,dumpEqSystem);
print("\n########################################\n" +& heading +& "\n########################################\n\n");
printBackendDAE(inBackendDAE);
print("\n");
dumpShared(shared);
end dumpBackendDAE;

public function dumpSparsityPattern "function dumpSparsityPattern
Expand Down
4 changes: 2 additions & 2 deletions Compiler/BackEnd/IndexReduction.mo
View file
Expand Up @@ -417,7 +417,7 @@ algorithm
eqns1 = List.uniqueIntN(eqns1,arrayLength(mapIncRowEqn));
Debug.fcall(Flags.BLT_DUMP, print, BackendDump.dumpMarkedEqns(isyst, eqns1));
syst = BackendDAEUtil.setEqSystemMatching(isyst,BackendDAE.MATCHING(inAssignments1,inAssignments2,{}));
Debug.fcall(Flags.BLT_DUMP, BackendDump.dumpBackendDAE, BackendDAE.DAE({syst},ishared));
Debug.fcall(Flags.BLT_DUMP, BackendDump.printBackendDAE, BackendDAE.DAE({syst},ishared));
Error.addMessage(Error.INTERNAL_ERROR, {"IndexReduction.pantelidesIndexReduction failed! Found empty set of continues equations. Use +d=bltdump to get more information."});
then
fail();
Expand All @@ -433,7 +433,7 @@ algorithm
varlst = List.map1r(unassignedStates,BackendVariable.getVarAt,BackendVariable.daeVars(isyst));
Debug.fcall(Flags.BLT_DUMP, BackendDump.printVarList,varlst);
syst = BackendDAEUtil.setEqSystemMatching(isyst,BackendDAE.MATCHING(inAssignments1,inAssignments2,{}));
Debug.fcall(Flags.BLT_DUMP, BackendDump.dumpBackendDAE, BackendDAE.DAE({syst},ishared));
Debug.fcall(Flags.BLT_DUMP, BackendDump.printBackendDAE, BackendDAE.DAE({syst},ishared));
Error.addMessage(Error.INTERNAL_ERROR, {"IndexReduction.pantelidesIndexReduction1 failed! System is structurally singulare and cannot handled because number of unassigned equations is larger than number of states. Use +d=bltdump to get more information."});
then
fail();
Expand Down
9 changes: 3 additions & 6 deletions Compiler/BackEnd/Initialization.mo
View file
Expand Up @@ -74,8 +74,7 @@ protected
algorithm
dae := inlineWhenForInitialization(inDAE);

Debug.fcall(Flags.DUMP_INITIAL_SYSTEM, print, "\ninlineWhenForInitialization\n########################################\n\n");
Debug.fcall(Flags.DUMP_INITIAL_SYSTEM, BackendDump.dumpBackendDAE, dae);
Debug.fcall2(Flags.DUMP_INITIAL_SYSTEM, BackendDump.dumpBackendDAE, dae, "inlineWhenForInitialization");

initVars := selectInitializationVariablesDAE(dae);
Debug.fcall2(Flags.DUMP_INITIAL_SYSTEM, BackendDump.dumpVariables, initVars, "selected initialization variables");
Expand Down Expand Up @@ -426,8 +425,7 @@ algorithm
{}));

// some debug prints
Debug.fcall(Flags.DUMP_INITIAL_SYSTEM, print, "\ninitial system\n########################################\n\n");
Debug.fcall(Flags.DUMP_INITIAL_SYSTEM, BackendDump.dumpBackendDAE, initdae);
Debug.fcall2(Flags.DUMP_INITIAL_SYSTEM, BackendDump.dumpBackendDAE, initdae, "initial system");

// now let's solve the system!
initdae = solveInitialSystem2(vars, eqns, inDAE, initdae);
Expand Down Expand Up @@ -482,8 +480,7 @@ algorithm
// simplify system
(isyst, Util.SUCCESS()) = BackendDAEUtil.pastoptimiseDAE(isyst, pastOptModules, matchingAlgorithm, daeHandler);

Debug.fcall(Flags.DUMP_INITIAL_SYSTEM, print, "\nsolved initial system\n########################################\n\n");
Debug.fcall(Flags.DUMP_INITIAL_SYSTEM, BackendDump.dumpBackendDAE, isyst);
Debug.fcall2(Flags.DUMP_INITIAL_SYSTEM, BackendDump.dumpBackendDAE, isyst, "solved initial system");
then isyst;

// under-determined system
Expand Down
6 changes: 2 additions & 4 deletions Compiler/BackEnd/InlineSolver.mo
View file
Expand Up @@ -76,13 +76,11 @@ algorithm
case dae equation
true = Flags.isSet(Flags.INLINE_SOLVER);

Debug.fcall(Flags.DUMP_INLINE_SOLVER, print, "\ninlineSolver (1)\n########################################\n\n");
Debug.fcall(Flags.DUMP_INLINE_SOLVER, BackendDump.dumpBackendDAE, dae);
Debug.fcall2(Flags.DUMP_INLINE_SOLVER, BackendDump.dumpBackendDAE, dae, "inlineSolver (1)");

dae = replaceStates(dae);

Debug.fcall(Flags.DUMP_INLINE_SOLVER, print, "\ninlineSolver (2)\n########################################\n\n");
Debug.fcall(Flags.DUMP_INLINE_SOLVER, BackendDump.dumpBackendDAE, dae);
Debug.fcall2(Flags.DUMP_INLINE_SOLVER, BackendDump.dumpBackendDAE, dae, "inlineSolver (2)");
then NONE();

else equation
Expand Down
3 changes: 1 addition & 2 deletions Compiler/BackEnd/SimCodeUtil.mo
View file
Expand Up @@ -6718,8 +6718,7 @@ algorithm
//mT = BackendDAEUtil.transposeMatrix(m);
v1 = listArray(lv1);
v2 = listArray(lv2);
Debug.fcall(Flags.PARAM_DLOW_DUMP, print, "Param DAE:\n");
Debug.fcall(Flags.PARAM_DLOW_DUMP, BackendDump.dumpBackendDAE, paramdlow);
Debug.fcall2(Flags.PARAM_DLOW_DUMP, BackendDump.dumpBackendDAE, paramdlow, "Param DAE");
Debug.fcall(Flags.PARAM_DLOW_DUMP, BackendDump.dumpIncidenceMatrix, m);
Debug.fcall(Flags.PARAM_DLOW_DUMP, BackendDump.dumpIncidenceMatrixT, mT);
Debug.fcall(Flags.PARAM_DLOW_DUMP, BackendDump.dumpMatching, v1);
Expand Down
6 changes: 4 additions & 2 deletions Compiler/Template/CodegenC.tpl
View file
Expand Up @@ -2301,9 +2301,11 @@ case SES_NONLINEAR(__) then
let namestr = cref(name)
<<
data->simulationInfo.nonlinearSystemData[<%nonlinindx%>].nlsx[<%i0%>] = <%namestr%>;
data->simulationInfo.nonlinearSystemData[<%nonlinindx%>].nlsxScaling[<%i0%>] = $P$ATTRIBUTE<%namestr%>.nominal;
data->simulationInfo.nonlinearSystemData[<%nonlinindx%>].nlsxExtrapolation[<%i0%>] = extraPolate(<%namestr%>, _<%namestr%>(1) /*old*/, _<%namestr%>(2) /*old2*/);
data->simulationInfo.nonlinearSystemData[<%nonlinindx%>].nlsxOld[<%i0%>] = _<%namestr%>(1) /*old*/;
data->simulationInfo.nonlinearSystemData[<%nonlinindx%>].nlsxExtrapolation[<%i0%>] = extraPolate(<%namestr%>, _<%namestr%>(1) /*old*/, _<%namestr%>(2) /*old2*/);
data->simulationInfo.nonlinearSystemData[<%nonlinindx%>].nominal[<%i0%>] = $P$ATTRIBUTE<%namestr%>.nominal;
data->simulationInfo.nonlinearSystemData[<%nonlinindx%>].min[<%i0%>] = $P$ATTRIBUTE<%namestr%>.min;
data->simulationInfo.nonlinearSystemData[<%nonlinindx%>].max[<%i0%>] = $P$ATTRIBUTE<%namestr%>.max;
>>
;separator="\n"%>
solve_nonlinear_system(data, <%indexNonLinear%>);
Expand Down
60 changes: 45 additions & 15 deletions SimulationRuntime/c/simulation/simulation_runtime.cpp
View file
Expand Up @@ -121,25 +121,30 @@ void setTermMsg(const char *msg)
{
size_t i;
size_t length = strlen(msg);
if(length > 0) {
if(TermMsg == NULL) {
if(length > 0)
{
if(TermMsg == NULL)
{
TermMsg = (char*)malloc((length+1)*sizeof(char));
}
else
{
if(strlen(msg) > strlen(TermMsg))
{
if(TermMsg != NULL)
{
free(TermMsg);
}
TermMsg = (char*)malloc((length+1)*sizeof(char));
} else {
if(strlen(msg) > strlen(TermMsg)) {
if(TermMsg != NULL) {
free(TermMsg);
}
TermMsg = (char*)malloc((length+1)*sizeof(char));
}
}
for(i=0;i<length;i++)
TermMsg[i] = msg[i];
/* set the terminating 0 */
TermMsg[i] = '\0';
}
for(i=0;i<length;i++)
TermMsg[i] = msg[i];
/* set the terminating 0 */
TermMsg[i] = '\0';
}
}


/* \brief determine verboselevel by investigating flag -lv flags
*
* Flags are or'ed to a returnvalue.
Expand Down Expand Up @@ -223,6 +228,27 @@ void setGlobalVerboseLevel(int argc, char**argv)
delete flags;
}

int getNonlinearSolverMethod(int argc, char**argv)
{
const string *method = getOption("nls", argc, argv);

if(!method)
return NS_HYBRID; /* default method */

if(*method == string("hybrid"))
return NS_HYBRID;
else if(*method == string("kinsol"))
return NS_KINSOL;

WARNING1(LOG_STDOUT, "unrecognized option -nls %s", method->c_str());
WARNING(LOG_STDOUT, "current options are:");
INDENT(LOG_STDOUT);
WARNING2(LOG_STDOUT, "%-18s [%s]", "hybrid", "default method");
WARNING2(LOG_STDOUT, "%-18s [%s]", "kinsol", "sundials/kinsol");
THROW("see last warning");
return NS_NONE;
}

/**
* Signals the type of the simulation
* retuns true for interactive and false for non-interactive
Expand Down Expand Up @@ -532,6 +558,8 @@ int initRuntimeAndSimulation(int argc, char**argv, DATA *data)
INFO(LOG_STDOUT, "\tspecify a new result file than the default Model_res.mat");
INFO(LOG_STDOUT, "<-m|s solver:{dassl,euler,rungekutta,inline-euler,inline-rungekutta,qss}>");
INFO(LOG_STDOUT, "\tspecify the solver");
INFO(LOG_STDOUT, "<-nls={hybrid|kinsol}>");
INFO(LOG_STDOUT, "\tspecify the nonlinear solver");
INFO(LOG_STDOUT, "<-interactive> <-port value>");
INFO(LOG_STDOUT, "\tspecify interactive simulation and port");
INFO(LOG_STDOUT, "<-iim initialization method:{none,numeric,symbolic}>");
Expand Down Expand Up @@ -573,6 +601,7 @@ int initRuntimeAndSimulation(int argc, char**argv, DATA *data)
EXIT(0);
}

setGlobalVerboseLevel(argc, argv);
initializeDataStruc(data);
if(!data)
{
Expand Down Expand Up @@ -619,7 +648,7 @@ int initRuntimeAndSimulation(int argc, char**argv, DATA *data)
}
#endif

setGlobalVerboseLevel(argc, argv);
data->simulationInfo.nlsMethod = getNonlinearSolverMethod(argc, argv);

return 0;
}
Expand Down Expand Up @@ -670,6 +699,7 @@ void communicateStatus(const char *phase, double completionPercent /*0.0 to 1.0*
int _main_SimulationRuntime(int argc, char**argv, DATA *data)
{
int retVal = -1;

if(!setjmp(globalJmpbuf))
{
if(initRuntimeAndSimulation(argc, argv, data)) //initRuntimeAndSimulation returns 1 if an error occurs
Expand Down
2 changes: 2 additions & 0 deletions SimulationRuntime/c/simulation/solver/CMakeLists.txt
View file
Expand Up @@ -10,13 +10,15 @@ SET(solver_sources dassl.c
model_help.c
nonlinearSystem.c
nonlinearSolverHybrd.c
kinsolSolver.c
solver_main.c)

SET(solver_headers dassl.h
delay.h
events.h
model_help.h
nonlinearSolverHybrd.h
kinsolSolver.h
nonlinearSystem.h
solver_main.h)

Expand Down

0 comments on commit 19fc2f3

Please sign in to comment.