Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 82.3% 200 / 0 / 243
Functions: 64.3% 9 / 0 / 14
Branches: 26.5% 161 / 0 / 607
OMCompiler/SimulationRuntime/cpp/Solver/Newton/Newton.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 solverNewton
29 *
30 * @{
31 */
32
33 #include <Core/ModelicaDefine.h>
34 #include <Core/Modelica.h>
35
36 #include <Solver/Newton/FactoryExport.h>
37 #include <Core/Utils/extension/logger.hpp>
38
39 #include <Solver/Newton/Newton.h>
40
41 #include <Core/Math/ILapack.h> // needed for solution of linear system with Lapack
42 #include <Core/Math/Constants.h> // definitializeion of constants like uround
43
44 52 Newton::Newton(INonLinSolverSettings* settings,shared_ptr<INonLinearAlgLoop> algLoop)
45 :AlgLoopSolverDefaultImplementation()
46 ,_algLoop (algLoop)
47 52 , _newtonSettings ((INonLinSolverSettings*)settings)
48 52 , _yNames (NULL)
49 52 , _yNominal (NULL)
50 52 , _yMin (NULL)
51 52 , _yMax (NULL)
52 52 , _y (NULL)
53 52 , _yHelp (NULL)
54 52 , _yTest (NULL)
55 52 , _fNominal (NULL)
56 52 , _f (NULL)
57 52 , _fHelp (NULL)
58 52 , _fTest (NULL)
59 52 , _iHelp (NULL)
60 52 , _jHelp (NULL)
61 52 , _jac (NULL)
62 52 , _firstCall (true)
63 52 , _iterationStatus (CONTINUE)
64
1/2
✓ Branch 2 taken 52 times.
✗ Branch 3 not taken.
 52 , _lc (LC_NLS)
65 {
66
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_algLoop)
67 {
68
3/3
✓ Branch 1 taken 52 times.
✓ Branch 4 taken 52 times.
✓ Branch 7 taken 52 times.
 52 AlgLoopSolverDefaultImplementation::initialize(_algLoop->getDimZeroFunc(),_algLoop->getDimReal());
69 }
70 else
71 {
72 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "solve for single instance is not supported");
73 }
74 52 }
75
76 104 Newton::~Newton()
77 {
78
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_yNames) delete [] _yNames;
79
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_yNominal) delete [] _yNominal;
80
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_yMin) delete [] _yMin;
81
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_yMax) delete [] _yMax;
82
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_y) delete [] _y;
83
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_yHelp) delete [] _yHelp;
84
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_yTest) delete [] _yTest;
85
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_fNominal) delete [] _fNominal;
86
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_f) delete [] _f;
87
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_fHelp) delete [] _fHelp;
88
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_fTest) delete [] _fTest;
89
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_iHelp) delete [] _iHelp;
90
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_jHelp) delete [] _jHelp;
91
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_jac) delete [] _jac;
92 104 }
93
94 52 void Newton::initialize()
95 {
96
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 _firstCall = false;
97
98 //(Re-) initializeialization of algebraic loop
99
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if(_algLoop)
100 52 _algLoop->initialize();
101 else
102 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "algloop system is not initialized");
103
104
105 52 _dimSys = _algLoop->getDimReal();
106
107
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 if (_dimSys > 0) {
108 // initialize of vectors of unknowns and residuals
109
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if (_yNames) delete [] _yNames;
110
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if (_yNominal) delete [] _yNominal;
111
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if (_yMin) delete [] _yMin;
112
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if (_yMax) delete [] _yMax;
113
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if (_y) delete [] _y;
114
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if (_fNominal) delete [] _fNominal;
115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if (_f) delete [] _f;
116
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if (_yHelp) delete [] _yHelp;
117
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if (_yTest) delete [] _yTest;
118
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if (_fHelp) delete [] _fHelp;
119
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if (_fTest) delete [] _fTest;
120
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if (_iHelp) delete [] _iHelp;
121
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if (_jHelp) delete [] _jHelp;
122
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if (_jac) delete [] _jac;
123
124
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 _yNames = new const char* [_dimSys];
125
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 _yNominal = new double[_dimSys];
126
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 _yMin = new double[_dimSys];
127
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 _yMax = new double[_dimSys];
128
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 _y = new double[_dimSys];
129
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 _fNominal = new double[_dimSys];
130
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 _f = new double[_dimSys];
131
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 _yHelp = new double[_dimSys];
132
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 _yTest = new double[_dimSys];
133
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 _fHelp = new double[_dimSys];
134
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 _fTest = new double[_dimSys];
135
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 _iHelp = new long int[_dimSys];
136
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 _jHelp = new long int[_dimSys];
137
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 _jac = new double[_dimSys*_dimSys];
138
139 52 _algLoop->getNamesReal(_yNames);
140 52 _algLoop->getNominalReal(_yNominal);
141 52 _algLoop->getMinReal(_yMin);
142 52 _algLoop->getMaxReal(_yMax);
143 }
144
145
146
147
2/16
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 52 LOGGER_WRITE_BEGIN("Newton: eq" + to_string(_algLoop->getEquationIndex()) +
148 " initialized", _lc, LL_DEBUG);
149
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 2 times.
 52 LOGGER_WRITE_VECTOR("yNames", _yNames, _dimSys, _lc, LL_DEBUG);
