23 Indx(
size_t nv,
size_t np) : m_nv(nv), m_np(np) {}
25 size_t operator()(
size_t m,
size_t j) {
35double bound_step(
const double* x,
const double* step, Domain1D& r,
int loglevel)
37 size_t np = r.nPoints();
38 size_t nv = r.nComponents();
41 bool wroteTitle =
false;
42 string separator = fmt::format(
"\n {:=>69}",
"");
44 for (
size_t m = 0; m < nv; m++) {
45 double above = r.upperBound(m);
46 double below = r.lowerBound(m);
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);
56 double newval = val + step[index(m,j)];
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));
65 if (loglevel > 1 && (newval > above || newval < below)) {
67 string header = fmt::format(
" {:=>10}",
"") +
68 "Undamped Newton step takes solution out of bounds" +
69 fmt::format(
"{:=>10}",
"");
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");
82 int domainLength = to_string(r.domainIndex()).length();
83 int gridLength = to_string(j).length();
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);
95 if (loglevel > 1 && wroteTitle) {
132double norm_square(
const double* x,
const double* step, Domain1D& r)
136 size_t nv = r.nComponents();
137 size_t np = r.nPoints();
139 for (
size_t n = 0; n < nv; n++) {
141 for (
size_t j = 0; j < np; j++) {
142 esum += fabs(x[nv*j + n]);
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;
148 f2max = std::max(f*f, f2max);
176 for (
size_t n = 0; n < nd; n++) {
187 for (
size_t n = 0; n < r.
size(); n++) {
191 auto jac = r.getJacobian();
195 if (jac->info() > 0) {
197 size_t row =
static_cast<size_t>(jac->info() - 1);
200 for (
size_t n = 0; n < r.
nDomains(); n++) {
203 if (row >= dom.
loc() && row < dom.
loc() + nComp * dom.
nPoints()) {
205 size_t pt =
offset / nComp;
206 size_t comp =
offset - pt * nComp;
208 "Jacobian is singular for domain {}, component {} at point {}\n"
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),
231 double* x1,
double* step1,
double& s1,
232 OneDim& r,
int loglevel,
bool writetitle)
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");
244 double s0 =
norm2(x0, step0, r);
247 double fbound =
boundStep(x0, step0, r, loglevel-1);
252 if (fbound < 1.e-10) {
253 debuglog(
"\n No damped step can be taken without violating solution component bounds.", loglevel);
261 double alpha = fbound*1.0;
263 auto jac = r.getJacobian();
267 for (
size_t j = 0; j <
m_n; j++) {
268 x1[j] = x0[j] + alpha*step0[j];
273 step(x1, step1, r, loglevel-1);
276 s1 =
norm2(x1, step1, r);
279 double ss = r.
ssnorm(x1,step1);
280 writelog(
"\n {:<4d} {:<9.3e} {:<9.3e} {:>6.3f} {:>6.3f} {:>6.3f} {:<5d} {:d}/{:d}",
283 jac->nEvals(), jac->age(),
m_maxAge);
290 if (s1 < 1.0 || s1 < s0) {
301 debuglog(
"\n Damping coefficient found (solution has not converged yet)", loglevel);
304 debuglog(
"\n Damping coefficient found (solution has converged)", loglevel);
308 debuglog(
"\n No damping coefficient found (max damping iterations reached)", loglevel);
315 clock_t t0 = clock();
317 bool forceNewJac =
false;
318 bool write_header =
true;
321 copy(x0, x0 +
m_n, &
m_x[0]);
323 double rdt = r.
rdt();
325 auto jac = r.getJacobian();
349 write_header =
false;
354 copy(x1, x1 +
m_n,
m_x.begin());
355 }
else if (status == 1) {
360 }
else if (status < 0) {
364 if (jac->age() > 1) {
366 if (nJacReeval > 3) {
371 writelog(
"\n Re-evaluating Jacobian (damping coefficient not found"
372 " with this Jacobian)");
385 copy(
m_x.begin(),
m_x.end(), x1);
387 m_elapsed += (clock() - t0)/(1.0*CLOCKS_PER_SEC);
Base class for exceptions thrown by Cantera classes.
Base class for one-dimensional domains.
size_t nComponents() const
Number of components at each grid point.
string id() const
Returns the identifying tag for this domain.
size_t nPoints() const
Number of grid points in this domain.
virtual string componentName(size_t n) const
Name of component n. May be overloaded.
virtual size_t loc(size_t j=0) const
Location of the start of the local solution vector in the global solution vector.
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.
size_t start(size_t i) const
The index of the start of domain i in the solution vector.
size_t size() const
Total solution vector length;.
void eval(size_t j, double *x, double *r, double rdt=-1.0, int count=1)
Evaluate the multi-domain residual function.
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.
double rdt() const
Reciprocal of the time step.
size_t nDomains() const
Number of domains.
void evalJacobian(double *x0)
Evaluates the Jacobian at x0 using finite differences.
Domain1D & domain(size_t i) const
Return a reference to domain i.
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 size_t npos
index returned by functions to indicate "no position"
const double SmallNumber
smallest number to compare to zero.
offset
Offsets of solution components in the 1D solution array.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...