Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_BackwardEulerStepper_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_BACKWARD_EULER_STEPPER_DEF_H
30#define Rythmos_BACKWARD_EULER_STEPPER_DEF_H
31
32#include "Rythmos_BackwardEulerStepper_decl.hpp"
33#include "Rythmos_DataStore.hpp"
34#include "Rythmos_LinearInterpolator.hpp"
35#include "Rythmos_InterpolatorBaseHelpers.hpp"
36#include "Rythmos_StepperHelpers.hpp"
37#include "Rythmos_FixedStepControlStrategy.hpp"
38#include "Rythmos_SimpleStepControlStrategy.hpp"
39#include "Rythmos_FirstOrderErrorStepControlStrategy.hpp"
40
41#include "Thyra_ModelEvaluatorHelpers.hpp"
42#include "Thyra_AssertOp.hpp"
43#include "Thyra_TestingTools.hpp"
44
45#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
46#include "Teuchos_as.hpp"
47
48namespace Rythmos {
49
50
51template<class Scalar>
52RCP<BackwardEulerStepper<Scalar> >
53backwardEulerStepper(
54 const RCP<Thyra::ModelEvaluator<Scalar> >& model,
55 const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
56 )
57{
58 return Teuchos::rcp(new BackwardEulerStepper<Scalar>(model, solver));
59}
60
61template<class Scalar>
62RCP<BackwardEulerStepper<Scalar> >
63backwardEulerStepper()
64{
65 return Teuchos::rcp(new BackwardEulerStepper<Scalar>());
66}
67
68
69// ////////////////////////////
70// Defintions
71
72
73// Constructors, intializers, Misc.
74
75
76template<class Scalar>
78{
79 typedef Teuchos::ScalarTraits<Scalar> ST;
80 this->defaultInitializeAll_();
81 numSteps_ = 0;
82 dt_ = ST::zero();
83}
84
85
86template<class Scalar>
88 const RCP<Thyra::ModelEvaluator<Scalar> >& model,
89 const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
90 )
91{
92 typedef Teuchos::ScalarTraits<Scalar> ST;
93 this->defaultInitializeAll_();
94 dt_ = ST::zero();
95 numSteps_ = 0;
96 setModel(model);
97 setSolver(solver);
98}
99
100template<class Scalar>
102{
103 typedef Teuchos::ScalarTraits<Scalar> ST;
104 isInitialized_ = false;
105 haveInitialCondition_ = false;
106 model_ = Teuchos::null;
107 solver_ = Teuchos::null;
108 x_old_ = Teuchos::null;
109 scaled_x_old_ = Teuchos::null;
110 x_dot_old_ = Teuchos::null;
111 // basePoint_;
112 x_ = Teuchos::null;
113 x_dot_ = Teuchos::null;
114 dx_ = Teuchos::null;
115 t_ = ST::nan();
116 t_old_ = ST::nan();
117 dt_ = ST::nan();
118 numSteps_ = -1;
119 neModel_ = Teuchos::null;
120 parameterList_ = Teuchos::null;
121 interpolator_ = Teuchos::null;
122 stepControl_ = Teuchos::null;
123 newtonConvergenceStatus_ = -1;
124}
125
126// Overridden from InterpolatorAcceptingObjectBase
127
128template<class Scalar>
130 const RCP<InterpolatorBase<Scalar> >& interpolator
131 )
132{
133#ifdef HAVE_RYTHMOS_DEBUG
134 TEUCHOS_TEST_FOR_EXCEPT(is_null(interpolator));
135#endif
136 interpolator_ = interpolator;
137 isInitialized_ = false;
138}
139
140template<class Scalar>
141RCP<InterpolatorBase<Scalar> >
143{
144 return interpolator_;
145}
146
147template<class Scalar>
148RCP<const InterpolatorBase<Scalar> >
150{
151 return interpolator_;
152}
153
154template<class Scalar>
155RCP<InterpolatorBase<Scalar> >
157{
158 RCP<InterpolatorBase<Scalar> > temp_interpolator = interpolator_;
159 interpolator_ = Teuchos::null;
160 //isInitialized_ = false;
161 return temp_interpolator; // Fix Rythmos_BackwardEulerStepper_def.hpp:162:1: error: control reaches end of non-void function [-Werror=return-type]
162}
163
164
165// Overridden from SolverAcceptingStepperBase
166
167
168template<class Scalar>
170 const RCP<Thyra::NonlinearSolverBase<Scalar> > &solver
171 )
172{
173 using Teuchos::as;
174
175 TEUCHOS_TEST_FOR_EXCEPTION(solver == Teuchos::null, std::logic_error,
176 "Error! Thyra::NonlinearSolverBase RCP passed in through "
177 "BackwardEulerStepper::setSolver is null!");
178
179 RCP<Teuchos::FancyOStream> out = this->getOStream();
180 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
181 Teuchos::OSTab ostab(out,1,"BES::setSolver");
182 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
183 *out << "solver = " << solver->description() << std::endl;
184 }
185
186 solver_ = solver;
187
188 isInitialized_ = false;
189
190}
191
192
193template<class Scalar>
194RCP<Thyra::NonlinearSolverBase<Scalar> >
199
200
201template<class Scalar>
202RCP<const Thyra::NonlinearSolverBase<Scalar> >
204{
205 return solver_;
206}
207
208
209// Overridden from StepperBase
210
211
212template<class Scalar>
214{
215 return true;
216}
217
218
219template<class Scalar>
220RCP<StepperBase<Scalar> >
222{
223 RCP<BackwardEulerStepper<Scalar> >
224 stepper = Teuchos::rcp(new BackwardEulerStepper<Scalar>);
225 stepper->isInitialized_ = isInitialized_;
226 stepper->model_ = model_; // Model is stateless so shallow copy is okay!
227 if (!is_null(solver_))
228 stepper->solver_ = solver_->cloneNonlinearSolver().assert_not_null();
229 if (!is_null(x_))
230 stepper->x_ = x_->clone_v().assert_not_null();
231 if (!is_null(scaled_x_old_))
232 stepper->scaled_x_old_ = scaled_x_old_->clone_v().assert_not_null();
233 stepper->t_ = t_;
234 stepper->t_old_ = t_old_;
235 stepper->dt_ = dt_;
236 stepper->numSteps_ = numSteps_;
237 if (!is_null(neModel_))
238 stepper->neModel_
239 = Teuchos::rcp(new Rythmos::SingleResidualModelEvaluator<Scalar>);
240 if (!is_null(parameterList_))
241 stepper->parameterList_ = parameterList(*parameterList_);
242 if (!is_null(interpolator_))
243 stepper->interpolator_
244 = interpolator_->cloneInterpolator().assert_not_null(); // ToDo: Implement cloneInterpolator()
245 if (!is_null(stepControl_)) {
246 if (stepControl_->supportsCloning())
247 stepper->setStepControlStrategy(
248 stepControl_->cloneStepControlStrategyAlgorithm().assert_not_null());
249 }
250 return stepper;
251}
252
253template<class Scalar>
255{
256 return true;
257}
258
259template<class Scalar>
261 const RCP<const Thyra::ModelEvaluator<Scalar> >& model)
262{
263 using Teuchos::as;
264
265 TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
266 assertValidModel( *this, *model );
267
268 RCP<Teuchos::FancyOStream> out = this->getOStream();
269 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
270 Teuchos::OSTab ostab(out,1,"BES::setModel");
271 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
272 *out << "model = " << model->description() << std::endl;
273 }
274 model_ = model;
275}
276
277template<class Scalar>
279 const RCP<Thyra::ModelEvaluator<Scalar> >& model
280 )
281{
282 this->setModel(model); // TODO: 09/09/09 tscoffe: Use ConstNonconstObjectContainer here!
283}
284
285template<class Scalar>
286RCP<const Thyra::ModelEvaluator<Scalar> >
288{
289 return model_;
290}
291
292
293template<class Scalar>
294RCP<Thyra::ModelEvaluator<Scalar> >
296{
297 return Teuchos::null;
298}
299
300
301template<class Scalar>
303 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition)
304{
305
306 typedef Teuchos::ScalarTraits<Scalar> ST;
307 // typedef Thyra::ModelEvaluatorBase MEB; // unused
308
309 basePoint_ = initialCondition;
310
311 // x
312 RCP<const Thyra::VectorBase<Scalar> > x_init = initialCondition.get_x();
313 TEUCHOS_TEST_FOR_EXCEPTION(
314 is_null(x_init), std::logic_error,
315 "Error, if the client passes in an intial condition to "
316 "setInitialCondition(...), then x can not be null!" );
317 x_ = x_init->clone_v();
318
319 // x_dot
320 RCP<const Thyra::VectorBase<Scalar> >
321 x_dot_init = initialCondition.get_x_dot();
322 if (!is_null(x_dot_init)) {
323 x_dot_ = x_dot_init->clone_v();
324 }
325 else {
326 x_dot_ = createMember(x_->space());
327 V_S(x_dot_.ptr(),ST::zero());
328 }
329
330 // dx
331 if (is_null(dx_)) dx_ = createMember(x_->space());
332 V_S(dx_.ptr(), ST::zero());
333
334 // t
335 t_ = initialCondition.get_t();
336 t_old_ = t_;
337
338 // x_old
339 if (is_null(x_old_)) x_old_ = createMember(x_->space());
340 x_old_ = x_init->clone_v();
341 scaled_x_old_ = x_->clone_v();
342
343 // x_dot_old
344 if (is_null(x_dot_old_)) x_dot_old_ = createMember(x_->space());
345 x_dot_old_ = x_dot_->clone_v();
346
347 haveInitialCondition_ = true;
348}
349
350
351template<class Scalar>
352Thyra::ModelEvaluatorBase::InArgs<Scalar>
354{
355 return basePoint_;
356}
357
358template<class Scalar>
360 const RCP<StepControlStrategyBase<Scalar> >& stepControl)
361{
362 TEUCHOS_TEST_FOR_EXCEPTION(stepControl == Teuchos::null,std::logic_error,
363 "Error, stepControl == Teuchos::null!\n");
364 stepControl_ = stepControl;
365}
366
367template<class Scalar>
368RCP<StepControlStrategyBase<Scalar> >
373
374template<class Scalar>
375RCP<const StepControlStrategyBase<Scalar> >
377{
378 return(stepControl_);
379}
380
381template<class Scalar>
383 StepSizeType stepSizeType)
384{
385
386 using Teuchos::as;
387 using Teuchos::incrVerbLevel;
388 typedef Teuchos::ScalarTraits<Scalar> ST;
389 typedef Thyra::NonlinearSolverBase<Scalar> NSB;
390 typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
391
392 RCP<Teuchos::FancyOStream> out = this->getOStream();
393 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
394 Teuchos::OSTab ostab(out,1,"BES::takeStep");
395 VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
396
397 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
398 *out << "\nEntering "
399 << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
400 << "::takeStep("<<dt<<","<<toString(stepSizeType)<<") ...\n";
401 }
402
403 initialize_();
404
405 if(!neModel_.get()) {
406 neModel_ =Teuchos::rcp(new Rythmos::SingleResidualModelEvaluator<Scalar>());
407 }
408
409 dt_ = dt;
410
411 if (dt <= ST::zero()) {
412 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
413 *out << "\nThe arguments to takeStep are not valid for "
414 << "BackwardEulerStepper at this time.\n"
415 << " dt = " << dt << "\n"
416 << "BackwardEulerStepper requires positive dt.\n" << std::endl;
417 return(Scalar(-ST::one()));
418 }
419 if ((stepSizeType == STEP_TYPE_VARIABLE) && (stepControl_ == Teuchos::null)) {
420 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
421 *out << "\nFor 'variable' time stepping (step size adjustable by "
422 << "BackwardEulerStepper), BackwardEulerStepper requires "
423 << "Step-Control Strategy.\n"
424 << " stepType = " << toString(stepSizeType) << "\n" << std::endl;
425 return(Scalar(-ST::one()));
426 }
427
428 stepControl_->setRequestedStepSize(*this,dt_,stepSizeType);
429 AttemptedStepStatusFlag status;
430 bool stepPass = false;
431 while (1) {
432
433 stepControl_->nextStepSize(*this,&dt_,&stepSizeType,NULL);
434 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
435 *out << "\nrequested dt = " << dt
436 << "\ncurrent dt = " << dt_ << "\n";
437 }
438
439 // Setup the nonlinear equations:
440 //
441 // f( (1/dt)* x + (-1/dt)*x_old), x, t ) = 0
442
443 V_StV( scaled_x_old_.ptr(), Scalar(-ST::one()/dt_), *x_old_ );
444 t_old_ = t_;
445 neModel_->initializeSingleResidualModel(
446 model_, basePoint_,
447 Scalar(ST::one()/dt_), scaled_x_old_,
448 ST::one(), Teuchos::null,
449 t_old_+dt_,
450 Teuchos::null
451 );
452 if( solver_->getModel().get() != neModel_.get() ) {
453 solver_->setModel(neModel_);
454 }
455 // 2007/05/18: rabartl: ToDo: Above, set the stream and the verbosity level
456 // on solver_ so that we an see what it is doing!
457
458 obtainPredictor_();
459
460 // Solve the implicit nonlinear system to a tolerance of ???
461
462 if ( as<int>(verbLevel) > as<int>(Teuchos::VERB_LOW) ) {
463 *out << "\nSolving the implicit backward-Euler timestep equation ...\n";
464 }
465
466 Thyra::SolveStatus<Scalar> neSolveStatus =
467 solver_->solve(&*x_, NULL, &*dx_);
468
469 // In the above solve, on input *x_ is the initial guess that comes from
470 // the predictor. On output, *x_ is the converged timestep solution and
471 // *dx_ is the difference computed from the intial guess in *x_ to the
472 // final solved value of *x_. This is needed for basic numerical stability.
473
474 if ( as<int>(verbLevel) > as<int>(Teuchos::VERB_LOW) ) {
475 *out << "\nOutput status of nonlinear solve:\n" << neSolveStatus;
476 }
477
478 // 2007/05/18: rabartl: ToDo: Above, get the solve status from the above
479 // solve and at least print warning message if the solve fails! Actually,
480 // you should most likely thrown an exception if the solve fails or return
481 // false if appropriate
482
483 if (neSolveStatus.solveStatus == Thyra::SOLVE_STATUS_CONVERGED) {
484 newtonConvergenceStatus_ = 0;
485 }
486 else {
487 newtonConvergenceStatus_ = -1;
488 }
489
490 stepControl_->setCorrection(*this,x_,dx_,newtonConvergenceStatus_);
491
492 stepPass = stepControl_->acceptStep(*this,NULL);
493
494 if (!stepPass) { // stepPass = false
495 status = stepControl_->rejectStep(*this);
496
497 if (status != PREDICT_AGAIN)
498 break;
499 } else { // stepPass = true
500 break;
501 }
502 }
503
504 // Update the step
505
506 if (stepPass) {
507 V_V( x_dot_old_.ptr(), *x_dot_ );
508 // x_dot = (1/dt)*x - (1/dt)*x_old
509 V_StVpStV( x_dot_.ptr(), Scalar(ST::one()/dt_), *x_,
510 Scalar(-ST::one()/dt_), *x_old_ );
511 V_V( x_old_.ptr(), *x_ );
512 t_ += dt_;
513 numSteps_++;
514 stepControl_->completeStep(*this);
515 } else {
516 // Complete failure. Return to Integrator with bad step size.
517 dt_ = Scalar(-ST::one());
518 }
519
520 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
521 *out << "\nt_old_ = " << t_old_ << std::endl;
522 *out << "\nt_ = " << t_ << std::endl;
523 }
524
525#ifdef HAVE_RYTHMOS_DEBUG
526 // 04/14/09 tscoffe: This code should be moved to StepperValidator
527
528 if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
529 *out << "\nChecking to make sure that solution and "
530 << "the interpolated solution are the same! ...\n";
531
532 {
533
534 typedef ScalarTraits<ScalarMag> SMT;
535
536 Teuchos::OSTab tab(out);
537
538 const StepStatus<Scalar> stepStatus = this->getStepStatus();
539
540 RCP<const Thyra::VectorBase<Scalar> >
541 x = stepStatus.solution,
542 xdot = stepStatus.solutionDot;
543
544 Array<Scalar> time_vec = Teuchos::tuple(stepStatus.time);
545 Array<RCP<const Thyra::VectorBase<Scalar> > > x_vec, xdot_vec;
546 this->getPoints(time_vec,&x_vec,&xdot_vec,0);
547
548 RCP<const Thyra::VectorBase<Scalar> >
549 x_interp = x_vec[0],
550 xdot_interp = xdot_vec[0];
551
552 TEUCHOS_TEST_FOR_EXCEPT(
553 !Thyra::testRelNormDiffErr(
554 "x", *x, "x_interp", *x_interp,
555 "2*epsilon", ScalarMag(100.0*SMT::eps()),
556 "2*epsilon", ScalarMag(100.0*SMT::eps()),
557 includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
558 )
559 );
560
561 TEUCHOS_TEST_FOR_EXCEPT(
562 !Thyra::testRelNormDiffErr(
563 "xdot", *xdot, "xdot_interp", *xdot_interp,
564 "2*epsilon", ScalarMag(100.0*SMT::eps()),
565 "2*epsilon", ScalarMag(100.0*SMT::eps()),
566 includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
567 )
568 );
569
570 }
571
572 // 2007/07/25: rabartl: ToDo: Move the above test into a helper function so
573 // that it can be used from lots of different places!
574
575#endif // HAVE_RYTHMOS_DEBUG
576
577 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
578 *out << "\nLeaving "
579 << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
580 << "::takeStep("<<dt_<<","<<toString(stepSizeType)<<") ...\n";
581 }
582
583 return(dt_);
584}
585
586
587template<class Scalar>
588const StepStatus<Scalar> BackwardEulerStepper<Scalar>::getStepStatus() const
589{
590
591 // typedef Teuchos::ScalarTraits<Scalar> ST; // unused
592
593 StepStatus<Scalar> stepStatus; // Defaults to unknown status
594
595 if (!isInitialized_) {
596 stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
597 }
598 else if (numSteps_ > 0) {
599 stepStatus.stepStatus = STEP_STATUS_CONVERGED;
600 }
601 // else unknown
602
603 stepStatus.stepSize = dt_;
604 stepStatus.order = 1;
605 stepStatus.time = t_;
606 stepStatus.solution = x_;
607 stepStatus.solutionDot = x_dot_;
608
609 return(stepStatus);
610
611}
612
613
614// Overridden from InterpolationBufferBase
615
616
617template<class Scalar>
618RCP<const Thyra::VectorSpaceBase<Scalar> >
620{
621 return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
622}
623
624
625template<class Scalar>
627 const Array<Scalar>& time_vec,
628 const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
629 const Array<RCP<const Thyra::VectorBase<Scalar> > >& /* xdot_vec */
630 )
631{
632
633 typedef Teuchos::ScalarTraits<Scalar> ST;
634 using Teuchos::as;
635
636 initialize_();
637
638#ifdef HAVE_RYTHMOS_DEBUG
639 TEUCHOS_TEST_FOR_EXCEPTION(
640 time_vec.size() == 0, std::logic_error,
641 "Error, addPoints called with an empty time_vec array!\n");
642#endif // HAVE_RYTHMOS_DEBUG
643
644 RCP<Teuchos::FancyOStream> out = this->getOStream();
645 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
646 Teuchos::OSTab ostab(out,1,"BES::setPoints");
647
648 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
649 *out << "time_vec = " << std::endl;
650 for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
651 *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
652 }
653 }
654 else if (time_vec.size() == 1) {
655 int n = 0;
656 t_ = time_vec[n];
657 t_old_ = t_;
658 Thyra::V_V(x_.ptr(),*x_vec[n]);
659 Thyra::V_V(scaled_x_old_.ptr(),*x_);
660 }
661 else {
662 int n = time_vec.size()-1;
663 int nm1 = time_vec.size()-2;
664 t_ = time_vec[n];
665 t_old_ = time_vec[nm1];
666 Thyra::V_V(x_.ptr(),*x_vec[n]);
667 Scalar dt = t_ - t_old_;
668 Thyra::V_StV(scaled_x_old_.ptr(),Scalar(-ST::one()/dt),*x_vec[nm1]);
669 }
670}
671
672
673template<class Scalar>
675{
676 if ( !isInitialized_ && haveInitialCondition_ )
677 return timeRange<Scalar>(t_,t_);
678 if ( !isInitialized_ && !haveInitialCondition_ )
679 return invalidTimeRange<Scalar>();
680 return timeRange<Scalar>(t_old_,t_);
681}
682
683
684template<class Scalar>
686 const Array<Scalar>& time_vec,
687 Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
688 Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
689 Array<ScalarMag>* accuracy_vec
690 ) const
691{
692 typedef Teuchos::ScalarTraits<Scalar> ST;
693 using Teuchos::constOptInArg;
694 using Teuchos::ptr;
695#ifdef HAVE_RYTHMOS_DEBUG
696 TEUCHOS_ASSERT(haveInitialCondition_);
697#endif // HAVE_RYTHMOS_DEBUG
698 RCP<Thyra::VectorBase<Scalar> > x_temp = x_;
699 if (compareTimeValues(t_old_,t_)!=0) {
700 Scalar dt = t_ - t_old_;
701 x_temp = scaled_x_old_->clone_v();
702 Thyra::Vt_S(x_temp.ptr(),Scalar(-ST::one()*dt)); // undo the scaling
703 }
704 defaultGetPoints<Scalar>(
705 t_old_, constOptInArg(*x_temp), constOptInArg(*x_dot_old_),
706 t_, constOptInArg(*x_), constOptInArg(*x_dot_),
707 time_vec, ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
708 ptr(interpolator_.get())
709 );
710
711 /*
712 using Teuchos::as;
713 typedef Teuchos::ScalarTraits<Scalar> ST;
714 typedef typename ST::magnitudeType ScalarMag;
715 typedef Teuchos::ScalarTraits<ScalarMag> SMT;
716 typename DataStore<Scalar>::DataStoreVector_t ds_nodes;
717 typename DataStore<Scalar>::DataStoreVector_t ds_out;
718
719#ifdef HAVE_RYTHMOS_DEBUG
720 TEUCHOS_TEST_FOR_EXCEPT(!haveInitialCondition_);
721 TEUCHOS_TEST_FOR_EXCEPT( 0 == x_vec );
722#endif
723
724 RCP<Teuchos::FancyOStream> out = this->getOStream();
725 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
726 Teuchos::OSTab ostab(out,1,"BES::getPoints");
727 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
728 *out
729 << "\nEntering " << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
730 << "::getPoints(...) ...\n";
731 }
732 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
733 for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
734 *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
735 }
736 *out << "I can interpolate in the interval [" << t_old_ << "," << t_ << "]." << std::endl;
737 }
738 if (t_old_ != t_) {
739 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
740 *out << "Passing two points to interpolator: " << t_old_ << " and " << t_ << std::endl;
741 }
742 DataStore<Scalar> ds_temp;
743 Scalar dt = t_ - t_old_;
744#ifdef HAVE_RYTHMOS_DEBUG
745 TEUCHOS_TEST_FOR_EXCEPT(
746 !Thyra::testRelErr(
747 "dt", dt, "dt_", dt_,
748 "1e+4*epsilon", ScalarMag(1e+4*SMT::eps()),
749 "1e+2*epsilon", ScalarMag(1e+2*SMT::eps()),
750 as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) ? out.get() : 0
751 )
752 );
753#endif
754 RCP<Thyra::VectorBase<Scalar> >
755 x_temp = scaled_x_old_->clone_v();
756 Thyra::Vt_S(&*x_temp,Scalar(-ST::one()*dt)); // undo the scaling
757 ds_temp.time = t_old_;
758 ds_temp.x = x_temp;
759 ds_temp.xdot = x_dot_old_;
760 ds_temp.accuracy = ScalarMag(dt);
761 ds_nodes.push_back(ds_temp);
762 ds_temp.time = t_;
763 ds_temp.x = x_;
764 ds_temp.xdot = x_dot_;
765 ds_temp.accuracy = ScalarMag(dt);
766 ds_nodes.push_back(ds_temp);
767 }
768 else {
769 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
770 *out << "Passing one point to interpolator: " << t_ << std::endl;
771 }
772 DataStore<Scalar> ds_temp;
773 ds_temp.time = t_;
774 ds_temp.x = x_;
775 ds_temp.xdot = x_dot_;
776 ds_temp.accuracy = ScalarMag(ST::zero());
777 ds_nodes.push_back(ds_temp);
778 }
779 interpolate<Scalar>(*interpolator_,rcp(&ds_nodes,false),time_vec,&ds_out);
780 Array<Scalar> time_out;
781 dataStoreVectorToVector(ds_out,&time_out,x_vec,xdot_vec,accuracy_vec);
782 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
783 *out << "Passing out the interpolated values:" << std::endl;
784 for (int i=0; i<Teuchos::as<int>(time_out.size()) ; ++i) {
785 *out << "time[" << i << "] = " << time_out[i] << std::endl;
786 *out << "x_vec[" << i << "] = " << std::endl;
787 (*x_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
788 if (xdot_vec) {
789 if ( (*xdot_vec)[i] == Teuchos::null) {
790 *out << "xdot_vec[" << i << "] = Teuchos::null" << std::endl;
791 }
792 else {
793 *out << "xdot_vec[" << i << "] = " << std::endl;
794 (*xdot_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
795 }
796 }
797 if(accuracy_vec) {
798 *out << "accuracy[" << i << "] = " << (*accuracy_vec)[i] << std::endl;
799 }
800 }
801 }
802 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
803 *out
804 << "Leaving " << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
805 << "::getPoints(...) ...\n";
806 }
807 */
808
809}
810
811
812template<class Scalar>
813void BackwardEulerStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
814{
815 using Teuchos::as;
816
817 TEUCHOS_ASSERT( time_vec != NULL );
818
819 time_vec->clear();
820
821 if (!haveInitialCondition_) {
822 return;
823 }
824
825 time_vec->push_back(t_old_);
826
827 if (numSteps_ > 0) {
828 time_vec->push_back(t_);
829 }
830
831 RCP<Teuchos::FancyOStream> out = this->getOStream();
832 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
833 Teuchos::OSTab ostab(out,1,"BES::getNodes");
834 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
835 *out << this->description() << std::endl;
836 for (int i=0 ; i<Teuchos::as<int>(time_vec->size()) ; ++i) {
837 *out << "time_vec[" << i << "] = " << (*time_vec)[i] << std::endl;
838 }
839 }
840}
841
842
843template<class Scalar>
844void BackwardEulerStepper<Scalar>::removeNodes(Array<Scalar>& time_vec)
845{
846 initialize_();
847 using Teuchos::as;
848 RCP<Teuchos::FancyOStream> out = this->getOStream();
849 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
850 Teuchos::OSTab ostab(out,1,"BES::removeNodes");
851 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
852 *out << "time_vec = " << std::endl;
853 for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
854 *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
855 }
856 }
857 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for BackwardEulerStepper at this time.\n");
858 // TODO:
859 // if any time in time_vec matches t_ or t_old_, then do the following:
860 // remove t_old_: set t_old_ = t_ and set scaled_x_old_ = x_
861 // remove t_: set t_ = t_old_ and set x_ = -dt*scaled_x_old_
862}
863
864
865template<class Scalar>
867{
868 return(1);
869}
870
871
872// Overridden from Teuchos::ParameterListAcceptor
873
874
875template <class Scalar>
877 RCP<Teuchos::ParameterList> const& paramList
878 )
879{
880 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
881 paramList->validateParametersAndSetDefaults(*this->getValidParameters());
882 parameterList_ = paramList;
883 Teuchos::readVerboseObjectSublist(&*parameterList_,this);
884}
885
886
887template <class Scalar>
888RCP<Teuchos::ParameterList>
890{
891 return(parameterList_);
892}
893
894
895template <class Scalar>
896RCP<Teuchos::ParameterList>
898{
899 RCP<Teuchos::ParameterList>
900 temp_param_list = parameterList_;
901 parameterList_ = Teuchos::null;
902 return(temp_param_list);
903}
904
905
906template<class Scalar>
907RCP<const Teuchos::ParameterList>
909{
910 using Teuchos::ParameterList;
911 static RCP<const ParameterList> validPL;
912 if (is_null(validPL)) {
913 RCP<ParameterList> pl = Teuchos::parameterList();
914 // This line is required to pass StepperValidator UnitTest!
915 pl->sublist(RythmosStepControlSettings_name);
916 Teuchos::setupVerboseObjectSublist(&*pl);
917 validPL = pl;
918 }
919 return validPL;
920}
921
922
923// Overridden from Teuchos::Describable
924
925
926template<class Scalar>
928 Teuchos::FancyOStream &out,
929 const Teuchos::EVerbosityLevel verbLevel
930 ) const
931{
932 using Teuchos::as;
933 Teuchos::OSTab tab(out);
934 if (!isInitialized_) {
935 out << this->description() << " : This stepper is not initialized yet"
936 << std::endl;
937 return;
938 }
939 if (
940 as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
941 ||
942 as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
943 )
944 {
945 out << this->description() << "::describe:" << std::endl;
946 out << "model = " << model_->description() << std::endl;
947 out << "solver = " << solver_->description() << std::endl;
948 if (neModel_ == Teuchos::null) {
949 out << "neModel = Teuchos::null" << std::endl;
950 } else {
951 out << "neModel = " << neModel_->description() << std::endl;
952 }
953 }
954 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
955 out << "t_ = " << t_ << std::endl;
956 out << "t_old_ = " << t_old_ << std::endl;
957 }
958 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
959 }
960 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
961 out << "model_ = " << std::endl;
962 model_->describe(out,verbLevel);
963 out << "solver_ = " << std::endl;
964 solver_->describe(out,verbLevel);
965 if (neModel_ == Teuchos::null) {
966 out << "neModel = Teuchos::null" << std::endl;
967 } else {
968 out << "neModel = " << std::endl;
969 neModel_->describe(out,verbLevel);
970 }
971 out << "x_ = " << std::endl;
972 x_->describe(out,verbLevel);
973 out << "scaled_x_old_ = " << std::endl;
974 scaled_x_old_->describe(out,verbLevel);
975 }
976}
977
978
979// private
980
981
982template <class Scalar>
984{
985
986 if (isInitialized_) return;
987
988 TEUCHOS_TEST_FOR_EXCEPT(is_null(model_));
989 TEUCHOS_TEST_FOR_EXCEPT(is_null(solver_));
990 TEUCHOS_TEST_FOR_EXCEPT(!haveInitialCondition_);
991
992 // Initialize Parameter List if none provided.
993 if (parameterList_ == Teuchos::null) {
994 RCP<Teuchos::ParameterList> emptyParameterList =
995 Teuchos::rcp(new Teuchos::ParameterList);
996 this->setParameterList(emptyParameterList);
997 }
998
999 // Initialize StepControl
1000 if (stepControl_ == Teuchos::null) {
1001 RCP<StepControlStrategyBase<Scalar> > stepControlStrategy =
1002 Teuchos::rcp(new FixedStepControlStrategy<Scalar>());
1003 RCP<Teuchos::ParameterList> stepControlPL =
1004 Teuchos::sublist(parameterList_, RythmosStepControlSettings_name);
1005
1006 stepControlStrategy->setParameterList(stepControlPL);
1007 this->setStepControlStrategy(stepControlStrategy);
1008 }
1009 stepControl_->initialize(*this);
1010 stepControl_->setOStream(this->getOStream());
1011 stepControl_->setVerbLevel(this->getVerbLevel());
1012
1013 //maxOrder_ = stepControl_->getMaxOrder(); // maximum order
1014
1015#ifdef HAVE_RYTHMOS_DEBUG
1016 THYRA_ASSERT_VEC_SPACES(
1017 "Rythmos::BackwardEulerStepper::initialize_(...)",
1018 *x_->space(), *model_->get_x_space() );
1019#endif // HAVE_RYTHMOS_DEBUG
1020
1021 if ( is_null(interpolator_) ) {
1022 // If an interpolator has not been explicitly set, then just create
1023 // a default linear interpolator.
1024 interpolator_ = linearInterpolator<Scalar>();
1025 // 2007/05/18: rabartl: ToDo: Replace this with a Hermite interplator
1026 // when it is implementated!
1027 }
1028 isInitialized_ = true;
1029}
1030
1031
1032template<class Scalar>
1033RCP<const MomentoBase<Scalar> >
1035{
1036 RCP<BackwardEulerStepperMomento<Scalar> > momento =
1037 Teuchos::rcp(new BackwardEulerStepperMomento<Scalar>());
1038 momento->set_scaled_x_old(scaled_x_old_);
1039 momento->set_x_old(x_old_);
1040 momento->set_x_dot_old(x_dot_old_);
1041 momento->set_x(x_);
1042 momento->set_x_dot(x_dot_);
1043 momento->set_dx(dx_);
1044 momento->set_t(t_);
1045 momento->set_t_old(t_old_);
1046 momento->set_dt(dt_);
1047 momento->set_numSteps(numSteps_);
1048 momento->set_newtonConvergenceStatus(newtonConvergenceStatus_);
1049 momento->set_isInitialized(isInitialized_);
1050 momento->set_haveInitialCondition(haveInitialCondition_);
1051 momento->set_parameterList(parameterList_);
1052 momento->set_basePoint(basePoint_);
1053 momento->set_neModel(neModel_);
1054 momento->set_interpolator(interpolator_);
1055 momento->set_stepControl(stepControl_);
1056 return momento;
1057}
1058
1059template<class Scalar>
1061 const Ptr<const MomentoBase<Scalar> >& momentoPtr,
1062 const RCP<Thyra::ModelEvaluator<Scalar> >& model,
1063 const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
1064 )
1065{
1066 Ptr<const BackwardEulerStepperMomento<Scalar> > beMomentoPtr =
1067 Teuchos::ptr_dynamic_cast<const BackwardEulerStepperMomento<Scalar> >
1068 (momentoPtr,true);
1069 const BackwardEulerStepperMomento<Scalar>& beMomento = *beMomentoPtr;
1070 model_ = model;
1071 solver_ = solver;
1072 scaled_x_old_ = beMomento.get_scaled_x_old();
1073 x_old_ = beMomento.get_x_old();
1074 x_dot_old_ = beMomento.get_x_dot_old();
1075 x_ = beMomento.get_x();
1076 x_dot_ = beMomento.get_x_dot();
1077 dx_ = beMomento.get_dx();
1078 t_ = beMomento.get_t();
1079 t_old_ = beMomento.get_t_old();
1080 dt_ = beMomento.get_dt();
1081 numSteps_ = beMomento.get_numSteps();
1082 newtonConvergenceStatus_ = beMomento.get_newtonConvergenceStatus();
1083 isInitialized_ = beMomento.get_isInitialized();
1084 haveInitialCondition_ = beMomento.get_haveInitialCondition();
1085 parameterList_ = beMomento.get_parameterList();
1086 basePoint_ = beMomento.get_basePoint();
1087 neModel_ = beMomento.get_neModel();
1088 interpolator_ = beMomento.get_interpolator();
1089 stepControl_ = beMomento.get_stepControl();
1090 this->checkConsistentState_();
1091}
1092
1093template<class Scalar>
1095{
1096 if (isInitialized_) {
1097 TEUCHOS_ASSERT( !Teuchos::is_null(model_) );
1098 TEUCHOS_ASSERT( !Teuchos::is_null(solver_) );
1099 TEUCHOS_ASSERT( haveInitialCondition_ );
1100 TEUCHOS_ASSERT( !Teuchos::is_null(interpolator_) );
1101 TEUCHOS_ASSERT( !Teuchos::is_null(stepControl_) );
1102 }
1103 if (haveInitialCondition_) {
1104 // basePoint_ should be defined
1105 typedef Teuchos::ScalarTraits<Scalar> ST;
1106 TEUCHOS_ASSERT( !ST::isnaninf(t_) );
1107 TEUCHOS_ASSERT( !ST::isnaninf(t_old_) );
1108 TEUCHOS_ASSERT( !Teuchos::is_null(scaled_x_old_) );
1109 TEUCHOS_ASSERT( !Teuchos::is_null(x_dot_old_) );
1110 TEUCHOS_ASSERT( !Teuchos::is_null(x_) );
1111 TEUCHOS_ASSERT( !Teuchos::is_null(x_dot_) );
1112 TEUCHOS_ASSERT( !Teuchos::is_null(x_old_) );
1113 TEUCHOS_ASSERT( !Teuchos::is_null(dx_) );
1114 TEUCHOS_ASSERT( t_ >= basePoint_.get_t() );
1115 TEUCHOS_ASSERT( t_old_ >= basePoint_.get_t() );
1116 }
1117 if (numSteps_ > 0) {
1118 TEUCHOS_ASSERT(isInitialized_);
1119 }
1120}
1121
1122template<class Scalar>
1123void BackwardEulerStepper<Scalar>::obtainPredictor_()
1124{
1125 using Teuchos::as;
1126 // typedef Teuchos::ScalarTraits<Scalar> ST; // unused
1127
1128 if (!isInitialized_) return;
1129
1130 RCP<Teuchos::FancyOStream> out = this->getOStream();
1131 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
1132 Teuchos::OSTab ostab(out,1,"obtainPredictor_");
1133
1134 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1135 *out << "Before predictor x_ = " << std::endl;
1136 x_->describe(*out,verbLevel);
1137 }
1138 // evaluate predictor -- basic Forward Euler
1139 // x_ = dt_*x_dot_old_ + x_old_
1140 V_StVpV(x_.ptr(), dt_, *x_dot_old_, *x_old_);
1141
1142 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1143 *out << "After predictor x_ = " << std::endl;
1144 x_->describe(*out,verbLevel);
1145 }
1146}
1147
1148//
1149// Explicit Instantiation macro
1150//
1151// Must be expanded from within the Rythmos namespace!
1152//
1153
1154#define RYTHMOS_BACKWARD_EULER_STEPPER_INSTANT(SCALAR) \
1155 \
1156 template class BackwardEulerStepper< SCALAR >; \
1157 \
1158 template RCP< BackwardEulerStepper< SCALAR > > \
1159 backwardEulerStepper( \
1160 const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \
1161 const RCP<Thyra::NonlinearSolverBase< SCALAR > >& solver \
1162 ); \
1163 template RCP< BackwardEulerStepper< SCALAR > > \
1164 backwardEulerStepper();
1165
1166
1167
1168
1169} // namespace Rythmos
1170
1171#endif //Rythmos_BACKWARD_EULER_STEPPER_DEF_H
Simple concrete stepper subclass implementing an implicit backward Euler method.
Teuchos::ScalarTraits< Scalar >::magnitudeType ScalarMag
RCP< const Thyra::VectorBase< Scalar > > get_x(const InterpolationBufferBase< Scalar > &interpBuffer, const Scalar &t)
Get a single point x(t) from an interpolation buffer.