150
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 2 times.
 52 LOGGER_WRITE_VECTOR("yNominal", _yNominal, _dimSys, _lc, LL_DEBUG);
151
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 2 times.
 52 LOGGER_WRITE_END(_lc, LL_DEBUG);
152 52 }
153 ✗ void Newton::solve( shared_ptr<INonLinearAlgLoop> algLoop,bool first_solve)
154 {
155 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "solve for single instance is not supported");
156 }
157
158
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 759169 times.
 759169 void Newton::solve( )
159 {
160
161
162
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 759169 times.
 759169 if(!_algLoop)
163 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "algloop system is not initialized");
164 long int
165 759169 dimRHS = 1, // Dimension of right hand side of linear system (=b)
166 759169 info = 0; // Retrun-flag of Fortran code
167 int
168 totSteps = 0; // Total number of steps taken
169 double
170 759169 atol = _newtonSettings->getAtol(),
171 759169 rtol = _newtonSettings->getRtol();
172
173 // If initialize() was not called yet
174
2/2
✓ Branch 0 taken 52 times.
✓ Branch 1 taken 759117 times.
 759169 if (_firstCall)
175 52 initialize();
176
177 // Get current values and residuals from system
178 759169 _algLoop->getReal(_y);
179 759169 _algLoop->evaluate();
180 759169 _algLoop->getRHS(_f);
181
182 // Reset status flag
183 759169 _iterationStatus = CONTINUE;
184
185
2/36
✓ Branch 0 taken 758007 times.
✓ Branch 1 taken 1162 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
 759169 LOGGER_WRITE_BEGIN("Newton: eq" + to_string(_algLoop->getEquationIndex()) +
186 " at time " + to_string(_algLoop->getSimTime()) + ":",
187 _lc, LL_DEBUG);
188
189
2/2
✓ Branch 0 taken 768026 times.
✓ Branch 1 taken 759169 times.
 1527195 while (_iterationStatus == CONTINUE) {
190
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 768026 times.
 768026 if (totSteps >= _newtonSettings->getNewtMax()) {
191 ✗ LOGGER_WRITE_END(_lc, LL_DEBUG);
192 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER,
193 ✗ "error solving nonlinear system (iteration limit: " + to_string(totSteps) + ")");
194 }
195
196 // Newton step for non-linear system
197
2/10
✓ Branch 0 taken 766864 times.
✓ Branch 1 taken 1162 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 768026 LOGGER_WRITE_VECTOR("y" + to_string(totSteps), _y, _dimSys, _lc, LL_DEBUG);
198
2/10
✓ Branch 0 taken 766864 times.
✓ Branch 1 taken 1162 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 768026 LOGGER_WRITE_VECTOR("f" + to_string(totSteps), _f, _dimSys, _lc, LL_DEBUG);
199
200 768026 calcJacobian(_jac, _fNominal);
201
202 // Initialize line search function
203 double phi = 0.0;
204
2/2
✓ Branch 0 taken 1134857 times.
✓ Branch 1 taken 768026 times.
 1902883 for (int i = 0; i < _dimSys; i++) {
205 1134857 _f[i] /= _fNominal[i];
206 1134857 phi += _f[i] * _f[i];
207 }
208
209 double det;
210
3/4
✓ Branch 0 taken 405273 times.
✓ Branch 1 taken 362753 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 405273 times.
 768026 if (_dimSys == 1 && (det = _jac[0]) != 0.0) {
211 det = _jac[0];
212 405273 _f[0] /= det;
213 405273 info = 0;
214 }
215
3/4
✓ Branch 0 taken 360927 times.
✓ Branch 1 taken 1826 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 360927 times.
 362753 else if (_dimSys == 2 && (det = _jac[0]*_jac[3] - _jac[1]*_jac[2]) != 0.0) {
216 det = _jac[0]*_jac[3] - _jac[1]*_jac[2];
217 360927 double f0 = (_f[0]*_jac[3] - _f[1]*_jac[2]) / det;
218 360927 _f[1] = (_jac[0]*_f[1] - _jac[1]*_f[0]) / det;
219 360927 _f[0] = f0;
220 360927 info = 0;
221 }
222 else {
223 // Solve linear system
224 1826 dgesv_(&_dimSys, &dimRHS, _jac, &_dimSys, _iHelp, _f, &_dimSys, &info);
225 }
226
227
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 768026 times.
 768026 if (info > 0) {
228 ✗ long int info2 = 0;
229 ✗ double scale = 0.0;
230 ✗ dgetc2_(&_dimSys, _jac, &_dimSys, _iHelp, _jHelp, &info2);
231 ✗ dgesc2_(&_dimSys, _jac, &_dimSys, _f, _iHelp, _jHelp, &scale);
232 // limit change of y by f to yNominal in case of singular Jacobian
233 ✗ for (int i = 0; i < _dimSys; i++) {
234 ✗ _f[i] = std::min(_yNominal[i], std::max(-_yNominal[i], _f[i]));
235 }
236 ✗ LOGGER_WRITE("total pivoting: dgesv/dgetc2 infos: " + to_string(info) + "/" + to_string(info2) +
237 ", dgesc2 scale: " + to_string(scale), _lc, LL_DEBUG);
238 }
239
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 768026 times.
 768026 else if (info < 0) {
240 ✗ LOGGER_WRITE_END(_lc, LL_DEBUG);
241 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER,
242 ✗ "error solving nonlinear system (iteration: " + to_string(totSteps)
243 ✗ + ", dgesv info: " + to_string(info) + ")");
244 }
245
246 // Increase counter
247 768026 ++ totSteps;
248
249 // New iterate
250 double lambda = 1.0; // step size
251 double alpha = 1e-4; // guard for sufficient decrease
252 // first find a feasible step
253 while (true) {
254
2/2
✓ Branch 0 taken 1134857 times.
✓ Branch 1 taken 768026 times.
 1902883 for (int i = 0; i < _dimSys; i++) {
255 1134857 _yHelp[i] = _y[i] - lambda * _f[i] /** _yNominal[i]*/;
256
2/4
✓ Branch 0 taken 1134857 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1134857 times.
✗ Branch 3 not taken.
 3404571 _yHelp[i] = std::min(_yMax[i], std::max(_yMin[i], _yHelp[i]));
257 }
258 // evaluate function
259 try {
260
1/1
✓ Branch 1 taken 768026 times.
 768026 calcFunction(_yHelp, _fHelp);
261 }
262 ✗ catch (ModelicaSimulationError& ex) {
263 ✗ if (lambda < 1e-10) {
264 ✗ LOGGER_WRITE_END(_lc, LL_DEBUG);
265 ✗ throw;
266 }
267 // reduce step size
268 ✗ lambda *= 0.5;
269 continue;
270 ✗ }
271 break;
272 }
273 // check stopping criterion
274 768026 _iterationStatus = DONE;
275
2/2
✓ Branch 0 taken 1134345 times.
✓ Branch 1 taken 759169 times.
 1893514 for (int i = 0; i < _dimSys; i++) {
276
2/2
✓ Branch 0 taken 8857 times.
✓ Branch 1 taken 1125488 times.
 1134345 if (std::abs(_fHelp[i]) > atol + rtol * _fNominal[i]) {
277 8857 _iterationStatus = CONTINUE;
278 8857 break;
279 }
280 }
281 // second do line search with quadratic approximation of phi(lambda)
282 // C.T.Kelley: Solving Nonlinear Equations with Newton's Method,
283 // no 1 in Fundamentals of Algorithms, SIAM 2003. ISBN 0-89871-546-6.
284 double phiHelp = 0.0;
285
2/2
✓ Branch 0 taken 1134857 times.
✓ Branch 1 taken 768026 times.
 1902883 for (int i = 0; i < _dimSys; i++) {
286 1134857 _fHelp[i] /= _fNominal[i];
287 1134857 phiHelp += _fHelp[i] * _fHelp[i];
288 }
289
2/2
✓ Branch 0 taken 8859 times.
✓ Branch 1 taken 759169 times.
 768028 while (_iterationStatus == CONTINUE) {
290 // test half step that also serves as max bound for step reduction
291 8859 double lambdaTest = 0.5*lambda;
292
2/2
✓ Branch 0 taken 9505 times.
✓ Branch 1 taken 8859 times.
 18364 for (int i = 0; i < _dimSys; i++) {
293 9505 _yTest[i] = _y[i] - lambdaTest * _f[i] /** _yNominal[i]*/;
294
2/4
✓ Branch 0 taken 9505 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 9505 times.
✗ Branch 3 not taken.
 28515 _yTest[i] = std::min(_yMax[i], std::max(_yMin[i], _yTest[i]));
295 }
296 8859 calcFunction(_yTest, _fTest);
297 double phiTest = 0.0;
298
2/2
✓ Branch 0 taken 9505 times.
✓ Branch 1 taken 8859 times.
 18364 for (int i = 0; i < _dimSys; i++) {
299 9505 _fTest[i] /= _fNominal[i];
300 9505 phiTest += _fTest[i] * _fTest[i];
301 }
302 // check for sufficient decrease of phiHelp
303 // and no further decrease with phiTest
304 // otherwise minimize quadratic approximation of phi(lambda)
305
2/2
✓ Branch 0 taken 8803 times.
✓ Branch 1 taken 56 times.
 8859 if (phiHelp > (1.0 - alpha * lambda) * phi ||
306
2/2
✓ Branch 0 taken 26 times.
✓ Branch 1 taken 8777 times.
 8803 phiTest < (1.0 - alpha * lambda) * phiHelp) {
307 82 long int n = 3;
308 long int ipiv[3];
309 82 double bx[] = {phi, phiTest, phiHelp};
310 82 double A[] = {1.0, 1.0, 1.0,
311 0.0, lambdaTest, lambda,
312 82 0.0, lambdaTest*lambdaTest, lambda*lambda};
313 82 dgesv_(&n, &dimRHS, A, &n, ipiv, bx, &n, &info);
314
2/2
✓ Branch 0 taken 80 times.
✓ Branch 1 taken 2 times.
 82 lambda = std::max(0.1*lambda, -0.5*bx[1]/bx[2]);
315
2/2
✓ Branch 0 taken 28 times.
✓ Branch 1 taken 54 times.
 82 if (lambda >= lambdaTest) {
316 // upper bound 0.5*lambda
317 lambda = lambdaTest;
318
1/2
✓ Branch 0 taken 28 times.
✗ Branch 1 not taken.
 28 std::copy(_yTest, _yTest + _dimSys, _yHelp);
319
1/2
✓ Branch 0 taken 28 times.
✗ Branch 1 not taken.
 28 std::copy(_fTest, _fTest + _dimSys, _fHelp);
320 phiHelp = phiTest;
321 }
322 else {
323
2/2
✓ Branch 0 taken 54 times.
✓ Branch 1 taken 54 times.
 108 for (int i = 0; i < _dimSys; i++) {
324 54 _yHelp[i] = _y[i] - lambda * _f[i] /** _yNominal[i]*/;
325
2/4
✓ Branch 0 taken 54 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 54 times.
✗ Branch 3 not taken.
 162 _yHelp[i] = std::min(_yMax[i], std::max(_yMin[i], _yHelp[i]));
326 }
327 54 calcFunction(_yHelp, _fHelp);
328 phiHelp = 0.0;
329
2/2
✓ Branch 0 taken 54 times.
✓ Branch 1 taken 54 times.
 108 for (int i = 0; i < _dimSys; i++) {
330 54 _fHelp[i] /= _fNominal[i];
331 54 phiHelp += _fHelp[i] * _fHelp[i];
332 }
333 }
334
1/42
✓ Branch 0 taken 82 times.
✗ Branch 1 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
 82 LOGGER_WRITE("lambda = " + to_string(lambda) +
335 ", phi = " + to_string(phi) +
336 " --> " + to_string(phiHelp),
337 _lc, LL_DEBUG);
338 }
339 // check for stalled step control
340
2/4
✓ Branch 0 taken 8859 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 8859 times.
✗ Branch 3 not taken.
 8859 if (!(lambda >= 1e-10) || !isfinite(phi)) {
341 ✗ LOGGER_WRITE_END(_lc, LL_DEBUG);
342 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER,
343 "can't get sufficient decrease of solution"
344 ✗ ", lambda = " + to_string(lambda) + ", phi = " + to_string(phi));
345 }
346 // check for sufficient decrease
347 // (also break for very small lambda and try a small step instead
348 // -- avoid "sufficient decrease" errors, still have iteration limit)
349
3/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8857 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 8859 if (phiHelp <= (1.0 - alpha * lambda) * phi || (lambda < alpha && isfinite(phiHelp)))
350 break;
351 }
352 // take iterate
353
1/2
✓ Branch 0 taken 768026 times.
✗ Branch 1 not taken.
 768026 std::copy(_yHelp, _yHelp + _dimSys, _y);
