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_Newton.H"
00043 #include "NOX_Common.H"
00044 #include "NOX_Abstract_Vector.H"
00045 #include "NOX_Abstract_Group.H"
00046 #include "NOX_Solver_Generic.H"
00047 #include "NOX_Solver_LineSearchBased.H"
00048 #include "NOX_Utils.H"
00049 #include "NOX_GlobalData.H"
00050
00051
00052 NOX::Direction::Newton::
00053 Newton(const Teuchos::RCP<NOX::GlobalData>& gd,
00054 Teuchos::ParameterList& p)
00055 {
00056 reset(gd, p);
00057 }
00058
00059 NOX::Direction::Newton::~Newton()
00060 {
00061 }
00062
00063 bool NOX::Direction::Newton::
00064 reset(const Teuchos::RCP<NOX::GlobalData>& gd,
00065 Teuchos::ParameterList& params)
00066 {
00067 globalDataPtr = gd;
00068 utils = gd->getUtils();
00069
00070 paramsPtr = ¶ms;
00071
00072 Teuchos::ParameterList& p = params.sublist("Newton");
00073
00074 doRescue = p.get("Rescue Bad Newton Solve", true);
00075 if (!p.sublist("Linear Solver").isParameter("Tolerance"))
00076 p.sublist("Linear Solver").get("Tolerance", 1.0e-10);
00077
00078
00079 if ( p.get("Forcing Term Method", "Constant") == "Constant" ) {
00080 useAdjustableForcingTerm = false;
00081 eta_k = p.sublist("Linear Solver").get("Tolerance", 1.0e-4);
00082 }
00083 else {
00084 useAdjustableForcingTerm = true;
00085 method = p.get("Forcing Term Method", "Type 1");
00086 eta_min = p.get("Forcing Term Minimum Tolerance", 1.0e-4);
00087 eta_max = p.get("Forcing Term Maximum Tolerance", 0.9);
00088 eta_initial = p.get("Forcing Term Initial Tolerance", 0.01);
00089 alpha = p.get("Forcing Term Alpha", 1.5);
00090 gamma = p.get("Forcing Term Gamma", 0.9);
00091 eta_k = eta_min;
00092 }
00093
00094 return true;
00095 }
00096
00097 bool NOX::Direction::Newton::compute(NOX::Abstract::Vector& dir,
00098 NOX::Abstract::Group& soln,
00099 const NOX::Solver::Generic& solver)
00100 {
00101 NOX::Abstract::Group::ReturnType status;
00102
00103
00104 status = soln.computeF();
00105 if (status != NOX::Abstract::Group::Ok)
00106 NOX::Direction::Newton::throwError("compute", "Unable to compute F");
00107
00108
00109 if (useAdjustableForcingTerm) {
00110 resetForcingTerm(soln, solver.getPreviousSolutionGroup(),
00111 solver.getNumIterations(), solver);
00112 }
00113 else {
00114 if (utils->isPrintType(Utils::Details)) {
00115 utils->out() << " CALCULATING FORCING TERM" << endl;
00116 utils->out() << " Method: Constant" << endl;
00117 utils->out() << " Forcing Term: " << eta_k << endl;
00118 }
00119 }
00120
00121
00122 status = soln.computeJacobian();
00123 if (status != NOX::Abstract::Group::Ok)
00124 NOX::Direction::Newton::throwError("compute", "Unable to compute Jacobian");
00125
00126
00127 status = soln.computeNewton(paramsPtr->sublist("Newton").sublist("Linear Solver"));
00128
00129
00130 if ((status != NOX::Abstract::Group::Ok) &&
00131 (doRescue == false)) {
00132 NOX::Direction::Newton::throwError("compute",
00133 "Unable to solve Newton system");
00134 }
00135 else if ((status != NOX::Abstract::Group::Ok) &&
00136 (doRescue == true)) {
00137 if (utils->isPrintType(NOX::Utils::Warning))
00138 utils->out() << "WARNING: NOX::Direction::Newton::compute() - Linear solve "
00139 << "failed to achieve convergence - using the step anyway "
00140 << "since \"Rescue Bad Newton Solve\" is true " << endl;
00141 }
00142
00143
00144 dir = soln.getNewton();
00145
00146 return true;
00147 }
00148
00149 bool NOX::Direction::Newton::compute(NOX::Abstract::Vector& dir,
00150 NOX::Abstract::Group& soln,
00151 const NOX::Solver::LineSearchBased& solver)
00152 {
00153 return NOX::Direction::Generic::compute( dir, soln, solver );
00154 }
00155
00156
00157 bool NOX::Direction::Newton::resetForcingTerm(const NOX::Abstract::Group& soln,
00158 const NOX::Abstract::Group& oldsoln,
00159 int niter,
00160 const NOX::Solver::Generic& solver)
00161 {
00162
00163
00164
00165 double eta_km1 = paramsPtr->sublist("Newton").sublist("Linear Solver")
00166 .get("Tolerance", 0.0);
00167
00168
00169
00170 const NOX::Solver::LineSearchBased* solverPtr = 0;
00171 solverPtr = dynamic_cast<const NOX::Solver::LineSearchBased*>(&solver);
00172 if (solverPtr != 0) {
00173 eta_km1 = 1.0 - solverPtr->getStepSize() * (1.0 - eta_km1);
00174 }
00175
00176 const string indent = " ";
00177
00178 if (utils->isPrintType(Utils::Details)) {
00179 utils->out() << indent << "CALCULATING FORCING TERM" << endl;
00180 utils->out() << indent << "Method: " << method << endl;
00181 }
00182
00183
00184 if (method == "Type 1") {
00185
00186 if (niter == 0) {
00187
00188 eta_k = eta_initial;
00189
00190 }
00191 else {
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201 if (Teuchos::is_null(predRhs)) {
00202 predRhs = oldsoln.getF().clone(ShapeCopy);
00203 }
00204 if (Teuchos::is_null(stepDir)) {
00205 stepDir = oldsoln.getF().clone(ShapeCopy);
00206 }
00207
00208
00209 stepDir->update(1.0, soln.getX(), -1.0, oldsoln.getX(), 0);
00210
00211
00212 if (!(oldsoln.isJacobian())) {
00213 if (utils->isPrintType(Utils::Details)) {
00214 utils->out() << "WARNING: NOX::Direction::Newton::resetForcingTerm() - "
00215 << "Jacobian is out of date! Recomputing Jacobian." << endl;
00216 }
00217 const_cast<NOX::Abstract::Group&>(oldsoln).computeJacobian();
00218 }
00219 oldsoln.applyJacobian(*stepDir, *predRhs);
00220
00221
00222 predRhs->update(1.0, oldsoln.getF(), 1.0);
00223
00224
00225 double normpredf = 0.0;
00226 double normf = 0.0;
00227 double normoldf = 0.0;
00228
00229 if (utils->isPrintType(Utils::Details)) {
00230 utils->out() << indent << "Forcing Term Norm: Using L-2 Norm."
00231 << endl;
00232 }
00233 normpredf = predRhs->norm();
00234 normf = soln.getNormF();
00235 normoldf = oldsoln.getNormF();
00236
00237
00238 eta_k = fabs(normf - normpredf) / normoldf;
00239
00240
00241 if (utils->isPrintType(Utils::Details)) {
00242 utils->out() << indent << "Residual Norm k-1 = " << normoldf << "\n";
00243 utils->out() << indent << "Residual Norm Linear Model k = " << normpredf << "\n";
00244 utils->out() << indent << "Residual Norm k = " << normf << "\n";
00245 utils->out() << indent << "Calculated eta_k (pre-bounds) = " << eta_k << endl;
00246 }
00247
00248
00249 const double tmp_alpha = (1.0 + sqrt(5.0)) / 2.0;
00250 const double eta_km1_alpha = pow(eta_km1, tmp_alpha);
00251 if (eta_km1_alpha > 0.1)
00252 eta_k = NOX_MAX(eta_k, eta_km1_alpha);
00253 eta_k = NOX_MAX(eta_k, eta_min);
00254 eta_k = NOX_MIN(eta_max, eta_k);
00255 }
00256 }
00257
00258 else if (method == "Type 2") {
00259
00260 if (niter == 0) {
00261
00262 eta_k = eta_initial;
00263
00264 }
00265 else {
00266
00267 double normf = 0.0;
00268 double normoldf = 0.0;
00269
00270 if (utils->isPrintType(Utils::Details)) {
00271 utils->out() << indent << "Forcing Term Norm: Using L-2 Norm."
00272 << endl;
00273 }
00274 normf = soln.getNormF();
00275 normoldf = oldsoln.getNormF();
00276
00277 const double residual_ratio = normf / normoldf;
00278
00279 eta_k = gamma * pow(residual_ratio, alpha);
00280
00281
00282 if (utils->isPrintType(Utils::Details)) {
00283 utils->out() << indent << "Residual Norm k-1 = " << normoldf << "\n";
00284 utils->out() << indent << "Residual Norm k = " << normf << "\n";
00285 utils->out() << indent << "Calculated eta_k (pre-bounds) = " << eta_k << endl;
00286 }
00287
00288
00289 const double eta_k_alpha = gamma * pow(eta_km1, alpha);
00290 if (eta_k_alpha > 0.1)
00291 eta_k = NOX_MAX(eta_k, eta_k_alpha);
00292 eta_k = NOX_MAX(eta_k, eta_min);
00293 eta_k = NOX_MIN(eta_max, eta_k);
00294 }
00295
00296 }
00297
00298 else {
00299
00300 if (utils->isPrintType(Utils::Warning))
00301 utils->out() << "NOX::Direction::Newton::resetForcingTerm "
00302 << "- invalid forcing term method (" << method << ")" << endl;
00303
00304 return false;
00305 }
00306
00307
00308 paramsPtr->sublist("Newton").sublist("Linear Solver")
00309 .set("Tolerance", eta_k);
00310
00311 if (utils->isPrintType(Utils::Details))
00312 utils->out() << indent << "Forcing Term: " << eta_k << endl;
00313
00314 return true;
00315 }
00316
00317 void NOX::Direction::Newton::throwError(const string& functionName, const string& errorMsg)
00318 {
00319 if (utils->isPrintType(Utils::Error))
00320 utils->err() << "NOX::Direction::Newton::" << functionName << " - " << errorMsg << endl;
00321 throw "NOX Error";
00322 }