Rythmos - Transient Integration for Differential Equations  Version of the Day
 All Classes Functions Variables Typedefs Pages
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 
48 namespace Rythmos {
49 
50 
51 template<class Scalar>
52 RCP<BackwardEulerStepper<Scalar> >
53 backwardEulerStepper(
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 
61 template<class Scalar>
62 RCP<BackwardEulerStepper<Scalar> >
63 backwardEulerStepper()
64 {
65  return Teuchos::rcp(new BackwardEulerStepper<Scalar>());
66 }
67 
68 
69 // ////////////////////////////
70 // Defintions
71 
72 
73 // Constructors, intializers, Misc.
74 
75 
76 template<class Scalar>
78 {
79  typedef Teuchos::ScalarTraits<Scalar> ST;
80  this->defaultInitializeAll_();
81  numSteps_ = 0;
82  dt_ = ST::zero();
83 }
84 
85 
86 template<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 
100 template<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 
128 template<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 
140 template<class Scalar>
141 RCP<InterpolatorBase<Scalar> >
143 {
144  return interpolator_;
145 }
146 
147 template<class Scalar>
148 RCP<const InterpolatorBase<Scalar> >
150 {
151  return interpolator_;
152 }
153 
154 template<class Scalar>
155 RCP<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 
168 template<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 
193 template<class Scalar>
194 RCP<Thyra::NonlinearSolverBase<Scalar> >
196 {
197  return solver_;
198 }
199 
200 
201 template<class Scalar>
202 RCP<const Thyra::NonlinearSolverBase<Scalar> >
204 {
205  return solver_;
206 }
207 
208 
209 // Overridden from StepperBase
210 
211 
212 template<class Scalar>
214 {
215  return true;
216 }
217 
218 
219 template<class Scalar>
220 RCP<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_
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 
253 template<class Scalar>
255 {
256  return true;
257 }
258 
259 template<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 
277 template<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 
285 template<class Scalar>
286 RCP<const Thyra::ModelEvaluator<Scalar> >
288 {
289  return model_;
290 }
291 
292 
293 template<class Scalar>
294 RCP<Thyra::ModelEvaluator<Scalar> >
296 {
297  return Teuchos::null;
298 }
299 
300 
301 template<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 
351 template<class Scalar>
352 Thyra::ModelEvaluatorBase::InArgs<Scalar>
354 {
355  return basePoint_;
356 }
357 
358 template<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 
367 template<class Scalar>
368 RCP<StepControlStrategyBase<Scalar> >
370 {
371  return(stepControl_);
372 }
373 
374 template<class Scalar>
375 RCP<const StepControlStrategyBase<Scalar> >
377 {
378  return(stepControl_);
379 }
380 
381 template<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 
587 template<class Scalar>
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 
617 template<class Scalar>
618 RCP<const Thyra::VectorSpaceBase<Scalar> >
620 {
621  return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
622 }
623 
624 
625 template<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 
673 template<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 
684 template<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 
812 template<class Scalar>
813 void 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 
843 template<class Scalar>
844 void 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 
865 template<class Scalar>
867 {
868  return(1);
869 }
870 
871 
872 // Overridden from Teuchos::ParameterListAcceptor
873 
874 
875 template <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 
887 template <class Scalar>
888 RCP<Teuchos::ParameterList>
890 {
891  return(parameterList_);
892 }
893 
894 
895 template <class Scalar>
896 RCP<Teuchos::ParameterList>
898 {
899  RCP<Teuchos::ParameterList>
900  temp_param_list = parameterList_;
901  parameterList_ = Teuchos::null;
902  return(temp_param_list);
903 }
904 
905 
906 template<class Scalar>
907 RCP<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 
926 template<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 
982 template <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 
1032 template<class Scalar>
1033 RCP<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 
1059 template<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 
1093 template<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 
1122 template<class Scalar>
1123 void 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.
Base strategy class for interpolation functionality.
RCP< const InterpolatorBase< Scalar > > getInterpolator() const
void setSolver(const RCP< Thyra::NonlinearSolverBase< Scalar > > &solver)
void setStepControlStrategy(const RCP< StepControlStrategyBase< Scalar > > &stepControlStrategy)
void getNodes(Array< Scalar > *time_vec) const
RCP< Teuchos::ParameterList > unsetParameterList()
void setInitialCondition(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &initialCondition)
void setModel(const RCP< const Thyra::ModelEvaluator< Scalar > > &model)
void removeNodes(Array< Scalar > &time_vec)
Scalar takeStep(Scalar dt, StepSizeType flag)
void addPoints(const Array< Scalar > &time_vec, const Array< RCP< const Thyra::VectorBase< Scalar > > > &x_vec, const Array< RCP< const Thyra::VectorBase< Scalar > > > &xdot_vec)
RCP< Teuchos::ParameterList > getNonconstParameterList()
RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
RCP< StepControlStrategyBase< Scalar > > getNonconstStepControlStrategy()
The member functions in the StepControlStrategyBase move you between these states in the following fa...
void setNonconstModel(const RCP< Thyra::ModelEvaluator< Scalar > > &model)
RCP< InterpolatorBase< Scalar > > unSetInterpolator()
Thyra::ModelEvaluatorBase::InArgs< Scalar > getInitialCondition() const
Teuchos::ScalarTraits< Scalar >::magnitudeType ScalarMag
void getPoints(const Array< Scalar > &time_vec, Array< RCP< const Thyra::VectorBase< Scalar > > > *x_vec, Array< RCP< const Thyra::VectorBase< Scalar > > > *xdot_vec, Array< ScalarMag > *accuracy_vec) const
RCP< const Thyra::VectorBase< Scalar > > solutionDot
void setMomento(const Ptr< const MomentoBase< Scalar > > &momentoPtr, const RCP< Thyra::ModelEvaluator< Scalar > > &model, const RCP< Thyra::NonlinearSolverBase< Scalar > > &solver)
Set momento object for use in restarts.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Concrete momento class for the BackwardEulerStepper.
RCP< Thyra::NonlinearSolverBase< Scalar > > getNonconstSolver()
RCP< const MomentoBase< Scalar > > getMomento() const
Get momento object for use in restarts.
Decorator subclass for a steady-state version of a DAE for single-residual time stepper methods...
RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
RCP< const StepControlStrategyBase< Scalar > > getStepControlStrategy() const
RCP< const Teuchos::ParameterList > getValidParameters() const
RCP< StepperBase< Scalar > > cloneStepperAlgorithm() const
Creates copies of all internal data (including the parameter list) except the model which is assumed ...
const StepStatus< Scalar > getStepStatus() const
RCP< const Thyra::VectorBase< Scalar > > solution
RCP< const Thyra::NonlinearSolverBase< Scalar > > getSolver() const
RCP< InterpolatorBase< Scalar > > getNonconstInterpolator()
void setParameterList(RCP< Teuchos::ParameterList > const &paramList)
RCP< Thyra::ModelEvaluator< Scalar > > getNonconstModel()
Base class for a momento object.
void setInterpolator(const RCP< InterpolatorBase< Scalar > > &interpolator)