Skip to content

Commit

Permalink
- fix non-linear system solver for iteration variables with bounded a…
Browse files Browse the repository at this point in the history 
…ttributes (#2841)

git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@22519 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
@lochel
lochel committed Sep 30, 2014
1 parent 7d61c46 commit 8f3b64b
Show file tree
Hide file tree
Showing 8 changed files with 148 additions and 65 deletions.
20 changes: 8 additions & 12 deletions SimulationRuntime/c/simulation/simulation_runtime.cpp
View file
Expand Up @@ -849,15 +849,11 @@ int initRuntimeAndSimulation(int argc, char**argv, DATA *data)
rt_tick(SIM_TIMER_INIT_XML);
read_input_xml(&(data->modelData), &(data->simulationInfo));
rt_accumulate(SIM_TIMER_INIT_XML);

/* allocate memory for mixed system solvers */
allocatemixedSystem(data);

/* allocate memory for linear system solvers */
allocatelinearSystem(data);

/* allocate memory for non-linear system solvers */
allocateNonlinearSystem(data);

/* initialize static data of mixed/linear/non-linear system solvers */
initializeMixedSystems(data);
initializeLinearSystems(data);
initializeNonlinearSystems(data);

// this sets the static variable that is in the file with the generated-model functions
if(data->modelData.nVariablesReal == 0 && data->modelData.nVariablesInteger && data->modelData.nVariablesBoolean)
Expand Down Expand Up @@ -978,9 +974,9 @@ int _main_SimulationRuntime(int argc, char**argv, DATA *data)
retVal = startNonInteractiveSimulation(argc, argv, data);
}

freemixedSystem(data); /* free mixed system data */
freelinearSystem(data); /* free linear system data */
freeNonlinearSystem(data); /* free nonlinear system data */
freeMixedSystems(data); /* free mixed system data */
freeLinearSystems(data); /* free linear system data */
freeNonlinearSystems(data); /* free nonlinear system data */

