30 #ifndef RYTHMOS_TIME_STEP_NONLINEAR_SOLVER_DEF_HPP 31 #define RYTHMOS_TIME_STEP_NONLINEAR_SOLVER_DEF_HPP 33 #include "Rythmos_TimeStepNonlinearSolver_decl.hpp" 35 #include "Thyra_TestingTools.hpp" 36 #include "Thyra_ModelEvaluatorHelpers.hpp" 37 #include "Teuchos_VerboseObjectParameterListHelpers.hpp" 38 #include "Teuchos_StandardParameterEntryValidators.hpp" 39 #include "Teuchos_as.hpp" 51 template<
class Scalar>
53 TimeStepNonlinearSolver<Scalar>::DefaultTol_name_ =
"Default Tol";
55 template<
class Scalar>
57 TimeStepNonlinearSolver<Scalar>::DefaultTol_default_ = 1e-2;
60 template<
class Scalar>
62 TimeStepNonlinearSolver<Scalar>::DefaultMaxIters_name_ =
"Default Max Iters";
64 template<
class Scalar>
66 TimeStepNonlinearSolver<Scalar>::DefaultMaxIters_default_ = 3;
69 template<
class Scalar>
71 TimeStepNonlinearSolver<Scalar>::NonlinearSafetyFactor_name_
72 =
"Nonlinear Safety Factor";
74 template<
class Scalar>
76 TimeStepNonlinearSolver<Scalar>::NonlinearSafetyFactor_default_ = 0.1;
79 template<
class Scalar>
81 TimeStepNonlinearSolver<Scalar>::LinearSafetyFactor_name_ =
"Linear Safety Factor";
83 template<
class Scalar>
85 TimeStepNonlinearSolver<Scalar>::LinearSafetyFactor_default_ = 0.05;
88 template<
class Scalar>
90 TimeStepNonlinearSolver<Scalar>::RMinFraction_name_ =
"R Min Fraction";
92 template<
class Scalar>
94 TimeStepNonlinearSolver<Scalar>::RMinFraction_default_ = 0.3;
97 template<
class Scalar>
99 TimeStepNonlinearSolver<Scalar>::ThrownOnLinearSolveFailure_name_
100 =
"Thrown on Linear Solve Failure";
102 template<
class Scalar>
104 TimeStepNonlinearSolver<Scalar>::ThrownOnLinearSolveFailure_default_ =
false;
110 template <
class Scalar>
112 :J_is_current_(false),
113 defaultTol_(DefaultTol_default_),
114 defaultMaxIters_(DefaultMaxIters_default_),
115 nonlinearSafetyFactor_(NonlinearSafetyFactor_default_),
116 linearSafetyFactor_(LinearSafetyFactor_default_),
117 RMinFraction_(RMinFraction_default_),
118 throwOnLinearSolveFailure_(ThrownOnLinearSolveFailure_default_)
125 template<
class Scalar>
127 RCP<ParameterList>
const& paramList
131 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
133 paramList_ = paramList;
134 defaultTol_ = get<double>(*paramList_,DefaultTol_name_);
135 defaultMaxIters_ = get<int>(*paramList_,DefaultMaxIters_name_);
136 nonlinearSafetyFactor_ = get<double>(*paramList_,NonlinearSafetyFactor_name_);
137 linearSafetyFactor_ = get<double>(*paramList_,LinearSafetyFactor_name_);
138 RMinFraction_ = get<double>(*paramList_,RMinFraction_name_);
139 throwOnLinearSolveFailure_ = get<bool>(
140 *paramList_,ThrownOnLinearSolveFailure_name_);
141 Teuchos::readVerboseObjectSublist(&*paramList_,
this);
142 #ifdef HAVE_RYTHMOS_DEBUG 144 #endif // HAVE_RYTHMOS_DEBUG 148 template<
class Scalar>
156 template<
class Scalar>
160 RCP<ParameterList> _paramList = paramList_;
161 paramList_ = Teuchos::null;
166 template<
class Scalar>
167 RCP<const ParameterList>
174 template<
class Scalar>
175 RCP<const ParameterList>
178 using Teuchos::setDoubleParameter;
using Teuchos::setIntParameter;
179 static RCP<const ParameterList> validPL;
180 if (is_null(validPL)) {
181 RCP<ParameterList> pl = Teuchos::parameterList();
183 DefaultTol_name_, DefaultTol_default_,
184 "The default base tolerance for the nonlinear timestep solve.\n" 185 "This tolerance can be overridden ???",
188 DefaultMaxIters_name_, DefaultMaxIters_default_,
189 "The default maximum number of Newton iterations to perform.\n" 190 "This default can be overridden ???",
193 NonlinearSafetyFactor_name_, NonlinearSafetyFactor_default_,
194 "The factor (< 1.0) to multiply tol to bound R*||dx|||.\n" 195 "The exact nonlinear convergence test is:\n" 196 " R*||dx|| <= \"" + NonlinearSafetyFactor_name_ +
"\" * tol.",
199 LinearSafetyFactor_name_, LinearSafetyFactor_default_,
200 "This factor multiplies the nonlinear safety factor which multiplies\n" 201 "tol when determining the linear solve tolerence.\n" 202 "The exact linear convergence tolerance is:\n" 203 " ||J*dx+f||/||f|| <= \"" + LinearSafetyFactor_name_ +
"\" * " 204 "\"" + NonlinearSafetyFactor_name_ +
"\" * tol.",
207 RMinFraction_name_, RMinFraction_default_,
208 "The faction below which the R factor is not allowed to drop\n" 209 "below each Newton iteration. The R factor is related to the\n" 210 "ratio of ||dx||/||dx_last|| between nonlinear iterations.",
213 ThrownOnLinearSolveFailure_name_, ThrownOnLinearSolveFailure_default_,
214 "If set to true (\"1\"), then an Thyra::CatastrophicSolveFailure\n" 215 "exception will be thrown when a linear solve fails to meet it's tolerance." 217 Teuchos::setupVerboseObjectSublist(&*pl);
227 template <
class Scalar>
229 const RCP<
const Thyra::ModelEvaluator<Scalar> > &model
232 TEUCHOS_TEST_FOR_EXCEPT(model.get()==NULL);
235 current_x_ = Teuchos::null;
236 J_is_current_ =
false;
240 template <
class Scalar>
241 RCP<const Thyra::ModelEvaluator<Scalar> >
247 template <
class Scalar>
248 Thyra::SolveStatus<Scalar>
250 Thyra::VectorBase<Scalar> *x,
251 const Thyra::SolveCriteria<Scalar> *solveCriteria,
252 Thyra::VectorBase<Scalar> *delta
256 RYTHMOS_FUNC_TIME_MONITOR(
"Rythmos:TimeStepNonlinearSolver::solve");
259 using Teuchos::incrVerbLevel;
260 using Teuchos::describe;
263 using Teuchos::OSTab;
264 using Teuchos::getFancyOStream;
265 typedef Thyra::ModelEvaluatorBase MEB;
266 typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
267 typedef Thyra::LinearOpWithSolveBase<Scalar> LOWSB;
268 typedef Teuchos::VerboseObjectTempState<LOWSB> VOTSLOWSB;
270 #ifdef HAVE_RYTHMOS_DEBUG 271 TEUCHOS_TEST_FOR_EXCEPT(0==x);
272 THYRA_ASSERT_VEC_SPACES(
273 "TimeStepNonlinearSolver<Scalar>::solve(...)",
274 *x->space(),*model_->get_x_space() );
275 TEUCHOS_TEST_FOR_EXCEPT(
276 0!=solveCriteria &&
"ToDo: Support passed in solve criteria!" );
279 const RCP<Teuchos::FancyOStream> out = this->getOStream();
280 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
281 const bool showNewtonDetails =
282 (!is_null(out) && (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)));
284 (!is_null(out) && (as<int>(verbLevel) == as<int>(Teuchos::VERB_EXTREME)));
286 VOTSME stateModel_outputTempState(model_,out,incrVerbLevel(verbLevel,-1));
288 if (showNewtonDetails)
290 <<
"\nEntering TimeStepNonlinearSolver::solve(...) ...\n" 291 <<
"\nmodel = " << Teuchos::describe(*model_,verbLevel);
294 *out <<
"\nInitial guess:\n";
295 *out <<
"\nx = " << *x;
299 if(!J_.get()) J_ = model_->create_W();
300 TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::is_null(J_), std::logic_error,
301 "Error! model->create_W() returned a null pointer!\n" 303 RCP<Thyra::VectorBase<Scalar> > f = createMember(model_->get_f_space());
304 RCP<Thyra::VectorBase<Scalar> > dx = createMember(model_->get_x_space());
305 RCP<Thyra::VectorBase<Scalar> > dx_last = createMember(model_->get_x_space());
306 RCP<Thyra::VectorBase<Scalar> > x_curr = createMember(model_->get_x_space());
308 Thyra::V_S(ptr(delta),ST::zero());
309 Thyra::assign(x_curr.ptr(),*x);
310 J_is_current_ =
false;
311 current_x_ = Teuchos::null;
315 ScalarMag linearTolSafety = linearSafetyFactor_ * nonlinearSafetyFactor_;
316 int maxIters = defaultMaxIters_;
321 bool converged =
false;
322 bool sawFailedLinearSolve =
false;
323 Thyra::SolveStatus<Scalar> failedLinearSolveStatus;
327 for( ; iter <= maxIters; ++iter ) {
328 if (showNewtonDetails)
329 *out <<
"\n*** newtonIter = " << iter << endl;
330 if (showNewtonDetails)
331 *out <<
"\nEvaluating the model f and W ...\n";
332 Thyra::eval_f_W( *model_, *x_curr, &*f, &*J_ );
333 if (showNewtonDetails)
334 *out <<
"\nSolving the system J*dx = -f ...\n";
335 Thyra::V_S(dx.ptr(),ST::zero());
336 Thyra::SolveCriteria<Scalar>
338 Thyra::SolveMeasureType(
339 Thyra::SOLVE_MEASURE_NORM_RESIDUAL, Thyra::SOLVE_MEASURE_NORM_RHS
343 VOTSLOWSB J_outputTempState(J_,out,incrVerbLevel(verbLevel,-1));
344 Thyra::SolveStatus<Scalar> linearSolveStatus
345 = J_->solve(Thyra::NOTRANS, *f, dx.ptr(), Teuchos::ptr(&linearSolveCriteria) );
346 if (showNewtonDetails)
347 *out <<
"\nLinear solve status:\n" << linearSolveStatus;
348 Thyra::Vt_S(dx.ptr(),Scalar(-ST::one()));
350 *out <<
"\ndx = " << Teuchos::describe(*dx,verbLevel);
352 Thyra::Vp_V(ptr(delta),*dx);
354 *out <<
"\ndelta = " << Teuchos::describe(*delta,verbLevel);
357 if(linearSolveStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) {
358 sawFailedLinearSolve =
true;
359 failedLinearSolveStatus = linearSolveStatus;
360 if (throwOnLinearSolveFailure_) {
361 TEUCHOS_TEST_FOR_EXCEPTION(
362 throwOnLinearSolveFailure_, Thyra::CatastrophicSolveFailure,
363 "Error, the linear solver did not converge!" 366 if (showNewtonDetails)
367 *out <<
"\nWarning, linear solve did not converge! Continuing anyway :-)\n";
370 Vp_V( x_curr.ptr(), *dx );
372 *out <<
"\nUpdated solution x = " << Teuchos::describe(*x_curr,verbLevel);
374 nrm_dx = Thyra::norm(*dx);
375 if ( R*nrm_dx <= nonlinearSafetyFactor_*tol )
377 if (showNewtonDetails)
379 <<
"\nConvergence test:\n" 380 <<
" R*||dx|| = " << R <<
"*" << nrm_dx
381 <<
" = " << (R*nrm_dx) <<
"\n" 382 <<
" <= nonlinearSafetyFactor*tol = " << nonlinearSafetyFactor_ <<
"*" << tol
383 <<
" = " << (nonlinearSafetyFactor_*tol)
384 <<
" : " << ( converged ?
"converged!" :
" unconverged" )
391 MinR = RMinFraction_*R,
392 nrm_dx_ratio = nrm_dx/nrm_dx_last;
393 R = std::max(MinR,nrm_dx_ratio);
394 if (showNewtonDetails)
397 <<
" = max(RMinFraction*R,||dx||/||dx_last||)\n" 398 <<
" = max("<<RMinFraction_<<
"*"<<R<<
","<<nrm_dx<<
"/"<<nrm_dx_last<<
")\n" 399 <<
" = max("<<MinR<<
","<<nrm_dx_ratio<<
")\n" 400 <<
" = " << R << endl;
403 std::swap(dx_last,dx);
404 nrm_dx_last = nrm_dx;
408 Thyra::assign(ptr(x),*x_curr);
411 *out <<
"\nFinal solution x = " << Teuchos::describe(*x,verbLevel);
415 Thyra::SolveStatus<Scalar> solveStatus;
417 std::ostringstream oss;
418 Teuchos::FancyOStream omsg(rcp(&oss,
false));
420 omsg <<
"Solver: " << this->description() << endl;
423 solveStatus.solveStatus = Thyra::SOLVE_STATUS_CONVERGED;
424 omsg <<
"CVODE status test converged!\n";
427 solveStatus.solveStatus = Thyra::SOLVE_STATUS_UNCONVERGED;
428 omsg <<
"CVODE status test failed!\n";
431 if (sawFailedLinearSolve) {
432 omsg <<
"Warning! A failed linear solve was encountered with status:\n";
434 omsg << failedLinearSolveStatus;
438 <<
"R*||dx|| = " << R <<
"*" << nrm_dx
439 <<
" <= nonlinearSafetyFactor*tol = " << nonlinearSafetyFactor_ <<
"*" << tol <<
" : " 440 << ( converged ?
"converged!" :
" unconverged" ) << endl;
443 <<
"Iterations = " << iter;
447 solveStatus.message = oss.str();
450 current_x_ = x->clone_v();
451 J_is_current_ =
false;
457 if (showNewtonDetails)
458 *out <<
"\nLeaving TimeStepNonlinearSolver::solve(...) ...\n";
465 template <
class Scalar>
472 template <
class Scalar>
473 RCP<Thyra::NonlinearSolverBase<Scalar> >
476 RCP<TimeStepNonlinearSolver<Scalar> >
478 nonlinearSolver->model_ = model_;
479 nonlinearSolver->defaultTol_ = defaultTol_;
480 nonlinearSolver->defaultMaxIters_ = defaultMaxIters_;
481 nonlinearSolver->nonlinearSafetyFactor_ = nonlinearSafetyFactor_;
482 nonlinearSolver->linearSafetyFactor_ = linearSafetyFactor_;
483 nonlinearSolver->RMinFraction_ = RMinFraction_;
484 nonlinearSolver->throwOnLinearSolveFailure_ = throwOnLinearSolveFailure_;
488 return nonlinearSolver;
492 template <
class Scalar>
493 RCP<const Thyra::VectorBase<Scalar> >
500 template <
class Scalar>
503 return J_is_current_;
507 template <
class Scalar>
508 RCP<Thyra::LinearOpWithSolveBase<Scalar> >
512 return Teuchos::null;
514 #ifdef HAVE_RYTHMOS_DEBUG 515 TEUCHOS_TEST_FOR_EXCEPT(is_null(current_x_));
517 Thyra::eval_f_W<Scalar>( *model_, *current_x_, 0, &*J_ );
518 J_is_current_ =
true;
524 template <
class Scalar>
525 RCP<const Thyra::LinearOpWithSolveBase<Scalar> >
532 template <
class Scalar>
535 J_is_current_ = W_is_current;
545 template <
class Scalar>
546 Teuchos::RCP<Rythmos::TimeStepNonlinearSolver<Scalar> >
547 Rythmos::timeStepNonlinearSolver()
553 template <
class Scalar>
554 Teuchos::RCP<Rythmos::TimeStepNonlinearSolver<Scalar> >
555 Rythmos::timeStepNonlinearSolver(
const RCP<ParameterList> &pl)
557 const RCP<Rythmos::TimeStepNonlinearSolver<Scalar> >
558 solver = timeStepNonlinearSolver<Scalar>();
559 solver->setParameterList(pl);
570 #define RYTHMOS_TIME_STEP_NONLINEAR_SOLVER_INSTANT(SCALAR) \ 572 template class TimeStepNonlinearSolver< SCALAR >; \ 574 template RCP<TimeStepNonlinearSolver< SCALAR > > timeStepNonlinearSolver(); \ 576 template RCP<TimeStepNonlinearSolver<SCALAR > > \ 577 timeStepNonlinearSolver(const RCP<ParameterList> &pl); 581 #endif // RYTHMOS_TIME_STEP_NONLINEAR_SOLVER_DEF_HPP Simple undampended Newton solver designed to solve time step equations in accurate times-tepping meth...
RCP< const Thyra::LinearOpWithSolveBase< Scalar > > get_W() const
RCP< ParameterList > getNonconstParameterList()
void setModel(const RCP< const Thyra::ModelEvaluator< Scalar > > &model)
RCP< const Thyra::VectorBase< Scalar > > get_current_x() const
void set_W_is_current(bool W_is_current)
bool is_W_current() const
RCP< ParameterList > unsetParameterList()
ST::magnitudeType ScalarMag
RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
RCP< Thyra::NonlinearSolverBase< Scalar > > cloneNonlinearSolver() const
Thyra::SolveStatus< Scalar > solve(Thyra::VectorBase< Scalar > *x, const Thyra::SolveCriteria< Scalar > *solveCriteria, Thyra::VectorBase< Scalar > *delta=NULL)
RCP< const ParameterList > getParameterList() const
RCP< Thyra::LinearOpWithSolveBase< Scalar > > get_nonconst_W(const bool forceUpToDate)
RCP< const ParameterList > getValidParameters() const
void setParameterList(RCP< ParameterList > const ¶mList)
bool supportsCloning() const
TimeStepNonlinearSolver()
Sets parameter defaults .