Cantera  3.2.0a1
Loading...
Searching...
No Matches
MultiNewton.cpp
Go to the documentation of this file.
1//! @file MultiNewton.cpp Damped Newton solver for 1D multi-domain problems
2
3// This file is part of Cantera. See License.txt in the top-level directory or
4// at https://cantera.org/license.txt for license and copyright information.
5
8
9#include <ctime>
10
11using namespace std;
12
13namespace Cantera
14{
15
16// unnamed-namespace for local helpers
17namespace
18{
19
20class Indx
21{
22public:
23 Indx(size_t nv, size_t np) : m_nv(nv), m_np(np) {}
24 size_t m_nv, m_np;
25 size_t operator()(size_t m, size_t j) {
26 return j*m_nv + m;
27 }
28};
29
30/**
31 * Return a damping coefficient that keeps the solution after taking one
32 * Newton step between specified lower and upper bounds. This function only
33 * considers one domain.
34 */
35double bound_step(const double* x, const double* step, Domain1D& r, int loglevel)
36{
37 size_t np = r.nPoints();
38 size_t nv = r.nComponents();
39 Indx index(nv, np);
40 double fbound = 1.0;
41 bool wroteTitle = false;
42 string separator = fmt::format("\n {:=>69}", ""); // equals sign separator
43
44 for (size_t m = 0; m < nv; m++) {
45 double above = r.upperBound(m);
46 double below = r.lowerBound(m);
47
48 for (size_t j = 0; j < np; j++) {
49 double val = x[index(m,j)];
50 if (loglevel > 0 && (val > above + 1.0e-12 || val < below - 1.0e-12)) {
51 writelog("\nERROR: solution out of bounds.\n");
52 writelog("domain {:d}: {:>20s}({:d}) = {:10.3e} ({:10.3e}, {:10.3e})\n",
53 r.domainIndex(), r.componentName(m), j, val, below, above);
54 }
55
56 double newval = val + step[index(m,j)];
57
58 if (newval > above) {
59 fbound = std::max(0.0, std::min(fbound,
60 (above - val)/(newval - val)));
61 } else if (newval < below) {
62 fbound = std::min(fbound, (val - below)/(val - newval));
63 }
64
65 if (loglevel > 1 && (newval > above || newval < below)) {
66 if (!wroteTitle){
67 string header = fmt::format(" {:=>10}","") +
68 "Undamped Newton step takes solution out of bounds" +
69 fmt::format("{:=>10}", "");
70 writelog("\n{}", header);
71 writelog(separator);
72 // Split header across 2 lines to shorten the line length
73 writelog("\n {:<7s} {:23s} {:<9s} {:<9s} {:<9s}",
74 "Domain/", "", "Value", "Min", "Max");
75 writelog("\n {:8s} {:<9s} {:<9s} {:6s} {:5s} {:5s}",
76 "Grid Loc", "Component", "Value", "Change", "Bound", "Bound");
77 writelog(separator);
78 wroteTitle = true;
79 }
80 // This create a dynamic spacing to allow for a nicely formatted output in the
81 // Domain/Grid Loc column
82 int domainLength = to_string(r.domainIndex()).length(); // For adjusting spacing
83 int gridLength = to_string(j).length(); // For adjusting spacing
84 int padding = 9; // Total spacing for first column
85 string formatString = fmt::format("{{:<{}d}} / {{:<{}d}}{:>{}}",
86 domainLength, gridLength, "",
87 padding-3-domainLength-gridLength);
88 writelog(fmt::format("\n {}", formatString), r.domainIndex(), j);
89 writelog(" {:<12s} {:>10.3e} {:>10.3e} {:>10.3e} {:>10.3e}",
90 r.componentName(m), val, step[index(m,j)], below, above);
91 }
92 }
93 }
94
95 if (loglevel > 1 && wroteTitle) { // If a title was written, close up the table
96 writelog(separator);
97 }
98 return fbound;
99}
100
101/**
102 * This function computes the square of a weighted norm of a step vector for one
103 * domain.
104 *
105 * @param x Solution vector for this domain.
106 * @param step Newton step vector for this domain.
107 * @param r Object representing the domain. Used to get tolerances,
108 * number of components, and number of points.
109 *
110 * The return value is
111 * @f[
112 * \sum_{n,j} \left(\frac{s_{n,j}}{w_n}\right)^2
113 * @f]
114 * where the error weight for solution component @f$ n @f$ is given by
115 * @f[
116 * w_n = \epsilon_{r,n} \frac{\sum_j |x_{n,j}|}{J} + \epsilon_{a,n}.
117 * @f]
118 * Here @f$ \epsilon_{r,n} @f$ is the relative error tolerance for component n,
119 * and multiplies the average magnitude of solution component n in the domain.
120 * The second term, @f$ \epsilon_{a,n} @f$, is the absolute error tolerance for
121 * component n.
122 *
123 * The purpose is to measure the "size" of the step vector @f$ s @f$ in a way that
124 * takes into account the relative importance or scale of different solution
125 * components.
126 *
127 * Normalization: Each component of the step vector is normalized by a weight that
128 * depends on both the current magnitude of the solution vector and specified
129 * tolerances. This makes the norm dimensionless and scaled appropriately, avoiding
130 * issues where some components dominate due to differences in their scales.
131 */
132double norm_square(const double* x, const double* step, Domain1D& r)
133{
134 double sum = 0.0;
135 double f2max = 0.0;
136 size_t nv = r.nComponents();
137 size_t np = r.nPoints();
138
139 for (size_t n = 0; n < nv; n++) {
140 double esum = 0.0;
141 for (size_t j = 0; j < np; j++) {
142 esum += fabs(x[nv*j + n]);
143 }
144 double ewt = r.rtol(n)*esum/np + r.atol(n);
145 for (size_t j = 0; j < np; j++) {
146 double f = step[nv*j + n]/ewt;
147 sum += f*f;
148 f2max = std::max(f*f, f2max);
149 }
150 }
151 return sum;
152}
153
154} // end unnamed-namespace
155
156
157// ---------------- MultiNewton methods ----------------
158
160 : m_n(sz)
161{
162}
163
164void MultiNewton::resize(size_t sz)
165{
166 m_n = sz;
167 m_x.resize(m_n);
168 m_stp.resize(m_n);
169 m_stp1.resize(m_n);
170}
171
172double MultiNewton::norm2(const double* x, const double* step, OneDim& r) const
173{
174 double sum = 0.0;
175 size_t nd = r.nDomains();
176 for (size_t n = 0; n < nd; n++) {
177 double f = norm_square(x + r.start(n), step + r.start(n), r.domain(n));
178 sum += f;
179 }
180 sum /= r.size();
181 return sqrt(sum);
182}
183
184void MultiNewton::step(double* x, double* step, OneDim& r, int loglevel)
185{
186 r.eval(npos, x, step);
187 for (size_t n = 0; n < r.size(); n++) {
188 step[n] = -step[n];
189 }
190
191 auto jac = r.getJacobian();
192 try {
193 jac->solve(r.size(), step, step);
194 } catch (CanteraError&) {
195 if (jac->info() > 0) {
196 // Positive value for "info" indicates the row where factorization failed
197 size_t row = static_cast<size_t>(jac->info() - 1);
198 // Find the domain, grid point, and solution component corresponding
199 // to this row
200 for (size_t n = 0; n < r.nDomains(); n++) {
201 Domain1D& dom = r.domain(n);
202 size_t nComp = dom.nComponents();
203 if (row >= dom.loc() && row < dom.loc() + nComp * dom.nPoints()) {
204 size_t offset = row - dom.loc();
205 size_t pt = offset / nComp;
206 size_t comp = offset - pt * nComp;
207 throw CanteraError("MultiNewton::step",
208 "Jacobian is singular for domain {}, component {} at point {}\n"
209 "(Matrix row {})",
210 dom.id(), dom.componentName(comp), pt, row);
211 }
212 }
213 }
214 throw;
215 }
216}
217
218double MultiNewton::boundStep(const double* x0, const double* step0, const OneDim& r,
219 int loglevel)
220{
221 double fbound = 1.0;
222 for (size_t i = 0; i < r.nDomains(); i++) {
223 fbound = std::min(fbound,
224 bound_step(x0 + r.start(i), step0 + r.start(i),
225 r.domain(i), loglevel));
226 }
227 return fbound;
228}
229
230int MultiNewton::dampStep(const double* x0, const double* step0,
231 double* x1, double* step1, double& s1,
232 OneDim& r, int loglevel, bool writetitle)
233{
234 // write header
235 if (loglevel > 0 && writetitle) {
236 writelog("\n\n {:-^70}", " Damped Newton iteration ");
237 writelog("\n {:<4s} {:<10s} {:<10s} {:<7s} {:<7s} {:<7s} {:<5s} {:<3s}\n",
238 "Iter", "F_damp", "F_bound", "log(ss)",
239 "log(s0)", "log(s1)", "N_jac", "Age");
240 writelog(" {:->70}", "");
241 }
242
243 // compute the weighted norm of the undamped step size step0
244 double s0 = norm2(x0, step0, r);
245
246 // compute the multiplier to keep all components in bounds
247 double fbound = boundStep(x0, step0, r, loglevel-1);
248
249 // if fbound is very small, then x0 is already close to the boundary and
250 // step0 points out of the allowed domain. In this case, the Newton
251 // algorithm fails, so return an error condition.
252 if (fbound < 1.e-10) {
253 debuglog("\n No damped step can be taken without violating solution component bounds.", loglevel);
254 return -3;
255 }
256
257 // ---------- Attempt damped step ----------
258
259 // damping coefficient starts at 1.0, but must be scaled by the
260 // fbound factor to ensure that the solution remains within bounds.
261 double alpha = fbound*1.0;
262 size_t m;
263 auto jac = r.getJacobian();
264 for (m = 0; m < m_maxDampIter; m++) {
265 // step the solution by the damped step size
266 // x_{k+1} = x_k + alpha_k*J(x_k)^-1 F(x_k)
267 for (size_t j = 0; j < m_n; j++) {
268 x1[j] = x0[j] + alpha*step0[j];
269 }
270
271 // compute the next undamped step that would result if x1 is accepted
272 // J(x_k)^-1 F(x_k+1)
273 step(x1, step1, r, loglevel-1);
274
275 // compute the weighted norm of step1
276 s1 = norm2(x1, step1, r);
277
278 if (loglevel > 0) {
279 double ss = r.ssnorm(x1,step1);
280 writelog("\n {:<4d} {:<9.3e} {:<9.3e} {:>6.3f} {:>6.3f} {:>6.3f} {:<5d} {:d}/{:d}",
281 m, alpha, fbound, log10(ss+SmallNumber),
282 log10(s0+SmallNumber), log10(s1+SmallNumber),
283 jac->nEvals(), jac->age(), m_maxAge);
284 }
285
286 // If the norm of s1 is less than the norm of s0, then accept this
287 // damping coefficient. Also accept it if this step would result in a
288 // converged solution. Otherwise, decrease the damping coefficient and
289 // try again.
290 if (s1 < 1.0 || s1 < s0) {
291 break;
292 }
293 alpha /= m_dampFactor;
294 }
295
296 // If a damping coefficient was found, return 1 if the solution after
297 // stepping by the damped step would represent a converged solution, and
298 // return 0 otherwise. If no damping coefficient could be found, return -2.
299 if (m < m_maxDampIter) {
300 if (s1 > 1.0) {
301 debuglog("\n Damping coefficient found (solution has not converged yet)", loglevel);
302 return 0;
303 } else {
304 debuglog("\n Damping coefficient found (solution has converged)", loglevel);
305 return 1;
306 }
307 } else {
308 debuglog("\n No damping coefficient found (max damping iterations reached)", loglevel);
309 return -2;
310 }
311}
312
313int MultiNewton::solve(double* x0, double* x1, OneDim& r, int loglevel)
314{
315 clock_t t0 = clock();
316 int status = 0;
317 bool forceNewJac = false;
318 bool write_header = true;
319 double s1=1.e30;
320
321 copy(x0, x0 + m_n, &m_x[0]);
322
323 double rdt = r.rdt();
324 int nJacReeval = 0;
325 auto jac = r.getJacobian();
326 while (true) {
327 // Check whether the Jacobian should be re-evaluated.
328 if (jac->age() > m_maxAge) {
329 if (loglevel > 1) {
330 writelog("\n Maximum Jacobian age reached ({}), updating it.", m_maxAge);
331 }
332 forceNewJac = true;
333 }
334
335 if (forceNewJac) {
336 r.evalJacobian(&m_x[0]);
337 jac->updateTransient(rdt, r.transientMask().data());
338 forceNewJac = false;
339 }
340
341 // compute the undamped Newton step
342 step(&m_x[0], &m_stp[0], r, loglevel-1);
343
344 // increment the Jacobian age
345 jac->incrementAge();
346
347 // damp the Newton step
348 status = dampStep(&m_x[0], &m_stp[0], x1, &m_stp1[0], s1, r, loglevel-1, write_header);
349 write_header = false;
350
351 // Successful step, but not converged yet. Take the damped step, and try
352 // again.
353 if (status == 0) {
354 copy(x1, x1 + m_n, m_x.begin());
355 } else if (status == 1) { // convergence
356 if (rdt == 0) {
357 jac->setAge(0); // for efficient sensitivity analysis
358 }
359 break;
360 } else if (status < 0) {
361 // If dampStep fails, first try a new Jacobian if an old one was
362 // being used. If it was a new Jacobian, then return -1 to signify
363 // failure.
364 if (jac->age() > 1) {
365 forceNewJac = true;
366 if (nJacReeval > 3) {
367 break;
368 }
369 nJacReeval++;
370 if (loglevel > 1) {
371 writelog("\n Re-evaluating Jacobian (damping coefficient not found"
372 " with this Jacobian)");
373 }
374 } else {
375 break;
376 }
377 }
378 }
379 // Close off the damped iteration table that is written by the dampedStep() method
380 if (loglevel > 1) {
381 writelog("\n {:->70}", "");
382 }
383
384 if (status < 0) { // Reset x1 to x0 if the solution failed
385 copy(m_x.begin(), m_x.end(), x1);
386 }
387 m_elapsed += (clock() - t0)/(1.0*CLOCKS_PER_SEC);
388 return status;
389}
390
391} // end namespace Cantera
Base class for exceptions thrown by Cantera classes.
Base class for one-dimensional domains.
Definition Domain1D.h:29
size_t nComponents() const
Number of components at each grid point.
Definition Domain1D.h:148
string id() const
Returns the identifying tag for this domain.
Definition Domain1D.h:471
size_t nPoints() const
Number of grid points in this domain.
Definition Domain1D.h:170
virtual string componentName(size_t n) const
Name of component n. May be overloaded.
Definition Domain1D.cpp:67
virtual size_t loc(size_t j=0) const
Location of the start of the local solution vector in the global solution vector.
Definition Domain1D.h:410
double m_dampFactor
Factor by which the damping coefficient is reduced in each iteration.
int solve(double *x0, double *x1, OneDim &r, int loglevel)
Find the solution to F(x) = 0 by damped Newton iteration.
void resize(size_t points)
Change the problem size.
vector< double > m_x
Work array holding the system state after the last successful step. Size m_n.
double norm2(const double *x, const double *step, OneDim &r) const
Compute the weighted 2-norm of step.
double m_elapsed
Elapsed CPU time spent computing the Jacobian.
vector< double > m_stp1
Work array holding the damped Newton step. Size m_n.
void step(double *x, double *step, OneDim &r, int loglevel)
Compute the undamped Newton step.
vector< double > m_stp
Work array holding the undamped Newton step or the system residual. Size m_n.
int dampStep(const double *x0, const double *step0, double *x1, double *step1, double &s1, OneDim &r, int loglevel, bool writetitle)
Performs a damped Newton step to solve the system of nonlinear equations.
size_t m_maxDampIter
Maximum number of damping iterations.
int m_maxAge
Maximum allowable Jacobian age before it is recomputed.
MultiNewton(int sz)
Constructor.
double boundStep(const double *x0, const double *step0, const OneDim &r, int loglevel)
Return the factor by which the undamped Newton step 'step0' must be multiplied in order to keep all s...
size_t m_n
number of variables
Container class for multiple-domain 1D problems.
Definition OneDim.h:28
size_t start(size_t i) const
The index of the start of domain i in the solution vector.
Definition OneDim.h:107
size_t size() const
Total solution vector length;.
Definition OneDim.h:118
void eval(size_t j, double *x, double *r, double rdt=-1.0, int count=1)
Evaluate the multi-domain residual function.
Definition OneDim.cpp:261
double ssnorm(double *x, double *r)
Steady-state max norm (infinity norm) of the residual evaluated using solution x.
Definition OneDim.cpp:339
vector< int > & transientMask()
Access the vector indicating which equations contain a transient term.
Definition OneDim.h:232
double rdt() const
Reciprocal of the time step.
Definition OneDim.h:171
size_t nDomains() const
Number of domains.
Definition OneDim.h:76
void evalJacobian(double *x0)
Evaluates the Jacobian at x0 using finite differences.
Definition OneDim.cpp:293
Domain1D & domain(size_t i) const
Return a reference to domain i.
Definition OneDim.h:81
void debuglog(const string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
Definition global.h:154
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:171
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
const double SmallNumber
smallest number to compare to zero.
Definition ct_defs.h:158
offset
Offsets of solution components in the 1D solution array.
Definition Flow1D.h:24
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...