Skip to content

Commit

Permalink
Make non-default nonlinear solvers ...
Browse files Browse the repository at this point in the history 
... compatible with adaptive homotopy.

kinsol still not compatible

Belonging to [master]:
  - OpenModelica/OMCompiler#2179
  - OpenModelica/OpenModelica-testsuite#846
  • Loading branch information
@ptaeuber @OpenModelica-Hudson
ptaeuber authored and OpenModelica-Hudson committed Feb 8, 2018
1 parent 10a2c8b commit f5acc12
Show file tree
Hide file tree
Showing 6 changed files with 58 additions and 20 deletions.
1 change: 1 addition & 0 deletions SimulationRuntime/c/simulation/solver/initialization/initialization.c
View file
Expand Up @@ -177,6 +177,7 @@ void dumpInitialSolution(DATA *simData)
* \param [ref] [threadData]
*
* \author lochel
* \author ptaeuber
*/
static int symbolic_initialization(DATA *data, threadData_t *threadData)
{
Expand Down
10 changes: 5 additions & 5 deletions SimulationRuntime/c/simulation/solver/newtonIteration.c
View file
Expand Up @@ -93,19 +93,19 @@ int allocateNewtonData(int size, void** voiddata)
data->fvecScaled = (double*) malloc(size*sizeof(double));

data->n = size;
data->x = (double*) malloc(size*sizeof(double));
data->x = (double*) malloc((size+1)*sizeof(double));
data->fvec = (double*) calloc(size,sizeof(double));
data->xtol = 1e-6;
data->ftol = 1e-6;
data->maxfev = size*100;
data->epsfcn = DBL_EPSILON;
data->fjac = (double*) malloc((size*size)*sizeof(double));
data->fjac = (double*) malloc((size*(size+1))*sizeof(double));

data->rwork = (double*) malloc((size)*sizeof(double));
data->iwork = (int*) malloc(size*sizeof(int));

/* damped newton */
data->x_new = (double*) malloc(size*sizeof(double));
data->x_new = (double*) malloc((size+1)*sizeof(double));
data->x_increment = (double*) malloc(size*sizeof(double));
data->f_old = (double*) calloc(size,sizeof(double));
data->fvec_minimum = (double*) calloc(size,sizeof(double));
Expand Down Expand Up @@ -166,7 +166,8 @@ int freeNewtonData(void **voiddata)
*/
int _omc_newton(int(*f)(int*, double*, double*, void*, int), DATA_NEWTON* solverData, void* userdata)
{

DATA_USER* uData = (DATA_USER*) userdata;
DATA* data = (DATA*) uData->data;
int i, j, k = 0, l = 0, nrsh = 1;
int *n = &(solverData->n);
double *x = solverData->x;
Expand All @@ -183,7 +184,6 @@ int _omc_newton(int(*f)(int*, double*, double*, void*, int), DATA_NEWTON* solver
double error_f = 1.0 + *eps, scaledError_f = 1.0 + *eps, delta_x = 1.0 + *eps, delta_f = 1.0 + *eps, delta_x_scaled = 1.0 + *eps, lambda = 1.0;
double current_fvec_enorm, enorm_new;


if(ACTIVE_STREAM(LOG_NLS_V))
{
infoStreamPrint(LOG_NLS_V, 1, "######### Start Newton maxfev: %d #########", (int)*maxfev);
Expand Down
16 changes: 10 additions & 6 deletions SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c
View file
Expand Up @@ -910,8 +910,6 @@ static int wrapper_fvec(DATA_HOMOTOPY* solverData, double* x, double* f)
void *dataAndThreadData[2] = {solverData->data, solverData->threadData};
int iflag = 0;

if ((solverData->data)->simulationInfo->nonlinearSystemData[solverData->sysNumber].homotopySupport && !solverData->initHomotopy && (solverData->data)->simulationInfo->nonlinearSystemData[solverData->sysNumber].size > solverData->n)
x[solverData->n] = 1.0;
/*TODO: change input to residualFunc from data to systemData */
(solverData->data)->simulationInfo->nonlinearSystemData[solverData->sysNumber].residualFunc(dataAndThreadData, x, f, &iflag);
solverData->numberOfFunctionEvaluations++;
Expand Down Expand Up @@ -2156,10 +2154,16 @@ int solveHomotopy(DATA *data, threadData_t *threadData, int sysNumber)
debugVectorDouble(LOG_NLS_V,"System extrapolation", solverData->xStart, solverData->n);
}
vecCopy(solverData->n, solverData->xStart, solverData->x0);
// Initialize lambda variable with 0
solverData->x0[solverData->n] = 0.0;
solverData->x[solverData->n] = 0.0;
solverData->x1[solverData->n] = 0.0;
// Initialize lambda variable
if ((solverData->data)->simulationInfo->nonlinearSystemData[solverData->sysNumber].homotopySupport && !solverData->initHomotopy && (solverData->data)->simulationInfo->nonlinearSystemData[solverData->sysNumber].size > solverData->n) {
solverData->x0[solverData->n] = 1.0;
solverData->x[solverData->n] = 1.0;
solverData->x1[solverData->n] = 1.0;
} else {
solverData->x0[solverData->n] = 0.0;
solverData->x[solverData->n] = 0.0;
solverData->x1[solverData->n] = 0.0;
}
/* Use actual working point for scaling */
for (i=0;i<solverData->n;i++){
solverData->xScaling[i] = fmax(systemData->nominal[i],fabs(solverData->x0[i]));
Expand Down
23 changes: 18 additions & 5 deletions SimulationRuntime/c/simulation/solver/nonlinearSolverHybrd.c
View file
Expand Up @@ -75,9 +75,9 @@ int allocateHybrdData(int size, void** voiddata)
data->xScalefactors = (double*) malloc(size*sizeof(double));

data->n = size;
data->x = (double*) malloc(size*sizeof(double));
data->xSave = (double*) malloc(size*sizeof(double));
data->xScaled = (double*) malloc(size*sizeof(double));
data->x = (double*) malloc((size+1)*sizeof(double));
data->xSave = (double*) malloc((size+1)*sizeof(double));
data->xScaled = (double*) malloc((size+1)*sizeof(double));
data->fvec = (double*) calloc(size, sizeof(double));
data->fvecSave = (double*) calloc(size, sizeof(double));
data->xtol = 1e-12;
Expand All @@ -93,8 +93,8 @@ int allocateHybrdData(int size, void** voiddata)
data->info = 0;
data->nfev = 0;
data->njev = 0;
data->fjac = (double*) calloc((size*size), sizeof(double));
data->fjacobian = (double*) calloc((size*size), sizeof(double));
data->fjac = (double*) calloc((size*(size+1)), sizeof(double));
data->fjacobian = (double*) calloc((size*(size+1)), sizeof(double));
data->ldfjac = size;
data->r__ = (double*) malloc(((size*(size+1))/2)*sizeof(double));
data->lr = (size*(size + 1)) / 2;
Expand Down Expand Up @@ -445,6 +445,19 @@ int solveHybrd(DATA *data, threadData_t *threadData, int sysNumber)
relationsPreBackup = (modelica_boolean*) malloc(data->modelData->nRelations*sizeof(modelica_boolean));

solverData->numberOfFunctionEvaluations = 0;

// Initialize lambda variable
if (data->simulationInfo->nonlinearSystemData[sysNumber].homotopySupport) {
solverData->x[solverData->n] = 1.0;
solverData->xSave[solverData->n] = 1.0;
solverData->xScaled[solverData->n] = 1.0;
}
else {
solverData->x[solverData->n] = 0.0;
solverData->xSave[solverData->n] = 0.0;
solverData->xScaled[solverData->n] = 0.0;
}

/* debug output */
if(ACTIVE_STREAM(LOG_NLS_V))
{
Expand Down
11 changes: 10 additions & 1 deletion SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.c
View file
Expand Up @@ -200,7 +200,6 @@ int solveNewton(DATA *data, threadData_t *threadData, int sysNumber)
int retries2 = 0;
int nonContinuousCase = 0;
modelica_boolean *relationsPreBackup = NULL;
// int casualTearingSet = systemData->strictTearingFunctionCall != NULL;
int casualTearingSet = data->simulationInfo->nonlinearSystemData[sysNumber].strictTearingFunctionCall != NULL;

DATA_USER* userdata = (DATA_USER*)malloc(sizeof(DATA_USER));
Expand All @@ -225,6 +224,16 @@ int solveNewton(DATA *data, threadData_t *threadData, int sysNumber)
/* try to calculate jacobian only once at the beginning of the iteration */
solverData->calculate_jacobian = 0;

// Initialize lambda variable
if (data->simulationInfo->nonlinearSystemData[sysNumber].homotopySupport) {
solverData->x[solverData->n] = 1.0;
solverData->x_new[solverData->n] = 1.0;
}
else {
solverData->x[solverData->n] = 0.0;
solverData->x_new[solverData->n] = 0.0;
}

/* debug output */
if(ACTIVE_STREAM(LOG_NLS_V))
{
Expand Down
17 changes: 14 additions & 3 deletions SimulationRuntime/c/simulation/solver/nonlinearSystem.c
View file
Expand Up @@ -448,9 +448,10 @@ int initializeNonlinearSystems(DATA *data, threadData_t *threadData)
case NLS_KINSOL:
solverData = (struct dataSolver*) malloc(sizeof(struct dataSolver));
if (nonlinsys[i].homotopySupport && (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3)) {
nlsKinsolAllocate(size-1, &nonlinsys[i], data->simulationInfo->nlsLinearSolver);
solverData->ordinaryData = nonlinsys[i].solverData;
allocateHomotopyData(size-1, &(solverData->initHomotopyData));
// nlsKinsolAllocate(size-1, &nonlinsys[i], data->simulationInfo->nlsLinearSolver);
// solverData->ordinaryData = nonlinsys[i].solverData;
// allocateHomotopyData(size-1, &(solverData->initHomotopyData));
throwStreamPrint(threadData, "kinsol solver not compatible with adaptive homotopy approaches");
} else {
nlsKinsolAllocate(size, &nonlinsys[i], data->simulationInfo->nlsLinearSolver);
solverData->ordinaryData = nonlinsys[i].solverData;
Expand Down Expand Up @@ -824,6 +825,14 @@ int solveNLS(DATA *data, threadData_t *threadData, int sysNumber)
return success;
}

/*! \fn solve system with homotopy solver
*
* \param [in] [data]
* \param [in] [threadData]
* \param [in] [sysNumber] index of corresponding non-linear system
*
* \author ptaeuber
*/
int solveWithInitHomotopy(DATA *data, threadData_t *threadData, int sysNumber)
{
int success = 0;
Expand Down Expand Up @@ -867,6 +876,8 @@ int solveWithInitHomotopy(DATA *data, threadData_t *threadData, int sysNumber)
* \param [in] [data]
* \param [in] [threadData]
* \param [in] [sysNumber] index of corresponding non-linear system
*
* \author ptaeuber
*/
int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
{
Expand Down

0 comments on commit f5acc12

Please sign in to comment.