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