Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_ImplicitBDFStepper_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#ifndef Rythmos_IMPLICITBDF_STEPPER_DEF_H
30#define Rythmos_IMPLICITBDF_STEPPER_DEF_H
31
32#include "Rythmos_ImplicitBDFStepper_decl.hpp"
33#include "Rythmos_StepperHelpers.hpp"
34#include "Rythmos_ImplicitBDFStepperStepControl.hpp"
35
36namespace Rythmos {
37
38// ////////////////////////////
39// Defintions
40
41// Nonmember constructor
42template<class Scalar>
43RCP<ImplicitBDFStepper<Scalar> > implicitBDFStepper() {
44 RCP<ImplicitBDFStepper<Scalar> >
45 stepper = rcp(new ImplicitBDFStepper<Scalar>() );
46 return stepper;
47}
48
49template<class Scalar>
50RCP<ImplicitBDFStepper<Scalar> > implicitBDFStepper(
51 const RCP<Thyra::ModelEvaluator<Scalar> >& model,
52 const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
53 )
54{
55 RCP<ImplicitBDFStepper<Scalar> >
56 stepper = Teuchos::rcp(new ImplicitBDFStepper<Scalar>(model, solver));
57 return stepper;
58}
59
60template<class Scalar>
61RCP<ImplicitBDFStepper<Scalar> > implicitBDFStepper(
62 const RCP<Thyra::ModelEvaluator<Scalar> >& model,
63 const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
64 const RCP<Teuchos::ParameterList>& parameterList
65 )
66{
67 RCP<ImplicitBDFStepper<Scalar> >
68 stepper = Teuchos::rcp(new ImplicitBDFStepper<Scalar>(model,
69 solver,
70 parameterList));
71 return stepper;
72}
73
74// Constructors, intializers, Misc.
75
76
77template<class Scalar>
79{
80 this->defaultInitializeAll_();
81 haveInitialCondition_ = false;
82 isInitialized_=false;
83}
84
85
86template<class Scalar>
88 const RCP<Thyra::ModelEvaluator<Scalar> >& model
89 ,const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
90 ,const RCP<Teuchos::ParameterList>& parameterList
91 )
92{
93 this->defaultInitializeAll_();
94 this->setParameterList(parameterList);
95 // Now we instantiate the model and the solver
96 setModel(model);
97 setSolver(solver);
98 haveInitialCondition_ = false;
99 isInitialized_=false;
100}
101
102
103template<class Scalar>
105 const RCP<Thyra::ModelEvaluator<Scalar> >& model
106 ,const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
107 )
108{
109 this->defaultInitializeAll_();
110 // Now we instantiate the model and the solver
111 setModel(model);
112 setSolver(solver);
113 haveInitialCondition_ = false;
114 isInitialized_=false;
115}
116
117template<class Scalar>
118const Thyra::VectorBase<Scalar>&
120{
121 TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,
122 "Error, attempting to call getxHistory before initialization!\n");
123 TEUCHOS_TEST_FOR_EXCEPT( !(( 0 <= index ) && ( index <= maxOrder_ )) );
124 TEUCHOS_TEST_FOR_EXCEPT( !( index <= usedOrder_+1 ) );
125 return(*(xHistory_[index]));
126}
127
128template<class Scalar>
129void ImplicitBDFStepper<Scalar>::setStepControlStrategy(const RCP<StepControlStrategyBase<Scalar> >& stepControl)
130{
131 TEUCHOS_TEST_FOR_EXCEPTION(stepControl == Teuchos::null,std::logic_error,"Error, stepControl == Teuchos::null!\n");
132 stepControl_ = stepControl;
133}
134
135template<class Scalar>
137{
138 return(stepControl_);
139}
140
141template<class Scalar>
142RCP<const StepControlStrategyBase<Scalar> > ImplicitBDFStepper<Scalar>::getStepControlStrategy() const
143{
144 return(stepControl_);
145}
146
147
148// Overridden from SolverAcceptingStepperBase
149
150
151template<class Scalar>
152void ImplicitBDFStepper<Scalar>::setSolver(const RCP<Thyra::NonlinearSolverBase<Scalar> > &solver)
153{
154 TEUCHOS_TEST_FOR_EXCEPT(solver == Teuchos::null)
155 solver_ = solver;
156}
157
158
159template<class Scalar>
160RCP<Thyra::NonlinearSolverBase<Scalar> >
162{
163 return (solver_);
164}
165
166
167template<class Scalar>
168RCP<const Thyra::NonlinearSolverBase<Scalar> >
170{
171 return (solver_);
172}
173
174
175// Overridden from StepperBase
176
177
178template<class Scalar>
180{
181 return true;
182}
183
184template<class Scalar>
186{
187 return true;
188}
189
190
191template<class Scalar>
192RCP<StepperBase<Scalar> >
194{
195
196 // Just use the interface to clone the algorithm in an basically
197 // uninitialized state
198
199 RCP<ImplicitBDFStepper<Scalar> >
200 stepper = Teuchos::rcp(new ImplicitBDFStepper<Scalar>());
201
202 if (!is_null(model_))
203 stepper->setModel(model_); // Shallow copy is okay!
204
205 if (!is_null(solver_))
206 stepper->setSolver(solver_->cloneNonlinearSolver().assert_not_null());
207
208 if (!is_null(parameterList_))
209 stepper->setParameterList(Teuchos::parameterList(*parameterList_));
210
211 if (!is_null(stepControl_)) {
212 if (stepControl_->supportsCloning())
213 stepper->setStepControlStrategy(
214 stepControl_->cloneStepControlStrategyAlgorithm().assert_not_null());
215 }
216
217 // At this point, we should have a valid algorithm state. What might be
218 // missing is the initial condition if it was not given in *model_ but was
219 // set explicitly. However, the specification for this function does not
220 // guarantee that the full state will be copied in any case!
221
222 return stepper;
223
224}
225
226
227template<class Scalar>
229 const RCP<const Thyra::ModelEvaluator<Scalar> >& model
230 )
231{
232 // typedef Teuchos::ScalarTraits<Scalar> ST; // unused
233 TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
234 assertValidModel( *this, *model );
235 model_ = model;
236}
237
238template<class Scalar>
240 const RCP<Thyra::ModelEvaluator<Scalar> >& model
241 )
242{
243 this->setModel(model); // TODO 09/09/09 tscoffe: use ConstNonconstObjectContainer!
244}
245
246
247template<class Scalar>
248RCP<const Thyra::ModelEvaluator<Scalar> >
250{
251 return model_;
252}
253
254
255template<class Scalar>
256RCP<Thyra::ModelEvaluator<Scalar> >
258{
259 return Teuchos::null;
260}
261
262
263template<class Scalar>
265 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
266 )
267{
268 typedef Teuchos::ScalarTraits<Scalar> ST;
269 // typedef Thyra::ModelEvaluatorBase MEB; // unused
270 TEUCHOS_TEST_FOR_EXCEPT(is_null(initialCondition.get_x()));
271 TEUCHOS_TEST_FOR_EXCEPT(is_null(initialCondition.get_x_dot()));
272 basePoint_ = initialCondition;
273 xn0_ = initialCondition.get_x()->clone_v();
274 xpn0_ = initialCondition.get_x_dot()->clone_v();
275 time_ = initialCondition.get_t();
276 // Generate vectors for use in calculations
277 x_dot_base_ = createMember(xpn0_->space());
278 V_S(x_dot_base_.ptr(),ST::zero());
279 ee_ = createMember(xn0_->space());
280 V_S(ee_.ptr(),ST::zero());
281 // x history
282 xHistory_.clear();
283 xHistory_.push_back(xn0_->clone_v());
284 xHistory_.push_back(xpn0_->clone_v());
285
286 haveInitialCondition_ = true;
287 if (isInitialized_) {
288 initialize_();
289 }
290}
291
292
293template<class Scalar>
294Thyra::ModelEvaluatorBase::InArgs<Scalar>
296{
297 return basePoint_;
298}
299
300
301template<class Scalar>
302Scalar ImplicitBDFStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepType)
303{
304
305 RYTHMOS_FUNC_TIME_MONITOR("Rythmos::ImplicitBDFStepper::takeStep");
306
307 using Teuchos::as;
308 using Teuchos::incrVerbLevel;
309 typedef Teuchos::ScalarTraits<Scalar> ST;
310 // typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag TScalarMag; // unused
311 typedef Thyra::NonlinearSolverBase<Scalar> NSB;
312 typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
313
314 RCP<Teuchos::FancyOStream> out = this->getOStream();
315 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
316 Teuchos::OSTab ostab(out,1,"takeStep");
317 VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
318
319 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
320 *out
321 << "\nEntering " << this->Teuchos::Describable::description()
322 << "::takeStep("<<dt<<","<<toString(stepType)<<") ...\n";
323 }
324
325 if (!isInitialized_) {
326 initialize_();
327 }
328
329 stepControl_->setOStream(out);
330 stepControl_->setVerbLevel(verbLevel);
331 stepControl_->setRequestedStepSize(*this,dt,stepType);
332
333 AttemptedStepStatusFlag status;
334 while (1) {
335 // Set up problem coefficients (and handle first step)
336 Scalar hh_old = hh_;
337 int desiredOrder;
338 stepControl_->nextStepSize(*this,&hh_,&stepType,&desiredOrder);
339 TEUCHOS_TEST_FOR_EXCEPT(!((1 <= desiredOrder) &&
340 (desiredOrder <= maxOrder_)));
341 TEUCHOS_TEST_FOR_EXCEPT(!(desiredOrder <= usedOrder_+1));
342 currentOrder_ = desiredOrder;
343 if (numberOfSteps_ == 0) {
344 psi_[0] = hh_;
345 if (nef_ == 0) {
346 Vt_S(xHistory_[1].ptr(),hh_);
347 } else {
348 Vt_S(xHistory_[1].ptr(),hh_/hh_old);
349 }
350 }
351 this->updateCoeffs_();
352 // compute predictor
353 obtainPredictor_();
354 // solve nonlinear problem (as follows)
355
356 //
357 // Setup the nonlinear equations:
358 //
359 // f_bar( x_dot_coeff * x_bar + x_dot_base,
360 // x_coeff * x_bar + x_base, t_base ) = 0
361 // x_dot_coeff = -alpha_s/dt
362 // x_dot_base = x_prime_pred + (alpha_s/dt) * x_pred
363 // x_coeff = 1
364 // x_base = 0
365 // t_base = tn+dt
366 //
367 Scalar coeff_x_dot = Scalar(-ST::one())*alpha_s_/hh_;
368 V_StVpStV( x_dot_base_.ptr(), ST::one(), *xpn0_, alpha_s_/hh_, *xn0_ );
369 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
370 *out << "model_ = " << std::endl;
371 model_->describe(*out,verbLevel);
372 *out << "basePoint_ = " << std::endl;
373 basePoint_.describe(*out,verbLevel);
374 *out << "coeff_x_dot = " << coeff_x_dot << std::endl;
375 *out << "x_dot_base_ = " << std::endl;
376 x_dot_base_->describe(*out,verbLevel);
377 *out << "time_+hh_ = " << time_+hh_ << std::endl;
378 *out << "xn0_ = " << std::endl;
379 xn0_->describe(*out,verbLevel);
380 }
381 neModel_.initializeSingleResidualModel(
382 model_, basePoint_,
383 coeff_x_dot, x_dot_base_,
384 ST::one(), Teuchos::null,
385 time_+hh_,
386 xn0_
387 );
388 //
389 // Solve the implicit nonlinear system to a tolerance of ???
390 //
391 if (solver_->getModel().get() != &neModel_) {
392 solver_->setModel( Teuchos::rcpFromRef(neModel_) );
393 }
394 // Rythmos::TimeStepNonlinearSolver uses a built in solveCriteria,
395 // so you can't pass one in.
396 // I believe this is the correct solveCriteria for IDA though.
397 /*
398 Thyra::SolveMeasureType nonlinear_solve_measure_type(
399 Thyra::SOLVE_MEASURE_NORM_RESIDUAL,Thyra::SOLVE_MEASURE_ONE);
400 // This should be changed to match the condition in IDA
401 TScalarMag tolerance = relErrTol_/TScalarMag(10.0);
402 Thyra::SolveCriteria<Scalar> nonlinearSolveCriteria(
403 nonlinear_solve_measure_type, tolerance);
404 Thyra::SolveStatus<Scalar> nonlinearSolveStatus =
405 solver_->solve( &*xn0_, &nonlinearSolveCriteria, &*delta_ );
406 */
407 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
408 *out << "xn0 = " << std::endl;
409 xn0_->describe(*out,verbLevel);
410 *out << "ee = " << std::endl;
411 ee_->describe(*out,verbLevel);
412 }
413 nonlinearSolveStatus_ = solver_->solve( &*xn0_, NULL, &*ee_ );
414 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
415 *out << "xn0 = " << std::endl;
416 xn0_->describe(*out,verbLevel);
417 *out << "ee = " << std::endl;
418 ee_->describe(*out,verbLevel);
419 }
420 // In the above solve, on input *xn0_ is the initial guess that comes from
421 // the predictor. On output, *xn0_ is the solved for solution and *ee_ is
422 // the difference computed from the intial guess in *xn0_ to the final
423 // solved value of *xn0_. This is needed for basic numerical stability.
424 if (nonlinearSolveStatus_.solveStatus == Thyra::SOLVE_STATUS_CONVERGED) {
425 newtonConvergenceStatus_ = 0;
426 }
427 else {
428 newtonConvergenceStatus_ = -1;
429 }
430
431 // check error and evaluate LTE
432 stepControl_->setCorrection(*this,xn0_,ee_,newtonConvergenceStatus_);
433 bool stepPass = stepControl_->acceptStep(*this,&LETvalue_);
434
435 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
436 *out << "xn0_ = " << std::endl;
437 xn0_->describe(*out,verbLevel);
438 *out << "xpn0_ = " << std::endl;
439 xpn0_->describe(*out,verbLevel);
440 *out << "ee_ = " << std::endl;
441 ee_->describe(*out,verbLevel);
442 for (int i=0; i<std::max(2,maxOrder_); ++i) {
443 *out << "xHistory_[" << i << "] = " << std::endl;
444 xHistory_[i]->describe(*out,verbLevel);
445 }
446 }
447
448 // Check LTE here:
449 if (!stepPass) { // stepPass = false
450 stepLETStatus_ = STEP_LET_STATUS_FAILED;
451 status = stepControl_->rejectStep(*this);
452 nef_++;
453 if (status == CONTINUE_ANYWAY) {
454 break;
455 } else {
456 restoreHistory_();
457 }
458 } else { // stepPass = true
459 stepLETStatus_ = STEP_LET_STATUS_PASSED;
460 break;
461 }
462 }
463
464 // 08/22/07 the history array must be updated before
465 // stepControl_->completeStep.
466 completeStep_();
467 stepControl_->completeStep(*this);
468
469 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
470 *out
471 << "\nLeaving " << this->Teuchos::Describable::description()
472 << "::takeStep("<<dt<<","<<toString(stepType)<<") ...\n";
473 }
474
475 return(usedStep_);
476}
477
478
479template<class Scalar>
480const StepStatus<Scalar> ImplicitBDFStepper<Scalar>::getStepStatus() const
481{
482
483 // 2007/08/24: rabartl: We agreed that getStepStatus() would be free
484 // so I have commented out removed all code that is not free
485
486 typedef Teuchos::ScalarTraits<Scalar> ST;
487 StepStatus<Scalar> stepStatus;
488 if (!isInitialized_) {
489 stepStatus.message = "This stepper is uninitialized.";
490 stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
491 stepStatus.stepSize = Scalar(-ST::one());
492 stepStatus.order = -1;
493 stepStatus.time = Scalar(-ST::one());
494// return(stepStatus);
495 }
496 else if (numberOfSteps_ > 0) {
497 stepStatus.stepStatus = STEP_STATUS_CONVERGED;
498 } else {
499 stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
500 }
501
502 stepStatus.stepLETStatus = stepLETStatus_;
503 stepStatus.stepSize = usedStep_;
504 stepStatus.order = usedOrder_;
505 stepStatus.time = time_;
506 stepStatus.stepLETValue = LETvalue_;
507 if(xHistory_.size() > 0)
508 stepStatus.solution = xHistory_[0];
509 else
510 stepStatus.solution = Teuchos::null;
511 stepStatus.solutionDot = Teuchos::null; // This is *not* free!
512 stepStatus.residual = Teuchos::null; // This is *not* free!
513
514 return(stepStatus);
515
516}
517
518
519// Overridden from InterpolationBufferBase
520
521
522template<class Scalar>
523RCP<const Thyra::VectorSpaceBase<Scalar> >
525{
526 //TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,"Error, attempting to call get_x_space before initialization!\n");
527 return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
528}
529
530
531template<class Scalar>
533 const Array<Scalar>& /* time_vec */,
534 const Array<RCP<const Thyra::VectorBase<Scalar> > >& /* x_vec */,
535 const Array<RCP<const Thyra::VectorBase<Scalar> > >& /* xdot_vec */
536 )
537{
538 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
539 "Error, addPoints is not implemented for ImplicitBDFStepper.\n");
540}
541
542
543template<class Scalar>
545{
546 if ( !isInitialized_ && haveInitialCondition_ )
547 return timeRange<Scalar>(time_,time_);
548 if ( !isInitialized_ && !haveInitialCondition_ )
549 return invalidTimeRange<Scalar>();
550 return timeRange<Scalar>(time_-usedStep_,time_);
551}
552
553
554template<class Scalar>
556 const Array<Scalar>& time_vec
557 ,Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec
558 ,Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
559 ,Array<ScalarMag>* accuracy_vec) const
560{
561 using Teuchos::as;
562 using Teuchos::constOptInArg;
563 using Teuchos::null;
564 using Teuchos::ptr;
565 typedef Teuchos::ScalarTraits<Scalar> ST;
566 typedef typename ST::magnitudeType mScalarMag;
567
568 TEUCHOS_ASSERT(haveInitialCondition_);
569 // Only do this if we're being called pre-initialization to get the IC.
570 if ( (numberOfSteps_ == -1) &&
571 (time_vec.length() == 1) &&
572 (compareTimeValues<Scalar>(time_vec[0],time_)==0) ) {
573 defaultGetPoints<Scalar>(
574 time_, constOptInArg(*xn0_), constOptInArg(*xpn0_),
575 time_, constOptInArg(*xn0_), constOptInArg(*xpn0_),
576 time_vec, ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
577 Ptr<InterpolatorBase<Scalar> >(null)
578 );
579 return;
580 }
581 TEUCHOS_ASSERT(isInitialized_);
582 RYTHMOS_FUNC_TIME_MONITOR("Rythmos::ImplicitBDFStepper::getPoints");
583 if (x_vec)
584 x_vec->clear();
585 if (xdot_vec)
586 xdot_vec->clear();
587 for (Teuchos::Ordinal i=0 ; i<time_vec.size() ; ++i) {
588 RCP<Thyra::VectorBase<Scalar> >
589 x_temp = createMember(xn0_->space());
590 RCP<Thyra::VectorBase<Scalar> >
591 xdot_temp = createMember(xn0_->space());
592 mScalarMag accuracy = -ST::zero();
593 interpolateSolution_(
594 time_vec[i], &*x_temp, &*xdot_temp,
595 accuracy_vec ? &accuracy : 0
596 );
597 if (x_vec)
598 x_vec->push_back(x_temp);
599 if (xdot_vec)
600 xdot_vec->push_back(xdot_temp);
601 if (accuracy_vec)
602 accuracy_vec->push_back(accuracy);
603 }
604 if ( as<int>(this->getVerbLevel()) >= as<int>(Teuchos::VERB_HIGH) ) {
605 RCP<Teuchos::FancyOStream> out = this->getOStream();
606 Teuchos::OSTab ostab(out,1,"getPoints");
607 *out << "Passing out the interpolated values:" << std::endl;
608 for (Teuchos::Ordinal i=0; i<time_vec.size() ; ++i) {
609 *out << "time_[" << i << "] = " << time_vec[i] << std::endl;
610 if (x_vec) {
611 *out << "x_vec[" << i << "] = " << std::endl;
612 (*x_vec)[i]->describe(*out,this->getVerbLevel());
613 }
614 if (xdot_vec) {
615 *out << "xdot_vec[" << i << "] = ";
616 if ( (*xdot_vec)[i] == Teuchos::null) {
617 *out << "Teuchos::null" << std::endl;
618 }
619 else {
620 *out << std::endl << Teuchos::describe(*(*xdot_vec)[i],this->getVerbLevel());
621 }
622 }
623 if (accuracy_vec)
624 *out << "accuracy[" << i << "] = " << (*accuracy_vec)[i] << std::endl;
625 }
626 }
627}
628
629
630template<class Scalar>
631void ImplicitBDFStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
632{
633 TEUCHOS_ASSERT( time_vec != NULL );
634 time_vec->clear();
635 if (!haveInitialCondition_) {
636 return;
637 }
638 if (numberOfSteps_ > 0) {
639 time_vec->push_back(time_-usedStep_);
640 }
641 time_vec->push_back(time_);
642}
643
644
645template<class Scalar>
646void ImplicitBDFStepper<Scalar>::removeNodes(Array<Scalar>& /* time_vec */)
647{
648 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
649 "Error, removeNodes is not implemented for ImplicitBDFStepper.\n");
650}
651
652
653template<class Scalar>
655{
656 if (!isInitialized_) {
657 return(-1);
658 }
659 return(usedOrder_);
660}
661
662
663// Overridden from Teuchos::ParameterListAcceptor
664
665
666template<class Scalar>
668 RCP<Teuchos::ParameterList> const& paramList
669 )
670{
671 TEUCHOS_TEST_FOR_EXCEPT(paramList == Teuchos::null);
672 paramList->validateParameters(*this->getValidParameters(),0);
673 parameterList_ = paramList;
674 Teuchos::readVerboseObjectSublist(&*parameterList_,this);
675}
676
677
678template<class Scalar>
680{
681 return(parameterList_);
682}
683
684
685template<class Scalar>
686RCP<Teuchos::ParameterList>
688{
689 RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
690 parameterList_ = Teuchos::null;
691 return(temp_param_list);
692}
693
694
695template<class Scalar>
696RCP<const Teuchos::ParameterList>
698{
699 static RCP<Teuchos::ParameterList> validPL;
700 if (is_null(validPL)) {
701 RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
702 // This line is required to pass StepperValidator UnitTest!
703 pl->sublist(RythmosStepControlSettings_name);
704 Teuchos::setupVerboseObjectSublist(&*pl);
705 validPL = pl;
706 }
707
708 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
709 if (Teuchos::as<int>(verbLevel) == Teuchos::VERB_HIGH) {
710 RCP<Teuchos::FancyOStream> out = this->getOStream();
711 Teuchos::OSTab ostab(out,1,"getValidParameters");
712 *out << "Setting up valid parameterlist." << std::endl;
713 validPL->print(*out);
714 }
715
716 return (validPL);
717
718}
719
720
721// Overridden from Teuchos::Describable
722
723
724template<class Scalar>
726{
727 std::ostringstream out;
728 out << this->Teuchos::Describable::description();
729 const TimeRange<Scalar> timeRange = this->getTimeRange();
730 if (timeRange.isValid())
731 out << " (timeRange="<<timeRange<<")";
732 else
733 out << " (This stepper is not initialized yet)";
734 out << std::endl;
735 return out.str();
736}
737
738
739template<class Scalar>
741 Teuchos::FancyOStream &out,
742 const Teuchos::EVerbosityLevel verbLevel
743 ) const
744{
745
746 using Teuchos::as;
747
748 if (!isInitialized_) {
749 out << this->description();
750 return;
751 }
752
753 if ( (as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT) ) ||
754 (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
755 )
756 {
757 out << this->description() << std::endl;
758 out << "model_ = " << Teuchos::describe(*model_,verbLevel);
759 out << "solver_ = " << Teuchos::describe(*solver_,verbLevel);
760 out << "neModel_ = " << Teuchos::describe(neModel_,verbLevel);
761 }
762 if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
763 out << "time_ = " << time_ << std::endl;
764 out << "hh_ = " << hh_ << std::endl;
765 out << "currentOrder_ = " << currentOrder_ << std::endl;
766 }
767 if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
768 out << "xn0_ = " << Teuchos::describe(*xn0_,verbLevel);
769 out << "xpn0_ = " << Teuchos::describe(*xpn0_,verbLevel);
770 out << "x_dot_base_ = " << Teuchos::describe(*x_dot_base_,verbLevel);
771 for (int i=0 ; i < std::max(2,maxOrder_) ; ++i) {
772 out << "xHistory_[" << i << "] = "
773 << Teuchos::describe(*xHistory_[i],verbLevel);
774 }
775 out << "ee_ = " << Teuchos::describe(*ee_,verbLevel);
776 }
777}
778
779
780// private
781
782
783// 2007/08/24: rabartl: Belos: We really should initialize all of this data in
784// a member initialization list but since there are like three constructors
785// this would mean that we would have to duplicate code (which is error prone)
786// or use a macro (which is not easy to debug). We really should remove all
787// but the default constructor which then would set this data once in the
788// initialization list.
789
790template<class Scalar>
792{
793 typedef Teuchos::ScalarTraits<Scalar> ST;
794 const Scalar nan = ST::nan(), one = ST::one(), zero = ST::zero();
795 // Initialize some data members to their rightful default values
796 haveInitialCondition_ = false;
797 isInitialized_ = false;
798 currentOrder_ = 1;
799 usedOrder_ = 1;
800 usedStep_ = zero;
801 // Initialize the rest of the private data members to invalid values to
802 // avoid uninitialed memory
803 time_ = nan;
804 hh_ = nan;
805 maxOrder_ = -1;
806 LETvalue_ = -one;
807 stepLETStatus_ = STEP_LET_STATUS_UNKNOWN;
808 alpha_s_ = -one;
809 numberOfSteps_ = -1;
810 nef_ = -1;
811 nscsco_ = -1;
812 newtonConvergenceStatus_ = -1;
813}
814
815
816template<class Scalar>
817void ImplicitBDFStepper<Scalar>::obtainPredictor_()
818{
819
820 using Teuchos::as;
821 typedef Teuchos::ScalarTraits<Scalar> ST;
822
823 if (!isInitialized_) {
824 return;
825 }
826
827 RCP<Teuchos::FancyOStream> out = this->getOStream();
828 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
829 Teuchos::OSTab ostab(out,1,"obtainPredictor_");
830 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
831 *out << "currentOrder_ = " << currentOrder_ << std::endl;
832 }
833
834 // prepare history array for prediction
835 for (int i=nscsco_;i<=currentOrder_;++i) {
836 Vt_S(xHistory_[i].ptr(),beta_[i]);
837 }
838
839 // evaluate predictor
840 V_V(xn0_.ptr(),*xHistory_[0]);
841 V_S(xpn0_.ptr(),ST::zero());
842 for (int i=1;i<=currentOrder_;++i) {
843 Vp_V(xn0_.ptr(),*xHistory_[i]);
844 Vp_StV(xpn0_.ptr(),gamma_[i],*xHistory_[i]);
845 }
846 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
847 *out << "xn0_ = " << std::endl;
848 xn0_->describe(*out,verbLevel);
849 *out << "xpn0_ = " << std::endl;
850 xpn0_->describe(*out,verbLevel);
851 }
852}
853
854
855template<class Scalar>
856void ImplicitBDFStepper<Scalar>::interpolateSolution_(
857 const Scalar& timepoint,
858 Thyra::VectorBase<Scalar>* x_ptr,
859 Thyra::VectorBase<Scalar>* xdot_ptr,
860 ScalarMag* accuracy_ptr
861 ) const
862{
863
864 // typedef std::numeric_limits<Scalar> NL; // unused
865 typedef Teuchos::ScalarTraits<Scalar> ST;
866
867#ifdef HAVE_RYTHMOS_DEBUG
868 TEUCHOS_TEST_FOR_EXCEPTION(
869 !isInitialized_,std::logic_error,
870 "Error, attempting to call interpolateSolution before initialization!\n");
871 const TimeRange<Scalar> currTimeRange = this->getTimeRange();
872 TEUCHOS_TEST_FOR_EXCEPTION(
873 !currTimeRange.isInRange(timepoint), std::logic_error,
874 "Error, timepoint = " << timepoint << " is not in the time range "
875 << currTimeRange << "!" );
876#endif
877
878 const Scalar tn = time_;
879 const int kused = usedOrder_;
880
881 // order of interpolation
882 int kord = kused;
883 if ( (kused == 0) || (timepoint == tn) ) {
884 kord = 1;
885 }
886
887 // Initialize interploation
888 Thyra::V_V(ptr(x_ptr),*xHistory_[0]);
889 Thyra::V_S(ptr(xdot_ptr),ST::zero());
890
891 // Add history array contributions
892 const Scalar delt = timepoint - tn;
893 Scalar c = ST::one(); // coefficient for interpolation of x
894 Scalar d = ST::zero(); // coefficient for interpolation of xdot
895 Scalar gam = delt/psi_[0]; // coefficient for interpolation
896 for (int j=1 ; j <= kord ; ++j) {
897 d = d*gam + c/psi_[j-1];
898 c = c*gam;
899 gam = (delt + psi_[j-1])/psi_[j];
900 Thyra::Vp_StV(ptr(x_ptr),c,*xHistory_[j]);
901 Thyra::Vp_StV(ptr(xdot_ptr),d,*xHistory_[j]);
902 }
903
904 // Set approximate accuracy
905 if (accuracy_ptr) {
906 *accuracy_ptr = Teuchos::ScalarTraits<Scalar>::pow(usedStep_,kord);
907 }
908
909}
910
911
912template<class Scalar>
913void ImplicitBDFStepper<Scalar>::updateHistory_()
914{
915
916 using Teuchos::as;
917
918 // Save Newton correction for potential order increase on next step.
919 if (usedOrder_ < maxOrder_) {
920 assign( xHistory_[usedOrder_+1].ptr(), *ee_ );
921 }
922 // Update history arrays
923 Vp_V( xHistory_[usedOrder_].ptr(), *ee_ );
924 for (int j=usedOrder_-1;j>=0;j--) {
925 Vp_V( xHistory_[j].ptr(), *xHistory_[j+1] );
926 }
927 RCP<Teuchos::FancyOStream> out = this->getOStream();
928 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
929 Teuchos::OSTab ostab(out,1,"updateHistory_");
930 if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
931 for (int i=0;i<std::max(2,maxOrder_);++i) {
932 *out << "xHistory_[" << i << "] = " << std::endl;
933 xHistory_[i]->describe(*out,verbLevel);
934 }
935 }
936
937}
938
939
940template<class Scalar>
941void ImplicitBDFStepper<Scalar>::restoreHistory_()
942{
943
944 using Teuchos::as;
945 typedef Teuchos::ScalarTraits<Scalar> ST;
946
947 // undo preparation of history array for prediction
948 for (int i=nscsco_;i<=currentOrder_;++i) {
949 Vt_S( xHistory_[i].ptr(), ST::one()/beta_[i] );
950 }
951 for (int i=1;i<=currentOrder_;++i) {
952 psi_[i-1] = psi_[i] - hh_;
953 }
954 RCP<Teuchos::FancyOStream> out = this->getOStream();
955 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
956 Teuchos::OSTab ostab(out,1,"restoreHistory_");
957 if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
958 for (int i=0;i<maxOrder_;++i) {
959 *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
960 }
961 for (int i=0;i<maxOrder_;++i) {
962 *out << "xHistory_[" << i << "] = " << std::endl;
963 xHistory_[i]->describe(*out,verbLevel);
964 }
965 }
966
967}
968
969
970template<class Scalar>
971void ImplicitBDFStepper<Scalar>::updateCoeffs_()
972{
973
974 using Teuchos::as;
975 typedef Teuchos::ScalarTraits<Scalar> ST;
976
977 // If the number of steps taken with constant order and constant stepsize is
978 // more than the current order + 1 then we don't bother to update the
979 // coefficients because we've reached a constant step-size formula. When
980 // this is is not true, then we update the coefficients for the variable
981 // step-sizes.
982 if ((hh_ != usedStep_) || (currentOrder_ != usedOrder_)) {
983 nscsco_ = 0;
984 }
985 nscsco_ = std::min(nscsco_+1,usedOrder_+2);
986 if (currentOrder_+1 >= nscsco_) {
987 beta_[0] = ST::one();
988 alpha_[0] = ST::one();
989 Scalar temp1 = hh_;
990 gamma_[0] = ST::zero();
991 for (int i=1;i<=currentOrder_;++i) {
992 Scalar temp2 = psi_[i-1];
993 psi_[i-1] = temp1;
994 beta_[i] = beta_[i-1]*psi_[i-1]/temp2;
995 temp1 = temp2 + hh_;
996 alpha_[i] = hh_/temp1;
997 gamma_[i] = gamma_[i-1]+alpha_[i-1]/hh_;
998 }
999 psi_[currentOrder_] = temp1;
1000 }
1001 alpha_s_ = ST::zero();
1002 for (int i=0;i<currentOrder_;++i) {
1003 alpha_s_ = alpha_s_ - Scalar(ST::one()/(i+ST::one()));
1004 }
1005 RCP<Teuchos::FancyOStream> out = this->getOStream();
1006 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
1007 Teuchos::OSTab ostab(out,1,"updateCoeffs_");
1008 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1009 for (int i=0;i<=maxOrder_;++i) {
1010 *out << "alpha_[" << i << "] = " << alpha_[i] << std::endl;
1011 *out << "beta_[" << i << "] = " << beta_[i] << std::endl;
1012 *out << "gamma_[" << i << "] = " << gamma_[i] << std::endl;
1013 *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
1014 *out << "alpha_s_ = " << alpha_s_ << std::endl;
1015 }
1016 //std::cout << "alpha_s_ = " << alpha_s_ << std::endl;
1017 //for (int i=0;i<=maxOrder_;++i) {
1018 // std::cout << " alpha_[" << i << "] = " << alpha_[i] << std::endl;
1019 // std::cout << " beta_[" << i << "] = " << beta_[i] << std::endl;
1020 // std::cout << " gamma_[" << i << "] = " << gamma_[i] << std::endl;
1021 // std::cout << " psi_[" << i << "] = " << psi_[i] << std::endl;
1022 //}
1023 }
1024}
1025
1026
1027template<class Scalar>
1028void ImplicitBDFStepper<Scalar>::initialize_()
1029{
1030
1031 using Teuchos::as;
1032 typedef Teuchos::ScalarTraits<Scalar> ST;
1033 using Thyra::createMember;
1034
1035 RCP<Teuchos::FancyOStream> out = this->getOStream();
1036 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
1037 const bool doTrace = (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH));
1038 Teuchos::OSTab ostab(out,1,"initialize_");
1039
1040 if (doTrace) {
1041 *out
1042 << "\nEntering " << this->Teuchos::Describable::description()
1043 << "::initialize_()...\n";
1044 }
1045
1046 TEUCHOS_TEST_FOR_EXCEPT(model_ == Teuchos::null);
1047 TEUCHOS_TEST_FOR_EXCEPT(solver_ == Teuchos::null);
1048 TEUCHOS_ASSERT(haveInitialCondition_);
1049
1050 // Initialize Parameter List if none provided.
1051 if (parameterList_ == Teuchos::null) {
1052 RCP<Teuchos::ParameterList> emptyParameterList =
1053 Teuchos::rcp(new Teuchos::ParameterList);
1054 this->setParameterList(emptyParameterList);
1055 }
1056
1057 // Initialize StepControl
1058 if (stepControl_ == Teuchos::null) {
1059 RCP<ImplicitBDFStepperStepControl<Scalar> > implicitBDFStepperStepControl =
1060 Teuchos::rcp(new ImplicitBDFStepperStepControl<Scalar>());
1061 RCP<Teuchos::ParameterList> stepControlPL =
1062 Teuchos::sublist(parameterList_, RythmosStepControlSettings_name);
1063 implicitBDFStepperStepControl->setParameterList(stepControlPL);
1064 this->setStepControlStrategy(implicitBDFStepperStepControl);
1065 }
1066 stepControl_->initialize(*this);
1067
1068 maxOrder_ = stepControl_->getMaxOrder(); // maximum order
1069 TEUCHOS_TEST_FOR_EXCEPTION(
1070 !((1 <= maxOrder_) && (maxOrder_ <= 5)), std::logic_error,
1071 "Error, maxOrder returned from stepControl_->getMaxOrder() = "
1072 << maxOrder_ << " is outside range of [1,5]!\n");
1073
1074 Scalar zero = ST::zero();
1075
1076 currentOrder_ = 1; // Current order of integration
1077 usedOrder_ = 1; // order used in current step (used after currentOrder_
1078 // is updated)
1079 usedStep_ = zero;
1080 nscsco_ = 0;
1081 LETvalue_ = zero;
1082
1083 alpha_.clear(); // $\alpha_j(n)=h_n/\psi_j(n)$ coefficient used in
1084 // local error test
1085 // note: $h_n$ = current step size, n = current time step
1086 gamma_.clear(); // calculate time derivative of history array for predictor
1087 beta_.clear(); // coefficients used to evaluate predictor from history array
1088 psi_.clear(); // $\psi_j(n) = t_n-t_{n-j}$ intermediary variable used to
1089 // compute $\beta_j(n):$
1090 for (int i=0 ; i<=maxOrder_ ; ++i) {
1091 alpha_.push_back(zero);
1092 beta_.push_back(zero);
1093 gamma_.push_back(zero);
1094 psi_.push_back(zero);
1095 }
1096 alpha_s_=Scalar(-ST::one()); // $\alpha_s$ fixed-leading coefficient of
1097 // this BDF method
1098 hh_=zero;
1099 numberOfSteps_=0; // number of total time integration steps taken
1100 nef_ = 0;
1101
1102 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1103 *out << "alpha_s_ = " << alpha_s_ << std::endl;
1104 for (int i=0 ; i<=maxOrder_ ; ++i) {
1105 *out << "alpha_[" << i << "] = " << alpha_[i] << std::endl;
1106 *out << "beta_[" << i << "] = " << beta_[i] << std::endl;
1107 *out << "gamma_[" << i << "] = " << gamma_[i] << std::endl;
1108 *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
1109 }
1110 *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
1111 }
1112
1113 // setInitialCondition initialized xHistory with xn0, xpn0.
1114 // Now we add the rest of the vectors. Store maxOrder_+1 vectors
1115 for (int i=2 ; i<=maxOrder_ ; ++i) {
1116 xHistory_.push_back(createMember(xn0_->space()));
1117 V_S(xHistory_[i].ptr(),zero);
1118 }
1119
1120 isInitialized_ = true;
1121
1122 if (doTrace) {
1123 *out
1124 << "\nLeaving " << this->Teuchos::Describable::description()
1125 << "::initialize_()...\n";
1126 }
1127
1128}
1129
1130
1131template<class Scalar>
1132void ImplicitBDFStepper<Scalar>::completeStep_()
1133{
1134
1135 using Teuchos::as;
1136
1137#ifdef HAVE_RYTHMOS_DEBUG
1138 typedef Teuchos::ScalarTraits<Scalar> ST;
1139 TEUCHOS_TEST_FOR_EXCEPT(ST::isnaninf(hh_));
1140#endif
1141
1142 numberOfSteps_ ++;
1143 nef_ = 0;
1144 usedStep_ = hh_;
1145 usedOrder_ = currentOrder_;
1146 time_ += hh_;
1147 RCP<Teuchos::FancyOStream> out = this->getOStream();
1148 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
1149 Teuchos::OSTab ostab(out,1,"completeStep_");
1150
1151 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1152 *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
1153 *out << "time_ = " << time_ << std::endl;
1154 }
1155
1156 // 03/22/04 tscoffe: Note that updating the history has nothing to do with
1157 // the step-size and everything to do with the newton correction vectors.
1158 updateHistory_();
1159
1160}
1161
1162template<class Scalar>
1164 const StepperBase<Scalar> & stepper)
1165{
1166 if (!isInitialized_) {
1167 initialize_();
1168 }
1169 stepControl_->setStepControlData(stepper);
1170}
1171
1172template<class Scalar>
1173const Thyra::SolveStatus<Scalar>& ImplicitBDFStepper<Scalar>::getNonlinearSolveStatus() const
1174{
1175 return nonlinearSolveStatus_;
1176}
1177
1178//
1179// Explicit Instantiation macro
1180//
1181// Must be expanded from within the Rythmos namespace!
1182//
1183
1184#define RYTHMOS_IMPLICITBDF_STEPPER_INSTANT(SCALAR) \
1185 \
1186 template class ImplicitBDFStepper< SCALAR >; \
1187 \
1188 template RCP< ImplicitBDFStepper< SCALAR > > \
1189 implicitBDFStepper(); \
1190 \
1191 template RCP< ImplicitBDFStepper< SCALAR > > \
1192 implicitBDFStepper( \
1193 const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \
1194 const RCP<Thyra::NonlinearSolverBase< SCALAR > >& solver, \
1195 const RCP<Teuchos::ParameterList>& parameterList \
1196 ); \
1197
1198
1199} // namespace Rythmos
1200
1201
1202#endif //Rythmos_IMPLICITBDF_STEPPER_DEF_H
ImplicitBDFStepper()
Constructors, intializers, Misc.