Skip to content

Commit

Permalink
adding context depending extrapolation
Browse files Browse the repository at this point in the history
  • Loading branch information
@OpenModelica-Hudson
Willi Braun authored and OpenModelica-Hudson committed Nov 19, 2015
1 parent 920eedf commit 57d05fd
Show file tree
Hide file tree
Showing 12 changed files with 254 additions and 113 deletions.
66 changes: 21 additions & 45 deletions SimulationRuntime/c/simulation/solver/dassl.c
View file
Expand Up @@ -128,34 +128,6 @@ dummy_precondition(int *neq, double *t, double *y, double *yprime, double *savr,

static int continue_DASSL(int* idid, double* tolarence);

enum EVAL_CONTEXT
{
CONTEXT_UNKNOWN = 0,

CONTEXT_ODE,
CONTEXT_JACOBIAN,
CONTEXT_EVENTS,

CONTEXT_MAX
};

const char *context_string[CONTEXT_MAX] = {
"context UNKNOWN",
"context ODE evaluation",
"context Jacobian",
"context Event Search"
};

void setDasslContext(DASSL_DATA* dasslData, double* currentTime, int currentContext){
dasslData->currentContextOld = dasslData->currentContext;
dasslData->currentContext = currentContext;
infoStreamPrint(LOG_DASSL, 0, "+++ Set DASSL %s +++ at time %f", context_string[dasslData->currentContext], *currentTime);
}
void unsetDasslContext(DASSL_DATA* dasslData){
infoStreamPrint(LOG_DASSL, 0, "--- Unset DASSL %s ---", context_string[dasslData->currentContext]);
dasslData->currentContext = dasslData->currentContextOld;
}

/* function for calculating state values on residual form */
static int functionODE_residual(double *t, double *x, double *xprime, double *cj, double *delta, int *ires, double *rpar, int* ipar);
/* function for calculating zeroCrossings */
Expand Down Expand Up @@ -200,7 +172,7 @@ int dassl_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo,
dasslData->newdelta = (double*) malloc(data->modelData.nStates*sizeof(double));
dasslData->stateDer = (double*) malloc(data->modelData.nStates*sizeof(double));

dasslData->currentContext = CONTEXT_UNKNOWN;
data->simulationInfo.currentContext = CONTEXT_ALGEBRAIC;

/* setup internal ring buffer for dassl */

Expand Down Expand Up @@ -390,9 +362,11 @@ int dassl_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo,
/* set up the appropriate function pointer */
switch (dasslData->dasslJacobian){
case DASSL_COLOREDNUMJAC:
data->simulationInfo.jacobianEvals = data->simulationInfo.analyticJacobians[data->callback->INDEX_JAC_A].sparsePattern.maxColors;
dasslData->jacobianFunction = JacobianOwnNumColored;
break;
case DASSL_COLOREDSYMJAC:
data->simulationInfo.jacobianEvals = data->simulationInfo.analyticJacobians[data->callback->INDEX_JAC_A].sparsePattern.maxColors;
dasslData->jacobianFunction = JacobianSymbolicColored;
break;
case DASSL_SYMJAC:
Expand Down Expand Up @@ -833,9 +807,9 @@ int functionODE_residual(double *t, double *y, double *yd, double* cj, double *d
int saveJumpState;
int success = 0;

if (dasslData->currentContext == CONTEXT_UNKNOWN)
if (data->simulationInfo.currentContext == CONTEXT_ALGEBRAIC)
{
setDasslContext(dasslData, t, CONTEXT_ODE);
setContext(data, t, CONTEXT_ODE);
}
printCurrentStatesVector(LOG_DASSL_STATES, y, data, *t);

Expand Down Expand Up @@ -876,8 +850,8 @@ int functionODE_residual(double *t, double *y, double *yd, double* cj, double *d

data->localData[0]->timeValue = timeBackup;

if (dasslData->currentContext == CONTEXT_ODE){
unsetDasslContext(dasslData);
if (data->simulationInfo.currentContext == CONTEXT_ODE){
unsetContext(data);
}

TRACE_POP
Expand All @@ -895,9 +869,9 @@ int function_ZeroCrossingsDASSL(int *neqm, double *t, double *y, double *yp,
double timeBackup;
int saveJumpState;

if (dasslData->currentContext == CONTEXT_UNKNOWN)
if (data->simulationInfo.currentContext == CONTEXT_ALGEBRAIC)
{
setDasslContext(dasslData, t, CONTEXT_EVENTS);
setContext(data, t, CONTEXT_EVENTS);
}

saveJumpState = threadData->currentErrorStage;
Expand All @@ -917,8 +891,8 @@ int function_ZeroCrossingsDASSL(int *neqm, double *t, double *y, double *yp,
threadData->currentErrorStage = saveJumpState;
data->localData[0]->timeValue = timeBackup;

if (dasslData->currentContext == CONTEXT_EVENTS){
unsetDasslContext(dasslData);
if (data->simulationInfo.currentContext == CONTEXT_EVENTS){
unsetContext(data);
}

TRACE_POP
Expand Down Expand Up @@ -1052,7 +1026,7 @@ static int JacobianSymbolicColored(double *t, double *y, double *yprime, double
int i;
int j;

setDasslContext(dasslData, t, CONTEXT_JACOBIAN);
setContext(data, t, CONTEXT_JACOBIAN);

backupStates = data->localData[0]->realVars;
timeBackup = data->localData[0]->timeValue;
Expand All @@ -1076,7 +1050,7 @@ static int JacobianSymbolicColored(double *t, double *y, double *yprime, double
data->localData[0]->realVars = backupStates;
data->localData[0]->timeValue = timeBackup;

unsetDasslContext(dasslData);
unsetContext(data);

TRACE_POP
return 0;
Expand All @@ -1097,7 +1071,7 @@ static int JacobianSymbolic(double *t, double *y, double *yprime, double *deltaD
int i;
int j;

setDasslContext(dasslData, t, CONTEXT_JACOBIAN);
setContext(data, t, CONTEXT_JACOBIAN);

backupStates = data->localData[0]->realVars;
timeBackup = data->localData[0]->timeValue;
Expand All @@ -1121,7 +1095,7 @@ static int JacobianSymbolic(double *t, double *y, double *yprime, double *deltaD
data->localData[0]->realVars = backupStates;
data->localData[0]->timeValue = timeBackup;

unsetDasslContext(dasslData);
unsetContext(data);

TRACE_POP
return 0;
Expand Down Expand Up @@ -1198,7 +1172,7 @@ static int JacobianOwnNum(double *t, double *y, double *yprime, double *deltaD,
DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];

setDasslContext(dasslData, t, CONTEXT_JACOBIAN);
setContext(data, t, CONTEXT_JACOBIAN);

if(jacA_num(data, t, y, yprime, deltaD, pd, cj, h, wt, rpar, ipar))
{
Expand All @@ -1213,7 +1187,7 @@ static int JacobianOwnNum(double *t, double *y, double *yprime, double *deltaD,
pd[j] -= (double) *cj;
j += data->modelData.nStates + 1;
}
unsetDasslContext(dasslData);
unsetContext(data);

TRACE_POP
return 0;
Expand Down Expand Up @@ -1257,6 +1231,8 @@ int jacA_numColored(DATA* data, double *t, double *y, double *yprime, double *de

functionODE_residual(t, y, yprime, cj, dasslData->newdelta, &ires, rpar, ipar);

increaseJacContext(data);

for(ii = 0; ii < data->simulationInfo.analyticJacobians[index].sizeCols; ii++)
{
if(data->simulationInfo.analyticJacobians[index].sparsePattern.colorCols[ii]-1 == i)
Expand Down Expand Up @@ -1308,7 +1284,7 @@ static int JacobianOwnNumColored(double *t, double *y, double *yprime, double *d
threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];
int i,j;

setDasslContext(dasslData, t, CONTEXT_JACOBIAN);
setContext(data, t, CONTEXT_JACOBIAN);

if(jacA_numColored(data, t, y, yprime, deltaD, pd, cj, h, wt, rpar, ipar))
{
Expand All @@ -1324,7 +1300,7 @@ static int JacobianOwnNumColored(double *t, double *y, double *yprime, double *d
pd[j] -= (double) *cj;
j += data->modelData.nStates + 1;
}
unsetDasslContext(dasslData);
unsetContext(data);

TRACE_POP
return 0;
Expand Down
4 changes: 0 additions & 4 deletions SimulationRuntime/c/simulation/solver/dassl.h
View file
Expand Up @@ -62,10 +62,6 @@ typedef struct DASSL_DATA{
unsigned int* dasslStatistics;
unsigned int* dasslStatisticsTmp;

/* current context evaulation */
int currentContext;
int currentContextOld;

int* info;

int idid;
Expand Down
51 changes: 51 additions & 0 deletions SimulationRuntime/c/simulation/solver/model_help.c
View file
Expand Up @@ -912,6 +912,8 @@ void initializeDataStruc(DATA *data, threadData_t *threadData)
data->simulationInfo.mixedMethod = MIXED_SEARCH;
data->simulationInfo.newtonStrategy = NEWTON_PURE;
data->simulationInfo.nlsCsvInfomation = 0;
data->simulationInfo.currentContext = CONTEXT_ALGEBRAIC;
data->simulationInfo.jacobianEvals = data->modelData.nStates;

data->simulationInfo.zeroCrossings = (modelica_real*) calloc(data->modelData.nZeroCrossings, sizeof(modelica_real));
data->simulationInfo.zeroCrossingsPre = (modelica_real*) calloc(data->modelData.nZeroCrossings, sizeof(modelica_real));
Expand Down Expand Up @@ -1321,3 +1323,52 @@ modelica_real _event_div_real(modelica_real x1, modelica_real x2, modelica_integ
}

int measure_time_flag=0;

const char *context_string[CONTEXT_MAX] = {
"context UNKNOWN",
"context ODE evaluation",
"context algebraic evaluation",
"context event search",
"context jacobian evaluation",
};

/*! \fn setContext
*
* \param [ref] [data]
* \param [in] [currentTime]
* \param [in] [currentContext]
*
* Set current context in simulation info object
*/
void setContext(DATA* data, double* currentTime, int currentContext){
data->simulationInfo.currentContextOld = data->simulationInfo.currentContext;
data->simulationInfo.currentContext = currentContext;
infoStreamPrint(LOG_SOLVER, 0, "+++ Set context %s +++ at time %f", context_string[data->simulationInfo.currentContext], *currentTime);
if (currentContext == CONTEXT_JACOBIAN){
data->simulationInfo.currentJacobianEval = 0;
}
}

/*! \fn increaseJacContext
*
* \param [ref] [data]
*
* Increase Jacobian column context in simulation info object
*/
void increaseJacContext(DATA* data){
if (data->simulationInfo.currentContext == CONTEXT_JACOBIAN){
data->simulationInfo.currentJacobianEval++;
infoStreamPrint(LOG_SOLVER, 0, "+++ Increase Jacobian column context %s +++ to %d", context_string[data->simulationInfo.currentContext], data->simulationInfo.currentJacobianEval);
}
}

/*! \fn unsetContext
*
* \param [ref] [data]
*
* Restores previous context in simulation info object
*/
void unsetContext(DATA* data){
infoStreamPrint(LOG_SOLVER, 0, "--- Unset context %s ---", context_string[data->simulationInfo.currentContext]);
data->simulationInfo.currentContext = data->simulationInfo.currentContextOld;
}
4 changes: 4 additions & 0 deletions SimulationRuntime/c/simulation/solver/model_help.h
View file
Expand Up @@ -153,6 +153,10 @@ modelica_boolean GreaterEqZC(double a, double b, modelica_boolean);

extern int measure_time_flag;

void setContext(DATA* data, double* currentTime, int currentContext);
void increaseJacContext(DATA* data);
void unsetContext(DATA* data);

#ifdef __cplusplus
}
#endif
Expand Down
46 changes: 26 additions & 20 deletions SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c
View file
Expand Up @@ -782,7 +782,6 @@ static int getNumericalJacobianHomotopy(DATA_HOMOTOPY* solverData, double *x, do
const double delta_h = sqrt(DBL_EPSILON*2e1);
double delta_hh;
double xsave;

int i,j,l;

/* solverData->f1 must be set outside this function based on x */
Expand Down Expand Up @@ -1253,6 +1252,7 @@ static int newtonAlgorithm(DATA_HOMOTOPY* solverData, double* x)
vecDivScaling(solverData->n, solverData->f1, solverData->resScaling, solverData->fvecScaled);
debugVectorDouble(LOG_NLS_V,"function values:",solverData->f1, n);
debugVectorDouble(LOG_NLS_V,"scaled function values:",solverData->fvecScaled, n);

vecDivScaling(solverData->n, solverData->dy0, solverData->xScaling, solverData->dxScaled);
delta_x = vecNorm2(solverData->n, solverData->dy0);
delta_x_scaled = vecNorm2(solverData->n, solverData->dxScaled);
Expand Down Expand Up @@ -1338,7 +1338,7 @@ static int newtonAlgorithm(DATA_HOMOTOPY* solverData, double* x)
debugString(LOG_NLS_V, "NEWTON SOLVER DID CONVERGE TO A SOLUTION WITH LESS ACCURACY!!!");
printUnknowns(LOG_NLS_V, solverData);
debugString(LOG_NLS_V, "******************************************************");
solverData->error_f = error_f;
solverData->error_f = 0;

} else
{
Expand All @@ -1352,12 +1352,11 @@ static int newtonAlgorithm(DATA_HOMOTOPY* solverData, double* x)
solverData->numberOfIterations += numberOfIterations;
break;
}

assert = 1;
#ifndef OMC_EMCC
MMC_TRY_INTERNAL(simulationJumpBuffer)
#endif
/* calculate jacobian and function values (both stored in fJac, last column is fvec)*/
/* calculate jacobian and function values (both stored in fJac, last column is fvec) */
solverData->fJac_f(solverData, x, solverData->fJac);
assert = 0;
#ifndef OMC_EMCC
Expand Down Expand Up @@ -1763,30 +1762,37 @@ int solveHomotopy(DATA *data, threadData_t *threadData, int sysNumber)
#ifndef OMC_EMCC
MMC_TRY_INTERNAL(simulationJumpBuffer)
#endif
if (mixedSystem)
memcpy(relationsPreBackup, data->simulationInfo.relations, sizeof(modelica_boolean)*data->modelData.nRelations);

solverData->f(solverData, solverData->x0, solverData->f1);
/* Try to get out of here!!! */
error_f = vecNorm2(solverData->n, solverData->f1);
if ((error_f - solverData->error_f)<=0)
{
success = 1;
/* debug information */
debugString(LOG_NLS_V, "NO ITERATION NECESSARY!!!");
debugString(LOG_NLS_V, "******************************************************");
debugString(LOG_NLS_V,"SYSTEM SOLVED");
debugInt(LOG_NLS_V, "number of function calls: ",solverData->numberOfFunctionEvaluations-numberOfFunctionEvaluationsOld);
debugString(LOG_NLS_V, "------------------------------------------------------");
/* take the solution */
vecCopy(solverData->n, solverData->x0, systemData->nlsx);
debugVectorDouble(LOG_NLS_V,"Solution", solverData->x0, solverData->n);
/* reset continous flag */
((DATA*)data)->simulationInfo.solveContinuous = 0;
//infoStreamPrint(LOG_STDOUT, 0, "No Iteration at time %g needed new f = %g and old f1 = %g", solverData->timeValue, error_f, solverData->error_f);
if (mixedSystem && data->simulationInfo.discreteCall && isNotEqualVectorInt(((DATA*)data)->modelData.nRelations, ((DATA*)data)->simulationInfo.relations, relationsPreBackup)){}
else
{
success = 1;

free(relationsPreBackup);
debugString(LOG_NLS_V, "NO ITERATION NECESSARY!!!");
debugString(LOG_NLS_V, "******************************************************");
debugString(LOG_NLS_V,"SYSTEM SOLVED");
debugInt(LOG_NLS_V, "number of function calls: ",solverData->numberOfFunctionEvaluations-numberOfFunctionEvaluationsOld);
debugString(LOG_NLS_V, "------------------------------------------------------");

/* write statistics */
systemData->numberOfFEval = solverData->numberOfFunctionEvaluations;
vecCopy(solverData->n, solverData->x0, systemData->nlsx);
debugVectorDouble(LOG_NLS_V,"Solution", solverData->x0, solverData->n);

return success;
((DATA*)data)->simulationInfo.solveContinuous = 0;

free(relationsPreBackup);

systemData->numberOfFEval = solverData->numberOfFunctionEvaluations;

return success;
}
}
solverData->fJac_f(solverData, solverData->x0, solverData->fJac);
vecCopy(solverData->n, solverData->f1, solverData->fJac + solverData->n*solverData->n);
Expand Down

0 comments on commit 57d05fd

Please sign in to comment.