Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
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
44namespace Rythmos {
45
50template<class Scalar>
51RCP<DefaultIntegrator<Scalar> >
53{
54 RCP<DefaultIntegrator<Scalar> >
55 integrator = Teuchos::rcp(new DefaultIntegrator<Scalar>());
56 return integrator;
57}
58
59
64template<class Scalar>
65RCP<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
83template<class Scalar>
84RCP<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
100template<class Scalar>
101RCP<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
121template<class Scalar>
122const std::string
124
125template<class Scalar>
126const int
128 std::numeric_limits<int>::max();
129
130
131
132// Constructors, Initializers, Misc
133
134
135template<class Scalar>
137 :landOnFinalTime_(true),
138 maxNumTimeSteps_(maxNumTimeSteps_default_),
139 currTimeStepIndex_(-1)
140{}
141
142
143template<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
154template<class Scalar>
155RCP<IntegrationControlStrategyBase<Scalar> >
157{
158 return integrationControlStrategy_;
159}
160
161template<class Scalar>
162RCP<const IntegrationControlStrategyBase<Scalar> >
164{
165 return integrationControlStrategy_;
166}
167
168
169template<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
181template<class Scalar>
183 const RCP<InterpolationBufferAppenderBase<Scalar> > &interpBufferAppender
184 )
185{
186 interpBufferAppender_ = interpBufferAppender.assert_not_null();
187}
188
189
190template<class Scalar>
191RCP<const InterpolationBufferAppenderBase<Scalar> >
193{
194 return interpBufferAppender_;
195}
196
197template<class Scalar>
198RCP<InterpolationBufferAppenderBase<Scalar> >
200{
201 return interpBufferAppender_;
202}
203
204template<class Scalar>
205RCP<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
217template<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
231template<class Scalar>
232RCP<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
250template<class Scalar>
251RCP<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
288template<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
317template<class Scalar>
319{
320 RCP<StepperBase<Scalar> > stepper_temp = stepper_;
321 stepper_ = Teuchos::null;
322 return(stepper_temp);
323}
324
325
326template<class Scalar>
327RCP<const StepperBase<Scalar> > DefaultIntegrator<Scalar>::getStepper() const
328{
329 return(stepper_);
330}
331
332
333template<class Scalar>
334RCP<StepperBase<Scalar> > DefaultIntegrator<Scalar>::getNonconstStepper() const
335{
336 return(stepper_);
337}
338
339
340template<class Scalar>
342 const RCP<InterpolationBufferBase<Scalar> > &trailingInterpBuffer
343 )
344{
345 trailingInterpBuffer_ = trailingInterpBuffer.assert_not_null();
346}
347
348
349template<class Scalar>
350RCP<InterpolationBufferBase<Scalar> >
352{
353 return trailingInterpBuffer_;
354}
355
356
357template<class Scalar>
358RCP<const InterpolationBufferBase<Scalar> >
360{
361 return trailingInterpBuffer_;
362}
363
364template<class Scalar>
365RCP<InterpolationBufferBase<Scalar> >
367{
368 RCP<InterpolationBufferBase<Scalar> > trailingInterpBuffer;
369 std::swap( trailingInterpBuffer, trailingInterpBuffer_ );
370 return trailingInterpBuffer;
371}
372
373
374template<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
508template<class Scalar>
511{
512 return timeRange(
513 stepper_->getTimeRange().lower(),
514 integrationTimeDomain_.upper()
515 );
516}
517
518
519// Overridden from InterpolationBufferBase
520
521
522template<class Scalar>
523RCP<const Thyra::VectorSpaceBase<Scalar> >
525{
526 return stepper_->get_x_space();
527}
528
529
530template<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
541template<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
566template<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
577template<class Scalar>
578void DefaultIntegrator<Scalar>::getNodes(Array<Scalar>* time_vec) const
579{
580 stepper_->getNodes(time_vec);
581}
582
583
584template<class Scalar>
585void DefaultIntegrator<Scalar>::removeNodes(Array<Scalar>& time_vec)
586{
587 stepper_->removeNodes(time_vec);
588}
589
590
591template<class Scalar>
593{
594 return stepper_->getOrder();
595}
596
597
598// private
599
600
601template<class Scalar>
603{
604 if (!is_null(trailingInterpBuffer_) && is_null(interpBufferAppender_))
605 interpBufferAppender_ = pointwiseInterpolationBufferAppender<Scalar>();
606 // ToDo: Do other setup?
607}
608
609
610template<class Scalar>
611bool 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
A concrete subclass for IntegratorBase that allows a good deal of customization.
RCP< DefaultIntegrator< Scalar > > defaultIntegrator(const RCP< IntegrationControlStrategyBase< Scalar > > &integrationControlStrategy, const RCP< IntegrationObserverBase< Scalar > > &integrationObserver)
void setStepper(const RCP< StepperBase< Scalar > > &stepper, const Scalar &finalTime, const bool landOnFinalTime=true)
RCP< const IntegrationControlStrategyBase< Scalar > > getIntegrationControlStrategy() const
RCP< InterpolationBufferBase< Scalar > > getNonconstTrailingInterpolationBuffer()
RCP< StepperBase< Scalar > > getNonconstStepper() const
RCP< InterpolationBufferAppenderBase< Scalar > > getNonconstInterpolationBufferAppender()
RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
void removeNodes(Array< Scalar > &time_vec)
void setTrailingInterpolationBuffer(const RCP< InterpolationBufferBase< Scalar > > &trailingInterpBuffer)
RCP< const StepperBase< Scalar > > getStepper() const
RCP< DefaultIntegrator< Scalar > > defaultIntegrator()
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< StepperBase< Scalar > > unSetStepper()
RCP< IntegratorBase< Scalar > > cloneIntegrator() const
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)
RCP< DefaultIntegrator< Scalar > > observedDefaultIntegrator(const RCP< IntegrationObserverBase< Scalar > > &integrationObserver)
void setIntegrationObserver(const RCP< IntegrationObserverBase< Scalar > > &integrationObserver)
void getNodes(Array< Scalar > *time_vec) const
void setInterpolationBufferAppender(const RCP< InterpolationBufferAppenderBase< Scalar > > &interpBufferAppender)
RCP< const ParameterList > getValidParameters() const
RCP< const InterpolationBufferAppenderBase< Scalar > > getInterpolationBufferAppender()
void setIntegrationControlStrategy(const RCP< IntegrationControlStrategyBase< Scalar > > &integrationControlStrategy)
TimeRange< Scalar > getTimeRange() const
void setParameterList(RCP< ParameterList > const &paramList)
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< InterpolationBufferBase< Scalar > > unSetTrailingInterpolationBuffer()
RCP< IntegrationControlStrategyBase< Scalar > > getNonconstIntegrationControlStrategy()
RCP< const InterpolationBufferBase< Scalar > > getTrailingInterpolationBuffer() const
TimeRange< Scalar > getFwdTimeRange() const
RCP< InterpolationBufferAppenderBase< Scalar > > unSetInterpolationBufferAppender()
RCP< DefaultIntegrator< Scalar > > controlledDefaultIntegrator(const RCP< IntegrationControlStrategyBase< Scalar > > &integrationControlStrategy)
Base class for strategy objects that control integration by selecting step sizes for a stepper.
Base class for strategy objects that observe and time integration by observing the stepper object.
Base class for strategy objects that append data from one InterplationBufferBase object to another.
Base class for an interpolation buffer.
Simple struct to aggregate integration/stepper control information.
Base class for defining stepper functionality.
Represent a time range.