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 / 400
Functions: 0.0% 0 / 0 / 18
Branches: 0.0% 0 / 0 / 475
OMCompiler/SimulationRuntime/cpp/Solver/Kinsol/Kinsol.cpp
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 /** @addtogroup solverKinsol
29 *
30 * @{
31 */
32
33 #include <Core/ModelicaDefine.h> // This has to be first to include
34 #include <Core/Modelica.h> // But... why?
35
36 #if defined(__vxworks)
37 #include<wvLib.h>
38 #endif
39
40 #include <kinsol/kinsol.h> // Main header file for KINSOL
41 #include <nvector/nvector_serial.h> // Default serial vectors
42 #include <sunlinsol/sunlinsol_dense.h> // Default dense linear solver
43 #include <sunlinsol/sunlinsol_spgmr.h> // Scaled, Preconditioned, Generalized Minimum Residual iterative linear solver
44 #include <sunlinsol/sunlinsol_spbcgs.h> // Scaled, Preconditioned, Bi-Conjugate Gradient, Stabilized iterative linear solver
45
46 #include <Solver/Kinsol/FactoryExport.h>
47 #include <Solver/Kinsol/Kinsol.h>
48 #include <Solver/Kinsol/KinsolSettings.h>
49
50 #include <Core/Math/ILapack.h>
51 #include <Core/Utils/extension/logger.hpp>
52
53 /**\Callback function for Kinsol to calculate right hand side, calls internal Kinsol member function
54 * \param [in] y variables vector
55 * \param [in] fval right hand side vecotre
56 * \param [in] user_data user data pointer is used to access Kinsol instance
57 * \return status value
58 */
59 ✗ int kin_fCallback(N_Vector y,N_Vector fval, void *user_data) {
60 Kinsol* myKinsol = (Kinsol*)(user_data);
61 ✗ return myKinsol->kin_f(y,fval,user_data);
62 }
63
64 /**
65 * @brief Construct a new Kinsol:: Kinsol object
66 *
67 * @param settings General parameters for a non-linear solver object
68 * @param algLoop Algebraic loop
69 */
70 ✗ Kinsol::Kinsol(INonLinSolverSettings* settings, shared_ptr<INonLinearAlgLoop> algLoop)
71 :AlgLoopSolverDefaultImplementation()
72 ,_algLoop (algLoop)
73 ✗ , _kinsolSettings ((INonLinSolverSettings*)settings)
74 ✗ , _y (NULL)
75 ✗ , _y0 (NULL)
76 ✗ , _yScale (NULL)
77 ✗ , _fScale (NULL)
78 ✗ , _f (NULL)
79 ✗ , _helpArray (NULL)
80 ✗ , _currentIterate (NULL)
81 ✗ , _jac (NULL)
82 ✗ , _fHelp (NULL)
83 ✗ , _yHelp (NULL)
84 ✗ , _fnorm (10.0)
85 ✗ , _currentIterateNorm (100.0)
86 ✗ , _firstCall (true)
87 ✗ , _usedCompletePivoting (false)
88 ✗ , _usedIterativeSolver (false)
89 ✗ , _iterationStatus (CONTINUE)
90 ✗ , _Kin_y (NULL)
91 ✗ , _Kin_y0 (NULL)
92 ✗ , _Kin_yScale (NULL)
93 ✗ , _Kin_fScale (NULL)
94 ✗ , _Kin_ySolver (NULL)
95 ✗ , _Kin_linSol (NULL)
96 ✗ , _Kin_J (NULL)
97 ✗ , _kinMem (NULL)
98 ✗ , _sunctx (NULL)
99 ✗ , _fValid (false)
100 ✗ , _y_old (NULL)
101 ✗ , _y_new (NULL)
102 ✗ , _solverErrorNotificationGiven(false)
103 {
104 ✗ _max_dimSys = 100;
105 ✗ _max_dimZeroFunc=50;
106 ✗ _data = ((void*)this);
107 ✗ if (_algLoop) {
108 ✗ _single_instance = false;
109 ✗ AlgLoopSolverDefaultImplementation::initialize(_algLoop->getDimZeroFunc(),_algLoop->getDimReal());
110 } else {
111 ✗ _single_instance = true;
112 ✗ AlgLoopSolverDefaultImplementation::initialize(_max_dimZeroFunc,_max_dimSys);
113 }
114 ✗ }
115
116 /**
117 * @brief Destroy the Kinsol:: Kinsol object
118 *
119 */
120 ✗ Kinsol::~Kinsol() {
121 ✗ if(_y) delete [] _y;
122 ✗ if(_y0) delete [] _y0;
123 ✗ if(_y_old) delete [] _y_old;
124 ✗ if(_y_new) delete [] _y_new;
125 ✗ if(_yScale) delete [] _yScale;
126 ✗ if(_fScale) delete [] _fScale;
127 ✗ if(_f) delete [] _f;
128 ✗ if(_helpArray) delete [] _helpArray;
129 ✗ if(_jac) delete [] _jac;
130 ✗ if(_fHelp) delete [] _fHelp;
131 ✗ if(_currentIterate) delete [] _currentIterate;
132 ✗ if(_yHelp) delete [] _yHelp;
133
134 ✗ N_VDestroy_Serial(_Kin_y);
135 ✗ N_VDestroy_Serial(_Kin_y0);
136 ✗ N_VDestroy_Serial(_Kin_yScale);
137 ✗ N_VDestroy_Serial(_Kin_fScale);
138 ✗ N_VDestroy_Serial(_Kin_ySolver);
139 ✗ SUNMatDestroy(_Kin_J);
140 ✗ SUNLinSolFree(_Kin_linSol);
141 ✗ KINFree(&_kinMem);
142 ✗ SUNContext_Free(&_sunctx);
143 ✗ }
144
145 /**
146 * @brief Initialize Kinsol solver
147 *
148 * If the solver was already initialized reinitialize the solver.
149 *
150 */
151 ✗ void Kinsol::initialize() {
152 int idid;
153 ✗ if(!_algLoop) {
154 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "algloop system is not initialized");
155 }
156 ✗ if(_firstCall) {
157 ✗ _algLoop->initialize();
158 }
159
160 ✗ _firstCall = false;
161 ✗ _sparse = _algLoop->getUseSparseFormat();
162 ✗ _dimSys =_algLoop->getDimReal();
163
164 // Free data, if it's not NULL allready
165 ✗ if(_y) delete [] _y;
166 ✗ if(_y0) delete [] _y0;
167 ✗ if(_yScale) delete [] _yScale;
168 ✗ if(_fScale) delete [] _fScale;
169 ✗ if(_f) delete [] _f;
170 ✗ if(_helpArray) delete [] _helpArray;
171 ✗ if(_jac) delete [] _jac;
172 ✗ if(_yHelp) delete [] _yHelp;
173 ✗ if(_fHelp) delete [] _fHelp;
174 ✗ if(_currentIterate) delete [] _currentIterate;
175 ✗ if(_y_old) delete [] _y_old;
176 ✗ if(_y_new) delete [] _y_new;
177 ✗ N_VDestroy_Serial(_Kin_y);
178 ✗ N_VDestroy_Serial(_Kin_y0);
179 ✗ N_VDestroy_Serial(_Kin_yScale);
180 ✗ N_VDestroy_Serial(_Kin_fScale);
181 ✗ N_VDestroy(_Kin_ySolver);
182 ✗ SUNMatDestroy(_Kin_J);
183 ✗ SUNLinSolFree(_Kin_linSol);
184 ✗ KINFree(&_kinMem);
185 ✗ SUNContext_Free(&_sunctx);
186
187 // Initialize vectors
188 ✗ _y = new double[_dimSys];
189 ✗ _y0 = new double[_dimSys];
190 ✗ _yScale = new double[_dimSys];
191 ✗ _fScale = new double[_dimSys];
192 ✗ _f = new double[_dimSys];
193 ✗ _helpArray = new double[_dimSys];
194 ✗ _currentIterate = new double[_dimSys];
195 ✗ _y_old = new double[_dimSys];
196 ✗ _y_new = new double[_dimSys];
197 ✗ _jac = new double[_dimSys*_dimSys];
198 ✗ _yHelp = new double[_dimSys];
199 ✗ _fHelp = new double[_dimSys];
200
201 ✗ _algLoop->getReal(_y);
202 ✗ _algLoop->getReal(_y0);
203 ✗ _algLoop->getReal(_y_new);
204 ✗ _algLoop->getReal(_y_old);
205 ✗ memset(_f, 0, _dimSys*sizeof(double));
206 ✗ memset(_helpArray, 0, _dimSys*sizeof(double));
207 ✗ memset(_yHelp, 0, _dimSys*sizeof(double));
208 ✗ memset(_fHelp, 0, _dimSys*sizeof(double));
209 ✗ memset(_jac, 0, _dimSys*_dimSys*sizeof(double));
210 ✗ memset(_currentIterate, 0, _dimSys*sizeof(double));
211
212 // Scale y
213 ✗ _algLoop->getNominalReal(_yScale);
214 ✗ for (int i=0; i<_dimSys; i++) {
215 ✗ if(_yScale[i] != 0) {
216 ✗ _yScale[i] = 1/_yScale[i];
217 } else {
218 ✗ _yScale[i] = 1;
219 }
220 }
221
222
223 // Create Kinsol memory
224 /* Create the SUNDIALS context first - every other SUNDIALS object below is
225 created with it. */
226 ✗ if (SUNContext_Create(SUN_COMM_NULL, &_sunctx) != SUN_SUCCESS)
227 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER,"SUNDIALS_ERROR: SUNContext_Create failed");
228 /* Mute SUNDIALS' own logger: package level messages go to stderr/stdout by default, and we report solver failures ourselves. */
229 {
230 ✗ SUNLogger _sunlogger = NULL;
231 ✗ if (SUNContext_GetLogger(_sunctx, &_sunlogger) == SUN_SUCCESS && _sunlogger != NULL) {
232 ✗ SUNLogger_SetErrorFilename(_sunlogger, "");
233 ✗ SUNLogger_SetWarningFilename(_sunlogger, "");
234 ✗ SUNLogger_SetInfoFilename(_sunlogger, "");
235 ✗ SUNLogger_SetDebugFilename(_sunlogger, "");
236 }
237 }
238
239 ✗ _Kin_y = N_VMake_Serial(_dimSys, _y, _sunctx);
240 ✗ _Kin_y0 = N_VMake_Serial(_dimSys, _y0, _sunctx);
241 ✗ _Kin_yScale = N_VMake_Serial(_dimSys, _yScale, _sunctx);
242 ✗ _Kin_fScale = N_VMake_Serial(_dimSys, _fScale, _sunctx);
243 ✗ _Kin_ySolver = N_VNew_Serial(_dimSys, _sunctx);
244 ✗ _kinMem = KINCreate(_sunctx);
245
246 // Set internal memory
247 ✗ idid = KINInit(_kinMem, kin_fCallback, _Kin_y);
248 ✗ if (check_flag(&idid, (char *)"KINInit", 1)) {
249 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER,"Kinsol::initialize()");
250 }
251
252 // Initialize dense linear solver
253 ✗ _Kin_J = SUNDenseMatrix(_dimSys, _dimSys, _sunctx);
254 ✗ _Kin_linSol = SUNLinSol_Dense(_Kin_ySolver, _Kin_J, _sunctx);
255 ✗ if (_Kin_linSol == NULL) {
256 ✗ fprintf(stderr,"\nSUNDIALS_ERROR: SUNLinSol_Dense() failed - returned NULL pointer\n\n");
257 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "Kinsol::solve()");
258 }
259 ✗ idid = KINSetLinearSolver(_kinMem, _Kin_linSol, _Kin_J);
260 ✗ if (check_flag(&idid, (char *)"KINSetLinearSolver", 1)) {
261 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "Kinsol::initialize()");
262 }
263
264 // Set optional inputs
265 ✗ idid = KINSetUserData(_kinMem, _data);
266 ✗ if (check_flag(&idid, (char *)"KINSetUserData", 1)) {
267 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER,"Kinsol::initialize()");
268 }
269
270 ✗ idid = KINSetNumMaxIters(_kinMem, 50);
271
272 ✗ _fnormtol = 1.e-13; /* function tolerance */
273 ✗ _scsteptol = 1.e-13; /* step tolerance */
274
275 ✗ idid = KINSetFuncNormTol(_kinMem, _fnormtol);
276 ✗ idid = KINSetScaledStepTol(_kinMem, _scsteptol);
277 ✗ idid = KINSetRelErrFunc(_kinMem, 1e-14);
278
279 ✗ _counter = 0;
280
281 ✗ LOGGER_WRITE("Kinsol: initialized",LC_NLS,LL_DEBUG);
282 ✗ }
283
284
285 /**
286 * @brief Wrapper for Kinsol::solve()
287 *
288 * Sets algLoop and first_solve
289 *
290 * @param algLoop
291 * @param first_solve
292 */
293 ✗ void Kinsol::solve(shared_ptr<INonLinearAlgLoop> algLoop, bool first_solve)
294 {
295 ✗ if (first_solve) {
296 _algLoop = algLoop;
297 ✗ _firstCall = true;
298 }
299 ✗ if (_algLoop != algLoop) {
300 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "algloop system is not initialized");
301 }
302 ✗ solve();
303 ✗ }
304
305
306 /**
307 * @brief Solve algebraic loop with former initialzed Kinsol solver
308 *
309 * _algLoop and _firstCall has to be set outside this function.
310 */
311 ✗ void Kinsol::solve() {
312 // Initialize at first call
313 ✗ if (_firstCall) {
314 ✗ initialize();
315 }
316
317 ✗ if(!_algLoop) {
318 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "algloop system is not initialized");
319 }
320
321 int idid;
322 ✗ _counter++;
323 ✗ _eventRetry = false;
324 ✗ _iterationStatus = CONTINUE;
325
326 //get variables vectors for last accepted step
327 ✗ _algLoop->getReal(_y);
328 ✗ _algLoop->getRealStartValues(_y0);
329
330 // Reinitialize dense linear solver if last call was with comlpete pivoting or iterative solver
331 ✗ if(_usedCompletePivoting || _usedIterativeSolver)
332 {
333 ✗ idid = SUNLinSolFree(_Kin_linSol);
334 ✗ if (check_flag(&idid, (char *)"SUNLinSolFree", 1)) {
335 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "Kinol::solve()");
336 }
337 ✗ _Kin_linSol = SUNLinSol_Dense(_Kin_ySolver, _Kin_J, _sunctx);
338 ✗ if (_Kin_linSol == NULL) {
339 ✗ fprintf(stderr,"\nSUNDIALS_ERROR: SUNLinSol_Dense() failed - returned NULL pointer\n\n");
340 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "Kinsol::solve()");
341 }
342 ✗ idid = KINSetLinearSolver(_kinMem, _Kin_linSol, _Kin_J);
343 ✗ if (check_flag(&idid, (char *)"KINSetUserData", 1)) {
344 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "Kinsol::initialize()");
345 }
346
347 ✗ _usedCompletePivoting = false;
348 ✗ _usedIterativeSolver = false;
349 }
350
351 // Reset Scaling
352 ✗ for(int i=0;i<_dimSys;i++) {
353 ✗ _fScale[i] = 1.0;
354 }
355
356 // Try Dense first
357 ////////////////////////////
358 ✗ solveNLS();
359 ✗ if(_iterationStatus == DONE)
360 {
361 ✗ _algLoop->setReal(_y);
362 ✗ _algLoop->evaluate();
363 ✗ return;
364 }
365
366 // Try Dense with scaling
367 ////////////////////////////
368 ✗ _iterationStatus = CONTINUE;
369 ✗ _algLoop->setReal(_y0);
370 ✗ _algLoop->evaluate();
371 ✗ _algLoop->getRHS(_fScale);
372 ✗ for(int i=0;i<_dimSys;i++) {
373 ✗ if(abs(_fScale[i]) >1.0) {
374 ✗ _fScale[i] = abs(1/_fScale[i]);
375 } else {
376 ✗ _fScale[i] = 1;
377 }
378 }
379 ✗ solveNLS();
380 ✗ if(_iterationStatus == DONE) {
381 ✗ _algLoop->setReal(_y);
382 ✗ _algLoop->evaluate();
383 ✗ return;
384 }
385
386 // Try SPGMR solver
387 // Scaled, Preconditioned, Generalized Minimum Residual iterative linear solver
388 /////////////////////////////////
389 ✗ _iterationStatus = CONTINUE;
390 ✗ _usedIterativeSolver = true;
391
392 // Reset Scaling
393 ✗ for(int i=0;i<_dimSys;i++) {
394 ✗ _fScale[i] = 1.0;
395 }
396
397 // Free linear solver and initialize linear solver
398 ✗ idid = SUNLinSolFree(_Kin_linSol);
399 ✗ if (check_flag(&idid, (char *)"SUNLinSolFree", 1)) {
400 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "Kinsol::solve()");
401 }
402 ✗ _Kin_linSol = SUNLinSol_SPGMR(_Kin_ySolver, SUN_PREC_NONE, _dimSys /* Krylov subspaces */, _sunctx);
403 ✗ if (_Kin_linSol == NULL) {
404 ✗ fprintf(stderr,"\nSUNDIALS_ERROR: SUNLinSol_SPGMR() failed - returned NULL pointer\n\n");
405 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "Kinsol::solve()");
406 }
407 ✗ idid = KINSetLinearSolver(_kinMem, _Kin_linSol, NULL);
408 ✗ if (check_flag(&idid, (char *)"KINSetLinearSolver", 1)) {
409 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "Kinsol::solve()");
410 }
411
412 // Solve
413 ✗ solveNLS();
414 ✗ if(_iterationStatus == DONE) {
415 ✗ _algLoop->setReal(_y);
416 ✗ _algLoop->evaluate();
417 ✗ return;
418 }
419
420 // Try SPGMR solver with scaling
421 /////////////////////////////////
422 ✗ _iterationStatus = CONTINUE;
423 ✗ _algLoop->setReal(_y0);
424 ✗ _algLoop->evaluate();
425 ✗ _algLoop->getRHS(_fScale);
426 ✗ for(int i=0;i<_dimSys;i++) {
427 ✗ if(abs(_fScale[i]) >1.0) {
428 ✗ _fScale[i] = abs(1/_fScale[i]);
429 } else {
430 ✗ _fScale[i] = 1;
431 }
432 }
433 ✗ solveNLS();
434 ✗ if(_iterationStatus == DONE) {
435 ✗ _algLoop->setReal(_y);
436 ✗ _algLoop->evaluate();
437 ✗ return;
438 }
439
440 // Try SPBCG solver
441 // Scaled, Preconditioned, Generalized Minimum Residual iterative linear solver
442 /////////////////////////////////
443 ✗ _iterationStatus = CONTINUE;
444
445 // Reset Scaling
446 ✗ for(int i=0;i<_dimSys;i++) {
447 ✗ _fScale[i] = 1.0;
448 }
449
450 // Free linear solver and initialize linear solver
451 ✗ idid = SUNLinSolFree(_Kin_linSol);
452 ✗ if (check_flag(&idid, (char *)"SUNLinSolFree", 1)) {
453 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "Kinsol::solve()");
454 }
455 ✗ _Kin_linSol = SUNLinSol_SPBCGS(_Kin_ySolver, SUN_PREC_NONE, _dimSys /* Krylov subspaces */, _sunctx);
456 ✗ if (_Kin_linSol == NULL) {
457 ✗ fprintf(stderr,"\nSUNDIALS_ERROR: SUNLinSol_SPGMR() failed - returned NULL pointer\n\n");
458 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "Kinsol::solve()");
459 }
460 ✗ idid = KINSetLinearSolver(_kinMem, _Kin_linSol, _Kin_J);
461 ✗ if (check_flag(&idid, (char *)"KINSetLinearSolver", 1)) {
462 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "Kinsol::solve()");
463 }
464
465 // Solve
466 ✗ solveNLS();
467 ✗ if(_iterationStatus == DONE) {
468 ✗ _algLoop->setReal(_y);
469 ✗ _algLoop->evaluate();
470 ✗ return;
471 }
472
473 // Try SPBCG solver with scaling
474 /////////////////////////////////
475 ✗ _iterationStatus = CONTINUE;
476 ✗ _algLoop->setReal(_y0);
477 ✗ _algLoop->evaluate();
478 ✗ _algLoop->getRHS(_fScale);
479 ✗ for(int i=0;i<_dimSys;i++) {
480 ✗ if(abs(_fScale[i]) >1.0) {
481 ✗ _fScale[i] = abs(1/_fScale[i]);
482 } else {
483 ✗ _fScale[i] = 1;
484 }
485 }
486 ✗ solveNLS();
487 ✗ if(_iterationStatus == DONE) {
488 ✗ _algLoop->setReal(_y);
489 ✗ _algLoop->evaluate();
490 ✗ return;
491 }
492
493 // TODO: Whats this event stuff doing????
494 ✗ if(_eventRetry) {
495 ✗ memcpy(_y, _helpArray ,_dimSys*sizeof(double));
496 ✗ _iterationStatus = CONTINUE;
497 ✗ return;
498 }
499
500 // Give up
501 /////////////////////////////////
502 ✗ if(_iterationStatus == SOLVERERROR && !_eventRetry) {
503 ✗ if(_kinsolSettings->getContinueOnError()) {
504 ✗ if(!_solverErrorNotificationGiven) {
505 ✗ LOGGER_WRITE("Kinsol: Solver error detected. The simulation will continue, but the results may be incorrect.",LC_NLS,LL_WARNING);
506 ✗ _solverErrorNotificationGiven = true;
507 }
508 } else {
509 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER,"Nonlinear solver failed!");
510 }
511 }
512 }
513
514 ✗ bool* Kinsol::getConditionsWorkArray() {
515 ✗ return AlgLoopSolverDefaultImplementation::getConditionsWorkArray();
516 }
517
518
519 ✗ bool* Kinsol::getConditions2WorkArray() {
520 ✗ return AlgLoopSolverDefaultImplementation::getConditions2WorkArray();
521 }
522
523
524 ✗ double* Kinsol::getVariableWorkArray() {
525 ✗ return AlgLoopSolverDefaultImplementation::getVariableWorkArray();
526 }
527
528 ✗ INonLinearAlgLoopSolver::ITERATIONSTATUS Kinsol::getIterationStatus() {
529 ✗ return _iterationStatus;
530 }
531
532
533 ✗ void Kinsol::calcFunction(const double *y, double *residual) {
534 ✗ if(!_algLoop) {
535 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "algloop system is not initialized");
536 }
537 ✗ _fValid = true;
538 ✗ _algLoop->setReal(y);
539 try {
540 ✗ _algLoop->evaluate();
541 ✗ } catch (std::exception & ex) {
542 ✗ _fValid = false;
543 ✗ }
544 ✗ _algLoop->getRHS(residual);
545
546 // Check for numerical overflow to avoid endless iterations
547 ✗ for(int i=0;i<_dimSys;i++) {
548 ✗ if(!isfinite(residual[i]) || !isfinite(y[i])) {
549 ✗ _fValid = false;
550 }
551 }
552 ✗ }
553
554 ✗ int Kinsol::kin_f(N_Vector y,N_Vector fval, void *user_data) {
555 ✗ ((Kinsol*) user_data)->calcFunction(NV_DATA_S(y),NV_DATA_S(fval));
556
557 ✗ if(((Kinsol*) user_data)->_fValid) {
558 return(0);
559 } else {
560 ✗ return(1);
561 }
562 }
563
564 ✗ void Kinsol::stepCompleted(double time) {
565 ✗ memcpy(_y0,_y,_dimSys*sizeof(double));
566 ✗ memcpy(_y_old,_y_new,_dimSys*sizeof(double));
567 ✗ memcpy(_y_new,_y,_dimSys*sizeof(double));
568 ✗ }
569
570 ✗ int Kinsol::check_flag(void *flagvalue, char *funcname, int opt) {
571 int *errflag;
572
573 // Check if SUNDIALS function returned NULL pointer - no memory allocated
574 ✗ if (opt == 0 && flagvalue == NULL) {
575 ✗ fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n", funcname);
576 ✗ return(1);
577 }
578
579 // Check if flag < 0
580 ✗ else if (opt == 1) {
581 errflag = (int *) flagvalue;
582 ✗ if (*errflag < 0) {
583 ✗ fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n", funcname, *errflag);
584 ✗ return(1);
585 }
586 }
587
588 // Check if function returned NULL pointer - no memory allocated
589 ✗ else if (opt == 2 && flagvalue == NULL) {
590 ✗ fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n", funcname);
591 ✗ return(1);
592 }
593
594 return(0);
595 }
596
597 ✗ void Kinsol::solveNLS() {
598 int
599 method = KIN_NONE,
600 iter = 0,
601 idid;
602 double
603 maxStepsStart = 0,
604 maxSteps = maxStepsStart,
605 maxStepsHigh =1e8,
606 locTol =5e-7,
607 delta = 1e-14;
608
609 ✗ _currentIterateNorm = 100.0;
610
611 ✗ while(_iterationStatus == CONTINUE) {
612 iter++;
613 //Increase max. Newton Step size
614 ✗ idid = KINSetMaxNewtonStep(_kinMem, maxSteps);
615
616 // Reset initial guess
617 ✗ memcpy(_y,_y0,_dimSys*sizeof(double));
618
619 // Call Kinsol
620 ✗ idid = KINSol(_kinMem, _Kin_y, method, _Kin_yScale, _Kin_fScale);
621
622 ✗ KINGetFuncNorm(_kinMem, &_fnorm);
623 ✗ if(!_fValid && (idid==KIN_FIRST_SYSFUNC_ERR)) {
624 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER,"Algloop could not be evaluated! Evaluation failed at the first call.");
625 }
626 ✗ if(idid != KIN_SUCCESS && _fnorm < locTol && _fnorm < _currentIterateNorm) {
627 ✗ _currentIterateNorm = _fnorm;
628 ✗ memcpy(_currentIterate,_y,_dimSys*sizeof(double));
629 }
630
631 // Check the return status for possible restarts
632 ✗ switch (idid) {
633 // Success
634 ✗ case KIN_SUCCESS:
635 ✗ _iterationStatus = DONE;
636 ✗ break;
637
638 // Fine initial guess
639 ✗ case KIN_INITIAL_GUESS_OK:
640 ✗ _iterationStatus = DONE;
641 ✗ break;
642
643 // Did we reach a saddle point
644 ✗ case KIN_STEP_LT_STPTOL:
645 ✗ KINGetFuncNorm(_kinMem, &_fnorm);
646 ✗ if(_fnorm < _fnormtol) {
647 ✗ _iterationStatus = DONE;
648 } else {
649 ✗ check4EventRetry(_y);
650 ✗ if(method==KIN_NONE) {
651 method = KIN_LINESEARCH;
652 maxSteps = maxStepsStart;
653 } else {
654 ✗ _iterationStatus = SOLVERERROR;
655 }
656 }
657 break;
658
659 // MaxStep too low
660 ✗ case KIN_MXNEWT_5X_EXCEEDED:
661 ✗ KINGetFuncNorm(_kinMem, &_fnorm);
662 ✗ if(_fnorm < _fnormtol) {
663 ✗ _iterationStatus = DONE;
664 } else {
665 ✗ check4EventRetry(_y);
666 ✗ if(method == KIN_NONE) {
667 ✗ if (maxSteps == maxStepsHigh) {
668 method = KIN_LINESEARCH;
669 maxSteps = maxStepsStart;
670 } else {
671 maxSteps = maxStepsHigh;
672 }
673 } else { // already trying Linesearch
674 ✗ _iterationStatus = SOLVERERROR;
675 }
676 }
677 break;
678
679 // Max Iterations exceeded
680 ✗ case KIN_MAXITER_REACHED:
681 ✗ KINGetFuncNorm(_kinMem, &_fnorm);
682 ✗ if(_fnorm < _fnormtol) {
683 ✗ _iterationStatus = DONE;
684 } else {
685 ✗ check4EventRetry(_y);
686 ✗ if(method == KIN_NONE) {
687 method = KIN_LINESEARCH;
688 maxSteps = maxStepsStart;
689 } else { // already trying Linesearch
690 ✗ if (maxSteps > 0 && maxSteps < 1) {
691 ✗ _iterationStatus = SOLVERERROR;
692 } else {
693 ✗ if (maxSteps==0) {
694 maxSteps = maxStepsHigh;
695 } else {
696 ✗ maxSteps /= 10;
697 }
698 }
699 }
700 }
701 break;
702
703 // Linesearch did not converge
704 ✗ case KIN_LINESEARCH_NONCONV:
705 ✗ KINGetFuncNorm(_kinMem, &_fnorm);
706 ✗ if(_fnorm < _fnormtol) {
707 ✗ _iterationStatus = DONE;
708 } else {
709 ✗ check4EventRetry(_y);
710 ✗ if(delta < 1e-16) {
711 ✗ _iterationStatus = SOLVERERROR;
712 } else {
713 ✗ delta /= 1e2;
714 ✗ idid = KINSetRelErrFunc(_kinMem, delta);
715 }
716 }
717 break;
718
719 // Other failures (setup etc) -> directly break
720 ✗ default:
721 ✗ KINGetFuncNorm(_kinMem, &_fnorm);
722 ✗ if(_fnorm < _fnormtol) { // Initial guess may be the solution
723 ✗ _iterationStatus = DONE;
724 } else {
725 ✗ _iterationStatus = SOLVERERROR;
726 }
727 break;
728 }
729 }
730
731 // Check if the best found solution suffices
732 ✗ if(_iterationStatus == SOLVERERROR && _currentIterateNorm < locTol) {
733 ✗ _iterationStatus = DONE;
734 ✗ for(int i=0;i<_dimSys;i++) {
735 ✗ _y[i] = _currentIterate[i];
736 }
737 }
738 ✗ }
739
740
741 /**
742 * \brief Restores all algloop variables for a output step
743 * \return Return_Description
744 * \details Details
745 */
746 ✗ void Kinsol::restoreOldValues() {
747 ✗ memcpy(_y,_y_old,_dimSys*sizeof(double));
748 ✗ }
749
750
751 /**
752 * \brief Restores all algloop variables for last output step
753 * \return Return_Description
754 * \details Details
755 */
756 ✗ void Kinsol::restoreNewValues() {
757 ✗ memcpy(_y,_y_new,_dimSys*sizeof(double));
758 ✗ }
759
760
761 ✗ void Kinsol::check4EventRetry(double* y) {
762 ✗ if(!_algLoop) {
763 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "algloop system is not initialized");
764 }
765 ✗ _algLoop->setReal(y);
766 ✗ if(!(_algLoop->isConsistent()) && !_eventRetry) {
767 ✗ memcpy(_helpArray, y,_dimSys*sizeof(double));
768 ✗ _eventRetry = true;
769 }
770 ✗ }
771 /** @} */ // end of solverKinsol
772