Rythmos - Transient Integration for Differential Equations  Version of the Day
 All Classes Functions Variables Typedefs Pages
Rythmos_DefaultIntegrator_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_DEFAULT_INTEGRATOR_DEF_HPP
30 #define RYTHMOS_DEFAULT_INTEGRATOR_DEF_HPP
31 
32 #include "Rythmos_DefaultIntegrator_decl.hpp"
33 #include "Rythmos_InterpolationBufferHelpers.hpp"
34 #include "Rythmos_IntegrationControlStrategyBase.hpp"
35 #include "Rythmos_IntegrationObserverBase.hpp"
36 #include "Rythmos_InterpolationBufferAppenderBase.hpp"
37 #include "Rythmos_PointwiseInterpolationBufferAppender.hpp"
38 #include "Rythmos_StepperHelpers.hpp"
39 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
40 #include "Teuchos_Assert.hpp"
41 #include "Teuchos_as.hpp"
42 #include <limits>
43 
44 namespace Rythmos {
45 
50 template<class Scalar>
51 RCP<DefaultIntegrator<Scalar> >
53 {
54  RCP<DefaultIntegrator<Scalar> >
55  integrator = Teuchos::rcp(new DefaultIntegrator<Scalar>());
56  return integrator;
57 }
58 
59 
64 template<class Scalar>
65 RCP<DefaultIntegrator<Scalar> >
67  const RCP<IntegrationControlStrategyBase<Scalar> > &integrationControlStrategy,
68  const RCP<IntegrationObserverBase<Scalar> > &integrationObserver
69  )
70 {
71  RCP<DefaultIntegrator<Scalar> >
72  integrator = Teuchos::rcp(new DefaultIntegrator<Scalar>());
73  integrator->setIntegrationControlStrategy(integrationControlStrategy);
74  integrator->setIntegrationObserver(integrationObserver);
75  return integrator;
76 }
77 
78 
83 template<class Scalar>
84 RCP<DefaultIntegrator<Scalar> >
86  const RCP<IntegrationControlStrategyBase<Scalar> > &integrationControlStrategy
87  )
88 {
89  RCP<DefaultIntegrator<Scalar> >
90  integrator = Teuchos::rcp(new DefaultIntegrator<Scalar>());
91  integrator->setIntegrationControlStrategy(integrationControlStrategy);
92  return integrator;
93 }
94 
95 
100 template<class Scalar>
101 RCP<DefaultIntegrator<Scalar> >
103  const RCP<IntegrationObserverBase<Scalar> > &integrationObserver
104  )
105 {
106  RCP<DefaultIntegrator<Scalar> >
107  integrator = Teuchos::rcp(new DefaultIntegrator<Scalar>());
108  integrator->setIntegrationObserver(integrationObserver);
109  return integrator;
110 }
111 
112 
113 
114 // ////////////////////////////
115 // Defintions
116 
117 
118 // Static data members
119 
120 
121 template<class Scalar>
122 const std::string
123 DefaultIntegrator<Scalar>::maxNumTimeSteps_name_ = "Max Number Time Steps";
124 
125 template<class Scalar>
126 const int
128  std::numeric_limits<int>::max();
129 
130 
131 
132 // Constructors, Initializers, Misc
133 
134 
135 template<class Scalar>
137  :landOnFinalTime_(true),
138  maxNumTimeSteps_(maxNumTimeSteps_default_),
139  currTimeStepIndex_(-1)
140 {}
141 
142 
143 template<class Scalar>
145  const RCP<IntegrationControlStrategyBase<Scalar> > &integrationControlStrategy
146  )
147 {
148 #ifdef HAVE_RYTHMOS_DEBUG
149  TEUCHOS_TEST_FOR_EXCEPT(is_null(integrationControlStrategy));
150 #endif
151  integrationControlStrategy_ = integrationControlStrategy;
152 }
153 
154 template<class Scalar>
155 RCP<IntegrationControlStrategyBase<Scalar> >
157 {
158  return integrationControlStrategy_;
159 }
160 
161 template<class Scalar>
162 RCP<const IntegrationControlStrategyBase<Scalar> >
164 {
165  return integrationControlStrategy_;
166 }
167 
168 
169 template<class Scalar>
171  const RCP<IntegrationObserverBase<Scalar> > &integrationObserver
172  )
173 {
174 #ifdef HAVE_RYTHMOS_DEBUG
175  TEUCHOS_TEST_FOR_EXCEPT(is_null(integrationObserver));
176 #endif
177  integrationObserver_ = integrationObserver;
178 }
179 
180 
181 template<class Scalar>
183  const RCP<InterpolationBufferAppenderBase<Scalar> > &interpBufferAppender
184  )
185 {
186  interpBufferAppender_ = interpBufferAppender.assert_not_null();
187 }
188 
189 
190 template<class Scalar>
191 RCP<const InterpolationBufferAppenderBase<Scalar> >
193 {
194  return interpBufferAppender_;
195 }
196 
197 template<class Scalar>
198 RCP<InterpolationBufferAppenderBase<Scalar> >
200 {
201  return interpBufferAppender_;
202 }
203 
204 template<class Scalar>
205 RCP<InterpolationBufferAppenderBase<Scalar> >
207 {
208  RCP<InterpolationBufferAppenderBase<Scalar> > InterpBufferAppender;
209  std::swap( InterpBufferAppender, interpBufferAppender_ );
210  return InterpBufferAppender;
211 }
212 
213 
214 // Overridden from ParameterListAcceptor
215 
216 
217 template<class Scalar>
219  RCP<ParameterList> const& paramList
220  )
221 {
222  TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
223  paramList->validateParameters(*getValidParameters());
224  this->setMyParamList(paramList);
225  maxNumTimeSteps_ = paramList->get(
226  maxNumTimeSteps_name_, maxNumTimeSteps_default_);
227  Teuchos::readVerboseObjectSublist(&*paramList,this);
228 }
229 
230 
231 template<class Scalar>
232 RCP<const ParameterList>
234 {
235  static RCP<const ParameterList> validPL;
236  if (is_null(validPL) ) {
237  RCP<ParameterList> pl = Teuchos::parameterList();
238  pl->set(maxNumTimeSteps_name_, maxNumTimeSteps_default_,
239  "Set the maximum number of integration time-steps allowed.");
240  Teuchos::setupVerboseObjectSublist(&*pl);
241  validPL = pl;
242  }
243  return validPL;
244 }
245 
246 
247 // Overridden from IntegratorBase
248 
249 
250 template<class Scalar>
251 RCP<IntegratorBase<Scalar> >
253 {
254 
255  using Teuchos::null;
256  RCP<DefaultIntegrator<Scalar> >
257  newIntegrator = Teuchos::rcp(new DefaultIntegrator<Scalar>());
258  // Only copy control information, not the stepper or the model it contains!
259  newIntegrator->stepper_ = Teuchos::null;
260  const RCP<const ParameterList> paramList = this->getParameterList();
261  if (!is_null(paramList)) {
262  newIntegrator->setParameterList(Teuchos::parameterList(*paramList));
263  }
264  if (!is_null(integrationControlStrategy_)) {
265  newIntegrator->integrationControlStrategy_ =
266  integrationControlStrategy_->cloneIntegrationControlStrategy().assert_not_null();
267  }
268  if (!is_null(integrationObserver_)) {
269  newIntegrator->integrationObserver_ =
270  integrationObserver_->cloneIntegrationObserver().assert_not_null();
271  }
272  if (!is_null(trailingInterpBuffer_)) {
273  // ToDo: implement the clone!
274  newIntegrator->trailingInterpBuffer_ = null;
275  //newIntegrator->trailingInterpBuffer_ =
276  // trailingInterpBuffer_->cloneInterpolationBuffer().assert_not_null();
277  }
278  if (!is_null(interpBufferAppender_)) {
279  // ToDo: implement the clone!
280  newIntegrator->interpBufferAppender_ = null;
281  //newIntegrator->interpBufferAppender_ =
282  // interpBufferAppender_->cloneInterpolationBufferAppender().assert_not_null();
283  }
284  return newIntegrator;
285 }
286 
287 
288 template<class Scalar>
290  const RCP<StepperBase<Scalar> > &stepper,
291  const Scalar &finalTime,
292  const bool landOnFinalTime
293  )
294 {
295  typedef Teuchos::ScalarTraits<Scalar> ST;
296  TEUCHOS_TEST_FOR_EXCEPT(is_null(stepper));
297  TEUCHOS_TEST_FOR_EXCEPT( finalTime < stepper->getTimeRange().lower() );
298  TEUCHOS_ASSERT( stepper->getTimeRange().length() == ST::zero() );
299  // 2007/07/25: rabartl: ToDo: Validate state of the stepper!
300  stepper_ = stepper;
301  integrationTimeDomain_ = timeRange(stepper_->getTimeRange().lower(),
302  finalTime);
303  landOnFinalTime_ = landOnFinalTime;
304  currTimeStepIndex_ = 0;
305  stepCtrlInfoLast_ = StepControlInfo<Scalar>();
306  if (!is_null(integrationControlStrategy_))
307  integrationControlStrategy_->resetIntegrationControlStrategy(
308  integrationTimeDomain_
309  );
310  if (!is_null(integrationObserver_))
311  integrationObserver_->resetIntegrationObserver(
312  integrationTimeDomain_
313  );
314 }
315 
316 
317 template<class Scalar>
318 RCP<StepperBase<Scalar> > DefaultIntegrator<Scalar>::unSetStepper()
319 {
320  RCP<StepperBase<Scalar> > stepper_temp = stepper_;
321  stepper_ = Teuchos::null;
322  return(stepper_temp);
323 }
324 
325 
326 template<class Scalar>
327 RCP<const StepperBase<Scalar> > DefaultIntegrator<Scalar>::getStepper() const
328 {
329  return(stepper_);
330 }
331 
332 
333 template<class Scalar>
334 RCP<StepperBase<Scalar> > DefaultIntegrator<Scalar>::getNonconstStepper() const
335 {
336  return(stepper_);
337 }
338 
339 
340 template<class Scalar>
342  const RCP<InterpolationBufferBase<Scalar> > &trailingInterpBuffer
343  )
344 {
345  trailingInterpBuffer_ = trailingInterpBuffer.assert_not_null();
346 }
347 
348 
349 template<class Scalar>
350 RCP<InterpolationBufferBase<Scalar> >
352 {
353  return trailingInterpBuffer_;
354 }
355 
356 
357 template<class Scalar>
358 RCP<const InterpolationBufferBase<Scalar> >
360 {
361  return trailingInterpBuffer_;
362 }
363 
364 template<class Scalar>
365 RCP<InterpolationBufferBase<Scalar> >
367 {
368  RCP<InterpolationBufferBase<Scalar> > trailingInterpBuffer;
369  std::swap( trailingInterpBuffer, trailingInterpBuffer_ );
370  return trailingInterpBuffer;
371 }
372 
373 
374 template<class Scalar>
376  const Array<Scalar>& time_vec,
377  Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
378  Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
379  Array<ScalarMag>* accuracy_vec)
380 {
381 
382  RYTHMOS_FUNC_TIME_MONITOR_DIFF("Rythmos:DefaultIntegrator::getFwdPoints",
383  TopLevel);
384 
385  using Teuchos::incrVerbLevel;
386 #ifndef _MSC_VER
387  using Teuchos::Describable;
388 #endif
389  // typedef Teuchos::ScalarTraits<Scalar> ST; // unused
391  typedef Teuchos::VerboseObjectTempState<IBB> VOTSIBB;
392 
393  finalizeSetup();
394 
395  RCP<Teuchos::FancyOStream> out = this->getOStream();
396  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
397  Teuchos::OSTab tab(out);
398  VOTSIBB stepper_outputTempState(stepper_,out,incrVerbLevel(verbLevel,-1));
399 
400  if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
401  *out << "\nEntering " << this->Describable::description()
402  << "::getFwdPoints(...) ...\n"
403  << "\nStepper: " << Teuchos::describe(*stepper_,verbLevel);
404 
405  if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
406  *out << "\nRequested time points: " << Teuchos::toString(time_vec) << "\n";
407 
408  // Observe start of a time integration
409  if (!is_null(integrationObserver_)) {
410  integrationObserver_->setOStream(out);
411  integrationObserver_->setVerbLevel(incrVerbLevel(verbLevel,-1));
412  integrationObserver_->observeStartTimeIntegration(*stepper_);
413  }
414 
415  //
416  // 0) Initial setup
417  //
418 
419  const int numTimePoints = time_vec.size();
420 
421  // Assert preconditions
422  assertTimePointsAreSorted(time_vec);
423  TEUCHOS_TEST_FOR_EXCEPT(accuracy_vec!=0); // ToDo: Remove accuracy_vec!
424 
425  // Resize the storage for the output arrays
426  if (x_vec)
427  x_vec->resize(numTimePoints);
428  if (xdot_vec)
429  xdot_vec->resize(numTimePoints);
430 
431  // This int records the next time point offset in time_vec[timePointIndex]
432  // that needs to be handled. This gets updated as the time points are
433  // filled below.
434  int nextTimePointIndex = 0;
435 
436  assertNoTimePointsBeforeCurrentTimeRange(*this,time_vec,nextTimePointIndex);
437 
438  //
439  // 1) First, get all time points that fall within the current time range
440  //
441 
442  {
443  RYTHMOS_FUNC_TIME_MONITOR(
444  "Rythmos:DefaultIntegrator::getFwdPoints: getPoints");
445  // 2007/10/05: rabartl: ToDo: Get points from trailingInterpBuffer_ first!
446  getCurrentPoints(*stepper_,time_vec,x_vec,xdot_vec,&nextTimePointIndex);
447  }
448 
449  //
450  // 2) Advance the stepper to satisfy time points in time_vec that fall
451  // before the current time.
452  //
453 
454  while ( nextTimePointIndex < numTimePoints ) {
455 
456  // Use the time stepping algorithm to step up to or past the next
457  // requested time but not so far as to step past the point entirely.
458  const Scalar t = time_vec[nextTimePointIndex];
459  bool advanceStepperToTimeSucceeded = false;
460 
461  // Make sure t is not beyond 'Final Time'.
462  if ( integrationTimeDomain_.upper() < t ) {
463  *out << "\nWARNING: Skipping time point, t = " << t
464  << ", because it is beyond 'Final Time' = "
465  << integrationTimeDomain_.upper() << "\n";
466  // Break out of the while loop and attempt to exit gracefully.
467  break;
468  } else {
469  RYTHMOS_FUNC_TIME_MONITOR(
470  "Rythmos:DefaultIntegrator::getFwdPoints: advanceStepperToTime");
471  advanceStepperToTimeSucceeded = advanceStepperToTime(t);
472  }
473  if (!advanceStepperToTimeSucceeded) {
474  bool reachedMaxNumTimeSteps = (currTimeStepIndex_ >= maxNumTimeSteps_);
475  if (reachedMaxNumTimeSteps) {
476  // Break out of the while loop and attempt to exit gracefully.
477  break;
478  }
479  TEUCHOS_TEST_FOR_EXCEPTION(
480  !advanceStepperToTimeSucceeded, Exceptions::GetFwdPointsFailed,
481  this->description() << "\n\n"
482  "Error: The integration failed to get to time " << t
483  << " and only achieved\ngetting to "
484  << stepper_->getTimeRange().upper() << "!"
485  );
486  }
487 
488  // Extract the next set of points (perhaps just one) from the stepper
489  {
490  RYTHMOS_FUNC_TIME_MONITOR(
491  "Rythmos:DefaultIntegrator::getFwdPoints: getPoints (fwd)");
492  getCurrentPoints(*stepper_,time_vec,x_vec,xdot_vec,&nextTimePointIndex);
493  }
494 
495  }
496 
497  // Observe end of a time integration
498  if (!is_null(integrationObserver_)) {
499  integrationObserver_->observeEndTimeIntegration(*stepper_);
500  }
501 
502  if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
503  *out << "\nLeaving " << this->Describable::description()
504  << "::getFwdPoints(...) ...\n";
505 }
506 
507 
508 template<class Scalar>
511 {
512  return timeRange(
513  stepper_->getTimeRange().lower(),
514  integrationTimeDomain_.upper()
515  );
516 }
517 
518 
519 // Overridden from InterpolationBufferBase
520 
521 
522 template<class Scalar>
523 RCP<const Thyra::VectorSpaceBase<Scalar> >
525 {
526  return stepper_->get_x_space();
527 }
528 
529 
530 template<class Scalar>
532  const Array<Scalar>& time_vec,
533  const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
534  const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
535  )
536 {
537  stepper_->addPoints(time_vec,x_vec,xdot_vec);
538 }
539 
540 
541 template<class Scalar>
543  const Array<Scalar>& time_vec,
544  Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
545  Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
546  Array<ScalarMag>* accuracy_vec
547  ) const
548 {
549 // if (nonnull(trailingInterpBuffer_)) {
550 // int nextTimePointIndex = 0;
551 // getCurrentPoints(*trailingInterpBuffer_, time_vec, x_vec, xdot_vec, &nextTimePointIndex);
552 // getCurrentPoints(*stepper_, time_vec, x_vec, xdot_vec, &nextTimePointIndex);
553 // TEUCHOS_TEST_FOR_EXCEPTION( nextTimePointIndex < Teuchos::as<int>(time_vec.size()),
554 // std::out_of_range,
555 // "Error, the time point time_vec["<<nextTimePointIndex<<"] = "
556 // << time_vec[nextTimePointIndex] << " falls outside of the time range "
557 // << getTimeRange() << "!"
558 // );
559 // }
560 // else {
561  stepper_->getPoints(time_vec, x_vec, xdot_vec, accuracy_vec);
562 // }
563 }
564 
565 
566 template<class Scalar>
568 {
569  if (nonnull(trailingInterpBuffer_)) {
570  return timeRange(trailingInterpBuffer_->getTimeRange().lower(),
571  stepper_->getTimeRange().upper());
572  }
573  return stepper_->getTimeRange();
574 }
575 
576 
577 template<class Scalar>
578 void DefaultIntegrator<Scalar>::getNodes(Array<Scalar>* time_vec) const
579 {
580  stepper_->getNodes(time_vec);
581 }
582 
583 
584 template<class Scalar>
585 void DefaultIntegrator<Scalar>::removeNodes(Array<Scalar>& time_vec)
586 {
587  stepper_->removeNodes(time_vec);
588 }
589 
590 
591 template<class Scalar>
593 {
594  return stepper_->getOrder();
595 }
596 
597 
598 // private
599 
600 
601 template<class Scalar>
603 {
604  if (!is_null(trailingInterpBuffer_) && is_null(interpBufferAppender_))
605  interpBufferAppender_ = pointwiseInterpolationBufferAppender<Scalar>();
606  // ToDo: Do other setup?
607 }
608 
609 
610 template<class Scalar>
611 bool DefaultIntegrator<Scalar>::advanceStepperToTime(const Scalar& advance_to_t)
612 {
613 
614  RYTHMOS_FUNC_TIME_MONITOR_DIFF(
615  "Rythmos:DefaultIntegrator::advanceStepperToTime", TopLevel);
616 
617  using std::endl;
618  typedef std::numeric_limits<Scalar> NL;
619  using Teuchos::incrVerbLevel;
620 #ifndef _MSC_VER
621  using Teuchos::Describable;
622 #endif
623  using Teuchos::OSTab;
624  typedef Teuchos::ScalarTraits<Scalar> ST;
625 
626  RCP<Teuchos::FancyOStream> out = this->getOStream();
627  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
628 
629  if (!is_null(integrationControlStrategy_)) {
630  integrationControlStrategy_->setOStream(out);
631  integrationControlStrategy_->setVerbLevel(incrVerbLevel(verbLevel,-1));
632  }
633 
634  if (!is_null(integrationObserver_)) {
635  integrationObserver_->setOStream(out);
636  integrationObserver_->setVerbLevel(incrVerbLevel(verbLevel,-1));
637  }
638 
639  if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
640  *out << "\nEntering " << this->Describable::description()
641  << "::advanceStepperToTime("<<advance_to_t<<") ...\n";
642 
643  // Remember what timestep index we are on so we can report it later
644  const int initCurrTimeStepIndex = currTimeStepIndex_;
645 
646  // Take steps until the requested time is reached (or passed)
647 
648  TimeRange<Scalar> currStepperTimeRange = stepper_->getTimeRange();
649 
650  // Start by assuming we can reach the time advance_to_t
651  bool return_val = true;
652 
653  while ( !currStepperTimeRange.isInRange(advance_to_t) ) {
654 
655  // Halt immediately if exceeded max iterations
656  if (currTimeStepIndex_ >= maxNumTimeSteps_) {
657  if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
658  *out
659  << "\n***"
660  << "\n*** NOTICE: currTimeStepIndex = "<<currTimeStepIndex_
661  << " >= maxNumTimeSteps = "<<maxNumTimeSteps_
662  << ", halting time integration!"
663  << "\n***\n";
664  return_val = false;
665  break; // Exit the loop immediately!
666  }
667 
668  if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) ) {
669  *out << "\nTake step: current_stepper_t = "
670  << currStepperTimeRange.upper()
671  << ", currTimeStepIndex = " << currTimeStepIndex_ << endl;
672  }
673 
674  OSTab tab(out);
675 
676  //
677  // A) Reinitialize if a hard breakpoint was reached on the last time step
678  //
679 
680  if (stepCtrlInfoLast_.limitedByBreakPoint) {
681  if ( stepCtrlInfoLast_.breakPointType == BREAK_POINT_TYPE_HARD ) {
682  RYTHMOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::restart");
683  if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
684  *out << "\nAt a hard-breakpoint, restarting time integrator ...\n";
685  restart(&*stepper_);
686  }
687  else {
688  if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
689  *out <<"\nAt a soft-breakpoint, NOT restarting time integrator ...\n";
690  }
691  }
692 
693  //
694  // B) Find an acceptable time step in a loop
695  //
696  // NOTE: Look for continue statements to iterate the loop!
697  //
698 
699  bool foundAcceptableTimeStep = false;
700  StepControlInfo<Scalar> stepCtrlInfo;
701 
702  // \todo Limit the maximum number of trial time steps to avoid an infinite
703  // loop!
704 
705  while (!foundAcceptableTimeStep) {
706 
707  //
708  // B.1) Get the trial step control info
709  //
710 
711  StepControlInfo<Scalar> trialStepCtrlInfo;
712  {
713  RYTHMOS_FUNC_TIME_MONITOR(
714  "Rythmos:DefaultIntegrator::advanceStepperToTime: getStepCtrl");
715  if (!is_null(integrationControlStrategy_)) {
716  // Let an external strategy object determine the step size and type.
717  // Note that any breakpoint info is also related through this call.
718  trialStepCtrlInfo =
719  integrationControlStrategy_->getNextStepControlInfo(
720  *stepper_, stepCtrlInfoLast_, currTimeStepIndex_);
721  }
722  else {
723  // Take a variable step if we have no control strategy
724  trialStepCtrlInfo.stepType = STEP_TYPE_VARIABLE;
725  trialStepCtrlInfo.stepSize = NL::max();
726  }
727  }
728 
729  // Print the initial trial step
730  if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
731  *out << "\nTrial step:\n";
732  OSTab tab2(out);
733  *out << trialStepCtrlInfo;
734  }
735 
736  // Halt immediately if we were told to do so
737  if (trialStepCtrlInfo.stepSize < ST::zero()) {
738  if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
739  *out
740  << "\n***"
741  << "\n*** NOTICE: The IntegrationControlStrategy object "
742  << "return stepSize < 0.0, halting time integration!"
743  << "\n***\n";
744  return_val = false;
745  break; // Exit the loop immediately!
746  }
747 
748  // Make sure we don't step past the final time if asked not to
749  bool updatedTrialStepCtrlInfo = false;
750  {
751  const Scalar finalTime = integrationTimeDomain_.upper();
752  if (landOnFinalTime_ && trialStepCtrlInfo.stepSize
753  + currStepperTimeRange.upper() > finalTime) {
754 
755  trialStepCtrlInfo.stepSize = finalTime - currStepperTimeRange.upper();
756  updatedTrialStepCtrlInfo = true;
757 
758  if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
759  *out << "\nCutting trial step to "<< trialStepCtrlInfo.stepSize
760  << " to avoid stepping past final time ...\n";
761  if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
762  *out << "\n"
763  << " integrationTimeDomain = " << integrationTimeDomain_<<"\n"
764  << " currStepperTimeRange = " << currStepperTimeRange <<"\n";
765  }
766  }
767  }
768 
769  // Print the modified trial step
770  if ( updatedTrialStepCtrlInfo
771  && includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
772  {
773  *out << "\nUpdated trial step:\n";
774  OSTab tab2(out);
775  *out << trialStepCtrlInfo;
776  }
777 
778  //
779  // B.2) Take the step
780  //
781 
782  // Output to the observer we are starting a step
783  if (!is_null(integrationObserver_))
784  integrationObserver_->observeStartTimeStep(
785  *stepper_, trialStepCtrlInfo, currTimeStepIndex_
786  );
787 
788  // Print step type and size
789  if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
790  if (trialStepCtrlInfo.stepType == STEP_TYPE_VARIABLE)
791  *out << "\nTaking a variable time step with max step size = "
792  << trialStepCtrlInfo.stepSize << " ....\n";
793  else
794  *out << "\nTaking a fixed time step of size = "
795  << trialStepCtrlInfo.stepSize << " ....\n";
796  }
797 
798  // Take step
799  Scalar stepSizeTaken;
800  {
801  RYTHMOS_FUNC_TIME_MONITOR(
802  "Rythmos:Defaultintegrator::advancesteppertotime: takestep");
803  stepSizeTaken = stepper_->takeStep(
804  trialStepCtrlInfo.stepSize, trialStepCtrlInfo.stepType
805  );
806  }
807 
808  // Update info about this step
809  currStepperTimeRange = stepper_->getTimeRange();
810  stepCtrlInfo = stepCtrlInfoTaken(trialStepCtrlInfo,stepSizeTaken);
811 
812  // Print the step actually taken
813  if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
814  *out << "\nStep actually taken:\n";
815  OSTab tab2(out);
816  *out << stepCtrlInfo;
817  }
818 
819  // Determine if the timestep failed
820  const bool timeStepFailed = (stepCtrlInfo.stepSize <= ST::zero());
821  if (timeStepFailed && includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM)) {
822  *out << "\nWARNING: timeStep = "<<trialStepCtrlInfo.stepSize
823  <<" failed!\n";
824  }
825 
826  // Notify observer of a failed time step
827  if (timeStepFailed) {
828  if (nonnull(integrationObserver_))
829  integrationObserver_->observeFailedTimeStep(
830  *stepper_, stepCtrlInfo, currTimeStepIndex_
831  );
832 
833  // Allow the IntegrationControlStrategy object to suggest another timestep
834  const bool handlesFailedTimeSteps =
835  nonnull(integrationControlStrategy_) &&
836  integrationControlStrategy_->handlesFailedTimeSteps();
837  if (handlesFailedTimeSteps)
838  {
839  // See if a new timestep can be suggested
840  if (integrationControlStrategy_->resetForFailedTimeStep(
841  *stepper_, stepCtrlInfoLast_,
842  currTimeStepIndex_, trialStepCtrlInfo) )
843  {
844  if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
845  *out << "\nThe IntegrationControlStrategy object indicated that"
846  << " it would like to suggest another timestep!\n";
847  }
848  // Skip the rest of the code in the loop and back to the top to try
849  // another timestep! Note: By doing this we skip the statement that
850  // sets
851  continue;
852  }
853  else
854  {
855  if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
856  *out << "\nThe IntegrationControlStrategy object could not "
857  << "suggest a better time step! Allowing to fail the "
858  << "time step!\n";
859  }
860  // Fall through to the failure checking!
861  }
862  }
863  }
864 
865  // Validate step taken
866  if (trialStepCtrlInfo.stepType == STEP_TYPE_VARIABLE) {
867  TEUCHOS_TEST_FOR_EXCEPTION(
868  stepSizeTaken < ST::zero(), std::logic_error,
869  "Error, stepper took negative step of dt = " << stepSizeTaken << "!\n"
870  );
871  TEUCHOS_TEST_FOR_EXCEPTION(
872  stepSizeTaken > trialStepCtrlInfo.stepSize, std::logic_error,
873  "Error, stepper took step of dt = " << stepSizeTaken
874  << " > max step size of = " << trialStepCtrlInfo.stepSize << "!\n"
875  );
876  }
877  else { // STEP_TYPE_FIXED
878  TEUCHOS_TEST_FOR_EXCEPTION(
879  stepSizeTaken != trialStepCtrlInfo.stepSize, std::logic_error,
880  "Error, stepper took step of dt = " << stepSizeTaken
881  << " when asked to take step of dt = " << trialStepCtrlInfo.stepSize
882  << "\n");
883  }
884 
885  // If we get here, the timestep is fine and is accepted!
886  foundAcceptableTimeStep = true;
887 
888  // Append the trailing interpolation buffer (if defined)
889  if (!is_null(trailingInterpBuffer_)) {
890  interpBufferAppender_->append(*stepper_,currStepperTimeRange,
891  trailingInterpBuffer_.ptr() );
892  }
893 
894  //
895  // B.3) Output info about step
896  //
897 
898  {
899 
900  RYTHMOS_FUNC_TIME_MONITOR(
901  "Rythmos:DefaultIntegrator::advanceStepperToTime: output");
902 
903  // Print our own brief output
904  if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
905  StepStatus<Scalar> stepStatus = stepper_->getStepStatus();
906  *out << "\nTime point reached = " << stepStatus.time << endl;
907  *out << "\nstepStatus:\n" << stepStatus;
908  if ( includesVerbLevel(verbLevel,Teuchos::VERB_EXTREME) ) {
909  RCP<const Thyra::VectorBase<Scalar> >
910  solution = stepStatus.solution,
911  solutionDot = stepStatus.solutionDot;
912  if (!is_null(solution))
913  *out << "\nsolution = \n"
914  << Teuchos::describe(*solution,verbLevel);
915  if (!is_null(solutionDot))
916  *out << "\nsolutionDot = \n"
917  << Teuchos::describe(*solutionDot,verbLevel);
918  }
919  }
920 
921  // Output to the observer
922  if (!is_null(integrationObserver_))
923  integrationObserver_->observeCompletedTimeStep(
924  *stepper_, stepCtrlInfo, currTimeStepIndex_
925  );
926 
927  }
928 
929  } // end loop to find a valid time step
930 
931  //
932  // C) Update info for next time step
933  //
934 
935  stepCtrlInfoLast_ = stepCtrlInfo;
936  ++currTimeStepIndex_;
937 
938  }
939 
940  if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
941  *out << "\nNumber of steps taken in this call to "
942  << "advanceStepperToTime(...) = "
943  << (currTimeStepIndex_ - initCurrTimeStepIndex) << endl
944  << "\nLeaving " << this->Describable::description()
945  << "::advanceStepperToTime("<<advance_to_t<<") ...\n";
946 
947  return return_val;
948 
949 }
950 
951 //
952 // Explicit Instantiation macro
953 //
954 // Must be expanded from within the Rythmos namespace!
955 //
956 
957 #define RYTHMOS_DEFAULT_INTEGRATOR_INSTANT(SCALAR) \
958  \
959  template class DefaultIntegrator< SCALAR >; \
960  \
961  template RCP<DefaultIntegrator< SCALAR > > \
962  defaultIntegrator(); \
963  \
964  template RCP<DefaultIntegrator< SCALAR > > \
965  defaultIntegrator( \
966  const RCP<IntegrationControlStrategyBase< SCALAR > > &integrationControlStrategy, \
967  const RCP<IntegrationObserverBase< SCALAR > > &integrationObserver \
968  ); \
969  \
970  template RCP<DefaultIntegrator< SCALAR > > \
971  controlledDefaultIntegrator( \
972  const RCP<IntegrationControlStrategyBase< SCALAR > > &integrationControlStrategy \
973  ); \
974  \
975  template RCP<DefaultIntegrator< SCALAR > > \
976  observedDefaultIntegrator( \
977  const RCP<IntegrationObserverBase< SCALAR > > &integrationObserver \
978  );
979 
980 } // namespace Rythmos
981 
982 
983 #endif //RYTHMOS_DEFAULT_INTEGRATOR_DEF_HPP
RCP< const ParameterList > getValidParameters() const
RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() 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
RCP< StepperBase< Scalar > > unSetStepper()
RCP< const InterpolationBufferAppenderBase< Scalar > > getInterpolationBufferAppender()
RCP< InterpolationBufferBase< Scalar > > unSetTrailingInterpolationBuffer()
Base class for defining stepper functionality.
RCP< IntegratorBase< Scalar > > cloneIntegrator() const
RCP< const StepperBase< Scalar > > getStepper() const
RCP< IntegrationControlStrategyBase< Scalar > > getNonconstIntegrationControlStrategy()
Simple struct to aggregate integration/stepper control information.
void getNodes(Array< Scalar > *time_vec) const
Base class for strategy objects that control integration by selecting step sizes for a stepper...
RCP< const IntegrationControlStrategyBase< Scalar > > getIntegrationControlStrategy() const
TimeRange< Scalar > getFwdTimeRange() const
Base class for strategy objects that observe and time integration by observing the stepper object...
RCP< DefaultIntegrator< Scalar > > defaultIntegrator(const RCP< IntegrationControlStrategyBase< Scalar > > &integrationControlStrategy, const RCP< IntegrationObserverBase< Scalar > > &integrationObserver)
RCP< InterpolationBufferAppenderBase< Scalar > > unSetInterpolationBufferAppender()
RCP< DefaultIntegrator< Scalar > > defaultIntegrator()
void setIntegrationControlStrategy(const RCP< IntegrationControlStrategyBase< Scalar > > &integrationControlStrategy)
A concrete subclass for IntegratorBase that allows a good deal of customization.
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 setTrailingInterpolationBuffer(const RCP< InterpolationBufferBase< Scalar > > &trailingInterpBuffer)
RCP< StepperBase< Scalar > > getNonconstStepper() const
Base class for an interpolation buffer.
RCP< const InterpolationBufferBase< Scalar > > getTrailingInterpolationBuffer() const
RCP< DefaultIntegrator< Scalar > > controlledDefaultIntegrator(const RCP< IntegrationControlStrategyBase< Scalar > > &integrationControlStrategy)
void getFwdPoints(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)
Base class for strategy objects that append data from one InterplationBufferBase object to another...
void removeNodes(Array< Scalar > &time_vec)
void setInterpolationBufferAppender(const RCP< InterpolationBufferAppenderBase< Scalar > > &interpBufferAppender)
RCP< InterpolationBufferBase< Scalar > > getNonconstTrailingInterpolationBuffer()
TimeRange< Scalar > getTimeRange() const
RCP< DefaultIntegrator< Scalar > > observedDefaultIntegrator(const RCP< IntegrationObserverBase< Scalar > > &integrationObserver)
RCP< InterpolationBufferAppenderBase< Scalar > > getNonconstInterpolationBufferAppender()
void setIntegrationObserver(const RCP< IntegrationObserverBase< Scalar > > &integrationObserver)
void setParameterList(RCP< ParameterList > const &paramList)
void setStepper(const RCP< StepperBase< Scalar > > &stepper, const Scalar &finalTime, const bool landOnFinalTime=true)