Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_ImplicitRKStepper_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_IMPLICIT_RK_STEPPER_DEF_H
30#define Rythmos_IMPLICIT_RK_STEPPER_DEF_H
31
32#include "Rythmos_ImplicitRKStepper_decl.hpp"
33
34#include "Rythmos_StepperHelpers.hpp"
35#include "Rythmos_SingleResidualModelEvaluator.hpp"
36#include "Rythmos_ImplicitRKModelEvaluator.hpp"
37#include "Rythmos_DiagonalImplicitRKModelEvaluator.hpp"
38#include "Rythmos_RKButcherTableau.hpp"
39#include "Rythmos_RKButcherTableauHelpers.hpp"
40
41#include "Thyra_ModelEvaluatorHelpers.hpp"
42#include "Thyra_ProductVectorSpaceBase.hpp"
43#include "Thyra_AssertOp.hpp"
44#include "Thyra_TestingTools.hpp"
45#include "Rythmos_ImplicitBDFStepperRampingStepControl.hpp"
46#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
47#include "Teuchos_as.hpp"
48
49
50namespace Rythmos {
51
56template<class Scalar>
57RCP<ImplicitRKStepper<Scalar> >
59{
60 RCP<ImplicitRKStepper<Scalar> > stepper(new ImplicitRKStepper<Scalar>());
61 return stepper;
62}
63
64template<class Scalar>
65RCP<ImplicitRKStepper<Scalar> >
66implicitRKStepper(
67 const RCP<const Thyra::ModelEvaluator<Scalar> >& model,
68 const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
69 const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >& irk_W_factory,
70 const RCP<const RKButcherTableauBase<Scalar> >& irkbt
71 )
72{
73 RCP<ImplicitRKStepper<Scalar> > stepper(new ImplicitRKStepper<Scalar>());
74
75 stepper->setModel(model);
76 stepper->setSolver(solver);
77 stepper->set_W_factory(irk_W_factory);
78 stepper->setRKButcherTableau(irkbt);
79 return stepper;
80}
81
82
83// ////////////////////////////
84// Defintions
85
86
87// Constructors, intializers, Misc.
88
89
90template<class Scalar>
92{
93 this->defaultInitializeAll_();
94 irkButcherTableau_ = rKButcherTableau<Scalar>();
95 numSteps_ = 0;
96}
97
98template<class Scalar>
100{
101 isInitialized_ = false;
102 model_ = Teuchos::null;
103 solver_ = Teuchos::null;
104 irk_W_factory_ = Teuchos::null;
105 paramList_ = Teuchos::null;
106 //basePoint_;
107 x_ = Teuchos::null;
108 x_old_ = Teuchos::null;
109 x_dot_ = Teuchos::null;
110 //timeRange_;
111 irkModel_ = Teuchos::null;
112 irkButcherTableau_ = Teuchos::null;
113 isDirk_ = false;
114 numSteps_ = -1;
115 haveInitialCondition_ = false;
116 x_stage_bar_ = Teuchos::null;
117}
118
119template<class Scalar>
120void ImplicitRKStepper<Scalar>::set_W_factory(
121 const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > &irk_W_factory
122 )
123{
124 TEUCHOS_ASSERT( !is_null(irk_W_factory) );
125 irk_W_factory_ = irk_W_factory;
126}
127
128template<class Scalar>
129RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > ImplicitRKStepper<Scalar>::get_W_factory() const
130{
131 return irk_W_factory_;
132}
133
134// Overridden from SolverAcceptingStepperBase
135
136
137template<class Scalar>
139 const RCP<Thyra::NonlinearSolverBase<Scalar> > &solver
140 )
141{
142 TEUCHOS_TEST_FOR_EXCEPT(is_null(solver));
143 solver_ = solver;
144}
145
146
147template<class Scalar>
148RCP<Thyra::NonlinearSolverBase<Scalar> >
150{
151 return solver_;
152}
153
154
155template<class Scalar>
156RCP<const Thyra::NonlinearSolverBase<Scalar> >
158{
159 return solver_;
160}
161
162
163// Overridden from StepperBase
164
165
166template<class Scalar>
168{
169 return true;
170}
171
172template<class Scalar>
174{
175 return true;
176}
177
178
179template<class Scalar>
180RCP<StepperBase<Scalar> >
182{
183 // Just use the interface to clone the algorithm in a basically
184 // uninitialized state
185 RCP<ImplicitRKStepper<Scalar> >
186 stepper = Teuchos::rcp(new ImplicitRKStepper<Scalar>());
187
188 if (!is_null(model_)) {
189 stepper->setModel(model_); // Shallow copy is okay!
190 }
191
192 if (!is_null(irkButcherTableau_)) {
193 // 06/16/09 tscoffe: should we clone the RKBT here?
194 stepper->setRKButcherTableau(irkButcherTableau_);
195 }
196
197 if (!is_null(solver_)) {
198 stepper->setSolver(solver_->cloneNonlinearSolver().assert_not_null());
199 }
200
201 if (!is_null(irk_W_factory_)) {
202 // 06/16/09 tscoffe: should we clone the W_factory here?
203 stepper->set_W_factory(irk_W_factory_);
204 }
205
206 if (!is_null(paramList_)) {
207 stepper->setParameterList(Teuchos::parameterList(*paramList_));
208 }
209
210 return stepper;
211}
212
213template<class Scalar>
215{
216 return(stepControl_);
217}
218
219template<class Scalar>
220RCP<const StepControlStrategyBase<Scalar> > ImplicitRKStepper<Scalar>::getStepControlStrategy() const
221{
222 return(stepControl_);
223}
224
225template<class Scalar>
227{
228 TEUCHOS_TEST_FOR_EXCEPTION(stepControl == Teuchos::null,std::logic_error,"Error, stepControl == Teuchos::null!\n");
229 stepControl_ = stepControl;
230}
231
232template<class Scalar>
234 const RCP<const Thyra::ModelEvaluator<Scalar> >& model
235 )
236{
237 TEUCHOS_TEST_FOR_EXCEPT(is_null(model));
238 assertValidModel( *this, *model );
239 model_ = model;
240}
241
242
243template<class Scalar>
245 const RCP<Thyra::ModelEvaluator<Scalar> >& model
246 )
247{
248 this->setModel(model); // TODO 09/09/09 tscoffe: use ConstNonconstObjectContainer!
249}
250
251
252template<class Scalar>
253RCP<const Thyra::ModelEvaluator<Scalar> >
255{
256 return model_;
257}
258
259
260template<class Scalar>
261RCP<Thyra::ModelEvaluator<Scalar> >
263{
264 return Teuchos::null;
265}
266
267
268template<class Scalar>
270 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
271 )
272{
273
274 typedef ScalarTraits<Scalar> ST;
275 typedef Thyra::ModelEvaluatorBase MEB;
276
277 basePoint_ = initialCondition;
278
279 // x
280
281 RCP<const Thyra::VectorBase<Scalar> >
282 x_init = initialCondition.get_x();
283
284#ifdef HAVE_RYTHMOS_DEBUG
285 TEUCHOS_TEST_FOR_EXCEPTION(
286 is_null(x_init), std::logic_error,
287 "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
288 "then x can not be null!" );
289#endif
290
291 x_ = x_init->clone_v();
292
293 // for the embedded RK method
294 xhat_ = x_init->clone_v();
295 ee_ = x_init->clone_v();
296
297 // x_dot
298
299 x_dot_ = createMember(x_->space());
300
301 RCP<const Thyra::VectorBase<Scalar> >
302 x_dot_init = initialCondition.get_x_dot();
303
304 if (!is_null(x_dot_init))
305 assign(x_dot_.ptr(),*x_dot_init);
306 else
307 assign(x_dot_.ptr(),ST::zero());
308
309 // t
310
311 const Scalar t =
312 (
313 initialCondition.supports(MEB::IN_ARG_t)
314 ? initialCondition.get_t()
315 : ST::zero()
316 );
317
318 timeRange_ = timeRange(t,t);
319
320 // x_old
321 x_old_ = x_->clone_v();
322
323 haveInitialCondition_ = true;
324
325}
326
327
328template<class Scalar>
329Thyra::ModelEvaluatorBase::InArgs<Scalar>
331{
332 return basePoint_;
333}
334
335
336template<class Scalar>
337Scalar ImplicitRKStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepSizeType)
338{
339 Scalar stepSizeTaken;
340 using Teuchos::as;
341 using Teuchos::incrVerbLevel;
342 typedef Thyra::NonlinearSolverBase<Scalar> NSB;
343 typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
344
345 RCP<FancyOStream> out = this->getOStream();
346 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
347 Teuchos::OSTab ostab(out,1,"takeStep");
348 VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
349
350 // not needed for this
351 int desiredOrder;
352
353 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
354 *out
355 << "\nEntering "
356 << Teuchos::TypeNameTraits<ImplicitRKStepper<Scalar> >::name()
357 << "::takeStep("<<dt<<","<<toString(stepSizeType)<<") ...\n";
358 }
359
360 if (!isInitialized_) {
361 initialize_();
362 }
363 if (stepSizeType == STEP_TYPE_FIXED) {
364 stepSizeTaken = takeFixedStep_(dt , stepSizeType);
365 return stepSizeTaken;
366 } else {
367 isVariableStep_ = true;
368 stepControl_->setOStream(out);
369 stepControl_->setVerbLevel(verbLevel);
370
371 rkNewtonConvergenceStatus_ = -1;
372
373 while (rkNewtonConvergenceStatus_ < 0){
374
375 stepControl_->setRequestedStepSize(*this, dt, stepSizeType);
376 stepControl_->nextStepSize(*this, &dt, &stepSizeType, &desiredOrder);
377
378 stepSizeTaken = takeVariableStep_(dt, stepSizeType);
379
380 }
381 return stepSizeTaken;
382 }
383
384}
385
386
387template<class Scalar>
388Scalar ImplicitRKStepper<Scalar>::takeFixedStep_(Scalar dt, StepSizeType stepSizeType)
389{
390 using Teuchos::as;
391 using Teuchos::incrVerbLevel;
392 typedef ScalarTraits<Scalar> ST;
393 typedef Thyra::NonlinearSolverBase<Scalar> NSB;
394 typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
395
396 RCP<FancyOStream> out = this->getOStream();
397 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
398 Teuchos::OSTab ostab(out,1,"takeStep");
399 VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
400
401 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
402 *out
403 << "\nEntering "
404 << Teuchos::TypeNameTraits<ImplicitRKStepper<Scalar> >::name()
405 << "::takeFixedStep_("<<dt<<","<<toString(stepSizeType)<<") ...\n";
406 }
407
408 if (!isInitialized_) {
409 initialize_();
410 }
411
412 TEUCHOS_TEST_FOR_EXCEPT( stepSizeType != STEP_TYPE_FIXED ); // ToDo: Handle variable case later
413
414 // A) Set up the IRK ModelEvaluator so that it can represent the time step
415 // equation to be solved.
416
417 // Set irkModel_ with x_old_, t_old_, and dt
418 V_V( x_old_.ptr(), *x_ );
419 Scalar current_dt = dt;
420 Scalar t = timeRange_.upper();
421
422 // B) Solve the timestep equation
423
424 // Set the guess for the stage derivatives to zero (unless we can think of
425 // something better)
426 V_S( Teuchos::rcp_dynamic_cast<Thyra::VectorBase<Scalar> >(x_stage_bar_).ptr(), ST::zero() );
427
428 if (!isDirk_) { // General Implicit RK Case:
429 RCP<ImplicitRKModelEvaluator<Scalar> > firkModel_ =
430 Teuchos::rcp_dynamic_cast<ImplicitRKModelEvaluator<Scalar> >(irkModel_,true);
431 firkModel_->setTimeStepPoint( x_old_, t, current_dt );
432
433 // Solve timestep equation
434 solver_->solve( &*x_stage_bar_ );
435
436 } else { // Diagonal Implicit RK Case:
437
438 RCP<DiagonalImplicitRKModelEvaluator<Scalar> > dirkModel_ =
439 Teuchos::rcp_dynamic_cast<DiagonalImplicitRKModelEvaluator<Scalar> >(irkModel_,true);
440 dirkModel_->setTimeStepPoint( x_old_, t, current_dt );
441 int numStages = irkButcherTableau_->numStages();
442 for (int stage=0 ; stage < numStages ; ++stage) {
443 dirkModel_->setCurrentStage(stage);
444 solver_->solve( &*(x_stage_bar_->getNonconstVectorBlock(stage)) );
445 dirkModel_->setStageSolution( stage, *(x_stage_bar_->getVectorBlock(stage)) );
446 }
447
448 }
449
450 // C) Complete the step ...
451
452 // Combine the stage derivatives with the Butcher tableau "b" vector to obtain the solution at the final time.
453 // x_{k+1} = x_k + dt*sum_{i}^{p}(b_i*x_stage_bar_[i])
454
455 assembleIRKSolution( irkButcherTableau_->b(), current_dt, *x_old_, *x_stage_bar_,
456 outArg(*x_)
457 );
458
459 // Update time range
460 timeRange_ = timeRange(t,t+current_dt);
461 numSteps_++;
462
463 return current_dt;
464
465}
466
467template<class Scalar>
468Scalar ImplicitRKStepper<Scalar>::takeVariableStep_(Scalar dt, StepSizeType stepSizeType)
469{
470 using Teuchos::as;
471 using Teuchos::incrVerbLevel;
472 typedef ScalarTraits<Scalar> ST;
473 typedef Thyra::NonlinearSolverBase<Scalar> NSB;
474 typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
475
476 RCP<FancyOStream> out = this->getOStream();
477 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
478 Teuchos::OSTab ostab(out,1,"takeStep");
479 VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
480
481 AttemptedStepStatusFlag status;
482 Scalar dt_old = dt;
483
484 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
485 *out
486 << "\nEntering "
487 << Teuchos::TypeNameTraits<ImplicitRKStepper<Scalar> >::name()
488 << "::takeVariableStep_("<<dt<<","<<toString(stepSizeType)<<") ...\n";
489 }
490
491 if (!isInitialized_) {
492 initialize_();
493 }
494
495 // A) Set up the IRK ModelEvaluator so that it can represent the time step
496 // equation to be solved.
497
498 // Set irkModel_ with x_old_, t_old_, and dt
499 V_V( x_old_.ptr(), *x_ );
500 Scalar current_dt = dt;
501 Scalar t = timeRange_.upper();
502 Scalar dt_to_return;
503
504 // B) Solve the timestep equation
505
506 // Set the guess for the stage derivatives to zero (unless we can think of
507 // something better)
508 V_S( Teuchos::rcp_dynamic_cast<Thyra::VectorBase<Scalar> >(x_stage_bar_).ptr(), ST::zero() );
509
510 if (!isDirk_) { // General Implicit RK Case:
511 RCP<ImplicitRKModelEvaluator<Scalar> > firkModel_ =
512 Teuchos::rcp_dynamic_cast<ImplicitRKModelEvaluator<Scalar> >(irkModel_,true);
513 firkModel_->setTimeStepPoint( x_old_, t, current_dt );
514
515 // Solve timestep equation
516 solver_->solve( &*x_stage_bar_ );
517
518 } else { // Diagonal Implicit RK Case:
519
520 RCP<DiagonalImplicitRKModelEvaluator<Scalar> > dirkModel_ =
521 Teuchos::rcp_dynamic_cast<DiagonalImplicitRKModelEvaluator<Scalar> >(irkModel_,true);
522 dirkModel_->setTimeStepPoint( x_old_, t, current_dt );
523 int numStages = irkButcherTableau_->numStages();
524 for (int stage=0 ; stage < numStages ; ++stage) {
525 dirkModel_->setCurrentStage(stage);
526 nonlinearSolveStatus_ = solver_->solve( &*(x_stage_bar_->getNonconstVectorBlock(stage)) );
527
528 if (nonlinearSolveStatus_.solveStatus == Thyra::SOLVE_STATUS_CONVERGED) {
529 rkNewtonConvergenceStatus_ = 0;
530 } else {
531 rkNewtonConvergenceStatus_ = -1;
532 }
533
534 // for now setCorrection just sets the rkNewtonConvergenceStatus_ in the stepControl
535 // and this is used by acceptStep method of the stepControl
536
537 stepControl_->setCorrection(*this, (x_stage_bar_->getNonconstVectorBlock(stage)), Teuchos::null , rkNewtonConvergenceStatus_);
538 bool stepPass = stepControl_->acceptStep(*this, &LETvalue_);
539
540 if (!stepPass) { // stepPass = false
541 stepLETStatus_ = STEP_LET_STATUS_FAILED;
542 rkNewtonConvergenceStatus_ = -1; // just making sure here
543 break; // leave the for loop
544 } else { // stepPass = true
545 stepLETStatus_ = STEP_LET_STATUS_PASSED;
546 dirkModel_->setStageSolution( stage, *(x_stage_bar_->getVectorBlock(stage)) );
547 rkNewtonConvergenceStatus_ = 0; // just making sure here
548 }
549 }
550 // if none of the stages failed, then I can complete the step
551 }
552
553 // check the nonlinearSolveStatus
554 if ( rkNewtonConvergenceStatus_ == 0) {
555
556 /*
557 * if the solver has converged, then I can go ahead and combine the stage solutions
558 * and get the new solution
559 */
560
561 // C) Complete the step ...
562
563 // Combine the stage derivatives with the Butcher tableau "b" vector to obtain the solution at the final time.
564 // x_{k+1} = x_k + dt*sum_{i}^{p}(b_i*x_stage_bar_[i])
565
566 assembleIRKSolution( irkButcherTableau_->b(), current_dt, *x_old_, *x_stage_bar_,
567 outArg(*x_)
568 );
569
570 //if using embedded method, estimate LTE
571 if (irkButcherTableau_->isEmbeddedMethod() ){
572
573 assembleIRKSolution( irkButcherTableau_->bhat(), current_dt, *x_old_, *x_stage_bar_,
574 outArg(*xhat_)
575 );
576
577 // ee_ = (x_ - xhat_)
578 Thyra::V_VmV(ee_.ptr(), *x_, *xhat_);
579 stepControl_->setCorrection(*this, x_, ee_ , rkNewtonConvergenceStatus_);
580
581 bool stepPass = stepControl_->acceptStep(*this, &LETvalue_);
582
583 if (!stepPass) { // stepPass = false
584 stepLETStatus_ = STEP_LET_STATUS_FAILED;
585 rkNewtonConvergenceStatus_ = -1; // just making sure here
586 } else { // stepPass = true
587 stepLETStatus_ = STEP_LET_STATUS_PASSED;
588 rkNewtonConvergenceStatus_ = 0; // just making sure here
589 }
590 }
591 }
592
593 if (rkNewtonConvergenceStatus_ == 0) {
594
595 // Update time range
596 timeRange_ = timeRange(t,t+current_dt);
597 numSteps_++;
598
599 // completeStep only if the none of the stage solution's failed to converged
600 stepControl_->completeStep(*this);
601
602 dt_to_return = current_dt;
603
604 } else {
605 rkNewtonConvergenceStatus_ = -1;
606 status = stepControl_-> rejectStep(*this); // reject the stage value
607 (void) status; // avoid "set but not used" build warning
608 dt_to_return = dt_old;
609 }
610
611 return dt_to_return;
612
613}
614
615
616template<class Scalar>
618{
619 StepStatus<Scalar> stepStatus;
620
621 if (!isInitialized_) {
622 stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
623 stepStatus.message = "This stepper is uninitialized.";
624// return stepStatus;
625 }
626 else if (numSteps_ > 0) {
627 stepStatus.stepStatus = STEP_STATUS_CONVERGED;
628 }
629 else {
630 stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
631 }
632 stepStatus.stepSize = timeRange_.length();
633 stepStatus.order = irkButcherTableau_->order();
634 stepStatus.time = timeRange_.upper();
635 if(Teuchos::nonnull(x_))
636 stepStatus.solution = x_;
637 else
638 stepStatus.solution = Teuchos::null;
639 stepStatus.solutionDot = Teuchos::null;
640 return(stepStatus);
641}
642
643
644// Overridden from InterpolationBufferBase
645
646
647template<class Scalar>
648RCP<const Thyra::VectorSpaceBase<Scalar> >
650{
651 return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
652}
653
654
655template<class Scalar>
657 const Array<Scalar>& /* time_vec */
658 ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& /* x_vec */
659 ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& /* xdot_vec */
660 )
661{
662 TEUCHOS_TEST_FOR_EXCEPT(true);
663}
664
665
666template<class Scalar>
668{
669 return timeRange_;
670}
671
672
673template<class Scalar>
675 const Array<Scalar>& time_vec
676 ,Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec
677 ,Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
678 ,Array<ScalarMag>* accuracy_vec) const
679{
680 using Teuchos::constOptInArg;
681 using Teuchos::null;
682 TEUCHOS_ASSERT(haveInitialCondition_);
683 defaultGetPoints<Scalar>(
684 timeRange_.lower(), constOptInArg(*x_old_),
685 Ptr<const VectorBase<Scalar> >(null), // Sun
686 timeRange_.upper(), constOptInArg(*x_),
687 Ptr<const VectorBase<Scalar> >(null), // Sun
688 time_vec,
689 ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
690 Ptr<InterpolatorBase<Scalar> >(null) // For Sun
691 );
692 // 04/17/09 tscoffe: Currently, we don't have x_dot to pass out (TODO)
693}
694
695
696template<class Scalar>
697void ImplicitRKStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
698{
699 TEUCHOS_ASSERT( time_vec != NULL );
700 time_vec->clear();
701 if (!haveInitialCondition_) {
702 return;
703 }
704 time_vec->push_back(timeRange_.lower());
705 if (numSteps_ > 0) {
706 time_vec->push_back(timeRange_.upper());
707 }
708}
709
710
711template<class Scalar>
712void ImplicitRKStepper<Scalar>::removeNodes(Array<Scalar>& /* time_vec */)
713{
714 TEUCHOS_TEST_FOR_EXCEPT(true);
715}
716
717
718template<class Scalar>
720{
721 return irkButcherTableau_->order();
722}
723
724
725// Overridden from Teuchos::ParameterListAcceptor
726
727
728template <class Scalar>
730 RCP<ParameterList> const& paramList
731 )
732{
733 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
734 paramList->validateParametersAndSetDefaults(*this->getValidParameters());
735 paramList_ = paramList;
736 Teuchos::readVerboseObjectSublist(&*paramList_,this);
737}
738
739
740template <class Scalar>
741RCP<ParameterList>
743{
744 return(paramList_);
745}
746
747
748template <class Scalar>
749RCP<ParameterList>
751{
752 RCP<ParameterList>
753 temp_param_list = paramList_;
754 paramList_ = Teuchos::null;
755 return(temp_param_list);
756}
757
758
759template<class Scalar>
760RCP<const ParameterList>
762{
763 static RCP<const ParameterList> validPL;
764 if (is_null(validPL)) {
765 RCP<ParameterList> pl = Teuchos::parameterList();
766 if (isVariableStep_){
767 pl->sublist(RythmosStepControlSettings_name);
768 }
769 Teuchos::setupVerboseObjectSublist(&*pl);
770 validPL = pl;
771 }
772 return validPL;
773}
774
775
776// Overridden from Teuchos::Describable
777
778
779template<class Scalar>
781 FancyOStream &out,
782 const Teuchos::EVerbosityLevel verbLevel
783 ) const
784{
785 using std::endl;
786 using Teuchos::as;
787 if (!isInitialized_) {
788 out << this->description() << " : This stepper is not initialized yet" << std::endl;
789 return;
790 }
791 if (
792 as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
793 ||
794 as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
795 )
796 {
797 out << this->description() << ":" << endl;
798 Teuchos::OSTab tab(out);
799 out << "model = " << Teuchos::describe(*model_,verbLevel);
800 out << "solver = " << Teuchos::describe(*solver_,verbLevel);
801 out << "irk_W_factory = " << Teuchos::describe(*irk_W_factory_,verbLevel);
802 }
803}
804
805
806// private
807
808
809template <class Scalar>
811{
812
813 // typedef ScalarTraits<Scalar> ST; // unused
814 using Teuchos::rcp_dynamic_cast;
815
816 TEUCHOS_TEST_FOR_EXCEPT(is_null(model_));
817 TEUCHOS_TEST_FOR_EXCEPT(is_null(solver_));
818 TEUCHOS_TEST_FOR_EXCEPT(irkButcherTableau_->numStages() == 0);
819 TEUCHOS_ASSERT(haveInitialCondition_);
820
821#ifdef HAVE_RYTHMOS_DEBUG
822 THYRA_ASSERT_VEC_SPACES(
823 "Rythmos::ImplicitRKStepper::initialize_(...)",
824 *x_->space(), *model_->get_x_space() );
825#endif
826
827 if (isVariableStep_ ) {
828 // Initialize StepControl
829
830 isEmbeddedRK_ = irkButcherTableau_->isEmbeddedMethod(); // determine if RK method is an embedded method
831 if (stepControl_ == Teuchos::null) {
832 RCP<ImplicitBDFStepperRampingStepControl<Scalar> > rkStepControl =
834 //RCP<StepControlStrategyBase<Scalar> > rkStepControl =
835 //Teuchos::rcp(new StepControlStrategyBase<Scalar>());
836 RCP<Teuchos::ParameterList> stepControlPL =
837 Teuchos::sublist(paramList_ , RythmosStepControlSettings_name);
838 rkStepControl->setParameterList(stepControlPL);
839 this->setStepControlStrategy(rkStepControl);
840 stepControl_->initialize(*this);
841 }
842 }
843
844 // Set up the IRK mdoel
845
846 if (!isDirk_) { // General Implicit RK
847 TEUCHOS_TEST_FOR_EXCEPT(is_null(irk_W_factory_));
848 irkModel_ = implicitRKModelEvaluator(
849 model_,basePoint_,irk_W_factory_,irkButcherTableau_);
850 } else { // Diagonal Implicit RK
851 irkModel_ = diagonalImplicitRKModelEvaluator(
852 model_,basePoint_,irk_W_factory_,irkButcherTableau_);
853 }
854
855 solver_->setModel(irkModel_);
856
857 // Set up the vector of stage derivatives ...
858 const int numStages = irkButcherTableau_->numStages();
859 RCP<const Thyra::ProductVectorSpaceBase<Scalar> > pvs = productVectorSpace(model_->get_x_space(),numStages);
860 RCP<const Thyra::VectorSpaceBase<Scalar> > vs = rcp_dynamic_cast<const Thyra::VectorSpaceBase<Scalar> >(pvs,true);
861 x_stage_bar_ = rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(createMember(vs),true);
862// x_stage_bar_ = rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(
863// createMember(irkModel_->get_x_space())
864// );
865
866 isInitialized_ = true;
867
868}
869
870template <class Scalar>
871void ImplicitRKStepper<Scalar>::setRKButcherTableau( const RCP<const RKButcherTableauBase<Scalar> > &rkButcherTableau )
872{
873 TEUCHOS_ASSERT( !is_null(rkButcherTableau) );
874 TEUCHOS_TEST_FOR_EXCEPTION( isInitialized_, std::logic_error,
875 "Error! The RK Butcher Tableau cannot be changed after internal initialization!"
876 );
877 validateIRKButcherTableau(*rkButcherTableau);
878 irkButcherTableau_ = rkButcherTableau;
879 E_RKButcherTableauTypes rkType = determineRKBTType<Scalar>(*irkButcherTableau_);
880 if (
881 (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_DIRK)
882 || (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_SDIRK)
883 || (irkButcherTableau_->numStages() == 1)
884 )
885 {
886 isDirk_ = true;
887 }
888}
889
890template <class Scalar>
891RCP<const RKButcherTableauBase<Scalar> > ImplicitRKStepper<Scalar>::getRKButcherTableau() const
892{
893 return irkButcherTableau_;
894}
895
896template<class Scalar>
898{
899 TEUCHOS_TEST_FOR_EXCEPTION(isInitialized_, std::logic_error,
900 "Error! Cannot change DIRK flag after internal initialization is completed\n"
901 );
902 if (isDirk == true) {
903 E_RKButcherTableauTypes rkType = determineRKBTType<Scalar>(*irkButcherTableau_);
904 bool RKBT_is_DIRK = (
905 (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_DIRK)
906 || (rkType == RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_SDIRK)
907 || (irkButcherTableau_->numStages() == 1)
908 );
909 TEUCHOS_TEST_FOR_EXCEPTION( !RKBT_is_DIRK, std::logic_error,
910 "Error! Cannot set DIRK flag on a non-DIRK RK Butcher Tableau\n"
911 );
912 } else { // isDirk = false;
913 isDirk_ = isDirk;
914 }
915}
916
917//
918// Explicit Instantiation macro
919//
920// Must be expanded from within the Rythmos namespace!
921//
922
923#define RYTHMOS_IMPLICIT_RK_STEPPER_INSTANT(SCALAR) \
924 \
925 template class ImplicitRKStepper< SCALAR >; \
926 \
927 template RCP< ImplicitRKStepper< SCALAR > > \
928 implicitRKStepper(); \
929 \
930 template RCP< ImplicitRKStepper< SCALAR > > \
931 implicitRKStepper( \
932 const RCP<const Thyra::ModelEvaluator< SCALAR > >& model, \
933 const RCP<Thyra::NonlinearSolverBase< SCALAR > >& solver, \
934 const RCP<Thyra::LinearOpWithSolveFactoryBase< SCALAR > >& irk_W_factory, \
935 const RCP<const RKButcherTableauBase< SCALAR > >& irkbt \
936 );
937
938} // namespace Rythmos
939
940#endif //Rythmos_IMPLICIT_RK_STEPPER_DEF_H
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< const RKButcherTableauBase< Scalar > > getRKButcherTableau() const
TimeRange< Scalar > getTimeRange() const
void getNodes(Array< Scalar > *time_vec) const
RCP< const Thyra::NonlinearSolverBase< Scalar > > getSolver() const
RCP< Thyra::ModelEvaluator< Scalar > > getNonconstModel()
void removeNodes(Array< Scalar > &time_vec)
void setNonconstModel(const RCP< Thyra::ModelEvaluator< Scalar > > &model)
void setSolver(const RCP< Thyra::NonlinearSolverBase< Scalar > > &solver)
void setParameterList(RCP< ParameterList > const &paramList)
RCP< StepperBase< Scalar > > cloneStepperAlgorithm() const
RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
void describe(FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
void setRKButcherTableau(const RCP< const RKButcherTableauBase< Scalar > > &rkButcherTableau)
RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
RCP< Thyra::NonlinearSolverBase< Scalar > > getNonconstSolver()
void setModel(const RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Thyra::ModelEvaluatorBase::InArgs< Scalar > getInitialCondition() const
void setInitialCondition(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &initialCondition)
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< const ParameterList > getValidParameters() const
void setStepControlStrategy(const RCP< StepControlStrategyBase< Scalar > > &stepControlStrategy)
Scalar takeStep(Scalar dt, StepSizeType flag)
const StepStatus< Scalar > getStepStatus() const
RCP< const StepControlStrategyBase< Scalar > > getStepControlStrategy() const
RCP< StepControlStrategyBase< Scalar > > getNonconstStepControlStrategy()
RCP< ImplicitRKStepper< Scalar > > implicitRKStepper()
Nonmember constructor.
Base strategy class for interpolation functionality.
The member functions in the StepControlStrategyBase move you between these states in the following fa...
Represent a time range.
RCP< const Thyra::VectorBase< Scalar > > solution
RCP< const Thyra::VectorBase< Scalar > > solutionDot