Rythmos - Transient Integration for Differential Equations  Version of the Day
Rythmos_ThetaStepper_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_THETA_STEPPER_DEF_H
30 #define Rythmos_THETA_STEPPER_DEF_H
31 
32 #include "Rythmos_ConfigDefs.h"
33 #ifdef HAVE_RYTHMOS_EXPERIMENTAL
34 
35 #include "Rythmos_ThetaStepper_decl.hpp"
36 
37 namespace Rythmos {
38 
39 
44 template<class Scalar>
45 RCP<ThetaStepper<Scalar> >
46 thetaStepper(
47  const RCP<Thyra::ModelEvaluator<Scalar> >& model,
48  const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
49  RCP<Teuchos::ParameterList>& parameterList
50  )
51 {
52  Teuchos::RCP<ThetaStepper<Scalar> > stepper =
53  Teuchos::rcp(new ThetaStepper<Scalar>());
54  stepper->setParameterList(parameterList);
55  stepper->setModel(model);
56  stepper->setSolver(solver);
57 
58  return stepper;
59 }
60 
61 // ////////////////////////////
62 // Defintions
63 
64 
65 // Constructors, intializers, Misc.
66 
67 
68 template<class Scalar>
70 {
71  typedef Teuchos::ScalarTraits<Scalar> ST;
72  this->defaultInitializeAll_();
73  Scalar zero = ST::zero();
74  t_ = -ST::one();
75  t_old_ = zero;
76  dt_ = zero;
77  dt_old_ = zero;
78  numSteps_ = 0;
79 }
80 
81 template<class Scalar>
82 void ThetaStepper<Scalar>::defaultInitializeAll_()
83 {
84  typedef Teuchos::ScalarTraits<Scalar> ST;
85  isInitialized_ = false;
86  haveInitialCondition_ = false;
87  model_ = Teuchos::null;
88  solver_ = Teuchos::null;
89  //basePoint_;
90  x_ = Teuchos::null;
91  x_old_ = Teuchos::null;
92  x_pre_ = Teuchos::null;
93  x_dot_ = Teuchos::null;
94  x_dot_old_ = Teuchos::null;
95  x_dot_really_old_ = Teuchos::null;
96  x_dot_base_ = Teuchos::null;
97  t_ = ST::nan();
98  t_old_ = ST::nan();
99  dt_ = ST::nan();
100  dt_old_ = ST::nan();
101  numSteps_ = -1;
102  thetaStepperType_ = INVALID_THETA_STEPPER_TYPE;
103  theta_ = ST::nan();
104  neModel_ = Teuchos::null;
105  parameterList_ = Teuchos::null;
106  interpolator_ = Teuchos::null;
107  predictor_corrector_begin_after_step_ = -1;
108  default_predictor_order_ = -1;
109 }
110 
111 template<class Scalar>
113 {
114  return true;
115 }
116 
117 
118 template<class Scalar>
120  const RCP<InterpolatorBase<Scalar> >& interpolator
121  )
122 {
123 #ifdef HAVE_RYTHMOS_DEBUG
124  TEUCHOS_TEST_FOR_EXCEPT(is_null(interpolator));
125 #endif
126  interpolator_ = interpolator;
127  isInitialized_ = false;
128 }
129 
130 template<class Scalar>
131 RCP<InterpolatorBase<Scalar> >
133 {
134  return interpolator_;
135 }
136 
137 template<class Scalar>
138 RCP<const InterpolatorBase<Scalar> >
140 {
141  return interpolator_;
142 }
143 
144 
145 template<class Scalar>
146 RCP<InterpolatorBase<Scalar> >
148 {
149  RCP<InterpolatorBase<Scalar> > temp_interpolator = interpolator_;
150  interpolator_ = Teuchos::null;
151  return(temp_interpolator);
152  isInitialized_ = false;
153 }
154 
155 
156 // Overridden from SolverAcceptingStepperBase
157 
158 
159 template<class Scalar>
161  const RCP<Thyra::NonlinearSolverBase<Scalar> > &solver
162  )
163 {
164  using Teuchos::as;
165 
166  TEUCHOS_TEST_FOR_EXCEPTION(solver == Teuchos::null, std::logic_error,
167  "Error! Thyra::NonlinearSolverBase RCP passed in through ThetaStepper::setSolver is null!"
168  );
169 
170  RCP<Teuchos::FancyOStream> out = this->getOStream();
171  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
172  Teuchos::OSTab ostab(out,1,"TS::setSolver");
173  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
174  *out << "solver = " << solver->description() << std::endl;
175  }
176 
177  solver_ = solver;
178 
179  isInitialized_ = false;
180 
181 }
182 
183 
184 template<class Scalar>
185 RCP<Thyra::NonlinearSolverBase<Scalar> >
187 {
188  return solver_;
189 }
190 
191 
192 template<class Scalar>
193 RCP<const Thyra::NonlinearSolverBase<Scalar> >
195 {
196  return solver_;
197 }
198 
199 
200 // Overridden from StepperBase
201 
202 
203 template<class Scalar>
205 {
206  return true;
207 }
208 
209 template<class Scalar>
210 RCP<StepperBase<Scalar> >
212 {
213  RCP<ThetaStepper<Scalar> >
214  stepper = Teuchos::rcp(new ThetaStepper<Scalar>);
215  stepper->isInitialized_ = isInitialized_;
216  stepper->model_ = model_; // Model is stateless so shallow copy is okay!
217 
218  if (!is_null(solver_))
219  stepper->solver_ = solver_->cloneNonlinearSolver().assert_not_null();
220 
221  stepper->basePoint_ = basePoint_;
222 
223  if (!is_null(x_))
224  stepper->x_ = x_->clone_v().assert_not_null();
225  if (!is_null(x_old_))
226  stepper->x_old_ = x_old_->clone_v().assert_not_null();
227  if (!is_null(x_pre_))
228  stepper->x_pre_ = x_pre_->clone_v().assert_not_null();
229 
230  if (!is_null(x_dot_))
231  stepper->x_dot_ = x_dot_->clone_v().assert_not_null();
232  if (!is_null(x_dot_old_))
233  stepper->x_dot_old_ = x_dot_old_->clone_v().assert_not_null();
234  if (!is_null(x_dot_really_old_))
235  stepper->x_dot_really_old_ = x_dot_really_old_->clone_v().assert_not_null();
236  if (!is_null(x_dot_base_))
237  stepper->x_dot_base_ = x_dot_base_->clone_v().assert_not_null();
238 
239  stepper->t_ = t_;
240  stepper->t_old_ = t_old_;
241 
242  stepper->dt_ = dt_;
243  stepper->dt_old_ = dt_old_;
244 
245  stepper->numSteps_ = numSteps_;
246 
247  stepper->thetaStepperType_ = thetaStepperType_;
248  stepper->theta_ = theta_;
249  stepper->predictor_corrector_begin_after_step_ = predictor_corrector_begin_after_step_;
250  stepper->default_predictor_order_ = default_predictor_order_;
251 
252  if (!is_null(neModel_))
253  stepper->neModel_
255 
256  if (!is_null(parameterList_))
257  stepper->parameterList_ = parameterList(*parameterList_);
258 
259  if (!is_null(interpolator_))
260  stepper->interpolator_
261  = interpolator_->cloneInterpolator().assert_not_null(); // ToDo: Implement cloneInterpolator()
262 
263  return stepper;
264 }
265 
266 template<class Scalar>
268  const RCP<const Thyra::ModelEvaluator<Scalar> >& model
269  )
270 {
271 
272  using Teuchos::as;
273 
274  TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
275  assertValidModel( *this, *model );
276 
277  RCP<Teuchos::FancyOStream> out = this->getOStream();
278  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
279  Teuchos::OSTab ostab(out,1,"TS::setModel");
280  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
281  *out << "model = " << model->description() << std::endl;
282  }
283  model_ = model;
284 
285  // Wipe out x. This will either be set thorugh setInitialCondition(...) or
286  // it will be taken from the model's nominal vlaues!
287  x_ = Teuchos::null;
288  x_old_ = Teuchos::null;
289  x_pre_ = Teuchos::null;
290 
291  x_dot_ = Teuchos::null;
292  x_dot_old_ = Teuchos::null;
293  x_dot_really_old_ = Teuchos::null;
294  x_dot_base_ = Teuchos::null;
295 
296  isInitialized_ = false;
297  haveInitialCondition_ = setDefaultInitialConditionFromNominalValues<Scalar>(
298  *model_, Teuchos::ptr(this) );
299 
300 }
301 
302 template<class Scalar>
304  const RCP<Thyra::ModelEvaluator<Scalar> >& model
305  )
306 {
307  this->setModel(model); // TODO 09/09/09 tscoffe: use ConstNonconstObjectContainer!
308 }
309 
310 
311 template<class Scalar>
312 RCP<const Thyra::ModelEvaluator<Scalar> >
314 {
315  return model_;
316 }
317 
318 
319 template<class Scalar>
320 RCP<Thyra::ModelEvaluator<Scalar> >
322 {
323  TEUCHOS_TEST_FOR_EXCEPT(true);
324  return Teuchos::null;
325 }
326 
327 
328 template<class Scalar>
330  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
331  )
332 {
333 
334  typedef Teuchos::ScalarTraits<Scalar> ST;
335  typedef Thyra::ModelEvaluatorBase MEB;
336 
337  basePoint_ = initialCondition;
338 
339  // x
340 
341  RCP<const Thyra::VectorBase<Scalar> >
342  x_init = initialCondition.get_x();
343 
344 #ifdef HAVE_RYTHMOS_DEBUG
345  TEUCHOS_TEST_FOR_EXCEPTION(
346  is_null(x_init), std::logic_error,
347  "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
348  "then x can not be null!" );
349 #endif
350 
351  x_ = x_init->clone_v();
352 
353  // x_dot
354 
355  RCP<const Thyra::VectorBase<Scalar> >
356  x_dot_init = initialCondition.get_x_dot();
357 
358  if (!is_null(x_dot_init)) {
359  x_dot_ = x_dot_init->clone_v();
360  }
361  else {
362  x_dot_ = createMember(x_->space());
363  assign(x_dot_.ptr(),ST::zero());
364  }
365 
366  // t
367 
368  t_ = initialCondition.get_t();
369 
370  t_old_ = t_;
371 
372  dt_old_ = 0.0;
373 
374  // x pre
375 
376  x_pre_ = x_->clone_v();
377 
378  // x old
379 
380  x_old_ = x_->clone_v();
381 
382  // x dot base
383 
384  x_dot_base_ = x_->clone_v();
385 
386  // x dot old
387 
388  x_dot_old_ = x_dot_->clone_v();
389 
390  // x dot really old
391 
392  x_dot_really_old_ = x_dot_->clone_v();
393 
394  haveInitialCondition_ = true;
395 
396 }
397 
398 
399 template<class Scalar>
400 Thyra::ModelEvaluatorBase::InArgs<Scalar>
402 {
403  return basePoint_;
404 }
405 
406 
407 template<class Scalar>
408 Scalar ThetaStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepSizeType)
409 {
410 
411  using Teuchos::as;
412  using Teuchos::incrVerbLevel;
413  typedef Teuchos::ScalarTraits<Scalar> ST;
414  typedef Thyra::NonlinearSolverBase<Scalar> NSB;
415  typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
416 
417  initialize_();
418 
419  // DEBUG
420  //this->setOverridingVerbLevel(Teuchos::VERB_EXTREME);
421  //this->setVerbLevel(Teuchos::VERB_EXTREME);
422 
423  RCP<Teuchos::FancyOStream> out = this->getOStream();
424  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
425  Teuchos::OSTab ostab(out,1,"TS::takeStep");
426  VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
427 
428  if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
429  *out
430  << "\nEntering " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
431  << "::takeStep("<<dt<<","<<toString(stepSizeType)<<") ...\n";
432  }
433 
434  dt_ = dt;
435  V_StV( x_old_.ptr(), Scalar(ST::one()), *x_ );
436  V_StV( x_dot_really_old_.ptr(), Scalar(ST::one()), *x_dot_old_ );
437  V_StV( x_dot_old_.ptr(), Scalar(ST::one()), *x_dot_ );
438 
439  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
440  *out << "\nSetting dt_ and old data ..." << std::endl;
441  *out << "\ndt_ = " << dt_;
442  *out << "\nx_old_ = " << *x_old_;
443  *out << "\nx_dot_old_ = " << *x_dot_old_;
444  *out << "\nx_dot_really_old_ = " << *x_dot_really_old_;
445  }
446 
447  if ((stepSizeType == STEP_TYPE_VARIABLE) || (dt == ST::zero())) {
448  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
449  *out << "\nThe arguments to takeStep are not valid for ThetaStepper at this time." << std::endl;
450  return(Scalar(-ST::one()));
451  }
452  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
453  *out << "\ndt = " << dt << std::endl;
454  *out << "\nstepSizeType = " << stepSizeType << std::endl;
455  }
456 
457  // compute predictor
458  obtainPredictor_();
459 
460  //
461  // Setup the nonlinear equations:
462  //
463  // substitute:
464  //
465  // x_dot = ( 1/(theta*dt) )*x + ( -1/(theta*dt) )*x_old + ( -(1-theta)/theta )*x_dot_old
466  //
467  // f( x_dot, x, t ) = 0
468  //
469 
470  const double theta = (numSteps_+1>=predictor_corrector_begin_after_step_) ? theta_ : 1.0;
471 
472  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
473  *out << "\nSetting x_dot_base_ ..." << std::endl;
474  *out << "\ntheta = " << theta;
475  *out << "\nx_ = " << *x_;
476  *out << "\nx_dot_old_ = " << *x_dot_old_;
477  }
478 
479  const Scalar coeff_x_dot = Scalar(ST::one()/(theta*dt));
480  const Scalar coeff_x = ST::one();
481 
482  const Scalar x_coeff = Scalar(-coeff_x_dot);
483  const Scalar x_dot_old_coeff = Scalar( -(ST::one()-theta)/theta);
484 
485  V_StV( x_dot_base_.ptr(), x_coeff, *x_old_ );
486  Vp_StV( x_dot_base_.ptr(), x_dot_old_coeff, *x_dot_old_);
487 
488  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
489  *out << "\nx_dot_base_ = " << *x_dot_base_;
490  }
491 
492  if(!neModel_.get()) {
493  neModel_ = Teuchos::rcp(new Rythmos::SingleResidualModelEvaluator<Scalar>());
494  }
495 
496  neModel_->initializeSingleResidualModel(
497  model_,
498  basePoint_,
499  coeff_x_dot,
500  x_dot_base_,
501  coeff_x,
502  Teuchos::null, // x_base
503  t_+dt, // t_base
504  Teuchos::null // x_bar_init
505  );
506  if( solver_->getModel().get() != neModel_.get() ) {
507  solver_->setModel(neModel_);
508  }
509 
510  solver_->setVerbLevel(this->getVerbLevel());
511 
512  //
513  // Solve the implicit nonlinear system to a tolerance of ???
514  //
515 
516  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
517  *out << "\nSolving the implicit theta-stepper timestep equation"
518  << " with theta = " << theta << "\n";
519  }
520 
521  Thyra::SolveStatus<Scalar>
522  neSolveStatus = solver_->solve(&*x_);
523 
524  // In the above solve, on input *x_ is the old value of x for the previous
525  // time step which is used as the initial guess for the solver. On output,
526  // *x_ is the converged timestep solution.
527 
528  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
529  *out << "\nOutput status of nonlinear solve:\n" << neSolveStatus;
530  }
531 
532  // 2007/05/18: rabartl: ToDo: Above, get the solve status from the above
533  // solve and at least print warning message if the solve fails! Actually,
534  // you should most likely thrown an exception if the solve fails or return
535  // false if appropriate
536 
537  //
538  // Update the step data
539  //
540 
541  // x_dot = ( 1/(theta*dt) )*x + ( -1/(theta*dt) )*x_old + ( -(1-theta)/theta )*x_dot_old
542 
543  V_StV( x_dot_.ptr(), Scalar( ST::one()/(theta*dt)), *x_ );
544  Vp_StV( x_dot_.ptr(), Scalar(-ST::one()/(theta*dt)), *x_old_ );
545  Vp_StV( x_dot_.ptr(), Scalar( -(ST::one()-theta)/theta), *x_dot_old_ );
546 
547  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
548  *out << "\nUpdating x_dot_ ...\n";
549  *out << "\nx_dot_ = " << *x_dot_ << std::endl;
550  }
551 
552  t_old_ = t_;
553  dt_old_ = dt;
554  t_ += dt;
555 
556  numSteps_++;
557 
558  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
559  *out << "\nt_old_ = " << t_old_ << std::endl;
560  *out << "\nt_ = " << t_ << std::endl;
561  }
562 
563 #ifdef HAVE_RYTHMOS_DEBUG
564 
565  if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
566  *out << "\nChecking to make sure that solution and the interpolated solution are the same! ...\n";
567 
568  {
569 
570  typedef ScalarTraits<Scalar> ST;
571  typedef typename ST::magnitudeType ScalarMag;
572  typedef ScalarTraits<ScalarMag> SMT;
573 
574  Teuchos::OSTab tab(out);
575 
576  const StepStatus<Scalar> stepStatus = this->getStepStatus();
577 
578  RCP<const Thyra::VectorBase<Scalar> >
579  x = stepStatus.solution,
580  xdot = stepStatus.solutionDot;
581 
582  Array<Scalar> time_vec = Teuchos::tuple(stepStatus.time);
583  Array<RCP<const Thyra::VectorBase<Scalar> > > x_vec, xdot_vec;
584  this->getPoints(time_vec,&x_vec,&xdot_vec,0);
585 
586  RCP<const Thyra::VectorBase<Scalar> >
587  x_interp = x_vec[0],
588  xdot_interp = xdot_vec[0];
589 
590  TEUCHOS_TEST_FOR_EXCEPT(
591  !Thyra::testRelNormDiffErr(
592  "x", *x, "x_interp", *x_interp,
593  "2*epsilon", ScalarMag(100.0*SMT::eps()),
594  "2*epsilon", ScalarMag(100.0*SMT::eps()),
595  includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
596  )
597  );
598 
599  TEUCHOS_TEST_FOR_EXCEPT(
600  !Thyra::testRelNormDiffErr(
601  "xdot", *xdot, "xdot_interp", *xdot_interp,
602  "2*epsilon", ScalarMag(100.0*SMT::eps()),
603  "2*epsilon", ScalarMag(100.0*SMT::eps()),
604  includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
605  )
606  );
607 
608  }
609 
610  // 2007/07/25: rabartl: ToDo: Move the above test into a helper function so
611  // that it can be used from lots of different places!
612 
613 #endif // HAVE_RYTHMOS_DEBUG
614 
615  if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
616  *out
617  << "\nLeaving " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
618  << "::takeStep(...) ...\n";
619  }
620 
621  return(dt);
622 
623 }
624 
625 
626 template<class Scalar>
627 const StepStatus<Scalar> ThetaStepper<Scalar>::getStepStatus() const
628 {
629 
630  typedef Teuchos::ScalarTraits<Scalar> ST;
631 
632  StepStatus<Scalar> stepStatus; // Defaults to unknown status
633 
634  if (!isInitialized_) {
635  stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
636  }
637  else if (numSteps_ > 0) {
638  stepStatus.stepStatus = STEP_STATUS_CONVERGED;
639  }
640  // else unknown
641 
642  stepStatus.stepSize = dt_;
643  stepStatus.order = 1;
644  stepStatus.time = t_;
645  stepStatus.solution = x_;
646  stepStatus.solutionDot = x_dot_;
647 
648  return(stepStatus);
649 
650 }
651 
652 
653 // Overridden from InterpolationBufferBase
654 
655 
656 template<class Scalar>
657 RCP<const Thyra::VectorSpaceBase<Scalar> >
659 {
660  return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
661 }
662 
663 
664 template<class Scalar>
666  const Array<Scalar>& time_vec,
667  const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
668  const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
669  )
670 {
671 
672  typedef Teuchos::ScalarTraits<Scalar> ST;
673  using Teuchos::as;
674 
675  initialize_();
676 
677 #ifdef HAVE_RYTHMOS_DEBUG
678  TEUCHOS_TEST_FOR_EXCEPTION(
679  time_vec.size() == 0, std::logic_error,
680  "Error, addPoints called with an empty time_vec array!\n");
681 #endif // HAVE_RYTHMOS_DEBUG
682 
683  RCP<Teuchos::FancyOStream> out = this->getOStream();
684  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
685  Teuchos::OSTab ostab(out,1,"TS::setPoints");
686 
687  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
688  *out << "time_vec = " << std::endl;
689  for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
690  *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
691  }
692  }
693  else if (time_vec.size() == 1) {
694  int n = 0;
695  t_ = time_vec[n];
696  t_old_ = t_;
697  Thyra::V_V(x_.ptr(),*x_vec[n]);
698  Thyra::V_V(x_dot_base_.ptr(),*x_);
699  }
700  else {
701  int n = time_vec.size()-1;
702  int nm1 = time_vec.size()-2;
703  t_ = time_vec[n];
704  t_old_ = time_vec[nm1];
705  Thyra::V_V(x_.ptr(),*x_vec[n]);
706  Scalar dt = t_ - t_old_;
707  Thyra::V_StV(x_dot_base_.ptr(),Scalar(-ST::one()/dt),*x_vec[nm1]);
708  }
709 }
710 
711 
712 template<class Scalar>
713 TimeRange<Scalar> ThetaStepper<Scalar>::getTimeRange() const
714 {
715  if ( !isInitialized_ && haveInitialCondition_ )
716  return timeRange<Scalar>(t_,t_);
717  if ( !isInitialized_ && !haveInitialCondition_ )
718  return invalidTimeRange<Scalar>();
719  return timeRange<Scalar>(t_old_,t_);
720 }
721 
722 
723 template<class Scalar>
725  const Array<Scalar>& time_vec,
726  Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
727  Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
728  Array<ScalarMag>* accuracy_vec
729  ) const
730 {
731  using Teuchos::constOptInArg;
732  using Teuchos::ptr;
733  typedef Teuchos::ScalarTraits<Scalar> ST;
734 
735  TEUCHOS_ASSERT(haveInitialCondition_);
736 
737  RCP<Thyra::VectorBase<Scalar> > x_temp = x_;
738  if (compareTimeValues(t_old_,t_)!=0) {
739  Scalar dt = t_ - t_old_;
740  x_temp = x_dot_base_->clone_v();
741  Thyra::Vt_S(x_temp.ptr(),Scalar(-ST::one()*dt)); // undo the scaling
742  }
743  defaultGetPoints<Scalar>(
744  t_old_, constOptInArg(*x_temp), constOptInArg(*x_dot_old_),
745  t_, constOptInArg(*x_), constOptInArg(*x_dot_),
746  time_vec, ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
747  ptr(interpolator_.get())
748  );
749 
750  /*
751  using Teuchos::as;
752  typedef Teuchos::ScalarTraits<Scalar> ST;
753  typedef typename ST::magnitudeType ScalarMag;
754  typedef Teuchos::ScalarTraits<ScalarMag> SMT;
755  typename DataStore<Scalar>::DataStoreVector_t ds_nodes;
756  typename DataStore<Scalar>::DataStoreVector_t ds_out;
757 
758 #ifdef HAVE_RYTHMOS_DEBUG
759  TEUCHOS_TEST_FOR_EXCEPT(!haveInitialCondition_);
760  TEUCHOS_TEST_FOR_EXCEPT( 0 == x_vec );
761 #endif
762 
763  RCP<Teuchos::FancyOStream> out = this->getOStream();
764  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
765  Teuchos::OSTab ostab(out,1,"TS::getPoints");
766  if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
767  *out
768  << "\nEntering " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
769  << "::getPoints(...) ...\n";
770  }
771  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
772  for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
773  *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
774  }
775  *out << "I can interpolate in the interval [" << t_old_ << "," << t_ << "]." << std::endl;
776  }
777 
778  if (t_old_ != t_) {
779  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
780  *out << "Passing two points to interpolator: " << t_old_ << " and " << t_ << std::endl;
781  }
782  DataStore<Scalar> ds_temp;
783  Scalar dt = t_ - t_old_;
784 #ifdef HAVE_RYTHMOS_DEBUG
785  TEUCHOS_TEST_FOR_EXCEPT(
786  !Thyra::testRelErr(
787  "dt", dt, "dt_", dt_,
788  "1e+4*epsilon", ScalarMag(1e+4*SMT::eps()),
789  "1e+2*epsilon", ScalarMag(1e+2*SMT::eps()),
790  as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) ? out.get() : 0
791  )
792  );
793 #endif
794  ds_temp.time = t_old_;
795  ds_temp.x = x_old_;
796  ds_temp.xdot = x_dot_old_;
797  ds_temp.accuracy = ScalarMag(dt);
798  ds_nodes.push_back(ds_temp);
799  ds_temp.time = t_;
800  ds_temp.x = x_;
801  ds_temp.xdot = x_dot_;
802  ds_temp.accuracy = ScalarMag(dt);
803  ds_nodes.push_back(ds_temp);
804  }
805  else {
806  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
807  *out << "Passing one point to interpolator: " << t_ << std::endl;
808  }
809  DataStore<Scalar> ds_temp;
810  ds_temp.time = t_;
811  ds_temp.x = x_;
812  ds_temp.xdot = x_dot_;
813  ds_temp.accuracy = ScalarMag(ST::zero());
814  ds_nodes.push_back(ds_temp);
815  }
816  interpolate<Scalar>(*interpolator_,rcp(&ds_nodes,false),time_vec,&ds_out);
817  Array<Scalar> time_out;
818  dataStoreVectorToVector(ds_out,&time_out,x_vec,xdot_vec,accuracy_vec);
819  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
820  *out << "Passing out the interpolated values:" << std::endl;
821  for (int i=0; i<Teuchos::as<int>(time_out.size()) ; ++i) {
822  if (x_vec) {
823  if ( (*x_vec)[i] == Teuchos::null) {
824  *out << "x_vec[" << i << "] = Teuchos::null" << std::endl;
825  }
826  else {
827  *out << "time[" << i << "] = " << time_out[i] << std::endl;
828  *out << "x_vec[" << i << "] = " << std::endl;
829  (*x_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
830  }
831  }
832  if (xdot_vec) {
833  if ( (*xdot_vec)[i] == Teuchos::null) {
834  *out << "xdot_vec[" << i << "] = Teuchos::null" << std::endl;
835  }
836  else {
837  *out << "xdot_vec[" << i << "] = " << std::endl;
838  (*xdot_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
839  }
840  }
841  if(accuracy_vec) {
842  *out << "accuracy[" << i << "] = " << (*accuracy_vec)[i] << std::endl;
843  }
844  }
845  }
846  if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
847  *out
848  << "Leaving " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
849  << "::getPoints(...) ...\n";
850  }
851  */
852 
853 }
854 
855 
856 template<class Scalar>
857 void ThetaStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
858 {
859  using Teuchos::as;
860 
861  TEUCHOS_ASSERT( time_vec != NULL );
862 
863  time_vec->clear();
864  if (!haveInitialCondition_) {
865  return;
866  }
867 
868  time_vec->push_back(t_old_);
869  if (numSteps_ > 0) {
870  time_vec->push_back(t_);
871  }
872 
873  RCP<Teuchos::FancyOStream> out = this->getOStream();
874  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
875  Teuchos::OSTab ostab(out,1,"TS::getNodes");
876  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
877  *out << this->description() << std::endl;
878  for (int i=0 ; i<Teuchos::as<int>(time_vec->size()) ; ++i) {
879  *out << "time_vec[" << i << "] = " << (*time_vec)[i] << std::endl;
880  }
881  }
882 }
883 
884 
885 template<class Scalar>
886 void ThetaStepper<Scalar>::removeNodes(Array<Scalar>& time_vec)
887 {
888  initialize_();
889  using Teuchos::as;
890  RCP<Teuchos::FancyOStream> out = this->getOStream();
891  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
892  Teuchos::OSTab ostab(out,1,"TS::removeNodes");
893  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
894  *out << "time_vec = " << std::endl;
895  for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
896  *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
897  }
898  }
899  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for ThetaStepper at this time.\n");
900  // TODO:
901  // if any time in time_vec matches t_ or t_old_, then do the following:
902  // remove t_old_: set t_old_ = t_ and set x_dot_base_ = x_
903  // remove t_: set t_ = t_old_ and set x_ = -dt*x_dot_base_
904 }
905 
906 
907 template<class Scalar>
909 {
910  return (thetaStepperType_==ImplicitEuler) ? 1 : 2;
911 }
912 
913 
914 // Overridden from Teuchos::ParameterListAcceptor
915 
916 
917 template <class Scalar>
919  RCP<Teuchos::ParameterList> const& paramList
920  )
921 {
922  TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
923  paramList->validateParametersAndSetDefaults(*this->getValidParameters());
924  parameterList_ = paramList;
925  Teuchos::readVerboseObjectSublist(&*parameterList_,this);
926 
927  RCP<ParameterList> pl_theta = Teuchos::sublist(parameterList_, RythmosStepControlSettings_name);
928 
929  std::string thetaStepperTypeString =
930  Teuchos::getParameter<std::string>(*pl_theta, ThetaStepperType_name);
931 
932  if (thetaStepperTypeString == "Implicit Euler")
933  thetaStepperType_ = ImplicitEuler;
934  else if (thetaStepperTypeString == "Trapezoid")
935  thetaStepperType_ = Trapezoid;
936  else
937  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
938  "Value of " << ThetaStepperType_name << " = " << thetaStepperTypeString
939  << " is invalid for Rythmos::ThetaStepper");
940 
941  default_predictor_order_ =
942  Teuchos::getParameter<int>(*pl_theta, PredictorOrder_name);
943 
944  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
945  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
946  if ( Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
947  *out << ThetaStepperType_name << " = " << thetaStepperTypeString << std::endl;
948  }
949 }
950 
951 
952 template <class Scalar>
953 RCP<Teuchos::ParameterList>
955 {
956  return(parameterList_);
957 }
958 
959 
960 template <class Scalar>
961 RCP<Teuchos::ParameterList>
963 {
964  RCP<Teuchos::ParameterList>
965  temp_param_list = parameterList_;
966  parameterList_ = Teuchos::null;
967  return(temp_param_list);
968 }
969 
970 
971 template<class Scalar>
972 RCP<const Teuchos::ParameterList>
974 {
975  using Teuchos::ParameterList;
976 
977  static RCP<const ParameterList> validPL;
978 
979  if (is_null(validPL)) {
980 
981  RCP<ParameterList> pl_top_level = Teuchos::parameterList();
982 
983  RCP<ParameterList> pl = Teuchos::sublist(pl_top_level, RythmosStepControlSettings_name);
984 
985  pl->set<std::string> ( ThetaStepperType_name, ThetaStepperType_default,
986  "Name of Stepper Type in Theta Stepper"
987  );
988 
989  pl->set<int> ( PredictorOrder_name, PredictorOrder_default,
990  "Order of Predictor in Theta Stepper, can be 0, 1, 2"
991  );
992 
993  Teuchos::setupVerboseObjectSublist(&*pl_top_level);
994  validPL = pl_top_level;
995  }
996  return validPL;
997 }
998 
999 
1000 // Overridden from Teuchos::Describable
1001 
1002 
1003 template<class Scalar>
1005  Teuchos::FancyOStream &out,
1006  const Teuchos::EVerbosityLevel verbLevel
1007  ) const
1008 {
1009  using Teuchos::as;
1010  Teuchos::OSTab tab(out);
1011  if (!isInitialized_) {
1012  out << this->description() << " : This stepper is not initialized yet" << std::endl;
1013  return;
1014  }
1015  if (
1016  as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
1017  ||
1018  as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
1019  )
1020  {
1021  out << this->description() << "::describe:" << std::endl;
1022  out << "model = " << model_->description() << std::endl;
1023  out << "solver = " << solver_->description() << std::endl;
1024  if (neModel_ == Teuchos::null) {
1025  out << "neModel = Teuchos::null" << std::endl;
1026  } else {
1027  out << "neModel = " << neModel_->description() << std::endl;
1028  }
1029  }
1030  else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
1031  out << "t_ = " << t_ << std::endl;
1032  out << "t_old_ = " << t_old_ << std::endl;
1033  }
1034  else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
1035  }
1036  else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
1037  out << "model_ = " << std::endl;
1038  model_->describe(out,verbLevel);
1039  out << "solver_ = " << std::endl;
1040  solver_->describe(out,verbLevel);
1041  if (neModel_ == Teuchos::null) {
1042  out << "neModel = Teuchos::null" << std::endl;
1043  } else {
1044  out << "neModel = " << std::endl;
1045  neModel_->describe(out,verbLevel);
1046  }
1047  out << "x_ = " << std::endl;
1048  x_->describe(out,verbLevel);
1049  out << "x_dot_base_ = " << std::endl;
1050  x_dot_base_->describe(out,verbLevel);
1051  }
1052 }
1053 
1054 
1055 // private
1056 
1057 
1058 template <class Scalar>
1059 void ThetaStepper<Scalar>::initialize_()
1060 {
1061 
1062  typedef Teuchos::ScalarTraits<Scalar> ST;
1063 
1064  if (isInitialized_)
1065  return;
1066 
1067  TEUCHOS_TEST_FOR_EXCEPT(is_null(model_));
1068  TEUCHOS_TEST_FOR_EXCEPT(is_null(solver_));
1069  TEUCHOS_TEST_FOR_EXCEPT(!haveInitialCondition_);
1070 
1071 #ifdef HAVE_RYTHMOS_DEBUG
1072  THYRA_ASSERT_VEC_SPACES(
1073  "Rythmos::ThetaStepper::setInitialCondition(...)",
1074  *x_->space(), *model_->get_x_space() );
1075 #endif // HAVE_RYTHMOS_DEBUG
1076 
1077  if ( is_null(interpolator_) ) {
1078  // If an interpolator has not been explicitly set, then just create
1079  // a default linear interpolator.
1080  interpolator_ = Teuchos::rcp(new LinearInterpolator<Scalar> );
1081  // 2007/05/18: rabartl: ToDo: Replace this with a Hermete interplator
1082  // when it is implementated!
1083  }
1084 
1085  if (thetaStepperType_ == ImplicitEuler)
1086  {
1087  theta_ = 1.0;
1088  predictor_corrector_begin_after_step_ = 2;
1089  }
1090  else
1091  {
1092  theta_ = 0.5;
1093  predictor_corrector_begin_after_step_ = 3;
1094  }
1095 
1096  isInitialized_ = true;
1097 
1098 }
1099 
1100 template<class Scalar>
1101 void ThetaStepper<Scalar>::obtainPredictor_()
1102 {
1103 
1104  using Teuchos::as;
1105  typedef Teuchos::ScalarTraits<Scalar> ST;
1106 
1107  if (!isInitialized_) {
1108  return;
1109  }
1110 
1111  RCP<Teuchos::FancyOStream> out = this->getOStream();
1112  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
1113  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1114  *out << "Obtaining predictor..." << std::endl;
1115  }
1116 
1117  const int preferred_predictor_order = std::min(default_predictor_order_, thetaStepperType_ + 1);
1118  const int max_predictor_order_at_this_timestep = std::max(0, numSteps_);
1119 
1120  const int predictor_order = std::min(preferred_predictor_order, max_predictor_order_at_this_timestep);
1121 
1122  switch (predictor_order)
1123  {
1124  case 0:
1125  V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
1126  break;
1127  case 1:
1128  {
1129  V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
1130 
1131  TEUCHOS_TEST_FOR_EXCEPT (dt_ <= 0.0);
1132 
1133  Vp_StV(x_pre_.ptr(), dt_, *x_dot_old_);
1134  }
1135  break;
1136  case 2:
1137  {
1138  V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
1139 
1140  TEUCHOS_TEST_FOR_EXCEPT (dt_ <= 0.0);
1141  TEUCHOS_TEST_FOR_EXCEPT (dt_old_ <= 0.0);
1142 
1143  const Scalar coeff_x_dot_old = (0.5 * dt_) * (2.0 + dt_/dt_old_);
1144  const Scalar coeff_x_dot_really_old = - (0.5 * dt_) * (dt_/dt_old_);
1145 
1146  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1147  *out << "x_dot_old_ = " << *x_dot_old_ << std::endl;
1148  *out << "x_dot_really_old_ = " << *x_dot_really_old_ << std::endl;
1149  }
1150 
1151  Vp_StV( x_pre_.ptr(), coeff_x_dot_old, *x_dot_old_);
1152  Vp_StV( x_pre_.ptr(), coeff_x_dot_really_old, *x_dot_really_old_);
1153  }
1154  break;
1155  default:
1156  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
1157  "Invalid predictor order " << predictor_order << ". Valid values are 0, 1, and 2.");
1158  }
1159 
1160  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1161  *out << "x_pre_ = " << *x_pre_ << std::endl;
1162  }
1163 
1164  // copy to current solution
1165  V_StV(x_.ptr(), Scalar(ST::one()), *x_pre_);
1166 }
1167 
1168 //
1169 // Explicit Instantiation macro
1170 //
1171 // Must be expanded from within the Rythmos namespace!
1172 //
1173 
1174 #define RYTHMOS_THETA_STEPPER_INSTANT(SCALAR) \
1175  \
1176  template class ThetaStepper< SCALAR >; \
1177  \
1178  template RCP< ThetaStepper< SCALAR > > \
1179  thetaStepper( \
1180  const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \
1181  const RCP<Thyra::NonlinearSolverBase< SCALAR > >& solver, \
1182  RCP<Teuchos::ParameterList>& parameterList \
1183  );
1184 
1185 
1186 } // namespace Rythmos
1187 
1188 #endif // HAVE_RYTHMOS_EXPERIMENTAL
1189 
1190 #endif //Rythmos_THETA_STEPPER_DEF_H
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
void getNodes(Array< Scalar > *time_vec) const
RCP< InterpolatorBase< Scalar > > getNonconstInterpolator()
void removeNodes(Array< Scalar > &time_vec)
RCP< Teuchos::ParameterList > unsetParameterList()
void setInterpolator(const RCP< InterpolatorBase< Scalar > > &interpolator)
Redefined from InterpolatorAcceptingObjectBase.
Thyra::ModelEvaluatorBase::InArgs< Scalar > getInitialCondition() const
RCP< Teuchos::ParameterList > getNonconstParameterList()
RCP< const Teuchos::ParameterList > getValidParameters() const
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)
void setModel(const RCP< const Thyra::ModelEvaluator< Scalar > > &model)
RCP< const InterpolatorBase< Scalar > > getInterpolator() const
RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
void setParameterList(RCP< Teuchos::ParameterList > const &paramList)
void setInitialCondition(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &initialCondition)
void setSolver(const RCP< Thyra::NonlinearSolverBase< Scalar > > &solver)
Decorator subclass for a steady-state version of a DAE for single-residual time stepper methods...
TimeRange< Scalar > getTimeRange() const
const StepStatus< Scalar > getStepStatus() const
RCP< Thyra::ModelEvaluator< Scalar > > getNonconstModel()
RCP< InterpolatorBase< Scalar > > unSetInterpolator()
bool supportsCloning() const
Returns true.
RCP< Thyra::NonlinearSolverBase< Scalar > > getNonconstSolver()
RCP< const Thyra::NonlinearSolverBase< Scalar > > getSolver() const
RCP< StepperBase< Scalar > > cloneStepperAlgorithm() const
Creates copies of all internal data (including the parameter list) except the model which is assumed ...
RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Scalar takeStep(Scalar dt, StepSizeType flag)
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
void setNonconstModel(const RCP< Thyra::ModelEvaluator< Scalar > > &model)
bool isImplicit() const
Return if this stepper is an implicit stepper.