00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "NOX_Direction_Utils_InexactNewton.H"
00043
00044 #include "NOX_Common.H"
00045 #include "NOX_Abstract_Vector.H"
00046 #include "NOX_Abstract_Group.H"
00047 #include "NOX_Solver_Generic.H"
00048 #include "NOX_Solver_LineSearchBased.H"
00049 #include "NOX_Utils.H"
00050 #include "NOX_GlobalData.H"
00051 #include "Teuchos_ParameterList.hpp"
00052
00053
00054
00055
00056 NOX::Direction::Utils::InexactNewton::
00057 InexactNewton(const Teuchos::RCP<NOX::GlobalData>& gd,
00058 Teuchos::ParameterList& directionSublist) :
00059 paramsPtr(0)
00060 {
00061 reset(gd, directionSublist);
00062 }
00063
00064
00065
00066
00067 NOX::Direction::Utils::InexactNewton::~InexactNewton()
00068 {
00069 }
00070
00071
00072
00073
00074 bool NOX::Direction::Utils::InexactNewton::
00075 reset(const Teuchos::RCP<NOX::GlobalData>& gd,
00076 Teuchos::ParameterList& directionSublist)
00077 {
00078 globalDataPtr = gd;
00079 printing = gd->getUtils();
00080 paramsPtr = &directionSublist;
00081
00082 directionMethod = paramsPtr->get("Method", "Newton");
00083
00084 Teuchos::ParameterList& p = paramsPtr->sublist(directionMethod);
00085
00086 setTolerance = p.get("Set Tolerance in Parameter List", true);
00087
00088 method = p.get("Forcing Term Method", "Constant");
00089
00090 if (method == "Constant") {
00091 forcingTermMethod = Constant;
00092 eta_k = p.sublist("Linear Solver").get("Tolerance", 1.0e-4);
00093 }
00094 else {
00095
00096 if (method == "Type 1") {
00097 forcingTermMethod = Type1;
00098 }
00099 else if (method == "Type 2") {
00100 forcingTermMethod = Type2;
00101 }
00102 else {
00103 throwError("reset", "\"Forcing Term Method\" is invalid!");
00104 }
00105
00106 eta_min = p.get("Forcing Term Minimum Tolerance", 1.0e-4);
00107 eta_max = p.get("Forcing Term Maximum Tolerance", 0.9);
00108 eta_initial = p.get("Forcing Term Initial Tolerance", 0.01);
00109 alpha = p.get("Forcing Term Alpha", 1.5);
00110 gamma = p.get("Forcing Term Gamma", 0.9);
00111 eta_k = eta_min;
00112
00113 }
00114
00115 return true;
00116 }
00117
00118
00119
00120
00121 double NOX::Direction::Utils::InexactNewton::
00122 computeForcingTerm(const NOX::Abstract::Group& soln,
00123 const NOX::Abstract::Group& oldsoln,
00124 int niter,
00125 const NOX::Solver::Generic& solver,
00126 double eta_last)
00127 {
00128 const string indent = " ";
00129
00130 if (forcingTermMethod == Constant) {
00131 if (printing->isPrintType(NOX::Utils::Details)) {
00132 printing->out() << indent << "CALCULATING FORCING TERM" << endl;
00133 printing->out() << indent << "Method: Constant" << endl;
00134 printing->out() << indent << "Forcing Term: " << eta_k << endl;
00135 }
00136 if (setTolerance)
00137 paramsPtr->sublist(directionMethod).sublist("Linear Solver").
00138 set("Tolerance", eta_k);
00139
00140 return eta_k;
00141 }
00142
00143
00144
00145
00146
00147 double eta_km1 = 0.0;
00148 if (eta_last < 0.0)
00149 eta_km1 = paramsPtr->sublist(directionMethod).
00150 sublist("Linear Solver").get("Tolerance", 0.0);
00151 else
00152 eta_km1 = eta_last;
00153
00154
00155
00156 const NOX::Solver::LineSearchBased* solverPtr = 0;
00157 solverPtr = dynamic_cast<const NOX::Solver::LineSearchBased*>(&solver);
00158 if (solverPtr != 0) {
00159 eta_km1 = 1.0 - solverPtr->getStepSize() * (1.0 - eta_km1);
00160 }
00161
00162 if (printing->isPrintType(NOX::Utils::Details)) {
00163 printing->out() << indent << "CALCULATING FORCING TERM" << endl;
00164 printing->out() << indent << "Method: " << method << endl;
00165 }
00166
00167
00168 if (forcingTermMethod == Type1) {
00169
00170 if (niter == 0) {
00171
00172 eta_k = eta_initial;
00173
00174 }
00175 else {
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185 if (Teuchos::is_null(predRhs)) {
00186 predRhs = oldsoln.getF().clone(ShapeCopy);
00187 }
00188 if (Teuchos::is_null(stepDir)) {
00189 stepDir = oldsoln.getF().clone(ShapeCopy);
00190 }
00191
00192
00193 stepDir->update(1.0, soln.getX(), -1.0, oldsoln.getX(), 0);
00194
00195
00196 if (!(oldsoln.isJacobian())) {
00197 if (printing->isPrintType(NOX::Utils::Details)) {
00198 printing->out() << "WARNING: NOX::InexactNewtonUtils::resetForcingTerm() - "
00199 << "Jacobian is out of date! Recomputing Jacobian." << endl;
00200 }
00201 const_cast<NOX::Abstract::Group&>(oldsoln).computeJacobian();
00202 }
00203 oldsoln.applyJacobian(*stepDir, *predRhs);
00204
00205
00206 predRhs->update(1.0, oldsoln.getF(), 1.0);
00207
00208
00209 double normpredf = predRhs->norm();
00210 double normf = soln.getNormF();
00211 double normoldf = oldsoln.getNormF();
00212
00213 if (printing->isPrintType(NOX::Utils::Details)) {
00214 printing->out() << indent << "Forcing Term Norm: Using L-2 Norm."
00215 << endl;
00216 }
00217
00218
00219 eta_k = fabs(normf - normpredf) / normoldf;
00220
00221
00222 if (printing->isPrintType(NOX::Utils::Details)) {
00223 printing->out() << indent << "Residual Norm k-1 = "
00224 << normoldf << "\n";
00225 printing->out() << indent << "Residual Norm Linear Model k = "
00226 << normpredf << "\n";
00227 printing->out() << indent << "Residual Norm k = " << normf << "\n";
00228 printing->out() << indent << "Calculated eta_k (pre-bounds) = " << eta_k << endl;
00229 }
00230
00231
00232 const double tmp = (1.0 + sqrt(5.0)) / 2.0;
00233 const double eta_km1_alpha = pow(eta_km1, tmp);
00234 if (eta_km1_alpha > 0.1)
00235 eta_k = NOX_MAX(eta_k, eta_km1_alpha);
00236 eta_k = NOX_MAX(eta_k, eta_min);
00237 eta_k = NOX_MIN(eta_max, eta_k);
00238 }
00239 }
00240
00241 else if (forcingTermMethod == Type2) {
00242
00243 if (niter == 0) {
00244
00245 eta_k = eta_initial;
00246
00247 }
00248 else {
00249
00250 double normf = soln.getNormF();
00251 double normoldf = oldsoln.getNormF();
00252
00253 if (printing->isPrintType(NOX::Utils::Details)) {
00254 printing->out() << indent << "Forcing Term Norm: Using L-2 Norm."
00255 << endl;
00256 }
00257
00258 const double residual_ratio = normf / normoldf;
00259
00260 eta_k = gamma * pow(residual_ratio, alpha);
00261
00262
00263 if (printing->isPrintType(NOX::Utils::Details)) {
00264 printing->out() << indent << "Residual Norm k-1 = "
00265 << normoldf << "\n";
00266 printing->out() << indent << "Residual Norm k = "
00267 << normf << "\n";
00268 printing->out() << indent << "Calculated eta_k (pre-bounds) = "
00269 << eta_k << endl;
00270 }
00271
00272
00273 const double eta_k_alpha = gamma * pow(eta_km1, alpha);
00274 if (eta_k_alpha > 0.1)
00275 eta_k = NOX_MAX(eta_k, eta_k_alpha);
00276 eta_k = NOX_MAX(eta_k, eta_min);
00277 eta_k = NOX_MIN(eta_max, eta_k);
00278 }
00279
00280 }
00281
00282
00283 if (setTolerance)
00284 paramsPtr->sublist(directionMethod).sublist("Linear Solver").
00285 set("Tolerance", eta_k);
00286
00287 if (printing->isPrintType(NOX::Utils::Details))
00288 printing->out() << indent << "Forcing Term: " << eta_k << endl;
00289
00290 return eta_k;
00291 }
00292
00293
00294
00295
00296 void NOX::Direction::Utils::InexactNewton::
00297 throwError(const string& functionName, const string& errorMsg)
00298 {
00299 if (printing->isPrintType(NOX::Utils::Error))
00300 printing->err() << "NOX::InexactNewtonUtils::" << functionName << " - "
00301 << errorMsg << endl;
00302 throw "NOX Error";
00303 }
00304
00305
00306
00307