Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_ImplicitBDFStepperRampingStepControl_def.hpp
1//@HEADER
2// ***********************************************************************
3//
4// Rythmos Package
5// Copyright (2006) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
25//
26// ***********************************************************************
27//@HEADER
28
29#ifndef Rythmos_IMPLICITBDF_STEPPER_RAMPING_STEP_CONTROL_DEF_H
30#define Rythmos_IMPLICITBDF_STEPPER_RAMPING_STEP_CONTROL_DEF_H
31
32#include "Rythmos_ImplicitBDFStepper.hpp"
33#include "Rythmos_ImplicitBDFStepperErrWtVecCalc.hpp"
34#include "Teuchos_StandardParameterEntryValidators.hpp"
35#include <string>
36#include <algorithm>
37
38namespace Rythmos {
39
40template<class Scalar>
41ImplicitBDFStepperRampingStepControl<Scalar>::
42ImplicitBDFStepperRampingStepControl() :
43 stepControlState_(UNINITIALIZED)
44{
45
46}
47
48template<class Scalar>
49void ImplicitBDFStepperRampingStepControl<Scalar>::setStepControlState_(
50 StepControlStrategyState newState)
51{
52 if (stepControlState_ == UNINITIALIZED) {
53 TEUCHOS_TEST_FOR_EXCEPT(newState != BEFORE_FIRST_STEP);
54 } else if (stepControlState_ == BEFORE_FIRST_STEP) {
55 TEUCHOS_TEST_FOR_EXCEPT(newState != MID_STEP);
56 } else if (stepControlState_ == MID_STEP) {
57 TEUCHOS_TEST_FOR_EXCEPT(newState != AFTER_CORRECTION);
58 } else if (stepControlState_ == AFTER_CORRECTION) {
59 TEUCHOS_TEST_FOR_EXCEPT(newState != READY_FOR_NEXT_STEP);
60 } else if (stepControlState_ == READY_FOR_NEXT_STEP) {
61 TEUCHOS_TEST_FOR_EXCEPT(newState != MID_STEP);
62 }
63 stepControlState_ = newState;
64}
65
66template<class Scalar>
67StepControlStrategyState
69{
70 return(stepControlState_);
71}
72
73template<class Scalar>
75{
76 TEUCHOS_TEST_FOR_EXCEPT(!((stepControlState_ == BEFORE_FIRST_STEP) ||
77 (stepControlState_ == READY_FOR_NEXT_STEP)));
78
79 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
80 "updateCoeffs_() is not implemented!");
81}
82
83template<class Scalar>
85 const StepperBase<Scalar>& stepper)
86{
87 // Initialize can be called from the stepper when setInitialCondition
88 // is called.
89 using Teuchos::as;
90 typedef Teuchos::ScalarTraits<Scalar> ST;
91 using Thyra::createMember;
92
93 // Set initial time:
94 TimeRange<Scalar> stepperRange = stepper.getTimeRange();
95 TEUCHOS_TEST_FOR_EXCEPTION(
96 !stepperRange.isValid(),
97 std::logic_error,
98 "Error, Stepper does not have valid time range for initialization "
99 "of ImplicitBDFStepperRampingStepControl!\n");
100
101 if (is_null(parameterList_)) {
102 RCP<Teuchos::ParameterList> emptyParameterList =
103 Teuchos::rcp(new Teuchos::ParameterList);
104 this->setParameterList(emptyParameterList);
105 }
106
107 if (is_null(errWtVecCalc_)) {
108 RCP<ImplicitBDFStepperErrWtVecCalc<Scalar> > IBDFErrWtVecCalc =
109 rcp(new ImplicitBDFStepperErrWtVecCalc<Scalar>());
110 errWtVecCalc_ = IBDFErrWtVecCalc;
111 }
112
113 stepControlState_ = UNINITIALIZED;
114
115 requestedStepSize_ = Scalar(-1.0);
116 currentStepSize_ = initialStepSize_;
117 currentOrder_ = 1;
118 nextStepSize_ = initialStepSize_;
119 nextOrder_ = 1;
120 numberOfSteps_ = 0;
121 totalNumberOfFailedSteps_ = 0;
122 countOfConstantStepsAfterFailure_ = 0;
123
124 if (is_null(delta_)) {
125 delta_ = createMember(stepper.get_x_space());
126 }
127 if (is_null(errWtVec_)) {
128 errWtVec_ = createMember(stepper.get_x_space());
129 }
130 V_S(delta_.ptr(),ST::zero());
131
132 if ( doOutput_(Teuchos::VERB_HIGH) ) {
133 RCP<Teuchos::FancyOStream> out = this->getOStream();
134 Teuchos::OSTab ostab(out,1,"initialize");
135 *out << "currentOrder_ = " << currentOrder_ << std::endl;
136 *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
137 }
138
139 if (breakPoints_.size() > 0) {
140 currentBreakPoints_.clear();
141 for (const auto& bp : breakPoints_) {
142 if (bp > stepperRange.lower())
143 currentBreakPoints_.push_back(bp);
144 }
145 }
146
147 setStepControlState_(BEFORE_FIRST_STEP);
148
149}
150
151template<class Scalar>
153 const StepperBase<Scalar>& stepper,
154 const Scalar& stepSize,
155 const StepSizeType& stepSizeType)
156{
157 typedef Teuchos::ScalarTraits<Scalar> ST;
158
159 TEUCHOS_TEST_FOR_EXCEPT(!((stepControlState_ == UNINITIALIZED) ||
160 (stepControlState_ == BEFORE_FIRST_STEP) ||
161 (stepControlState_ == READY_FOR_NEXT_STEP) ||
162 (stepControlState_ == MID_STEP)));
163
164 TEUCHOS_TEST_FOR_EXCEPTION(
165 ((stepSizeType == STEP_TYPE_FIXED) && (stepSize == ST::zero())),
166 std::logic_error,
167 "Error, step size type == STEP_TYPE_FIXED, "
168 "but requested step size == 0!\n");
169
170 bool didInitialization = false;
171 if (stepControlState_ == UNINITIALIZED) {
172 initialize(stepper);
173 didInitialization = true;
174 }
175
176 // errWtVecSet_ is called during initialize
177 if (!didInitialization) {
178 const ImplicitBDFStepper<Scalar>& implicitBDFStepper =
179 Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
180 const Thyra::VectorBase<Scalar>& xHistory =
181 implicitBDFStepper.getxHistory(0);
182 errWtVecCalc_->errWtVecSet(&*errWtVec_,xHistory,relErrTol_,absErrTol_);
183 }
184
185 requestedStepSize_ = stepSize;
186 stepSizeType_ = stepSizeType;
187}
188
189template<class Scalar>
191 const StepperBase<Scalar>& stepper, Scalar* stepSize,
192 StepSizeType* stepSizeType, int* order)
193{
194 TEUCHOS_TEST_FOR_EXCEPT(!((stepControlState_ == BEFORE_FIRST_STEP) ||
195 (stepControlState_ == MID_STEP) ||
196 (stepControlState_ == READY_FOR_NEXT_STEP) )
197 );
198
199 if (stepControlState_ == BEFORE_FIRST_STEP) {
200 nextStepSize_ = initialStepSize_;
201 nextOrder_ = 1;
202 }
203
204 // Now starting a step - rotate next values into current values
205 if (stepSizeType_ == STEP_TYPE_FIXED)
206 currentStepSize_ = requestedStepSize_;
207 else
208 currentStepSize_ = nextStepSize_;
209
210 currentOrder_ = nextOrder_;
211
212 // Limit the step size to the requested step size
213 currentStepSize_ = std::min(requestedStepSize_, currentStepSize_);
214
215 // Cut if a break point is in range
216 bool hitBreakPoint = false;
217 if (currentBreakPoints_.size() > 0) {
218 // Break points are sorted and in range. Remove as we hit them.
219 const Scalar time = stepper.getStepStatus().time;
220 if (time < *currentBreakPoints_.begin() && (time + currentStepSize_) >= *currentBreakPoints_.begin()) {
221 currentStepSize_ = *currentBreakPoints_.begin() - time;
222 currentBreakPoints_.pop_front();
223 hitBreakPoint = true;
224 }
225 }
226
227 *stepSize = currentStepSize_;
228 *stepSizeType = stepSizeType_;
229 *order = currentOrder_;
230
231 if (stepControlState_ != MID_STEP) {
232 setStepControlState_(MID_STEP);
233 }
234
235 // Output
236 if (doOutput_(Teuchos::VERB_MEDIUM)){
237 Teuchos::FancyOStream& out = *this->getOStream();
238 Teuchos::OSTab ostab1(out,2,"** nextStepSize_ **");
239 out << "Values returned to stepper:" << std::endl;
240 Teuchos::OSTab ostab2(out,2,"** nextStepSize_ **");
241 out << "currentStepSize_ = " << currentStepSize_ << std::endl;
242 out << "currentOrder_ = " << currentOrder_ << std::endl;
243 out << "requestedStepSize_ = " << requestedStepSize_ << std::endl;
244 if (breakPoints_.size() > 0) {
245 if (hitBreakPoint)
246 out << "On break point = true" << std::endl;
247 else
248 out << "On break point = false" << std::endl;
249 }
250 if (currentBreakPoints_.size() > 0)
251 out << "Next break point = " << *currentBreakPoints_.begin() << std::endl;
252 }
253
254}
255
256template<class Scalar>
258 const StepperBase<Scalar>& /* stepper */
259 ,const RCP<const Thyra::VectorBase<Scalar> >& /* soln */
260 ,const RCP<const Thyra::VectorBase<Scalar> >& ee
261 ,int solveStatus)
262{
263 TEUCHOS_TEST_FOR_EXCEPT(stepControlState_ != MID_STEP);
264
265 TEUCHOS_TEST_FOR_EXCEPTION(is_null(ee), std::logic_error,
266 "Error, ee == Teuchos::null!\n");
267
268 ee_ = ee;
269
270 newtonConvergenceStatus_ = solveStatus;
271
272 if ( doOutput_(Teuchos::VERB_MEDIUM) && newtonConvergenceStatus_ < 0) {
273 RCP<Teuchos::FancyOStream> out = this->getOStream();
274 Teuchos::OSTab ostab(out,1,"setCorrection");
275 *out << "\nImplicitBDFStepperRampingStepControl::setCorrection(): "
276 << "Nonlinear Solver Failed!\n";
277 }
278
279 setStepControlState_(AFTER_CORRECTION);
280}
281
282template<class Scalar>
284 const StepperBase<Scalar>& /* stepper */, Scalar* LETValue)
285{
286 using Teuchos::as;
287 typedef Teuchos::ScalarTraits<Scalar> ST;
288
289 TEUCHOS_TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
290
291
292 if ( doOutput_(Teuchos::VERB_HIGH) ) {
293 RCP<Teuchos::FancyOStream> out = this->getOStream();
294 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
295 Teuchos::OSTab ostab(out,1,"acceptStep");
296 *out << "ee_ = " << std::endl;
297 ee_->describe(*out,verbLevel);
298 *out << "errWtVec_ = " << std::endl;
299 errWtVec_->describe(*out,verbLevel);
300 }
301
302 Scalar enorm = wRMSNorm_(*errWtVec_,*ee_);
303
304 Scalar LET = ck_ * enorm;
305
306 if (LETValue) {
307 *LETValue = LET;
308 *LETValue = Scalar(0.0);
309 }
310
311 if (newtonConvergenceStatus_ < 0)
312 return false;
313
314 bool return_status = false;
315
316 if (LET < ST::one() || !useLETToDetermineConvergence_)
317 return_status = true;
318
319 if ( doOutput_(Teuchos::VERB_HIGH) ) {
320 RCP<Teuchos::FancyOStream> out = this->getOStream();
321 Teuchos::OSTab ostab(out,1,"acceptStep");
322 *out << "return_status = " << return_status << std::endl;
323 *out << "Local Truncation Error Check: (ck*enorm) < 1: (" << LET
324 << ") <?= 1" << std::endl;
325 if ( doOutput_(Teuchos::VERB_EXTREME) ) {
326 *out << "ck_ = " << ck_ << std::endl;
327 *out << "enorm = " << enorm << std::endl;
328 }
329 }
330
331 return(return_status);
332}
333
334template<class Scalar>
336 const StepperBase<Scalar>& stepper)
337{
338 TEUCHOS_TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
339 using Teuchos::as;
340 // typedef Teuchos::ScalarTraits<Scalar> ST; // unused
341
342 if ( doOutput_(Teuchos::VERB_HIGH) ) {
343 RCP<Teuchos::FancyOStream> out = this->getOStream();
344
345 Teuchos::OSTab ostab1(out,2,"completeStep_");
346 *out << "\n** Begin completeStep() **" << std::endl;
347 Teuchos::OSTab ostab2(out,2,"** Begin completeStep_ **");
348 *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
349 *out << "numConstantSteps_ = " << numConstantSteps_ << std::endl;
350 *out << "currentStepSize_ = " << currentStepSize_ << std::endl;
351 *out << "nextStepSize_ = " << nextStepSize_ << std::endl;
352 *out << "currentOrder_ = " << currentOrder_ << std::endl;
353 *out << "nextOrder_ = " << nextOrder_ << std::endl;
354 *out << "stepSizeIncreaseFactor_ = " << stepSizeIncreaseFactor_ <<std::endl;
355 *out << "countOfConstantStepsAfterFailure_ = "
356 << countOfConstantStepsAfterFailure_ << std::endl;
357 }
358
359 numberOfSteps_ ++;
360
361 if (countOfConstantStepsAfterFailure_ > 0) {
362 // We track the number of consecutive time step failures so that
363 // if we have a bunch of nonlinear failures, lets keep the time
364 // step constant for a while before we start to ramp again. This
365 // keeps us from oscillating between ramping and cutting step
366 // sizes and wasting resources.
367
368 nextStepSize_ = currentStepSize_;
369 nextOrder_ = currentOrder_;
370
371 // Decrement failure counter
372 countOfConstantStepsAfterFailure_ =
373 std::max( (countOfConstantStepsAfterFailure_ - 1), 0);
374
375 if ( doOutput_(Teuchos::VERB_HIGH) ) {
376 RCP<Teuchos::FancyOStream> out = this->getOStream();
377 Teuchos::OSTab ostab(out,1,"completeStep_");
378 *out << "\nNext Step Size held constant due to previous failed steps!\n";
379 *out << "countOfConstantStepsAfterFailure_ = "
380 << countOfConstantStepsAfterFailure_ << std::endl;
381 }
382 }
383 else {
384
385 // Phase 1: Constant step size at 1st order
386 if (numberOfSteps_ < numConstantSteps_) {
387 if (currentStepSize_ < initialStepSize_)
388 nextStepSize_ = std::min(initialStepSize_,
389 currentStepSize_ * stepSizeIncreaseFactor_);
390 nextOrder_ = 1;
391 }
392 // Phase 2: Constant step size, ramping the order
393 else if (currentOrder_ < maxOrder_) {
394 if (currentStepSize_ < initialStepSize_)
395 nextStepSize_ = std::min(initialStepSize_,
396 currentStepSize_ * stepSizeIncreaseFactor_);
397 else
398 nextStepSize_ = currentStepSize_;
399
400 nextOrder_ = currentOrder_ + 1;
401 }
402 // Phase 3: Ramping dt to max step size, highest order
403 else if ( (numberOfSteps_ >= numConstantSteps_) &&
404 (currentOrder_ == maxOrder_) ) {
405 nextStepSize_ = std::min(maxStepSize_,
406 currentStepSize_ * stepSizeIncreaseFactor_);
407 nextOrder_ = maxOrder_;
408 }
409 else {
410 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
411 "RampingStepControlStrategy logic is broken. Please contact "
412 "developers. Aborting run!");
413 }
414
415 if (restrictStepSizeByNumberOfNonlinearIterations_) {
416 const Rythmos::ImplicitBDFStepper<Scalar>* ibdfStepper =
417 dynamic_cast<const Rythmos::ImplicitBDFStepper<Scalar>* >(&stepper);
418 TEUCHOS_ASSERT(ibdfStepper != NULL);
419 TEUCHOS_ASSERT(nonnull(ibdfStepper->getNonlinearSolveStatus().extraParameters));
420 int numberOfNonlinearIterations = ibdfStepper->getNonlinearSolveStatus().extraParameters->template get<int>("Number of Iterations");
421 if (numberOfNonlinearIterations >= numberOfNonlinearIterationsForStepSizeRestriction_) {
422 nextStepSize_ = currentStepSize_;
423 }
424 }
425
426
427 } // if (countOfConstantStepsAfterFailure_ > 0)
428
429 setStepControlState_(READY_FOR_NEXT_STEP);
430
431 if ( doOutput_(Teuchos::VERB_HIGH) ) {
432 RCP<Teuchos::FancyOStream> out = this->getOStream();
433 Teuchos::OSTab ostab1(out,2,"** completeStep_ **");
434 *out << "** End of completeStep() **" << std::endl;
435 Teuchos::OSTab ostab2(out,2,"** End completeStep_ **");
436 *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
437 *out << "numConstantSteps_ = " << numConstantSteps_ << std::endl;
438 *out << "currentStepSize_ = " << currentStepSize_ << std::endl;
439 *out << "nextStepSize_ = " << nextStepSize_ << std::endl;
440 *out << "currentOrder_ = " << currentOrder_ << std::endl;
441 *out << "nextOrder_ = " << nextOrder_ << std::endl;
442 *out << "stepSizeIncreaseFactor_ = " << stepSizeIncreaseFactor_ <<std::endl;
443 *out << "countOfConstantStepsAfterFailure_ = "
444 << countOfConstantStepsAfterFailure_ << std::endl;
445 }
446}
447
448template<class Scalar>
449AttemptedStepStatusFlag
451 const StepperBase<Scalar>& /* stepper */)
452{
453 TEUCHOS_TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
454
455 using Teuchos::as;
456
457 ++totalNumberOfFailedSteps_;
458 ++countOfConstantStepsAfterFailure_;
459
460 // If time step size is already at the min time step, then quit
461 if (currentStepSize_ <= minStepSize_)
462 return (REP_ERR_FAIL);
463
464 // Otherwise, cut the time step and keep order the same
465 nextStepSize_ = std::max(minStepSize_,
466 (currentStepSize_ * stepSizeDecreaseFactor_) );
467
468 setStepControlState_(READY_FOR_NEXT_STEP);
469
470 return(PREDICT_AGAIN);
471}
472
473template<class Scalar>
475 Teuchos::FancyOStream &out,
476 const Teuchos::EVerbosityLevel verbLevel
477 ) const
478{
479
480 using Teuchos::as;
481
482 if ( (as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT) ) ||
483 (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
484 ) {
485 out << this->description() << "::describe" << std::endl;
486 }
487 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
488 out << "currentStepSize_ = " << currentStepSize_ << std::endl;
489 out << "currentOrder_ = " << currentOrder_ << std::endl;
490 }
491 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
492 }
493 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
494 out << "ee_ = ";
495 if (ee_ == Teuchos::null) {
496 out << "Teuchos::null" << std::endl;
497 } else {
498 ee_->describe(out,verbLevel);
499 }
500 out << "delta_ = ";
501 if (delta_ == Teuchos::null) {
502 out << "Teuchos::null" << std::endl;
503 } else {
504 delta_->describe(out,verbLevel);
505 }
506 out << "errWtVec_ = ";
507 if (errWtVec_ == Teuchos::null) {
508 out << "Teuchos::null" << std::endl;
509 } else {
510 errWtVec_->describe(out,verbLevel);
511 }
512 }
513}
514
515template<class Scalar>
517 RCP<Teuchos::ParameterList> const& paramList
518 )
519{
520
521 using Teuchos::as;
522 // typedef Teuchos::ScalarTraits<Scalar> ST; // unused
523
524 TEUCHOS_TEST_FOR_EXCEPT(paramList == Teuchos::null);
525
526 parameterList_ = Teuchos::parameterList(*paramList);
527
528 parameterList_->validateParametersAndSetDefaults(*this->getValidParameters());
529
530 Teuchos::ParameterList& p = *parameterList_;
531
532 numConstantSteps_ = p.get<int>("Number of Constant First Order Steps");
533 initialStepSize_ = p.get<Scalar>("Initial Step Size");
534 minStepSize_ = p.get<Scalar>("Min Step Size");
535 maxStepSize_ = p.get<Scalar>("Max Step Size");
536 stepSizeIncreaseFactor_ = p.get<Scalar>("Step Size Increase Factor");
537 stepSizeDecreaseFactor_ = p.get<Scalar>("Step Size Decrease Factor");
538
539 minOrder_ = p.get<int>("Min Order");
540 TEUCHOS_TEST_FOR_EXCEPTION(
541 !((1 <= minOrder_) && (minOrder_ <= 5)), std::logic_error,
542 "Error, minOrder_ = " << minOrder_ << " is not in range [1,5]!\n"
543 );
544 maxOrder_ = p.get<int>("Max Order");
545 TEUCHOS_TEST_FOR_EXCEPTION(
546 !((1 <= maxOrder_) && (maxOrder_ <= 5)), std::logic_error,
547 "Error, maxOrder_ = " << maxOrder_ << " is not in range [1,5]!\n"
548 );
549
550 absErrTol_ = p.get<Scalar>("Absolute Error Tolerance");
551 relErrTol_ = p.get<Scalar>("Relative Error Tolerance");
552
553 {
554 std::string let_acceptance =
555 p.get<std::string>("Use LET To Determine Step Acceptance");
556 useLETToDetermineConvergence_ = (let_acceptance == "TRUE");
557
558 // Currently the using LET for step acceptance is not supported
559 // since we can't calculate the LETValue. Once this is
560 // implemented, delete the line below.
561 TEUCHOS_TEST_FOR_EXCEPTION(useLETToDetermineConvergence_, std::logic_error,
562 "Error - the flag \"Use LET To Determine Step Acceptance\" is set to "
563 "\"TRUE\" but the local error computation is currently not supported. "
564 "Please set this flag to \"FALSE\" for now.");
565 }
566
567 if (p.get<std::string>("Restrict Step Size Increase by Number of Nonlinear Iterations") == "TRUE")
568 restrictStepSizeByNumberOfNonlinearIterations_ = true;
569 else if (p.get<std::string>("Restrict Step Size Increase by Number of Nonlinear Iterations") == "FALSE")
570 restrictStepSizeByNumberOfNonlinearIterations_ = false;
571
572 numberOfNonlinearIterationsForStepSizeRestriction_ =
573 p.get<int>("Number of Nonlinear Iterations for Step Size Restriction");
574
575 // Parse the break points
576 {
577 breakPoints_.clear();
578 std::string str = p.get<std::string>("Break Points");
579 std::string delimiters(",");
580 std::string::size_type lastPos = str.find_first_not_of(delimiters, 0);
581 std::string::size_type pos = str.find_first_of(delimiters, lastPos);
582 while ( (pos != std::string::npos) || (lastPos != std::string::npos) ) {
583 std::string token = str.substr(lastPos,pos-lastPos);
584 breakPoints_.push_back(Scalar(std::stod(token)));
585 if(pos==std::string::npos)
586 break;
587
588 lastPos = str.find_first_not_of(delimiters, pos);
589 pos = str.find_first_of(delimiters, lastPos);
590 }
591
592 // order the break points
593 std::sort(breakPoints_.begin(),breakPoints_.end());
594
595 // copy into current
596 currentBreakPoints_.resize(breakPoints_.size());
597 std::copy(breakPoints_.begin(),breakPoints_.end(),currentBreakPoints_.begin());
598 }
599
600 if ( doOutput_(Teuchos::VERB_HIGH) ) {
601 RCP<Teuchos::FancyOStream> out = this->getOStream();
602 Teuchos::OSTab ostab(out,1,"setParameterList");
603 out->precision(15);
604 *out << "minOrder_ = " << minOrder_ << std::endl;
605 *out << "maxOrder_ = " << maxOrder_ << std::endl;
606 *out << "relErrTol_ = " << relErrTol_ << std::endl;
607 *out << "absErrTol_ = " << absErrTol_ << std::endl;
608 *out << "stepSizeType = " << stepSizeType_ << std::endl;
609 *out << "stopTime_ = " << stopTime_ << std::endl;
610 }
611
612}
613
614template<class Scalar>
615RCP<const Teuchos::ParameterList>
617{
618 using Teuchos::RCP;
619 using Teuchos::rcp;
620 using Teuchos::ParameterList;
621
622 static RCP<ParameterList> p;
623
624 if (is_null(p)) {
625
626 p = rcp(new ParameterList);
627
628 p->set<int>("Number of Constant First Order Steps", 10,
629 "Number of constant steps to take before handing control to "
630 "variable stepper.");
631 p->set<Scalar>("Initial Step Size", Scalar(1.0e-3),
632 "Initial time step size and target step size to take during the "
633 "initial constant step phase (could be reduced due to step failures).");
634 p->set<Scalar>("Min Step Size", Scalar(1.0e-7), "Minimum time step size.");
635 p->set<Scalar>("Max Step Size", Scalar(1.0), "Maximum time step size.");
636 p->set<Scalar>("Step Size Increase Factor", Scalar(1.2),
637 "Time step growth factor used after a successful time step. dt_{n+1} = "
638 "(increase factor) * dt_n");
639 p->set<Scalar>("Step Size Decrease Factor", Scalar(0.5),
640 "Time step reduction factor used for a failed time step. dt_{n+1} = "
641 "(decrease factor) * dt_n");
642 p->set<int>("Min Order", 1, "Minimum order to run at.");
643 p->set<int>("Max Order", 5, "Maximum order to run at.");
644 p->set<Scalar>("Absolute Error Tolerance", Scalar(1.0e-5),
645 "abstol value used in WRMS calculation.");
646 p->set<Scalar>("Relative Error Tolerance", Scalar(1.0e-3),
647 "reltol value used in WRMS calculation.");
648 Teuchos::setStringToIntegralParameter<int>(
649 "Use LET To Determine Step Acceptance",
650 "FALSE",
651 "If set to TRUE, then acceptance of step dependes on LET in addition "
652 "to Nonlinear solver converging.",
653 Teuchos::tuple<std::string>("TRUE","FALSE"),
654 p.get());
655 Teuchos::setStringToIntegralParameter<int>(
656 "Restrict Step Size Increase by Number of Nonlinear Iterations",
657 "FALSE",
658 "If set to TRUE, then the step size will not be allowed to increase "
659 "if the number of nonlinear iterations was greater than or equal to the "
660 "specified value.",
661 Teuchos::tuple<std::string>("TRUE","FALSE"),
662 p.get());
663 p->set<int>("Number of Nonlinear Iterations for Step Size Restriction",
664 2,
665 "If \" Restrct Step Size Increase by Number of Nonlinear Iterations\" is "
666 "true, the step size will not be allowed to increase if the number of nonlinear "
667 "iterations was greater than or equal to the specified value.");
668 p->set<std::string>("Break Points","");
669 }
670
671 return (p);
672}
673
674template<class Scalar>
675RCP<Teuchos::ParameterList>
677{
678 RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
679 parameterList_ = Teuchos::null;
680 return(temp_param_list);
681}
682
683template<class Scalar>
684RCP<Teuchos::ParameterList>
686{
687 return(parameterList_);
688}
689
690template<class Scalar>
692 const StepperBase<Scalar>& stepper)
693{
694 if (stepControlState_ == UNINITIALIZED) {
695 initialize(stepper);
696 }
697 const ImplicitBDFStepper<Scalar>& bdfstepper =
698 Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
699 int desiredOrder = bdfstepper.getOrder();
700 TEUCHOS_TEST_FOR_EXCEPT(!((1 <= desiredOrder) &&
701 (desiredOrder <= maxOrder_)));
702 if (stepControlState_ == BEFORE_FIRST_STEP) {
703 TEUCHOS_TEST_FOR_EXCEPTION(
704 desiredOrder > 1,
705 std::logic_error,
706 "Error, this ImplicitBDF stepper has not taken a step yet, so it "
707 "cannot take a step of order " << desiredOrder << " > 1!\n");
708 }
709 TEUCHOS_TEST_FOR_EXCEPT(!(desiredOrder <= currentOrder_+1));
710 currentOrder_ = desiredOrder;
711
712 if ( doOutput_(Teuchos::VERB_EXTREME) ) {
713 RCP<Teuchos::FancyOStream> out = this->getOStream();
714 Teuchos::OSTab ostab(out,1,"setStepControlData");
715 *out << "currentOrder_ = " << currentOrder_ << std::endl;
716 }
717}
718
719template<class Scalar>
721{
722 return true;
723}
724
725
726template<class Scalar>
727RCP<StepControlStrategyBase<Scalar> >
729{
730
731 RCP<ImplicitBDFStepperRampingStepControl<Scalar> > stepControl =
733
734 if (!is_null(parameterList_)) {
735 stepControl->setParameterList(parameterList_);
736 }
737
738 return stepControl;
739}
740
741template<class Scalar>
743 const RCP<ErrWtVecCalcBase<Scalar> >& errWtVecCalc)
744{
745 TEUCHOS_TEST_FOR_EXCEPT(is_null(errWtVecCalc));
746 errWtVecCalc_ = errWtVecCalc;
747}
748
749template<class Scalar>
750RCP<const ErrWtVecCalcBase<Scalar> >
752{
753 return(errWtVecCalc_);
754}
755
756template<class Scalar>
758 const Thyra::VectorBase<Scalar>& weight,
759 const Thyra::VectorBase<Scalar>& vector) const
760{
761 return(norm_2(weight,vector));
762}
763
764template<class Scalar>
766{
767 TEUCHOS_TEST_FOR_EXCEPTION(
768 stepControlState_ == UNINITIALIZED, std::logic_error,
769 "Error, attempting to call getMinOrder before intiialization!\n"
770 );
771 return(minOrder_);
772}
773
774template<class Scalar>
776{
777 TEUCHOS_TEST_FOR_EXCEPTION(
778 stepControlState_ == UNINITIALIZED, std::logic_error,
779 "Error, attempting to call getMaxOrder before initialization!\n"
780 );
781 return(maxOrder_);
782}
783
784template<class Scalar>
786 Teuchos::EVerbosityLevel verbLevel)
787{
788 Teuchos::EVerbosityLevel currentObjectVerbLevel = this->getVerbLevel();
789
790 if ( Teuchos::as<int>(currentObjectVerbLevel) >= Teuchos::as<int>(verbLevel) )
791 return true;
792
793 return false;
794}
795
796template<class Scalar>
797int ImplicitBDFStepperRampingStepControl<Scalar>::numberOfSteps() const
798{
799 return numberOfSteps_;
800}
801
802template<class Scalar>
803int ImplicitBDFStepperRampingStepControl<Scalar>::numberOfFailedSteps() const
804{
805 return totalNumberOfFailedSteps_;
806}
807
808template<class Scalar>
809Scalar ImplicitBDFStepperRampingStepControl<Scalar>::currentStepSize() const
810{
811 return currentStepSize_;
812}
813
814template<class Scalar>
815int ImplicitBDFStepperRampingStepControl<Scalar>::currentOrder() const
816{
817 return currentOrder_;
818}
819
820
821//
822// Explicit Instantiation macro
823//
824// Must be expanded from within the Rythmos namespace!
825//
826
827#define RYTHMOS_IMPLICITBDF_STEPPER_RAMPING_STEPCONTROL_INSTANT(SCALAR) \
828 template class ImplicitBDFStepperRampingStepControl< SCALAR >;
829
830
831} // namespace Rythmos
832#endif // Rythmos_IMPLICITBDF_STEPPER_RAMPING_STEP_CONTROL_DEF_H
833
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
const Thyra::SolveStatus< Scalar > & getNonlinearSolveStatus() const
Returns the Thyra::SolveStatus object from the last nonlinear solve.
virtual RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const =0
Return the space for x and x_dot.
virtual TimeRange< Scalar > getTimeRange() const =0
Return the range of time values where interpolation calls can be performed.
Base class for defining stepper functionality.
virtual const StepStatus< Scalar > getStepStatus() const =0
Get current stepper status after a step has been taken.
Represent a time range.