Rythmos - Transient Integration for Differential Equations  Version of the Day
Rythmos_ImplicitRKStepper_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_IMPLICIT_RK_STEPPER_DEF_H
30 #define Rythmos_IMPLICIT_RK_STEPPER_DEF_H
31 
32 #include "Rythmos_ImplicitRKStepper_decl.hpp"
33 
34 #include "Rythmos_StepperHelpers.hpp"
35 #include "Rythmos_SingleResidualModelEvaluator.hpp"
36 #include "Rythmos_ImplicitRKModelEvaluator.hpp"
37 #include "Rythmos_DiagonalImplicitRKModelEvaluator.hpp"
38 #include "Rythmos_RKButcherTableau.hpp"
39 #include "Rythmos_RKButcherTableauHelpers.hpp"
40 
41 #include "Thyra_ModelEvaluatorHelpers.hpp"
42 #include "Thyra_ProductVectorSpaceBase.hpp"
43 #include "Thyra_AssertOp.hpp"
44 #include "Thyra_TestingTools.hpp"
45 #include "Rythmos_ImplicitBDFStepperRampingStepControl.hpp"
46 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
47 #include "Teuchos_as.hpp"
48 
49 
50 namespace Rythmos {
51 
56 template<class Scalar>
57 RCP<ImplicitRKStepper<Scalar> >
59 {
60  RCP<ImplicitRKStepper<Scalar> > stepper(new ImplicitRKStepper<Scalar>());
61  return stepper;
62 }
63 
64 template<class Scalar>
65 RCP<ImplicitRKStepper<Scalar> >
66 implicitRKStepper(
67  const RCP<const Thyra::ModelEvaluator<Scalar> >& model,
68  const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
69  const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >& irk_W_factory,
70  const RCP<const RKButcherTableauBase<Scalar> >& irkbt
71  )
72 {
73  RCP<ImplicitRKStepper<Scalar> > stepper(new ImplicitRKStepper<Scalar>());
74 
75  stepper->setModel(model);
76  stepper->setSolver(solver);
77  stepper->set_W_factory(irk_W_factory);
78  stepper->setRKButcherTableau(irkbt);
79  return stepper;
80 }
81 
82 
83 // ////////////////////////////
84 // Defintions
85 
86 
87 // Constructors, intializers, Misc.
88 
89 
90 template<class Scalar>
92 {
93  this->defaultInitializeAll_();
94  irkButcherTableau_ = rKButcherTableau<Scalar>();
95  numSteps_ = 0;
96 }
97 
98 template<class Scalar>
100 {
101  isInitialized_ = false;
102  model_ = Teuchos::null;
103  solver_ = Teuchos::null;
104  irk_W_factory_ = Teuchos::null;
105  paramList_ = Teuchos::null;
106  //basePoint_;
107  x_ = Teuchos::null;
108  x_old_ = Teuchos::null;
109  x_dot_ = Teuchos::null;
110  //timeRange_;
111  irkModel_ = Teuchos::null;
112  irkButcherTableau_ = Teuchos::null;
113  isDirk_ = false;
114  numSteps_ = -1;
115  haveInitialCondition_ = false;
116  x_stage_bar_ = Teuchos::null;
117 }
118 
119 template<class Scalar>
121  const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > &irk_W_factory
122  )
123 {
124  TEUCHOS_ASSERT( !is_null(irk_W_factory) );
125  irk_W_factory_ = irk_W_factory;
126 }
127 
128 template<class Scalar>
129 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > ImplicitRKStepper<Scalar>::get_W_factory() const
130 {
131  return irk_W_factory_;
132 }
133 
134 // Overridden from SolverAcceptingStepperBase
135 
136 
137 template<class Scalar>
139  const RCP<Thyra::NonlinearSolverBase<Scalar> > &solver
140  )
141 {
142  TEUCHOS_TEST_FOR_EXCEPT(is_null(solver));
143  solver_ = solver;
144 }
145 
146 
147 template<class Scalar>
148 RCP<Thyra::NonlinearSolverBase<Scalar> >
150 {
151  return solver_;
152 }
153 
154 
155 template<class Scalar>
156 RCP<const Thyra::NonlinearSolverBase<Scalar> >
158 {
159  return solver_;
160 }
161 
162 
163 // Overridden from StepperBase
164 
165 
166 template<class Scalar>
168 {
169  return true;
170 }
171 
172 template<class Scalar>
174 {
175  return true;
176 }
177 
178 
179 template<class Scalar>
180 RCP<StepperBase<Scalar> >
182 {
183  // Just use the interface to clone the algorithm in a basically
184  // uninitialized state
185  RCP<ImplicitRKStepper<Scalar> >
186  stepper = Teuchos::rcp(new ImplicitRKStepper<Scalar>());
187 
188  if (!is_null(model_)) {
189  stepper->setModel(model_); // Shallow copy is okay!
190  }
191 
192  if (!is_null(irkButcherTableau_)) {
193  // 06/16/09 tscoffe: should we clone the RKBT here?
194  stepper->setRKButcherTableau(irkButcherTableau_);
195  }
196 
197  if (!is_null(solver_)) {
198  stepper->setSolver(solver_->cloneNonlinearSolver().assert_not_null());
199  }
200 
201  if (!is_null(irk_W_factory_)) {
202  // 06/16/09 tscoffe: should we clone the W_factory here?
203  stepper->set_W_factory(irk_W_factory_);
204  }
205 
206  if (!is_null(paramList_)) {
207  stepper->setParameterList(Teuchos::parameterList(*paramList_));
208  }
209 
210  return stepper;
211 }
212 
213 template<class Scalar>
214 RCP<StepControlStrategyBase<Scalar> > ImplicitRKStepper<Scalar>::getNonconstStepControlStrategy()
215 {
216  return(stepControl_);
217 }
218 
219 template<class Scalar>
220 RCP<const StepControlStrategyBase<Scalar> > ImplicitRKStepper<Scalar>::getStepControlStrategy() const
221 {
222  return(stepControl_);
223 }
224 
225 template<class Scalar>
227 {
228  TEUCHOS_TEST_FOR_EXCEPTION(stepControl == Teuchos::null,std::logic_error,"Error, stepControl == Teuchos::null!\n");
229  stepControl_ = stepControl;
230 }
231 
232 template<class Scalar>
234  const RCP<const Thyra::ModelEvaluator<Scalar> >& model
235  )
236 {
237  TEUCHOS_TEST_FOR_EXCEPT(is_null(model));
238  assertValidModel( *this, *model );
239  model_ = model;
240 }
241 
242 
243 template<class Scalar>
245  const RCP<Thyra::ModelEvaluator<Scalar> >& model
246  )
247 {
248  this->setModel(model); // TODO 09/09/09 tscoffe: use ConstNonconstObjectContainer!
249 }
250 
251 
252 template<class Scalar>
253 RCP<const Thyra::ModelEvaluator<Scalar> >
255 {
256  return model_;
257 }
258 
259 
260 template<class Scalar>
261 RCP<Thyra::ModelEvaluator<Scalar> >
263 {
264  return Teuchos::null;
265 }
266 
267 
268 template<class Scalar>
270  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
271  )
272 {
273 
274  typedef ScalarTraits<Scalar> ST;
275  typedef Thyra::ModelEvaluatorBase MEB;
276 
277  basePoint_ = initialCondition;
278 
279  // x
280 
281  RCP<const Thyra::VectorBase<Scalar> >
282  x_init = initialCondition.get_x();
283 
284 #ifdef HAVE_RYTHMOS_DEBUG
285  TEUCHOS_TEST_FOR_EXCEPTION(
286  is_null(x_init), std::logic_error,
287  "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
288  "then x can not be null!" );
289 #endif
290 
291  x_ = x_init->clone_v();
292 
293  // x_dot
294 
295  x_dot_ = createMember(x_->space());
296 
297  RCP<const Thyra::VectorBase<Scalar> >
298  x_dot_init = initialCondition.get_x_dot();
299 
300  if (!is_null(x_dot_init))
301  assign(x_dot_.ptr(),*x_dot_init);
302  else
303  assign(x_dot_.ptr(),ST::zero());
304 
305  // t
306 
307  const Scalar t =
308  (
309  initialCondition.supports(MEB::IN_ARG_t)
310  ? initialCondition.get_t()
311  : ST::zero()
312  );
313 
314  timeRange_ = timeRange(t,t);
315 
316  // x_old
317  x_old_ = x_->clone_v();
318 
319  haveInitialCondition_ = true;
320 
321 }
322 
323 
324 template<class Scalar>
325 Thyra::ModelEvaluatorBase::InArgs<Scalar>
327 {
328  return basePoint_;
329 }
330 
331 
332 template<class Scalar>
333 Scalar ImplicitRKStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepSizeType)
334 {
335  Scalar stepSizeTaken;
336  using Teuchos::as;
337  using Teuchos::incrVerbLevel;
338  typedef Thyra::NonlinearSolverBase<Scalar> NSB;
339  typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
340 
341  RCP<FancyOStream> out = this->getOStream();
342  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
343  Teuchos::OSTab ostab(out,1,"takeStep");
344  VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
345 
346  // not needed for this
347  int desiredOrder;
348 
349  if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
350  *out
351  << "\nEntering "
352  << Teuchos::TypeNameTraits<ImplicitRKStepper<Scalar> >::name()
353  << "::takeStep("<<dt<<","<<toString(stepSizeType)<<") ...\n";
354  }
355 
356  if (!isInitialized_) {
357  initialize_();
358  }
359  if (stepSizeType == STEP_TYPE_FIXED) {
360  stepSizeTaken = takeFixedStep_(dt , stepSizeType);
361  return stepSizeTaken;
362  } else {
363  isVariableStep_ = true;
364  stepControl_->setOStream(out);
365  stepControl_->setVerbLevel(verbLevel);
366 
367  rkNewtonConvergenceStatus_ = -1;
368 
369  while (rkNewtonConvergenceStatus_ < 0){
370 
371  stepControl_->setRequestedStepSize(*this, dt, stepSizeType);
372  stepControl_->nextStepSize(*this, &dt, &stepSizeType, &desiredOrder);
373 
374  stepSizeTaken = takeVariableStep_(dt, stepSizeType);
375 
376  }
377  return stepSizeTaken;
378  }
379 
380 }
381 
382 
383 template<class Scalar>
384 Scalar ImplicitRKStepper<Scalar>::takeFixedStep_(Scalar dt, StepSizeType stepSizeType)
385 {
386  using Teuchos::as;
387  using Teuchos::incrVerbLevel;
388  typedef ScalarTraits<Scalar> ST;
389  typedef Thyra::NonlinearSolverBase<Scalar> NSB;
390  typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
391 
392  RCP<FancyOStream> out = this->getOStream();
393  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
394  Teuchos::OSTab ostab(out,1,"takeStep");
395  VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
396 
397  if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
398  *out
399  << "\nEntering "
400  << Teuchos::TypeNameTraits<ImplicitRKStepper<Scalar> >::name()
401  << "::takeFixedStep_("<<dt<<","<<toString(stepSizeType)<<") ...\n";
402  }
403 
404  if (!isInitialized_) {
405  initialize_();
406  }
407 
408  TEUCHOS_TEST_FOR_EXCEPT( stepSizeType != STEP_TYPE_FIXED ); // ToDo: Handle variable case later
409 
410  // A) Set up the IRK ModelEvaluator so that it can represent the time step
411  // equation to be solved.
412 
413  // Set irkModel_ with x_old_, t_old_, and dt
414  V_V( x_old_.ptr(), *x_ );
415  Scalar current_dt = dt;
416  Scalar t = timeRange_.upper();
417 
418  // B) Solve the timestep equation
419 
420  // Set the guess for the stage derivatives to zero (unless we can think of
421  // something better)
422  V_S( Teuchos::rcp_dynamic_cast<Thyra::VectorBase<Scalar> >(x_stage_bar_).ptr(), ST::zero() );
423 
424  if (!isDirk_) { // General Implicit RK Case:
425  RCP<ImplicitRKModelEvaluator<Scalar> > firkModel_ =
426  Teuchos::rcp_dynamic_cast<ImplicitRKModelEvaluator<Scalar> >(irkModel_,true);
427  firkModel_->setTimeStepPoint( x_old_, t, current_dt );
428 
429  // Solve timestep equation
430  solver_->solve( &*x_stage_bar_ );
431 
432  } else { // Diagonal Implicit RK Case:
433 
434  RCP<DiagonalImplicitRKModelEvaluator<Scalar> > dirkModel_ =
435  Teuchos::rcp_dynamic_cast<DiagonalImplicitRKModelEvaluator<Scalar> >(irkModel_,true);
436  dirkModel_->setTimeStepPoint( x_old_, t, current_dt );
437  int numStages = irkButcherTableau_->numStages();
438  for (int stage=0 ; stage < numStages ; ++stage) {
439  dirkModel_->setCurrentStage(stage);
440  solver_->solve( &*(x_stage_bar_->getNonconstVectorBlock(stage)) );
441  dirkModel_->setStageSolution( stage, *(x_stage_bar_->getVectorBlock(stage)) );
442  }
443 
444  }
445 
446  // C) Complete the step ...
447 
448  // Combine the stage derivatives with the Butcher tableau "b" vector to obtain the solution at the final time.
449  // x_{k+1} = x_k + dt*sum_{i}^{p}(b_i*x_stage_bar_[i])
450 
451  assembleIRKSolution( irkButcherTableau_->b(), current_dt, *x_old_, *x_stage_bar_,
452  outArg(*x_)
453  );
454 
455  // Update time range
456  timeRange_ = timeRange(t,t+current_dt);
457  numSteps_++;
458 
459  return current_dt;
460 
461 }
462 
463 template<class Scalar>
464 Scalar ImplicitRKStepper<Scalar>::takeVariableStep_(Scalar dt, StepSizeType stepSizeType)
465 {
466  using Teuchos::as;
467  using Teuchos::incrVerbLevel;
468  typedef ScalarTraits<Scalar> ST;
469  typedef Thyra::NonlinearSolverBase<Scalar> NSB;
470  typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
471 
472  RCP<FancyOStream> out = this->getOStream();
473  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
474  Teuchos::OSTab ostab(out,1,"takeStep");
475  VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
476 
477  AttemptedStepStatusFlag status;
478  Scalar dt_old = dt;
479 
480  if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
481  *out
482  << "\nEntering "
483  << Teuchos::TypeNameTraits<ImplicitRKStepper<Scalar> >::name()
484  << "::takeVariableStep_("<<dt<<","<<toString(stepSizeType)<<") ...\n";
485  }
486 
487  if (!isInitialized_) {
488  initialize_();
489  }
490 
491  // A) Set up the IRK ModelEvaluator so that it can represent the time step
492  // equation to be solved.
493 
494  // Set irkModel_ with x_old_, t_old_, and dt
495  V_V( x_old_.ptr(), *x_ );
496  Scalar current_dt = dt;
497  Scalar t = timeRange_.upper();
498  Scalar dt_to_return;
499 
500  // B) Solve the timestep equation
501 
502  // Set the guess for the stage derivatives to zero (unless we can think of
503  // something better)
504  V_S( Teuchos::rcp_dynamic_cast<Thyra::VectorBase<Scalar> >(x_stage_bar_).ptr(), ST::zero() );
505 
506  if (!isDirk_) { // General Implicit RK Case:
507  RCP<ImplicitRKModelEvaluator<Scalar> > firkModel_ =
508  Teuchos::rcp_dynamic_cast<ImplicitRKModelEvaluator<Scalar> >(irkModel_,true);
509  firkModel_->setTimeStepPoint( x_old_, t, current_dt );
510 
511  // Solve timestep equation
512  solver_->solve( &*x_stage_bar_ );
513 
514  } else { // Diagonal Implicit RK Case:
515 
516  RCP<DiagonalImplicitRKModelEvaluator<Scalar> > dirkModel_ =
517  Teuchos::rcp_dynamic_cast<DiagonalImplicitRKModelEvaluator<Scalar> >(irkModel_,true);
518  dirkModel_->setTimeStepPoint( x_old_, t, current_dt );
519  int numStages = irkButcherTableau_->numStages();
520  for (int stage=0 ; stage < numStages ; ++stage) {
521  dirkModel_->setCurrentStage(stage);
522  nonlinearSolveStatus_ = solver_->solve( &*(x_stage_bar_->getNonconstVectorBlock(stage)) );
523 
524  if (nonlinearSolveStatus_.solveStatus == Thyra::SOLVE_STATUS_CONVERGED) {
525  rkNewtonConvergenceStatus_ = 0;
526  } else {
527  rkNewtonConvergenceStatus_ = -1;
528  }
529 
530  // for now setCorrection just sets the rkNewtonConvergenceStatus_ in the stepControl
531  // and this is used by acceptStep method of the stepControl
532 
533  stepControl_->setCorrection(*this, (x_stage_bar_->getNonconstVectorBlock(stage)), (x_stage_bar_->getNonconstVectorBlock(stage)), rkNewtonConvergenceStatus_);
534  bool stepPass = stepControl_->acceptStep(*this, &LETvalue_);
535 
536  if (!stepPass) { // stepPass = false
537  stepLETStatus_ = STEP_LET_STATUS_FAILED;
538  rkNewtonConvergenceStatus_ = -1; // just making sure here
539  break; // leave the for loop
540  } else { // stepPass = true
541  stepLETStatus_ = STEP_LET_STATUS_PASSED;
542  dirkModel_->setStageSolution( stage, *(x_stage_bar_->getVectorBlock(stage)) );
543  rkNewtonConvergenceStatus_ = 0; // just making sure here
544  }
545  }
546  // if none of the stages failed, then I can complete the step
547  }
548 
549  // check the nonlinearSolveStatus
550  if ( rkNewtonConvergenceStatus_ == 0) {
551 
552  /*
553  * if the solver has converged, then I can go ahead and combine the stage solutions
554  * and get the new solution
555  */
556 
557  // C) Complete the step ...
558 
559  // Combine the stage derivatives with the Butcher tableau "b" vector to obtain the solution at the final time.
560  // x_{k+1} = x_k + dt*sum_{i}^{p}(b_i*x_stage_bar_[i])
561 
562  assembleIRKSolution( irkButcherTableau_->b(), current_dt, *x_old_, *x_stage_bar_,
563  outArg(*x_)
564  );
565 
566  // Update time range
567  timeRange_ = timeRange(t,t+current_dt);
568  numSteps_++;
569 
570  // completeStep only if the none of the stage solution's failed to converged
571  stepControl_->completeStep(*this);
572 
573  dt_to_return = current_dt;
574 
575  } else {
576  rkNewtonConvergenceStatus_ = -1;
577  status = stepControl_-> rejectStep(*this); // reject the stage value
578  (void) status; // avoid "set but not used" build warning
579  dt_to_return = dt_old;
580  }
581 
582  return dt_to_return;
583 
584 }
585 
586 
587 template<class Scalar>
589 {
590  StepStatus<Scalar> stepStatus;
591 
592  if (!isInitialized_) {
593  stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
594  stepStatus.message = "This stepper is uninitialized.";
595 // return stepStatus;
596  }
597  else if (numSteps_ > 0) {
598  stepStatus.stepStatus = STEP_STATUS_CONVERGED;
599  }
600  else {
601  stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
602  }
603  stepStatus.stepSize = timeRange_.length();
604  stepStatus.order = irkButcherTableau_->order();
605  stepStatus.time = timeRange_.upper();
606  if(Teuchos::nonnull(x_))
607  stepStatus.solution = x_;
608  else
609  stepStatus.solution = Teuchos::null;
610  stepStatus.solutionDot = Teuchos::null;
611  return(stepStatus);
612 }
613 
614 
615 // Overridden from InterpolationBufferBase
616 
617 
618 template<class Scalar>
619 RCP<const Thyra::VectorSpaceBase<Scalar> >
621 {
622  return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
623 }
624 
625 
626 template<class Scalar>
628  const Array<Scalar>& time_vec
629  ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec
630  ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
631  )
632 {
633  TEUCHOS_TEST_FOR_EXCEPT(true);
634 }
635 
636 
637 template<class Scalar>
639 {
640  return timeRange_;
641 }
642 
643 
644 template<class Scalar>
646  const Array<Scalar>& time_vec
647  ,Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec
648  ,Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
649  ,Array<ScalarMag>* accuracy_vec) const
650 {
651  using Teuchos::constOptInArg;
652  using Teuchos::null;
653  TEUCHOS_ASSERT(haveInitialCondition_);
654  defaultGetPoints<Scalar>(
655  timeRange_.lower(), constOptInArg(*x_old_),
656  Ptr<const VectorBase<Scalar> >(null), // Sun
657  timeRange_.upper(), constOptInArg(*x_),
658  Ptr<const VectorBase<Scalar> >(null), // Sun
659  time_vec,
660  ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
661  Ptr<InterpolatorBase<Scalar> >(null) // For Sun
662  );
663  // 04/17/09 tscoffe: Currently, we don't have x_dot to pass out (TODO)
664 }
665 
666 
667 template<class Scalar>
668 void ImplicitRKStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
669 {
670  TEUCHOS_ASSERT( time_vec != NULL );
671  time_vec->clear();
672  if (!haveInitialCondition_) {
673  return;
674  }
675  time_vec->push_back(timeRange_.lower());
676  if (numSteps_ > 0) {
677  time_vec->push_back(timeRange_.upper());
678  }
679 }
680 
681 
682 template<class Scalar>
683 void ImplicitRKStepper<Scalar>::removeNodes(Array<Scalar>& time_vec)
684 {
685  TEUCHOS_TEST_FOR_EXCEPT(true);
686 }
687 
688 
689 template<class Scalar>
691 {
692  return irkButcherTableau_->order();
693 }
694 
695 
696 // Overridden from Teuchos::ParameterListAcceptor
697 
698 
699 template <class Scalar>
701  RCP<ParameterList> const& paramList
702  )
703 {
704  TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
705  paramList->validateParametersAndSetDefaults(*this->getValidParameters());
706  paramList_ = paramList;
707  Teuchos::readVerboseObjectSublist(&*paramList_,this);
708 }
709 
710 
711 template <class Scalar>
712 RCP<ParameterList>
714 {
715  return(paramList_);
716 }
717 
718 
719 template <class Scalar>
720 RCP<ParameterList>
722 {
723  RCP<ParameterList>
724  temp_param_list = paramList_;
725  paramList_ = Teuchos::null;
726  return(temp_param_list);
727 }
728 
729 
730 template<class Scalar>
731 RCP<const ParameterList>
733 {
734  static RCP<const ParameterList> validPL;
735  if (is_null(validPL)) {
736  RCP<ParameterList> pl = Teuchos::parameterList();
737  if (isVariableStep_){
738  pl->sublist(RythmosStepControlSettings_name);
739  }
740  Teuchos::setupVerboseObjectSublist(&*pl);
741  validPL = pl;
742  }
743  return validPL;
744 }
745 
746 
747 // Overridden from Teuchos::Describable
748 
749 
750 template<class Scalar>
752  FancyOStream &out,
753  const Teuchos::EVerbosityLevel verbLevel
754  ) const
755 {
756  using std::endl;
757  using Teuchos::as;
758  if (!isInitialized_) {
759  out << this->description() << " : This stepper is not initialized yet" << std::endl;
760  return;
761  }
762  if (
763  as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
764  ||
765  as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
766  )
767  {
768  out << this->description() << ":" << endl;
769  Teuchos::OSTab tab(out);
770  out << "model = " << Teuchos::describe(*model_,verbLevel);
771  out << "solver = " << Teuchos::describe(*solver_,verbLevel);
772  out << "irk_W_factory = " << Teuchos::describe(*irk_W_factory_,verbLevel);
773  }
774 }
775 
776 
777 // private
778 
779 
780 template <class Scalar>
782 {
783 
784  // typedef ScalarTraits<Scalar> ST; // unused
785  using Teuchos::rcp_dynamic_cast;
786 
787  TEUCHOS_TEST_FOR_EXCEPT(is_null(model_));
788  TEUCHOS_TEST_FOR_EXCEPT(is_null(solver_));
789  TEUCHOS_TEST_FOR_EXCEPT(irkButcherTableau_->numStages() == 0);
790  TEUCHOS_ASSERT(haveInitialCondition_);
791 
792 #ifdef HAVE_RYTHMOS_DEBUG
793  THYRA_ASSERT_VEC_SPACES(
794  "Rythmos::ImplicitRKStepper::initialize_(...)",
795  *x_->space(), *model_->get_x_space() );
796 #endif
797 
798  if (isVariableStep_ ) {
799  // Initialize StepControl
800  if (stepControl_ == Teuchos::null) {
801  RCP<ImplicitBDFStepperRampingStepControl<Scalar> > rkStepControl =
803  //RCP<StepControlStrategyBase<Scalar> > rkStepControl =
804  //Teuchos::rcp(new StepControlStrategyBase<Scalar>());
805  RCP<Teuchos::ParameterList> stepControlPL =
806  Teuchos::sublist(paramList_ , RythmosStepControlSettings_name);
807  rkStepControl->setParameterList(stepControlPL);
808  this->setStepControlStrategy(rkStepControl);
809  stepControl_->initialize(*this);
810  }
811  }
812 
813  // Set up the IRK mdoel
814 
815  if (!isDirk_) { // General Implicit RK
816  TEUCHOS_TEST_FOR_EXCEPT(is_null(irk_W_factory_));
817  irkModel_ = implicitRKModelEvaluator(
818  model_,basePoint_,irk_W_factory_,irkButcherTableau_);
819  } else { // Diagonal Implicit RK
820  irkModel_ = diagonalImplicitRKModelEvaluator(
821  model_,basePoint_,irk_W_factory_,irkButcherTableau_);
822  }
823 
824  solver_->setModel(irkModel_);
825 
826  // Set up the vector of stage derivatives ...
827  const int numStages = irkButcherTableau_->numStages();
828  RCP<const Thyra::ProductVectorSpaceBase<Scalar> > pvs = productVectorSpace(model_->get_x_space(),numStages);
829  RCP<const Thyra::VectorSpaceBase<Scalar> > vs = rcp_dynamic_cast<const Thyra::VectorSpaceBase<Scalar> >(pvs,true);
830  x_stage_bar_ = rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(createMember(vs),true);
831 // x_stage_bar_ = rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(
832 // createMember(irkModel_->get_x_space())
833 // );
834 
835  isInitialized_ = true;
836 
837 }
838 
839 template <class Scalar>
840 void ImplicitRKStepper<Scalar>::setRKButcherTableau( const RCP<const RKButcherTableauBase<Scalar> > &rkButcherTableau )
841 {
842  TEUCHOS_ASSERT( !is_null(rkButcherTableau) );
843  TEUCHOS_TEST_FOR_EXCEPTION( isInitialized_, std::logic_error,
844  "Error! The RK Butcher Tableau cannot be changed after internal initialization!"
845  );
846  validateIRKButcherTableau(*rkButcherTableau);
847  irkButcherTableau_ = rkButcherTableau;
848  E_RKButcherTableauTypes rkType = determineRKBTType<Scalar>(*irkButcherTableau_);
849  if (
850  (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_DIRK)
851  || (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_SDIRK)
852  || (irkButcherTableau_->numStages() == 1)
853  )
854  {
855  isDirk_ = true;
856  }
857 }
858 
859 template <class Scalar>
860 RCP<const RKButcherTableauBase<Scalar> > ImplicitRKStepper<Scalar>::getRKButcherTableau() const
861 {
862  return irkButcherTableau_;
863 }
864 
865 template<class Scalar>
867 {
868  TEUCHOS_TEST_FOR_EXCEPTION(isInitialized_, std::logic_error,
869  "Error! Cannot change DIRK flag after internal initialization is completed\n"
870  );
871  if (isDirk == true) {
872  E_RKButcherTableauTypes rkType = determineRKBTType<Scalar>(*irkButcherTableau_);
873  bool RKBT_is_DIRK = (
874  (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_DIRK)
875  || (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_SDIRK)
876  || (irkButcherTableau_->numStages() == 1)
877  );
878  TEUCHOS_TEST_FOR_EXCEPTION( !RKBT_is_DIRK, std::logic_error,
879  "Error! Cannot set DIRK flag on a non-DIRK RK Butcher Tableau\n"
880  );
881  } else { // isDirk = false;
882  isDirk_ = isDirk;
883  }
884 }
885 
886 //
887 // Explicit Instantiation macro
888 //
889 // Must be expanded from within the Rythmos namespace!
890 //
891 
892 #define RYTHMOS_IMPLICIT_RK_STEPPER_INSTANT(SCALAR) \
893  \
894  template class ImplicitRKStepper< SCALAR >; \
895  \
896  template RCP< ImplicitRKStepper< SCALAR > > \
897  implicitRKStepper(); \
898  \
899  template RCP< ImplicitRKStepper< SCALAR > > \
900  implicitRKStepper( \
901  const RCP<const Thyra::ModelEvaluator< SCALAR > >& model, \
902  const RCP<Thyra::NonlinearSolverBase< SCALAR > >& solver, \
903  const RCP<Thyra::LinearOpWithSolveFactoryBase< SCALAR > >& irk_W_factory, \
904  const RCP<const RKButcherTableauBase< SCALAR > >& irkbt \
905  );
906 
907 } // namespace Rythmos
908 
909 #endif //Rythmos_IMPLICIT_RK_STEPPER_DEF_H
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)
bool isImplicit() const
Returns true.
RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
RCP< const RKButcherTableauBase< Scalar > > getRKButcherTableau() const
void describe(FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
void setTimeStepPoint(const RCP< const Thyra::VectorBase< Scalar > > &x_old, const Scalar &t_old, const Scalar &delta_t)
Scalar takeStep(Scalar dt, StepSizeType flag)
RCP< const ParameterList > getValidParameters() const
void setRKButcherTableau(const RCP< const RKButcherTableauBase< Scalar > > &rkButcherTableau)
RCP< StepperBase< Scalar > > cloneStepperAlgorithm() const
void setSolver(const RCP< Thyra::NonlinearSolverBase< Scalar > > &solver)
RCP< StepControlStrategyBase< Scalar > > getNonconstStepControlStrategy()
RCP< const StepControlStrategyBase< Scalar > > getStepControlStrategy() const
const StepStatus< Scalar > getStepStatus() const
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
The member functions in the StepControlStrategyBase move you between these states in the following fa...
RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
RCP< const Thyra::VectorBase< Scalar > > solutionDot
void removeNodes(Array< Scalar > &time_vec)
void setInitialCondition(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &initialCondition)
RCP< ImplicitRKStepper< Scalar > > implicitRKStepper()
Nonmember constructor.
void setModel(const RCP< const Thyra::ModelEvaluator< Scalar > > &model)
void setStepControlStrategy(const RCP< StepControlStrategyBase< Scalar > > &stepControlStrategy)
void setNonconstModel(const RCP< Thyra::ModelEvaluator< Scalar > > &model)
void getNodes(Array< Scalar > *time_vec) const
bool supportsCloning() const
Returns true.
RCP< const Thyra::VectorBase< Scalar > > solution
RCP< Thyra::ModelEvaluator< Scalar > > getNonconstModel()
RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
void setTimeStepPoint(const RCP< const Thyra::VectorBase< Scalar > > &x_old, const Scalar &t_old, const Scalar &delta_t)
TimeRange< Scalar > getTimeRange() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getInitialCondition() const
RCP< const Thyra::NonlinearSolverBase< Scalar > > getSolver() const
RCP< Thyra::NonlinearSolverBase< Scalar > > getNonconstSolver()
void setParameterList(RCP< ParameterList > const &paramList)