25ba149aab
Moved the external libraries to separate library files so that libcantera.a just contains its own namespace externals. Fixed several errors in the equilibrium program that occurred during the port. (int to size_t issues). Moved some equilibrium program headers to the include file system, so that it can link with equilibrium program. Worked on Cantera.mak. Needs more work. Fixed an issue with the Residual virtual base classes within numerics. They didn't inherit due to int to size_t migration. This caused numerous test problems to fail (issue with backwards compatibility - do we want it and how much do we want?). Added csvdiff back so that it's available for shell environment runtests.
4197 lines
161 KiB
C++
4197 lines
161 KiB
C++
/**
|
|
*
|
|
* @file NonlinearSolver.cpp
|
|
*
|
|
* Damped Newton solver for 0D and 1D problems
|
|
*/
|
|
|
|
/*
|
|
* $Date$
|
|
* $Revision$
|
|
*/
|
|
/*
|
|
* Copyright 2004 Sandia Corporation. Under the terms of Contract
|
|
* DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
|
|
* retains certain rights in this software.
|
|
* See file License.txt for licensing information.
|
|
*/
|
|
|
|
#include <limits>
|
|
|
|
#include "cantera/numerics/SquareMatrix.h"
|
|
#include "cantera/numerics/GeneralMatrix.h"
|
|
#include "cantera/numerics/NonlinearSolver.h"
|
|
#include "cantera/numerics/ctlapack.h"
|
|
|
|
#include "cantera/base/clockWC.h"
|
|
#include "cantera/base/vec_functions.h"
|
|
#include "cantera/base/mdp_allo.h"
|
|
#include "cantera/base/stringUtils.h"
|
|
|
|
#include <cfloat>
|
|
#include <ctime>
|
|
#include <vector>
|
|
#include <cstdio>
|
|
#include <cmath>
|
|
|
|
//@{
|
|
|
|
#ifndef CONSTD_DATA_PTR
|
|
#define CONSTD_DATA_PTR(x) (( const doublereal *) (&x[0]))
|
|
#endif
|
|
|
|
//@}
|
|
using namespace std;
|
|
|
|
namespace Cantera
|
|
{
|
|
|
|
|
|
//====================================================================================================================
|
|
//-----------------------------------------------------------
|
|
// Constants
|
|
//-----------------------------------------------------------
|
|
//! Dampfactor is the factor by which the damping factor is reduced by when a reduction in step length is warranted
|
|
const doublereal DampFactor = 4.0;
|
|
//! Number of damping steps that are carried out before the solution is deemed a failure
|
|
const int NDAMP = 7;
|
|
//====================================================================================================================
|
|
//! Print a line of a single repeated character string
|
|
/*!
|
|
* @param str Character string
|
|
* @param n Iteration length
|
|
*/
|
|
static void print_line(const char* str, int n)
|
|
{
|
|
for (int i = 0; i < n; i++) {
|
|
printf("%s", str);
|
|
}
|
|
printf("\n");
|
|
}
|
|
|
|
bool NonlinearSolver::s_TurnOffTiming(false);
|
|
|
|
#ifdef DEBUG_NUMJAC
|
|
bool NonlinearSolver::s_print_NumJac(true);
|
|
#else
|
|
bool NonlinearSolver::s_print_NumJac(false);
|
|
#endif
|
|
|
|
// Turn off printing of dogleg information
|
|
bool NonlinearSolver::s_print_DogLeg(false);
|
|
|
|
// Turn off solving the system twice and comparing the answer.
|
|
/*
|
|
* Turn this on if you want to compare the Hessian and Newton solve results.
|
|
*/
|
|
bool NonlinearSolver::s_doBothSolvesAndCompare(false);
|
|
|
|
// This toggle turns off the use of the Hessian when it is warranted by the condition number.
|
|
/*
|
|
* This is a debugging option.
|
|
*/
|
|
bool NonlinearSolver::s_alwaysAssumeNewtonGood(false);
|
|
|
|
//====================================================================================================================
|
|
// Default constructor
|
|
/*
|
|
* @param func Residual and jacobian evaluator function object
|
|
*/
|
|
NonlinearSolver::NonlinearSolver(ResidJacEval* func) :
|
|
m_func(func),
|
|
solnType_(NSOLN_TYPE_STEADY_STATE),
|
|
neq_(0),
|
|
m_ewt(0),
|
|
m_manualDeltaStepSet(0),
|
|
m_deltaStepMinimum(0),
|
|
m_y_n_curr(0),
|
|
m_ydot_n_curr(0),
|
|
m_y_nm1(0),
|
|
m_y_n_1(0),
|
|
m_ydot_n_1(0),
|
|
m_colScales(0),
|
|
m_rowScales(0),
|
|
m_rowWtScales(0),
|
|
m_resid(0),
|
|
m_wksp(0),
|
|
m_wksp_2(0),
|
|
m_residWts(0),
|
|
m_normResid_0(0.0),
|
|
m_normResid_Bound(0.0),
|
|
m_normResid_1(0.0),
|
|
m_normDeltaSoln_Newton(0.0),
|
|
m_normDeltaSoln_CP(0.0),
|
|
m_normResidTrial(0.0),
|
|
m_resid_scaled(false),
|
|
m_y_high_bounds(0),
|
|
m_y_low_bounds(0),
|
|
m_dampBound(1.0),
|
|
m_dampRes(1.0),
|
|
delta_t_n(-1.0),
|
|
m_nfe(0),
|
|
m_colScaling(0),
|
|
m_rowScaling(0),
|
|
m_numTotalLinearSolves(0),
|
|
m_numTotalNewtIts(0),
|
|
m_min_newt_its(0),
|
|
maxNewtIts_(100),
|
|
m_jacFormMethod(NSOLN_JAC_NUM),
|
|
m_nJacEval(0),
|
|
time_n(0.0),
|
|
m_matrixConditioning(0),
|
|
m_order(1),
|
|
rtol_(1.0E-3),
|
|
atolBase_(1.0E-10),
|
|
m_ydot_nm1(0),
|
|
atolk_(0),
|
|
userResidAtol_(0),
|
|
userResidRtol_(1.0E-3),
|
|
checkUserResidualTols_(0),
|
|
m_print_flag(0),
|
|
m_ScaleSolnNormToResNorm(0.001),
|
|
jacCopyPtr_(0),
|
|
HessianPtr_(0),
|
|
deltaX_CP_(0),
|
|
deltaX_Newton_(0),
|
|
residNorm2Cauchy_(0.0),
|
|
dogLegID_(0),
|
|
dogLegAlpha_(1.0),
|
|
RJd_norm_(0.0),
|
|
lambdaStar_(0.0),
|
|
Jd_(0),
|
|
deltaX_trust_(0),
|
|
norm_deltaX_trust_(0.0),
|
|
trustDelta_(1.0),
|
|
trustRegionInitializationMethod_(2),
|
|
trustRegionInitializationFactor_(1.0),
|
|
Nuu_(0.0),
|
|
dist_R0_(0.0),
|
|
dist_R1_(0.0),
|
|
dist_R2_(0.0),
|
|
dist_Total_(0.0),
|
|
JdJd_norm_(0.0),
|
|
normTrust_Newton_(0.0),
|
|
normTrust_CP_(0.0),
|
|
doDogLeg_(0),
|
|
doAffineSolve_(0) ,
|
|
CurrentTrustFactor_(1.0),
|
|
NextTrustFactor_(1.0),
|
|
ResidWtsReevaluated_(false),
|
|
ResidDecreaseSDExp_(0.0),
|
|
ResidDecreaseSD_(0.0),
|
|
ResidDecreaseNewtExp_(0.0),
|
|
ResidDecreaseNewt_(0.0)
|
|
{
|
|
neq_ = m_func->nEquations();
|
|
|
|
m_ewt.resize(neq_, rtol_);
|
|
m_deltaStepMinimum.resize(neq_, 0.001);
|
|
m_deltaStepMaximum.resize(neq_, 1.0E10);
|
|
m_y_n_curr.resize(neq_, 0.0);
|
|
m_ydot_n_curr.resize(neq_, 0.0);
|
|
m_y_nm1.resize(neq_, 0.0);
|
|
m_y_n_1.resize(neq_, 0.0);
|
|
m_ydot_n_1.resize(neq_, 0.0);
|
|
m_colScales.resize(neq_, 1.0);
|
|
m_rowScales.resize(neq_, 1.0);
|
|
m_rowWtScales.resize(neq_, 1.0);
|
|
m_resid.resize(neq_, 0.0);
|
|
m_wksp.resize(neq_, 0.0);
|
|
m_wksp_2.resize(neq_, 0.0);
|
|
m_residWts.resize(neq_, 0.0);
|
|
atolk_.resize(neq_, atolBase_);
|
|
deltaX_Newton_.resize(neq_, 0.0);
|
|
m_step_1.resize(neq_, 0.0);
|
|
m_y_n_1.resize(neq_, 0.0);
|
|
doublereal hb = std::numeric_limits<double>::max();
|
|
m_y_high_bounds.resize(neq_, hb);
|
|
m_y_low_bounds.resize(neq_, -hb);
|
|
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
atolk_[i] = atolBase_;
|
|
m_ewt[i] = atolk_[i];
|
|
}
|
|
|
|
|
|
// jacCopyPtr_->resize(neq_, 0.0);
|
|
deltaX_CP_.resize(neq_, 0.0);
|
|
Jd_.resize(neq_, 0.0);
|
|
deltaX_trust_.resize(neq_, 1.0);
|
|
|
|
}
|
|
//====================================================================================================================
|
|
NonlinearSolver::NonlinearSolver(const NonlinearSolver& right) :
|
|
m_func(right.m_func),
|
|
solnType_(NSOLN_TYPE_STEADY_STATE),
|
|
neq_(0),
|
|
m_ewt(0),
|
|
m_manualDeltaStepSet(0),
|
|
m_deltaStepMinimum(0),
|
|
m_y_n_curr(0),
|
|
m_ydot_n_curr(0),
|
|
m_y_nm1(0),
|
|
m_y_n_1(0),
|
|
m_ydot_n_1(0),
|
|
m_step_1(0),
|
|
m_colScales(0),
|
|
m_rowScales(0),
|
|
m_rowWtScales(0),
|
|
m_resid(0),
|
|
m_wksp(0),
|
|
m_wksp_2(0),
|
|
m_residWts(0),
|
|
m_normResid_0(0.0),
|
|
m_normResid_Bound(0.0),
|
|
m_normResid_1(0.0),
|
|
m_normDeltaSoln_Newton(0.0),
|
|
m_normDeltaSoln_CP(0.0),
|
|
m_normResidTrial(0.0),
|
|
m_resid_scaled(false),
|
|
m_y_high_bounds(0),
|
|
m_y_low_bounds(0),
|
|
m_dampBound(1.0),
|
|
m_dampRes(1.0),
|
|
delta_t_n(-1.0),
|
|
m_nfe(0),
|
|
m_colScaling(0),
|
|
m_rowScaling(0),
|
|
m_numTotalLinearSolves(0),
|
|
m_numTotalNewtIts(0),
|
|
m_min_newt_its(0),
|
|
maxNewtIts_(100),
|
|
m_jacFormMethod(NSOLN_JAC_NUM),
|
|
m_nJacEval(0),
|
|
time_n(0.0),
|
|
m_matrixConditioning(0),
|
|
m_order(1),
|
|
rtol_(1.0E-3),
|
|
atolBase_(1.0E-10),
|
|
m_ydot_nm1(0),
|
|
atolk_(0),
|
|
userResidAtol_(0),
|
|
userResidRtol_(1.0E-3),
|
|
checkUserResidualTols_(0),
|
|
m_print_flag(0),
|
|
m_ScaleSolnNormToResNorm(0.001),
|
|
jacCopyPtr_(0),
|
|
HessianPtr_(0),
|
|
deltaX_CP_(0),
|
|
deltaX_Newton_(0),
|
|
residNorm2Cauchy_(0.0),
|
|
dogLegID_(0),
|
|
dogLegAlpha_(1.0),
|
|
RJd_norm_(0.0),
|
|
lambdaStar_(0.0),
|
|
Jd_(0),
|
|
deltaX_trust_(0),
|
|
norm_deltaX_trust_(0.0),
|
|
trustDelta_(1.0),
|
|
trustRegionInitializationMethod_(2),
|
|
trustRegionInitializationFactor_(1.0),
|
|
Nuu_(0.0),
|
|
dist_R0_(0.0),
|
|
dist_R1_(0.0),
|
|
dist_R2_(0.0),
|
|
dist_Total_(0.0),
|
|
JdJd_norm_(0.0),
|
|
normTrust_Newton_(0.0),
|
|
normTrust_CP_(0.0),
|
|
doDogLeg_(0),
|
|
doAffineSolve_(0),
|
|
CurrentTrustFactor_(1.0),
|
|
NextTrustFactor_(1.0),
|
|
ResidWtsReevaluated_(false),
|
|
ResidDecreaseSDExp_(0.0),
|
|
ResidDecreaseSD_(0.0),
|
|
ResidDecreaseNewtExp_(0.0),
|
|
ResidDecreaseNewt_(0.0)
|
|
{
|
|
*this =operator=(right);
|
|
}
|
|
|
|
//====================================================================================================================
|
|
NonlinearSolver::~NonlinearSolver()
|
|
{
|
|
if (jacCopyPtr_) {
|
|
delete jacCopyPtr_;
|
|
}
|
|
if (HessianPtr_) {
|
|
delete HessianPtr_;
|
|
}
|
|
}
|
|
//====================================================================================================================
|
|
NonlinearSolver& NonlinearSolver::operator=(const NonlinearSolver& right)
|
|
{
|
|
if (this == &right) {
|
|
return *this;
|
|
}
|
|
// rely on the ResidJacEval duplMyselfAsresidJacEval() function to
|
|
// create a deep copy
|
|
m_func = right.m_func->duplMyselfAsResidJacEval();
|
|
|
|
solnType_ = right.solnType_;
|
|
neq_ = right.neq_;
|
|
m_ewt = right.m_ewt;
|
|
m_manualDeltaStepSet = right.m_manualDeltaStepSet;
|
|
m_deltaStepMinimum = right.m_deltaStepMinimum;
|
|
m_y_n_curr = right.m_y_n_curr;
|
|
m_ydot_n_curr = right.m_ydot_n_curr;
|
|
m_y_nm1 = right.m_y_nm1;
|
|
m_y_n_1 = right.m_y_n_1;
|
|
m_ydot_n_1 = right.m_ydot_n_1;
|
|
m_step_1 = right.m_step_1;
|
|
m_colScales = right.m_colScales;
|
|
m_rowScales = right.m_rowScales;
|
|
m_rowWtScales = right.m_rowWtScales;
|
|
m_resid = right.m_resid;
|
|
m_wksp = right.m_wksp;
|
|
m_wksp_2 = right.m_wksp_2;
|
|
m_residWts = right.m_residWts;
|
|
m_normResid_0 = right.m_normResid_0;
|
|
m_normResid_Bound = right.m_normResid_Bound;
|
|
m_normResid_1 = right.m_normResid_1;
|
|
m_normDeltaSoln_Newton = right.m_normDeltaSoln_Newton;
|
|
m_normDeltaSoln_CP = right.m_normDeltaSoln_CP;
|
|
m_normResidTrial = right.m_normResidTrial;
|
|
m_resid_scaled = right.m_resid_scaled;
|
|
m_y_high_bounds = right.m_y_high_bounds;
|
|
m_y_low_bounds = right.m_y_low_bounds;
|
|
m_dampBound = right.m_dampBound;
|
|
m_dampRes = right.m_dampRes;
|
|
delta_t_n = right.delta_t_n;
|
|
m_nfe = right.m_nfe;
|
|
m_colScaling = right.m_colScaling;
|
|
m_rowScaling = right.m_rowScaling;
|
|
m_numTotalLinearSolves = right.m_numTotalLinearSolves;
|
|
m_numTotalNewtIts = right.m_numTotalNewtIts;
|
|
m_min_newt_its = right.m_min_newt_its;
|
|
maxNewtIts_ = right.maxNewtIts_;
|
|
m_jacFormMethod = right.m_jacFormMethod;
|
|
m_nJacEval = right.m_nJacEval;
|
|
time_n = right.time_n;
|
|
m_matrixConditioning = right.m_matrixConditioning;
|
|
m_order = right.m_order;
|
|
rtol_ = right.rtol_;
|
|
atolBase_ = right.atolBase_;
|
|
atolk_ = right.atolk_;
|
|
userResidAtol_ = right.userResidAtol_;
|
|
userResidRtol_ = right.userResidRtol_;
|
|
checkUserResidualTols_ = right.checkUserResidualTols_;
|
|
m_print_flag = right.m_print_flag;
|
|
m_ScaleSolnNormToResNorm = right.m_ScaleSolnNormToResNorm;
|
|
|
|
if (jacCopyPtr_) {
|
|
delete(jacCopyPtr_);
|
|
}
|
|
jacCopyPtr_ = (right.jacCopyPtr_)->duplMyselfAsGeneralMatrix();
|
|
if (HessianPtr_) {
|
|
delete(HessianPtr_);
|
|
}
|
|
HessianPtr_ = (right.HessianPtr_)->duplMyselfAsGeneralMatrix();
|
|
|
|
deltaX_CP_ = right.deltaX_CP_;
|
|
deltaX_Newton_ = right.deltaX_Newton_;
|
|
residNorm2Cauchy_ = right.residNorm2Cauchy_;
|
|
dogLegID_ = right.dogLegID_;
|
|
dogLegAlpha_ = right.dogLegAlpha_;
|
|
RJd_norm_ = right.RJd_norm_;
|
|
lambdaStar_ = right.lambdaStar_;
|
|
Jd_ = right.Jd_;
|
|
deltaX_trust_ = right.deltaX_trust_;
|
|
norm_deltaX_trust_ = right.norm_deltaX_trust_;
|
|
trustDelta_ = right.trustDelta_;
|
|
trustRegionInitializationMethod_ = right.trustRegionInitializationMethod_;
|
|
trustRegionInitializationFactor_ = right.trustRegionInitializationFactor_;
|
|
Nuu_ = right.Nuu_;
|
|
dist_R0_ = right.dist_R0_;
|
|
dist_R1_ = right.dist_R1_;
|
|
dist_R2_ = right.dist_R2_;
|
|
dist_Total_ = right.dist_Total_;
|
|
JdJd_norm_ = right.JdJd_norm_;
|
|
normTrust_Newton_ = right.normTrust_Newton_;
|
|
normTrust_CP_ = right.normTrust_CP_;
|
|
doDogLeg_ = right.doDogLeg_;
|
|
doAffineSolve_ = right.doAffineSolve_;
|
|
CurrentTrustFactor_ = right.CurrentTrustFactor_;
|
|
NextTrustFactor_ = right.NextTrustFactor_;
|
|
|
|
ResidWtsReevaluated_ = right.ResidWtsReevaluated_;
|
|
ResidDecreaseSDExp_ = right.ResidDecreaseSDExp_;
|
|
ResidDecreaseSD_ = right.ResidDecreaseSD_;
|
|
ResidDecreaseNewtExp_ = right.ResidDecreaseNewtExp_;
|
|
ResidDecreaseNewt_ = right.ResidDecreaseNewt_;
|
|
|
|
return *this;
|
|
}
|
|
//====================================================================================================================
|
|
// Create solution weights for convergence criteria
|
|
/*
|
|
* We create soln weights from the following formula
|
|
*
|
|
* wt[i] = rtol * abs(y[i]) + atol[i]
|
|
*
|
|
* The program always assumes that atol is specific
|
|
* to the solution component
|
|
*
|
|
* @param y vector of the current solution values
|
|
*/
|
|
void NonlinearSolver::createSolnWeights(const doublereal* const y)
|
|
{
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
m_ewt[i] = rtol_ * fabs(y[i]) + atolk_[i];
|
|
}
|
|
}
|
|
//====================================================================================================================
|
|
// set bounds constraints for all variables in the problem
|
|
/*
|
|
*
|
|
* @param y_low_bounds Vector of lower bounds
|
|
* @param y_high_bounds Vector of high bounds
|
|
*/
|
|
void NonlinearSolver::setBoundsConstraints(const doublereal* const y_low_bounds,
|
|
const doublereal* const y_high_bounds)
|
|
{
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
m_y_low_bounds[i] = y_low_bounds[i];
|
|
m_y_high_bounds[i] = y_high_bounds[i];
|
|
}
|
|
}
|
|
//====================================================================================================================
|
|
void NonlinearSolver::setSolverScheme(int doDogLeg, int doAffineSolve)
|
|
{
|
|
doDogLeg_ = doDogLeg;
|
|
doAffineSolve_ = doAffineSolve;
|
|
}
|
|
//====================================================================================================================
|
|
std::vector<doublereal> & NonlinearSolver::lowBoundsConstraintVector()
|
|
{
|
|
return m_y_low_bounds;
|
|
}
|
|
//====================================================================================================================
|
|
std::vector<doublereal> & NonlinearSolver::highBoundsConstraintVector()
|
|
{
|
|
return m_y_high_bounds;
|
|
}
|
|
//====================================================================================================================
|
|
// L2 norm of the delta of the solution vector
|
|
/*
|
|
* calculate the norm of the solution vector. This will
|
|
* involve the column scaling of the matrix
|
|
*
|
|
* The third argument has a default of false. However,
|
|
* if true, then a table of the largest values is printed
|
|
* out to standard output.
|
|
*
|
|
* @param delta_y Vector to take the norm of
|
|
* @param title Optional title to be printed out
|
|
* @param printLargest int indicating how many specific lines should be printed out
|
|
* @param dampFactor Current value of the damping factor. Defaults to 1.
|
|
* only used for printout out a table.
|
|
*/
|
|
doublereal NonlinearSolver::solnErrorNorm(const doublereal* const delta_y, const char* title, int printLargest,
|
|
const doublereal dampFactor) const
|
|
{
|
|
doublereal sum_norm = 0.0, error;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
error = delta_y[i] / m_ewt[i];
|
|
sum_norm += (error * error);
|
|
}
|
|
sum_norm = sqrt(sum_norm / neq_);
|
|
if (printLargest) {
|
|
if ((printLargest == 1) || (m_print_flag >= 4 && m_print_flag <= 5)) {
|
|
|
|
printf("\t\t solnErrorNorm(): ");
|
|
if (title) {
|
|
printf("%s", title);
|
|
} else {
|
|
printf(" Delta soln norm ");
|
|
}
|
|
printf(" = %-11.4E\n", sum_norm);
|
|
} else if (m_print_flag >= 6) {
|
|
const int num_entries = printLargest;
|
|
printf("\t\t ");
|
|
print_line("-", 90);
|
|
printf("\t\t solnErrorNorm(): ");
|
|
if (title) {
|
|
printf("%s", title);
|
|
} else {
|
|
printf(" Delta soln norm ");
|
|
}
|
|
printf(" = %-11.4E\n", sum_norm);
|
|
|
|
doublereal dmax1, normContrib;
|
|
int j;
|
|
std::vector<size_t> imax(num_entries, npos);
|
|
printf("\t\t Printout of Largest Contributors: (damp = %g)\n", dampFactor);
|
|
printf("\t\t I weightdeltaY/sqtN| deltaY "
|
|
"ysolnOld ysolnNew Soln_Weights\n");
|
|
printf("\t\t ");
|
|
print_line("-", 88);
|
|
|
|
for (int jnum = 0; jnum < num_entries; jnum++) {
|
|
dmax1 = -1.0;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
bool used = false;
|
|
for (j = 0; j < jnum; j++) {
|
|
if (imax[j] == i) {
|
|
used = true;
|
|
}
|
|
}
|
|
if (!used) {
|
|
error = delta_y[i] / m_ewt[i];
|
|
normContrib = sqrt(error * error);
|
|
if (normContrib > dmax1) {
|
|
imax[jnum] = i;
|
|
dmax1 = normContrib;
|
|
}
|
|
}
|
|
}
|
|
size_t i = imax[jnum];
|
|
if (i != npos) {
|
|
error = delta_y[i] / m_ewt[i];
|
|
normContrib = sqrt(error * error);
|
|
printf("\t\t %4s %12.4e | %12.4e %12.4e %12.4e %12.4e\n",
|
|
int2str(i).c_str(), normContrib/sqrt((double)neq_), delta_y[i],
|
|
m_y_n_curr[i], m_y_n_curr[i] + dampFactor * delta_y[i], m_ewt[i]);
|
|
|
|
}
|
|
}
|
|
printf("\t\t ");
|
|
print_line("-", 90);
|
|
}
|
|
}
|
|
return sum_norm;
|
|
}
|
|
//====================================================================================================================
|
|
/*
|
|
* L2 Norm of the residual
|
|
*
|
|
* The second argument has a default of false. However,
|
|
* if true, then a table of the largest values is printed
|
|
* out to standard output.
|
|
*/
|
|
doublereal NonlinearSolver::residErrorNorm(const doublereal* const resid, const char* title, const int printLargest,
|
|
const doublereal* const y) const
|
|
{
|
|
doublereal sum_norm = 0.0, error;
|
|
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
#ifdef DEBUG_HKM
|
|
mdp::checkFinite(resid[i]);
|
|
#endif
|
|
error = resid[i] / m_residWts[i];
|
|
#ifdef DEBUG_HKM
|
|
mdp::checkFinite(error);
|
|
#endif
|
|
sum_norm += (error * error);
|
|
}
|
|
sum_norm = sqrt(sum_norm / neq_);
|
|
#ifdef DEBUG_HKM
|
|
mdp::checkFinite(sum_norm);
|
|
#endif
|
|
if (printLargest) {
|
|
const int num_entries = printLargest;
|
|
doublereal dmax1, normContrib;
|
|
int j;
|
|
std::vector<size_t> imax(num_entries, npos);
|
|
|
|
if (m_print_flag >= 4 && m_print_flag <= 5) {
|
|
printf("\t\t residErrorNorm():");
|
|
if (title) {
|
|
printf(" %s ", title);
|
|
} else {
|
|
printf(" residual L2 norm ");
|
|
}
|
|
printf("= %12.4E\n", sum_norm);
|
|
}
|
|
if (m_print_flag >= 6) {
|
|
printf("\t\t ");
|
|
print_line("-", 90);
|
|
printf("\t\t residErrorNorm(): ");
|
|
if (title) {
|
|
printf(" %s ", title);
|
|
} else {
|
|
printf(" residual L2 norm ");
|
|
}
|
|
printf("= %12.4E\n", sum_norm);
|
|
printf("\t\t Printout of Largest Contributors to norm:\n");
|
|
printf("\t\t I |Resid/ResWt| UnsclRes ResWt | y_curr\n");
|
|
printf("\t\t ");
|
|
print_line("-", 88);
|
|
for (int jnum = 0; jnum < num_entries; jnum++) {
|
|
dmax1 = -1.0;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
bool used = false;
|
|
for (j = 0; j < jnum; j++) {
|
|
if (imax[j] == i) {
|
|
used = true;
|
|
}
|
|
}
|
|
if (!used) {
|
|
error = resid[i] / m_residWts[i];
|
|
normContrib = sqrt(error * error);
|
|
if (normContrib > dmax1) {
|
|
imax[jnum] = i;
|
|
dmax1 = normContrib;
|
|
}
|
|
}
|
|
}
|
|
size_t i = imax[jnum];
|
|
if (i != npos) {
|
|
error = resid[i] / m_residWts[i];
|
|
normContrib = sqrt(error * error);
|
|
printf("\t\t %4s %12.4e %12.4e %12.4e | %12.4e\n",
|
|
int2str(i).c_str(), normContrib, resid[i], m_residWts[i], y[i]);
|
|
}
|
|
}
|
|
|
|
printf("\t\t ");
|
|
print_line("-", 90);
|
|
}
|
|
}
|
|
return sum_norm;
|
|
}
|
|
//====================================================================================================================
|
|
// Set the column scaling that are used for the inversion of the matrix
|
|
/*
|
|
* There are three ways to do this.
|
|
*
|
|
* The first method is to set the bool useColScaling to true, leaving the scaling factors unset.
|
|
* Then, the column scales will be set to the solution error weighting factors. This has the
|
|
* effect of ensuring that all delta variables will have the same order of magnitude at convergence
|
|
* end.
|
|
*
|
|
* The second way is the explicity set the column factors in the second argument of this function call.
|
|
*
|
|
* The final way to input the scales is to override the ResidJacEval member function call,
|
|
*
|
|
* calcSolnScales(double time_n, const double *m_y_n_curr, const double *m_y_nm1, double *m_colScales)
|
|
*
|
|
* Overriding this function call will trump all other ways to specify the column scaling factors.
|
|
*
|
|
* @param useColScaling Turn this on if you want to use column scaling in the calculations
|
|
* @param scaleFactors A vector of doubles that specifies the column factors.
|
|
*/
|
|
void NonlinearSolver::setColumnScaling(bool useColScaling, const double* const scaleFactors)
|
|
{
|
|
if (useColScaling) {
|
|
if (scaleFactors) {
|
|
m_colScaling = 2;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
m_colScales[i] = scaleFactors[i];
|
|
if (m_colScales[i] <= 1.0E-200) {
|
|
throw CanteraError("NonlinearSolver::setColumnScaling() ERROR", "Bad column scale factor");
|
|
}
|
|
}
|
|
} else {
|
|
m_colScaling = 1;
|
|
}
|
|
} else {
|
|
m_colScaling = 0;
|
|
}
|
|
}
|
|
//====================================================================================================================
|
|
// Set the rowscaling that are used for the inversion of the matrix
|
|
/*
|
|
* Row scaling is set here. Right now the row scaling is set internally in the code.
|
|
*
|
|
* @param useRowScaling Turn row scaling on or off.
|
|
*/
|
|
void NonlinearSolver::setRowScaling(bool useRowScaling)
|
|
{
|
|
m_rowScaling = useRowScaling;
|
|
}
|
|
//====================================================================================================================
|
|
/*
|
|
* calcColumnScales():
|
|
*
|
|
* Set the column scaling vector at the current time
|
|
*/
|
|
void NonlinearSolver::calcColumnScales()
|
|
{
|
|
if (m_colScaling == 1) {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
m_colScales[i] = m_ewt[i];
|
|
}
|
|
} else {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
m_colScales[i] = 1.0;
|
|
}
|
|
}
|
|
if (m_colScaling) {
|
|
m_func->calcSolnScales(time_n, DATA_PTR(m_y_n_curr), DATA_PTR(m_y_nm1), DATA_PTR(m_colScales));
|
|
}
|
|
}
|
|
//====================================================================================================================
|
|
// Compute the current residual
|
|
/*
|
|
* @param time_curr Value of the time
|
|
* @param typeCalc Type of the calculation
|
|
* @param y_curr Current value of the solution vector
|
|
* @param ydot_curr Current value of the time derivative of the solution vector
|
|
*
|
|
* @return Returns a flag to indicate that operation is successful.
|
|
* 1 Means a successful operation
|
|
* -0 or neg value Means an unsuccessful operation
|
|
*/
|
|
int NonlinearSolver::doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal* const y_curr,
|
|
const doublereal* const ydot_curr, const ResidEval_Type_Enum evalType) const
|
|
{
|
|
int retn = m_func->evalResidNJ(time_curr, delta_t_n, y_curr, ydot_curr, DATA_PTR(m_resid), evalType);
|
|
m_nfe++;
|
|
m_resid_scaled = false;
|
|
return retn;
|
|
}
|
|
//====================================================================================================================
|
|
// Scale the matrix
|
|
/*
|
|
* @param jac Jacobian
|
|
* @param y_comm Current value of the solution vector
|
|
* @param ydot_comm Current value of the time derivative of the solution vector
|
|
* @param time_curr current value of the time
|
|
*/
|
|
void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm, doublereal* const ydot_comm,
|
|
doublereal time_curr, int num_newt_its)
|
|
{
|
|
size_t irow, jcol;
|
|
size_t ku, kl;
|
|
size_t ivec[2];
|
|
jac.nRowsAndStruct(ivec);
|
|
double* colP_j;
|
|
|
|
/*
|
|
* Column scaling -> We scale the columns of the Jacobian
|
|
* by the nominal important change in the solution vector
|
|
*/
|
|
if (m_colScaling) {
|
|
if (!jac.factored()) {
|
|
if (jac.matrixType_ == 0) {
|
|
/*
|
|
* Go get new scales -> Took this out of this inner loop.
|
|
* Needs to be done at a larger scale.
|
|
*/
|
|
// setColumnScales();
|
|
|
|
/*
|
|
* Scale the new Jacobian
|
|
*/
|
|
doublereal* jptr = &(*(jac.begin()));
|
|
for (jcol = 0; jcol < neq_; jcol++) {
|
|
for (irow = 0; irow < neq_; irow++) {
|
|
*jptr *= m_colScales[jcol];
|
|
jptr++;
|
|
}
|
|
}
|
|
} else if (jac.matrixType_ == 1) {
|
|
kl = ivec[0];
|
|
ku = ivec[1];
|
|
for (jcol = 0; jcol < neq_; jcol++) {
|
|
colP_j = (doublereal*) jac.ptrColumn(jcol);
|
|
for (irow = jcol - ku; irow <= jcol + kl; irow++) {
|
|
if (irow < neq_) {
|
|
colP_j[kl + ku + irow - jcol] *= m_colScales[jcol];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* row sum scaling -> Note, this is an unequivocal success
|
|
* at keeping the small numbers well balanced and nonnegative.
|
|
*/
|
|
if (! jac.factored()) {
|
|
/*
|
|
* Ok, this is ugly. jac.begin() returns an vector<double> iterator
|
|
* to the first data location.
|
|
* Then &(*()) reverts it to a doublereal *.
|
|
*/
|
|
doublereal* jptr = &(*(jac.begin()));
|
|
for (irow = 0; irow < neq_; irow++) {
|
|
m_rowScales[irow] = 0.0;
|
|
m_rowWtScales[irow] = 0.0;
|
|
}
|
|
if (jac.matrixType_ == 0) {
|
|
for (jcol = 0; jcol < neq_; jcol++) {
|
|
for (irow = 0; irow < neq_; irow++) {
|
|
if (m_rowScaling) {
|
|
m_rowScales[irow] += fabs(*jptr);
|
|
}
|
|
if (m_colScaling) {
|
|
// This is needed in order to mitigate the change in J_ij carried out just above this loop.
|
|
// Alternatively, we could move this loop up to the top
|
|
m_rowWtScales[irow] += fabs(*jptr) * m_ewt[jcol] / m_colScales[jcol];
|
|
} else {
|
|
m_rowWtScales[irow] += fabs(*jptr) * m_ewt[jcol];
|
|
}
|
|
jptr++;
|
|
}
|
|
}
|
|
} else if (jac.matrixType_ == 1) {
|
|
kl = ivec[0];
|
|
ku = ivec[1];
|
|
for (jcol = 0; jcol < neq_; jcol++) {
|
|
colP_j = (doublereal*) jac.ptrColumn(jcol);
|
|
for (irow = jcol - ku; irow <= jcol + kl; irow++) {
|
|
if (irow < neq_) {
|
|
double vv = fabs(colP_j[kl + ku + irow - jcol]);
|
|
if (m_rowScaling) {
|
|
m_rowScales[irow] += vv;
|
|
}
|
|
if (m_colScaling) {
|
|
// This is needed in order to mitigate the change in J_ij carried out just above this loop.
|
|
// Alternatively, we could move this loop up to the top
|
|
m_rowWtScales[irow] += vv * m_ewt[jcol] / m_colScales[jcol];
|
|
} else {
|
|
m_rowWtScales[irow] += vv * m_ewt[jcol];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if (m_rowScaling) {
|
|
for (irow = 0; irow < neq_; irow++) {
|
|
m_rowScales[irow] = 1.0/m_rowScales[irow];
|
|
}
|
|
} else {
|
|
for (irow = 0; irow < neq_; irow++) {
|
|
m_rowScales[irow] = 1.0;
|
|
}
|
|
}
|
|
// What we have defined is a maximum value that the residual can be and still pass.
|
|
// This isn't sufficient.
|
|
|
|
if (m_rowScaling) {
|
|
if (jac.matrixType_ == 0) {
|
|
jptr = &(*(jac.begin()));
|
|
for (jcol = 0; jcol < neq_; jcol++) {
|
|
for (irow = 0; irow < neq_; irow++) {
|
|
*jptr *= m_rowScales[irow];
|
|
jptr++;
|
|
}
|
|
}
|
|
} else if (jac.matrixType_ == 1) {
|
|
kl = ivec[0];
|
|
ku = ivec[1];
|
|
for (jcol = 0; jcol < neq_; jcol++) {
|
|
colP_j = (doublereal*) jac.ptrColumn(jcol);
|
|
for (irow = jcol - ku; irow <= jcol + kl; irow++) {
|
|
if (irow < neq_) {
|
|
colP_j[kl + ku + irow - jcol] *= m_rowScales[irow];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if (num_newt_its % 5 == 1) {
|
|
computeResidWts();
|
|
}
|
|
|
|
}
|
|
}
|
|
//====================================================================================================================
|
|
// Calculate the scaling factor for translating residual norms into solution norms.
|
|
/*
|
|
* This routine calls computeResidWts() a couple of times in the calculation of m_ScaleSolnNormToResNorm.
|
|
* A more sophisticated routine may do more with signs to get a better value. Perhaps, a series of calculations
|
|
* with different signs attached may be in order. Then, m_ScaleSolnNormToResNorm would be calculated
|
|
* as the minimum of a series of calculations.
|
|
*/
|
|
void NonlinearSolver::calcSolnToResNormVector()
|
|
{
|
|
if (! jacCopyPtr_->factored()) {
|
|
|
|
if (checkUserResidualTols_ != 1) {
|
|
doublereal sum = 0.0;
|
|
for (size_t irow = 0; irow < neq_; irow++) {
|
|
m_residWts[irow] = m_rowWtScales[irow] / neq_;
|
|
sum += m_residWts[irow];
|
|
}
|
|
sum /= neq_;
|
|
for (size_t irow = 0; irow < neq_; irow++) {
|
|
m_residWts[irow] = (m_residWts[irow] + atolBase_ * atolBase_ * sum);
|
|
}
|
|
if (checkUserResidualTols_ == 2) {
|
|
for (size_t irow = 0; irow < neq_; irow++) {
|
|
m_residWts[irow] = std::min(m_residWts[irow], userResidAtol_[irow] + userResidRtol_ * m_rowWtScales[irow] / neq_);
|
|
}
|
|
}
|
|
} else {
|
|
for (size_t irow = 0; irow < neq_; irow++) {
|
|
m_residWts[irow] = userResidAtol_[irow] + userResidRtol_ * m_rowWtScales[irow] / neq_;
|
|
}
|
|
}
|
|
|
|
|
|
for (size_t irow = 0; irow < neq_; irow++) {
|
|
m_wksp[irow] = 0.0;
|
|
}
|
|
doublereal* jptr = &(jacCopyPtr_->operator()(0,0));
|
|
for (size_t jcol = 0; jcol < neq_; jcol++) {
|
|
for (size_t irow = 0; irow < neq_; irow++) {
|
|
m_wksp[irow] += (*jptr) * m_ewt[jcol];
|
|
jptr++;
|
|
}
|
|
}
|
|
doublereal resNormOld = 0.0;
|
|
doublereal error;
|
|
|
|
for (size_t irow = 0; irow < neq_; irow++) {
|
|
error = m_wksp[irow] / m_residWts[irow];
|
|
resNormOld += error * error;
|
|
}
|
|
resNormOld = sqrt(resNormOld / neq_);
|
|
|
|
if (resNormOld > 0.0) {
|
|
m_ScaleSolnNormToResNorm = resNormOld;
|
|
}
|
|
if (m_ScaleSolnNormToResNorm < 1.0E-8) {
|
|
m_ScaleSolnNormToResNorm = 1.0E-8;
|
|
}
|
|
|
|
// Recalculate the residual weights now that we know the value of m_ScaleSolnNormToResNorm
|
|
computeResidWts();
|
|
} else {
|
|
throw CanteraError("NonlinearSolver::calcSolnToResNormVector()" , "Logic error");
|
|
}
|
|
}
|
|
//====================================================================================================================
|
|
// Compute the undamped Newton step based on the current jacobian and an input rhs
|
|
/*
|
|
* Compute the undamped Newton step. The residual function is
|
|
* evaluated at the current time, t_n, at the current values of the
|
|
* solution vector, m_y_n_curr, and the solution time derivative, m_ydot_n.
|
|
* The Jacobian is not recomputed.
|
|
*
|
|
* A factored jacobian is reused, if available. If a factored jacobian
|
|
* is not available, then the jacobian is factored. Before factoring,
|
|
* the jacobian is row and column-scaled. Column scaling is not
|
|
* recomputed. The row scales are recomputed here, after column
|
|
* scaling has been implemented.
|
|
*/
|
|
int NonlinearSolver::doNewtonSolve(const doublereal time_curr, const doublereal* const y_curr,
|
|
const doublereal* const ydot_curr, doublereal* const delta_y,
|
|
GeneralMatrix& jac)
|
|
{
|
|
// multiply the residual by -1
|
|
if (m_rowScaling && !m_resid_scaled) {
|
|
for (size_t n = 0; n < neq_; n++) {
|
|
delta_y[n] = -m_rowScales[n] * m_resid[n];
|
|
}
|
|
m_resid_scaled = true;
|
|
} else {
|
|
for (size_t n = 0; n < neq_; n++) {
|
|
delta_y[n] = -m_resid[n];
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Solve the system -> This also involves inverting the
|
|
* matrix
|
|
*/
|
|
int info = jac.solve(DATA_PTR(delta_y));
|
|
|
|
/*
|
|
* reverse the column scaling if there was any.
|
|
*/
|
|
if (m_colScaling) {
|
|
for (size_t irow = 0; irow < neq_; irow++) {
|
|
delta_y[irow] = delta_y[irow] * m_colScales[irow];
|
|
}
|
|
}
|
|
|
|
#ifdef DEBUG_JAC
|
|
if (printJacContributions) {
|
|
for (size_t iNum = 0; iNum < numRows; iNum++) {
|
|
if (iNum > 0) {
|
|
focusRow++;
|
|
}
|
|
doublereal dsum = 0.0;
|
|
vector_fp& Jdata = jacBack.data();
|
|
doublereal dRow = Jdata[neq_ * focusRow + focusRow];
|
|
printf("\n Details on delta_Y for row %d \n", focusRow);
|
|
printf(" Value before = %15.5e, delta = %15.5e,"
|
|
"value after = %15.5e\n", y_curr[focusRow],
|
|
delta_y[focusRow], y_curr[focusRow] + delta_y[focusRow]);
|
|
if (!freshJac) {
|
|
printf(" Old Jacobian\n");
|
|
}
|
|
printf(" col delta_y aij "
|
|
"contrib \n");
|
|
printf("--------------------------------------------------"
|
|
"---------------------------------------------\n");
|
|
printf(" Res(%d) %15.5e %15.5e %15.5e (Res = %g)\n",
|
|
focusRow, delta_y[focusRow],
|
|
dRow, RRow[iNum] / dRow, RRow[iNum]);
|
|
dsum += RRow[iNum] / dRow;
|
|
for (size_t ii = 0; ii < neq_; ii++) {
|
|
if (ii != focusRow) {
|
|
doublereal aij = Jdata[neq_ * ii + focusRow];
|
|
doublereal contrib = aij * delta_y[ii] * (-1.0) / dRow;
|
|
dsum += contrib;
|
|
if (fabs(contrib) > Pcutoff) {
|
|
printf("%6d %15.5e %15.5e %15.5e\n", ii,
|
|
delta_y[ii] , aij, contrib);
|
|
}
|
|
}
|
|
}
|
|
printf("--------------------------------------------------"
|
|
"---------------------------------------------\n");
|
|
printf(" %15.5e %15.5e\n",
|
|
delta_y[focusRow], dsum);
|
|
}
|
|
}
|
|
|
|
#endif
|
|
|
|
m_numTotalLinearSolves++;
|
|
m_numLocalLinearSolves++;
|
|
return info;
|
|
}
|
|
//====================================================================================================================
|
|
// Compute the newton step, either by direct newton's or by solving a close problem that is represented
|
|
// by a Hessian (
|
|
/*
|
|
* This is algorith A.6.5.1 in Dennis / Schnabel
|
|
*
|
|
* Compute the QR decomposition
|
|
*
|
|
* Notes on banded Hessian solve:
|
|
* The matrix for jT j has a larger band width. Both the top and bottom band widths
|
|
* are doubled, going from KU to KU+KL and KL to KU+KL in size. This is not an impossible increase in cost, but
|
|
* has to be considered.
|
|
*/
|
|
int NonlinearSolver::doAffineNewtonSolve(const doublereal* const y_curr, const doublereal* const ydot_curr,
|
|
doublereal* const delta_y, GeneralMatrix& jac)
|
|
{
|
|
bool newtonGood = true;
|
|
doublereal* delyNewton = 0;
|
|
// We can default to QR here ( or not )
|
|
jac.useFactorAlgorithm(1);
|
|
int useQR = jac.factorAlgorithm();
|
|
// multiply the residual by -1
|
|
// Scale the residual if there is row scaling. Note, the matrix has already been scaled
|
|
if (m_rowScaling && !m_resid_scaled) {
|
|
for (size_t n = 0; n < neq_; n++) {
|
|
delta_y[n] = -m_rowScales[n] * m_resid[n];
|
|
}
|
|
m_resid_scaled = true;
|
|
} else {
|
|
for (size_t n = 0; n < neq_; n++) {
|
|
delta_y[n] = -m_resid[n];
|
|
}
|
|
}
|
|
|
|
// Factor the matrix using a standard Newton solve
|
|
m_conditionNumber = 1.0E300;
|
|
int info = 0;
|
|
if (!jac.factored()) {
|
|
if (useQR) {
|
|
info = jac.factorQR();
|
|
} else {
|
|
info = jac.factor();
|
|
}
|
|
}
|
|
/*
|
|
* Find the condition number of the matrix
|
|
* If we have failed to factor, we will fall back to calculating and factoring a modified Hessian
|
|
*/
|
|
if (info == 0) {
|
|
doublereal rcond = 0.0;
|
|
if (useQR) {
|
|
rcond = jac.rcondQR();
|
|
} else {
|
|
doublereal a1norm = jac.oneNorm();
|
|
rcond = jac.rcond(a1norm);
|
|
}
|
|
if (rcond > 0.0) {
|
|
m_conditionNumber = 1.0 / rcond;
|
|
}
|
|
} else {
|
|
m_conditionNumber = 1.0E300;
|
|
newtonGood = false;
|
|
if (m_print_flag >= 1) {
|
|
printf("\t\t doAffineNewtonSolve: ");
|
|
if (useQR) {
|
|
printf("factorQR()");
|
|
} else {
|
|
printf("factor()");
|
|
}
|
|
printf(" returned with info = %d, indicating a zero row or column\n", info);
|
|
}
|
|
}
|
|
bool doHessian = false;
|
|
if (s_doBothSolvesAndCompare) {
|
|
doHessian = true;
|
|
}
|
|
if (m_conditionNumber < 1.0E7) {
|
|
if (m_print_flag >= 4) {
|
|
printf("\t\t doAffineNewtonSolve: Condition number = %g during regular solve\n", m_conditionNumber);
|
|
}
|
|
|
|
/*
|
|
* Solve the system -> This also involves inverting the matrix
|
|
*/
|
|
int info = jac.solve(DATA_PTR(delta_y));
|
|
if (info) {
|
|
if (m_print_flag >= 2) {
|
|
printf("\t\t doAffineNewtonSolve() ERROR: QRSolve returned INFO = %d. Switching to Hessian solve\n", info);
|
|
}
|
|
doHessian = true;
|
|
newtonGood = false;
|
|
}
|
|
/*
|
|
* reverse the column scaling if there was any on a successful solve
|
|
*/
|
|
if (m_colScaling) {
|
|
for (size_t irow = 0; irow < neq_; irow++) {
|
|
delta_y[irow] = delta_y[irow] * m_colScales[irow];
|
|
}
|
|
}
|
|
|
|
} else {
|
|
if (jac.matrixType_ == 1) {
|
|
newtonGood = true;
|
|
if (m_print_flag >= 3) {
|
|
printf("\t\t doAffineNewtonSolve() WARNING: Condition number too large, %g, But Banded Hessian solve "
|
|
"not implemented yet \n", m_conditionNumber);
|
|
}
|
|
} else {
|
|
doHessian = true;
|
|
newtonGood = false;
|
|
if (m_print_flag >= 3) {
|
|
printf("\t\t doAffineNewtonSolve() WARNING: Condition number too large, %g. Doing a Hessian solve \n", m_conditionNumber);
|
|
}
|
|
}
|
|
}
|
|
|
|
if (doHessian) {
|
|
// Store the old value for later comparison
|
|
|
|
delyNewton = mdp::mdp_alloc_dbl_1((int) neq_, MDP_DBL_NOINIT);
|
|
for (size_t irow = 0; irow < neq_; irow++) {
|
|
delyNewton[irow] = delta_y[irow];
|
|
}
|
|
|
|
// Get memory if not done before
|
|
if (HessianPtr_ == 0) {
|
|
HessianPtr_ = jac.duplMyselfAsGeneralMatrix();
|
|
}
|
|
|
|
/*
|
|
* Calculate the symmetric Hessian
|
|
*/
|
|
GeneralMatrix& hessian = *HessianPtr_;
|
|
GeneralMatrix& jacCopy = *jacCopyPtr_;
|
|
hessian.zero();
|
|
if (m_rowScaling) {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
for (size_t j = i; j < neq_; j++) {
|
|
for (size_t k = 0; k < neq_; k++) {
|
|
hessian(i,j) += jacCopy(k,i) * jacCopy(k,j) * m_rowScales[k] * m_rowScales[k];
|
|
}
|
|
hessian(j,i) = hessian(i,j);
|
|
}
|
|
}
|
|
} else {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
for (size_t j = i; j < neq_; j++) {
|
|
for (size_t k = 0; k < neq_; k++) {
|
|
hessian(i,j) += jacCopy(k,i) * jacCopy(k,j);
|
|
}
|
|
hessian(j,i) = hessian(i,j);
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Calculate the matrix norm of the Hessian
|
|
*/
|
|
doublereal hnorm = 0.0;
|
|
doublereal hcol = 0.0;
|
|
if (m_colScaling) {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
for (size_t j = i; j < neq_; j++) {
|
|
hcol += fabs(hessian(j,i)) * m_colScales[j];
|
|
}
|
|
for (size_t j = i+1; j < neq_; j++) {
|
|
hcol += fabs(hessian(i,j)) * m_colScales[j];
|
|
}
|
|
hcol *= m_colScales[i];
|
|
if (hcol > hnorm) {
|
|
hnorm = hcol;
|
|
}
|
|
}
|
|
} else {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
for (size_t j = i; j < neq_; j++) {
|
|
hcol += fabs(hessian(j,i));
|
|
}
|
|
for (size_t j = i+1; j < neq_; j++) {
|
|
hcol += fabs(hessian(i,j));
|
|
}
|
|
if (hcol > hnorm) {
|
|
hnorm = hcol;
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* Add junk to the Hessian diagonal
|
|
* -> Note, testing indicates that this will get too big for ill-conditioned systems.
|
|
*/
|
|
hcol = sqrt(1.0*neq_) * 1.0E-7 * hnorm;
|
|
#ifdef DEBUG_HKM_NOT
|
|
if (hcol > 1.0) {
|
|
hcol = 1.0E1;
|
|
}
|
|
#endif
|
|
if (m_colScaling) {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
hessian(i,i) += hcol / (m_colScales[i] * m_colScales[i]);
|
|
}
|
|
} else {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
hessian(i,i) += hcol;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Factor the Hessian
|
|
*/
|
|
int info = 0;
|
|
ct_dpotrf(ctlapack::UpperTriangular, neq_, &(*(HessianPtr_->begin())), neq_, info);
|
|
if (info) {
|
|
if (m_print_flag >= 2) {
|
|
printf("\t\t doAffineNewtonSolve() ERROR: Hessian isn't positive definate DPOTRF returned INFO = %d\n", info);
|
|
}
|
|
return info;
|
|
}
|
|
|
|
// doublereal *JTF = delta_y;
|
|
doublereal* delyH = mdp::mdp_alloc_dbl_1((int) neq_, MDP_DBL_NOINIT);
|
|
// First recalculate the scaled residual. It got wiped out doing the newton solve
|
|
if (m_rowScaling) {
|
|
for (size_t n = 0; n < neq_; n++) {
|
|
delyH[n] = -m_rowScales[n] * m_resid[n];
|
|
}
|
|
} else {
|
|
for (size_t n = 0; n < neq_; n++) {
|
|
delyH[n] = -m_resid[n];
|
|
}
|
|
}
|
|
|
|
if (m_rowScaling) {
|
|
for (size_t j = 0; j < neq_; j++) {
|
|
delta_y[j] = 0.0;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
delta_y[j] += delyH[i] * jacCopy(i,j) * m_rowScales[i];
|
|
}
|
|
}
|
|
} else {
|
|
for (size_t j = 0; j < neq_; j++) {
|
|
delta_y[j] = 0.0;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
delta_y[j] += delyH[i] * jacCopy(i,j);
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Solve the factored Hessian System
|
|
*/
|
|
ct_dpotrs(ctlapack::UpperTriangular, neq_, 1,&(*(hessian.begin())), neq_, delta_y, neq_, info);
|
|
if (info) {
|
|
if (m_print_flag >= 2) {
|
|
printf("\t\t NonlinearSolver::doAffineNewtonSolve() ERROR: DPOTRS returned INFO = %d\n", info);
|
|
}
|
|
return info;
|
|
}
|
|
/*
|
|
* reverse the column scaling if there was any.
|
|
*/
|
|
if (m_colScaling) {
|
|
for (size_t irow = 0; irow < neq_; irow++) {
|
|
delta_y[irow] = delta_y[irow] * m_colScales[irow];
|
|
}
|
|
}
|
|
|
|
|
|
if (doDogLeg_ && m_print_flag > 7) {
|
|
double normNewt = solnErrorNorm(CONSTD_DATA_PTR(delyNewton));
|
|
double normHess = solnErrorNorm(CONSTD_DATA_PTR(delta_y));
|
|
printf("\t\t doAffineNewtonSolve(): Printout Comparison between Hessian deltaX and Newton deltaX\n");
|
|
|
|
printf("\t\t I Hessian+Junk Newton");
|
|
if (newtonGood || s_alwaysAssumeNewtonGood) {
|
|
printf(" (USING NEWTON DIRECTION)\n");
|
|
} else {
|
|
printf(" (USING HESSIAN DIRECTION)\n");
|
|
}
|
|
printf("\t\t Norm: %12.4E %12.4E\n", normHess, normNewt);
|
|
|
|
printf("\t\t --------------------------------------------------------\n");
|
|
for (size_t i =0; i < neq_; i++) {
|
|
printf("\t\t %3s %13.5E %13.5E\n",
|
|
int2str(i).c_str(), delta_y[i], delyNewton[i]);
|
|
}
|
|
printf("\t\t --------------------------------------------------------\n");
|
|
} else if (doDogLeg_ && m_print_flag >= 4) {
|
|
double normNewt = solnErrorNorm(CONSTD_DATA_PTR(delyNewton));
|
|
double normHess = solnErrorNorm(CONSTD_DATA_PTR(delta_y));
|
|
printf("\t\t doAffineNewtonSolve(): Hessian update norm = %12.4E \n"
|
|
"\t\t Newton update norm = %12.4E \n", normHess, normNewt);
|
|
if (newtonGood || s_alwaysAssumeNewtonGood) {
|
|
printf("\t\t (USING NEWTON DIRECTION)\n");
|
|
} else {
|
|
printf("\t\t (USING HESSIAN DIRECTION)\n");
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Choose the delta_y to use
|
|
*/
|
|
if (newtonGood || s_alwaysAssumeNewtonGood) {
|
|
mdp::mdp_copy_dbl_1(DATA_PTR(delta_y), CONSTD_DATA_PTR(delyNewton), (int) neq_);
|
|
}
|
|
mdp::mdp_safe_free((void**) &delyH);
|
|
mdp::mdp_safe_free((void**) &delyNewton);
|
|
}
|
|
|
|
#ifdef DEBUG_JAC
|
|
if (printJacContributions) {
|
|
for (int iNum = 0; iNum < numRows; iNum++) {
|
|
if (iNum > 0) {
|
|
focusRow++;
|
|
}
|
|
doublereal dsum = 0.0;
|
|
vector_fp& Jdata = jacBack.data();
|
|
doublereal dRow = Jdata[neq_ * focusRow + focusRow];
|
|
printf("\n Details on delta_Y for row %d \n", focusRow);
|
|
printf(" Value before = %15.5e, delta = %15.5e,"
|
|
"value after = %15.5e\n", y_curr[focusRow],
|
|
delta_y[focusRow], y_curr[focusRow] + delta_y[focusRow]);
|
|
if (!freshJac) {
|
|
printf(" Old Jacobian\n");
|
|
}
|
|
printf(" col delta_y aij "
|
|
"contrib \n");
|
|
printf("-----------------------------------------------------------------------------------------------\n");
|
|
printf(" Res(%d) %15.5e %15.5e %15.5e (Res = %g)\n",
|
|
focusRow, delta_y[focusRow],
|
|
dRow, RRow[iNum] / dRow, RRow[iNum]);
|
|
dsum += RRow[iNum] / dRow;
|
|
for (int ii = 0; ii < neq_; ii++) {
|
|
if (ii != focusRow) {
|
|
doublereal aij = Jdata[neq_ * ii + focusRow];
|
|
doublereal contrib = aij * delta_y[ii] * (-1.0) / dRow;
|
|
dsum += contrib;
|
|
if (fabs(contrib) > Pcutoff) {
|
|
printf("%6d %15.5e %15.5e %15.5e\n", ii,
|
|
delta_y[ii] , aij, contrib);
|
|
}
|
|
}
|
|
}
|
|
printf("-----------------------------------------------------------------------------------------------\n");
|
|
printf(" %15.5e %15.5e\n",
|
|
delta_y[focusRow], dsum);
|
|
}
|
|
}
|
|
|
|
#endif
|
|
|
|
m_numTotalLinearSolves++;
|
|
m_numLocalLinearSolves++;
|
|
return info;
|
|
|
|
}
|
|
//====================================================================================================================
|
|
// Do a steepest descent calculation
|
|
/*
|
|
* This call must be made on the unfactored jacobian!
|
|
*/
|
|
doublereal NonlinearSolver::doCauchyPointSolve(GeneralMatrix& jac)
|
|
{
|
|
doublereal rowFac = 1.0;
|
|
doublereal colFac = 1.0;
|
|
doublereal normSoln;
|
|
// Calculate the descent direction
|
|
/*
|
|
* For confirmation of the scaling factors, see Dennis and Schnabel p, 152, p, 156 and my notes
|
|
*
|
|
* The colFac and rowFac values are used to eliminate the scaling of the matrix from the
|
|
* actual equation
|
|
*
|
|
* Here we calculate the steepest descent direction. This is equation (11) in the notes. It is
|
|
* stored in deltaX_CP_[].The value corresponds to d_descent[].
|
|
*/
|
|
for (size_t j = 0; j < neq_; j++) {
|
|
deltaX_CP_[j] = 0.0;
|
|
if (m_colScaling) {
|
|
colFac = 1.0 / m_colScales[j];
|
|
}
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
if (m_rowScaling) {
|
|
rowFac = 1.0 / m_rowScales[i];
|
|
}
|
|
deltaX_CP_[j] -= m_resid[i] * jac(i,j) * colFac * rowFac * m_ewt[j] * m_ewt[j]
|
|
/ (m_residWts[i] * m_residWts[i]);
|
|
#ifdef DEBUG_MODE
|
|
mdp::checkFinite(deltaX_CP_[j]);
|
|
#endif
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Calculate J_hat d_y_descent. This is formula 18 in the notes.
|
|
*/
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
Jd_[i] = 0.0;
|
|
if (m_rowScaling) {
|
|
rowFac = 1.0 / m_rowScales[i];
|
|
} else {
|
|
rowFac = 1.0;
|
|
}
|
|
for (size_t j = 0; j < neq_; j++) {
|
|
if (m_colScaling) {
|
|
colFac = 1.0 / m_colScales[j];
|
|
}
|
|
Jd_[i] += deltaX_CP_[j] * jac(i,j) * rowFac * colFac / m_residWts[i];
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Calculate the distance along the steepest descent until the Cauchy point
|
|
* This is Eqn. 17 in the notes.
|
|
*/
|
|
RJd_norm_ = 0.0;
|
|
JdJd_norm_ = 0.0;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
RJd_norm_ += m_resid[i] * Jd_[i] / m_residWts[i];
|
|
JdJd_norm_ += Jd_[i] * Jd_[i];
|
|
}
|
|
//if (RJd_norm_ > -1.0E-300) {
|
|
// printf("we are here: zero residual\n");
|
|
//}
|
|
if (fabs(JdJd_norm_) < 1.0E-290) {
|
|
if (fabs(RJd_norm_) < 1.0E-300) {
|
|
lambdaStar_ = 0.0;
|
|
} else {
|
|
throw CanteraError("NonlinearSolver::doCauchyPointSolve()", "Unexpected condition: norms are zero");
|
|
}
|
|
} else {
|
|
lambdaStar_ = - RJd_norm_ / (JdJd_norm_);
|
|
}
|
|
|
|
/*
|
|
* Now we modify the steepest descent vector such that its length is equal to the
|
|
* Cauchy distance. From now on, if we want to recreate the descent vector, we have
|
|
* to unnormalize it by dividing by lambdaStar_.
|
|
*/
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
deltaX_CP_[i] *= lambdaStar_;
|
|
}
|
|
|
|
|
|
doublereal normResid02 = m_normResid_0 * m_normResid_0 * neq_;
|
|
|
|
/*
|
|
* Calculate the expected square of the risdual at the Cauchy point if the linear model is correct
|
|
*/
|
|
if (fabs(JdJd_norm_) < 1.0E-290) {
|
|
residNorm2Cauchy_ = normResid02;
|
|
} else {
|
|
residNorm2Cauchy_ = normResid02 - RJd_norm_ * RJd_norm_ / (JdJd_norm_);
|
|
}
|
|
|
|
|
|
// Extra printout section
|
|
if (m_print_flag > 2) {
|
|
// Calculate the expected residual at the Cauchy point if the linear model is correct
|
|
doublereal residCauchy = 0.0;
|
|
if (residNorm2Cauchy_ > 0.0) {
|
|
residCauchy = sqrt(residNorm2Cauchy_ / neq_);
|
|
} else {
|
|
if (fabs(JdJd_norm_) < 1.0E-290) {
|
|
residCauchy = m_normResid_0;
|
|
} else {
|
|
residCauchy = m_normResid_0 - sqrt(RJd_norm_ * RJd_norm_ / (JdJd_norm_));
|
|
}
|
|
}
|
|
// Compute the weighted norm of the undamped step size descentDir_[]
|
|
if ((s_print_DogLeg || doDogLeg_) && m_print_flag >= 6) {
|
|
normSoln = solnErrorNorm(DATA_PTR(deltaX_CP_), "SteepestDescentDir", 10);
|
|
} else {
|
|
normSoln = solnErrorNorm(DATA_PTR(deltaX_CP_), "SteepestDescentDir", 0);
|
|
}
|
|
if ((s_print_DogLeg || doDogLeg_) && m_print_flag >= 5) {
|
|
printf("\t\t doCauchyPointSolve: Steepest descent to Cauchy point: \n");
|
|
printf("\t\t\t R0 = %g \n", m_normResid_0);
|
|
printf("\t\t\t Rpred = %g\n", residCauchy);
|
|
printf("\t\t\t Rjd = %g\n", RJd_norm_);
|
|
printf("\t\t\t JdJd = %g\n", JdJd_norm_);
|
|
printf("\t\t\t deltaX = %g\n", normSoln);
|
|
printf("\t\t\t lambda = %g\n", lambdaStar_);
|
|
}
|
|
} else {
|
|
// Calculate the norm of the Cauchy solution update in any case
|
|
normSoln = solnErrorNorm(DATA_PTR(deltaX_CP_), "SteepestDescentDir", 0);
|
|
}
|
|
return normSoln;
|
|
}
|
|
//===================================================================================================================
|
|
void NonlinearSolver::descentComparison(doublereal time_curr, doublereal* ydot0, doublereal* ydot1, int& numTrials)
|
|
{
|
|
doublereal ff = 1.0E-5;
|
|
doublereal ffNewt = 1.0E-5;
|
|
doublereal* y_n_1 = DATA_PTR(m_wksp);
|
|
doublereal cauchyDistanceNorm = solnErrorNorm(DATA_PTR(deltaX_CP_));
|
|
if (cauchyDistanceNorm < 1.0E-2) {
|
|
ff = 1.0E-9 / cauchyDistanceNorm;
|
|
if (ff > 1.0E-2) {
|
|
ff = 1.0E-2;
|
|
}
|
|
}
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
y_n_1[i] = m_y_n_curr[i] + ff * deltaX_CP_[i];
|
|
}
|
|
/*
|
|
* Calculate the residual that would result if y1[] were the new solution vector
|
|
* -> m_resid[] contains the result of the residual calculation
|
|
*/
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
doResidualCalc(time_curr, solnType_, y_n_1, ydot1, Base_LaggedSolutionComponents);
|
|
} else {
|
|
doResidualCalc(time_curr, solnType_, y_n_1, ydot0, Base_LaggedSolutionComponents);
|
|
}
|
|
|
|
doublereal normResid02 = m_normResid_0 * m_normResid_0 * neq_;
|
|
doublereal residSteep = residErrorNorm(DATA_PTR(m_resid));
|
|
doublereal residSteep2 = residSteep * residSteep * neq_;
|
|
doublereal funcDecreaseSD = 0.5 * (residSteep2 - normResid02) / (ff * cauchyDistanceNorm);
|
|
|
|
doublereal sNewt = solnErrorNorm(DATA_PTR(deltaX_Newton_));
|
|
if (sNewt > 1.0) {
|
|
ffNewt = ffNewt / sNewt;
|
|
}
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
y_n_1[i] = m_y_n_curr[i] + ffNewt * deltaX_Newton_[i];
|
|
}
|
|
/*
|
|
* Calculate the residual that would result if y1[] were the new solution vector.
|
|
* Here we use the lagged solution components in the residual calculation as well. We are
|
|
* interested in the linear model and its agreement with the nonlinear model.
|
|
*
|
|
* -> m_resid[] contains the result of the residual calculation
|
|
*/
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
doResidualCalc(time_curr, solnType_, y_n_1, ydot1, Base_LaggedSolutionComponents);
|
|
} else {
|
|
doResidualCalc(time_curr, solnType_, y_n_1, ydot0, Base_LaggedSolutionComponents);
|
|
}
|
|
doublereal residNewt = residErrorNorm(DATA_PTR(m_resid));
|
|
doublereal residNewt2 = residNewt * residNewt * neq_;
|
|
|
|
doublereal funcDecreaseNewt2 = 0.5 * (residNewt2 - normResid02) / (ffNewt * sNewt);
|
|
|
|
// This is the expected initial rate of decrease in the Cauchy direction.
|
|
// -> This is Eqn. 29 = Rhat dot Jhat dy / || d ||
|
|
doublereal funcDecreaseSDExp = RJd_norm_ / cauchyDistanceNorm * lambdaStar_;
|
|
|
|
doublereal funcDecreaseNewtExp2 = - normResid02 / sNewt;
|
|
|
|
if (m_normResid_0 > 1.0E-100) {
|
|
ResidDecreaseSDExp_ = funcDecreaseSDExp / neq_ / m_normResid_0;
|
|
ResidDecreaseSD_ = funcDecreaseSD / neq_ / m_normResid_0;
|
|
ResidDecreaseNewtExp_ = funcDecreaseNewtExp2 / neq_ / m_normResid_0;
|
|
ResidDecreaseNewt_ = funcDecreaseNewt2 / neq_ / m_normResid_0;
|
|
} else {
|
|
ResidDecreaseSDExp_ = 0.0;
|
|
ResidDecreaseSD_ = funcDecreaseSD / neq_;
|
|
ResidDecreaseNewtExp_ = 0.0;
|
|
ResidDecreaseNewt_ = funcDecreaseNewt2 / neq_;
|
|
}
|
|
numTrials += 2;
|
|
|
|
/*
|
|
* HKM These have been shown to exactly match up.
|
|
* The steepest direction is always largest even when there are variable solution weights
|
|
*
|
|
* HKM When a hessian is used with junk on the diagonal, funcDecreaseNewtExp2 is no longer accurate as the
|
|
* direction gets significantly shorter with increasing condition number. This suggests an algorithm where the
|
|
* newton step from the Hessian should be increased so as to match funcDecreaseNewtExp2 = funcDecreaseNewt2.
|
|
* This roughly equals the ratio of the norms of the hessian and newton steps. This increased Newton step can
|
|
* then be used with the trust region double dogleg algorithm.
|
|
*/
|
|
if ((s_print_DogLeg && m_print_flag >= 3) || (doDogLeg_ && m_print_flag >= 5)) {
|
|
printf("\t\t descentComparison: initial rate of decrease of func in cauchy dir (expected) = %g\n", funcDecreaseSDExp);
|
|
printf("\t\t descentComparison: initial rate of decrease of func in cauchy dir = %g\n", funcDecreaseSD);
|
|
printf("\t\t descentComparison: initial rate of decrease of func in newton dir (expected) = %g\n", funcDecreaseNewtExp2);
|
|
printf("\t\t descentComparison: initial rate of decrease of func in newton dir = %g\n", funcDecreaseNewt2);
|
|
}
|
|
if ((s_print_DogLeg && m_print_flag >= 3) || (doDogLeg_ && m_print_flag >= 4)) {
|
|
printf("\t\t descentComparison: initial rate of decrease of Resid in cauchy dir (expected) = %g\n", ResidDecreaseSDExp_);
|
|
printf("\t\t descentComparison: initial rate of decrease of Resid in cauchy dir = %g\n", ResidDecreaseSD_);
|
|
printf("\t\t descentComparison: initial rate of decrease of Resid in newton dir (expected) = %g\n", ResidDecreaseNewtExp_);
|
|
printf("\t\t descentComparison: initial rate of decrease of Resid in newton dir = %g\n", ResidDecreaseNewt_);
|
|
}
|
|
|
|
if ((s_print_DogLeg && m_print_flag >= 5) || (doDogLeg_ && m_print_flag >= 5)) {
|
|
if (funcDecreaseNewt2 >= 0.0) {
|
|
printf("\t\t %13.5E %22.16E\n", funcDecreaseNewtExp2, m_normResid_0);
|
|
double ff = ffNewt * 1.0E-5;
|
|
for (int ii = 0; ii < 13; ii++) {
|
|
ff *= 10.;
|
|
if (ii == 12) {
|
|
ff = ffNewt;
|
|
}
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
y_n_1[i] = m_y_n_curr[i] + ff * deltaX_Newton_[i];
|
|
}
|
|
numTrials += 1;
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
doResidualCalc(time_curr, solnType_, y_n_1, ydot1, Base_LaggedSolutionComponents);
|
|
} else {
|
|
doResidualCalc(time_curr, solnType_, y_n_1, ydot0, Base_LaggedSolutionComponents);
|
|
}
|
|
residNewt = residErrorNorm(DATA_PTR(m_resid));
|
|
residNewt2 = residNewt * residNewt * neq_;
|
|
funcDecreaseNewt2 = 0.5 * (residNewt2 - normResid02) / (ff * sNewt);
|
|
printf("\t\t %10.3E %13.5E %22.16E\n", ff, funcDecreaseNewt2, residNewt);
|
|
}
|
|
|
|
}
|
|
|
|
|
|
}
|
|
|
|
}
|
|
|
|
//====================================================================================================================
|
|
// Setup the parameters for the double dog leg
|
|
/*
|
|
* The calls to the doCauchySolve() and doNewtonSolve() routines are done at the main level. This routine comes
|
|
* after those calls. We calculate the point Nuu_ here, the distances of the dog-legs,
|
|
* and the norms of the CP and Newton points in terms of the trust vectors.
|
|
*/
|
|
void NonlinearSolver::setupDoubleDogleg()
|
|
{
|
|
/*
|
|
* Gamma = ||grad f ||**4
|
|
* ---------------------------------------------
|
|
* (grad f)T H (grad f) (grad f)T H-1 (grad f)
|
|
*/
|
|
// doublereal sumG = 0.0;
|
|
// doublereal sumH = 0.0;
|
|
// for (int i = 0; i < neq_; i++) {
|
|
// sumG = deltax_cp_[i] * deltax_cp_[i];
|
|
// sumH = deltax_cp_[i] * newtDir[i];
|
|
// }
|
|
// double fac1 = sumG / lambdaStar_;
|
|
// double fac2 = sumH / lambdaStar_;
|
|
// double gamma = fac1 / fac2;
|
|
// doublereal gamma = m_normDeltaSoln_CP / m_normDeltaSoln_Newton;
|
|
/*
|
|
* This hasn't worked. so will do it heuristically. One issue is that the newton
|
|
* direction is not the inverse of the Hessian times the gradient. The Hession
|
|
* is the matrix squared. Until I have the inverse of the Hessian from QR factorization
|
|
* I may not be able to do it this way.
|
|
*/
|
|
|
|
/*
|
|
* Heuristic algorithm - Find out where on the Newton line the residual is the same
|
|
* as the residual at the cauchy point. Then, go halfway to
|
|
* the newton point and call that Nuu.
|
|
* Maybe we need to check that the linearized residual is
|
|
* monotonic along that line. However, we haven't needed to yet.
|
|
*/
|
|
doublereal residSteepLin = expectedResidLeg(0, 1.0);
|
|
doublereal Nres2CP = residSteepLin * residSteepLin * neq_;
|
|
doublereal Nres2_o = m_normResid_0 * m_normResid_0 * neq_;
|
|
doublereal a = Nres2CP / Nres2_o;
|
|
doublereal betaEqual = (2.0 - sqrt(4.0 - 4 * (1.0 - a))) / 2.0;
|
|
doublereal beta = (1.0 + betaEqual) / 2.0;
|
|
|
|
|
|
Nuu_ = beta;
|
|
|
|
dist_R0_ = m_normDeltaSoln_CP;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
m_wksp[i] = Nuu_ * deltaX_Newton_[i] - deltaX_CP_[i];
|
|
}
|
|
dist_R1_ = solnErrorNorm(DATA_PTR(m_wksp));
|
|
dist_R2_ = (1.0 - Nuu_) * m_normDeltaSoln_Newton;
|
|
dist_Total_ = dist_R0_ + dist_R1_ + dist_R2_;
|
|
|
|
/*
|
|
* Calculate the trust distances
|
|
*/
|
|
normTrust_Newton_ = calcTrustDistance(deltaX_Newton_);
|
|
normTrust_CP_ = calcTrustDistance(deltaX_CP_);
|
|
|
|
}
|
|
//====================================================================================================================
|
|
// Change the global lambda coordinate into the (leg,alpha) coordinate for the double dogleg
|
|
/*
|
|
* @param lambda Global value of the distance along the double dogleg
|
|
* @param alpha relative value along the particular leg
|
|
*
|
|
* @return Returns the leg number ( 0, 1, or 2).
|
|
*/
|
|
int NonlinearSolver::lambdaToLeg(const doublereal lambda, doublereal& alpha) const
|
|
{
|
|
|
|
if (lambda < dist_R0_ / dist_Total_) {
|
|
alpha = lambda * dist_Total_ / dist_R0_;
|
|
return 0;
|
|
} else if (lambda < ((dist_R0_ + dist_R1_)/ dist_Total_)) {
|
|
alpha = (lambda * dist_Total_ - dist_R0_) / dist_R1_;
|
|
return 1;
|
|
}
|
|
alpha = (lambda * dist_Total_ - dist_R0_ - dist_R1_) / dist_R2_;
|
|
return 2;
|
|
}
|
|
//====================================================================================================================
|
|
// Calculated the expected residual along the double dogleg curve.
|
|
/*
|
|
* @param leg 0, 1, or 2 representing the curves of the dogleg
|
|
* @param alpha Relative distance along the particular curve.
|
|
*
|
|
* @return Returns the expected value of the residual at that point according to the quadratic model.
|
|
* The residual at the newton point will always be zero.
|
|
*/
|
|
doublereal NonlinearSolver::expectedResidLeg(int leg, doublereal alpha) const
|
|
{
|
|
|
|
doublereal resD2, res2, resNorm;
|
|
doublereal normResid02 = m_normResid_0 * m_normResid_0 * neq_;
|
|
|
|
if (leg == 0) {
|
|
/*
|
|
* We are on the steepest descent line
|
|
* along that line
|
|
* R2 = R2 + 2 lambda R dot Jd + lambda**2 Jd dot Jd
|
|
*/
|
|
|
|
doublereal tmp = - 2.0 * alpha + alpha * alpha;
|
|
doublereal tmp2 = - RJd_norm_ * lambdaStar_;
|
|
resD2 = tmp2 * tmp;
|
|
|
|
} else if (leg == 1) {
|
|
|
|
/*
|
|
* Same formula as above for lambda=1.
|
|
*/
|
|
doublereal tmp2 = - RJd_norm_ * lambdaStar_;
|
|
doublereal RdotJS = - tmp2;
|
|
doublereal JsJs = tmp2;
|
|
|
|
|
|
doublereal res0_2 = m_normResid_0 * m_normResid_0 * neq_;
|
|
|
|
res2 = res0_2 + (1.0 - alpha) * 2 * RdotJS - 2 * alpha * Nuu_ * res0_2
|
|
+ (1.0 - alpha) * (1.0 - alpha) * JsJs
|
|
+ alpha * alpha * Nuu_ * Nuu_ * res0_2
|
|
- 2 * alpha * Nuu_ * (1.0 - alpha) * RdotJS;
|
|
|
|
resNorm = sqrt(res2 / neq_);
|
|
return resNorm;
|
|
|
|
} else {
|
|
doublereal beta = Nuu_ + alpha * (1.0 - Nuu_);
|
|
doublereal tmp2 = normResid02;
|
|
doublereal tmp = 1.0 - 2.0 * beta + 1.0 * beta * beta - 1.0;
|
|
resD2 = tmp * tmp2;
|
|
}
|
|
|
|
res2 = m_normResid_0 * m_normResid_0 * neq_ + resD2;
|
|
if (res2 < 0.0) {
|
|
resNorm = m_normResid_0 - sqrt(resD2/neq_);
|
|
} else {
|
|
resNorm = sqrt(res2 / neq_);
|
|
}
|
|
|
|
return resNorm;
|
|
|
|
}
|
|
//====================================================================================================================
|
|
// Here we print out the residual at various points along the double dogleg, comparing against the quadratic model
|
|
// in a table format
|
|
/*
|
|
* @param time_curr INPUT current time
|
|
* @param ydot0 INPUT Current value of the derivative of the solution vector for non-time dependent
|
|
* determinations
|
|
* @param legBest OUTPUT leg of the dogleg that gives the lowest residual
|
|
* @param alphaBest OUTPUT distance along dogleg for best result.
|
|
*/
|
|
void NonlinearSolver::residualComparisonLeg(const doublereal time_curr, const doublereal* const ydot0, int& legBest,
|
|
doublereal& alphaBest) const
|
|
{
|
|
doublereal* y1 = DATA_PTR(m_wksp);
|
|
doublereal* ydot1 = DATA_PTR(m_wksp_2);
|
|
doublereal sLen;
|
|
doublereal alpha = 0;
|
|
|
|
doublereal residSteepBest = 1.0E300;
|
|
doublereal residSteepLinBest = 0.0;
|
|
if (s_print_DogLeg || (doDogLeg_ && m_print_flag > 6)) {
|
|
printf("\t\t residualComparisonLeg() \n");
|
|
printf("\t\t Point StepLen Residual_Actual Residual_Linear RelativeMatch\n");
|
|
}
|
|
// First compare at 1/4 along SD curve
|
|
std::vector<doublereal> alphaT;
|
|
alphaT.push_back(0.00);
|
|
alphaT.push_back(0.01);
|
|
alphaT.push_back(0.1);
|
|
alphaT.push_back(0.25);
|
|
alphaT.push_back(0.50);
|
|
alphaT.push_back(0.75);
|
|
alphaT.push_back(1.0);
|
|
for (size_t iteration = 0; iteration < alphaT.size(); iteration++) {
|
|
alpha = alphaT[iteration];
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
y1[i] = m_y_n_curr[i] + alpha * deltaX_CP_[i];
|
|
}
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
calc_ydot(m_order, y1, ydot1);
|
|
}
|
|
sLen = alpha * solnErrorNorm(DATA_PTR(deltaX_CP_));
|
|
/*
|
|
* Calculate the residual that would result if y1[] were the new solution vector
|
|
* -> m_resid[] contains the result of the residual calculation
|
|
*/
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
doResidualCalc(time_curr, solnType_, y1, ydot1, Base_LaggedSolutionComponents);
|
|
} else {
|
|
doResidualCalc(time_curr, solnType_, y1, ydot0, Base_LaggedSolutionComponents);
|
|
}
|
|
|
|
|
|
doublereal residSteep = residErrorNorm(DATA_PTR(m_resid));
|
|
doublereal residSteepLin = expectedResidLeg(0, alpha);
|
|
if (residSteep < residSteepBest) {
|
|
legBest = 0;
|
|
alphaBest = alpha;
|
|
residSteepBest = residSteep;
|
|
residSteepLinBest = residSteepLin;
|
|
}
|
|
|
|
doublereal relFit = (residSteep - residSteepLin) / (fabs(residSteepLin) + 1.0E-10);
|
|
if (s_print_DogLeg || (doDogLeg_ && m_print_flag > 6)) {
|
|
printf("\t\t (%2d - % 10.3g) % 15.8E % 15.8E % 15.8E % 15.8E\n", 0, alpha, sLen, residSteep, residSteepLin , relFit);
|
|
}
|
|
}
|
|
|
|
for (size_t iteration = 0; iteration < alphaT.size(); iteration++) {
|
|
doublereal alpha = alphaT[iteration];
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
y1[i] = m_y_n_curr[i] + (1.0 - alpha) * deltaX_CP_[i];
|
|
y1[i] += alpha * Nuu_ * deltaX_Newton_[i];
|
|
}
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
calc_ydot(m_order, y1, ydot1);
|
|
}
|
|
/*
|
|
* Calculate the residual that would result if y1[] were the new solution vector
|
|
* -> m_resid[] contains the result of the residual calculation
|
|
*/
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
doResidualCalc(time_curr, solnType_, y1, ydot1, Base_LaggedSolutionComponents);
|
|
} else {
|
|
doResidualCalc(time_curr, solnType_, y1, ydot0, Base_LaggedSolutionComponents);
|
|
}
|
|
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
y1[i] -= m_y_n_curr[i];
|
|
}
|
|
sLen = solnErrorNorm(DATA_PTR(y1));
|
|
|
|
doublereal residSteep = residErrorNorm(DATA_PTR(m_resid));
|
|
doublereal residSteepLin = expectedResidLeg(1, alpha);
|
|
if (residSteep < residSteepBest) {
|
|
legBest = 1;
|
|
alphaBest = alpha;
|
|
residSteepBest = residSteep;
|
|
residSteepLinBest = residSteepLin;
|
|
}
|
|
|
|
doublereal relFit = (residSteep - residSteepLin) / (fabs(residSteepLin) + 1.0E-10);
|
|
if (s_print_DogLeg || (doDogLeg_ && m_print_flag > 6)) {
|
|
printf("\t\t (%2d - % 10.3g) % 15.8E % 15.8E % 15.8E % 15.8E\n", 1, alpha, sLen, residSteep, residSteepLin , relFit);
|
|
}
|
|
}
|
|
|
|
for (size_t iteration = 0; iteration < alphaT.size(); iteration++) {
|
|
doublereal alpha = alphaT[iteration];
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
y1[i] = m_y_n_curr[i] + (Nuu_ + alpha * (1.0 - Nuu_))* deltaX_Newton_[i];
|
|
}
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
calc_ydot(m_order, y1, ydot1);
|
|
}
|
|
sLen = (Nuu_ + alpha * (1.0 - Nuu_)) * solnErrorNorm(DATA_PTR(deltaX_Newton_));
|
|
/*
|
|
* Calculate the residual that would result if y1[] were the new solution vector
|
|
* -> m_resid[] contains the result of the residual calculation
|
|
*/
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
doResidualCalc(time_curr, solnType_, y1, ydot1, Base_LaggedSolutionComponents);
|
|
} else {
|
|
doResidualCalc(time_curr, solnType_, y1, ydot0, Base_LaggedSolutionComponents);
|
|
}
|
|
|
|
|
|
|
|
doublereal residSteep = residErrorNorm(DATA_PTR(m_resid));
|
|
doublereal residSteepLin = expectedResidLeg(2, alpha);
|
|
if (residSteep < residSteepBest) {
|
|
legBest = 2;
|
|
alphaBest = alpha;
|
|
residSteepBest = residSteep;
|
|
residSteepLinBest = residSteepLin;
|
|
}
|
|
doublereal relFit = (residSteep - residSteepLin) / (fabs(residSteepLin) + 1.0E-10);
|
|
if (s_print_DogLeg || (doDogLeg_ && m_print_flag > 6)) {
|
|
printf("\t\t (%2d - % 10.3g) % 15.8E % 15.8E % 15.8E % 15.8E\n", 2, alpha, sLen, residSteep, residSteepLin , relFit);
|
|
}
|
|
}
|
|
if (s_print_DogLeg || (doDogLeg_ && m_print_flag > 6)) {
|
|
printf("\t\t Best Result: \n");
|
|
doublereal relFit = (residSteepBest - residSteepLinBest) / (fabs(residSteepLinBest) + 1.0E-10);
|
|
if (m_print_flag <= 6) {
|
|
printf("\t\t Leg %2d alpha %5g: NonlinResid = %g LinResid = %g, relfit = %g\n",
|
|
legBest, alphaBest, residSteepBest, residSteepLinBest, relFit);
|
|
} else {
|
|
if (legBest == 0) {
|
|
sLen = alpha * solnErrorNorm(DATA_PTR(deltaX_CP_));
|
|
} else if (legBest == 1) {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
y1[i] = (1.0 - alphaBest) * deltaX_CP_[i];
|
|
y1[i] += alphaBest * Nuu_ * deltaX_Newton_[i];
|
|
}
|
|
sLen = solnErrorNorm(DATA_PTR(y1));
|
|
} else {
|
|
sLen = (Nuu_ + alpha * (1.0 - Nuu_)) * solnErrorNorm(DATA_PTR(deltaX_Newton_));
|
|
}
|
|
printf("\t\t (%2d - % 10.3g) % 15.8E % 15.8E % 15.8E % 15.8E\n", legBest, alphaBest, sLen,
|
|
residSteepBest, residSteepLinBest , relFit);
|
|
}
|
|
}
|
|
|
|
}
|
|
//====================================================================================================================
|
|
// Calculate the length of the current trust region in terms of the solution error norm
|
|
/*
|
|
* We carry out a norm of deltaX_trust_ first. Then, we multiply that value
|
|
* by trustDelta_
|
|
*/
|
|
doublereal NonlinearSolver::trustRegionLength() const
|
|
{
|
|
norm_deltaX_trust_ = solnErrorNorm(DATA_PTR(deltaX_trust_));
|
|
return (trustDelta_ * norm_deltaX_trust_);
|
|
}
|
|
//====================================================================================================================
|
|
void NonlinearSolver::setDefaultDeltaBoundsMagnitudes()
|
|
{
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
m_deltaStepMinimum[i] = 1000. * atolk_[i];
|
|
m_deltaStepMinimum[i] = std::max(m_deltaStepMinimum[i], 0.1 * fabs(m_y_n_curr[i]));
|
|
}
|
|
}
|
|
//====================================================================================================================
|
|
void NonlinearSolver::adjustUpStepMinimums()
|
|
{
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
doublereal goodVal = deltaX_trust_[i] * trustDelta_;
|
|
if (deltaX_trust_[i] * trustDelta_ > m_deltaStepMinimum[i]) {
|
|
m_deltaStepMinimum[i] = 1.1 * goodVal;
|
|
}
|
|
|
|
}
|
|
}
|
|
//====================================================================================================================
|
|
void NonlinearSolver::setDeltaBoundsMagnitudes(const doublereal* const deltaStepMinimum)
|
|
{
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
m_deltaStepMinimum[i] = deltaStepMinimum[i];
|
|
}
|
|
m_manualDeltaStepSet = 1;
|
|
}
|
|
//====================================================================================================================
|
|
/*
|
|
*
|
|
* Return the factor by which the undamped Newton step 'step0'
|
|
* must be multiplied in order to keep the update within the bounds of an accurate jacobian.
|
|
*
|
|
* The idea behind these is that the Jacobian couldn't possibly be representative, if the
|
|
* variable is changed by a lot. (true for nonlinear systems, false for linear systems)
|
|
* Maximum increase in variable in any one newton iteration:
|
|
* factor of 1.5
|
|
* Maximum decrease in variable in any one newton iteration:
|
|
* factor of 2
|
|
*
|
|
* @param y_n_curr Initial value of the solution vector
|
|
* @param step_1 initial proposed step size
|
|
*
|
|
* @return returns the damping factor
|
|
*/
|
|
double
|
|
NonlinearSolver::deltaBoundStep(const doublereal* const y_n_curr, const doublereal* const step_1)
|
|
{
|
|
|
|
size_t i_fbounds = 0;
|
|
int ifbd = 0;
|
|
int i_fbd = 0;
|
|
doublereal UPFAC = 2.0;
|
|
|
|
doublereal sameSign = 0.0;
|
|
doublereal ff;
|
|
doublereal f_delta_bounds = 1.0;
|
|
doublereal ff_alt;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
doublereal y_new = y_n_curr[i] + step_1[i];
|
|
sameSign = y_new * y_n_curr[i];
|
|
|
|
/*
|
|
* Now do a delta bounds
|
|
* Increase variables by a factor of UPFAC only
|
|
* decrease variables by a factor of 2 only
|
|
*/
|
|
ff = 1.0;
|
|
|
|
|
|
if (sameSign >= 0.0) {
|
|
if ((fabs(y_new) > UPFAC * fabs(y_n_curr[i])) &&
|
|
(fabs(y_new - y_n_curr[i]) > m_deltaStepMinimum[i])) {
|
|
ff = (UPFAC - 1.0) * fabs(y_n_curr[i]/(y_new - y_n_curr[i]));
|
|
ff_alt = fabs(m_deltaStepMinimum[i] / (y_new - y_n_curr[i]));
|
|
ff = std::max(ff, ff_alt);
|
|
ifbd = 1;
|
|
}
|
|
if ((fabs(2.0 * y_new) < fabs(y_n_curr[i])) &&
|
|
(fabs(y_new - y_n_curr[i]) > m_deltaStepMinimum[i])) {
|
|
ff = y_n_curr[i]/(y_new - y_n_curr[i]) * (1.0 - 2.0)/2.0;
|
|
ff_alt = fabs(m_deltaStepMinimum[i] / (y_new - y_n_curr[i]));
|
|
ff = std::max(ff, ff_alt);
|
|
ifbd = 0;
|
|
}
|
|
} else {
|
|
/*
|
|
* This handles the case where the value crosses the origin.
|
|
* - First we don't let it cross the origin until it's shrunk to the size of m_deltaStepMinimum[i]
|
|
*/
|
|
if (fabs(y_n_curr[i]) > m_deltaStepMinimum[i]) {
|
|
ff = y_n_curr[i]/(y_new - y_n_curr[i]) * (1.0 - 2.0)/2.0;
|
|
ff_alt = fabs(m_deltaStepMinimum[i] / (y_new - y_n_curr[i]));
|
|
ff = std::max(ff, ff_alt);
|
|
if (y_n_curr[i] >= 0.0) {
|
|
ifbd = 0;
|
|
} else {
|
|
ifbd = 1;
|
|
}
|
|
}
|
|
/*
|
|
* Second when it does cross the origin, we make sure that its magnitude is only 50% of the previous value.
|
|
*/
|
|
else if (fabs(y_new) > 0.5 * fabs(y_n_curr[i])) {
|
|
ff = y_n_curr[i]/(y_new - y_n_curr[i]) * (-1.5);
|
|
ff_alt = fabs(m_deltaStepMinimum[i] / (y_new - y_n_curr[i]));
|
|
ff = std::max(ff, ff_alt);
|
|
ifbd = 0;
|
|
}
|
|
}
|
|
|
|
if (ff < f_delta_bounds) {
|
|
f_delta_bounds = ff;
|
|
i_fbounds = i;
|
|
i_fbd = ifbd;
|
|
}
|
|
|
|
|
|
}
|
|
|
|
|
|
/*
|
|
* Report on any corrections
|
|
*/
|
|
if (m_print_flag >= 3) {
|
|
if (f_delta_bounds < 1.0) {
|
|
if (i_fbd) {
|
|
printf("\t\tdeltaBoundStep: Increase of Variable %s causing "
|
|
"delta damping of %g: origVal = %10.3g, undampedNew = %10.3g, dampedNew = %10.3g\n",
|
|
int2str(i_fbounds).c_str(), f_delta_bounds, y_n_curr[i_fbounds], y_n_curr[i_fbounds] + step_1[i_fbounds],
|
|
y_n_curr[i_fbounds] + f_delta_bounds * step_1[i_fbounds]);
|
|
} else {
|
|
printf("\t\tdeltaBoundStep: Decrease of variable %s causing"
|
|
"delta damping of %g: origVal = %10.3g, undampedNew = %10.3g, dampedNew = %10.3g\n",
|
|
int2str(i_fbounds).c_str(), f_delta_bounds, y_n_curr[i_fbounds], y_n_curr[i_fbounds] + step_1[i_fbounds],
|
|
y_n_curr[i_fbounds] + f_delta_bounds * step_1[i_fbounds]);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
return f_delta_bounds;
|
|
}
|
|
//====================================================================================================================
|
|
// Readjust the trust region vectors
|
|
/*
|
|
* The trust region is made up of the trust region vector calculation and the trustDelta_ value
|
|
* We periodically recalculate the trustVector_ values so that they renormalize to the
|
|
* correct length.
|
|
*/
|
|
void NonlinearSolver::readjustTrustVector()
|
|
{
|
|
doublereal trustDeltaOld = trustDelta_;
|
|
doublereal wtSum = 0.0;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
wtSum += m_ewt[i];
|
|
}
|
|
wtSum /= neq_;
|
|
doublereal trustNorm = solnErrorNorm(DATA_PTR(deltaX_trust_));
|
|
doublereal deltaXSizeOld = trustNorm;
|
|
doublereal trustNormGoal = trustNorm * trustDelta_;
|
|
|
|
// This is the size of each component.
|
|
// doublereal trustDeltaEach = trustDelta_ * trustNorm / neq_;
|
|
doublereal oldVal;
|
|
doublereal fabsy;
|
|
// we use the old value of the trust region as an indicator
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
oldVal = deltaX_trust_[i];
|
|
fabsy = fabs(m_y_n_curr[i]);
|
|
// First off make sure that each trust region vector is 1/2 the size of each variable or smaller
|
|
// unless overridden by the deltaStepMininum value.
|
|
// doublereal newValue = trustDeltaEach * m_ewt[i] / wtSum;
|
|
doublereal newValue = trustNormGoal * m_ewt[i];
|
|
if (newValue > 0.5 * fabsy) {
|
|
if (fabsy * 0.5 > m_deltaStepMinimum[i]) {
|
|
deltaX_trust_[i] = 0.5 * fabsy;
|
|
} else {
|
|
deltaX_trust_[i] = m_deltaStepMinimum[i];
|
|
}
|
|
} else {
|
|
if (newValue > 4.0 * oldVal) {
|
|
newValue = 4.0 * oldVal;
|
|
} else if (newValue < 0.25 * oldVal) {
|
|
newValue = 0.25 * oldVal;
|
|
if (deltaX_trust_[i] < m_deltaStepMinimum[i]) {
|
|
newValue = m_deltaStepMinimum[i];
|
|
}
|
|
}
|
|
deltaX_trust_[i] = newValue;
|
|
if (deltaX_trust_[i] > 0.75 * m_deltaStepMaximum[i]) {
|
|
deltaX_trust_[i] = 0.75 * m_deltaStepMaximum[i];
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Final renormalization.
|
|
norm_deltaX_trust_ = solnErrorNorm(DATA_PTR(deltaX_trust_));
|
|
doublereal sum = trustNormGoal / trustNorm;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
deltaX_trust_[i] = deltaX_trust_[i] * sum;
|
|
}
|
|
norm_deltaX_trust_ = solnErrorNorm(DATA_PTR(deltaX_trust_));
|
|
trustDelta_ = trustNormGoal / norm_deltaX_trust_;
|
|
|
|
if (doDogLeg_ && m_print_flag >= 4) {
|
|
printf("\t\t reajustTrustVector(): Trust size = %11.3E: Old deltaX size = %11.3E trustDelta_ = %11.3E\n"
|
|
"\t\t new deltaX size = %11.3E trustdelta_ = %11.3E\n",
|
|
trustNormGoal, deltaXSizeOld, trustDeltaOld, norm_deltaX_trust_, trustDelta_);
|
|
}
|
|
}
|
|
//====================================================================================================================
|
|
//! Initialize the size of the trust vector.
|
|
/*!
|
|
* The algorithm we use is to set it equal to the length of the Distance to the Cauchy point.
|
|
*/
|
|
void NonlinearSolver::initializeTrustRegion()
|
|
{
|
|
if (trustRegionInitializationMethod_ == 0) {
|
|
return;
|
|
}
|
|
if (trustRegionInitializationMethod_ == 1) {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
deltaX_trust_[i] = m_ewt[i] * trustRegionInitializationFactor_;
|
|
}
|
|
trustDelta_ = 1.0;
|
|
}
|
|
if (trustRegionInitializationMethod_ == 2) {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
deltaX_trust_[i] = m_ewt[i] * m_normDeltaSoln_CP * trustRegionInitializationFactor_;
|
|
}
|
|
doublereal cpd = calcTrustDistance(deltaX_CP_);
|
|
if ((doDogLeg_ && m_print_flag >= 4)) {
|
|
printf("\t\t initializeTrustRegion(): Relative Distance of Cauchy Vector wrt Trust Vector = %g\n", cpd);
|
|
}
|
|
trustDelta_ = trustDelta_ * cpd * trustRegionInitializationFactor_;
|
|
readjustTrustVector();
|
|
cpd = calcTrustDistance(deltaX_CP_);
|
|
if ((doDogLeg_ && m_print_flag >= 4)) {
|
|
printf("\t\t initializeTrustRegion(): Relative Distance of Cauchy Vector wrt Trust Vector = %g\n", cpd);
|
|
}
|
|
}
|
|
if (trustRegionInitializationMethod_ == 3) {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
deltaX_trust_[i] = m_ewt[i] * m_normDeltaSoln_Newton * trustRegionInitializationFactor_;
|
|
}
|
|
doublereal cpd = calcTrustDistance(deltaX_Newton_);
|
|
if ((doDogLeg_ && m_print_flag >= 4)) {
|
|
printf("\t\t initializeTrustRegion(): Relative Distance of Newton Vector wrt Trust Vector = %g\n", cpd);
|
|
}
|
|
trustDelta_ = trustDelta_ * cpd;
|
|
readjustTrustVector();
|
|
cpd = calcTrustDistance(deltaX_Newton_);
|
|
if ((doDogLeg_ && m_print_flag >= 4)) {
|
|
printf("\t\t initializeTrustRegion(): Relative Distance of Newton Vector wrt Trust Vector = %g\n", cpd);
|
|
}
|
|
}
|
|
}
|
|
|
|
//====================================================================================================================
|
|
// Fill a dogleg solution step vector
|
|
/*
|
|
* Previously, we have filled up deltaX_Newton_[], deltaX_CP_[], and Nuu_, so that
|
|
* this routine is straightforward.
|
|
*
|
|
* @param leg Leg of the dog leg you are on (0, 1, or 2)
|
|
* @param alpha Relative length along the dog length that you are on.
|
|
* @param deltaX Vector to be filled up
|
|
*/
|
|
void NonlinearSolver::fillDogLegStep(int leg, doublereal alpha, std::vector<doublereal> & deltaX) const
|
|
{
|
|
if (leg == 0) {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
deltaX[i] = alpha * deltaX_CP_[i];
|
|
}
|
|
} else if (leg == 2) {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
deltaX[i] = (alpha + (1.0 - alpha) * Nuu_) * deltaX_Newton_[i];
|
|
}
|
|
} else {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
deltaX[i] = deltaX_CP_[i] * (1.0 - alpha) + alpha * Nuu_ * deltaX_Newton_[i];
|
|
}
|
|
}
|
|
}
|
|
//====================================================================================================================
|
|
// Calculate the trust distance of a step in the solution variables
|
|
/*
|
|
* The trust distance is defined as the length of the step according to the norm wrt to the trust region.
|
|
* We calculate the trust distance by the following method
|
|
*
|
|
* trustDist = || delta_x dot 1/trustDeltaX_ || / trustDelta_
|
|
*
|
|
* @param deltaX Current value of deltaX
|
|
*/
|
|
doublereal NonlinearSolver::calcTrustDistance(std::vector<doublereal> const& deltaX) const
|
|
{
|
|
doublereal sum = 0.0;
|
|
doublereal tmp = 0.0;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
tmp = deltaX[i] / deltaX_trust_[i];
|
|
sum += tmp * tmp;
|
|
}
|
|
sum = sqrt(sum / neq_) / trustDelta_;
|
|
return sum;
|
|
}
|
|
//====================================================================================================================
|
|
// Given a trust distance, this routine calculates the intersection of the this distance with the
|
|
// double dogleg curve
|
|
/*
|
|
* @param trustDelta (INPUT) Value of the trust distance
|
|
* @param lambda (OUTPUT) Returns the internal coordinate of the double dogleg
|
|
* @param alpha (OUTPUT) Returns the relative distance along the appropriate leg
|
|
* @return leg (OUTPUT) Returns the leg ID (0, 1, or 2)
|
|
*/
|
|
int NonlinearSolver::calcTrustIntersection(doublereal trustDelta, doublereal& lambda, doublereal& alpha) const
|
|
{
|
|
doublereal dist;
|
|
if (normTrust_Newton_ < trustDelta) {
|
|
lambda = 1.0;
|
|
alpha = 1.0;
|
|
return 2;
|
|
}
|
|
|
|
if (normTrust_Newton_ * Nuu_ < trustDelta) {
|
|
alpha = (trustDelta - normTrust_Newton_ * Nuu_) / (normTrust_Newton_ - normTrust_Newton_ * Nuu_);
|
|
dist = dist_R0_ + dist_R1_ + alpha * dist_R2_;
|
|
lambda = dist / dist_Total_;
|
|
return 2;
|
|
}
|
|
if (normTrust_CP_ > trustDelta) {
|
|
lambda = 1.0;
|
|
dist = dist_R0_ * trustDelta / normTrust_CP_;
|
|
lambda = dist / dist_Total_;
|
|
alpha = trustDelta / normTrust_CP_;
|
|
return 0;
|
|
}
|
|
doublereal sumv = 0.0;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
sumv += (deltaX_Newton_[i] / deltaX_trust_[i]) * (deltaX_CP_[i] / deltaX_trust_[i]);
|
|
}
|
|
|
|
doublereal a = normTrust_Newton_ * normTrust_Newton_ * Nuu_ * Nuu_;
|
|
doublereal b = 2.0 * Nuu_ * sumv;
|
|
doublereal c = normTrust_CP_ * normTrust_CP_ - trustDelta * trustDelta;
|
|
|
|
alpha =(-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
|
|
|
|
|
|
dist = dist_R0_ + alpha * dist_R1_;
|
|
lambda = dist / dist_Total_;
|
|
return 1;
|
|
}
|
|
//====================================================================================================================
|
|
/*
|
|
*
|
|
* boundStep():
|
|
*
|
|
* Return the factor by which the undamped Newton step 'step0'
|
|
* must be multiplied in order to keep all solution components in
|
|
* all domains between their specified lower and upper bounds.
|
|
* Other bounds may be applied here as well.
|
|
*
|
|
* Currently the bounds are hard coded into this routine:
|
|
*
|
|
* Minimum value for all variables: - 0.01 * m_ewt[i]
|
|
* Maximum value = none.
|
|
*
|
|
* Thus, this means that all solution components are expected
|
|
* to be numerical greater than zero in the limit of time step
|
|
* truncation errors going to zero.
|
|
*
|
|
* Delta bounds: The idea behind these is that the Jacobian
|
|
* couldn't possibly be representative if the
|
|
* variable is changed by a lot. (true for
|
|
* nonlinear systems, false for linear systems)
|
|
* Maximum increase in variable in any one newton iteration:
|
|
* factor of 2
|
|
* Maximum decrease in variable in any one newton iteration:
|
|
* factor of 5
|
|
*/
|
|
doublereal NonlinearSolver::boundStep(const doublereal* const y, const doublereal* const step0)
|
|
{
|
|
size_t i_lower = npos;
|
|
doublereal fbound = 1.0, f_bounds = 1.0;
|
|
doublereal ff, y_new;
|
|
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
y_new = y[i] + step0[i];
|
|
/*
|
|
* Force the step to only take 80% a step towards the lower bounds
|
|
*/
|
|
if (step0[i] < 0.0) {
|
|
if (y_new < (y[i] + 0.8 * (m_y_low_bounds[i] - y[i]))) {
|
|
doublereal legalDelta = 0.8*(m_y_low_bounds[i] - y[i]);
|
|
ff = legalDelta / step0[i];
|
|
if (ff < f_bounds) {
|
|
f_bounds = ff;
|
|
i_lower = i;
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* Force the step to only take 80% a step towards the high bounds
|
|
*/
|
|
if (step0[i] > 0.0) {
|
|
if (y_new > (y[i] + 0.8 * (m_y_high_bounds[i] - y[i]))) {
|
|
doublereal legalDelta = 0.8*(m_y_high_bounds[i] - y[i]);
|
|
ff = legalDelta / step0[i];
|
|
if (ff < f_bounds) {
|
|
f_bounds = ff;
|
|
i_lower = i;
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
/*
|
|
* Report on any corrections
|
|
*/
|
|
if (m_print_flag >= 3) {
|
|
if (f_bounds != 1.0) {
|
|
printf("\t\tboundStep: Variable %s causing bounds damping of %g\n",
|
|
int2str(i_lower).c_str(), f_bounds);
|
|
}
|
|
}
|
|
|
|
doublereal f_delta_bounds = deltaBoundStep(y, step0);
|
|
fbound = std::min(f_bounds, f_delta_bounds);
|
|
|
|
return fbound;
|
|
}
|
|
//===================================================================================================================
|
|
// Find a damping coefficient through a look-ahead mechanism
|
|
/*
|
|
*
|
|
* On entry, step0 must contain an undamped Newton step to the
|
|
* current solution y0. This method attempts to find a damping coefficient
|
|
* such that the next undamped step would have a norm smaller than
|
|
* that of step0. If successful, the new solution after taking the
|
|
* damped step is returned in y1, and the undamped step at y1 is
|
|
* returned in step1.
|
|
*
|
|
*
|
|
* @return 1 Successful step was taken: Next step was less than previous step.
|
|
* s1 is calculated
|
|
* 2 Successful step: Next step's norm is less than 0.8
|
|
* 3 Success: The final residual is less than 1.0
|
|
* A predicted deltaSoln1 is not produced however. s1 is estimated.
|
|
* 4 Success: The final residual is less than the residual
|
|
* from the previous step.
|
|
* A predicted deltaSoln1 is not produced however. s1 is estimated.
|
|
* 0 Uncertain Success: s1 is about the same as s0
|
|
* NSOLN_RETN_FAIL_DAMPSTEP
|
|
* Unsuccessful step. We can not find a damping factor that is suitable.
|
|
*/
|
|
int NonlinearSolver::dampStep(const doublereal time_curr, const doublereal* const y_n_curr,
|
|
const doublereal* const ydot_n_curr, doublereal* const step_1,
|
|
doublereal* const y_n_1, doublereal* const ydot_n_1, doublereal* const step_2,
|
|
doublereal& stepNorm_2, GeneralMatrix& jac, bool writetitle, int& num_backtracks)
|
|
{
|
|
int m;
|
|
int info = 0;
|
|
int retnTrial = NSOLN_RETN_FAIL_DAMPSTEP;
|
|
// Compute the weighted norm of the undamped step size step_1
|
|
doublereal stepNorm_1 = solnErrorNorm(step_1);
|
|
|
|
doublereal* step_1_orig = DATA_PTR(m_wksp);
|
|
for (size_t j = 0; j < neq_; j++) {
|
|
step_1_orig[j] = step_1[j];
|
|
}
|
|
|
|
|
|
// Compute the multiplier to keep all components in bounds.A value of one indicates that there is no limitation
|
|
// on the current step size in the nonlinear method due to bounds constraints (either negative values of delta
|
|
// bounds constraints.
|
|
m_dampBound = boundStep(y_n_curr, step_1);
|
|
|
|
// If fbound is very small, then y0 is already close to the boundary and step0 points out of the allowed domain. In
|
|
// this case, the Newton algorithm fails, so return an error condition.
|
|
if (m_dampBound < 1.e-30) {
|
|
if (m_print_flag > 1) {
|
|
printf("\t\t\tdampStep(): At limits.\n");
|
|
}
|
|
return -3;
|
|
}
|
|
|
|
//--------------------------------------------
|
|
// Attempt damped step
|
|
//--------------------------------------------
|
|
|
|
// damping coefficient starts at 1.0
|
|
m_dampRes = 1.0;
|
|
|
|
doublereal ff = m_dampBound;
|
|
num_backtracks = 0;
|
|
for (m = 0; m < NDAMP; m++) {
|
|
|
|
ff = m_dampBound * m_dampRes;
|
|
|
|
// step the solution by the damped step size
|
|
/*
|
|
* Whenever we update the solution, we must also always
|
|
* update the time derivative.
|
|
*/
|
|
for (size_t j = 0; j < neq_; j++) {
|
|
step_1[j] = ff * step_1_orig[j];
|
|
y_n_1[j] = y_n_curr[j] + step_1[j];
|
|
}
|
|
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
calc_ydot(m_order, y_n_1, ydot_n_1);
|
|
}
|
|
/*
|
|
* Calculate the residual that would result if y1[] were the new solution vector
|
|
* -> m_resid[] contains the result of the residual calculation
|
|
*/
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
info = doResidualCalc(time_curr, solnType_, y_n_1, ydot_n_1, Base_LaggedSolutionComponents);
|
|
} else {
|
|
info = doResidualCalc(time_curr, solnType_, y_n_1, ydot_n_curr, Base_LaggedSolutionComponents);
|
|
}
|
|
if (info != 1) {
|
|
if (m_print_flag > 0) {
|
|
printf("\t\t\tdampStep(): current trial step and damping led to Residual Calc ERROR %d. Bailing\n", info);
|
|
}
|
|
return -1;
|
|
}
|
|
m_normResidTrial = residErrorNorm(DATA_PTR(m_resid));
|
|
m_normResid_1 = m_normResidTrial;
|
|
if (m == 0) {
|
|
m_normResid_Bound = m_normResidTrial;
|
|
}
|
|
|
|
bool steepEnough = (m_normResidTrial < m_normResid_0 * (0.9 * (1.0 - ff) * (1.0 - ff)* (1.0 - ff) + 0.1));
|
|
|
|
if (m_normResidTrial < 1.0 || steepEnough) {
|
|
if (m_print_flag >= 5) {
|
|
if (m_normResidTrial < 1.0) {
|
|
printf("\t dampStep(): Current trial step and damping"
|
|
" coefficient accepted because residTrial test step < 1:\n");
|
|
printf("\t resid0 = %g, residTrial = %g\n", m_normResid_0, m_normResidTrial);
|
|
} else if (steepEnough) {
|
|
printf("\t dampStep(): Current trial step and damping"
|
|
" coefficient accepted because resid0 > residTrial and steep enough:\n");
|
|
printf("\t resid0 = %g, residTrial = %g\n", m_normResid_0, m_normResidTrial);
|
|
} else {
|
|
printf("\t dampStep(): Current trial step and damping"
|
|
" coefficient accepted because residual solution damping is turned off:\n");
|
|
printf("\t resid0 = %g, residTrial = %g\n", m_normResid_0, m_normResidTrial);
|
|
}
|
|
}
|
|
/*
|
|
* We aren't going to solve the system if we don't need to. Therefore, return an estimate
|
|
* of the next solution update based on the ratio of the residual reduction.
|
|
*/
|
|
if (m_normResid_0 > 0.0) {
|
|
stepNorm_2 = stepNorm_1 * m_normResidTrial / m_normResid_0;
|
|
} else {
|
|
stepNorm_2 = 0;
|
|
}
|
|
if (m_normResidTrial < 1.0) {
|
|
retnTrial = 3;
|
|
} else {
|
|
retnTrial = 4;
|
|
}
|
|
break;
|
|
}
|
|
|
|
// Compute the next undamped step, step1[], that would result if y1[] were accepted.
|
|
// We now have two steps that we have calculated step0[] and step1[]
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
info = doNewtonSolve(time_curr, y_n_1, ydot_n_1, step_2, jac);
|
|
} else {
|
|
info = doNewtonSolve(time_curr, y_n_1, ydot_n_curr, step_2, jac);
|
|
}
|
|
if (info) {
|
|
if (m_print_flag > 0) {
|
|
printf("\t\t\tdampStep: current trial step and damping led to LAPACK ERROR %d. Bailing\n", info);
|
|
}
|
|
return -1;
|
|
}
|
|
|
|
// compute the weighted norm of step1
|
|
stepNorm_2 = solnErrorNorm(step_2);
|
|
|
|
// write log information
|
|
if (m_print_flag >= 5) {
|
|
print_solnDelta_norm_contrib((const doublereal*) step_1_orig, "DeltaSoln",
|
|
(const doublereal*) step_2, "DeltaSolnTrial",
|
|
"dampNewt: Important Entries for Weighted Soln Updates:",
|
|
y_n_curr, y_n_1, ff, 5);
|
|
}
|
|
if (m_print_flag >= 4) {
|
|
printf("\t\t\tdampStep(): s1 = %g, s2 = %g, dampBound = %g,"
|
|
"dampRes = %g\n", stepNorm_1, stepNorm_2, m_dampBound, m_dampRes);
|
|
}
|
|
|
|
|
|
// if the norm of s1 is less than the norm of s0, then
|
|
// accept this damping coefficient. Also accept it if this
|
|
// step would result in a converged solution. Otherwise,
|
|
// decrease the damping coefficient and try again.
|
|
|
|
if (stepNorm_2 < 0.8 || stepNorm_2 < stepNorm_1) {
|
|
if (stepNorm_2 < 1.0) {
|
|
if (m_print_flag >= 3) {
|
|
if (stepNorm_2 < 1.0) {
|
|
printf("\t\t\tdampStep: current trial step and damping coefficient accepted because test step < 1\n");
|
|
printf("\t\t\t s2 = %g, s1 = %g\n", stepNorm_2, stepNorm_1);
|
|
}
|
|
}
|
|
retnTrial = 2;
|
|
} else {
|
|
retnTrial = 1;
|
|
}
|
|
break;
|
|
} else {
|
|
if (m_print_flag > 1) {
|
|
printf("\t\t\tdampStep: current step rejected: (s1 = %g > "
|
|
"s0 = %g)", stepNorm_2, stepNorm_1);
|
|
if (m < (NDAMP-1)) {
|
|
printf(" Decreasing damping factor and retrying");
|
|
} else {
|
|
printf(" Giving up!!!");
|
|
}
|
|
printf("\n");
|
|
}
|
|
}
|
|
num_backtracks++;
|
|
m_dampRes /= DampFactor;
|
|
}
|
|
|
|
// If a damping coefficient was found, return 1 if the
|
|
// solution after stepping by the damped step would represent
|
|
// a converged solution, and return 0 otherwise. If no damping
|
|
// coefficient could be found, return NSOLN_RETN_FAIL_DAMPSTEP.
|
|
if (m < NDAMP) {
|
|
if (m_print_flag >= 4) {
|
|
printf("\t dampStep(): current trial step accepted retnTrial = %d, its = %d, damp = %g\n", retnTrial, m+1, ff);
|
|
}
|
|
return retnTrial;
|
|
} else {
|
|
if (stepNorm_2 < 0.5 && (stepNorm_1 < 0.5)) {
|
|
if (m_print_flag >= 4) {
|
|
printf("\t dampStep(): current trial step accepted kindof retnTrial = %d, its = %d, damp = %g\n", 2, m+1, ff);
|
|
}
|
|
return 2;
|
|
}
|
|
if (stepNorm_2 < 1.0) {
|
|
if (m_print_flag >= 4) {
|
|
printf("\t dampStep(): current trial step accepted and soln converged retnTrial ="
|
|
"%d, its = %d, damp = %g\n", 0, m+1, ff);
|
|
}
|
|
return 0;
|
|
}
|
|
}
|
|
if (m_print_flag >= 4) {
|
|
printf("\t dampStep(): current direction is rejected! retnTrial = %d, its = %d, damp = %g\n",
|
|
NSOLN_RETN_FAIL_DAMPSTEP, m+1, ff);
|
|
}
|
|
return NSOLN_RETN_FAIL_DAMPSTEP;
|
|
}
|
|
//====================================================================================================================
|
|
// Damp using the dog leg approach
|
|
/*
|
|
*
|
|
* @param time_curr INPUT Current value of the time
|
|
* @param y_n_curr INPUT Current value of the solution vector
|
|
* @param ydot_n_curr INPUT Current value of the derivative of the solution vector
|
|
* @param step_1 INPUT First trial step for the first iteration
|
|
* @param y_n_1 INPUT First trial value of the solution vector
|
|
* @param ydot_n_1 INPUT First trial value of the derivative of the solution vector
|
|
* @param s1 OUTPUT Norm of the vector step_1
|
|
* @param jac INPUT jacobian
|
|
* @param numTrials OUTPUT number of trials taken in the current damping step
|
|
*
|
|
*
|
|
* @return 1 Success: Good step was taken. The predicted residual norm is less than one
|
|
* 2 Success: Good step: Next step's norm is less than 0.8
|
|
* 3 Success: The final residual is less than 1.0
|
|
* A predicted deltaSoln1 is not produced however. s1 is estimated.
|
|
* 4 Success: The final residual is less than the residual from the previous step.
|
|
* A predicted deltaSoln1 is not produced however. s1 is estimated.
|
|
* 0 Unknown Uncertain Success: s1 is about the same as s0
|
|
* NSOLN_RETN_FAIL_DAMPSTEP
|
|
* Unsuccessful step. Can not find a damping coefficient that is suitable
|
|
*/
|
|
int NonlinearSolver::dampDogLeg(const doublereal time_curr, const doublereal* y_n_curr,
|
|
const doublereal* ydot_n_curr, std::vector<doublereal> & step_1,
|
|
doublereal* const y_n_1, doublereal* const ydot_n_1,
|
|
doublereal& stepNorm_1, doublereal& stepNorm_2, GeneralMatrix& jac, int& numTrials)
|
|
{
|
|
doublereal lambda;
|
|
int info;
|
|
|
|
bool success = false;
|
|
bool haveASuccess = false;
|
|
doublereal trustDeltaOld = trustDelta_;
|
|
doublereal* stepLastGood = DATA_PTR(m_wksp);
|
|
//--------------------------------------------
|
|
// Attempt damped step
|
|
//--------------------------------------------
|
|
|
|
// damping coefficient starts at 1.0
|
|
m_dampRes = 1.0;
|
|
int m;
|
|
doublereal tlen;
|
|
|
|
|
|
for (m = 0; m < NDAMP; m++) {
|
|
numTrials++;
|
|
/*
|
|
* Find the initial value of lambda that satisfies the trust distance, trustDelta_
|
|
*/
|
|
dogLegID_ = calcTrustIntersection(trustDelta_, lambda, dogLegAlpha_);
|
|
if (m_print_flag >= 4) {
|
|
tlen = trustRegionLength();
|
|
printf("\t\t dampDogLeg: trust region with length %13.5E has intersection at leg = %d, alpha = %g, lambda = %g\n",
|
|
tlen, dogLegID_, dogLegAlpha_, lambda);
|
|
}
|
|
/*
|
|
* Figure out the new step vector, step_1, based on (leg, alpha). Here we are using the
|
|
* intersection of the trust oval with the dog-leg curve.
|
|
*/
|
|
fillDogLegStep(dogLegID_, dogLegAlpha_, step_1);
|
|
|
|
/*
|
|
* OK, now that we have step0, Bound the step
|
|
*/
|
|
m_dampBound = boundStep(y_n_curr, DATA_PTR(step_1));
|
|
/*
|
|
* Decrease the step length if we are bound
|
|
*/
|
|
if (m_dampBound < 1.0) {
|
|
for (size_t j = 0; j < neq_; j++) {
|
|
step_1[j] = step_1[j] * m_dampBound;
|
|
}
|
|
}
|
|
/*
|
|
* Calculate the new solution value y1[] given the step size
|
|
*/
|
|
for (size_t j = 0; j < neq_; j++) {
|
|
y_n_1[j] = y_n_curr[j] + step_1[j];
|
|
}
|
|
/*
|
|
* Calculate the new solution time derivative given the step size
|
|
*/
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
calc_ydot(m_order, y_n_1, ydot_n_1);
|
|
}
|
|
/*
|
|
* OK, we have the step0. Now, ask the question whether it satisfies the acceptance criteria
|
|
* as a good step. The overall outcome is returned in the variable info.
|
|
*/
|
|
info = decideStep(time_curr, dogLegID_, dogLegAlpha_, y_n_curr, ydot_n_curr, step_1,
|
|
y_n_1, ydot_n_1, trustDeltaOld);
|
|
m_normResid_Bound = m_normResid_1;
|
|
|
|
/*
|
|
* The algorithm failed to find a solution vector sufficiently different than the current point
|
|
*/
|
|
if (info == -1) {
|
|
|
|
if (m_print_flag >= 1) {
|
|
doublereal stepNorm = solnErrorNorm(DATA_PTR(step_1));
|
|
printf("\t\t dampDogLeg: Current direction rejected, update became too small %g\n", stepNorm);
|
|
success = false;
|
|
break;
|
|
}
|
|
}
|
|
if (info == -2) {
|
|
if (m_print_flag >= 1) {
|
|
printf("\t\t dampDogLeg: current trial step and damping led to LAPACK ERROR %d. Bailing\n", info);
|
|
success = false;
|
|
break;
|
|
}
|
|
}
|
|
if (info == 0) {
|
|
success = true;
|
|
break;
|
|
}
|
|
if (info == 3) {
|
|
|
|
haveASuccess = true;
|
|
// Store the good results in stepLastGood
|
|
mdp::mdp_copy_dbl_1(DATA_PTR(stepLastGood), CONSTD_DATA_PTR(step_1), (int) neq_);
|
|
// Within the program decideStep(), we have already increased the value of trustDelta_. We store the
|
|
// value of step0 in step1, recalculate a larger step0 in the next fillDogLegStep(),
|
|
// and then attempt to see if the larger step works in the next iteration
|
|
}
|
|
if (info == 2) {
|
|
// Step was a failure. If we had a previous success with a smaller stepsize, haveASuccess is true
|
|
// and we execute the next block and break. If we didn't have a previous success, trustDelta_ has
|
|
// already been decreased in the decideStep() routine. We go back and try another iteration with
|
|
// a smaller trust region.
|
|
if (haveASuccess) {
|
|
mdp::mdp_copy_dbl_1(DATA_PTR(step_1), CONSTD_DATA_PTR(stepLastGood), (int) neq_);
|
|
for (size_t j = 0; j < neq_; j++) {
|
|
y_n_1[j] = y_n_curr[j] + step_1[j];
|
|
}
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
calc_ydot(m_order, y_n_1, ydot_n_1);
|
|
}
|
|
success = true;
|
|
break;
|
|
} else {
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Estimate s1, the norm after the next step
|
|
*/
|
|
stepNorm_1 = solnErrorNorm(DATA_PTR(step_1));
|
|
stepNorm_2 = stepNorm_1;
|
|
if (m_dampBound < 1.0) {
|
|
stepNorm_2 /= m_dampBound;
|
|
}
|
|
stepNorm_2 /= lambda;
|
|
stepNorm_2 *= m_normResidTrial / m_normResid_0;
|
|
|
|
|
|
if (success) {
|
|
if (m_normResidTrial < 1.0) {
|
|
if (normTrust_Newton_ < trustDelta_ && m_dampBound == 1.0) {
|
|
return 1;
|
|
} else {
|
|
return 0;
|
|
}
|
|
}
|
|
return 0;
|
|
}
|
|
return NSOLN_RETN_FAIL_DAMPSTEP;
|
|
}
|
|
//====================================================================================================================
|
|
// Decide whether the current step is acceptable and adjust the trust region size
|
|
/*
|
|
* This is an extension of algorithm 6.4.5 of Dennis and Schnabel.
|
|
*
|
|
* Here we decide whether to accept the current step
|
|
* At the end of the calculation a new estimate of the trust region is calculated
|
|
*
|
|
* @param time_curr INPUT Current value of the time
|
|
* @param leg INPUT Leg of the dogleg that we are on
|
|
* @param alpha INPUT Distance down that leg that we are on
|
|
* @param y0 INPUT Current value of the solution vector
|
|
* @param ydot0 INPUT Current value of the derivative of the solution vector
|
|
* @param step0 INPUT Trial step
|
|
* @param y1 OUTPUT Solution values at the conditions which are evaluated for success
|
|
* @param ydot1 OUTPUT Time derivates of solution at the conditions which are evaluated for success
|
|
* @param trustDeltaOld INPUT Value of the trust length at the old conditions
|
|
*
|
|
*
|
|
* @return This function returns a code which indicates whether the step will be accepted or not.
|
|
* 3 Step passed with flying colors. Try redoing the calculation with a bigger trust region.
|
|
* 2 Step didn't pass deltaF requirement. Decrease the size of the next trust region for a retry and return
|
|
* 0 The step passed.
|
|
* -1 The step size is now too small (||d || < 0.1). A really small step isn't decreasing the function.
|
|
* This is an error condition.
|
|
* -2 Current value of the solution vector caused a residual error in its evaluation.
|
|
* Step is a failure, and the step size must be reduced in order to proceed further.
|
|
*/
|
|
int NonlinearSolver::decideStep(const doublereal time_curr, int leg, doublereal alpha,
|
|
const doublereal* const y_n_curr,
|
|
const doublereal* const ydot_n_curr, const std::vector<doublereal> & step_1,
|
|
const doublereal* const y_n_1, const doublereal* const ydot_n_1,
|
|
doublereal trustDeltaOld)
|
|
{
|
|
int retn = 2;
|
|
int info;
|
|
doublereal ll;
|
|
// Calculate the solution step length
|
|
doublereal stepNorm = solnErrorNorm(DATA_PTR(step_1));
|
|
|
|
// Calculate the initial (R**2 * neq) value for the old function
|
|
doublereal normResid0_2 = m_normResid_0 * m_normResid_0 * neq_;
|
|
|
|
// Calculate the distance to the cauchy point
|
|
doublereal cauchyDistanceNorm = solnErrorNorm(DATA_PTR(deltaX_CP_));
|
|
|
|
// This is the expected initial rate of decrease in the cauchy direction.
|
|
// -> This is Eqn. 29 = Rhat dot Jhat dy / || d ||
|
|
doublereal funcDecreaseSDExp = RJd_norm_ / cauchyDistanceNorm * lambdaStar_;
|
|
if (funcDecreaseSDExp > 0.0) {
|
|
if (m_print_flag >= 5) {
|
|
printf("\t\tdecideStep(): Unexpected condition -> cauchy slope is positive\n");
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Calculate the residual that would result if y1[] were the new solution vector.
|
|
* The Lagged solution components are kept lagged here. Unfortunately, it just doesn't work in some cases to use a
|
|
* Jacobian from a lagged state and then use a residual from an unlagged condition. The linear model doesn't
|
|
* agree with the nonlinear model.
|
|
* -> m_resid[] contains the result of the residual calculation
|
|
*/
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
info = doResidualCalc(time_curr, solnType_, y_n_1, ydot_n_1, Base_LaggedSolutionComponents);
|
|
} else {
|
|
info = doResidualCalc(time_curr, solnType_, y_n_1, ydot_n_curr, Base_LaggedSolutionComponents);
|
|
}
|
|
|
|
if (info != 1) {
|
|
if (m_print_flag >= 2) {
|
|
printf("\t\tdecideStep: current trial step and damping led to Residual Calc ERROR %d. Bailing\n", info);
|
|
}
|
|
return -2;
|
|
}
|
|
/*
|
|
* Ok we have a successful new residual. Calculate the normalized residual value and store it in
|
|
* m_normResidTrial
|
|
*/
|
|
m_normResidTrial = residErrorNorm(DATA_PTR(m_resid));
|
|
doublereal normResidTrial_2 = neq_ * m_normResidTrial * m_normResidTrial;
|
|
|
|
/*
|
|
* We have a minimal acceptance test for passage. deltaf < 1.0E-4 (CauchySlope) (deltS)
|
|
* This is the condition that D&S use in 6.4.5
|
|
*/
|
|
doublereal funcDecrease = 0.5 * (normResidTrial_2 - normResid0_2);
|
|
doublereal acceptableDelF = funcDecreaseSDExp * stepNorm * 1.0E-4;
|
|
if (funcDecrease < acceptableDelF) {
|
|
m_normResid_1 = m_normResidTrial;
|
|
m_normResid_1 = m_normResidTrial;
|
|
retn = 0;
|
|
if (m_print_flag >= 4) {
|
|
printf("\t\t decideStep: Norm Residual(leg=%1d, alpha=%10.2E) = %11.4E passes\n",
|
|
dogLegID_, dogLegAlpha_, m_normResidTrial);
|
|
}
|
|
} else {
|
|
if (m_print_flag >= 4) {
|
|
printf("\t\t decideStep: Norm Residual(leg=%1d, alpha=%10.2E) = %11.4E failes\n",
|
|
dogLegID_, dogLegAlpha_, m_normResidTrial);
|
|
}
|
|
trustDelta_ *= 0.33;
|
|
CurrentTrustFactor_ *= 0.33;
|
|
retn = 2;
|
|
// error condition if step is getting too small
|
|
if (rtol_ * stepNorm < 1.0E-6) {
|
|
retn = -1;
|
|
}
|
|
return retn;
|
|
}
|
|
/*
|
|
* Figure out the next trust region. We are here iff retn = 0
|
|
*
|
|
* If we had to bounds delta the update, decrease the trust region
|
|
*/
|
|
if (m_dampBound < 1.0) {
|
|
// trustDelta_ *= 0.5;
|
|
// NextTrustFactor_ *= 0.5;
|
|
// ll = trustRegionLength();
|
|
// if (m_print_flag >= 5) {
|
|
// printf("\t\tdecideStep(): Trust region decreased from %g to %g due to bounds constraint\n", ll*2, ll);
|
|
//}
|
|
} else {
|
|
retn = 0;
|
|
/*
|
|
* Calculate the expected residual from the quadratic model
|
|
*/
|
|
doublereal expectedNormRes = expectedResidLeg(leg, alpha);
|
|
doublereal expectedFuncDecrease = 0.5 * (neq_ * expectedNormRes * expectedNormRes - normResid0_2);
|
|
if (funcDecrease > 0.1 * expectedFuncDecrease) {
|
|
if ((m_normResidTrial > 0.5 * m_normResid_0) && (m_normResidTrial > 0.1)) {
|
|
trustDelta_ *= 0.5;
|
|
NextTrustFactor_ *= 0.5;
|
|
ll = trustRegionLength();
|
|
if (m_print_flag >= 4) {
|
|
printf("\t\t decideStep: Trust region decreased from %g to %g due to bad quad approximation\n",
|
|
ll*2, ll);
|
|
}
|
|
}
|
|
} else {
|
|
/*
|
|
* If we are doing well, consider increasing the trust region and recalculating
|
|
*/
|
|
if (funcDecrease < 0.8 * expectedFuncDecrease || (m_normResidTrial < 0.33 * m_normResid_0)) {
|
|
if (trustDelta_ <= trustDeltaOld && (leg != 2 || alpha < 0.75)) {
|
|
trustDelta_ *= 2.0;
|
|
CurrentTrustFactor_ *= 2;
|
|
adjustUpStepMinimums();
|
|
ll = trustRegionLength();
|
|
if (m_print_flag >= 4) {
|
|
if (m_normResidTrial < 0.33 * m_normResid_0) {
|
|
printf("\t\t decideStep: Redo line search with trust region increased from %g to %g due to good nonlinear behavior\n",
|
|
ll*0.5, ll);
|
|
} else {
|
|
printf("\t\t decideStep: Redi line search with trust region increased from %g to %g due to good linear model approximation\n",
|
|
ll*0.5, ll);
|
|
}
|
|
}
|
|
retn = 3;
|
|
} else {
|
|
/*
|
|
* Increase the size of the trust region for the next calculation
|
|
*/
|
|
if (m_normResidTrial < 0.99 * expectedNormRes || (m_normResidTrial < 0.20 * m_normResid_0) ||
|
|
(funcDecrease < -1.0E-50 && (funcDecrease < 0.9 *expectedFuncDecrease))) {
|
|
if (leg == 2 && alpha == 1.0) {
|
|
ll = trustRegionLength();
|
|
if (ll < 2.0 * m_normDeltaSoln_Newton) {
|
|
trustDelta_ *= 2.0;
|
|
NextTrustFactor_ *= 2.0;
|
|
adjustUpStepMinimums();
|
|
ll = trustRegionLength();
|
|
if (m_print_flag >= 4) {
|
|
printf("\t\t decideStep: Trust region further increased from %g to %g next step due to good linear model behavior\n",
|
|
ll*0.5, ll);
|
|
}
|
|
}
|
|
} else {
|
|
ll = trustRegionLength();
|
|
trustDelta_ *= 2.0;
|
|
NextTrustFactor_ *= 2.0;
|
|
adjustUpStepMinimums();
|
|
ll = trustRegionLength();
|
|
if (m_print_flag >= 4) {
|
|
printf("\t\t decideStep: Trust region further increased from %g to %g next step due to good linear model behavior\n",
|
|
ll*0.5, ll);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return retn;
|
|
}
|
|
//====================================================================================================================
|
|
/*
|
|
* solve_nonlinear_problem():
|
|
*
|
|
* Find the solution to F(X) = 0 by damped Newton iteration. On
|
|
* entry, x0 contains an initial estimate of the solution. On
|
|
* successful return, x1 contains the converged solution.
|
|
*
|
|
* SolnType = TRANSIENT -> we will assume we are relaxing a transient
|
|
* equation system for now. Will make it more general later,
|
|
* if an application comes up.
|
|
*
|
|
* @return A positive value indicates a successful convergence
|
|
* -1 Failed convergence
|
|
*/
|
|
int NonlinearSolver::solve_nonlinear_problem(int SolnType, doublereal* const y_comm, doublereal* const ydot_comm,
|
|
doublereal CJ, doublereal time_curr, GeneralMatrix& jac,
|
|
int& num_newt_its, int& num_linear_solves,
|
|
int& num_backtracks, int loglevelInput)
|
|
{
|
|
clockWC wc;
|
|
int convRes = 0;
|
|
solnType_ = SolnType;
|
|
int info = 0;
|
|
|
|
num_linear_solves -= m_numTotalLinearSolves;
|
|
int retnDamp = 0;
|
|
int retnCode = 0;
|
|
bool forceNewJac = false;
|
|
|
|
if (jacCopyPtr_) {
|
|
delete jacCopyPtr_;
|
|
}
|
|
jacCopyPtr_ = jac.duplMyselfAsGeneralMatrix();
|
|
|
|
doublereal stepNorm_1;
|
|
doublereal stepNorm_2;
|
|
#ifdef DEBUG_MODE
|
|
int legBest;
|
|
doublereal alphaBest;
|
|
#endif
|
|
bool trInit = false;
|
|
|
|
|
|
mdp::mdp_copy_dbl_1(DATA_PTR(m_y_n_curr), DATA_PTR(y_comm), (int) neq_);
|
|
|
|
if (SolnType != NSOLN_TYPE_STEADY_STATE || ydot_comm) {
|
|
mdp::mdp_copy_dbl_1(DATA_PTR(m_ydot_n_curr), ydot_comm, (int) neq_);
|
|
mdp::mdp_copy_dbl_1(DATA_PTR(m_ydot_n_1), ydot_comm, (int) neq_);
|
|
}
|
|
// Redo the solution weights every time we enter the function
|
|
createSolnWeights(DATA_PTR(m_y_n_curr));
|
|
m_normDeltaSoln_Newton = 1.0E1;
|
|
bool frst = true;
|
|
num_newt_its = 0;
|
|
num_backtracks = 0;
|
|
int i_numTrials;
|
|
m_print_flag = loglevelInput;
|
|
|
|
if (trustRegionInitializationMethod_ == 0) {
|
|
trInit = true;
|
|
} else if (trustRegionInitializationMethod_ == 1) {
|
|
trInit = true;
|
|
initializeTrustRegion();
|
|
} else {
|
|
mdp::mdp_init_dbl_1(DATA_PTR(deltaX_trust_), 1.0, (int) neq_);
|
|
trustDelta_ = 1.0;
|
|
}
|
|
|
|
if (m_print_flag == 2 || m_print_flag == 3) {
|
|
printf("\tsolve_nonlinear_problem():\n\n");
|
|
if (doDogLeg_) {
|
|
printf("\tWt Iter Resid NewJac log(CN)| dRdS_CDexp dRdS_CD dRdS_Newtexp dRdS_Newt |"
|
|
"DS_Cauchy DS_Newton DS_Trust | legID legAlpha Fbound | CTF NTF | nTr|"
|
|
"DS_Final ResidLag ResidFull\n");
|
|
printf("\t---------------------------------------------------------------------------------------------------"
|
|
"--------------------------------------------------------------------------------\n");
|
|
} else {
|
|
printf("\t Wt Iter Resid NewJac | Fbound ResidBound | DampIts Fdamp DS_Step1 DS_Step2"
|
|
"ResidLag | DS_Damp DS_Newton ResidFull\n");
|
|
printf("\t--------------------------------------------------------------------------------------------------"
|
|
"----------------------------------\n");
|
|
}
|
|
}
|
|
|
|
while (1 > 0) {
|
|
|
|
CurrentTrustFactor_ = 1.0;
|
|
NextTrustFactor_ = 1.0;
|
|
ResidWtsReevaluated_ = false;
|
|
i_numTrials = 0;
|
|
/*
|
|
* Increment Newton Solve counter
|
|
*/
|
|
m_numTotalNewtIts++;
|
|
num_newt_its++;
|
|
m_numLocalLinearSolves = 0;
|
|
|
|
if (m_print_flag > 3) {
|
|
printf("\t");
|
|
print_line("=", 119);
|
|
printf("\tsolve_nonlinear_problem(): iteration %d:\n",
|
|
num_newt_its);
|
|
}
|
|
/*
|
|
* If we are far enough away from the solution, redo the solution weights and the trust vectors.
|
|
*/
|
|
if (m_normDeltaSoln_Newton > 1.0E2) {
|
|
createSolnWeights(DATA_PTR(m_y_n_curr));
|
|
#ifdef DEBUG_MODE
|
|
if (trInit) {
|
|
readjustTrustVector();
|
|
}
|
|
#else
|
|
if (doDogLeg_ && trInit) {
|
|
readjustTrustVector();
|
|
}
|
|
#endif
|
|
} else {
|
|
// Do this stuff every 5 iterations
|
|
if ((num_newt_its % 5) == 1) {
|
|
createSolnWeights(DATA_PTR(m_y_n_curr));
|
|
#ifdef DEBUG_MODE
|
|
if (trInit) {
|
|
readjustTrustVector();
|
|
}
|
|
#else
|
|
if (doDogLeg_ && trInit) {
|
|
readjustTrustVector();
|
|
}
|
|
#endif
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Set default values of Delta bounds constraints
|
|
*/
|
|
if (!m_manualDeltaStepSet) {
|
|
setDefaultDeltaBoundsMagnitudes();
|
|
}
|
|
|
|
// Check whether the Jacobian should be re-evaluated.
|
|
|
|
forceNewJac = true;
|
|
|
|
if (forceNewJac) {
|
|
if (m_print_flag > 3) {
|
|
printf("\t solve_nonlinear_problem(): Getting a new Jacobian\n");
|
|
}
|
|
info = beuler_jac(jac, DATA_PTR(m_resid), time_curr, CJ, DATA_PTR(m_y_n_curr),
|
|
DATA_PTR(m_ydot_n_curr), num_newt_its);
|
|
if (info != 1) {
|
|
if (m_print_flag > 0) {
|
|
printf("\t solve_nonlinear_problem(): Jacobian Formation Error: %d Bailing\n", info);
|
|
}
|
|
retnDamp = NSOLN_RETN_JACOBIANFORMATIONERROR ;
|
|
goto done;
|
|
}
|
|
} else {
|
|
if (m_print_flag > 1) {
|
|
printf("\t solve_nonlinear_problem(): Solving system with old jacobian\n");
|
|
}
|
|
}
|
|
/*
|
|
* Go get new scales
|
|
*/
|
|
calcColumnScales();
|
|
|
|
|
|
/*
|
|
* Calculate the base residual
|
|
*/
|
|
if (m_print_flag >= 6) {
|
|
printf("\t solve_nonlinear_problem(): Calculate the base residual\n");
|
|
}
|
|
info = doResidualCalc(time_curr, NSOLN_TYPE_STEADY_STATE, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr));
|
|
if (info != 1) {
|
|
if (m_print_flag > 0) {
|
|
printf("\t solve_nonlinear_problem(): Residual Calc ERROR %d. Bailing\n", info);
|
|
}
|
|
retnDamp = NSOLN_RETN_RESIDUALFORMATIONERROR;
|
|
goto done;
|
|
}
|
|
|
|
/*
|
|
* Scale the matrix and the rhs, if they aren't already scaled
|
|
* Figure out and store the residual scaling factors.
|
|
*/
|
|
scaleMatrix(jac, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr), time_curr, num_newt_its);
|
|
|
|
|
|
/*
|
|
* Optional print out the initial residual
|
|
*/
|
|
if (m_print_flag >= 6) {
|
|
m_normResid_0 = residErrorNorm(DATA_PTR(m_resid), "Initial norm of the residual", 10, DATA_PTR(m_y_n_curr));
|
|
} else {
|
|
m_normResid_0 = residErrorNorm(DATA_PTR(m_resid), "Initial norm of the residual", 0, DATA_PTR(m_y_n_curr));
|
|
if (m_print_flag == 4 || m_print_flag == 5) {
|
|
printf("\t solve_nonlinear_problem(): Initial Residual Norm = %13.4E\n", m_normResid_0);
|
|
}
|
|
}
|
|
|
|
|
|
#ifdef DEBUG_MODE
|
|
if (m_print_flag > 3) {
|
|
printf("\t solve_nonlinear_problem(): Calculate the steepest descent direction and Cauchy Point\n");
|
|
}
|
|
m_normDeltaSoln_CP = doCauchyPointSolve(jac);
|
|
|
|
#else
|
|
if (doDogLeg_) {
|
|
if (m_print_flag > 3) {
|
|
printf("\t solve_nonlinear_problem(): Calculate the steepest descent direction and Cauchy Point\n");
|
|
}
|
|
m_normDeltaSoln_CP = doCauchyPointSolve(jac);
|
|
}
|
|
#endif
|
|
|
|
// compute the undamped Newton step
|
|
if (doAffineSolve_) {
|
|
if (m_print_flag >= 4) {
|
|
printf("\t solve_nonlinear_problem(): Calculate the Newton direction via an Affine solve\n");
|
|
}
|
|
info = doAffineNewtonSolve(DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr), DATA_PTR(deltaX_Newton_), jac);
|
|
} else {
|
|
if (m_print_flag >= 4) {
|
|
printf("\t solve_nonlinear_problem(): Calculate the Newton direction via a Newton solve\n");
|
|
}
|
|
info = doNewtonSolve(time_curr, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr), DATA_PTR(deltaX_Newton_), jac);
|
|
}
|
|
|
|
if (info) {
|
|
retnDamp = NSOLN_RETN_MATRIXINVERSIONERROR;
|
|
if (m_print_flag > 0) {
|
|
printf("\t solve_nonlinear_problem(): Matrix Inversion Error: %d Bailing\n", info);
|
|
}
|
|
goto done;
|
|
}
|
|
mdp::mdp_copy_dbl_1(DATA_PTR(m_step_1), CONSTD_DATA_PTR(deltaX_Newton_), (int) neq_);
|
|
|
|
if (m_print_flag >= 6) {
|
|
m_normDeltaSoln_Newton = solnErrorNorm(DATA_PTR(deltaX_Newton_), "Initial Undamped Newton Step of the iteration", 10);
|
|
} else {
|
|
m_normDeltaSoln_Newton = solnErrorNorm(DATA_PTR(deltaX_Newton_), "Initial Undamped Newton Step of the iteration", 0);
|
|
}
|
|
|
|
if (m_numTotalNewtIts == 1) {
|
|
if (trustRegionInitializationMethod_ == 2 || trustRegionInitializationMethod_ == 3) {
|
|
if (m_print_flag > 3) {
|
|
if (trustRegionInitializationMethod_ == 2) {
|
|
printf("\t solve_nonlinear_problem(): Initialize the trust region size as the length of the Cauchy Vector times %f\n",
|
|
trustRegionInitializationFactor_);
|
|
} else {
|
|
printf("\t solve_nonlinear_problem(): Initialize the trust region size as the length of the Newton Vector times %f\n",
|
|
trustRegionInitializationFactor_);
|
|
}
|
|
}
|
|
initializeTrustRegion();
|
|
trInit = true;
|
|
}
|
|
}
|
|
|
|
|
|
if (doDogLeg_) {
|
|
|
|
|
|
|
|
#ifdef DEBUG_MODE
|
|
doublereal trustD = calcTrustDistance(m_step_1);
|
|
if (m_print_flag >= 4) {
|
|
if (trustD > trustDelta_) {
|
|
printf("\t\t Newton's method step size, %g trustVectorUnits, larger than trust region, %g trustVectorUnits\n",
|
|
trustD, trustDelta_);
|
|
printf("\t\t Newton's method step size, %g trustVectorUnits, larger than trust region, %g trustVectorUnits\n",
|
|
trustD, trustDelta_);
|
|
} else {
|
|
printf("\t\t Newton's method step size, %g trustVectorUnits, smaller than trust region, %g trustVectorUnits\n",
|
|
trustD, trustDelta_);
|
|
}
|
|
}
|
|
#endif
|
|
}
|
|
|
|
/*
|
|
* Filter out bad directions
|
|
*/
|
|
filterNewStep(time_curr, DATA_PTR(m_y_n_curr), DATA_PTR(m_step_1));
|
|
|
|
|
|
|
|
if (s_print_DogLeg && m_print_flag >= 4) {
|
|
printf("\t solve_nonlinear_problem(): Compare descent rates for Cauchy and Newton directions\n");
|
|
descentComparison(time_curr, DATA_PTR(m_ydot_n_curr), DATA_PTR(m_ydot_n_1), i_numTrials);
|
|
} else {
|
|
if (doDogLeg_) {
|
|
descentComparison(time_curr, DATA_PTR(m_ydot_n_curr), DATA_PTR(m_ydot_n_1), i_numTrials);
|
|
}
|
|
}
|
|
|
|
|
|
|
|
if (doDogLeg_) {
|
|
setupDoubleDogleg();
|
|
#ifdef DEBUG_MODE
|
|
if (s_print_DogLeg && m_print_flag >= 5) {
|
|
printf("\t solve_nonlinear_problem(): Compare Linear and nonlinear residuals along double dog-leg path\n");
|
|
residualComparisonLeg(time_curr, DATA_PTR(m_ydot_n_curr), legBest, alphaBest);
|
|
}
|
|
#endif
|
|
if (m_print_flag >= 4) {
|
|
printf("\t solve_nonlinear_problem(): Calculate damping along dog-leg path to ensure residual decrease\n");
|
|
}
|
|
retnDamp = dampDogLeg(time_curr, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr),
|
|
m_step_1, DATA_PTR(m_y_n_1), DATA_PTR(m_ydot_n_1), stepNorm_1, stepNorm_2, jac, i_numTrials);
|
|
}
|
|
#ifdef DEBUG_MODE
|
|
else {
|
|
if (s_print_DogLeg && m_print_flag >= 5) {
|
|
printf("\t solve_nonlinear_problem(): Compare Linear and nonlinear residuals along double dog-leg path\n");
|
|
residualComparisonLeg(time_curr, DATA_PTR(m_ydot_n_curr), legBest, alphaBest);
|
|
}
|
|
}
|
|
#endif
|
|
|
|
// Damp the Newton step
|
|
/*
|
|
* On return the recommended new solution and derivatisve is located in:
|
|
* y_new
|
|
* y_dot_new
|
|
* The update delta vector is located in
|
|
* stp1
|
|
* The estimate of the solution update norm for the next step is located in
|
|
* s1
|
|
*/
|
|
if (!doDogLeg_) {
|
|
retnDamp = dampStep(time_curr, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr),
|
|
DATA_PTR(m_step_1), DATA_PTR(m_y_n_1), DATA_PTR(m_ydot_n_1),
|
|
DATA_PTR(m_wksp_2), stepNorm_2, jac, frst, i_numTrials);
|
|
frst = false;
|
|
num_backtracks += i_numTrials;
|
|
stepNorm_1 = solnErrorNorm(DATA_PTR(m_step_1));
|
|
}
|
|
|
|
|
|
/*
|
|
* Impose the minimum number of newton iterations critera
|
|
*/
|
|
if (num_newt_its < m_min_newt_its) {
|
|
if (retnDamp > NSOLN_RETN_CONTINUE) {
|
|
if (m_print_flag > 2) {
|
|
printf("\t solve_nonlinear_problem(): Damped Newton successful (m=%d) but minimum newton"
|
|
"iterations not attained. Resolving ...\n", retnDamp);
|
|
}
|
|
retnDamp = NSOLN_RETN_CONTINUE;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Impose max newton iteration
|
|
*/
|
|
if (num_newt_its > maxNewtIts_) {
|
|
retnDamp = NSOLN_RETN_MAXIMUMITERATIONSEXCEEDED;
|
|
if (m_print_flag > 1) {
|
|
printf("\t solve_nonlinear_problem(): Damped newton unsuccessful (max newts exceeded) sfinal = %g\n",
|
|
stepNorm_1);
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Do a full residual calculation with the unlagged solution components.
|
|
* Then get the norm of the residual
|
|
*/
|
|
info = doResidualCalc(time_curr, NSOLN_TYPE_STEADY_STATE, DATA_PTR(m_y_n_1), DATA_PTR(m_ydot_n_1));
|
|
if (info != 1) {
|
|
if (m_print_flag > 0) {
|
|
printf("\t solve_nonlinear_problem(): current trial step and damping led to Residual Calc "
|
|
"ERROR %d. Bailing\n", info);
|
|
}
|
|
retnDamp = NSOLN_RETN_RESIDUALFORMATIONERROR;
|
|
goto done;
|
|
}
|
|
if (m_print_flag >= 4) {
|
|
m_normResid_full = residErrorNorm(DATA_PTR(m_resid), " Resulting full residual norm", 10, DATA_PTR(m_y_n_1));
|
|
if (fabs(m_normResid_full - m_normResid_1) > 1.0E-3 * (m_normResid_1 + m_normResid_full + 1.0E-4)) {
|
|
if (m_print_flag >= 4) {
|
|
printf("\t solve_nonlinear_problem(): Full residual norm changed from %g to %g due to "
|
|
"lagging of components\n", m_normResid_1, m_normResid_full);
|
|
}
|
|
}
|
|
} else {
|
|
m_normResid_full = residErrorNorm(DATA_PTR(m_resid));
|
|
}
|
|
|
|
/*
|
|
* Check the convergence criteria
|
|
*/
|
|
convRes = 0;
|
|
if (retnDamp > NSOLN_RETN_CONTINUE) {
|
|
convRes = convergenceCheck(retnDamp, stepNorm_1);
|
|
}
|
|
|
|
|
|
|
|
|
|
bool m_filterIntermediate = false;
|
|
if (m_filterIntermediate) {
|
|
if (retnDamp == NSOLN_RETN_CONTINUE) {
|
|
(void) filterNewSolution(time_n, DATA_PTR(m_y_n_1), DATA_PTR(m_ydot_n_1));
|
|
}
|
|
}
|
|
|
|
// Exchange new for curr solutions
|
|
if (retnDamp >= NSOLN_RETN_CONTINUE) {
|
|
mdp::mdp_copy_dbl_1(DATA_PTR(m_y_n_curr), CONSTD_DATA_PTR(m_y_n_1), (int) neq_);
|
|
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
calc_ydot(m_order, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr));
|
|
}
|
|
}
|
|
|
|
if (m_print_flag == 2 || m_print_flag == 3) {
|
|
// printf("\t Iter Resid NewJac | Fbound | ResidBound | Fdamp DampIts | DeltaSolnNewton ResidFinal \n");
|
|
if (ResidWtsReevaluated_) {
|
|
printf("\t*");
|
|
} else {
|
|
printf("\t ");
|
|
}
|
|
printf(" %3d %11.3E", num_newt_its, m_normResid_0);
|
|
bool m_jacAge = false;
|
|
if (!m_jacAge) {
|
|
printf(" Y ");
|
|
} else {
|
|
printf(" N ");
|
|
}
|
|
if (doDogLeg_) {
|
|
printf("%5.1f |", log10(m_conditionNumber));
|
|
// printf("\t Iter Resid NewJac | DS_Cauchy DS_Newton DS_Trust | legID legAlpha Fbound | | DS_F ResidFinal \n");
|
|
printf("%10.3E %10.3E %10.3E %10.3E|", ResidDecreaseSDExp_, ResidDecreaseSD_,
|
|
ResidDecreaseNewtExp_, ResidDecreaseNewt_);
|
|
printf("%10.3E %10.3E %10.3E|", m_normDeltaSoln_CP , m_normDeltaSoln_Newton, norm_deltaX_trust_ * trustDelta_);
|
|
printf("%2d %10.2E %10.2E", dogLegID_ , dogLegAlpha_, m_dampBound);
|
|
printf("| %3.2f %3.2f |", CurrentTrustFactor_, NextTrustFactor_);
|
|
printf(" %2d ", i_numTrials);
|
|
printf("| %10.3E %10.3E %10.3E", stepNorm_1, m_normResid_1, m_normResid_full);
|
|
} else {
|
|
printf(" |");
|
|
printf("%10.2E %10.3E |", m_dampBound, m_normResid_Bound);
|
|
printf("%2d %10.2E %10.3E %10.3E %10.3E", i_numTrials + 1, m_dampRes,
|
|
stepNorm_1 / (m_dampRes * m_dampBound), stepNorm_2, m_normResid_1);
|
|
printf("| %10.3E %10.3E %10.3E", stepNorm_1, m_normDeltaSoln_Newton, m_normResid_full);
|
|
}
|
|
printf("\n");
|
|
|
|
}
|
|
if (m_print_flag >= 4) {
|
|
if (doDogLeg_) {
|
|
if (convRes > 0) {
|
|
printf("\t solve_nonlinear_problem(): Problem Converged, stepNorm = %11.3E, reduction of res from %11.3E to %11.3E\n",
|
|
stepNorm_1, m_normResid_0, m_normResid_full);
|
|
printf("\t");
|
|
print_line("=", 119);
|
|
} else {
|
|
printf("\t solve_nonlinear_problem(): Successfull step taken with stepNorm = %11.3E, reduction of res from %11.3E to %11.3E\n",
|
|
stepNorm_1, m_normResid_0, m_normResid_full);
|
|
}
|
|
} else {
|
|
if (convRes > 0) {
|
|
printf("\t solve_nonlinear_problem(): Damped Newton iteration successful, nonlin "
|
|
"converged, final estimate of the next solution update norm = %-12.4E\n", stepNorm_2);
|
|
printf("\t");
|
|
print_line("=", 119);
|
|
} else if (retnDamp >= NSOLN_RETN_CONTINUE) {
|
|
printf("\t solve_nonlinear_problem(): Damped Newton iteration successful, "
|
|
"estimate of the next solution update norm = %-12.4E\n", stepNorm_2);
|
|
} else {
|
|
printf("\t solve_nonlinear_problem(): Damped Newton unsuccessful, final estimate "
|
|
"of the next solution update norm = %-12.4E\n", stepNorm_2);
|
|
}
|
|
}
|
|
}
|
|
// convergence
|
|
if (convRes) {
|
|
goto done;
|
|
}
|
|
|
|
// If dampStep fails, first try a new Jacobian if an old
|
|
// one was being used. If it was a new Jacobian, then
|
|
// return -1 to signify failure.
|
|
else if (retnDamp < NSOLN_RETN_CONTINUE) {
|
|
goto done;
|
|
}
|
|
}
|
|
|
|
done:
|
|
|
|
|
|
if (m_print_flag == 2 || m_print_flag == 3) {
|
|
if (convRes > 0) {
|
|
if (doDogLeg_) {
|
|
if (convRes == 3) {
|
|
printf("\t | | "
|
|
" | | converged = 3 |(%11.3E) \n", stepNorm_2);
|
|
} else {
|
|
printf("\t | | "
|
|
" | | converged = %1d | %10.3E %10.3E\n", convRes,
|
|
stepNorm_2, m_normResidTrial);
|
|
}
|
|
printf("\t-----------------------------------------------------------------------------------------------------"
|
|
"------------------------------------------------------------------------------\n");
|
|
} else {
|
|
if (convRes == 3) {
|
|
printf("\t | "
|
|
" | converged = 3 | (%11.3E) \n", stepNorm_2);
|
|
} else {
|
|
printf("\t | "
|
|
" | converged = %1d | %10.3E %10.3E\n", convRes,
|
|
stepNorm_2, m_normResidTrial);
|
|
}
|
|
printf("\t------------------------------------------------------------------------------------"
|
|
"-----------------------------------------------\n");
|
|
}
|
|
}
|
|
|
|
|
|
|
|
}
|
|
|
|
mdp::mdp_copy_dbl_1(y_comm, CONSTD_DATA_PTR(m_y_n_curr), (int) neq_);
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
mdp::mdp_copy_dbl_1(ydot_comm, CONSTD_DATA_PTR(m_ydot_n_curr), (int) neq_);
|
|
}
|
|
|
|
num_linear_solves += m_numTotalLinearSolves;
|
|
|
|
doublereal time_elapsed = wc.secondsWC();
|
|
if (m_print_flag > 1) {
|
|
if (retnDamp > 0) {
|
|
if (NonlinearSolver::s_TurnOffTiming) {
|
|
printf("\tNonlinear problem solved successfully in %d its\n",
|
|
num_newt_its);
|
|
} else {
|
|
printf("\tNonlinear problem solved successfully in %d its, time elapsed = %g sec\n",
|
|
num_newt_its, time_elapsed);
|
|
}
|
|
} else {
|
|
printf("\tNonlinear problem failed to solve after %d its\n", num_newt_its);
|
|
}
|
|
}
|
|
retnCode = retnDamp;
|
|
if (retnDamp > 0) {
|
|
retnCode = NSOLN_RETN_SUCCESS;
|
|
}
|
|
|
|
|
|
return retnCode;
|
|
}
|
|
//====================================================================================================================
|
|
// Print solution norm contribution
|
|
/*
|
|
* Prints out the most important entries to the update to the solution vector for the current step
|
|
*
|
|
* @param step_1 Raw update vector for the current nonlinear step
|
|
* @param stepNorm_1 Norm of the vector step_1
|
|
* @param step_2 Raw update vector for the next solution value based on the old matrix
|
|
* @param stepNorm_2 Norm of the vector step_2
|
|
* @param title title of the printout
|
|
* @param y_n_curr Old value of the solution
|
|
* @param y_n_1 New value of the solution after damping corrections
|
|
* @param damp Value of the damping factor
|
|
* @param num_entries Number of entries to print out
|
|
*/
|
|
void NonlinearSolver::
|
|
print_solnDelta_norm_contrib(const doublereal* const step_1,
|
|
const char* const stepNorm_1,
|
|
const doublereal* const step_2,
|
|
const char* const stepNorm_2,
|
|
const char* const title,
|
|
const doublereal* const y_n_curr,
|
|
const doublereal* const y_n_1,
|
|
doublereal damp,
|
|
size_t num_entries)
|
|
{
|
|
bool used;
|
|
doublereal dmax0, dmax1, error, rel_norm;
|
|
printf("\t\t%s currentDamp = %g\n", title, damp);
|
|
printf("\t\t I ysolnOld %13s ysolnNewRaw | ysolnNewTrial "
|
|
"%10s ysolnNewTrialRaw | solnWeight wtDelSoln wtDelSolnTrial\n", stepNorm_1, stepNorm_2);
|
|
std::vector<size_t> imax(num_entries, npos);
|
|
printf("\t\t ");
|
|
print_line("-", 125);
|
|
for (size_t jnum = 0; jnum < num_entries; jnum++) {
|
|
dmax1 = -1.0;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
used = false;
|
|
for (size_t j = 0; j < jnum; j++) {
|
|
if (imax[j] == i) {
|
|
used = true;
|
|
}
|
|
}
|
|
if (!used) {
|
|
error = step_1[i] / m_ewt[i];
|
|
rel_norm = sqrt(error * error);
|
|
error = step_2[i] / m_ewt[i];
|
|
rel_norm += sqrt(error * error);
|
|
if (rel_norm > dmax1) {
|
|
imax[jnum] = i;
|
|
dmax1 = rel_norm;
|
|
}
|
|
}
|
|
}
|
|
if (imax[jnum] != npos) {
|
|
size_t i = imax[jnum];
|
|
error = step_1[i] / m_ewt[i];
|
|
dmax0 = sqrt(error * error);
|
|
error = step_2[i] / m_ewt[i];
|
|
dmax1 = sqrt(error * error);
|
|
printf("\t\t %4s %12.4e %12.4e %12.4e | %12.4e %12.4e %12.4e |%12.4e %12.4e %12.4e\n",
|
|
int2str(i).c_str(), y_n_curr[i], step_1[i], y_n_curr[i] + step_1[i], y_n_1[i],
|
|
step_2[i], y_n_1[i]+ step_2[i], m_ewt[i], dmax0, dmax1);
|
|
}
|
|
}
|
|
printf("\t\t ");
|
|
print_line("-", 125);
|
|
}
|
|
//====================================================================================================================
|
|
//! This routine subtracts two numbers for one another
|
|
/*!
|
|
* This routine subtracts 2 numbers. If the difference is less
|
|
* than 1.0E-14 times the magnitude of the smallest number, then diff returns an exact zero.
|
|
* It also returns an exact zero if the difference is less than
|
|
* 1.0E-300.
|
|
*
|
|
* returns: a - b
|
|
*
|
|
* This routine is used in numerical differencing schemes in order
|
|
* to avoid roundoff errors resulting in creating Jacobian terms.
|
|
* Note: This is a slow routine. However, jacobian errors may cause
|
|
* loss of convergence. Therefore, in practice this routine has proved cost-effective.
|
|
*
|
|
* @param a Value of a
|
|
* @param b value of b
|
|
*
|
|
* @return returns the difference between a and b
|
|
*/
|
|
static inline doublereal subtractRD(doublereal a, doublereal b)
|
|
{
|
|
doublereal diff = a - b;
|
|
doublereal d = std::min(fabs(a), fabs(b));
|
|
d *= 1.0E-14;
|
|
doublereal ad = fabs(diff);
|
|
if (ad < 1.0E-300) {
|
|
diff = 0.0;
|
|
}
|
|
if (ad < d) {
|
|
diff = 0.0;
|
|
}
|
|
return diff;
|
|
}
|
|
//====================================================================================================================
|
|
/*
|
|
*
|
|
* Function called by BEuler to evaluate the Jacobian matrix and the
|
|
* current residual at the current time step.
|
|
* @param N = The size of the equation system
|
|
* @param J = Jacobian matrix to be filled in
|
|
* @param f = Right hand side. This routine returns the current
|
|
* value of the rhs (output), so that it does
|
|
* not have to be computed again.
|
|
*
|
|
* @return Returns a flag to indicate that operation is successful.
|
|
* 1 Means a successful operation
|
|
* 0 Means an unsuccessful operation
|
|
*/
|
|
int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f,
|
|
doublereal time_curr, doublereal CJ,
|
|
doublereal* const y, doublereal* const ydot,
|
|
int num_newt_its)
|
|
{
|
|
double* col_j;
|
|
int info;
|
|
doublereal ysave, ydotsave, dy;
|
|
int retn = 1;
|
|
|
|
/*
|
|
* Clear the factor flag
|
|
*/
|
|
J.clearFactorFlag();
|
|
if (m_jacFormMethod == NSOLN_JAC_ANAL) {
|
|
/********************************************************************
|
|
* Call the function to get a jacobian.
|
|
*/
|
|
info = m_func->evalJacobian(time_curr, delta_t_n, CJ, y, ydot, J, f);
|
|
m_nJacEval++;
|
|
m_nfe++;
|
|
if (info != 1) {
|
|
return info;
|
|
}
|
|
} else {
|
|
if (J.matrixType_ == 0) {
|
|
/*******************************************************************
|
|
* Generic algorithm to calculate a numerical Jacobian
|
|
*/
|
|
/*
|
|
* Calculate the current value of the rhs given the
|
|
* current conditions.
|
|
*/
|
|
|
|
info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, f, JacBase_ResidEval);
|
|
m_nfe++;
|
|
if (info != 1) {
|
|
return info;
|
|
}
|
|
m_nJacEval++;
|
|
|
|
/*
|
|
* Malloc a vector and call the function object to return a set of
|
|
* deltaY's that are appropriate for calculating the numerical
|
|
* derivative.
|
|
*/
|
|
doublereal* dyVector = mdp::mdp_alloc_dbl_1((int) neq_, MDP_DBL_NOINIT);
|
|
retn = m_func->calcDeltaSolnVariables(time_curr, y, ydot, dyVector, DATA_PTR(m_ewt));
|
|
if (s_print_NumJac) {
|
|
if (m_print_flag >= 7) {
|
|
if (retn != 1) {
|
|
printf("\t\tbeuler_jac ERROR: calcDeltaSolnVariables() returned an error condition.\n");
|
|
printf("\t\t We will bail after calculating the Jacobian\n");
|
|
}
|
|
if (neq_ < 20) {
|
|
printf("\t\tUnk m_ewt y dyVector ResN\n");
|
|
for (size_t iii = 0; iii < neq_; iii++) {
|
|
printf("\t\t %4s %16.8e %16.8e %16.8e %16.8e \n",
|
|
int2str(iii).c_str(), m_ewt[iii], y[iii], dyVector[iii], f[iii]);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Loop over the variables, formulating a numerical derivative
|
|
* of the dense matrix.
|
|
* For the delta in the variable, we will use a variety of approaches
|
|
* The original approach was to use the error tolerance amount.
|
|
* This may not be the best approach, as it could be overly large in
|
|
* some instances and overly small in others.
|
|
* We will first protect from being overly small, by using the usual
|
|
* sqrt of machine precision approach, i.e., 1.0E-7,
|
|
* to bound the lower limit of the delta.
|
|
*/
|
|
for (size_t j = 0; j < neq_; j++) {
|
|
|
|
|
|
/*
|
|
* Get a pointer into the column of the matrix
|
|
*/
|
|
|
|
|
|
col_j = (doublereal*) J.ptrColumn(j);
|
|
ysave = y[j];
|
|
dy = dyVector[j];
|
|
//dy = fmaxx(1.0E-6 * m_ewt[j], fabs(ysave)*1.0E-7);
|
|
|
|
y[j] = ysave + dy;
|
|
dy = y[j] - ysave;
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
ydotsave = ydot[j];
|
|
ydot[j] += dy * CJ;
|
|
}
|
|
/*
|
|
* Call the function
|
|
*/
|
|
|
|
|
|
info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, DATA_PTR(m_wksp),
|
|
JacDelta_ResidEval, j, dy);
|
|
m_nfe++;
|
|
if (info != 1) {
|
|
mdp::mdp_safe_free((void**) &dyVector);
|
|
return info;
|
|
}
|
|
|
|
doublereal diff;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
diff = subtractRD(m_wksp[i], f[i]);
|
|
col_j[i] = diff / dy;
|
|
}
|
|
y[j] = ysave;
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
ydot[j] = ydotsave;
|
|
}
|
|
|
|
}
|
|
/*
|
|
* Release memory
|
|
*/
|
|
mdp::mdp_safe_free((void**) &dyVector);
|
|
} else if (J.matrixType_ == 1) {
|
|
size_t ku, kl;
|
|
size_t ivec[2];
|
|
size_t n = J.nRowsAndStruct(ivec);
|
|
kl = ivec[0];
|
|
ku = ivec[1];
|
|
if (n != neq_) {
|
|
printf("we have probs\n");
|
|
exit(-1);
|
|
}
|
|
|
|
// --------------------------------- BANDED MATRIX BRAIN DEAD ---------------------------------------------------
|
|
info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, f, JacBase_ResidEval);
|
|
m_nfe++;
|
|
if (info != 1) {
|
|
return info;
|
|
}
|
|
m_nJacEval++;
|
|
|
|
|
|
doublereal* dyVector = mdp::mdp_alloc_dbl_1((int) neq_, MDP_DBL_NOINIT);
|
|
retn = m_func->calcDeltaSolnVariables(time_curr, y, ydot, dyVector, DATA_PTR(m_ewt));
|
|
if (s_print_NumJac) {
|
|
if (m_print_flag >= 7) {
|
|
if (retn != 1) {
|
|
printf("\t\tbeuler_jac ERROR: calcDeltaSolnVariables() returned an error condition.\n");
|
|
printf("\t\t We will bail after calculating the Jacobian\n");
|
|
}
|
|
if (neq_ < 20) {
|
|
printf("\t\tUnk m_ewt y dyVector ResN\n");
|
|
for (size_t iii = 0; iii < neq_; iii++) {
|
|
printf("\t\t %4s %16.8e %16.8e %16.8e %16.8e \n",
|
|
int2str(iii).c_str(), m_ewt[iii], y[iii], dyVector[iii], f[iii]);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
for (size_t j = 0; j < neq_; j++) {
|
|
|
|
|
|
col_j = (doublereal*) J.ptrColumn(j);
|
|
ysave = y[j];
|
|
dy = dyVector[j];
|
|
|
|
|
|
y[j] = ysave + dy;
|
|
dy = y[j] - ysave;
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
ydotsave = ydot[j];
|
|
ydot[j] += dy * CJ;
|
|
}
|
|
|
|
info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, DATA_PTR(m_wksp), JacDelta_ResidEval, j, dy);
|
|
m_nfe++;
|
|
if (info != 1) {
|
|
mdp::mdp_safe_free((void**) &dyVector);
|
|
return info;
|
|
}
|
|
|
|
doublereal diff;
|
|
|
|
|
|
|
|
for (size_t i = j - ku; i <= j + kl; i++) {
|
|
if (i < neq_) {
|
|
diff = subtractRD(m_wksp[i], f[i]);
|
|
col_j[kl + ku + i - j] = diff / dy;
|
|
}
|
|
}
|
|
y[j] = ysave;
|
|
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
|
|
ydot[j] = ydotsave;
|
|
}
|
|
|
|
}
|
|
|
|
mdp::mdp_safe_free((void**) &dyVector);
|
|
double vSmall;
|
|
size_t ismall = J.checkRows(vSmall);
|
|
if (vSmall < 1.0E-100) {
|
|
printf("WE have a zero row, %s\n", int2str(ismall).c_str());
|
|
exit(-1);
|
|
}
|
|
ismall = J.checkColumns(vSmall);
|
|
if (vSmall < 1.0E-100) {
|
|
printf("WE have a zero column, %s\n", int2str(ismall).c_str());
|
|
exit(-1);
|
|
}
|
|
|
|
// ---------------------BANDED MATRIX BRAIN DEAD -----------------------
|
|
}
|
|
}
|
|
|
|
if (m_print_flag >= 7 && s_print_NumJac) {
|
|
if (neq_ < 30) {
|
|
printf("\t\tCurrent Matrix and Residual:\n");
|
|
printf("\t\t I,J | ");
|
|
for (size_t j = 0; j < neq_; j++) {
|
|
printf(" %5s ", int2str(j).c_str());
|
|
}
|
|
printf("| Residual \n");
|
|
printf("\t\t --");
|
|
for (size_t j = 0; j < neq_; j++) {
|
|
printf("------------");
|
|
}
|
|
printf("| -----------\n");
|
|
|
|
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
printf("\t\t %4s |", int2str(i).c_str());
|
|
for (size_t j = 0; j < neq_; j++) {
|
|
printf(" % 11.4E", J(i,j));
|
|
}
|
|
printf(" | % 11.4E\n", f[i]);
|
|
}
|
|
|
|
printf("\t\t --");
|
|
for (size_t j = 0; j < neq_; j++) {
|
|
printf("------------");
|
|
}
|
|
printf("--------------\n");
|
|
}
|
|
}
|
|
/*
|
|
* Make a copy of the data. Note, this jacobian copy occurs before any matrix scaling operations.
|
|
* It's the raw matrix producted by this routine.
|
|
*/
|
|
jacCopyPtr_->copyData(J);
|
|
|
|
return retn;
|
|
}
|
|
//====================================================================================================================
|
|
// Internal function to calculate the time derivative of the solution at the new step
|
|
/*
|
|
* Previously, the user must have supplied information about the previous time step for this routine to
|
|
* work as intended.
|
|
*
|
|
* @param order of the BDF method
|
|
* @param y_curr current value of the solution
|
|
* @param ydot_curr Calculated value of the solution derivative that is consistent with y_curr
|
|
*/
|
|
void NonlinearSolver::
|
|
calc_ydot(const int order, const doublereal* const y_curr, doublereal* const ydot_curr) const
|
|
{
|
|
if (!ydot_curr) {
|
|
return;
|
|
}
|
|
doublereal c1;
|
|
switch (order) {
|
|
case 0:
|
|
case 1: /* First order forward Euler/backward Euler */
|
|
c1 = 1.0 / delta_t_n;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
ydot_curr[i] = c1 * (y_curr[i] - m_y_nm1[i]);
|
|
}
|
|
return;
|
|
case 2: /* Second order Adams-Bashforth / Trapezoidal Rule */
|
|
c1 = 2.0 / delta_t_n;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
ydot_curr[i] = c1 * (y_curr[i] - m_y_nm1[i]) - m_ydot_nm1[i];
|
|
}
|
|
|
|
return;
|
|
default:
|
|
throw CanteraError("calc_ydot()", "Case not covered");
|
|
}
|
|
}
|
|
//====================================================================================================================
|
|
// Apply a filtering process to the new step
|
|
/*
|
|
* @param timeCurrent Current value of the time
|
|
* @param y_current current value of the solution
|
|
* @param ydot_current Current value of the solution derivative.
|
|
*
|
|
* @return Returns the norm of the value of the amount filtered
|
|
*/
|
|
doublereal NonlinearSolver::filterNewStep(const doublereal timeCurrent,
|
|
const doublereal* const ybase, doublereal* const step0)
|
|
{
|
|
doublereal tmp = m_func->filterNewStep(timeCurrent, ybase, step0);
|
|
return tmp;
|
|
}
|
|
//====================================================================================================================
|
|
// Apply a filtering process to the new solution
|
|
/*
|
|
* @param timeCurrent Current value of the time
|
|
* @param y_current current value of the solution
|
|
* @param ydot_current Current value of the solution derivative.
|
|
*
|
|
* @return Returns the norm of the value of the amount filtered
|
|
*/
|
|
doublereal NonlinearSolver::filterNewSolution(const doublereal timeCurrent,
|
|
doublereal* const y_current, doublereal* const ydot_current)
|
|
{
|
|
doublereal tmp = m_func->filterSolnPrediction(timeCurrent, y_current);
|
|
return tmp;
|
|
}
|
|
//====================================================================================================================
|
|
// Compute the Residual Weights
|
|
/*
|
|
* The residual weights are defined here to be equal to the inverse of the row scaling factors used to
|
|
* row scale the matrix, after column scaling is used. They are multiplied by rtol and an atol factor
|
|
* is added as well so that if the residual is less than 1, then the calculation is deemed to be converged.
|
|
*
|
|
* The basic idea is that a change in the solution vector on the order of the convergence tolerance
|
|
* multiplied by [RJC] which is of order one after row scaling should give you the relative weight
|
|
* of the row. Values of the residual for that row can then be normalized by the value of this weight.
|
|
* When the tolerance in delta x is achieved, the tolerance in the residual should also be achieved
|
|
* and should be checked.
|
|
*/
|
|
void
|
|
NonlinearSolver::computeResidWts()
|
|
{
|
|
ResidWtsReevaluated_ = true;
|
|
if (checkUserResidualTols_ == 1) {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
m_residWts[i] = userResidAtol_[i] + userResidRtol_ * m_rowWtScales[i] / neq_;
|
|
}
|
|
} else {
|
|
doublereal sum = 0.0;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
m_residWts[i] = m_rowWtScales[i] / neq_;
|
|
sum += m_residWts[i];
|
|
}
|
|
sum /= neq_;
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
m_residWts[i] = m_ScaleSolnNormToResNorm * (m_residWts[i] + atolBase_ * atolBase_ * sum);
|
|
}
|
|
if (checkUserResidualTols_ == 2) {
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
double uR = userResidAtol_[i] + userResidRtol_ * m_rowWtScales[i] / neq_;
|
|
m_residWts[i] = std::min(m_residWts[i], uR);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
//=====================================================================================================================
|
|
// return the residual weights
|
|
/*
|
|
* @param residWts Vector of length neq_
|
|
*/
|
|
void
|
|
NonlinearSolver::getResidWts(doublereal* const residWts) const
|
|
{
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
residWts[i] = (m_residWts)[i];
|
|
}
|
|
}
|
|
//=====================================================================================================================
|
|
// Check to see if the nonlinear problem has converged
|
|
/*
|
|
*
|
|
* @return integer is returned. If positive, then the problem has converged
|
|
* 1 Successful step was taken: Next step's norm is less than 1.0.
|
|
* The final residual norm is less than 1.0.
|
|
* 2 Successful step: Next step's norm is less than 0.8.
|
|
* This step's norm is less than 1.0.
|
|
* The residual norm can be anything.
|
|
* 3 Success: The final residual is less than 1.0E-2
|
|
* The predicted deltaSoln is below 1.0E-2.
|
|
* 0 Not converged yet
|
|
*/
|
|
int
|
|
NonlinearSolver::convergenceCheck(int dampCode, doublereal s1)
|
|
{
|
|
int retn = 0;
|
|
if (m_dampBound < 0.9999) {
|
|
return retn;
|
|
}
|
|
if (m_dampRes < 0.9999) {
|
|
return retn;
|
|
}
|
|
if (dampCode <= 0) {
|
|
return retn;
|
|
}
|
|
if (dampCode == 3) {
|
|
if (s1 < 1.0E-2) {
|
|
if (m_normResidTrial < 1.0E-6) {
|
|
return 3;
|
|
}
|
|
}
|
|
if (s1 < 0.8) {
|
|
if (m_normDeltaSoln_Newton < 1.0) {
|
|
return 2;
|
|
}
|
|
}
|
|
}
|
|
if (dampCode == 4) {
|
|
if (s1 < 1.0E-2) {
|
|
if (m_normResidTrial < 1.0E-6) {
|
|
return 3;
|
|
}
|
|
}
|
|
}
|
|
|
|
if (s1 < 0.8) {
|
|
if (m_normDeltaSoln_Newton < 1.0) {
|
|
return 2;
|
|
}
|
|
}
|
|
if (dampCode == 1 || dampCode == 2) {
|
|
if (s1 < 1.0) {
|
|
if (m_normResidTrial < 1.0) {
|
|
return 1;
|
|
}
|
|
}
|
|
}
|
|
return retn;
|
|
}
|
|
//=====================================================================================================================
|
|
// Set the absolute tolerances for the solution variables
|
|
/*
|
|
* Set the absolute tolerances used in the calculation
|
|
*
|
|
* @param atol Vector of length neq_ that contains the tolerances to be used for the solution variables
|
|
*/
|
|
void NonlinearSolver::setAtol(const doublereal* const atol)
|
|
{
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
atolk_[i]= atol[i];
|
|
}
|
|
}
|
|
//=====================================================================================================================
|
|
// Set the relative tolerances for the solution variables
|
|
/*
|
|
* Set the relative tolerances used in the calculation for the solution variables.
|
|
*
|
|
* @param rtol single double
|
|
*/
|
|
void NonlinearSolver::setRtol(const doublereal rtol)
|
|
{
|
|
rtol_ = rtol;
|
|
}
|
|
//=====================================================================================================================
|
|
// Set the relative and absolute tolerances for the Residual norm comparisons, if used
|
|
/*
|
|
*
|
|
* residWeightNorm[i] = residAtol[i] + residRtol * m_rowWtScales[i] / neq
|
|
*
|
|
* @param residNormHandling Parameter that sets the default handling of the residual norms
|
|
* 0 The residual weighting vector is calculated to make sure that the solution
|
|
* norms are roughly 1 when the residual norm is roughly 1.
|
|
* This is the default if this routine is not called.
|
|
* 1 Use the user residual norm specified by the parameters in this routine
|
|
* 2 Use the minimum value of the residual weights calculcated by method 1 and 2.
|
|
* This is the default if this routine is called and this parameter isn't specified.
|
|
*/
|
|
void NonlinearSolver::setResidualTols(double residRtol, double* residATol, int residNormHandling)
|
|
{
|
|
if (residNormHandling < 0 || residNormHandling > 2) {
|
|
throw CanteraError("NonlinearSolver::setResidualTols()",
|
|
"Unknown int for residNormHandling");
|
|
}
|
|
checkUserResidualTols_ = residNormHandling;
|
|
userResidRtol_ = residRtol;
|
|
if (residATol) {
|
|
userResidAtol_.resize(neq_);
|
|
for (size_t i = 0; i < neq_; i++) {
|
|
userResidAtol_[i] = residATol[i];
|
|
}
|
|
} else {
|
|
if (residNormHandling ==1 || residNormHandling == 2) {
|
|
throw CanteraError("NonlinearSolver::setResidualTols()",
|
|
"Must set residATol vector");
|
|
}
|
|
}
|
|
}
|
|
//=====================================================================================================================
|
|
void NonlinearSolver::setPrintLvl(int printLvl)
|
|
{
|
|
m_print_flag = printLvl;
|
|
}
|
|
//=====================================================================================================================
|
|
}
|
|
|