Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_ImplicitBDFStepperStepControl_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_STEP_CONTROL_DEF_H
30#define Rythmos_IMPLICITBDF_STEPPER_STEP_CONTROL_DEF_H
31
32// disable clang warnings
33#ifdef __clang__
34#pragma clang system_header
35#endif
36
37#include "Rythmos_ImplicitBDFStepperStepControl_decl.hpp"
38#include "Rythmos_ImplicitBDFStepper.hpp"
39#include "Rythmos_ImplicitBDFStepperErrWtVecCalc.hpp"
40
41namespace Rythmos {
42
43template<class Scalar>
44void ImplicitBDFStepperStepControl<Scalar>::setStepControlState_(StepControlStrategyState newState)
45{
46 if (stepControlState_ == UNINITIALIZED) {
47 TEUCHOS_TEST_FOR_EXCEPT(newState != BEFORE_FIRST_STEP);
48 } else if (stepControlState_ == BEFORE_FIRST_STEP) {
49 TEUCHOS_TEST_FOR_EXCEPT(newState != MID_STEP);
50 } else if (stepControlState_ == MID_STEP) {
51 TEUCHOS_TEST_FOR_EXCEPT(newState != AFTER_CORRECTION);
52 } else if (stepControlState_ == AFTER_CORRECTION) {
53 TEUCHOS_TEST_FOR_EXCEPT(newState != READY_FOR_NEXT_STEP);
54 checkReduceOrderCalled_ = false;
55 } else if (stepControlState_ == READY_FOR_NEXT_STEP) {
56 TEUCHOS_TEST_FOR_EXCEPT(newState != MID_STEP);
57 }
58 stepControlState_ = newState;
59}
60
61template<class Scalar>
62StepControlStrategyState ImplicitBDFStepperStepControl<Scalar>::getCurrentState()
63{
64 return(stepControlState_);
65}
66
67template<class Scalar>
68void ImplicitBDFStepperStepControl<Scalar>::updateCoeffs_()
69{
70 TEUCHOS_TEST_FOR_EXCEPT(!((stepControlState_ == BEFORE_FIRST_STEP) || (stepControlState_ == READY_FOR_NEXT_STEP)));
71 using Teuchos::as;
72 typedef Teuchos::ScalarTraits<Scalar> ST;
73
74 // If the number of steps taken with constant order and constant stepsize is
75 // more than the current order + 1 then we don't bother to update the
76 // coefficients because we've reached a constant step-size formula. When
77 // this is is not true, then we update the coefficients for the variable
78 // step-sizes.
79 if ((hh_ != usedStep_) || (currentOrder_ != usedOrder_)) {
80 nscsco_ = 0;
81 }
82 nscsco_ = std::min(nscsco_+1,usedOrder_+2);
83 if (currentOrder_+1 >= nscsco_) {
84 alpha_[0] = ST::one();
85 Scalar temp1 = hh_;
86 sigma_[0] = ST::one();
87 gamma_[0] = ST::zero();
88 for (int i=1;i<=currentOrder_;++i) {
89 Scalar temp2 = psi_[i-1];
90 psi_[i-1] = temp1;
91 temp1 = temp2 + hh_;
92 alpha_[i] = hh_/temp1;
93 sigma_[i] = Scalar(i+1)*sigma_[i-1]*alpha_[i];
94 gamma_[i] = gamma_[i-1]+alpha_[i-1]/hh_;
95 }
96 psi_[currentOrder_] = temp1;
97 alpha_s_ = ST::zero();
98 alpha_0_ = ST::zero();
99 for (int i=0;i<currentOrder_;++i) {
100 alpha_s_ = alpha_s_ - Scalar(ST::one()/(i+ST::one()));
101 alpha_0_ = alpha_0_ - alpha_[i];
102 }
103 cj_ = -alpha_s_/hh_;
104 ck_ = std::abs(alpha_[currentOrder_]+alpha_s_-alpha_0_);
105 ck_ = std::max(ck_,alpha_[currentOrder_]);
106 }
107 RCP<Teuchos::FancyOStream> out = this->getOStream();
108 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
109 Teuchos::OSTab ostab(out,1,"updateCoeffs_");
110 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
111 for (int i=0;i<=maxOrder_;++i) {
112 *out << "alpha_[" << i << "] = " << alpha_[i] << std::endl;
113 *out << "sigma_[" << i << "] = " << sigma_[i] << std::endl;
114 *out << "gamma_[" << i << "] = " << gamma_[i] << std::endl;
115 *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
116 *out << "alpha_s_ = " << alpha_s_ << std::endl;
117 *out << "alpha_0_ = " << alpha_0_ << std::endl;
118 *out << "cj_ = " << cj_ << std::endl;
119 *out << "ck_ = " << ck_ << std::endl;
120 }
121 }
122}
123
124template<class Scalar>
125ImplicitBDFStepperStepControl<Scalar>::ImplicitBDFStepperStepControl()
126{
127 defaultInitializeAllData_();
128}
129
130template<class Scalar>
131void ImplicitBDFStepperStepControl<Scalar>::initialize(const StepperBase<Scalar>& stepper)
132{
133 using Teuchos::as;
134 typedef Teuchos::ScalarTraits<Scalar> ST;
135 using Thyra::createMember;
136
137 RCP<Teuchos::FancyOStream> out = this->getOStream();
138 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
139 const bool doTrace = (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH));
140 Teuchos::OSTab ostab(out,1,"initialize");
141
142 if (doTrace) {
143 *out
144 << "\nEntering " << this->Teuchos::Describable::description()
145 << "::initialize()...\n";
146 }
147
148 // Set initial time:
149 TimeRange<Scalar> stepperRange = stepper.getTimeRange();
150 TEUCHOS_TEST_FOR_EXCEPTION(
151 !stepperRange.isValid(),
152 std::logic_error,
153 "Error, Stepper does not have valid time range for initialization of ImplicitBDFStepperStepControl!\n"
154 );
155 time_ = stepperRange.upper();
156
157 if (parameterList_ == Teuchos::null) {
158 RCP<Teuchos::ParameterList> emptyParameterList = Teuchos::rcp(new Teuchos::ParameterList);
159 this->setParameterList(emptyParameterList);
160 }
161
162 if (is_null(errWtVecCalc_)) {
163 RCP<ImplicitBDFStepperErrWtVecCalc<Scalar> > IBDFErrWtVecCalc = rcp(new ImplicitBDFStepperErrWtVecCalc<Scalar>());
164 errWtVecCalc_ = IBDFErrWtVecCalc;
165 }
166
167 // 08/22/07 initialize can be called from the stepper when setInitialCondition is called.
168 checkReduceOrderCalled_ = false;
169 stepControlState_ = UNINITIALIZED;
170
171 currentOrder_=1; // Current order of integration
172 oldOrder_=1; // previous order of integration
173 usedOrder_=1; // order used in current step (used after currentOrder_ is updated)
174 alpha_s_=Scalar(-ST::one()); // $\alpha_s$ fixed-leading coefficient of this BDF method
175 alpha_.clear(); // $\alpha_j(n)=h_n/\psi_j(n)$ coefficient used in local error test
176 // note: $h_n$ = current step size, n = current time step
177 gamma_.clear(); // calculate time derivative of history array for predictor
178 psi_.clear(); // $\psi_j(n) = t_n-t_{n-j}$ intermediary variable used to
179 sigma_.clear(); // $\sigma_j(n) = \frac{h_n^j(j-1)!}{\psi_1(n)*\cdots *\psi_j(n)}$
180 Scalar zero = ST::zero();
181 for (int i=0 ; i<=maxOrder_ ; ++i) {
182 alpha_.push_back(zero);
183 gamma_.push_back(zero);
184 psi_.push_back(zero);
185 sigma_.push_back(zero);
186 }
187 alpha_0_=zero; // $-\sum_{j=1}^k \alpha_j(n)$ coefficient used in local error test
188 cj_=zero; // $-\alpha_s/h_n$ coefficient used in local error test
189 ck_=zero; // local error coefficient gamma_[0] = 0; // $\gamma_j(n)=\sum_{l=1}^{j-1}1/\psi_l(n)$ coefficient used to
190 hh_=zero;
191 numberOfSteps_=0; // number of total time integration steps taken
192 nef_=0;
193 usedStep_=zero;
194 nscsco_=0;
195 Ek_=zero;
196 Ekm1_=zero;
197 Ekm2_=zero;
198 Ekp1_=zero;
199 Est_=zero;
200 Tk_=zero;
201 Tkm1_=zero;
202 Tkm2_=zero;
203 Tkp1_=zero;
204 newOrder_=currentOrder_;
205 initialPhase_=true;
206
207 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
208 *out << "currentOrder_ = " << currentOrder_ << std::endl;
209 *out << "oldOrder_ = " << oldOrder_ << std::endl;
210 *out << "usedOrder_ = " << usedOrder_ << std::endl;
211 *out << "alpha_s_ = " << alpha_s_ << std::endl;
212 for (int i=0 ; i<=maxOrder_ ; ++i) {
213 *out << "alpha_[" << i << "] = " << alpha_[i] << std::endl;
214 *out << "gamma_[" << i << "] = " << gamma_[i] << std::endl;
215 *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
216 *out << "sigma_[" << i << "] = " << sigma_[i] << std::endl;
217 }
218 *out << "alpha_0_ = " << alpha_0_ << std::endl;
219 *out << "cj_ = " << cj_ << std::endl;
220 *out << "ck_ = " << ck_ << std::endl;
221 *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
222 *out << "nef_ = " << nef_ << std::endl;
223 *out << "usedStep_ = " << usedStep_ << std::endl;
224 *out << "nscsco_ = " << nscsco_ << std::endl;
225 *out << "Ek_ = " << Ek_ << std::endl;
226 *out << "Ekm1_ = " << Ekm1_ << std::endl;
227 *out << "Ekm2_ = " << Ekm2_ << std::endl;
228 *out << "Ekp1_ = " << Ekp1_ << std::endl;
229 *out << "Est_ = " << Est_ << std::endl;
230 *out << "Tk_ = " << Tk_ << std::endl;
231 *out << "Tkm1_ = " << Tkm1_ << std::endl;
232 *out << "Tkm2_ = " << Tkm2_ << std::endl;
233 *out << "Tkp1_ = " << Tkp1_ << std::endl;
234 *out << "newOrder_ = " << newOrder_ << std::endl;
235 *out << "initialPhase_ = " << initialPhase_ << std::endl;
236 }
237
238
239 if (is_null(delta_)) {
240 delta_ = createMember(stepper.get_x_space());
241 }
242 if (is_null(errWtVec_)) {
243 errWtVec_ = createMember(stepper.get_x_space());
244 }
245 V_S(delta_.ptr(),zero);
246
247 setStepControlState_(BEFORE_FIRST_STEP);
248
249 if (doTrace) {
250 *out
251 << "\nLeaving " << this->Teuchos::Describable::description()
252 << "::initialize()...\n";
253 }
254}
255
256template<class Scalar>
257void ImplicitBDFStepperStepControl<Scalar>::getFirstTimeStep_(const StepperBase<Scalar>& stepper)
258{
259
260 TEUCHOS_TEST_FOR_EXCEPT(!(stepControlState_ == BEFORE_FIRST_STEP));
261
262 using Teuchos::as;
263 typedef Teuchos::ScalarTraits<Scalar> ST;
264 Scalar zero = ST::zero();
265
266 RCP<Teuchos::FancyOStream> out = this->getOStream();
267 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
268 const bool doTrace = (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH));
269 Teuchos::OSTab ostab(out,1,"getFirstTimeStep_");
270
271 const ImplicitBDFStepper<Scalar>& implicitBDFStepper
272 = Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
273 const Thyra::VectorBase<Scalar>& xHistory0
274 = implicitBDFStepper.getxHistory(0);
275 errWtVecCalc_->errWtVecSet(&*errWtVec_,xHistory0,relErrTol_,absErrTol_);
276
277 // Choose initial step-size
278 Scalar time_to_stop = stopTime_ - time_;
279 Scalar currentTimeStep = ST::nan();
280 if (stepSizeType_ == STEP_TYPE_FIXED) {
281 currentTimeStep = hh_;
282 //currentTimeStep = 0.1 * time_to_stop;
283 //currentTimeStep = std::min(hh_, currentTimeStep);
284 } else {
285 // compute an initial step-size based on rate of change
286 // in the solution initially
287 const Thyra::VectorBase<Scalar>& xHistory1
288 = implicitBDFStepper.getxHistory(1);
289 Scalar ypnorm = wRMSNorm_(*errWtVec_,xHistory1);
290 if (ypnorm > zero) { // time-dependent DAE
291 currentTimeStep = std::min( h0_max_factor_*std::abs(time_to_stop),
292 std::sqrt(2.0)/(h0_safety_*ypnorm) );
293 } else { // non-time-dependent DAE
294 currentTimeStep = h0_max_factor_*std::abs(time_to_stop);
295 }
296 // choose std::min of user specified value and our value:
297 if (hh_ > zero) {
298 currentTimeStep = std::min(hh_, currentTimeStep);
299 }
300 // check for maximum step-size:
301#ifdef HAVE_RYTHMOS_DEBUG
302 TEUCHOS_TEST_FOR_EXCEPT(ST::isnaninf(currentTimeStep));
303#endif // HAVE_RYTHMOS_DEBUG
304 Scalar rh = std::abs(currentTimeStep)*h_max_inv_;
305 if (rh>1.0) {
306 currentTimeStep = currentTimeStep/rh;
307 }
308 }
309 hh_ = currentTimeStep;
310
311 // Coefficient initialization
312 psi_[0] = hh_;
313 cj_ = 1/psi_[0];
314
315 if (doTrace) {
316 *out << "\nhh_ = " << hh_ << std::endl;
317 *out << "psi_[0] = " << psi_[0] << std::endl;
318 *out << "cj_ = " << cj_ << std::endl;
319 }
320
321}
322
323template<class Scalar>
324void ImplicitBDFStepperStepControl<Scalar>::setRequestedStepSize(
325 const StepperBase<Scalar>& stepper
326 ,const Scalar& stepSize
327 ,const StepSizeType& stepSizeType
328 )
329{
330 typedef Teuchos::ScalarTraits<Scalar> ST;
331 TEUCHOS_TEST_FOR_EXCEPT(!((stepControlState_ == UNINITIALIZED) ||
332 (stepControlState_ == BEFORE_FIRST_STEP) ||
333 (stepControlState_ == READY_FOR_NEXT_STEP) ||
334 (stepControlState_ == MID_STEP)));
335 TEUCHOS_TEST_FOR_EXCEPTION(
336 ((stepSizeType == STEP_TYPE_FIXED) && (stepSize == ST::zero())),
337 std::logic_error,
338 "Error, step size type == STEP_TYPE_FIXED, but requested step size == 0!\n"
339 );
340 bool didInitialization = false;
341 if (stepControlState_ == UNINITIALIZED) {
342 initialize(stepper);
343 didInitialization = true;
344 }
345
346 // errWtVecSet_ is called during initialize
347 if (!didInitialization) {
348 const ImplicitBDFStepper<Scalar>& implicitBDFStepper = Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
349 const Thyra::VectorBase<Scalar>& xHistory = implicitBDFStepper.getxHistory(0);
350 errWtVecCalc_->errWtVecSet(&*errWtVec_,xHistory,relErrTol_,absErrTol_);
351 }
352
353 stepSizeType_ = stepSizeType;
354 if (stepSizeType_ == STEP_TYPE_FIXED) {
355 hh_ = stepSize;
356 if (numberOfSteps_ == 0) {
357 psi_[0] = hh_;
358 cj_ = 1/psi_[0];
359 }
360 } else { // STEP_TYPE_VARIABLE
361 if (stepSize != ST::zero()) {
362 maxTimeStep_ = stepSize;
363 h_max_inv_ = Scalar(ST::one()/maxTimeStep_);
364 }
365 }
366}
367
368template<class Scalar>
369void ImplicitBDFStepperStepControl<Scalar>::nextStepSize(const StepperBase<Scalar>& stepper, Scalar* stepSize, StepSizeType* stepSizeType, int* order)
370{
371 TEUCHOS_TEST_FOR_EXCEPT(!((stepControlState_ == BEFORE_FIRST_STEP) ||
372 (stepControlState_ == MID_STEP) ||
373 (stepControlState_ == READY_FOR_NEXT_STEP) )
374 );
375 if (stepControlState_ == BEFORE_FIRST_STEP) {
376 getFirstTimeStep_(stepper);
377 }
378 if (stepControlState_ != MID_STEP) {
379 this->updateCoeffs_();
380 }
381 if (stepSizeType_ == STEP_TYPE_VARIABLE) {
382 if (hh_ > maxTimeStep_) {
383 hh_ = maxTimeStep_;
384 }
385 }
386 *stepSize = hh_;
387 *stepSizeType = stepSizeType_;
388 *order = currentOrder_;
389 if (stepControlState_ != MID_STEP) {
390 setStepControlState_(MID_STEP);
391 }
392}
393
394template<class Scalar>
395void ImplicitBDFStepperStepControl<Scalar>::setCorrection(
396 const StepperBase<Scalar>& /* stepper */
397 ,const RCP<const Thyra::VectorBase<Scalar> >& /* soln */
398 ,const RCP<const Thyra::VectorBase<Scalar> >& ee
399 ,int solveStatus)
400{
401 TEUCHOS_TEST_FOR_EXCEPT(stepControlState_ != MID_STEP);
402 TEUCHOS_TEST_FOR_EXCEPTION(is_null(ee), std::logic_error, "Error, ee == Teuchos::null!\n");
403 ee_ = ee;
404 newtonConvergenceStatus_ = solveStatus;
405 setStepControlState_(AFTER_CORRECTION);
406}
407
408template<class Scalar>
409void ImplicitBDFStepperStepControl<Scalar>::completeStep(const StepperBase<Scalar>& stepper)
410{
411 TEUCHOS_TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
412 using Teuchos::as;
413 typedef Teuchos::ScalarTraits<Scalar> ST;
414
415 numberOfSteps_ ++;
416 nef_ = 0;
417 time_ += hh_;
418 RCP<Teuchos::FancyOStream> out = this->getOStream();
419 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
420 Teuchos::OSTab ostab(out,1,"completeStep_");
421
422 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
423 *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
424 *out << "nef_ = " << nef_ << std::endl;
425 *out << "time_ = " << time_ << std::endl;
426 }
427
428 // Only update the time step if we are NOT running constant stepsize.
429 bool adjustStep = (stepSizeType_ == STEP_TYPE_VARIABLE);
430
431 Scalar newTimeStep = hh_;
432 Scalar rr = ST::one(); // step size ratio = new step / old step
433 // 03/11/04 tscoffe: Here is the block for choosing order & step-size when
434 // the local error test PASSES (and Newton succeeded).
435 int orderDiff = currentOrder_ - usedOrder_;
436 usedOrder_ = currentOrder_;
437 usedStep_ = hh_;
438 if ((newOrder_ == currentOrder_-1) || (currentOrder_ == maxOrder_)) {
439 // If we reduced our order or reached std::max order then move to the next phase
440 // of integration where we don't automatically double the step-size and
441 // increase the order.
442 initialPhase_ = false;
443 }
444 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
445 *out << "initialPhase_ = " << initialPhase_ << std::endl;
446 }
447 if (initialPhase_) {
448 currentOrder_++;
449 newTimeStep = h_phase0_incr_ * hh_;
450 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
451 *out << "currentOrder_ = " << currentOrder_ << std::endl;
452 *out << "newTimeStep = " << newTimeStep << std::endl;
453 }
454 } else { // not in the initial phase of integration
455 BDFactionFlag action = ACTION_UNSET;
456 if (newOrder_ == currentOrder_-1) {
457 action = ACTION_LOWER;
458 } else if (newOrder_ == maxOrder_) {
459 action = ACTION_MAINTAIN;
460 } else if ((currentOrder_+1>=nscsco_) || (orderDiff == 1)) {
461 // If we just raised the order last time then we won't raise it again
462 // until we've taken currentOrder_+1 steps at order currentOrder_.
463 action = ACTION_MAINTAIN;
464 } else { // consider changing the order
465 const ImplicitBDFStepper<Scalar>& implicitBDFStepper = Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
466 const Thyra::VectorBase<Scalar>& xHistory = implicitBDFStepper.getxHistory(currentOrder_+1);
467 V_StVpStV(delta_.ptr(),ST::one(),*ee_,Scalar(-ST::one()),xHistory);
468 Tkp1_ = wRMSNorm_(*errWtVec_,*delta_);
469 Ekp1_ = Tkp1_/(currentOrder_+2);
470 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
471 *out << "delta_ = " << std::endl;
472 delta_->describe(*out,verbLevel);
473 *out << "Tkp1_ = ||delta_||_WRMS = " << Tkp1_ << std::endl;
474 *out << "Ekp1_ = " << Ekp1_ << std::endl;
475 }
476 if (currentOrder_ == 1) {
477 if (Tkp1_ >= Tkp1_Tk_safety_ * Tk_) {
478 action = ACTION_MAINTAIN;
479 } else {
480 action = ACTION_RAISE;
481 }
482 } else {
483 if (Tkm1_ <= std::min(Tk_,Tkp1_)) {
484 action = ACTION_LOWER;
485 } else if (Tkp1_ >= Tk_) {
486 action = ACTION_MAINTAIN;
487 } else {
488 action = ACTION_RAISE;
489 }
490 }
491 }
492 if (currentOrder_ < minOrder_) {
493 action = ACTION_RAISE;
494 } else if ( (currentOrder_ == minOrder_) && (action == ACTION_LOWER) ) {
495 action = ACTION_MAINTAIN;
496 }
497 if (action == ACTION_RAISE) {
498 currentOrder_++;
499 Est_ = Ekp1_;
500 } else if (action == ACTION_LOWER) {
501 currentOrder_--;
502 Est_ = Ekm1_;
503 }
504 newTimeStep = hh_;
505 rr = pow(r_safety_*Est_+r_fudge_,-1.0/(currentOrder_+1.0));
506 if (rr >= r_hincr_test_) {
507 rr = r_hincr_;
508 newTimeStep = rr*hh_;
509 } else if (rr <= 1) {
510 rr = std::max(r_min_,std::min(r_max_,rr));
511 newTimeStep = rr*hh_;
512 }
513 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
514 *out << "Est_ = " << Est_ << std::endl;
515 *out << "rr = " << rr << std::endl;
516 *out << "newTimeStep = " << newTimeStep << std::endl;
517 }
518 }
519
520 if (time_ < stopTime_) {
521 // If the step needs to be adjusted:
522 if (adjustStep) {
523 newTimeStep = std::max(newTimeStep, minTimeStep_);
524 newTimeStep = std::min(newTimeStep, maxTimeStep_);
525
526 Scalar nextTimePt = time_ + newTimeStep;
527
528 if (nextTimePt > stopTime_) {
529 nextTimePt = stopTime_;
530 newTimeStep = stopTime_ - time_;
531 }
532
533 hh_ = newTimeStep;
534
535 } else { // if time step is constant for this step:
536 Scalar nextTimePt = time_ + hh_;
537
538 if (nextTimePt > stopTime_) {
539 nextTimePt = stopTime_;
540 hh_ = stopTime_ - time_;
541 }
542 }
543 }
544 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
545 *out << "hh_ = " << hh_ << std::endl;
546 *out << "currentOrder_ = " << currentOrder_ << std::endl;
547 }
548 setStepControlState_(READY_FOR_NEXT_STEP);
549}
550
551template<class Scalar>
552AttemptedStepStatusFlag ImplicitBDFStepperStepControl<Scalar>::rejectStep(const StepperBase<Scalar>& /* stepper */)
553{
554 TEUCHOS_TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
555
556 using Teuchos::as;
557
558 // This routine puts its output in newTimeStep and newOrder
559
560 // This routine changes the following variables:
561 // initialPhase, nef, psi, newTimeStep,
562 // newOrder, currentOrder_, currentTimeStep, dsDae.xHistory,
563 // dsDae.qHistory, nextTimePt,
564 // currentTimeStepSum, nextTimePt
565
566 // This routine reads but does not change the following variables:
567 // r_factor, r_safety, Est_, r_fudge_, r_min_, r_max_,
568 // minTimeStep_, maxTimeStep_, time, stopTime_
569
570 // Only update the time step if we are NOT running constant stepsize.
571 bool adjustStep = (stepSizeType_ == STEP_TYPE_VARIABLE);
572
573 Scalar newTimeStep = hh_;
574 Scalar rr = 1.0; // step size ratio = new step / old step
575 nef_++;
576 RCP<Teuchos::FancyOStream> out = this->getOStream();
577 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
578 Teuchos::OSTab ostab(out,1,"rejectStep_");
579 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
580 *out << "adjustStep = " << adjustStep << std::endl;
581 *out << "nef_ = " << nef_ << std::endl;
582 }
583 if (nef_ >= max_LET_fail_) {
584 TEUCHOS_TEST_FOR_EXCEPTION(nef_ >= max_LET_fail_, std::logic_error, "Error, maximum number of local error test failures.\n");
585 }
586 initialPhase_ = false;
587 if (adjustStep) {
588 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
589 *out << "initialPhase_ = " << initialPhase_ << std::endl;
590 }
591 for (int i=1;i<=currentOrder_;++i) {
592 psi_[i-1] = psi_[i] - hh_;
593 }
594
595 if ((newtonConvergenceStatus_ < 0)) {
597 // rely on the error estimate - it may be full of Nan's.
598 rr = r_min_;
599 newTimeStep = rr * hh_;
600
601 if (nef_ > 2) {
602 newOrder_ = 1;//consistent with block below.
603 }
604 } else {
605 // 03/11/04 tscoffe: Here is the block for choosing order &
606 // step-size when the local error test FAILS (but Newton
607 // succeeded).
608 if (nef_ == 1) { // first local error test failure
609 rr = r_factor_*pow(r_safety_*Est_+r_fudge_,-1.0/(newOrder_+1.0));
610 rr = std::max(r_min_,std::min(r_max_,rr));
611 newTimeStep = rr * hh_;
612 } else if (nef_ == 2) { // second failure
613 rr = r_min_;
614 newTimeStep = rr * hh_;
615 } else if (nef_ > 2) { // third and later failures
616 newOrder_ = 1;
617 rr = r_min_;
618 newTimeStep = rr * hh_;
619 }
620 }
621 if (newOrder_ >= minOrder_) {
622 currentOrder_ = newOrder_;
623 }
624 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
625 *out << "newTimeStep = " << newTimeStep << std::endl;
626 *out << "rr = " << rr << std::endl;
627 *out << "newOrder_ = " << newOrder_ << std::endl;
628 *out << "currentOrder_ = " << currentOrder_ << std::endl;
629 }
630 if (numberOfSteps_ == 0) { // still first step
631 psi_[0] = newTimeStep;
632 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
633 *out << "numberOfSteps_ == 0:" << std::endl;
634 *out << "psi_[0] = " << psi_[0] << std::endl;
635 }
636 }
637 } else if (!adjustStep) {
638 if ( as<int>(verbLevel) != as<int>(Teuchos::VERB_NONE) ) {
639 *out << "Rythmos_ImplicitBDFStepperStepControl::rejectStep(...): "
640 << "Warning: Local error test failed with constant step-size."
641 << std::endl;
642 }
643 }
644
645 AttemptedStepStatusFlag return_status = PREDICT_AGAIN;
646
647 // If the step needs to be adjusted:
648 if (adjustStep) {
649 newTimeStep = std::max(newTimeStep, minTimeStep_);
650 newTimeStep = std::min(newTimeStep, maxTimeStep_);
651
652 Scalar nextTimePt = time_ + newTimeStep;
653
654 if (nextTimePt > stopTime_) {
655 nextTimePt = stopTime_;
656 newTimeStep = stopTime_ - time_;
657 }
658
659 hh_ = newTimeStep;
660
661 } else { // if time step is constant for this step:
662 Scalar nextTimePt = time_ + hh_;
663
664 if (nextTimePt > stopTime_) {
665 nextTimePt = stopTime_;
666 hh_ = stopTime_ - time_;
667 }
668 return_status = CONTINUE_ANYWAY;
669 }
670 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
671 *out << "hh_ = " << hh_ << std::endl;
672 }
673
674 if (return_status == PREDICT_AGAIN) {
675 setStepControlState_(READY_FOR_NEXT_STEP);
676 } else if (return_status == CONTINUE_ANYWAY) {
677 // do nothing, as we'll call completeStep which must be in AFTER_CORRECTION state.
678 }
679
680 return(return_status);
681}
682
683template<class Scalar>
684Scalar ImplicitBDFStepperStepControl<Scalar>::checkReduceOrder_(
685 const StepperBase<Scalar>& stepper)
686{
687 TEUCHOS_TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
688 TEUCHOS_TEST_FOR_EXCEPT(checkReduceOrderCalled_ == true);
689
690 using Teuchos::as;
691
692 const ImplicitBDFStepper<Scalar>& implicitBDFStepper =
693 Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
694
695 // This routine puts its output in newOrder_
696
697 RCP<Teuchos::FancyOStream> out = this->getOStream();
698 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
699 Teuchos::OSTab ostab(out,1,"checkReduceOrder_");
700 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
701 *out << "sigma_[" << currentOrder_ << "] = "
702 << sigma_[currentOrder_] << std::endl;
703 *out << "ee_ = " << std::endl;
704 ee_->describe(*out,verbLevel);
705 *out << "errWtVec_ = " << std::endl;
706 errWtVec_->describe(*out,verbLevel);
707 }
708
709 Scalar enorm = wRMSNorm_(*errWtVec_,*ee_);
710 Ek_ = sigma_[currentOrder_]*enorm;
711 Tk_ = Scalar(currentOrder_+1)*Ek_;
712 Est_ = Ek_;
713 newOrder_ = currentOrder_;
714 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
715 *out << "currentOrder_ = " << currentOrder_ << std::endl;
716 *out << "Ek_ = " << Ek_ << std::endl;
717 *out << "Tk_ = " << Tk_ << std::endl;
718 *out << "enorm = " << enorm << std::endl;
719 }
720 if (currentOrder_>1) {
721 const Thyra::VectorBase<Scalar>& xHistoryCur =
722 implicitBDFStepper.getxHistory(currentOrder_);
723 V_VpV(delta_.ptr(),xHistoryCur,*ee_);
724 Ekm1_ = sigma_[currentOrder_-1]*wRMSNorm_(*errWtVec_,*delta_);
725 Tkm1_ = currentOrder_*Ekm1_;
726 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
727 *out << "Ekm1_ = " << Ekm1_ << std::endl;
728 *out << "Tkm1_ = " << Tkm1_ << std::endl;
729 }
730 if (currentOrder_>2) {
731 const Thyra::VectorBase<Scalar>& xHistoryPrev =
732 implicitBDFStepper.getxHistory(currentOrder_-1);
733 Vp_V(delta_.ptr(),xHistoryPrev);
734 Ekm2_ = sigma_[currentOrder_-2]*wRMSNorm_(*errWtVec_,*delta_);
735 Tkm2_ = (currentOrder_-1)*Ekm2_;
736 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
737 *out << "Ekm2_ = " << Ekm2_ << std::endl;
738 *out << "Tkm2_ = " << Tkm2_ << std::endl;
739 }
740 if (std::max(Tkm1_,Tkm2_)<=Tk_) {
741 newOrder_--;
742 Est_ = Ekm1_;
743 }
744 }
745 else if (Tkm1_ <= Tkm1_Tk_safety_ * Tk_) {
746 newOrder_--;
747 Est_ = Ekm1_;
748 }
749 }
750 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
751 *out << "Est_ = " << Est_ << std::endl;
752 *out << "newOrder_= " << newOrder_ << std::endl;
753 }
754 checkReduceOrderCalled_ = true;
755 return(enorm);
756}
757
758template<class Scalar>
759bool ImplicitBDFStepperStepControl<Scalar>::acceptStep(const StepperBase<Scalar>& stepper, Scalar* LETValue)
760{
761 TEUCHOS_TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
762 typedef Teuchos::ScalarTraits<Scalar> ST;
763 RCP<Teuchos::FancyOStream> out = this->getOStream();
764 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
765 Teuchos::OSTab ostab(out,1,"acceptStep");
766
767 bool return_status = false;
768 Scalar enorm = checkReduceOrder_(stepper);
769 Scalar LET = ck_*enorm;
770
771 if (failStepIfNonlinearSolveFails_ && (newtonConvergenceStatus_ < 0) )
772 return false;
773
774 if (LETValue) {
775 *LETValue = LET;
776 }
777 if (LET < ST::one()) {
778 return_status = true;
779 }
780 if ( Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
781 *out << "ck_ = " << ck_ << std::endl;
782 *out << "enorm = " << enorm << std::endl;
783 *out << "Local Truncation Error Check: (ck*enorm) < 1: (" << LET << ") <?= 1" << std::endl;
784 }
785 return(return_status);
786}
787
788template<class Scalar>
789void ImplicitBDFStepperStepControl<Scalar>::describe(
790 Teuchos::FancyOStream &out,
791 const Teuchos::EVerbosityLevel verbLevel
792 ) const
793{
794
795 using Teuchos::as;
796
797 if ( (as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT) ) ||
798 (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
799 ) {
800 out << this->description() << "::describe" << std::endl;
801 }
802 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
803 out << "time_ = " << time_ << std::endl;
804 out << "hh_ = " << hh_ << std::endl;
805 out << "currentOrder_ = " << currentOrder_ << std::endl;
806 }
807 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
808 }
809 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
810 out << "ee_ = ";
811 if (ee_ == Teuchos::null) {
812 out << "Teuchos::null" << std::endl;
813 } else {
814 ee_->describe(out,verbLevel);
815 }
816 out << "delta_ = ";
817 if (delta_ == Teuchos::null) {
818 out << "Teuchos::null" << std::endl;
819 } else {
820 delta_->describe(out,verbLevel);
821 }
822 out << "errWtVec_ = ";
823 if (errWtVec_ == Teuchos::null) {
824 out << "Teuchos::null" << std::endl;
825 } else {
826 errWtVec_->describe(out,verbLevel);
827 }
828 }
829}
830
831template<class Scalar>
832void ImplicitBDFStepperStepControl<Scalar>::setParameterList(
833 RCP<Teuchos::ParameterList> const& paramList)
834{
835
836 using Teuchos::as;
837 // typedef Teuchos::ScalarTraits<Scalar> ST; // unused
838
839 TEUCHOS_TEST_FOR_EXCEPT(paramList == Teuchos::null);
840 paramList->validateParameters(*this->getValidParameters(),0);
841 parameterList_ = paramList;
842 Teuchos::readVerboseObjectSublist(&*parameterList_,this);
843
844 minOrder_ = parameterList_->get("minOrder",int(1)); // minimum order
845 TEUCHOS_TEST_FOR_EXCEPTION(
846 !((1 <= minOrder_) && (minOrder_ <= 5)), std::logic_error,
847 "Error, minOrder_ = " << minOrder_ << " is not in range [1,5]!\n"
848 );
849 maxOrder_ = parameterList_->get("maxOrder",int(5)); // maximum order
850 TEUCHOS_TEST_FOR_EXCEPTION(
851 !((1 <= maxOrder_) && (maxOrder_ <= 5)), std::logic_error,
852 "Error, maxOrder_ = " << maxOrder_ << " is not in range [1,5]!\n"
853 );
854
855 relErrTol_ = parameterList_->get( "relErrTol", Scalar(1.0e-4) );
856 absErrTol_ = parameterList_->get( "absErrTol", Scalar(1.0e-6) );
857 bool constantStepSize = parameterList_->get( "constantStepSize", false );
858 stopTime_ = parameterList_->get( "stopTime", Scalar(1.0) );
859
860 if (constantStepSize == true) {
861 stepSizeType_ = STEP_TYPE_FIXED;
862 } else {
863 stepSizeType_ = STEP_TYPE_VARIABLE;
864 }
865
866 failStepIfNonlinearSolveFails_ =
867 parameterList_->get( "failStepIfNonlinearSolveFails", false );
868
869 RCP<Teuchos::FancyOStream> out = this->getOStream();
870 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
871 Teuchos::OSTab ostab(out,1,"setParameterList");
872 out->precision(15);
873
874 setDefaultMagicNumbers_(parameterList_->sublist("magicNumbers"));
875
876 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
877 *out << "minOrder_ = " << minOrder_ << std::endl;
878 *out << "maxOrder_ = " << maxOrder_ << std::endl;
879 *out << "relErrTol = " << relErrTol_ << std::endl;
880 *out << "absErrTol = " << absErrTol_ << std::endl;
881 *out << "stepSizeType = " << stepSizeType_ << std::endl;
882 *out << "stopTime_ = " << stopTime_ << std::endl;
883 *out << "failStepIfNonlinearSolveFails_ = "
884 << failStepIfNonlinearSolveFails_ << std::endl;
885 }
886
887}
888
889template<class Scalar>
890void ImplicitBDFStepperStepControl<Scalar>::setDefaultMagicNumbers_(
891 Teuchos::ParameterList &magicNumberList)
892{
893
894 using Teuchos::as;
895
896 // Magic Number Defaults:
897 h0_safety_ = magicNumberList.get( "h0_safety", Scalar(2.0) );
898 h0_max_factor_ = magicNumberList.get( "h0_max_factor", Scalar(0.001) );
899 h_phase0_incr_ = magicNumberList.get( "h_phase0_incr", Scalar(2.0) );
900 h_max_inv_ = magicNumberList.get( "h_max_inv", Scalar(0.0) );
901 Tkm1_Tk_safety_ = magicNumberList.get( "Tkm1_Tk_safety", Scalar(2.0) );
902 Tkp1_Tk_safety_ = magicNumberList.get( "Tkp1_Tk_safety", Scalar(0.5) );
903 r_factor_ = magicNumberList.get( "r_factor", Scalar(0.9) );
904 r_safety_ = magicNumberList.get( "r_safety", Scalar(2.0) );
905 r_fudge_ = magicNumberList.get( "r_fudge", Scalar(0.0001) );
906 r_min_ = magicNumberList.get( "r_min", Scalar(0.125) );
907 r_max_ = magicNumberList.get( "r_max", Scalar(0.9) );
908 r_hincr_test_ = magicNumberList.get( "r_hincr_test", Scalar(2.0) );
909 r_hincr_ = magicNumberList.get( "r_hincr", Scalar(2.0) );
910 max_LET_fail_ = magicNumberList.get( "max_LET_fail", int(15) );
911 minTimeStep_ = magicNumberList.get( "minTimeStep", Scalar(0.0) );
912 maxTimeStep_ = magicNumberList.get( "maxTimeStep", Scalar(10.0) );
913
914 RCP<Teuchos::FancyOStream> out = this->getOStream();
915 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
916 Teuchos::OSTab ostab(out,1,"setDefaultMagicNumbers_");
917 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
918 *out << "h0_safety_ = " << h0_safety_ << std::endl;
919 *out << "h0_max_factor_ = " << h0_max_factor_ << std::endl;
920 *out << "h_phase0_incr_ = " << h_phase0_incr_ << std::endl;
921 *out << "h_max_inv_ = " << h_max_inv_ << std::endl;
922 *out << "Tkm1_Tk_safety_ = " << Tkm1_Tk_safety_ << std::endl;
923 *out << "Tkp1_Tk_safety_ = " << Tkp1_Tk_safety_ << std::endl;
924 *out << "r_factor_ = " << r_factor_ << std::endl;
925 *out << "r_safety_ = " << r_safety_ << std::endl;
926 *out << "r_fudge_ = " << r_fudge_ << std::endl;
927 *out << "r_min_ = " << r_min_ << std::endl;
928 *out << "r_max_ = " << r_max_ << std::endl;
929 *out << "r_hincr_test_ = " << r_hincr_test_ << std::endl;
930 *out << "r_hincr_ = " << r_hincr_ << std::endl;
931 *out << "max_LET_fail_ = " << max_LET_fail_ << std::endl;
932 *out << "minTimeStep_ = " << minTimeStep_ << std::endl;
933 *out << "maxTimeStep_ = " << maxTimeStep_ << std::endl;
934 }
935
936}
937
938template<class Scalar>
939RCP<const Teuchos::ParameterList>
940ImplicitBDFStepperStepControl<Scalar>::getValidParameters() const
941{
942
943 static RCP<Teuchos::ParameterList> validPL;
944
945 if (is_null(validPL)) {
946
947 RCP<Teuchos::ParameterList>
948 pl = Teuchos::parameterList();
949
950 pl->set<int> ( "minOrder", 1,
951 "lower limit of order selection, guaranteed");
952 pl->set<int> ( "maxOrder", 5,
953 "upper limit of order selection, does not guarantee this order");
954 pl->set<Scalar>( "relErrTol", Scalar(1.0e-4),
955 "Relative tolerance value used in WRMS calculation.");
956 pl->set<Scalar>( "absErrTol", Scalar(1.0e-6),
957 "Absolute tolerance value used in WRMS calculation.");
958 pl->set<bool> ( "constantStepSize", false,
959 "Use constant (=0) or variable (=1) step sizes during BDF steps.");
960 pl->set<Scalar>( "stopTime", Scalar(10.0),
961 "The final time for time integration. This may be in conflict with "
962 "the Integrator's final time.");
963 pl->set<bool>("failStepIfNonlinearSolveFails", false,
964 "Power user command. Will force the function acceptStep() to return "
965 "false even if the LET is acceptable. Used to run with loose "
966 "tolerances but enforce a correct nonlinear solution to the step.");
967
968 Teuchos::ParameterList
969 &magicNumberList = pl->sublist("magicNumbers", false,
970 "These are knobs in the algorithm that have been set to reasonable "
971 "values using lots of testing and heuristics and some theory.");
972 magicNumberList.set<Scalar>( "h0_safety", Scalar(2.0),
973 "Safety factor on the initial step size for time-dependent DAEs. "
974 "Larger values will tend to reduce the initial step size.");
975 magicNumberList.set<Scalar>( "h0_max_factor", Scalar(0.001),
976 "Factor multipler for the initial step size for non-time-dependent "
977 "DAEs.");
978 magicNumberList.set<Scalar>( "h_phase0_incr", Scalar(2.0),
979 "Initial ramp-up in variable mode (stepSize multiplier).");
980 magicNumberList.set<Scalar>( "h_max_inv", Scalar(0.0),
981 "The inverse of the maximum time step (maxTimeStep). (Not used?)");
982 magicNumberList.set<Scalar>( "Tkm1_Tk_safety", Scalar(2.0),
983 "Used to help control the decrement of the order when the order is "
984 "less than or equal to 2. Larger values will tend to reduce the "
985 "order from one step to the next.");
986 magicNumberList.set<Scalar>( "Tkp1_Tk_safety", Scalar(0.5),
987 "Used to control the increment of the order when the order is one. "
988 "Larger values tend to increase the order for the next step.");
989 magicNumberList.set<Scalar>( "r_factor", Scalar(0.9),
990 "Used in rejectStep: time step ratio multiplier.");
991 magicNumberList.set<Scalar>( "r_safety", Scalar(2.0),
992 "Local error multiplier as part of time step ratio calculation.");
993 magicNumberList.set<Scalar>( "r_fudge", Scalar(0.0001),
994 "Local error addition as part of time step ratio calculation.");
995 magicNumberList.set<Scalar>( "r_min", Scalar(0.125),
996 "Used in rejectStep: How much to cut step and lower bound for time "
997 "step ratio.");
998 magicNumberList.set<Scalar>( "r_max", Scalar(0.9),
999 "Upper bound for time step ratio.");
1000 magicNumberList.set<Scalar>( "r_hincr_test", Scalar(2.0),
1001 "Used in completeStep: If time step ratio > this then set time step "
1002 "ratio to r_hincr.");
1003 magicNumberList.set<Scalar>( "r_hincr", Scalar(2.0),
1004 "Used in completeStep: Limit on time step ratio increases, not "
1005 "checked by r_max.");
1006 magicNumberList.set<int> ( "max_LET_fail", int(15),
1007 "Max number of rejected steps");
1008 magicNumberList.set<Scalar>( "minTimeStep", Scalar(0.0),
1009 "Bound on smallest time step in variable mode.");
1010 magicNumberList.set<Scalar>( "maxTimeStep", Scalar(10.0),
1011 "bound on largest time step in variable mode.");
1012
1013 Teuchos::setupVerboseObjectSublist(&*pl);
1014
1015 validPL = pl;
1016
1017 }
1018
1019 return (validPL);
1020
1021}
1022
1023template<class Scalar>
1024RCP<Teuchos::ParameterList>
1025ImplicitBDFStepperStepControl<Scalar>::unsetParameterList()
1026{
1027 RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
1028 parameterList_ = Teuchos::null;
1029 return(temp_param_list);
1030}
1031
1032template<class Scalar>
1033RCP<Teuchos::ParameterList>
1034ImplicitBDFStepperStepControl<Scalar>::getNonconstParameterList()
1035{
1036 return(parameterList_);
1037}
1038
1039template<class Scalar>
1040void ImplicitBDFStepperStepControl<Scalar>::setStepControlData(const StepperBase<Scalar>& stepper)
1041{
1042 if (stepControlState_ == UNINITIALIZED) {
1043 initialize(stepper);
1044 }
1045 const ImplicitBDFStepper<Scalar>& bdfstepper = Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
1046 int desiredOrder = bdfstepper.getOrder();
1047 TEUCHOS_TEST_FOR_EXCEPT(!((1 <= desiredOrder) && (desiredOrder <= maxOrder_)));
1048 if (stepControlState_ == BEFORE_FIRST_STEP) {
1049 TEUCHOS_TEST_FOR_EXCEPTION(
1050 desiredOrder > 1,
1051 std::logic_error,
1052 "Error, this ImplicitBDF stepper has not taken a step yet, so it cannot take a step of order " << desiredOrder << " > 1!\n"
1053 );
1054 }
1055 TEUCHOS_TEST_FOR_EXCEPT(!(desiredOrder <= currentOrder_+1));
1056 currentOrder_ = desiredOrder;
1057
1058 using Teuchos::as;
1059 RCP<Teuchos::FancyOStream> out = this->getOStream();
1060 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
1061 Teuchos::OSTab ostab(out,1,"setStepControlData");
1062
1063 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
1064 *out << "currentOrder_ = " << currentOrder_ << std::endl;
1065 }
1066}
1067
1068template<class Scalar>
1069bool ImplicitBDFStepperStepControl<Scalar>::supportsCloning() const
1070{
1071 return true;
1072}
1073
1074
1075template<class Scalar>
1076RCP<StepControlStrategyBase<Scalar> >
1077ImplicitBDFStepperStepControl<Scalar>::cloneStepControlStrategyAlgorithm() const
1078{
1079
1080 RCP<ImplicitBDFStepperStepControl<Scalar> > stepControl = rcp(new ImplicitBDFStepperStepControl<Scalar>());
1081
1082 if (!is_null(parameterList_)) {
1083 stepControl->setParameterList(parameterList_);
1084 }
1085
1086 return stepControl;
1087}
1088
1089template<class Scalar>
1090void ImplicitBDFStepperStepControl<Scalar>::setErrWtVecCalc(const RCP<ErrWtVecCalcBase<Scalar> >& errWtVecCalc)
1091{
1092 TEUCHOS_TEST_FOR_EXCEPT(is_null(errWtVecCalc));
1093 errWtVecCalc_ = errWtVecCalc;
1094}
1095
1096template<class Scalar>
1097RCP<const ErrWtVecCalcBase<Scalar> > ImplicitBDFStepperStepControl<Scalar>::getErrWtVecCalc() const
1098{
1099 return(errWtVecCalc_);
1100}
1101
1102template<class Scalar>
1103Scalar ImplicitBDFStepperStepControl<Scalar>::wRMSNorm_(
1104 const Thyra::VectorBase<Scalar>& weight,
1105 const Thyra::VectorBase<Scalar>& vector) const
1106{
1107 return(norm_2(weight,vector));
1108}
1109
1110template<class Scalar>
1111void ImplicitBDFStepperStepControl<Scalar>::defaultInitializeAllData_()
1112{
1113 typedef Teuchos::ScalarTraits<Scalar> ST;
1114 Scalar zero = ST::zero();
1115 Scalar mone = Scalar(-ST::one());
1116
1117 stepControlState_ = UNINITIALIZED;
1118 hh_ = zero;
1119 numberOfSteps_ = 0;
1120 stepSizeType_ = STEP_TYPE_VARIABLE;
1121
1122 minOrder_ = -1;
1123 maxOrder_ = -1;
1124 nef_ = 0;
1125 midStep_ = false;
1126 checkReduceOrderCalled_ = false;
1127 time_ = -std::numeric_limits<Scalar>::max();
1128 relErrTol_ = mone;
1129 absErrTol_ = mone;
1130 usedStep_ = mone;
1131 currentOrder_ = 1;
1132 usedOrder_ = -1;
1133 nscsco_ = -1;
1134 alpha_s_ = mone;
1135 alpha_0_ = mone;
1136 cj_ = mone;
1137 ck_ = mone;
1138 ck_enorm_ = mone;
1139 constantStepSize_ = false;
1140 Ek_ = mone;
1141 Ekm1_ = mone;
1142 Ekm2_ = mone;
1143 Ekp1_ = mone;
1144 Est_ = mone;
1145 Tk_ = mone;
1146 Tkm1_ = mone;
1147 Tkm2_ = mone;
1148 Tkp1_ = mone;
1149 newOrder_ = -1;
1150 oldOrder_ = -1;
1151 initialPhase_ = false;
1152 stopTime_ = mone;
1153 h0_safety_ = mone;
1154 h0_max_factor_ = mone;
1155 h_phase0_incr_ = mone;
1156 h_max_inv_ = mone;
1157 Tkm1_Tk_safety_ = mone;
1158 Tkp1_Tk_safety_ = mone;
1159 r_factor_ = mone;
1160 r_safety_ = mone;
1161 r_fudge_ = mone;
1162 r_min_ = mone;
1163 r_max_ = mone;
1164 r_hincr_test_ = mone;
1165 r_hincr_ = mone;
1166 max_LET_fail_ = -1;
1167 minTimeStep_ = mone;
1168 maxTimeStep_ = mone;
1169 newtonConvergenceStatus_ = -1;
1170}
1171
1172template<class Scalar>
1173int ImplicitBDFStepperStepControl<Scalar>::getMinOrder() const
1174{
1175 TEUCHOS_TEST_FOR_EXCEPTION(
1176 stepControlState_ == UNINITIALIZED, std::logic_error,
1177 "Error, attempting to call getMinOrder before intiialization!\n"
1178 );
1179 return(minOrder_);
1180}
1181
1182template<class Scalar>
1183int ImplicitBDFStepperStepControl<Scalar>::getMaxOrder() const
1184{
1185 TEUCHOS_TEST_FOR_EXCEPTION(
1186 stepControlState_ == UNINITIALIZED, std::logic_error,
1187 "Error, attempting to call getMaxOrder before initialization!\n"
1188 );
1189 return(maxOrder_);
1190}
1191
1192//
1193// Explicit Instantiation macro
1194//
1195// Must be expanded from within the Rythmos namespace!
1196//
1197
1198#define RYTHMOS_IMPLICITBDF_STEPPER_STEPCONTROL_INSTANT(SCALAR) \
1199 template class ImplicitBDFStepperStepControl< SCALAR >;
1200
1201
1202} // namespace Rythmos
1203#endif // Rythmos_IMPLICITBDF_STEPPER_STEP_CONTROL_DEF_H
1204