Rythmos - Transient Integration for Differential Equations  Version of the Day
 All Classes Functions Variables Typedefs Pages
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 
336  basePoint_ = initialCondition;
337 
338  // x
339 
340  RCP<const Thyra::VectorBase<Scalar> >
341  x_init = initialCondition.get_x();
342 
343 #ifdef HAVE_RYTHMOS_DEBUG
344  TEUCHOS_TEST_FOR_EXCEPTION(
345  is_null(x_init), std::logic_error,
346  "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
347  "then x can not be null!" );
348 #endif
349 
350  x_ = x_init->clone_v();
351 
352  // x_dot
353 
354  RCP<const Thyra::VectorBase<Scalar> >
355  x_dot_init = initialCondition.get_x_dot();
356 
357  if (!is_null(x_dot_init)) {
358  x_dot_ = x_dot_init->clone_v();
359  }
360  else {
361  x_dot_ = createMember(x_->space());
362  assign(x_dot_.ptr(),ST::zero());
363  }
364 
365  // t
366 
367  t_ = initialCondition.get_t();
368 
369  t_old_ = t_;
370 
371  dt_old_ = 0.0;
372 
373  // x pre
374 
375  x_pre_ = x_->clone_v();
376 
377  // x old
378 
379  x_old_ = x_->clone_v();
380 
381  // x dot base
382 
383  x_dot_base_ = x_->clone_v();
384 
385  // x dot old
386 
387  x_dot_old_ = x_dot_->clone_v();
388 
389  // x dot really old
390 
391  x_dot_really_old_ = x_dot_->clone_v();
392 
393  haveInitialCondition_ = true;
394 
395 }
396 
397 
398 template<class Scalar>
399 Thyra::ModelEvaluatorBase::InArgs<Scalar>
401 {
402  return basePoint_;
403 }
404 
405 
406 template<class Scalar>
407 Scalar ThetaStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepSizeType)
408 {
409 
410  using Teuchos::as;
411  using Teuchos::incrVerbLevel;
412  typedef Teuchos::ScalarTraits<Scalar> ST;
413  typedef Thyra::NonlinearSolverBase<Scalar> NSB;
414  typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
415 
416  initialize_();
417 
418  // DEBUG
419  //this->setOverridingVerbLevel(Teuchos::VERB_EXTREME);
420  //this->setVerbLevel(Teuchos::VERB_EXTREME);
421 
422  RCP<Teuchos::FancyOStream> out = this->getOStream();
423  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
424  Teuchos::OSTab ostab(out,1,"TS::takeStep");
425  VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
426 
427  if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
428  *out
429  << "\nEntering " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
430  << "::takeStep("<<dt<<","<<toString(stepSizeType)<<") ...\n";
431  }
432 
433  dt_ = dt;
434  V_StV( x_old_.ptr(), Scalar(ST::one()), *x_ );
435  V_StV( x_dot_really_old_.ptr(), Scalar(ST::one()), *x_dot_old_ );
436  V_StV( x_dot_old_.ptr(), Scalar(ST::one()), *x_dot_ );
437 
438  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
439  *out << "\nSetting dt_ and old data ..." << std::endl;
440  *out << "\ndt_ = " << dt_;
441  *out << "\nx_old_ = " << *x_old_;
442  *out << "\nx_dot_old_ = " << *x_dot_old_;
443  *out << "\nx_dot_really_old_ = " << *x_dot_really_old_;
444  }
445 
446  if ((stepSizeType == STEP_TYPE_VARIABLE) || (dt == ST::zero())) {
447  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
448  *out << "\nThe arguments to takeStep are not valid for ThetaStepper at this time." << std::endl;
449  return(Scalar(-ST::one()));
450  }
451  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
452  *out << "\ndt = " << dt << std::endl;
453  *out << "\nstepSizeType = " << stepSizeType << std::endl;
454  }
455 
456  // compute predictor
457  obtainPredictor_();
458 
459  //
460  // Setup the nonlinear equations:
461  //
462  // substitute:
463  //
464  // x_dot = ( 1/(theta*dt) )*x + ( -1/(theta*dt) )*x_old + ( -(1-theta)/theta )*x_dot_old
465  //
466  // f( x_dot, x, t ) = 0
467  //
468 
469  const double theta = (numSteps_+1>=predictor_corrector_begin_after_step_) ? theta_ : 1.0;
470 
471  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
472  *out << "\nSetting x_dot_base_ ..." << std::endl;
473  *out << "\ntheta = " << theta;
474  *out << "\nx_ = " << *x_;
475  *out << "\nx_dot_old_ = " << *x_dot_old_;
476  }
477 
478  const Scalar coeff_x_dot = Scalar(ST::one()/(theta*dt));
479  const Scalar coeff_x = ST::one();
480 
481  const Scalar x_coeff = Scalar(-coeff_x_dot);
482  const Scalar x_dot_old_coeff = Scalar( -(ST::one()-theta)/theta);
483 
484  V_StV( x_dot_base_.ptr(), x_coeff, *x_old_ );
485  Vp_StV( x_dot_base_.ptr(), x_dot_old_coeff, *x_dot_old_);
486 
487  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
488  *out << "\nx_dot_base_ = " << *x_dot_base_;
489  }
490 
491  if(!neModel_.get()) {
492  neModel_ = Teuchos::rcp(new Rythmos::SingleResidualModelEvaluator<Scalar>());
493  }
494 
495  neModel_->initializeSingleResidualModel(
496  model_,
497  basePoint_,
498  coeff_x_dot,
499  x_dot_base_,
500  coeff_x,
501  Teuchos::null, // x_base
502  t_+dt, // t_base
503  Teuchos::null // x_bar_init
504  );
505  if( solver_->getModel().get() != neModel_.get() ) {
506  solver_->setModel(neModel_);
507  }
508 
509  solver_->setVerbLevel(this->getVerbLevel());
510 
511  //
512  // Solve the implicit nonlinear system to a tolerance of ???
513  //
514 
515  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
516  *out << "\nSolving the implicit theta-stepper timestep equation"
517  << " with theta = " << theta << "\n";
518  }
519 
520  Thyra::SolveStatus<Scalar>
521  neSolveStatus = solver_->solve(&*x_);
522 
523  // In the above solve, on input *x_ is the old value of x for the previous
524  // time step which is used as the initial guess for the solver. On output,
525  // *x_ is the converged timestep solution.
526 
527  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
528  *out << "\nOutput status of nonlinear solve:\n" << neSolveStatus;
529  }
530 
531  // 2007/05/18: rabartl: ToDo: Above, get the solve status from the above
532  // solve and at least print warning message if the solve fails! Actually,
533  // you should most likely thrown an exception if the solve fails or return
534  // false if appropriate
535 
536  //
537  // Update the step data
538  //
539 
540  // x_dot = ( 1/(theta*dt) )*x + ( -1/(theta*dt) )*x_old + ( -(1-theta)/theta )*x_dot_old
541 
542  V_StV( x_dot_.ptr(), Scalar( ST::one()/(theta*dt)), *x_ );
543  Vp_StV( x_dot_.ptr(), Scalar(-ST::one()/(theta*dt)), *x_old_ );
544  Vp_StV( x_dot_.ptr(), Scalar( -(ST::one()-theta)/theta), *x_dot_old_ );
545 
546  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
547  *out << "\nUpdating x_dot_ ...\n";
548  *out << "\nx_dot_ = " << *x_dot_ << std::endl;
549  }
550 
551  t_old_ = t_;
552  dt_old_ = dt;
553  t_ += dt;
554 
555  numSteps_++;
556 
557  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
558  *out << "\nt_old_ = " << t_old_ << std::endl;
559  *out << "\nt_ = " << t_ << std::endl;
560  }
561 
562 #ifdef HAVE_RYTHMOS_DEBUG
563 
564  if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
565  *out << "\nChecking to make sure that solution and the interpolated solution are the same! ...\n";
566 
567  {
568 
569  using SMT = ScalarTraits<ScalarMag>;
570 
571  Teuchos::OSTab tab(out);
572 
573  const StepStatus<Scalar> stepStatus = this->getStepStatus();
574 
575  RCP<const Thyra::VectorBase<Scalar> >
576  x = stepStatus.solution,
577  xdot = stepStatus.solutionDot;
578 
579  Array<Scalar> time_vec = Teuchos::tuple(stepStatus.time);
580  Array<RCP<const Thyra::VectorBase<Scalar> > > x_vec, xdot_vec;
581  this->getPoints(time_vec,&x_vec,&xdot_vec,0);
582 
583  RCP<const Thyra::VectorBase<Scalar> >
584  x_interp = x_vec[0],
585  xdot_interp = xdot_vec[0];
586 
587  TEUCHOS_TEST_FOR_EXCEPT(
588  !Thyra::testRelNormDiffErr(
589  "x", *x, "x_interp", *x_interp,
590  "2*epsilon", ScalarMag(100.0*SMT::eps()),
591  "2*epsilon", ScalarMag(100.0*SMT::eps()),
592  includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
593  )
594  );
595 
596  TEUCHOS_TEST_FOR_EXCEPT(
597  !Thyra::testRelNormDiffErr(
598  "xdot", *xdot, "xdot_interp", *xdot_interp,
599  "2*epsilon", ScalarMag(100.0*SMT::eps()),
600  "2*epsilon", ScalarMag(100.0*SMT::eps()),
601  includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
602  )
603  );
604 
605  }
606 
607  // 2007/07/25: rabartl: ToDo: Move the above test into a helper function so
608  // that it can be used from lots of different places!
609 
610 #endif // HAVE_RYTHMOS_DEBUG
611 
612  if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
613  *out
614  << "\nLeaving " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
615  << "::takeStep(...) ...\n";
616  }
617 
618  return(dt);
619 
620 }
621 
622 
623 template<class Scalar>
624 const StepStatus<Scalar> ThetaStepper<Scalar>::getStepStatus() const
625 {
626 
627  StepStatus<Scalar> stepStatus; // Defaults to unknown status
628 
629  if (!isInitialized_) {
630  stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
631  }
632  else if (numSteps_ > 0) {
633  stepStatus.stepStatus = STEP_STATUS_CONVERGED;
634  }
635  // else unknown
636 
637  stepStatus.stepSize = dt_;
638  stepStatus.order = 1;
639  stepStatus.time = t_;
640  stepStatus.solution = x_;
641  stepStatus.solutionDot = x_dot_;
642 
643  return(stepStatus);
644 
645 }
646 
647 
648 // Overridden from InterpolationBufferBase
649 
650 
651 template<class Scalar>
652 RCP<const Thyra::VectorSpaceBase<Scalar> >
654 {
655  return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
656 }
657 
658 
659 template<class Scalar>
661  const Array<Scalar>& time_vec,
662  const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
663  const Array<RCP<const Thyra::VectorBase<Scalar> > >& /* xdot_vec */
664  )
665 {
666 
667  typedef Teuchos::ScalarTraits<Scalar> ST;
668  using Teuchos::as;
669 
670  initialize_();
671 
672 #ifdef HAVE_RYTHMOS_DEBUG
673  TEUCHOS_TEST_FOR_EXCEPTION(
674  time_vec.size() == 0, std::logic_error,
675  "Error, addPoints called with an empty time_vec array!\n");
676 #endif // HAVE_RYTHMOS_DEBUG
677 
678  RCP<Teuchos::FancyOStream> out = this->getOStream();
679  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
680  Teuchos::OSTab ostab(out,1,"TS::setPoints");
681 
682  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
683  *out << "time_vec = " << std::endl;
684  for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
685  *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
686  }
687  }
688  else if (time_vec.size() == 1) {
689  int n = 0;
690  t_ = time_vec[n];
691  t_old_ = t_;
692  Thyra::V_V(x_.ptr(),*x_vec[n]);
693  Thyra::V_V(x_dot_base_.ptr(),*x_);
694  }
695  else {
696  int n = time_vec.size()-1;
697  int nm1 = time_vec.size()-2;
698  t_ = time_vec[n];
699  t_old_ = time_vec[nm1];
700  Thyra::V_V(x_.ptr(),*x_vec[n]);
701  Scalar dt = t_ - t_old_;
702  Thyra::V_StV(x_dot_base_.ptr(),Scalar(-ST::one()/dt),*x_vec[nm1]);
703  }
704 }
705 
706 
707 template<class Scalar>
708 TimeRange<Scalar> ThetaStepper<Scalar>::getTimeRange() const
709 {
710  if ( !isInitialized_ && haveInitialCondition_ )
711  return timeRange<Scalar>(t_,t_);
712  if ( !isInitialized_ && !haveInitialCondition_ )
713  return invalidTimeRange<Scalar>();
714  return timeRange<Scalar>(t_old_,t_);
715 }
716 
717 
718 template<class Scalar>
720  const Array<Scalar>& time_vec,
721  Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
722  Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
723  Array<ScalarMag>* accuracy_vec
724  ) const
725 {
726  using Teuchos::constOptInArg;
727  using Teuchos::ptr;
728  typedef Teuchos::ScalarTraits<Scalar> ST;
729 
730  TEUCHOS_ASSERT(haveInitialCondition_);
731 
732  RCP<Thyra::VectorBase<Scalar> > x_temp = x_;
733  if (compareTimeValues(t_old_,t_)!=0) {
734  Scalar dt = t_ - t_old_;
735  x_temp = x_dot_base_->clone_v();
736  Thyra::Vt_S(x_temp.ptr(),Scalar(-ST::one()*dt)); // undo the scaling
737  }
738  defaultGetPoints<Scalar>(
739  t_old_, constOptInArg(*x_temp), constOptInArg(*x_dot_old_),
740  t_, constOptInArg(*x_), constOptInArg(*x_dot_),
741  time_vec, ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
742  ptr(interpolator_.get())
743  );
744 
745  /*
746  using Teuchos::as;
747  typedef Teuchos::ScalarTraits<Scalar> ST;
748  typedef typename ST::magnitudeType ScalarMag;
749  typedef Teuchos::ScalarTraits<ScalarMag> SMT;
750  typename DataStore<Scalar>::DataStoreVector_t ds_nodes;
751  typename DataStore<Scalar>::DataStoreVector_t ds_out;
752 
753 #ifdef HAVE_RYTHMOS_DEBUG
754  TEUCHOS_TEST_FOR_EXCEPT(!haveInitialCondition_);
755  TEUCHOS_TEST_FOR_EXCEPT( 0 == x_vec );
756 #endif
757 
758  RCP<Teuchos::FancyOStream> out = this->getOStream();
759  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
760  Teuchos::OSTab ostab(out,1,"TS::getPoints");
761  if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
762  *out
763  << "\nEntering " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
764  << "::getPoints(...) ...\n";
765  }
766  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
767  for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
768  *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
769  }
770  *out << "I can interpolate in the interval [" << t_old_ << "," << t_ << "]." << std::endl;
771  }
772 
773  if (t_old_ != t_) {
774  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
775  *out << "Passing two points to interpolator: " << t_old_ << " and " << t_ << std::endl;
776  }
777  DataStore<Scalar> ds_temp;
778  Scalar dt = t_ - t_old_;
779 #ifdef HAVE_RYTHMOS_DEBUG
780  TEUCHOS_TEST_FOR_EXCEPT(
781  !Thyra::testRelErr(
782  "dt", dt, "dt_", dt_,
783  "1e+4*epsilon", ScalarMag(1e+4*SMT::eps()),
784  "1e+2*epsilon", ScalarMag(1e+2*SMT::eps()),
785  as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) ? out.get() : 0
786  )
787  );
788 #endif
789  ds_temp.time = t_old_;
790  ds_temp.x = x_old_;
791  ds_temp.xdot = x_dot_old_;
792  ds_temp.accuracy = ScalarMag(dt);
793  ds_nodes.push_back(ds_temp);
794  ds_temp.time = t_;
795  ds_temp.x = x_;
796  ds_temp.xdot = x_dot_;
797  ds_temp.accuracy = ScalarMag(dt);
798  ds_nodes.push_back(ds_temp);
799  }
800  else {
801  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
802  *out << "Passing one point to interpolator: " << t_ << std::endl;
803  }
804  DataStore<Scalar> ds_temp;
805  ds_temp.time = t_;
806  ds_temp.x = x_;
807  ds_temp.xdot = x_dot_;
808  ds_temp.accuracy = ScalarMag(ST::zero());
809  ds_nodes.push_back(ds_temp);
810  }
811  interpolate<Scalar>(*interpolator_,rcp(&ds_nodes,false),time_vec,&ds_out);
812  Array<Scalar> time_out;
813  dataStoreVectorToVector(ds_out,&time_out,x_vec,xdot_vec,accuracy_vec);
814  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
815  *out << "Passing out the interpolated values:" << std::endl;
816  for (int i=0; i<Teuchos::as<int>(time_out.size()) ; ++i) {
817  if (x_vec) {
818  if ( (*x_vec)[i] == Teuchos::null) {
819  *out << "x_vec[" << i << "] = Teuchos::null" << std::endl;
820  }
821  else {
822  *out << "time[" << i << "] = " << time_out[i] << std::endl;
823  *out << "x_vec[" << i << "] = " << std::endl;
824  (*x_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
825  }
826  }
827  if (xdot_vec) {
828  if ( (*xdot_vec)[i] == Teuchos::null) {
829  *out << "xdot_vec[" << i << "] = Teuchos::null" << std::endl;
830  }
831  else {
832  *out << "xdot_vec[" << i << "] = " << std::endl;
833  (*xdot_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
834  }
835  }
836  if(accuracy_vec) {
837  *out << "accuracy[" << i << "] = " << (*accuracy_vec)[i] << std::endl;
838  }
839  }
840  }
841  if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
842  *out
843  << "Leaving " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
844  << "::getPoints(...) ...\n";
845  }
846  */
847 
848 }
849 
850 
851 template<class Scalar>
852 void ThetaStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
853 {
854  using Teuchos::as;
855 
856  TEUCHOS_ASSERT( time_vec != NULL );
857 
858  time_vec->clear();
859  if (!haveInitialCondition_) {
860  return;
861  }
862 
863  time_vec->push_back(t_old_);
864  if (numSteps_ > 0) {
865  time_vec->push_back(t_);
866  }
867 
868  RCP<Teuchos::FancyOStream> out = this->getOStream();
869  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
870  Teuchos::OSTab ostab(out,1,"TS::getNodes");
871  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
872  *out << this->description() << std::endl;
873  for (int i=0 ; i<Teuchos::as<int>(time_vec->size()) ; ++i) {
874  *out << "time_vec[" << i << "] = " << (*time_vec)[i] << std::endl;
875  }
876  }
877 }
878 
879 
880 template<class Scalar>
881 void ThetaStepper<Scalar>::removeNodes(Array<Scalar>& time_vec)
882 {
883  initialize_();
884  using Teuchos::as;
885  RCP<Teuchos::FancyOStream> out = this->getOStream();
886  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
887  Teuchos::OSTab ostab(out,1,"TS::removeNodes");
888  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
889  *out << "time_vec = " << std::endl;
890  for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
891  *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
892  }
893  }
894  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for ThetaStepper at this time.\n");
895  // TODO:
896  // if any time in time_vec matches t_ or t_old_, then do the following:
897  // remove t_old_: set t_old_ = t_ and set x_dot_base_ = x_
898  // remove t_: set t_ = t_old_ and set x_ = -dt*x_dot_base_
899 }
900 
901 
902 template<class Scalar>
904 {
905  return (thetaStepperType_==ImplicitEuler) ? 1 : 2;
906 }
907 
908 
909 // Overridden from Teuchos::ParameterListAcceptor
910 
911 
912 template <class Scalar>
914  RCP<Teuchos::ParameterList> const& paramList
915  )
916 {
917  TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
918  paramList->validateParametersAndSetDefaults(*this->getValidParameters());
919  parameterList_ = paramList;
920  Teuchos::readVerboseObjectSublist(&*parameterList_,this);
921 
922  RCP<ParameterList> pl_theta = Teuchos::sublist(parameterList_, RythmosStepControlSettings_name);
923 
924  std::string thetaStepperTypeString =
925  Teuchos::getParameter<std::string>(*pl_theta, ThetaStepperType_name);
926 
927  if (thetaStepperTypeString == "Implicit Euler")
928  thetaStepperType_ = ImplicitEuler;
929  else if (thetaStepperTypeString == "Trapezoid")
930  thetaStepperType_ = Trapezoid;
931  else
932  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
933  "Value of " << ThetaStepperType_name << " = " << thetaStepperTypeString
934  << " is invalid for Rythmos::ThetaStepper");
935 
936  default_predictor_order_ =
937  Teuchos::getParameter<int>(*pl_theta, PredictorOrder_name);
938 
939  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
940  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
941  if ( Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
942  *out << ThetaStepperType_name << " = " << thetaStepperTypeString << std::endl;
943  }
944 }
945 
946 
947 template <class Scalar>
948 RCP<Teuchos::ParameterList>
950 {
951  return(parameterList_);
952 }
953 
954 
955 template <class Scalar>
956 RCP<Teuchos::ParameterList>
958 {
959  RCP<Teuchos::ParameterList>
960  temp_param_list = parameterList_;
961  parameterList_ = Teuchos::null;
962  return(temp_param_list);
963 }
964 
965 
966 template<class Scalar>
967 RCP<const Teuchos::ParameterList>
969 {
970  using Teuchos::ParameterList;
971 
972  static RCP<const ParameterList> validPL;
973 
974  if (is_null(validPL)) {
975 
976  RCP<ParameterList> pl_top_level = Teuchos::parameterList();
977 
978  RCP<ParameterList> pl = Teuchos::sublist(pl_top_level, RythmosStepControlSettings_name);
979 
980  pl->set<std::string> ( ThetaStepperType_name, ThetaStepperType_default,
981  "Name of Stepper Type in Theta Stepper"
982  );
983 
984  pl->set<int> ( PredictorOrder_name, PredictorOrder_default,
985  "Order of Predictor in Theta Stepper, can be 0, 1, 2"
986  );
987 
988  Teuchos::setupVerboseObjectSublist(&*pl_top_level);
989  validPL = pl_top_level;
990  }
991  return validPL;
992 }
993 
994 
995 // Overridden from Teuchos::Describable
996 
997 
998 template<class Scalar>
1000  Teuchos::FancyOStream &out,
1001  const Teuchos::EVerbosityLevel verbLevel
1002  ) const
1003 {
1004  using Teuchos::as;
1005  Teuchos::OSTab tab(out);
1006  if (!isInitialized_) {
1007  out << this->description() << " : This stepper is not initialized yet" << std::endl;
1008  return;
1009  }
1010  if (
1011  as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
1012  ||
1013  as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
1014  )
1015  {
1016  out << this->description() << "::describe:" << std::endl;
1017  out << "model = " << model_->description() << std::endl;
1018  out << "solver = " << solver_->description() << std::endl;
1019  if (neModel_ == Teuchos::null) {
1020  out << "neModel = Teuchos::null" << std::endl;
1021  } else {
1022  out << "neModel = " << neModel_->description() << std::endl;
1023  }
1024  }
1025  else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
1026  out << "t_ = " << t_ << std::endl;
1027  out << "t_old_ = " << t_old_ << std::endl;
1028  }
1029  else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
1030  }
1031  else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
1032  out << "model_ = " << std::endl;
1033  model_->describe(out,verbLevel);
1034  out << "solver_ = " << std::endl;
1035  solver_->describe(out,verbLevel);
1036  if (neModel_ == Teuchos::null) {
1037  out << "neModel = Teuchos::null" << std::endl;
1038  } else {
1039  out << "neModel = " << std::endl;
1040  neModel_->describe(out,verbLevel);
1041  }
1042  out << "x_ = " << std::endl;
1043  x_->describe(out,verbLevel);
1044  out << "x_dot_base_ = " << std::endl;
1045  x_dot_base_->describe(out,verbLevel);
1046  }
1047 }
1048 
1049 
1050 // private
1051 
1052 
1053 template <class Scalar>
1054 void ThetaStepper<Scalar>::initialize_()
1055 {
1056 
1057  if (isInitialized_)
1058  return;
1059 
1060  TEUCHOS_TEST_FOR_EXCEPT(is_null(model_));
1061  TEUCHOS_TEST_FOR_EXCEPT(is_null(solver_));
1062  TEUCHOS_TEST_FOR_EXCEPT(!haveInitialCondition_);
1063 
1064 #ifdef HAVE_RYTHMOS_DEBUG
1065  THYRA_ASSERT_VEC_SPACES(
1066  "Rythmos::ThetaStepper::setInitialCondition(...)",
1067  *x_->space(), *model_->get_x_space() );
1068 #endif // HAVE_RYTHMOS_DEBUG
1069 
1070  if ( is_null(interpolator_) ) {
1071  // If an interpolator has not been explicitly set, then just create
1072  // a default linear interpolator.
1073  interpolator_ = Teuchos::rcp(new LinearInterpolator<Scalar> );
1074  // 2007/05/18: rabartl: ToDo: Replace this with a Hermete interplator
1075  // when it is implementated!
1076  }
1077 
1078  if (thetaStepperType_ == ImplicitEuler)
1079  {
1080  theta_ = 1.0;
1081  predictor_corrector_begin_after_step_ = 2;
1082  }
1083  else
1084  {
1085  theta_ = 0.5;
1086  predictor_corrector_begin_after_step_ = 3;
1087  }
1088 
1089  isInitialized_ = true;
1090 
1091 }
1092 
1093 template<class Scalar>
1094 void ThetaStepper<Scalar>::obtainPredictor_()
1095 {
1096 
1097  using Teuchos::as;
1098  typedef Teuchos::ScalarTraits<Scalar> ST;
1099 
1100  if (!isInitialized_) {
1101  return;
1102  }
1103 
1104  RCP<Teuchos::FancyOStream> out = this->getOStream();
1105  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
1106  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1107  *out << "Obtaining predictor..." << std::endl;
1108  }
1109 
1110  const int preferred_predictor_order = std::min(default_predictor_order_, thetaStepperType_ + 1);
1111  const int max_predictor_order_at_this_timestep = std::max(0, numSteps_);
1112 
1113  const int predictor_order = std::min(preferred_predictor_order, max_predictor_order_at_this_timestep);
1114 
1115  switch (predictor_order)
1116  {
1117  case 0:
1118  V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
1119  break;
1120  case 1:
1121  {
1122  V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
1123 
1124  TEUCHOS_TEST_FOR_EXCEPT (dt_ <= 0.0);
1125 
1126  Vp_StV(x_pre_.ptr(), dt_, *x_dot_old_);
1127  }
1128  break;
1129  case 2:
1130  {
1131  V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
1132 
1133  TEUCHOS_TEST_FOR_EXCEPT (dt_ <= 0.0);
1134  TEUCHOS_TEST_FOR_EXCEPT (dt_old_ <= 0.0);
1135 
1136  const Scalar coeff_x_dot_old = (0.5 * dt_) * (2.0 + dt_/dt_old_);
1137  const Scalar coeff_x_dot_really_old = - (0.5 * dt_) * (dt_/dt_old_);
1138 
1139  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1140  *out << "x_dot_old_ = " << *x_dot_old_ << std::endl;
1141  *out << "x_dot_really_old_ = " << *x_dot_really_old_ << std::endl;
1142  }
1143 
1144  Vp_StV( x_pre_.ptr(), coeff_x_dot_old, *x_dot_old_);
1145  Vp_StV( x_pre_.ptr(), coeff_x_dot_really_old, *x_dot_really_old_);
1146  }
1147  break;
1148  default:
1149  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
1150  "Invalid predictor order " << predictor_order << ". Valid values are 0, 1, and 2.");
1151  }
1152 
1153  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1154  *out << "x_pre_ = " << *x_pre_ << std::endl;
1155  }
1156 
1157  // copy to current solution
1158  V_StV(x_.ptr(), Scalar(ST::one()), *x_pre_);
1159 }
1160 
1161 //
1162 // Explicit Instantiation macro
1163 //
1164 // Must be expanded from within the Rythmos namespace!
1165 //
1166 
1167 #define RYTHMOS_THETA_STEPPER_INSTANT(SCALAR) \
1168  \
1169  template class ThetaStepper< SCALAR >; \
1170  \
1171  template RCP< ThetaStepper< SCALAR > > \
1172  thetaStepper( \
1173  const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \
1174  const RCP<Thyra::NonlinearSolverBase< SCALAR > >& solver, \
1175  RCP<Teuchos::ParameterList>& parameterList \
1176  );
1177 
1178 
1179 } // namespace Rythmos
1180 
1181 #endif // HAVE_RYTHMOS_EXPERIMENTAL
1182 
1183 #endif //Rythmos_THETA_STEPPER_DEF_H
bool supportsCloning() const
Returns true.
RCP< const Teuchos::ParameterList > getValidParameters() const
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) 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.
RCP< Teuchos::ParameterList > getNonconstParameterList()
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::NonlinearSolverBase< Scalar > > getSolver() const
RCP< const InterpolatorBase< Scalar > > getInterpolator() 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 Thyra::VectorSpaceBase< Scalar > > get_x_space() const
void setParameterList(RCP< Teuchos::ParameterList > const &paramList)
void setInitialCondition(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &initialCondition)
void setSolver(const RCP< Thyra::NonlinearSolverBase< Scalar > > &solver)
void getNodes(Array< Scalar > *time_vec) const
Decorator subclass for a steady-state version of a DAE for single-residual time stepper methods...
bool isImplicit() const
Return if this stepper is an implicit stepper.
int getOrder() const
RCP< Thyra::ModelEvaluator< Scalar > > getNonconstModel()
RCP< InterpolatorBase< Scalar > > unSetInterpolator()
TimeRange< Scalar > getTimeRange() const
RCP< Thyra::NonlinearSolverBase< Scalar > > getNonconstSolver()
Thyra::ModelEvaluatorBase::InArgs< Scalar > getInitialCondition() const
const StepStatus< Scalar > getStepStatus() const
Scalar takeStep(Scalar dt, StepSizeType flag)
RCP< StepperBase< Scalar > > cloneStepperAlgorithm() const
Creates copies of all internal data (including the parameter list) except the model which is assumed ...
RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
void setNonconstModel(const RCP< Thyra::ModelEvaluator< Scalar > > &model)