85 const StepperBase<Scalar>& stepper)
90 typedef Teuchos::ScalarTraits<Scalar> ST;
91 using Thyra::createMember;
94 TimeRange<Scalar> stepperRange = stepper.getTimeRange();
95 TEUCHOS_TEST_FOR_EXCEPTION(
96 !stepperRange.isValid(),
98 "Error, Stepper does not have valid time range for initialization "
99 "of ImplicitBDFStepperRampingStepControl!\n");
101 if (is_null(parameterList_)) {
102 RCP<Teuchos::ParameterList> emptyParameterList =
103 Teuchos::rcp(
new Teuchos::ParameterList);
104 this->setParameterList(emptyParameterList);
107 if (is_null(errWtVecCalc_)) {
108 RCP<ImplicitBDFStepperErrWtVecCalc<Scalar> > IBDFErrWtVecCalc =
109 rcp(
new ImplicitBDFStepperErrWtVecCalc<Scalar>());
110 errWtVecCalc_ = IBDFErrWtVecCalc;
113 stepControlState_ = UNINITIALIZED;
115 requestedStepSize_ = Scalar(-1.0);
116 currentStepSize_ = initialStepSize_;
118 nextStepSize_ = initialStepSize_;
121 totalNumberOfFailedSteps_ = 0;
122 countOfConstantStepsAfterFailure_ = 0;
124 if (is_null(delta_)) {
125 delta_ = createMember(stepper.get_x_space());
127 if (is_null(errWtVec_)) {
128 errWtVec_ = createMember(stepper.get_x_space());
130 V_S(delta_.ptr(),ST::zero());
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;
139 if (breakPoints_.size() > 0) {
140 currentBreakPoints_.clear();
141 for (
const auto& bp : breakPoints_) {
142 if (bp > stepperRange.lower())
143 currentBreakPoints_.push_back(bp);
147 setStepControlState_(BEFORE_FIRST_STEP);
153 const StepperBase<Scalar>& stepper,
154 const Scalar& stepSize,
155 const StepSizeType& stepSizeType)
157 typedef Teuchos::ScalarTraits<Scalar> ST;
159 TEUCHOS_TEST_FOR_EXCEPT(!((stepControlState_ == UNINITIALIZED) ||
160 (stepControlState_ == BEFORE_FIRST_STEP) ||
161 (stepControlState_ == READY_FOR_NEXT_STEP) ||
162 (stepControlState_ == MID_STEP)));
164 TEUCHOS_TEST_FOR_EXCEPTION(
165 ((stepSizeType == STEP_TYPE_FIXED) && (stepSize == ST::zero())),
167 "Error, step size type == STEP_TYPE_FIXED, "
168 "but requested step size == 0!\n");
170 bool didInitialization =
false;
171 if (stepControlState_ == UNINITIALIZED) {
173 didInitialization =
true;
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_);
185 requestedStepSize_ = stepSize;
186 stepSizeType_ = stepSizeType;
191 const StepperBase<Scalar>& stepper, Scalar* stepSize,
192 StepSizeType* stepSizeType,
int* order)
194 TEUCHOS_TEST_FOR_EXCEPT(!((stepControlState_ == BEFORE_FIRST_STEP) ||
195 (stepControlState_ == MID_STEP) ||
196 (stepControlState_ == READY_FOR_NEXT_STEP) )
199 if (stepControlState_ == BEFORE_FIRST_STEP) {
200 nextStepSize_ = initialStepSize_;
205 if (stepSizeType_ == STEP_TYPE_FIXED)
206 currentStepSize_ = requestedStepSize_;
208 currentStepSize_ = nextStepSize_;
210 currentOrder_ = nextOrder_;
213 currentStepSize_ = std::min(requestedStepSize_, currentStepSize_);
216 bool hitBreakPoint =
false;
217 if (currentBreakPoints_.size() > 0) {
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;
227 *stepSize = currentStepSize_;
228 *stepSizeType = stepSizeType_;
229 *order = currentOrder_;
231 if (stepControlState_ != MID_STEP) {
232 setStepControlState_(MID_STEP);
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) {
246 out <<
"On break point = true" << std::endl;
248 out <<
"On break point = false" << std::endl;
250 if (currentBreakPoints_.size() > 0)
251 out <<
"Next break point = " << *currentBreakPoints_.begin() << std::endl;
258 const StepperBase<Scalar>&
259 ,
const RCP<
const Thyra::VectorBase<Scalar> >&
260 ,
const RCP<
const Thyra::VectorBase<Scalar> >& ee
263 TEUCHOS_TEST_FOR_EXCEPT(stepControlState_ != MID_STEP);
265 TEUCHOS_TEST_FOR_EXCEPTION(is_null(ee), std::logic_error,
266 "Error, ee == Teuchos::null!\n");
270 newtonConvergenceStatus_ = solveStatus;
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";
279 setStepControlState_(AFTER_CORRECTION);
284 const StepperBase<Scalar>& , Scalar* LETValue)
287 typedef Teuchos::ScalarTraits<Scalar> ST;
289 TEUCHOS_TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
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);
302 Scalar enorm = wRMSNorm_(*errWtVec_,*ee_);
304 Scalar LET = ck_ * enorm;
308 *LETValue = Scalar(0.0);
311 if (newtonConvergenceStatus_ < 0)
314 bool return_status =
false;
316 if (LET < ST::one() || !useLETToDetermineConvergence_)
317 return_status =
true;
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;
331 return(return_status);
336 const StepperBase<Scalar>& stepper)
338 TEUCHOS_TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
342 if ( doOutput_(Teuchos::VERB_HIGH) ) {
343 RCP<Teuchos::FancyOStream> out = this->getOStream();
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;
361 if (countOfConstantStepsAfterFailure_ > 0) {
368 nextStepSize_ = currentStepSize_;
369 nextOrder_ = currentOrder_;
372 countOfConstantStepsAfterFailure_ =
373 std::max( (countOfConstantStepsAfterFailure_ - 1), 0);
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;
386 if (numberOfSteps_ < numConstantSteps_) {
387 if (currentStepSize_ < initialStepSize_)
388 nextStepSize_ = std::min(initialStepSize_,
389 currentStepSize_ * stepSizeIncreaseFactor_);
393 else if (currentOrder_ < maxOrder_) {
394 if (currentStepSize_ < initialStepSize_)
395 nextStepSize_ = std::min(initialStepSize_,
396 currentStepSize_ * stepSizeIncreaseFactor_);
398 nextStepSize_ = currentStepSize_;
400 nextOrder_ = currentOrder_ + 1;
403 else if ( (numberOfSteps_ >= numConstantSteps_) &&
404 (currentOrder_ == maxOrder_) ) {
405 nextStepSize_ = std::min(maxStepSize_,
406 currentStepSize_ * stepSizeIncreaseFactor_);
407 nextOrder_ = maxOrder_;
410 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
411 "RampingStepControlStrategy logic is broken. Please contact "
412 "developers. Aborting run!");
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_;
429 setStepControlState_(READY_FOR_NEXT_STEP);
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;
475 Teuchos::FancyOStream &out,
476 const Teuchos::EVerbosityLevel verbLevel
482 if ( (as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT) ) ||
483 (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
485 out << this->description() <<
"::describe" << std::endl;
487 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
488 out <<
"currentStepSize_ = " << currentStepSize_ << std::endl;
489 out <<
"currentOrder_ = " << currentOrder_ << std::endl;
491 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
493 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
495 if (ee_ == Teuchos::null) {
496 out <<
"Teuchos::null" << std::endl;
498 ee_->describe(out,verbLevel);
501 if (delta_ == Teuchos::null) {
502 out <<
"Teuchos::null" << std::endl;
504 delta_->describe(out,verbLevel);
506 out <<
"errWtVec_ = ";
507 if (errWtVec_ == Teuchos::null) {
508 out <<
"Teuchos::null" << std::endl;
510 errWtVec_->describe(out,verbLevel);
517 RCP<Teuchos::ParameterList>
const& paramList
524 TEUCHOS_TEST_FOR_EXCEPT(paramList == Teuchos::null);
526 parameterList_ = Teuchos::parameterList(*paramList);
528 parameterList_->validateParametersAndSetDefaults(*this->getValidParameters());
530 Teuchos::ParameterList& p = *parameterList_;
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");
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"
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"
550 absErrTol_ = p.get<Scalar>(
"Absolute Error Tolerance");
551 relErrTol_ = p.get<Scalar>(
"Relative Error Tolerance");
554 std::string let_acceptance =
555 p.get<std::string>(
"Use LET To Determine Step Acceptance");
556 useLETToDetermineConvergence_ = (let_acceptance ==
"TRUE");
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.");
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;
572 numberOfNonlinearIterationsForStepSizeRestriction_ =
573 p.get<
int>(
"Number of Nonlinear Iterations for Step Size Restriction");
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)
588 lastPos = str.find_first_not_of(delimiters, pos);
589 pos = str.find_first_of(delimiters, lastPos);
593 std::sort(breakPoints_.begin(),breakPoints_.end());
596 currentBreakPoints_.resize(breakPoints_.size());
597 std::copy(breakPoints_.begin(),breakPoints_.end(),currentBreakPoints_.begin());
600 if ( doOutput_(Teuchos::VERB_HIGH) ) {
601 RCP<Teuchos::FancyOStream> out = this->getOStream();
602 Teuchos::OSTab ostab(out,1,
"setParameterList");
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;
692 const StepperBase<Scalar>& stepper)
694 if (stepControlState_ == UNINITIALIZED) {
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(
706 "Error, this ImplicitBDF stepper has not taken a step yet, so it "
707 "cannot take a step of order " << desiredOrder <<
" > 1!\n");
709 TEUCHOS_TEST_FOR_EXCEPT(!(desiredOrder <= currentOrder_+1));
710 currentOrder_ = desiredOrder;
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;