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
390 typedef InterpolationBufferBase<Scalar> IBB;
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>
509TimeRange<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)
RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
void removeNodes(Array< Scalar > &time_vec)
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< DefaultIntegrator< Scalar > > observedDefaultIntegrator(const RCP< IntegrationObserverBase< Scalar > > &integrationObserver)
void getNodes(Array< Scalar > *time_vec) const
void setIntegrationControlStrategy(const RCP< IntegrationControlStrategyBase< Scalar > > &integrationControlStrategy)
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< DefaultIntegrator< Scalar > > controlledDefaultIntegrator(const RCP< IntegrationControlStrategyBase< Scalar > > &integrationControlStrategy)