Skip to content
This repository was archived by the owner on May 18, 2019. It is now read-only.

Commit f5acc12

Browse files
ptaeuberOpenModelica-Hudson
ptaeuber
authored and
OpenModelica-Hudson
committed
Make non-default nonlinear solvers ...
... compatible with adaptive homotopy. kinsol still not compatible Belonging to [master]: - #2179 - OpenModelica/OpenModelica-testsuite#846
1 parent 10a2c8b commit f5acc12

File tree

6 files changed

+58
-20
lines changed

6 files changed

+58
-20
lines changed

‎SimulationRuntime/c/simulation/solver/initialization/initialization.c‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -177,6 +177,7 @@ void dumpInitialSolution(DATA *simData)
177177
* \param [ref] [threadData]
178178
*
179179
* \author lochel
180+
* \author ptaeuber
180181
*/
181182
static int symbolic_initialization(DATA *data, threadData_t *threadData)
182183
{

‎SimulationRuntime/c/simulation/solver/newtonIteration.c‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -93,19 +93,19 @@ int allocateNewtonData(int size, void** voiddata)
9393
data->fvecScaled = (double*) malloc(size*sizeof(double));
9494

9595
data->n = size;
96-
data->x = (double*) malloc(size*sizeof(double));
96+
data->x = (double*) malloc((size+1)*sizeof(double));
9797
data->fvec = (double*) calloc(size,sizeof(double));
9898
data->xtol = 1e-6;
9999
data->ftol = 1e-6;
100100
data->maxfev = size*100;
101101
data->epsfcn = DBL_EPSILON;
102-
data->fjac = (double*) malloc((size*size)*sizeof(double));
102+
data->fjac = (double*) malloc((size*(size+1))*sizeof(double));
103103

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

107107
/* damped newton */
108-
data->x_new = (double*) malloc(size*sizeof(double));
108+
data->x_new = (double*) malloc((size+1)*sizeof(double));
109109
data->x_increment = (double*) malloc(size*sizeof(double));
110110
data->f_old = (double*) calloc(size,sizeof(double));
111111
data->fvec_minimum = (double*) calloc(size,sizeof(double));
@@ -166,7 +166,8 @@ int freeNewtonData(void **voiddata)
166166
*/
167167
int _omc_newton(int(*f)(int*, double*, double*, void*, int), DATA_NEWTON* solverData, void* userdata)
168168
{
169-
169+
DATA_USER* uData = (DATA_USER*) userdata;
170+
DATA* data = (DATA*) uData->data;
170171
int i, j, k = 0, l = 0, nrsh = 1;
171172
int *n = &(solverData->n);
172173
double *x = solverData->x;
@@ -183,7 +184,6 @@ int _omc_newton(int(*f)(int*, double*, double*, void*, int), DATA_NEWTON* solver
183184
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;
184185
double current_fvec_enorm, enorm_new;
185186

186-
187187
if(ACTIVE_STREAM(LOG_NLS_V))
188188
{
189189
infoStreamPrint(LOG_NLS_V, 1, "######### Start Newton maxfev: %d #########", (int)*maxfev);

‎SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c‎

Lines changed: 10 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -910,8 +910,6 @@ static int wrapper_fvec(DATA_HOMOTOPY* solverData, double* x, double* f)
910910
void *dataAndThreadData[2] = {solverData->data, solverData->threadData};
911911
int iflag = 0;
912912

913-
if ((solverData->data)->simulationInfo->nonlinearSystemData[solverData->sysNumber].homotopySupport && !solverData->initHomotopy && (solverData->data)->simulationInfo->nonlinearSystemData[solverData->sysNumber].size > solverData->n)
914-
x[solverData->n] = 1.0;
915913
/*TODO: change input to residualFunc from data to systemData */
916914
(solverData->data)->simulationInfo->nonlinearSystemData[solverData->sysNumber].residualFunc(dataAndThreadData, x, f, &iflag);
917915
solverData->numberOfFunctionEvaluations++;
@@ -2156,10 +2154,16 @@ int solveHomotopy(DATA *data, threadData_t *threadData, int sysNumber)
21562154
debugVectorDouble(LOG_NLS_V,"System extrapolation", solverData->xStart, solverData->n);
21572155
}
21582156
vecCopy(solverData->n, solverData->xStart, solverData->x0);
2159-
// Initialize lambda variable with 0
2160-
solverData->x0[solverData->n] = 0.0;
2161-
solverData->x[solverData->n] = 0.0;
2162-
solverData->x1[solverData->n] = 0.0;
2157+
// Initialize lambda variable
2158+
if ((solverData->data)->simulationInfo->nonlinearSystemData[solverData->sysNumber].homotopySupport && !solverData->initHomotopy && (solverData->data)->simulationInfo->nonlinearSystemData[solverData->sysNumber].size > solverData->n) {
2159+
solverData->x0[solverData->n] = 1.0;
2160+
solverData->x[solverData->n] = 1.0;
2161+
solverData->x1[solverData->n] = 1.0;
2162+
} else {
2163+
solverData->x0[solverData->n] = 0.0;
2164+
solverData->x[solverData->n] = 0.0;
2165+
solverData->x1[solverData->n] = 0.0;
2166+
}
21632167
/* Use actual working point for scaling */
21642168
for (i=0;i<solverData->n;i++){
21652169
solverData->xScaling[i] = fmax(systemData->nominal[i],fabs(solverData->x0[i]));

‎SimulationRuntime/c/simulation/solver/nonlinearSolverHybrd.c‎

Lines changed: 18 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -75,9 +75,9 @@ int allocateHybrdData(int size, void** voiddata)
7575
data->xScalefactors = (double*) malloc(size*sizeof(double));
7676

7777
data->n = size;
78-
data->x = (double*) malloc(size*sizeof(double));
79-
data->xSave = (double*) malloc(size*sizeof(double));
80-
data->xScaled = (double*) malloc(size*sizeof(double));
78+
data->x = (double*) malloc((size+1)*sizeof(double));
79+
data->xSave = (double*) malloc((size+1)*sizeof(double));
80+
data->xScaled = (double*) malloc((size+1)*sizeof(double));
8181
data->fvec = (double*) calloc(size, sizeof(double));
8282
data->fvecSave = (double*) calloc(size, sizeof(double));
8383
data->xtol = 1e-12;
@@ -93,8 +93,8 @@ int allocateHybrdData(int size, void** voiddata)
9393
data->info = 0;
9494
data->nfev = 0;
9595
data->njev = 0;
96-
data->fjac = (double*) calloc((size*size), sizeof(double));
97-
data->fjacobian = (double*) calloc((size*size), sizeof(double));
96+
data->fjac = (double*) calloc((size*(size+1)), sizeof(double));
97+
data->fjacobian = (double*) calloc((size*(size+1)), sizeof(double));
9898
data->ldfjac = size;
9999
data->r__ = (double*) malloc(((size*(size+1))/2)*sizeof(double));
100100
data->lr = (size*(size + 1)) / 2;
@@ -445,6 +445,19 @@ int solveHybrd(DATA *data, threadData_t *threadData, int sysNumber)
445445
relationsPreBackup = (modelica_boolean*) malloc(data->modelData->nRelations*sizeof(modelica_boolean));
446446

447447
solverData->numberOfFunctionEvaluations = 0;
448+
449+
// Initialize lambda variable
450+
if (data->simulationInfo->nonlinearSystemData[sysNumber].homotopySupport) {
451+
solverData->x[solverData->n] = 1.0;
452+
solverData->xSave[solverData->n] = 1.0;
453+
solverData->xScaled[solverData->n] = 1.0;
454+
}
455+
else {
456+
solverData->x[solverData->n] = 0.0;
457+
solverData->xSave[solverData->n] = 0.0;
458+
solverData->xScaled[solverData->n] = 0.0;
459+
}
460+
448461
/* debug output */
449462
if(ACTIVE_STREAM(LOG_NLS_V))
450463
{

‎SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.c‎

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -200,7 +200,6 @@ int solveNewton(DATA *data, threadData_t *threadData, int sysNumber)
200200
int retries2 = 0;
201201
int nonContinuousCase = 0;
202202
modelica_boolean *relationsPreBackup = NULL;
203-
// int casualTearingSet = systemData->strictTearingFunctionCall != NULL;
204203
int casualTearingSet = data->simulationInfo->nonlinearSystemData[sysNumber].strictTearingFunctionCall != NULL;
205204

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

227+
// Initialize lambda variable
228+
if (data->simulationInfo->nonlinearSystemData[sysNumber].homotopySupport) {
229+
solverData->x[solverData->n] = 1.0;
230+
solverData->x_new[solverData->n] = 1.0;
231+
}
232+
else {
233+
solverData->x[solverData->n] = 0.0;
234+
solverData->x_new[solverData->n] = 0.0;
235+
}
236+
228237
/* debug output */
229238
if(ACTIVE_STREAM(LOG_NLS_V))
230239
{

‎SimulationRuntime/c/simulation/solver/nonlinearSystem.c‎

Lines changed: 14 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -448,9 +448,10 @@ int initializeNonlinearSystems(DATA *data, threadData_t *threadData)
448448
case NLS_KINSOL:
449449
solverData = (struct dataSolver*) malloc(sizeof(struct dataSolver));
450450
if (nonlinsys[i].homotopySupport && (data->callback->useHomotopy == 2 || data->callback->useHomotopy == 3)) {
451-
nlsKinsolAllocate(size-1, &nonlinsys[i], data->simulationInfo->nlsLinearSolver);
452-
solverData->ordinaryData = nonlinsys[i].solverData;
453-
allocateHomotopyData(size-1, &(solverData->initHomotopyData));
451+
// nlsKinsolAllocate(size-1, &nonlinsys[i], data->simulationInfo->nlsLinearSolver);
452+
// solverData->ordinaryData = nonlinsys[i].solverData;
453+
// allocateHomotopyData(size-1, &(solverData->initHomotopyData));
454+
throwStreamPrint(threadData, "kinsol solver not compatible with adaptive homotopy approaches");
454455
} else {
455456
nlsKinsolAllocate(size, &nonlinsys[i], data->simulationInfo->nlsLinearSolver);
456457
solverData->ordinaryData = nonlinsys[i].solverData;
@@ -824,6 +825,14 @@ int solveNLS(DATA *data, threadData_t *threadData, int sysNumber)
824825
return success;
825826
}
826827

828+
/*! \fn solve system with homotopy solver
829+
*
830+
* \param [in] [data]
831+
* \param [in] [threadData]
832+
* \param [in] [sysNumber] index of corresponding non-linear system
833+
*
834+
* \author ptaeuber
835+
*/
827836
int solveWithInitHomotopy(DATA *data, threadData_t *threadData, int sysNumber)
828837
{
829838
int success = 0;
@@ -867,6 +876,8 @@ int solveWithInitHomotopy(DATA *data, threadData_t *threadData, int sysNumber)
867876
* \param [in] [data]
868877
* \param [in] [threadData]
869878
* \param [in] [sysNumber] index of corresponding non-linear system
879+
*
880+
* \author ptaeuber
870881
*/
871882
int solve_nonlinear_system(DATA *data, threadData_t *threadData, int sysNumber)
872883
{

0 commit comments

Comments
 (0)