354
2/2
✓ Branch 0 taken 1134857 times.
✓ Branch 1 taken 768026 times.
 1902883 for (int i = 0; i < _dimSys; i++)
355 1134857 _f[i] = _fHelp[i] * _fNominal[i];
356 phi = phiHelp;
357 } // end while
358
359
2/2
✓ Branch 0 taken 758007 times.
✓ Branch 1 taken 1162 times.
 759169 LOGGER_WRITE_VECTOR("y*", _y, _dimSys, _lc, LL_DEBUG);
360
2/2
✓ Branch 0 taken 758007 times.
✓ Branch 1 taken 1162 times.
 759169 LOGGER_WRITE_END(_lc, LL_DEBUG);
361 759169 }
362
363 ✗ INonLinearAlgLoopSolver::ITERATIONSTATUS Newton::getIterationStatus()
364 {
365 ✗ return _iterationStatus;
366 }
367
368
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1903552 times.
 1903552 void Newton::calcFunction(const double *y, double *residual)
369 {
370
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1903552 times.
 1903552 if(!_algLoop)
371 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "algloop system is not initialized");
372 1903552 _algLoop->setReal(y);
373 1903552 _algLoop->evaluate();
374 1903552 _algLoop->getRHS(residual);
375 1903552 }
376
377 ✗ void Newton::stepCompleted(double time)
378 {
379
380 ✗ }
381
382
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 768026 times.
 768026 void Newton::calcJacobian(double *jac, double *fNominal)
