Cantera  3.2.0a2
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
17 : m_n(sz)
18{
19}
20
21void MultiNewton::resize(size_t sz)
22{
23 m_n = sz;
24 m_x.resize(m_n);
25 m_stp.resize(m_n);
26 m_stp1.resize(m_n);
27}
28
29void MultiNewton::step(double* x, double* step, SteadyStateSystem& r, int loglevel)
30{
31 r.eval(x, step);
32 for (size_t n = 0; n < r.size(); n++) {
33 step[n] = -step[n];
34 }
35
36 auto jac = r.linearSolver();
37 try {
38 jac->solve(r.size(), step, step);
39 } catch (CanteraError&) {
40 if (jac->info() > 0) {
41 // Positive value for "info" indicates the row where factorization failed
42 size_t row = static_cast<size_t>(jac->info() - 1);
43 throw CanteraError("MultiNewton::step",
44 "Jacobian is singular for matrix row {}:\n{}",
45 row, r.componentName(row));
46 }
47 throw;
48 }
49}
50
51double MultiNewton::boundStep(const double* x0, const double* step0, const SteadyStateSystem& r,
52 int loglevel)
53{
54 const static string separator = fmt::format("\n {:=>71}", ""); // equals sign separator
55 double fbound = 1.0;
56 bool wroteTitle = false;
57
58 for (size_t i = 0; i < size(); i++) {
59 double above = r.upperBound(i);
60 double below = r.lowerBound(i);
61
62 double val = x0[i];
63 if (loglevel > 0 && (val > above + 1.0e-12 || val < below - 1.0e-12)) {
64 writelog("\nERROR: solution component {} out of bounds.\n", i);
65 writelog("{}: value = {:10.3e} (lower = {:10.3e}, upper = {:10.3e})\n",
66 r.componentName(i), val, below, above);
67 }
68
69 double newval = val + step0[i];
70
71 if (newval > above) {
72 fbound = std::max(0.0, std::min(fbound,
73 (above - val)/(newval - val)));
74 } else if (newval < below) {
75 fbound = std::min(fbound, (val - below)/(val - newval));
76 }
77
78 if (loglevel > 1 && (newval > above || newval < below)) {
79 if (!wroteTitle){
80 string header = fmt::format(" {:=>10}", "") +
81 " Undamped Newton step takes solution out of bounds " +
82 fmt::format("{:=>10}", "");
83 writelog("\n{}", header);
84 writelog(separator);
85 const auto& [custom1, custom2] = r.componentTableHeader();
86 // Split header across 2 lines to shorten the line length
87 writelog("\n {:<24s} {:<10s} {:<10s} {:<9s} {:<9s}",
88 custom1, "", "Value", "Min", "Max");
89 writelog("\n {:<24s} {:<10s} {:<10s} {:<9s} {:<9s}",
90 custom2, "Value", "Change", "Bound", "Bound");
91 writelog(separator);
92 wroteTitle = true;
93 }
94 string comp_info = r.componentTableLabel(i);
95 writelog("\n {:<24s} {:>10.3e} {:>10.3e} {:>9.2e} {:>9.2e}",
96 comp_info, val, step0[i], below, above);
97 }
98 }
99 if (loglevel > 1 && wroteTitle) { // If a title was written, close up the table
100 writelog(separator);
101 }
102return fbound;
103}
104
105int MultiNewton::dampStep(const double* x0, const double* step0,
106 double* x1, double* step1, double& s1,
107 SteadyStateSystem& r, int loglevel, bool writetitle)
108{
109 // write header
110 if (loglevel > 0 && writetitle) {
111 writelog("\n\n {:-^70}", " Damped Newton iteration ");
112 writelog("\n {:<4s} {:<10s} {:<10s} {:<7s} {:<7s} {:<7s} {:<5s} {:<3s}\n",
113 "Iter", "F_damp", "F_bound", "log(ss)",
114 "log(s0)", "log(s1)", "N_jac", "Age");
115 writelog(" {:->70}", "");
116 }
117
118 // compute the weighted norm of the undamped step size step0
119 double s0 = r.weightedNorm(step0);
120
121 // compute the multiplier to keep all components in bounds
122 double fbound = boundStep(x0, step0, r, loglevel-1);
123
124 // if fbound is very small, then x0 is already close to the boundary and
125 // step0 points out of the allowed domain. In this case, the Newton
126 // algorithm fails, so return an error condition.
127 if (fbound < 1.e-10) {
128 debuglog("\n No damped step can be taken without violating solution component bounds.", loglevel);
129 return -3;
130 }
131
132 // ---------- Attempt damped step ----------
133
134 // damping coefficient starts at 1.0, but must be scaled by the
135 // fbound factor to ensure that the solution remains within bounds.
136 double alpha = fbound*1.0;
137 size_t m;
138 auto jac = r.linearSolver();
139 for (m = 0; m < m_maxDampIter; m++) {
140 // step the solution by the damped step size
141 // x_{k+1} = x_k + alpha_k*J(x_k)^-1 F(x_k)
142 for (size_t j = 0; j < m_n; j++) {
143 x1[j] = x0[j] + alpha*step0[j];
144 }
145
146 // compute the next undamped step that would result if x1 is accepted
147 // J(x_k)^-1 F(x_k+1)
148 try {
149 // Allow solver to continue after step failure (which can involve problems
150 // such as errors setting the thermodynamic state) by returning to time
151 // stepping mode.
152 step(x1, step1, r, loglevel-1);
153 } catch (CanteraError& err) {
154 if (r.rdt() == 0) {
155 if (loglevel > 1) {
156 writelog("\n Steady-state Newton step failed:"
157 "\n {}:\n {}", err.getMethod(), err.getMessage());
158 }
159 return -5;
160 } else {
161 throw;
162 }
163 }
164
165 // compute the weighted norm of step1
166 s1 = r.weightedNorm(step1);
167
168 if (loglevel > 0) {
169 double ss = r.ssnorm(x1,step1);
170 writelog("\n {:<4d} {:<9.3e} {:<9.3e} {:>6.3f} {:>6.3f} {:>6.3f} {:<5d} {:d}/{:d}",
171 m, alpha, fbound, log10(ss+SmallNumber),
172 log10(s0+SmallNumber), log10(s1+SmallNumber),
173 jac->nEvals(), jac->age(), m_maxAge);
174 }
175
176 // If the norm of s1 is less than the norm of s0, then accept this
177 // damping coefficient. Also accept it if this step would result in a
178 // converged solution. Otherwise, decrease the damping coefficient and
179 // try again.
180 if (s1 < 1.0 || s1 < s0) {
181 break;
182 }
183 alpha /= m_dampFactor;
184 }
185
186 // If a damping coefficient was found, return 1 if the solution after
187 // stepping by the damped step would represent a converged solution, and
188 // return 0 otherwise. If no damping coefficient could be found, return -2.
189 if (m < m_maxDampIter) {
190 if (s1 > 1.0) {
191 debuglog("\n Damping coefficient found (solution has not converged yet)", loglevel);
192 return 0;
193 } else {
194 debuglog("\n Damping coefficient found (solution has converged)", loglevel);
195 return 1;
196 }
197 } else {
198 debuglog("\n No damping coefficient found (max damping iterations reached)", loglevel);
199 return -2;
200 }
201}
202
203int MultiNewton::solve(double* x0, double* x1, SteadyStateSystem& r, int loglevel)
204{
205 clock_t t0 = clock();
206 int status = 0;
207 bool forceNewJac = false;
208 bool write_header = true;
209 double s1=1.e30;
210
211 copy(x0, x0 + m_n, &m_x[0]);
212
213 double rdt = r.rdt();
214 int nJacReeval = 0;
215 auto jac = r.linearSolver();
216 while (true) {
217 // Check whether the Jacobian should be re-evaluated.
218 if (jac->age() > m_maxAge) {
219 if (loglevel > 1) {
220 writelog("\n Maximum Jacobian age reached ({}), updating it.", m_maxAge);
221 }
222 forceNewJac = true;
223 }
224
225 if (forceNewJac) {
226 r.evalJacobian(&m_x[0]);
227 try {
228 jac->updateTransient(rdt, r.transientMask().data());
229 } catch (CanteraError& err) {
230 // Allow solver to continue after failure to factorize the steady-state
231 // Jacobian by returning to time stepping mode
232 if (rdt == 0.0) {
233 if (loglevel > 1) {
234 writelog("\n Steady Jacobian factorization failed:"
235 "\n {}:\n {}",
236 err.getMethod(), err.getMessage());
237 }
238 return -4;
239 } else {
240 throw;
241 }
242 }
243 forceNewJac = false;
244 }
245
246 // compute the undamped Newton step
247 step(&m_x[0], &m_stp[0], r, loglevel-1);
248
249 // increment the Jacobian age
250 jac->incrementAge();
251
252 // damp the Newton step
253 status = dampStep(&m_x[0], &m_stp[0], x1, &m_stp1[0], s1, r, loglevel-1, write_header);
254 write_header = false;
255
256 // Successful step, but not converged yet. Take the damped step, and try
257 // again.
258 if (status == 0) {
259 copy(x1, x1 + m_n, m_x.begin());
260 } else if (status == 1) { // convergence
261 if (rdt == 0) {
262 jac->setAge(0); // for efficient sensitivity analysis
263 }
264 break;
265 } else if (status < 0) {
266 // If dampStep fails, first try a new Jacobian if an old one was
267 // being used. If it was a new Jacobian, then return -1 to signify
268 // failure.
269 if (jac->age() > 1) {
270 forceNewJac = true;
271 if (nJacReeval > 3) {
272 break;
273 }
274 nJacReeval++;
275 if (loglevel > 1) {
276 writelog("\n Re-evaluating Jacobian (damping coefficient not found"
277 " with this Jacobian)");
278 }
279 } else {
280 break;
281 }
282 }
283 }
284 // Close off the damped iteration table that is written by the dampedStep() method
285 if (loglevel > 1) {
286 writelog("\n {:->70}", "");
287 }
288
289 if (status < 0) { // Reset x1 to x0 if the solution failed
290 copy(m_x.begin(), m_x.end(), x1);
291 }
292 m_elapsed += (clock() - t0)/(1.0*CLOCKS_PER_SEC);
293 return status;
294}
295
296} // end namespace Cantera
Base class for exceptions thrown by Cantera classes.
virtual string getMessage() const
Method overridden by derived classes to format the error message.
virtual string getMethod() const
Get the name of the method that threw the exception.
double m_dampFactor
Factor by which the damping coefficient is reduced in each iteration.
size_t size()
Get the number of variables in the system.
Definition MultiNewton.h:35
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 m_elapsed
Elapsed CPU time spent computing the Jacobian.
vector< double > m_stp1
Work array holding the damped Newton step. Size m_n.
vector< double > m_stp
Work array holding the undamped Newton step or the system residual. Size m_n.
size_t m_maxDampIter
Maximum number of damping iterations.
int dampStep(const double *x0, const double *step0, double *x1, double *step1, double &s1, SteadyStateSystem &r, int loglevel, bool writetitle)
Performs a damped Newton step to solve the system of nonlinear equations.
int m_maxAge
Maximum allowable Jacobian age before it is recomputed.
MultiNewton(int sz)
Constructor.
int solve(double *x0, double *x1, SteadyStateSystem &r, int loglevel)
Find the solution to F(x) = 0 by damped Newton iteration.
double boundStep(const double *x0, const double *step0, const SteadyStateSystem &r, int loglevel)
Return the factor by which the undamped Newton step 'step0' must be multiplied in order to keep all s...
void step(double *x, double *step, SteadyStateSystem &r, int loglevel)
Compute the undamped Newton step.
size_t m_n
number of variables
Base class for representing a system of differential-algebraic equations and solving for its steady-s...
shared_ptr< SystemJacobian > linearSolver() const
Get the the linear solver being used to hold the Jacobian matrix and solve linear systems as part of ...
size_t size() const
Total solution vector length;.
virtual double weightedNorm(const double *step) const =0
Compute the weighted norm of step.
double ssnorm(double *x, double *r)
Steady-state max norm (infinity norm) of the residual evaluated using solution x.
vector< int > & transientMask()
Access the vector indicating which equations contain a transient term.
virtual double upperBound(size_t i) const =0
Get the upper bound for global component i in the state vector.
double rdt() const
Reciprocal of the time step.
virtual string componentName(size_t i) const
Get the name of the i-th component of the state vector.
virtual string componentTableLabel(size_t i) const
Get elements of the component name, aligned with the column headings given by componentTableHeader().
virtual double lowerBound(size_t i) const =0
Get the lower bound for global component i in the state vector.
virtual void eval(double *x, double *r, double rdt=-1.0, int count=1)=0
Evaluate the residual function.
virtual pair< string, string > componentTableHeader() const
Get header lines describing the column names included in a component label.
virtual void evalJacobian(double *x0)=0
Evaluates the Jacobian at x0 using finite differences.
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 double SmallNumber
smallest number to compare to zero.
Definition ct_defs.h:158
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...