research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[1D] Merged newton_utils.cpp with MultiNewton.cpp

Browse source
This commit is contained in:
Ray Speth 2012-12-30 00:17:02 +00:00
parent d3c0411f3e
commit 8da49ed999
2 changed files with 124 additions and 148 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

134
src/oneD/MultiNewton.cpp
View file

@ -26,17 +26,133 @@ using namespace std;
namespace Cantera
{
//----------------------------------------------------------
// function declarations
//----------------------------------------------------------
// unnamed-namespace for local helpers
namespace {
// declarations for functions in newton_utils.h
doublereal bound_step(const doublereal* x,
const doublereal* step, Domain1D& r, int loglevel=0);
class Indx
{
public:
Indx(size_t nv, size_t np) : m_nv(nv), m_np(np) {}
size_t m_nv, m_np;
size_t operator()(size_t m, size_t j) {
return j*m_nv + m;
}
};
/**
* Return a damping coefficient that keeps the solution after taking one
* Newton step between specified lower and upper bounds. This function only
* considers one domain.
*/
doublereal bound_step(const doublereal* x, const doublereal* step,
Domain1D& r, int loglevel)
{
char buf[100];
size_t np = r.nPoints();
size_t nv = r.nComponents();
Indx index(nv, np);
doublereal above, below, val, newval;
size_t m, j;
doublereal fbound = 1.0;
bool wroteTitle = false;
for (m = 0; m < nv; m++) {
above = r.upperBound(m);
below = r.lowerBound(m);
for (j = 0; j < np; j++) {
val = x[index(m,j)];
if (loglevel > 0) {
if (val > above + 1.0e-12 || val < below - 1.0e-12) {
sprintf(buf, "domain %s: %20s(%s) = %10.3e (%10.3e, %10.3e)\n",
int2str(r.domainIndex()).c_str(),
r.componentName(m).c_str(), int2str(j).c_str(),
val, below, above);
writelog(string("\nERROR: solution out of bounds.\n")+buf);
}
}
newval = val + step[index(m,j)];
if (newval > above) {
fbound = std::max(0.0, std::min(fbound,
(above - val)/(newval - val)));
} else if (newval < below) {
fbound = std::min(fbound, (val - below)/(val - newval));
}
if (loglevel > 1 && (newval > above || newval < below)) {
if (!wroteTitle) {
writelog("\nNewton step takes solution out of bounds.\n\n");
sprintf(buf," %12s %12s %4s %10s %10s %10s %10s\n",
"domain","component","pt","value","step","min","max");
wroteTitle = true;
writelog(buf);
}
sprintf(buf, " %4s %12s %4s %10.3e %10.3e %10.3e %10.3e\n",
int2str(r.domainIndex()).c_str(),
r.componentName(m).c_str(), int2str(j).c_str(),
val, step[index(m,j)], below, above);
writelog(buf);
}
}
}
return fbound;
}
/**
* This function computes the square of a weighted norm of a step
* vector for one domain.
*
* @param x Solution vector for this domain.
* @param step Newton step vector for this domain.
* @param r Object representing the domain. Used to get tolerances,
* number of components, and number of points.
*
* The return value is
* \f[
* \sum_{n,j} \left(\frac{s_{n,j}}{w_n}\right)^2
* \f]
* where the error weight for solution component \f$n\f$ is given by
* \f[
* w_n = \epsilon_{r,n} \frac{\sum_j |x_{n,j}|}{J} + \epsilon_{a,n}.
* \f]
* Here \f$\epsilon_{r,n} \f$ is the relative error tolerance for
* component n, and multiplies the average magnitude of
* solution component n in the domain. The second term,
* \f$\epsilon_{a,n}\f$, is the absolute error tolerance for component
* n.
*
*/
doublereal norm_square(const doublereal* x,
const doublereal* step, Domain1D& r);
const doublereal* step, Domain1D& r)
{
doublereal f, ewt, esum, sum = 0.0;
size_t n, j;
doublereal f2max = 0.0;
size_t nv = r.nComponents();
size_t np = r.nPoints();
for (n = 0; n < nv; n++) {
esum = 0.0;
for (j = 0; j < np; j++) {
esum += fabs(x[nv*j + n]);
}
ewt = r.rtol(n)*esum/np + r.atol(n);
for (j = 0; j < np; j++) {
f = step[nv*j + n]/ewt;
sum += f*f;
if (f*f > f2max) {
f2max = f*f;
}
}
}
return sum;
}
} // end unnamed-namespace
//-----------------------------------------------------------
// constants
@ -424,7 +540,5 @@ void MultiNewton::releaseWorkArray(doublereal* work)
{
m_workarrays.push_back(work);
}
}
// $Log: Newton.cpp,v
} // end namespace Cantera