383 {
384
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 768026 times.
 768026 if(!_algLoop)
385 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "algloop system is not initialized");
386 const double *Adata = NULL;
387 768026 std::fill(fNominal, fNominal + _dimSys, 1e2 * _newtonSettings->getAtol());
388
389 // Use analytic Jacobian if available
390 try {
391
1/1
✓ Branch 1 taken 768026 times.
 768026 const matrix_t& A = _algLoop->getSystemMatrix();
392
3/4
✓ Branch 0 taken 2302 times.
✓ Branch 1 taken 765724 times.
✓ Branch 2 taken 2302 times.
✗ Branch 3 not taken.
 768026 if (A.size1() == _dimSys && A.size2() == _dimSys) {
393 Adata = A.data().begin();
394
1/2
✓ Branch 0 taken 2302 times.
✗ Branch 1 not taken.
 2302 std::copy(Adata, Adata + _dimSys * _dimSys, jac);
395
2/2
✓ Branch 0 taken 8244 times.
✓ Branch 1 taken 2302 times.
 10546 for (int j = 0, idx = 0; j < _dimSys; j++)
396
2/2
✓ Branch 0 taken 31048 times.
✓ Branch 1 taken 8244 times.
 39292 for (int i = 0; i < _dimSys; i++, idx++) {
397
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 31048 times.
 31048 if (!isfinite(jac[idx]))
398 ✗ jac[idx] = 0.0; // remove infinite element in favor of potential singularity
399
2/2
✓ Branch 0 taken 22798 times.
✓ Branch 1 taken 8250 times.
 53846 fNominal[i] = std::max(std::abs(jac[idx]) /** _yNominal[j]*/, fNominal[i]);
400 }
401 }
402 }
403 ✗ catch (ModelicaSimulationError& ex) {
404 ✗ LOGGER_WRITE("Analytic Jacobian failed for eq" +
405 to_string(_algLoop->getEquationIndex()) + " at time " +
406 to_string(_algLoop->getSimTime()) + ": " + ex.what(),
407 _lc, LL_WARNING);
408 ✗ }
409
410 // Alternatively apply finite differences
411
1/2
✓ Branch 0 taken 2302 times.
✗ Branch 1 not taken.
 2302 if (Adata == NULL) {
412
2/2
✓ Branch 0 taken 1126613 times.
✓ Branch 1 taken 765724 times.
 1892337 for (int j = 0; j < _dimSys; j++) {
413 // Reset variables for every column
414
1/2
✓ Branch 0 taken 1126613 times.
✗ Branch 1 not taken.
 1126613 std::copy(_y, _y + _dimSys, _yHelp);
415 1126613 double stepsize = 1e2 * _newtonSettings->getRtol() * _yNominal[j];
416
417 // Finite differences
418 1126613 _yHelp[j] += stepsize;
419
420 1126613 calcFunction(_yHelp, _fHelp);
421
422 // Build Jacobian in Fortran format
423
2/2
✓ Branch 0 taken 1880803 times.
✓ Branch 1 taken 1126613 times.
 3007416 for (int i = 0, idx = j * _dimSys; i < _dimSys; i++, idx++) {
424
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1880803 times.
 1880803 jac[idx] = (_fHelp[i] - _f[i]) / stepsize;
425
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1880803 times.
 1880803 if (!isfinite(jac[idx]))
426 ✗ jac[idx] = 0.0; // remove infinite element in favor of potential singularity
427
2/2
✓ Branch 0 taken 515458 times.
✓ Branch 1 taken 1365345 times.
 2396261 fNominal[i] = std::max(std::abs(jac[idx]) /** _yNominal[j]*/, fNominal[i]);
428 }
429
430 1126613 _yHelp[j] -= stepsize;
431 }
432 }
433
434 // Scale Jacobian
435
2/2
✓ Branch 0 taken 766864 times.
✓ Branch 1 taken 1162 times.
 768026 LOGGER_WRITE_VECTOR("fNominal", fNominal, _dimSys, _lc, LL_DEBUG);
436
2/2
✓ Branch 0 taken 1134857 times.
✓ Branch 1 taken 768026 times.
 1902883 for (int j = 0, idx = 0; j < _dimSys; j++)
437
2/2
✓ Branch 0 taken 1911851 times.
✓ Branch 1 taken 1134857 times.
 3046708 for (int i = 0; i < _dimSys; i++, idx++)
438 //jac[idx] *= _yNominal[j] / fNominal[i];
439 1911851 jac[idx] /= fNominal[i];
440 768026 }
441 759165 bool* Newton::getConditionsWorkArray()
442 {
443 759165 return AlgLoopSolverDefaultImplementation::getConditionsWorkArray();
444
445 }
446 759165 bool* Newton::getConditions2WorkArray()
447 {
448
449 759165 return AlgLoopSolverDefaultImplementation::getConditions2WorkArray();
450 }
451
452
453 759165 double* Newton::getVariableWorkArray()
454 {
455
456 759165 return AlgLoopSolverDefaultImplementation::getVariableWorkArray();
457
458 }
459 ✗ void Newton::restoreOldValues()
460 {
461 ✗ }
462
463 ✗ void Newton::restoreNewValues()
464 {
465 ✗ }
466
467 /** @} */ // end of solverNewton
468