Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 0.0% 0 / 0 / 116
Functions: 0.0% 0 / 0 / 2
Branches: 0.0% 0 / 0 / 66
OMCompiler/SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.c
Line Branch Exec Source
1 /*
2 * This file belongs to the OpenModelica Run-Time System
3 *
4 * Copyright (c) 1998-2026, Open Source Modelica Consortium (OSMC), c/o Linköpings
5 * universitet, Department of Computer and Information Science, SE-58183 Linköping, Sweden. All rights
6 * reserved.
7 *
8 * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE
9 * AGPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8. ANY
10 * USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES RECIPIENT'S
11 * ACCEPTANCE OF THE BSD NEW LICENSE OR THE OSMC PUBLIC LICENSE OR THE AGPL
12 * VERSION 3, ACCORDING TO RECIPIENTS CHOICE.
13 *
14 * The OpenModelica software and the OSMC (Open Source Modelica Consortium) Public License
15 * (OSMC-PL) are obtained from OSMC, either from the above address, from the URLs:
16 * http://www.openmodelica.org or https://github.com/OpenModelica/ or
17 * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica distribution. GNU
18 * AGPL version 3 is obtained from: https://www.gnu.org/licenses/licenses.html#GPL. The BSD NEW
19 * License is obtained from: http://www.opensource.org/licenses/BSD-3-Clause.
20 *
21 * This program is distributed WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY
23 * SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF
24 * OSMC-PL.
25 *
26 */
27
28 /*! \file nonlinearSolverNewton.c
29 */
30
31 #ifdef __cplusplus
32 extern "C" {
33 #endif
34
35 #include <math.h>
36 #include <stdlib.h>
37 #include <string.h> /* memcpy */
38
39 #include "../simulation_info_json.h"
40 #include "../jacobian_util.h"
41 #include "util/omc_error.h"
42
43 #include "util/varinfo.h"
44 #include "model_help.h"
45
46 #include "nonlinearSystem.h"
47 #include "nonlinearSolverNewton.h"
48 #include "newtonIteration.h"
49
50 #include "external_input.h"
51
52 /* Private function prototypes */
53
54 int wrapper_fvec_newton(int n, double* x, double* fvec, NLS_USERDATA* userData, int fj);
55
56 /* External function prototypes */
57
58 extern double enorm_(int *n, double *x);
59 extern int dgesv_(int *n, int *nrhs, doublereal *a, int *lda, int *ipiv, doublereal *b, int *ldb, int *info);
60
61
62 /**
63 * @brief Calculate residual f(x) or Jacobian J(x).
64 *
65 * @param n Size of vector x.
66 * @param x Input vector x.
67 * Also used as work array, but will be reverted before function exits.
68 * @param fvec Value of f(x).
69 * Will be computed if fj = 1.
70 * Will be used to compute Jacobian if fj = 0.
71 * @param userData Pointer to Newton user data.
72 * @param fj Decides whether the function values or the jacobian matrix shall be calculated.
73 * fj = 1: calculate function values
74 * fj = 0: calculate jacobian matrix
75 * @return int Returns 1 on success (probably)
76 */
77 ✗ int wrapper_fvec_newton(int n, double* x, double* fvec, NLS_USERDATA* userData, int fj)
78 {
79 ✗ DATA* data = userData->data;
80 ✗ threadData_t *threadData = userData->threadData;
81 int sysNumber = userData->sysNumber;
82 ✗ NONLINEAR_SYSTEM_DATA* nlsData = userData->nlsData;
83 ✗ JACOBIAN* jacobian = userData->analyticJacobian;
84
85 ✗ DATA_NEWTON* solverData = (DATA_NEWTON*)(nlsData->solverData);
86 ✗ RESIDUAL_USERDATA resUserData = {.data=data, .threadData=threadData, .solverData=userData->solverData};
87 ✗ int flag = 1;
88
89 ✗ if (fj) {
90 ✗ nlsData->residualFunc(&resUserData, x, fvec, &flag);
91 } else {
92 /* performance measurement */
93 ✗ rt_ext_tp_tick(&nlsData->jacobianTimeClock);
94
95 ✗ if(nlsData->jacobianIndex != -1 && jacobian != NULL ) {
96 /* call generic dense Jacobian */
97 ✗ evalJacobian(data, threadData, jacobian, NULL, solverData->fjac, TRUE);
98 } else {
99 ✗ double delta_h = sqrt(solverData->epsfcn);
100 double delta_hh;
101 double xsave;
102
103 int i,j,l, linear=0;
104
105 ✗ for(i = 0; i < n; i++) {
106 ✗ delta_hh = fmax(delta_h * fmax(fabs(x[i]), fabs(fvec[i])), delta_h);
107 ✗ delta_hh = ((fvec[i] >= 0) ? delta_hh : -delta_hh);
108 ✗ delta_hh = x[i] + delta_hh - x[i];
109 xsave = x[i];
110 ✗ x[i] += delta_hh;
111 ✗ delta_hh = 1. / delta_hh;
112
113 ✗ wrapper_fvec_newton(n, x, solverData->rwork, userData, 1);
114 ✗ solverData->nfev++;
115
116 ✗ for(j = 0; j < n; j++) {
117 ✗ l = i * n + j;
118 ✗ solverData->fjac[l] = (solverData->rwork[j] - fvec[j]) * delta_hh;
119 }
120 ✗ x[i] = xsave;
121 }
122 }
123 /* performance measurement and statistics */
124 ✗ nlsData->jacobianTime += rt_ext_tp_tock(&(nlsData->jacobianTimeClock));
125 ✗ nlsData->numberOfJEval++;
126 }
127 ✗ return flag;
128 }
129
130 /**
131 * @brief Solve non-linear system with Newton method.
132 *
133 * @param data Runtime data struct.
134 * @param threadData Thread data for error handling.
135 * @param nlsData Pointer to non-linear system data.
136 * @return NLS_SOLVER_STATUS Return NLS_SOLVED on success and NLS_FAILED otherwise.
137 */
138 ✗ NLS_SOLVER_STATUS solveNewton(DATA *data, threadData_t *threadData, NONLINEAR_SYSTEM_DATA* nlsData)
139 {
140 ✗ DATA_NEWTON* solverData = (DATA_NEWTON*)(nlsData->solverData);
141
142 int eqSystemNumber = 0;
143 int i;
144 double xerror = -1, xerror_scaled = -1;
145 NLS_SOLVER_STATUS success = NLS_FAILED;
146 int nfunc_evals = 0;
147 ✗ double local_tol = solverData->ftol;
148
149 int giveUp = 0;
150 int retries = 0;
151 int retries2 = 0;
152 int nonContinuousCase = 0;
153 modelica_boolean *relationsPreBackup = NULL;
154 ✗ int casualTearingSet = nlsData->strictTearingFunctionCall != NULL;
155
156 /*
157 * We are given the number of the non-linear system.
158 * We want to look it up among all equations.
159 */
160 ✗ eqSystemNumber = nlsData->equationIndex;
161
162 ✗ relationsPreBackup = (modelica_boolean*) malloc(data->modelData->nRelations*sizeof(modelica_boolean));
163
164 ✗ solverData->nfev = 0;
165
166 /* try to calculate jacobian only once at the beginning of the iteration */
167 ✗ solverData->calculate_jacobian = 0;
168
169 // Initialize lambda variable
170 ✗ if (nlsData->homotopySupport) {
171 ✗ solverData->x[solverData->n] = 1.0;
172 ✗ solverData->x_new[solverData->n] = 1.0;
173 }
174 else {
175 ✗ solverData->x[solverData->n] = 0.0;
176 ✗ solverData->x_new[solverData->n] = 0.0;
177 }
178
179 /* debug output */
180 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
181 {
182 ✗ int indexes[2] = {1, eqSystemNumber};
183 ✗ infoStreamPrintWithEquationIndexes(OMC_LOG_NLS_V, omc_dummyFileInfo, 1, indexes,
184 "Start solving Non-Linear System %d (size %d) at time %g with Newton Solver",
185 ✗ eqSystemNumber, (int) nlsData->size, data->localData[0]->timeValue);
186
187 ✗ for(i = 0; i < solverData->n; i++) {
188 ✗ infoStreamPrint(OMC_LOG_NLS_V, 1, "x[%d] = %.15e", i, data->simulationInfo->discreteCall ? nlsData->nlsx[i] : nlsData->nlsxExtrapolation[i]);
189 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "nominal = %g +++ nlsx = %g +++ old = %g +++ extrapolated = %g",
190 ✗ nlsData->nominal[i], nlsData->nlsx[i], nlsData->nlsxOld[i], nlsData->nlsxExtrapolation[i]);
191 ✗ messageClose(OMC_LOG_NLS_V);
192 }
193 ✗ messageClose(OMC_LOG_NLS_V);
194 }
195
196 /* set x vector */
197 ✗ if(data->simulationInfo->discreteCall) {
198 ✗ memcpy(solverData->x, nlsData->nlsx, solverData->n*(sizeof(double)));
199 } else {
200 ✗ memcpy(solverData->x, nlsData->nlsxExtrapolation, solverData->n*(sizeof(double)));
201 }
202 ✗ solverData->time = data->localData[0]->timeValue;
203 ✗ solverData->initial = data->simulationInfo->initial;
204
205 /* start solving loop */
206 ✗ while(!giveUp && success != NLS_SOLVED)
207 {
208
209 giveUp = 1;
210 ✗ solverData->newtonStrategy = data->simulationInfo->newtonStrategy;
211 ✗ _omc_newton((genericResidualFunc*)wrapper_fvec_newton, solverData, solverData->userData);
212
213 /* check for proper inputs */
214 ✗ if(solverData->info == 0)
215 ✗ printErrorEqSyst(IMPROPER_INPUT, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber), data->localData[0]->timeValue);
216
217 /* reset non-contunuousCase */
218 ✗ if(nonContinuousCase && xerror > local_tol && xerror_scaled > local_tol)
219 {
220 ✗ memcpy(data->simulationInfo->relationsPre, relationsPreBackup, sizeof(modelica_boolean)*data->modelData->nRelations);
221 nonContinuousCase = 0;
222 }
223
224 /* check for error */
225 ✗ xerror_scaled = enorm_(&solverData->n, solverData->fvecScaled);
226 ✗ xerror = enorm_(&solverData->n, solverData->fvec);
227
228 /* solution found */
229 ✗ if((xerror <= local_tol || xerror_scaled <= local_tol) && solverData->info > 0)
230 {
231 success = NLS_SOLVED;
232 ✗ nfunc_evals += solverData->nfev;
233 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
234 {
235 ✗ infoStreamPrint(OMC_LOG_NLS_V, 1, "System solved");
236 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "%d restarts", retries);
237 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "nfunc = %d +++ error = %.15e +++ error_scaled = %.15e", nfunc_evals, xerror, xerror_scaled);
238 ✗ for(i = 0; i < solverData->n; i++)
239 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "x[%d] = %.15e\n\tresidual = %e", i, solverData->x[i], solverData->fvec[i]);
240 ✗ messageClose(OMC_LOG_NLS_V);
241 }
242
243 /* take the solution */
244 ✗ memcpy(nlsData->nlsx, solverData->x, solverData->n*(sizeof(double)));
245
246 /* Then try with old values (instead of extrapolating )*/
247 }
248 // If this is the casual tearing set (only exists for dynamic tearing), break after first try
249 ✗ else if(retries < 1 && casualTearingSet)
250 {
251 giveUp = 1;
252 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "### No Solution for the casual tearing set at the first try! ###");
253 }
254 ✗ else if(retries < 1)
255 {
256 ✗ memcpy(solverData->x, nlsData->nlsxOld, solverData->n*(sizeof(double)));
257
258 ✗ retries++;
259 giveUp = 0;
260 ✗ nfunc_evals += solverData->nfev;
261 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, " - iteration making no progress:\t try old values.");
262 /* try to vary the initial values */
263
264 /* evaluate jacobian in every step now */
265 ✗ solverData->calculate_jacobian = 1;
266 }
267 ✗ else if(retries < 2)
268 {
269 ✗ for(i = 0; i < solverData->n; i++)
270 ✗ solverData->x[i] += nlsData->nominal[i] * 0.01;
271 retries++;
272 giveUp = 0;
273 ✗ nfunc_evals += solverData->nfev;
274 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, " - iteration making no progress:\t vary solution point by 1%%.");
275 /* try to vary the initial values */
276 }
277 ✗ else if(retries < 3)
278 {
279 ✗ for(i = 0; i < solverData->n; i++)
280 ✗ solverData->x[i] = nlsData->nominal[i];
281 retries++;
282 giveUp = 0;
283 ✗ nfunc_evals += solverData->nfev;
284 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, " - iteration making no progress:\t try nominal values as initial solution.");
285 }
286 ✗ else if(retries < 4 && data->simulationInfo->discreteCall)
287 {
288 /* try to solve non-continuous
289 * work-a-round: since other wise some model does
290 * stuck in event iteration. e.g.: Modelica.Mechanics.Rotational.Examples.HeatLosses
291 */
292
293 ✗ memcpy(solverData->x, nlsData->nlsxOld, solverData->n*(sizeof(double)));
294 retries++;
295
296 /* try to solve a discontinuous system */
297 nonContinuousCase = 1;
298 ✗ memcpy(relationsPreBackup, data->simulationInfo->relationsPre, sizeof(modelica_boolean)*data->modelData->nRelations);
299
300 giveUp = 0;
301 ✗ nfunc_evals += solverData->nfev;
302 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, " - iteration making no progress:\t try to solve a discontinuous system.");
303 }
304 ✗ else if(retries2 < 4)
305 {
306 ✗ memcpy(solverData->x, nlsData->nlsxOld, solverData->n*(sizeof(double)));
307 /* reduce tolarance */
308 ✗ local_tol = local_tol*10;
309
310 retries = 0;
311 ✗ retries2++;
312 giveUp = 0;
313 ✗ nfunc_evals += solverData->nfev;
314 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, " - iteration making no progress:\t reduce the tolerance slightly to %e.", local_tol);
315 }
316 else
317 {
318 ✗ printErrorEqSyst(ERROR_AT_TIME, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber), data->localData[0]->timeValue);
319 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
320 {
321 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "### No Solution! ###\n after %d restarts", retries);
322 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "nfunc = %d +++ error = %.15e +++ error_scaled = %.15e", nfunc_evals, xerror, xerror_scaled);
323 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
324 ✗ for(i = 0; i < solverData->n; i++)
325 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "x[%d] = %.15e\n\tresidual = %e", i, solverData->x[i], solverData->fvec[i]);
326 }
327 }
328 }
329
330 ✗ free(relationsPreBackup);
331
332 /* write statistics */
333 ✗ nlsData->numberOfFEval = solverData->numberOfFunctionEvaluations;
334 ✗ nlsData->numberOfIterations = solverData->numberOfIterations;
335
336 ✗ return success;
337 }
338