138
src/oneD/newton_utils.cpp
View file

@ -1,138 +0,0 @@
/**
* @file newton_utils.cpp
*/
#include "cantera/base/ct_defs.h"
#include "cantera/oneD/Domain1D.h"
#include <cstdio>
using namespace std;
namespace Cantera
{
class Indx
{
public:
Indx(size_t nv, size_t np) : m_nv(nv), m_np(np) {}
size_t m_nv, m_np;
size_t operator()(size_t m, size_t j) {
return j*m_nv + m;
}
};
/**
* Return a damping coefficient that keeps the solution after taking one
* Newton step between specified lower and upper bounds. This function only
* considers one domain.
*/
doublereal bound_step(const doublereal* x, const doublereal* step,
Domain1D& r, int loglevel)
{
char buf[100];
size_t np = r.nPoints();
size_t nv = r.nComponents();
Indx index(nv, np);
doublereal above, below, val, newval;
size_t m, j;
doublereal fbound = 1.0;
bool wroteTitle = false;
for (m = 0; m < nv; m++) {
above = r.upperBound(m);
below = r.lowerBound(m);
for (j = 0; j < np; j++) {
val = x[index(m,j)];
if (loglevel > 0) {
if (val > above + 1.0e-12 || val < below - 1.0e-12) {
sprintf(buf, "domain %s: %20s(%s) = %10.3e (%10.3e, %10.3e)\n",
int2str(r.domainIndex()).c_str(),
r.componentName(m).c_str(), int2str(j).c_str(),
val, below, above);
writelog(string("\nERROR: solution out of bounds.\n")+buf);
}
}
newval = val + step[index(m,j)];
if (newval > above) {
fbound = std::max(0.0, std::min(fbound,
(above - val)/(newval - val)));
} else if (newval < below) {
fbound = std::min(fbound, (val - below)/(val - newval));
}
if (loglevel > 1 && (newval > above || newval < below)) {
if (!wroteTitle) {
writelog("\nNewton step takes solution out of bounds.\n\n");
sprintf(buf," %12s %12s %4s %10s %10s %10s %10s\n",
"domain","component","pt","value","step","min","max");
wroteTitle = true;
writelog(buf);
}
sprintf(buf, " %4s %12s %4s %10.3e %10.3e %10.3e %10.3e\n",
int2str(r.domainIndex()).c_str(),
r.componentName(m).c_str(), int2str(j).c_str(),
val, step[index(m,j)], below, above);
writelog(buf);
}
}
}
return fbound;
}
/**
* This function computes the square of a weighted norm of a step
* vector for one domain.
*
* @param x Solution vector for this domain.
* @param step Newton step vector for this domain.
* @param r Object representing the domain. Used to get tolerances,
* number of components, and number of points.
*
* The return value is
* \f[
* \sum_{n,j} \left(\frac{s_{n,j}}{w_n}\right)^2
* \f]
* where the error weight for solution component \f$n\f$ is given by
* \f[
* w_n = \epsilon_{r,n} \frac{\sum_j |x_{n,j}|}{J} + \epsilon_{a,n}.
* \f]
* Here \f$\epsilon_{r,n} \f$ is the relative error tolerance for
* component n, and multiplies the average magnitude of
* solution component n in the domain. The second term,
* \f$\epsilon_{a,n}\f$, is the absolute error tolerance for component
* n.
*
*/
doublereal norm_square(const doublereal* x,
const doublereal* step, Domain1D& r)
{
doublereal f, ewt, esum, sum = 0.0;
size_t n, j;
doublereal f2max = 0.0;
size_t nv = r.nComponents();
size_t np = r.nPoints();
for (n = 0; n < nv; n++) {
esum = 0.0;
for (j = 0; j < np; j++) {
esum += fabs(x[nv*j + n]);
}
ewt = r.rtol(n)*esum/np + r.atol(n);
for (j = 0; j < np; j++) {
f = step[nv*j + n]/ewt;
sum += f*f;
if (f*f > f2max) {
f2max = f*f;
}
}
}
return sum;
}
}