32 for (
size_t n = 0; n < r.
size(); n++) {
40 if (jac->info() > 0) {
42 size_t row =
static_cast<size_t>(jac->info() - 1);
44 "Jacobian is singular for matrix row {}:\n{}",
54 const static string separator = fmt::format(
"\n {:=>71}",
"");
56 bool wroteTitle =
false;
58 for (
size_t i = 0; i <
size(); 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",
69 double newval = val + step0[i];
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));
78 if (loglevel > 1 && (newval > above || newval < below)) {
80 string header = fmt::format(
" {:=>10}",
"") +
81 " Undamped Newton step takes solution out of bounds " +
82 fmt::format(
"{:=>10}",
"");
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");
95 writelog(
"\n {:<24s} {:>10.3e} {:>10.3e} {:>9.2e} {:>9.2e}",
96 comp_info, val, step0[i], below, above);
99 if (loglevel > 1 && wroteTitle) {
106 double* x1,
double* step1,
double& s1,
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");
122 double fbound =
boundStep(x0, step0, r, loglevel-1);
127 if (fbound < 1.e-10) {
128 debuglog(
"\n No damped step can be taken without violating solution component bounds.", loglevel);
136 double alpha = fbound*1.0;
142 for (
size_t j = 0; j <
m_n; j++) {
143 x1[j] = x0[j] + alpha*step0[j];
152 step(x1, step1, r, loglevel-1);
156 writelog(
"\n Steady-state Newton step failed:"
169 double ss = r.
ssnorm(x1,step1);
170 writelog(
"\n {:<4d} {:<9.3e} {:<9.3e} {:>6.3f} {:>6.3f} {:>6.3f} {:<5d} {:d}/{:d}",
173 jac->nEvals(), jac->age(),
m_maxAge);
180 if (s1 < 1.0 || s1 < s0) {
191 debuglog(
"\n Damping coefficient found (solution has not converged yet)", loglevel);
194 debuglog(
"\n Damping coefficient found (solution has converged)", loglevel);
198 debuglog(
"\n No damping coefficient found (max damping iterations reached)", loglevel);
205 clock_t t0 = clock();
207 bool forceNewJac =
false;
208 bool write_header =
true;
211 copy(x0, x0 +
m_n, &
m_x[0]);
213 double rdt = r.
rdt();
234 writelog(
"\n Steady Jacobian factorization failed:"
254 write_header =
false;
259 copy(x1, x1 +
m_n,
m_x.begin());
260 }
else if (status == 1) {
265 }
else if (status < 0) {
269 if (jac->age() > 1) {
271 if (nJacReeval > 3) {
276 writelog(
"\n Re-evaluating Jacobian (damping coefficient not found"
277 " with this Jacobian)");
290 copy(
m_x.begin(),
m_x.end(), x1);
292 m_elapsed += (clock() - t0)/(1.0*CLOCKS_PER_SEC);
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.
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.
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Namespace for the Cantera kernel.
const double SmallNumber
smallest number to compare to zero.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...