Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_ExplicitRKStepper_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_ExplicitRK_STEPPER_DEF_H
30#define Rythmos_ExplicitRK_STEPPER_DEF_H
31
32#include "Rythmos_ExplicitRKStepper_decl.hpp"
33
34#include "Rythmos_RKButcherTableau.hpp"
35#include "Rythmos_RKButcherTableauHelpers.hpp"
36#include "Rythmos_RKButcherTableauBuilder.hpp"
37#include "Rythmos_StepperHelpers.hpp"
38#include "Rythmos_LinearInterpolator.hpp"
39#include "Rythmos_InterpolatorBaseHelpers.hpp"
40
41#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
42
43#include "Thyra_ModelEvaluatorHelpers.hpp"
44#include "Thyra_MultiVectorStdOps.hpp"
45#include "Thyra_VectorStdOps.hpp"
46
47
48namespace Rythmos {
49
50// Non-member constructors
51template<class Scalar>
52RCP<ExplicitRKStepper<Scalar> > explicitRKStepper()
53{
54 RCP<ExplicitRKStepper<Scalar> > stepper = rcp(new ExplicitRKStepper<Scalar>());
55 return stepper;
56}
57
58template<class Scalar>
59RCP<ExplicitRKStepper<Scalar> > explicitRKStepper(
60 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model
61 )
62{
63 RCP<RKButcherTableauBase<Scalar> > rkbt = createRKBT<Scalar>("Explicit 4 Stage");
64 //RCP<RKButcherTableauBase<Scalar> > rkbt = rcp(new Explicit4Stage4thOrder_RKBT<Scalar>());
65 RCP<ExplicitRKStepper<Scalar> > stepper = rcp(new ExplicitRKStepper<Scalar>());
66 stepper->setModel(model);
67 stepper->setRKButcherTableau(rkbt);
68 return stepper;
69}
70
71template<class Scalar>
72RCP<ExplicitRKStepper<Scalar> > explicitRKStepper(
73 const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
74 const RCP<const RKButcherTableauBase<Scalar> >& rkbt
75 )
76{
77 RCP<ExplicitRKStepper<Scalar> > stepper = rcp(new ExplicitRKStepper<Scalar>());
78 stepper->setModel(model);
79 stepper->setRKButcherTableau(rkbt);
80 return stepper;
81}
82
83template<class Scalar>
85{
86 this->defaultInitializeAll_();
87 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
88 out->precision(15);
89 erkButcherTableau_ = rKButcherTableau<Scalar>();
90 numSteps_ = 0;
91}
92
93template<class Scalar>
95{
96 model_ = Teuchos::null;
97 solution_vector_ = Teuchos::null;
98 solution_vector_old_ = Teuchos::null;
99 //k_vector_;
100 ktemp_vector_ = Teuchos::null;
101 //basePoint_;
102 erkButcherTableau_ = Teuchos::null;
103 t_ = ST::nan();
104 t_old_ = ST::nan();
105 dt_ = ST::nan();
106 numSteps_ = -1;
107 parameterList_ = Teuchos::null;
108 isInitialized_ = false;
109 haveInitialCondition_ = false;
110}
111
112template<class Scalar>
113void ExplicitRKStepper<Scalar>::setRKButcherTableau(const RCP<const RKButcherTableauBase<Scalar> > &rkbt)
114{
115 TEUCHOS_ASSERT( !is_null(rkbt) );
116 validateERKButcherTableau(*rkbt);
117 int numStages_old = erkButcherTableau_->numStages();
118 int numStages_new = rkbt->numStages();
119 TEUCHOS_TEST_FOR_EXCEPTION( numStages_new == 0, std::logic_error,
120 "Error! The Runge-Kutta Butcher tableau has no stages!"
121 );
122 if (!is_null(model_)) {
123 int numNewStages = numStages_new - numStages_old;
124 if ( numNewStages > 0 ) {
125 k_vector_.reserve(numStages_new);
126 for (int i=0 ; i<numNewStages ; ++i) {
127 k_vector_.push_back(Thyra::createMember(model_->get_f_space()));
128 }
129 }
130 }
131 erkButcherTableau_ = rkbt;
132}
133
134template<class Scalar>
135RCP<const RKButcherTableauBase<Scalar> > ExplicitRKStepper<Scalar>::getRKButcherTableau() const
136{
137 return erkButcherTableau_;
138}
139
140template<class Scalar>
142{
143 if (!isInitialized_) {
144 TEUCHOS_ASSERT( !is_null(model_) );
145 TEUCHOS_ASSERT( !is_null(erkButcherTableau_) );
146 TEUCHOS_ASSERT( haveInitialCondition_ );
147 TEUCHOS_TEST_FOR_EXCEPTION( erkButcherTableau_->numStages() == 0, std::logic_error,
148 "Error! The Runge-Kutta Butcher tableau has no stages!"
149 );
150 ktemp_vector_ = Thyra::createMember(model_->get_f_space());
151 // Initialize the stage vectors
152 int numStages = erkButcherTableau_->numStages();
153 k_vector_.reserve(numStages);
154 for (int i=0 ; i<numStages ; ++i) {
155 k_vector_.push_back(Thyra::createMember(model_->get_f_space()));
156 }
157 }
158#ifdef HAVE_RYTHMOS_DEBUG
159 THYRA_ASSERT_VEC_SPACES(
160 "Rythmos::ExplicitRKStepper::initialize_(...)",
161 *solution_vector_->space(), *model_->get_x_space() );
162#endif // HAVE_RYTHMOS_DEBUG
163 isInitialized_ = true;
164}
165
166
167template<class Scalar>
169{
170}
171
172template<class Scalar>
173Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > ExplicitRKStepper<Scalar>::get_x_space() const
174{
175 TEUCHOS_ASSERT( !is_null(model_) );
176 return(model_->get_x_space());
177}
178
179template<class Scalar>
180Scalar ExplicitRKStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepSizeType)
181{
182
183 RCP<FancyOStream> out = this->getOStream();
184 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
185 Teuchos::OSTab ostab(out,1,"takeStep");
186 Scalar stepSizeTaken;
187
188
189 // not needed for this
190 int desiredOrder;
191
192 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
193 *out
194 << "\nEntering "
195 << Teuchos::TypeNameTraits<ExplicitRKStepper<Scalar> >::name()
196 << "::takeStep("<<dt<<","<<toString(stepSizeType)<<") ...\n";
197 }
198
199
200 if (stepSizeType == STEP_TYPE_FIXED) {
201 stepSizeTaken = takeFixedStep_(dt , stepSizeType);
202 return stepSizeTaken;
203 } else {
204 isVariableStep_ = true;
205
206 stepControl_->setOStream(out);
207 stepControl_->setVerbLevel(verbLevel);
208
209 rkNewtonConvergenceStatus_ = -1;
210
211 while (rkNewtonConvergenceStatus_ < 0){
212
213 stepControl_->setRequestedStepSize(*this, dt, stepSizeType);
214 stepControl_->nextStepSize(*this, &dt, &stepSizeType, &desiredOrder);
215
216 stepSizeTaken = takeVariableStep_(dt, stepSizeType);
217
218 }
219 return stepSizeTaken;
220 }
221
222
223}
224
225
226template<class Scalar>
227Scalar ExplicitRKStepper<Scalar>::takeVariableStep_(Scalar dt, StepSizeType /* stepSizeType */)
228{
229 typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag TScalarMag;
230 this->initialize_();
231
232 // Store old solution & old time
233 V_V(solution_vector_old_.ptr(), *solution_vector_);
234 t_old_ = t_;
235 Scalar current_dt = dt;
236 Scalar t = timeRange_.upper();
237 Scalar dt_to_return;
238 AttemptedStepStatusFlag status;
239 Scalar dt_old = dt;
240
241 dt_ = dt;
242
243 int stages = erkButcherTableau_->numStages();
244 Teuchos::SerialDenseMatrix<int,Scalar> A = erkButcherTableau_->A();
245 Teuchos::SerialDenseVector<int,Scalar> b = erkButcherTableau_->b();
246 Teuchos::SerialDenseVector<int,Scalar> c = erkButcherTableau_->c();
247
248
249 // Compute stage solutions
250 for (int s=0 ; s < stages ; ++s) {
251 Thyra::assign(ktemp_vector_.ptr(), *solution_vector_); // ktemp = solution_vector
252 for (int j=0 ; j < s ; ++j) { // assuming Butcher matix is strictly lower triangular
253 if (A(s,j) != ST::zero()) {
254 Thyra::Vp_StV(ktemp_vector_.ptr(), A(s,j), *k_vector_[j]); // ktemp = ktemp + a_{s+1,j+1}*k_{j+1}
255 }
256 }
257 TScalarMag ts = t_ + c(s)*dt;
258 TScalarMag scaled_dt = (s == 0)? Scalar(dt/stages) : c(s)*dt;
259
260 // need to check here the status of the solver (the linear solve)
261 //eval_model_explicit<Scalar>(*model_,basePoint_,*ktemp_vector_,ts,Teuchos::outArg(*k_vector_[s]));
262 eval_model_explicit<Scalar>(*model_,basePoint_,*ktemp_vector_,ts,Teuchos::outArg(*k_vector_[s]), scaled_dt, c(s));
263 Thyra::Vt_S(k_vector_[s].ptr(),dt); // k_s = k_s*dt
264 }
265 // Sum for solution:
266 for (int s=0 ; s < stages ; ++s) {
267 if (b(s) != ST::zero()) {
268 Thyra::Vp_StV(solution_vector_.ptr(), b(s), *k_vector_[s]); // solution_vector += b_{s+1}*k_{s+1}
269 }
270 }
271
272 // cheat and say that the solver converged ( although no solver is needed for explicit method )
273 rkNewtonConvergenceStatus_ = 0;
274 stepControl_->setCorrection(*this, solution_vector_, Teuchos::null , rkNewtonConvergenceStatus_);
275 bool stepPass = stepControl_->acceptStep(*this, &LETvalue_);
276
277 if (erkButcherTableau_->isEmbeddedMethod() ){
278
279 Teuchos::SerialDenseVector<int,Scalar> bhat = erkButcherTableau_->bhat();
280
281 // Sum for Embedded solution:
282 for (int s=0 ; s < stages ; ++s) {
283 if (bhat(s) != ST::zero()) {
284
285 // solution_hat_vector += bhat_{s+1}*k_{s+1}
286 Thyra::Vp_StV(solution_hat_vector_.ptr(), bhat(s), *k_vector_[s]);
287 }
288 }
289
290 // ee_ = (solution_vector__ - solution_hat_vector_)
291 Thyra::V_VmV(ee_.ptr(), *solution_vector_, *solution_hat_vector_);
292
293 stepControl_->setCorrection(*this, solution_vector_, ee_ , rkNewtonConvergenceStatus_);
294 }
295 else {
296 stepControl_->setCorrection(*this, solution_vector_, Teuchos::null, rkNewtonConvergenceStatus_);
297 }
298
299 stepPass = stepControl_->acceptStep(*this, &LETvalue_);
300
301 if (!stepPass) { // stepPass = false
302 stepLETStatus_ = STEP_LET_STATUS_FAILED;
303 rkNewtonConvergenceStatus_ = -1; // just making sure here
304 } else { // stepPass = true
305 stepLETStatus_ = STEP_LET_STATUS_PASSED;
306 rkNewtonConvergenceStatus_ = 0; // just making sure here
307 }
308
309 if (rkNewtonConvergenceStatus_ == 0) {
310
311 // Update time range
312 timeRange_ = timeRange(t,t+current_dt);
313 numSteps_++;
314
315 // completeStep only if the none of the stage solution's failed to converged
316 stepControl_->completeStep(*this);
317
318 dt_to_return = current_dt;
319
320 } else {
321
322 rkNewtonConvergenceStatus_ = -1;
323 status = stepControl_-> rejectStep(*this); // reject the stage value
324 (void) status; // avoid "set but not used" build warning
325 dt_to_return = dt_old;
326 }
327
328 // update current time:
329 t_ = t_ + dt;
330 numSteps_++;
331 return( dt_to_return );
332}
333
334
335template<class Scalar>
336Scalar ExplicitRKStepper<Scalar>::takeFixedStep_(Scalar dt, StepSizeType flag)
337{
338 typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag TScalarMag;
339 this->initialize_();
340#ifdef HAVE_RYTHMOS_DEBUG
341 TEUCHOS_TEST_FOR_EXCEPTION( flag == STEP_TYPE_VARIABLE, std::logic_error,
342 "Error! ExplicitRKStepper does not support variable time steps at this time."
343 );
344#endif // HAVE_RYTHMOS_DEBUG
345 if ((flag == STEP_TYPE_VARIABLE) || (dt == ST::zero())) {
346 return(Scalar(-ST::one()));
347 }
348 // Store old solution & old time
349 V_V(solution_vector_old_.ptr(), *solution_vector_);
350 t_old_ = t_;
351
352 dt_ = dt;
353
354 int stages = erkButcherTableau_->numStages();
355 Teuchos::SerialDenseMatrix<int,Scalar> A = erkButcherTableau_->A();
356 Teuchos::SerialDenseVector<int,Scalar> b = erkButcherTableau_->b();
357 Teuchos::SerialDenseVector<int,Scalar> c = erkButcherTableau_->c();
358 // Compute stage solutions
359 for (int s=0 ; s < stages ; ++s) {
360 Thyra::assign(ktemp_vector_.ptr(), *solution_vector_); // ktemp = solution_vector
361 for (int j=0 ; j < s ; ++j) { // assuming Butcher matix is strictly lower triangular
362 if (A(s,j) != ST::zero()) {
363 Thyra::Vp_StV(ktemp_vector_.ptr(), A(s,j), *k_vector_[j]); // ktemp = ktemp + a_{s+1,j+1}*k_{j+1}
364 }
365 }
366 TScalarMag ts = t_ + c(s)*dt;
367 TScalarMag scaled_dt = (s == 0)? Scalar(dt/stages) : c(s)*dt;
368
369 eval_model_explicit<Scalar>(*model_,basePoint_,*ktemp_vector_,ts,Teuchos::outArg(*k_vector_[s]), scaled_dt, c(s));
370 Thyra::Vt_S(k_vector_[s].ptr(),dt); // k_s = k_s*dt
371 }
372 // Sum for solution:
373 for (int s=0 ; s < stages ; ++s) {
374 if (b(s) != ST::zero()) {
375 Thyra::Vp_StV(solution_vector_.ptr(), b(s), *k_vector_[s]); // solution_vector += b_{s+1}*k_{s+1}
376 }
377 }
378
379 // update current time:
380 t_ = t_ + dt;
381
382 numSteps_++;
383
384 return(dt);
385}
386
387template<class Scalar>
389{
390 StepStatus<Scalar> stepStatus;
391
392 if (!haveInitialCondition_) {
393 stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
394 } else if (numSteps_ == 0) {
395 stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
396 stepStatus.order = erkButcherTableau_->order();
397 stepStatus.time = t_;
398 stepStatus.solution = solution_vector_;
399 } else {
400 stepStatus.stepStatus = STEP_STATUS_CONVERGED;
401 stepStatus.stepSize = dt_;
402 stepStatus.order = erkButcherTableau_->order();
403 stepStatus.time = t_;
404 stepStatus.solution = solution_vector_;
405 }
406
407 return(stepStatus);
408}
409
410template<class Scalar>
412 Teuchos::FancyOStream &out,
413 const Teuchos::EVerbosityLevel verbLevel
414 ) const
415{
416 if ( (static_cast<int>(verbLevel) == static_cast<int>(Teuchos::VERB_DEFAULT) ) ||
417 (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW) )
418 ) {
419 out << this->description() << "::describe" << std::endl;
420 out << "model = " << model_->description() << std::endl;
421 out << erkButcherTableau_->numStages() << " stage Explicit RK method" << std::endl;
422 } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)) {
423 out << "solution_vector = " << std::endl;
424 out << Teuchos::describe(*solution_vector_,verbLevel);
425 } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM)) {
426 } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_HIGH)) {
427 out << "model = " << std::endl;
428 out << Teuchos::describe(*model_,verbLevel);
429 int stages = erkButcherTableau_->numStages();
430 for (int i=0 ; i<stages ; ++i) {
431 out << "k_vector[" << i << "] = " << std::endl;
432 out << Teuchos::describe(*k_vector_[i],verbLevel);
433 }
434 out << "ktemp_vector = " << std::endl;
435 out << Teuchos::describe(*ktemp_vector_,verbLevel);
436 out << "ERK Butcher Tableau A matrix: " << printMat(erkButcherTableau_->A()) << std::endl;
437 out << "ERK Butcher Tableau b vector: " << printMat(erkButcherTableau_->b()) << std::endl;
438 out << "ERK Butcher Tableau c vector: " << printMat(erkButcherTableau_->c()) << std::endl;
439 out << "t = " << t_ << std::endl;
440 }
441}
442
443template<class Scalar>
445 const Array<Scalar>& /* time_vec */
446 ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& /* x_vec */
447 ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& /* xdot_vec */
448 )
449{
450 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, addPoints is not implemented for ExplicitRKStepper at this time.\n");
451}
452
453template<class Scalar>
455{
456 if (!haveInitialCondition_) {
457 return(invalidTimeRange<Scalar>());
458 } else {
459 return(TimeRange<Scalar>(t_old_,t_));
460 }
461}
462
463template<class Scalar>
465 const Array<Scalar>& time_vec,
466 Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
467 Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
468 Array<ScalarMag>* accuracy_vec
469 ) const
470{
471 TEUCHOS_ASSERT( haveInitialCondition_ );
472 using Teuchos::constOptInArg;
473 using Teuchos::null;
474 defaultGetPoints<Scalar>(
475 t_old_, constOptInArg(*solution_vector_old_),
476 Ptr<const VectorBase<Scalar> >(null),
477 t_, constOptInArg(*solution_vector_),
478 Ptr<const VectorBase<Scalar> >(null),
479 time_vec,ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
480 Ptr<InterpolatorBase<Scalar> >(null)
481 );
482}
483
484template<class Scalar>
485void ExplicitRKStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
486{
487 TEUCHOS_ASSERT( time_vec != NULL );
488 time_vec->clear();
489 if (!haveInitialCondition_) {
490 return;
491 }
492 time_vec->push_back(t_old_);
493 if (t_ != t_old_) {
494 time_vec->push_back(t_);
495 }
496}
497
498template<class Scalar>
499void ExplicitRKStepper<Scalar>::removeNodes(Array<Scalar>& /* time_vec */)
500{
501 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for ExplicitRKStepper at this time.\n");
502}
503
504template<class Scalar>
506{
507 return(erkButcherTableau_->order());
508}
509
510template <class Scalar>
511void ExplicitRKStepper<Scalar>::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
512{
513 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
514 paramList->validateParametersAndSetDefaults(*this->getValidParameters());
515 parameterList_ = paramList;
516 Teuchos::readVerboseObjectSublist(&*parameterList_,this);
517}
518
519template <class Scalar>
520Teuchos::RCP<Teuchos::ParameterList> ExplicitRKStepper<Scalar>::getNonconstParameterList()
521{
522 return(parameterList_);
523}
524
525template <class Scalar>
526Teuchos::RCP<Teuchos::ParameterList> ExplicitRKStepper<Scalar>::unsetParameterList()
527{
528 Teuchos::RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
529 parameterList_ = Teuchos::null;
530 return(temp_param_list);
531}
532
533template<class Scalar>
534RCP<const Teuchos::ParameterList>
536{
537 using Teuchos::ParameterList;
538 static RCP<const ParameterList> validPL;
539 if (is_null(validPL)) {
540 RCP<ParameterList> pl = Teuchos::parameterList();
541 Teuchos::setupVerboseObjectSublist(&*pl);
542 validPL = pl;
543 }
544 return validPL;
545}
546
547template<class Scalar>
548void ExplicitRKStepper<Scalar>::setModel(const RCP<const Thyra::ModelEvaluator<Scalar> >& model)
549{
550 TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
551 TEUCHOS_TEST_FOR_EXCEPT( !is_null(model_) ); // For now you can only call this once.
552 assertValidModel( *this, *model );
553 model_ = model;
554}
555
556
557template<class Scalar>
558void ExplicitRKStepper<Scalar>::setNonconstModel(const RCP<Thyra::ModelEvaluator<Scalar> >& model)
559{
560 this->setModel(model); // TODO 09/09/09 tscoffe: use ConstNonconstObjectContainer!
561}
562
563
564template<class Scalar>
565Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
567{
568 return model_;
569}
570
571
572template<class Scalar>
573Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >
575{
576 return Teuchos::null;
577}
578
579
580template<class Scalar>
582 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
583 )
584{
585
586 // typedef Thyra::ModelEvaluatorBase MEB; // unused
587
588 basePoint_ = initialCondition;
589
590 // x
591
592 RCP<const Thyra::VectorBase<Scalar> >
593 x_init = initialCondition.get_x();
594
595#ifdef HAVE_RYTHMOS_DEBUG
596 TEUCHOS_TEST_FOR_EXCEPTION(
597 is_null(x_init), std::logic_error,
598 "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
599 "then x can not be null!" );
600#endif
601
602 solution_vector_ = x_init->clone_v();
603 solution_vector_old_ = x_init->clone_v();
604
605 // for the embedded RK method
606 solution_hat_vector_ = x_init->clone_v();
607 ee_ = x_init->clone_v();
608
609 // t
610
611 t_ = initialCondition.get_t();
612 t_old_ = t_;
613
614 haveInitialCondition_ = true;
615
616}
617
618
619template<class Scalar>
620Thyra::ModelEvaluatorBase::InArgs<Scalar>
622{
623 return basePoint_;
624}
625
626template<class Scalar>
628{
629 return true;
630}
631
632template<class Scalar>
634{
635 // Just use the interface to clone the algorithm in a basically
636 // uninitialized state
637 RCP<ExplicitRKStepper<Scalar> >
638 stepper = Teuchos::rcp(new ExplicitRKStepper<Scalar>());
639
640 if (!is_null(model_)) {
641 stepper->setModel(model_); // Shallow copy is okay!
642 }
643
644 if (!is_null(erkButcherTableau_)) {
645 // 06/16/09 tscoffe: should we clone the RKBT here?
646 stepper->setRKButcherTableau(erkButcherTableau_);
647 }
648
649 if (!is_null(parameterList_)) {
650 stepper->setParameterList(Teuchos::parameterList(*parameterList_));
651 }
652
653 return stepper;
654
655}
656
657template<class Scalar>
659{
660 return(stepControl_);
661}
662
663template<class Scalar>
664RCP<const StepControlStrategyBase<Scalar> > ExplicitRKStepper<Scalar>::getStepControlStrategy() const
665{
666 return(stepControl_);
667}
668
669template<class Scalar>
671{
672 TEUCHOS_TEST_FOR_EXCEPTION(stepControl == Teuchos::null,std::logic_error,"Error, stepControl == Teuchos::null!\n");
673 stepControl_ = stepControl;
674}
675//
676// Explicit Instantiation macro
677//
678// Must be expanded from within the Rythmos namespace!
679//
680
681#define RYTHMOS_EXPLICIT_RK_STEPPER_INSTANT(SCALAR) \
682 \
683 template class ExplicitRKStepper< SCALAR >; \
684 \
685 template RCP< ExplicitRKStepper< SCALAR > > \
686 explicitRKStepper(); \
687 \
688 template RCP< ExplicitRKStepper< SCALAR > > \
689 explicitRKStepper( \
690 const RCP<Thyra::ModelEvaluator< SCALAR > >& model \
691 ); \
692 \
693 template RCP< ExplicitRKStepper< SCALAR > > \
694 explicitRKStepper( \
695 const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \
696 const RCP<const RKButcherTableauBase< SCALAR > >& rkbt \
697 ); \
698
699} // namespace Rythmos
700
701#endif //Rythmos_ExplicitRK_STEPPER_DEF_H
702
void setNonconstModel(const Teuchos::RCP< Thyra::ModelEvaluator< Scalar > > &model)
void setStepControlStrategy(const RCP< StepControlStrategyBase< Scalar > > &stepControlStrategy)
void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
const StepStatus< Scalar > getStepStatus() const
void getNodes(Array< Scalar > *time_vec) const
Get interpolation nodes.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
RCP< const StepControlStrategyBase< Scalar > > getStepControlStrategy() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
void removeNodes(Array< Scalar > &time_vec)
Remove interpolation nodes.
TimeRange< Scalar > getTimeRange() const
void getPoints(const Array< Scalar > &time_vec, Array< RCP< const VectorBase< Scalar > > > *x_vec, Array< RCP< const VectorBase< Scalar > > > *xdot_vec, Array< ScalarMag > *accuracy_vec) const
Get values from buffer.
void setRKButcherTableau(const RCP< const RKButcherTableauBase< Scalar > > &rkbt)
RCP< StepperBase< Scalar > > cloneStepperAlgorithm() const
RCP< Thyra::ModelEvaluator< Scalar > > getNonconstModel()
Scalar takeStep(Scalar dt, StepSizeType flag)
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
int getOrder() const
Get order of interpolation.
void setInitialCondition(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &initialCondition)
void addPoints(const Array< Scalar > &time_vec, const Array< Teuchos::RCP< const Thyra::VectorBase< Scalar > > > &x_vec, const Array< Teuchos::RCP< const Thyra::VectorBase< Scalar > > > &xdot_vec)
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
Redefined from Teuchos::ParameterListAcceptor.
Thyra::ModelEvaluatorBase::InArgs< Scalar > getInitialCondition() const
RCP< const RKButcherTableauBase< Scalar > > getRKButcherTableau() const
RCP< const Teuchos::ParameterList > getValidParameters() const
RCP< StepControlStrategyBase< Scalar > > getNonconstStepControlStrategy()
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
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