23 Indx(
size_t nv,
size_t np) : m_nv(nv), m_np(np) {}
25 size_t operator()(
size_t m,
size_t j) {
35 doublereal bound_step(
const doublereal* x,
const doublereal* step,
36 Domain1D& r,
int loglevel)
38 size_t np = r.nPoints();
39 size_t nv = r.nComponents();
41 doublereal fbound = 1.0;
42 bool wroteTitle =
false;
43 for (
size_t m = 0; m < nv; m++) {
44 double above = r.upperBound(m);
45 double below = r.lowerBound(m);
47 for (
size_t j = 0; j < np; j++) {
48 double val = x[index(m,j)];
49 if (loglevel > 0 && (val > above + 1.0e-12 || val < below - 1.0e-12)) {
50 writelog(
"\nERROR: solution out of bounds.\n");
51 writelog(
"domain {:d}: {:>20s}({:d}) = {:10.3e} ({:10.3e}, {:10.3e})\n",
52 r.domainIndex(), r.componentName(m), j, val, below, above);
55 double newval = val + step[index(m,j)];
58 fbound = std::max(0.0, std::min(fbound,
59 (above - val)/(newval - val)));
60 }
else if (newval < below) {
61 fbound = std::min(fbound, (val - below)/(val - newval));
64 if (loglevel > 1 && (newval > above || newval < below)) {
66 writelog(
"\nNewton step takes solution out of bounds.\n\n");
67 writelog(
" {:>12s} {:>12s} {:>4s} {:>10s} {:>10s} {:>10s} {:>10s}\n",
68 "domain",
"component",
"pt",
"value",
"step",
"min",
"max");
71 writelog(
" {:4d} {:>12s} {:4d} {:10.3e} {:10.3e} {:10.3e} {:10.3e}\n",
72 r.domainIndex(), r.componentName(m), j,
73 val, step[index(m,j)], below, above);
102 doublereal norm_square(
const doublereal* x,
103 const doublereal* step, Domain1D& r)
106 doublereal f2max = 0.0;
107 size_t nv = r.nComponents();
108 size_t np = r.nPoints();
110 for (
size_t n = 0; n < nv; n++) {
112 for (
size_t j = 0; j < np; j++) {
113 esum += fabs(x[nv*j + n]);
115 double ewt = r.rtol(n)*esum/np + r.atol(n);
116 for (
size_t j = 0; j < np; j++) {
117 double f = step[nv*j + n]/ewt;
119 f2max = std::max(f*f, f2max);
129 const doublereal DampFactor = sqrt(2.0);
130 const size_t NDAMP = 7;
134 MultiNewton::MultiNewton(
int sz)
150 const doublereal* step,
OneDim& r)
const 154 for (
size_t n = 0; n < nd; n++) {
166 for (
size_t n = 0; n < r.
size(); n++) {
173 int iok = jac.
info() - 1;
177 for (n = nd-1; n !=
npos; n--) {
178 if (iok >= static_cast<int>(r.
start(n))) {
183 size_t offset = iok - r.
start(n);
187 "Jacobian is singular for domain {}, component {} at point {}\n" 197 const doublereal* step0,
const OneDim& r,
int loglevel)
199 doublereal fbound = 1.0;
200 for (
size_t i = 0; i < r.
nDomains(); i++) {
201 fbound = std::min(fbound,
202 bound_step(x0 + r.
start(i), step0 + r.
start(i),
209 doublereal* x1, doublereal* step1, doublereal& s1,
213 if (loglevel > 0 && writetitle) {
214 writelog(
"\n\nDamped Newton iteration:\n");
215 writeline(
'-', 65,
false);
217 writelog(
"\n{} {:>9s} {:>9s} {:>9s} {:>9s} {:>9s} {:>5s} {:>5s}\n",
218 "m",
"F_damp",
"F_bound",
"log10(ss)",
219 "log10(s0)",
"log10(s1)",
"N_jac",
"Age");
224 doublereal s0 =
norm2(x0, step0, r);
227 doublereal fbound =
boundStep(x0, step0, r, loglevel-1);
232 if (fbound < 1.e-10) {
233 debuglog(
"\nAt limits.\n", loglevel);
240 doublereal damp = 1.0;
242 for (m = 0; m < NDAMP; m++) {
243 double ff = fbound*damp;
246 for (
size_t j = 0; j <
m_n; j++) {
247 x1[j] = ff*step0[j] + x0[j];
251 step(x1, step1, r, jac, loglevel-1);
254 s1 =
norm2(x1, step1, r);
258 doublereal ss = r.
ssnorm(x1,step1);
259 writelog(
"\n{:d} {:9.5f} {:9.5f} {:9.5f} {:9.5f} {:9.5f} {:4d} {:d}/{:d}",
269 if (s1 < 1.0 || s1 < s0) {
292 clock_t t0 = clock();
294 bool forceNewJac =
false;
297 copy(x0, x0 +
m_n, &
m_x[0]);
300 doublereal rdt = r.
rdt();
306 if (jac.
age() > m_maxAge) {
308 writelog(
"\nMaximum Jacobian age reached ({})\n", m_maxAge);
315 jac.
eval(&
m_x[0], &m_stp[0], 0.0);
316 jac.updateTransient(rdt, r.transientMask().data());
321 step(&
m_x[0], &m_stp[0], r, jac, loglevel-1);
327 m =
dampStep(&
m_x[0], &m_stp[0], x1, &m_stp1[0], s1, r, jac, loglevel-1, frst);
328 if (loglevel == 1 && m >= 0) {
330 writelog(
"\n\n {:>10s} {:>10s} {:>5s}",
331 "log10(ss)",
"log10(s1)",
"N_jac");
332 writelog(
"\n ------------------------------------");
334 doublereal ss = r.
ssnorm(&
m_x[0], &m_stp[0]);
335 writelog(
"\n {:10.4f} {:10.4f} {:d}",
336 log10(ss),log10(s1),jac.
nEvals());
343 copy(x1, x1 +
m_n,
m_x.begin());
356 if (nJacReeval > 3) {
360 debuglog(
"\nRe-evaluating Jacobian, since no damping " 361 "coefficient\ncould be found with this Jacobian.\n",
370 copy(
m_x.begin(),
m_x.end(), x1);
372 if (m > 0 && jac.
nEvals() == j0) {
375 m_elapsed += (clock() - t0)/(1.0*CLOCKS_PER_SEC);
Container class for multiple-domain 1D problems.
size_t m_n
number of variables
doublereal norm2(const doublereal *x, const doublereal *step, OneDim &r) const
Compute the weighted 2-norm of step.
doublereal ssnorm(doublereal *x, doublereal *r)
Steady-state max norm (infinity norm) of the residual evaluated using solution x. ...
void eval(size_t j, double *x, double *r, doublereal rdt=-1.0, int count=1)
Evaluate the multi-domain residual function.
Various templated functions that carry out common vector operations (see Templated Utility Functions)...
size_t size() const
Total solution vector length;.
const size_t npos
index returned by functions to indicate "no position"
void step(doublereal *x, doublereal *step, OneDim &r, MultiJac &jac, int loglevel)
Compute the undamped Newton step.
int solve(const doublereal *const b, doublereal *const x)
Solve the matrix problem Ax = b.
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
int info() const
LAPACK "info" flag after last factor/solve operation.
void incrementAge()
Increment the Jacobian age.
size_t m_n
Number of rows and columns of the matrix.
size_t nDomains() const
Number of domains.
Domain1D & domain(size_t i) const
Return a reference to domain i.
Base class for one-dimensional domains.
void setAge(int age)
Set the Jacobian age.
int solve(doublereal *x0, doublereal *x1, OneDim &r, MultiJac &jac, int loglevel)
Find the solution to F(X) = 0 by damped Newton iteration.
int age() const
Number of times 'incrementAge' has been called since the last evaluation.
int nEvals() const
Number of Jacobian evaluations.
Base class for exceptions thrown by Cantera classes.
void resize(size_t points)
Change the problem size.
size_t start(size_t i) const
The index of the start of domain i in the solution vector.
void eval(doublereal *x0, doublereal *resid0, double rdt)
Evaluate the Jacobian at x0.
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
void debuglog(const std::string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
Class MultiJac evaluates the Jacobian of a system of equations defined by a residual function supplie...
size_t nComponents() const
Number of components at each grid point.
doublereal boundStep(const doublereal *x0, const doublereal *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...
const doublereal SmallNumber
smallest number to compare to zero.
int dampStep(const doublereal *x0, const doublereal *step0, doublereal *x1, doublereal *step1, doublereal &s1, OneDim &r, MultiJac &jac, int loglevel, bool writetitle)
On entry, step0 must contain an undamped Newton step for the solution x0.
doublereal rdt() const
Reciprocal of the time step.
Namespace for the Cantera kernel.
vector_fp m_x
Work arrays of size m_n used in solve().