Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_TimeStepNonlinearSolver_def.hpp
1//@HEADER
2// ***********************************************************************
3//
4// Rythmos Package
5// Copyright (2006) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
25//
26// ***********************************************************************
27//@HEADER
28
29
30#ifndef RYTHMOS_TIME_STEP_NONLINEAR_SOLVER_DEF_HPP
31#define RYTHMOS_TIME_STEP_NONLINEAR_SOLVER_DEF_HPP
32
33#include "Rythmos_TimeStepNonlinearSolver_decl.hpp"
34
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"
40
41namespace Rythmos {
42
43
44// ////////////////////////
45// Defintions
46
47
48// Static members
49
50
51template<class Scalar>
52const std::string
53TimeStepNonlinearSolver<Scalar>::DefaultTol_name_ = "Default Tol";
54
55template<class Scalar>
56const double
57TimeStepNonlinearSolver<Scalar>::DefaultTol_default_ = 1e-2;
58
59
60template<class Scalar>
61const std::string
62TimeStepNonlinearSolver<Scalar>::DefaultMaxIters_name_ = "Default Max Iters";
63
64template<class Scalar>
65const int
66TimeStepNonlinearSolver<Scalar>::DefaultMaxIters_default_ = 3;
67
68
69template<class Scalar>
70const std::string
71TimeStepNonlinearSolver<Scalar>::NonlinearSafetyFactor_name_
72= "Nonlinear Safety Factor";
73
74template<class Scalar>
75const double
76TimeStepNonlinearSolver<Scalar>::NonlinearSafetyFactor_default_ = 0.1;
77
78
79template<class Scalar>
80const std::string
81TimeStepNonlinearSolver<Scalar>::LinearSafetyFactor_name_ = "Linear Safety Factor";
82
83template<class Scalar>
84const double
85TimeStepNonlinearSolver<Scalar>::LinearSafetyFactor_default_ = 0.05;
86
87
88template<class Scalar>
89const std::string
90TimeStepNonlinearSolver<Scalar>::RMinFraction_name_ = "R Min Fraction";
91
92template<class Scalar>
93const double
94TimeStepNonlinearSolver<Scalar>::RMinFraction_default_ = 0.3;
95
96
97template<class Scalar>
98const std::string
99TimeStepNonlinearSolver<Scalar>::ThrownOnLinearSolveFailure_name_
100= "Thrown on Linear Solve Failure";
101
102template<class Scalar>
103const bool
104TimeStepNonlinearSolver<Scalar>::ThrownOnLinearSolveFailure_default_ = false;
105
106
107// Constructors/Intializers/Misc
108
109
110template <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_)
119{}
120
121
122// Overridden from ParameterListAcceptor
123
124
125template<class Scalar>
127 RCP<ParameterList> const& paramList
128 )
129{
130 using Teuchos::get;
131 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
132 paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
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
143 paramList_->validateParameters(*getValidParameters(),0);
144#endif // HAVE_RYTHMOS_DEBUG
145}
146
147
148template<class Scalar>
149RCP<ParameterList>
151{
152 return paramList_;
153}
154
155
156template<class Scalar>
157RCP<ParameterList>
159{
160 RCP<ParameterList> _paramList = paramList_;
161 paramList_ = Teuchos::null;
162 return _paramList;
163}
164
165
166template<class Scalar>
167RCP<const ParameterList>
169{
170 return paramList_;
171}
172
173
174template<class Scalar>
175RCP<const ParameterList>
177{
178 using Teuchos::setDoubleParameter; using Teuchos::setIntParameter;
179 static RCP<const ParameterList> validPL;
180 if (is_null(validPL)) {
181 RCP<ParameterList> pl = Teuchos::parameterList();
182 setDoubleParameter(
183 DefaultTol_name_, DefaultTol_default_,
184 "The default base tolerance for the nonlinear timestep solve.\n"
185 "This tolerance can be overridden ???",
186 &*pl );
187 setIntParameter(
188 DefaultMaxIters_name_, DefaultMaxIters_default_,
189 "The default maximum number of Newton iterations to perform.\n"
190 "This default can be overridden ???",
191 &*pl );
192 setDoubleParameter(
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.",
197 &*pl );
198 setDoubleParameter(
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.",
205 &*pl );
206 setDoubleParameter(
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.",
211 &*pl );
212 pl->set(
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."
216 );
217 Teuchos::setupVerboseObjectSublist(&*pl);
218 validPL = pl;
219 }
220 return validPL;
221}
222
223
224// Overridden from NonlinearSolverBase
225
226
227template <class Scalar>
229 const RCP<const Thyra::ModelEvaluator<Scalar> > &model
230 )
231{
232 TEUCHOS_TEST_FOR_EXCEPT(model.get()==NULL);
233 model_ = model;
234 J_ = Teuchos::null;
235 current_x_ = Teuchos::null;
236 J_is_current_ = false;
237}
238
239
240template <class Scalar>
241RCP<const Thyra::ModelEvaluator<Scalar> >
243{
244 return model_;
245}
246
247template <class Scalar>
248Thyra::SolveStatus<Scalar>
250 Thyra::VectorBase<Scalar> *x,
251 const Thyra::SolveCriteria<Scalar> *solveCriteria,
252 Thyra::VectorBase<Scalar> *delta
253 )
254{
255
256 RYTHMOS_FUNC_TIME_MONITOR("Rythmos:TimeStepNonlinearSolver::solve");
257
258 using std::endl;
259 using Teuchos::incrVerbLevel;
260 using Teuchos::describe;
261 using Teuchos::as;
262 using Teuchos::rcp;
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;
269
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!" );
277#else
278 (void)solveCriteria;
279#endif
280
281 const RCP<Teuchos::FancyOStream> out = this->getOStream();
282 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
283 const bool showNewtonDetails =
284 (!is_null(out) && (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)));
285 const bool dumpAll =
286 (!is_null(out) && (as<int>(verbLevel) == as<int>(Teuchos::VERB_EXTREME)));
287 TEUCHOS_OSTAB;
288 VOTSME stateModel_outputTempState(model_,out,incrVerbLevel(verbLevel,-1));
289
290 if (showNewtonDetails)
291 *out
292 << "\nEntering TimeStepNonlinearSolver::solve(...) ...\n"
293 << "\nmodel = " << Teuchos::describe(*model_,verbLevel);
294
295 if(dumpAll) {
296 *out << "\nInitial guess:\n";
297 *out << "\nx = " << *x;
298 }
299
300 // Initialize storage for algorithm
301 if(!J_.get()) J_ = model_->create_W();
302 TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::is_null(J_), std::logic_error,
303 "Error! model->create_W() returned a null pointer!\n"
304 );
305 RCP<Thyra::VectorBase<Scalar> > f = createMember(model_->get_f_space());
306 RCP<Thyra::VectorBase<Scalar> > dx = createMember(model_->get_x_space());
307 RCP<Thyra::VectorBase<Scalar> > dx_last = createMember(model_->get_x_space());
308 RCP<Thyra::VectorBase<Scalar> > x_curr = createMember(model_->get_x_space());
309 if (delta != NULL)
310 Thyra::V_S(ptr(delta),ST::zero()); // delta stores the cumulative update to x over the whole Newton solve.
311 Thyra::assign(x_curr.ptr(),*x);
312 J_is_current_ = false;
313 current_x_ = Teuchos::null;
314
315 // Initialize convergence criteria
316 ScalarMag R = SMT::one();
317 ScalarMag linearTolSafety = linearSafetyFactor_ * nonlinearSafetyFactor_;
318 int maxIters = defaultMaxIters_;
319 ScalarMag tol = defaultTol_;
320 // ToDo: Get above from solveCriteria!
321
322 // Do the undampened Newton iterations
323 bool converged = false;
324 bool sawFailedLinearSolve = false;
325 Thyra::SolveStatus<Scalar> failedLinearSolveStatus;
326 ScalarMag nrm_dx = SMT::nan();
327 ScalarMag nrm_dx_last = SMT::nan();
328 int iter = 1;
329 for( ; iter <= maxIters; ++iter ) {
330 if (showNewtonDetails)
331 *out << "\n*** newtonIter = " << iter << endl;
332 if (showNewtonDetails)
333 *out << "\nEvaluating the model f and W ...\n";
334 Thyra::eval_f_W( *model_, *x_curr, &*f, &*J_ );
335 if (showNewtonDetails)
336 *out << "\nSolving the system J*dx = -f ...\n";
337 Thyra::V_S(dx.ptr(),ST::zero()); // Initial guess is needed!
338 Thyra::SolveCriteria<Scalar>
339 linearSolveCriteria(
340 Thyra::SolveMeasureType(
341 Thyra::SOLVE_MEASURE_NORM_RESIDUAL, Thyra::SOLVE_MEASURE_NORM_RHS
342 ),
343 linearTolSafety*tol
344 );
345 VOTSLOWSB J_outputTempState(J_,out,incrVerbLevel(verbLevel,-1));
346 Thyra::SolveStatus<Scalar> linearSolveStatus
347 = J_->solve(Thyra::NOTRANS, *f, dx.ptr(), Teuchos::ptr(&linearSolveCriteria) );
348 if (showNewtonDetails)
349 *out << "\nLinear solve status:\n" << linearSolveStatus;
350 Thyra::Vt_S(dx.ptr(),Scalar(-ST::one()));
351 if (dumpAll)
352 *out << "\ndx = " << Teuchos::describe(*dx,verbLevel);
353 if (delta != NULL) {
354 Thyra::Vp_V(ptr(delta),*dx);
355 if (dumpAll)
356 *out << "\ndelta = " << Teuchos::describe(*delta,verbLevel);
357 }
358 // Check the linear solve
359 if(linearSolveStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) {
360 sawFailedLinearSolve = true;
361 failedLinearSolveStatus = linearSolveStatus;
362 if (throwOnLinearSolveFailure_) {
363 TEUCHOS_TEST_FOR_EXCEPTION(
364 throwOnLinearSolveFailure_, Thyra::CatastrophicSolveFailure,
365 "Error, the linear solver did not converge!"
366 );
367 }
368 if (showNewtonDetails)
369 *out << "\nWarning, linear solve did not converge! Continuing anyway :-)\n";
370 }
371 // Update the solution: x_curr = x_curr + dx
372 Vp_V( x_curr.ptr(), *dx );
373 if (dumpAll)
374 *out << "\nUpdated solution x = " << Teuchos::describe(*x_curr,verbLevel);
375 // Convergence test
376 nrm_dx = Thyra::norm(*dx);
377 if ( R*nrm_dx <= nonlinearSafetyFactor_*tol )
378 converged = true;
379 if (showNewtonDetails)
380 *out
381 << "\nConvergence test:\n"
382 << " R*||dx|| = " << R << "*" << nrm_dx
383 << " = " << (R*nrm_dx) << "\n"
384 << " <= nonlinearSafetyFactor*tol = " << nonlinearSafetyFactor_ << "*" << tol
385 << " = " << (nonlinearSafetyFactor_*tol)
386 << " : " << ( converged ? "converged!" : " unconverged" )
387 << endl;
388 if(converged)
389 break; // We have converged!!!
390 // Update convergence criteria for the next iteration ...
391 if(iter > 1) {
392 const Scalar
393 MinR = RMinFraction_*R,
394 nrm_dx_ratio = nrm_dx/nrm_dx_last;
395 R = std::max(MinR,nrm_dx_ratio);
396 if (showNewtonDetails)
397 *out
398 << "\nUpdated R\n"
399 << " = max(RMinFraction*R,||dx||/||dx_last||)\n"
400 << " = max("<<RMinFraction_<<"*"<<R<<","<<nrm_dx<<"/"<<nrm_dx_last<<")\n"
401 << " = max("<<MinR<<","<<nrm_dx_ratio<<")\n"
402 << " = " << R << endl;
403 }
404 // Save to old
405 std::swap(dx_last,dx);
406 nrm_dx_last = nrm_dx;
407 }
408
409 // Set the solution
410 Thyra::assign(ptr(x),*x_curr);
411
412 if (dumpAll)
413 *out << "\nFinal solution x = " << Teuchos::describe(*x,verbLevel);
414
415 // Check the status
416
417 Thyra::SolveStatus<Scalar> solveStatus;
418
419 std::ostringstream oss;
420 Teuchos::FancyOStream omsg(rcp(&oss,false));
421
422 omsg << "Solver: " << this->description() << endl;
423
424 if(converged) {
425 solveStatus.solveStatus = Thyra::SOLVE_STATUS_CONVERGED;
426 omsg << "CVODE status test converged!\n";
427 }
428 else {
429 solveStatus.solveStatus = Thyra::SOLVE_STATUS_UNCONVERGED;
430 omsg << "CVODE status test failed!\n";
431 }
432
433 if (sawFailedLinearSolve) {
434 omsg << "Warning! A failed linear solve was encountered with status:\n";
435 OSTab tab(omsg);
436 omsg << failedLinearSolveStatus;
437 }
438
439 omsg
440 << "R*||dx|| = " << R << "*" << nrm_dx
441 << " <= nonlinearSafetyFactor*tol = " << nonlinearSafetyFactor_ << "*" << tol << " : "
442 << ( converged ? "converged!" : " unconverged" ) << endl;
443
444 omsg
445 << "Iterations = " << iter;
446 // Above, we leave off the last newline since this is the convention for the
447 // SolveStatus::message string!
448
449 solveStatus.message = oss.str();
450
451 // Update the solution state for external clients
452 current_x_ = x->clone_v();
453 J_is_current_ = false;
454 // 2007/09/04: rabartl: Note, above the Jacobian J is always going to be out
455 // of date since this algorithm computes x_curr = x_curr + dx for at least
456 // one solve for dx = -inv(J)*f. Therefore, J is never at the updated
457 // x_curr, only the old x_curr!
458
459 if (showNewtonDetails)
460 *out << "\nLeaving TimeStepNonlinearSolver::solve(...) ...\n";
461
462 return solveStatus;
463
464}
465
466
467template <class Scalar>
469{
470 return true;
471}
472
473
474template <class Scalar>
475RCP<Thyra::NonlinearSolverBase<Scalar> >
477{
478 RCP<TimeStepNonlinearSolver<Scalar> >
479 nonlinearSolver = Teuchos::rcp(new TimeStepNonlinearSolver<Scalar>);
480 nonlinearSolver->model_ = model_; // Shallow copy is okay, model is stateless
481 nonlinearSolver->defaultTol_ = defaultTol_;
482 nonlinearSolver->defaultMaxIters_ = defaultMaxIters_;
483 nonlinearSolver->nonlinearSafetyFactor_ = nonlinearSafetyFactor_;
484 nonlinearSolver->linearSafetyFactor_ = linearSafetyFactor_;
485 nonlinearSolver->RMinFraction_ = RMinFraction_;
486 nonlinearSolver->throwOnLinearSolveFailure_ = throwOnLinearSolveFailure_;
487 // Note: The specification of this virtual function in the interface class
488 // allows us to just copy the algorithm, not the entire state so we are
489 // done!
490 return nonlinearSolver;
491}
492
493
494template <class Scalar>
495RCP<const Thyra::VectorBase<Scalar> >
497{
498 return current_x_;
499}
500
501
502template <class Scalar>
504{
505 return J_is_current_;
506}
507
508
509template <class Scalar>
510RCP<Thyra::LinearOpWithSolveBase<Scalar> >
512{
513 if (is_null(J_))
514 return Teuchos::null;
515 if (forceUpToDate) {
516#ifdef HAVE_RYTHMOS_DEBUG
517 TEUCHOS_TEST_FOR_EXCEPT(is_null(current_x_));
518#endif
519 Thyra::eval_f_W<Scalar>( *model_, *current_x_, 0, &*J_ );
520 J_is_current_ = true;
521 }
522 return J_;
523}
524
525
526template <class Scalar>
527RCP<const Thyra::LinearOpWithSolveBase<Scalar> >
529{
530 return J_;
531}
532
533
534template <class Scalar>
536{
537 J_is_current_ = W_is_current;
538}
539
540
541} // namespace Rythmos
542
543
544// Nonmember constructors
545
546
547template <class Scalar>
548Teuchos::RCP<Rythmos::TimeStepNonlinearSolver<Scalar> >
549Rythmos::timeStepNonlinearSolver()
550{
551 return Teuchos::rcp(new TimeStepNonlinearSolver<Scalar>);
552}
553
554
555template <class Scalar>
556Teuchos::RCP<Rythmos::TimeStepNonlinearSolver<Scalar> >
557Rythmos::timeStepNonlinearSolver(const RCP<ParameterList> &pl)
558{
559 const RCP<Rythmos::TimeStepNonlinearSolver<Scalar> >
560 solver = timeStepNonlinearSolver<Scalar>();
561 solver->setParameterList(pl);
562 return solver;
563}
564
565
566//
567// Explicit Instantiation macro
568//
569// Must be expanded from within the Rythmos namespace!
570//
571
572#define RYTHMOS_TIME_STEP_NONLINEAR_SOLVER_INSTANT(SCALAR) \
573 \
574 template class TimeStepNonlinearSolver< SCALAR >; \
575 \
576 template RCP<TimeStepNonlinearSolver< SCALAR > > timeStepNonlinearSolver(); \
577 \
578 template RCP<TimeStepNonlinearSolver<SCALAR > > \
579 timeStepNonlinearSolver(const RCP<ParameterList> &pl);
580
581
582
583#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 ParameterList > getParameterList() const
Thyra::SolveStatus< Scalar > solve(Thyra::VectorBase< Scalar > *x, const Thyra::SolveCriteria< Scalar > *solveCriteria, Thyra::VectorBase< Scalar > *delta=NULL)
RCP< const Thyra::VectorBase< Scalar > > get_current_x() const
RCP< Thyra::LinearOpWithSolveBase< Scalar > > get_nonconst_W(const bool forceUpToDate)
RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
void setParameterList(RCP< ParameterList > const &paramList)
RCP< const ParameterList > getValidParameters() const
void setModel(const RCP< const Thyra::ModelEvaluator< Scalar > > &model)
RCP< Thyra::NonlinearSolverBase< Scalar > > cloneNonlinearSolver() const
RCP< const Thyra::LinearOpWithSolveBase< Scalar > > get_W() const