data->callback->callExternalObjectDestructors(data);
deInitializeDataStruc(data);
Expand Down
75 changes: 52 additions & 23 deletions SimulationRuntime/c/simulation/solver/initialization/initialization.c
View file
Expand Up @@ -977,8 +977,8 @@ static int importStartValues(DATA *data, const char *pInitFile, double initTime)
*/
int initialization(DATA *data, const char* pInitMethod, const char* pOptiMethod, const char* pInitFile, double initTime, int lambda_steps)
{
int initMethod = data->callback->useSymbolicInitialization ? IIM_SYMBOLIC : IIM_NUMERIC; /* default method */
int optiMethod = IOM_NELDER_MEAD_EX; /* default method */
int initMethod = data->callback->useSymbolicInitialization ? IIM_SYMBOLIC : IIM_NUMERIC; /* default method */
int optiMethod = IOM_NELDER_MEAD_EX; /* default method */
int retVal = -1;
int i;

Expand All @@ -1003,47 +1003,61 @@ int initialization(DATA *data, const char* pInitMethod, const char* pOptiMethod,
/* set up all variables with their start-values */
setAllVarsToStart(data);

/* update static data of non-linear system solvers */
updateStaticDataOfNonlinearSystems(data);

/* if there are user-specified options, use them! */
if(pInitMethod && strcmp(pInitMethod, "")) {
if(pInitMethod && strcmp(pInitMethod, ""))
{
initMethod = IIM_UNKNOWN;

for(i=1; i<IIM_MAX; ++i) {
if(!strcmp(pInitMethod, INIT_METHOD_NAME[i])) {
for(i=1; i<IIM_MAX; ++i)
{
if(!strcmp(pInitMethod, INIT_METHOD_NAME[i]))
{
initMethod = i;
}
}

if(initMethod == IIM_UNKNOWN) {
if(initMethod == IIM_UNKNOWN)
{
warningStreamPrint(LOG_STDOUT, 0, "unrecognized option -iim %s", pInitMethod);
warningStreamPrint(LOG_STDOUT, 0, "current options are:");
for(i=1; i<IIM_MAX; ++i) {
for(i=1; i<IIM_MAX; ++i)
{
warningStreamPrint(LOG_STDOUT, 0, "| %-15s [%s]", INIT_METHOD_NAME[i], INIT_METHOD_DESC[i]);
}
throwStreamPrint(data->threadData, "see last warning");
}
}

if(pOptiMethod && strcmp(pOptiMethod, "")) {
if(pOptiMethod && strcmp(pOptiMethod, ""))
{
optiMethod = IOM_UNKNOWN;

for(i=1; i<IOM_MAX; ++i) {
if(!strcmp(pOptiMethod, OPTI_METHOD_NAME[i])) {
for(i=1; i<IOM_MAX; ++i)
{
if(!strcmp(pOptiMethod, OPTI_METHOD_NAME[i]))
{
optiMethod = i;
}
}

if(optiMethod == IOM_UNKNOWN) {
if(optiMethod == IOM_UNKNOWN)
{
warningStreamPrint(LOG_STDOUT, 0, "unrecognized option -iom %s", pOptiMethod);
warningStreamPrint(LOG_STDOUT, 0, "current options are:");
for(i=1; i<IOM_MAX; ++i) {
for(i=1; i<IOM_MAX; ++i)
{
warningStreamPrint(LOG_STDOUT, 0, "| %-15s [%s]", OPTI_METHOD_NAME[i], OPTI_METHOD_DESC[i]);
}
throwStreamPrint(data->threadData, "see last warning");
}
}

infoStreamPrint(LOG_INIT, 0, "initialization method: %-15s [%s]", INIT_METHOD_NAME[initMethod], INIT_METHOD_DESC[initMethod]);
if(initMethod == IIM_NUMERIC) {
if(initMethod == IIM_NUMERIC)
{
infoStreamPrint(LOG_INIT, 0, "optimization method: %-15s [%s]", OPTI_METHOD_NAME[optiMethod], OPTI_METHOD_DESC[optiMethod]);
}

Expand All @@ -1053,36 +1067,51 @@ int initialization(DATA *data, const char* pInitMethod, const char* pOptiMethod,
/* initialize all (nonlinear|linear|mixed) systems
* This is a workaround and should be removed as soon as possible.
*/
for(i=0; i<data->modelData.nNonLinearSystems; ++i) {
for(i=0; i<data->modelData.nNonLinearSystems; ++i)
{
data->simulationInfo.nonlinearSystemData[i].solved = 1;
}
for(i=0; i<data->modelData.nLinearSystems; ++i) {
for(i=0; i<data->modelData.nLinearSystems; ++i)
{
data->simulationInfo.linearSystemData[i].solved = 1;
}
for(i=0; i<data->modelData.nMixedSystems; ++i) {
for(i=0; i<data->modelData.nMixedSystems; ++i)
{
data->simulationInfo.mixedSystemData[i].solved = 1;
}
/* end workaround */

/* select the right initialization-method */
if(initMethod == IIM_NONE) {
if(IIM_NONE == initMethod)
{
retVal = 0;
} else if(initMethod == IIM_NUMERIC) {
}
else if(IIM_NUMERIC == initMethod)
{
retVal = numeric_initialization(data, optiMethod, lambda_steps);
} else if(initMethod == IIM_SYMBOLIC) {
}
else if(IIM_SYMBOLIC == initMethod)
{
retVal = symbolic_initialization(data, lambda_steps);
} else {
}
else
{
throwStreamPrint(data->threadData, "unsupported option -iim");
}

/* check for unsolved (nonlinear|linear|mixed) systems
* This is a workaround and should be removed as soon as possible.
*/
if(check_nonlinear_solutions(data, 1)) {
if(check_nonlinear_solutions(data, 1))
{
retVal = -2;
} else if(check_linear_solutions(data, 1)) {
}
else if(check_linear_solutions(data, 1))
{
retVal = -3;
} else if(check_mixed_solutions(data, 1)) {
}
else if(check_mixed_solutions(data, 1))
{
retVal = -4;
}
/* end workaround */
Expand Down
52 changes: 36 additions & 16 deletions SimulationRuntime/c/simulation/solver/linearSystem.c
View file
Expand Up @@ -41,18 +41,20 @@
#include "blaswrap.h"
#include "f2c.h"

/*! \fn int allocatelinearSystem(DATA *data)
/*! \fn int initializeLinearSystems(DATA *data)
*
* This function allocates memory for all linear systems.
*
* \param [ref] [data]
*/
int allocatelinearSystem(DATA *data)
int initializeLinearSystems(DATA *data)
{
int i,nnz;
int i, nnz;
int size;
LINEAR_SYSTEM_DATA *linsys = data->simulationInfo.linearSystemData;

infoStreamPrint(LOG_LS, 1, "initialize linear system solvers");

for(i=0; i<data->modelData.nLinearSystems; ++i)
{
size = linsys[i].size;
Expand All @@ -63,8 +65,10 @@ int allocatelinearSystem(DATA *data)
linsys[i].b = (double*) malloc(size*sizeof(double));

/* check if analytical jacobian is created */
if (1 == linsys[i].method){
if(linsys[i].jacobianIndex != -1){
if (1 == linsys[i].method)
{
if(linsys[i].jacobianIndex != -1)
{
assertStreamPrint(data->threadData, 0 != linsys[i].analyticalJacobianColumn, "jacobian function pointer is invalid" );
}
if(linsys[i].initialAnalyticalJacobian(data))
Expand All @@ -74,62 +78,73 @@ int allocatelinearSystem(DATA *data)
}
/* allocate solver data */
/* the implementation of matrix A is solver-specific */
switch(data->simulationInfo.lsMethod){
switch(data->simulationInfo.lsMethod)
{
case LS_LAPACK:
linsys[i].A = (double*) malloc(size*size*sizeof(double));
linsys[i].setAElement = setAElementLAPACK;
allocateLapackData(size, &linsys[i].solverData);
break;

case LS_LIS:
linsys[i].setAElement = setAElementLis;
allocateLisData(size, size, nnz, &linsys[i].solverData);
break;

default:
throwStreamPrint(data->threadData, "unrecognized linear solver");
}
}

messageClose(LOG_LS);
return 0;
}

/*! \fn freelinearSystem
/*! \fn freeLinearSystems
*
* This function frees memory of linear systems.
*
* \param [ref] [data]
*/
int freelinearSystem(DATA *data)
int freeLinearSystems(DATA *data)
{
int i;
LINEAR_SYSTEM_DATA* linsys = data->simulationInfo.linearSystemData;

for(i=0;i<data->modelData.nLinearSystems;++i)
infoStreamPrint(LOG_LS, 1, "free linear system solvers");

for(i=0; i<data->modelData.nLinearSystems; ++i)
{
/* free system and solver data */
free(linsys[i].x);
free(linsys[i].b);

switch(data->simulationInfo.lsMethod){
switch(data->simulationInfo.lsMethod)
{
case LS_LAPACK:
freeLapackData(&linsys[i].solverData);
free(linsys[i].A);
break;

case LS_LIS:
freeLisData(&linsys[i].solverData);
break;

default:
throwStreamPrint(data->threadData, "unrecognized linear solver");
}

free(linsys[i].solverData);
}

messageClose(LOG_LS);
return 0;
}

/*! \fn solve non-linear systems
*
* \param [in] [data]
* [sysNumber] index of corresponding non-linear System
* [in] [sysNumber] index of corresponding non-linear System
*
* \author wbraun
*/
Expand All @@ -138,13 +153,16 @@ int solve_linear_system(DATA *data, int sysNumber)
int success;
LINEAR_SYSTEM_DATA* linsys = data->simulationInfo.linearSystemData;

switch(data->simulationInfo.lsMethod){
switch(data->simulationInfo.lsMethod)
{
case LS_LAPACK:
success = solveLapack(data, sysNumber);
break;

case LS_LIS:
success = solveLis(data, sysNumber);
break;

default:
throwStreamPrint(data->threadData, "unrecognized linear solver");
}
Expand All @@ -168,27 +186,29 @@ int check_linear_solutions(DATA *data, int printFailingSystems)
int i, j, retVal=0;

for(i=0; i<data->modelData.nLinearSystems; ++i)
if(linsys[i].solved == 0)
{
if(0 == linsys[i].solved)
{
retVal = 1;
if(printFailingSystems && ACTIVE_WARNING_STREAM(LOG_LS))
{
int indexes[2] = {1,modelInfoXmlGetEquation(&data->modelData.modelDataXml, linsys->equationIndex).id};
int indexes[2] = {1, modelInfoXmlGetEquation(&data->modelData.modelDataXml, linsys->equationIndex).id};
warningStreamPrintWithEquationIndexes(LOG_LS, 1, indexes, "linear system %d fails at t=%g", indexes[1], data->localData[0]->timeValue);
messageClose(LOG_LS);
}
}
}

return retVal;
}

void setAElementLAPACK(int row, int col, double value, int nth, void *data )
void setAElementLAPACK(int row, int col, double value, int nth, void *data)
{
LINEAR_SYSTEM_DATA* linsys = (LINEAR_SYSTEM_DATA*) data;
linsys->A[row + col * linsys->size] = value;
}

void setAElementLis(int row, int col, double value, int nth, void *data )
void setAElementLis(int row, int col, double value, int nth, void *data)
{
LINEAR_SYSTEM_DATA* linSys = (LINEAR_SYSTEM_DATA*) data;
DATA_LIS* sData = (DATA_LIS*) linSys->solverData;
Expand Down
4 changes: 2 additions & 2 deletions SimulationRuntime/c/simulation/solver/linearSystem.h
View file
Expand Up @@ -55,8 +55,8 @@ enum LINEAR_SOLVER

typedef void* LS_SOLVER_DATA;

int allocatelinearSystem(DATA *data);
int freelinearSystem(DATA *data);
int initializeLinearSystems(DATA *data);
int freeLinearSystems(DATA *data);
int solve_linear_system(DATA *data, int sysNumber);
int check_linear_solutions(DATA *data, int printFailingSystems);

Expand Down

0 comments on commit 8f3b64b

Please sign in to comment.