Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_ForwardEulerStepper_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_FORWARDEULER_STEPPER_DEF_H
30#define Rythmos_FORWARDEULER_STEPPER_DEF_H
31
32#include "Rythmos_ForwardEulerStepper_decl.hpp"
33#include "Rythmos_StepperHelpers.hpp"
34#include "Thyra_ModelEvaluatorHelpers.hpp"
35#include "Thyra_MultiVectorStdOps.hpp"
36#include "Thyra_VectorStdOps.hpp"
37#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
38
39
40namespace Rythmos {
41
42
43// --------------------------------------------------------------------
44// ForwardEulerStepperMomento definitions:
45// --------------------------------------------------------------------
46
47template<class Scalar>
48ForwardEulerStepperMomento<Scalar>::ForwardEulerStepperMomento()
49{}
50
51template<class Scalar>
52ForwardEulerStepperMomento<Scalar>::~ForwardEulerStepperMomento()
53{}
54
55template<class Scalar>
56void ForwardEulerStepperMomento<Scalar>::serialize(
57 const StateSerializerStrategy<Scalar>& stateSerializer,
58 std::ostream& oStream
59 ) const
60{
61 using Teuchos::is_null;
62 TEUCHOS_ASSERT( !is_null(model_) );
63 RCP<VectorBase<Scalar> > sol_vec = solution_vector_;
64 if (is_null(sol_vec)) {
65 sol_vec = Thyra::createMember(model_->get_x_space());
66 }
67 RCP<VectorBase<Scalar> > res_vec = residual_vector_;
68 if (is_null(res_vec)) {
69 res_vec = Thyra::createMember(model_->get_f_space());
70 }
71 RCP<VectorBase<Scalar> > sol_vec_old = solution_vector_old_;
72 if (is_null(sol_vec_old)) {
73 sol_vec_old = Thyra::createMember(model_->get_x_space());
74 }
75 stateSerializer.serializeVectorBase(*sol_vec,oStream);
76 stateSerializer.serializeVectorBase(*res_vec,oStream);
77 stateSerializer.serializeVectorBase(*sol_vec_old,oStream);
78 stateSerializer.serializeScalar(t_,oStream);
79 stateSerializer.serializeScalar(t_old_,oStream);
80 stateSerializer.serializeScalar(dt_,oStream);
81 stateSerializer.serializeInt(numSteps_,oStream);
82 stateSerializer.serializeBool(isInitialized_,oStream);
83 stateSerializer.serializeBool(haveInitialCondition_,oStream);
84 RCP<ParameterList> pl = parameterList_;
85 if (Teuchos::is_null(pl)) {
86 pl = Teuchos::parameterList();
87 }
88 stateSerializer.serializeParameterList(*pl,oStream);
89}
90
91template<class Scalar>
92void ForwardEulerStepperMomento<Scalar>::deSerialize(
93 const StateSerializerStrategy<Scalar>& stateSerializer,
94 std::istream& iStream
95 )
96{
97 using Teuchos::outArg;
98 using Teuchos::is_null;
99 TEUCHOS_ASSERT( !is_null(model_) );
100 if (is_null(solution_vector_)) {
101 solution_vector_ = Thyra::createMember(*model_->get_x_space());
102 }
103 if (is_null(residual_vector_)) {
104 residual_vector_ = Thyra::createMember(*model_->get_f_space());
105 }
106 if (is_null(solution_vector_old_)) {
107 solution_vector_old_ = Thyra::createMember(*model_->get_x_space());
108 }
109 stateSerializer.deSerializeVectorBase(outArg(*solution_vector_),iStream);
110 stateSerializer.deSerializeVectorBase(outArg(*residual_vector_),iStream);
111 stateSerializer.deSerializeVectorBase(outArg(*solution_vector_old_),iStream);
112 stateSerializer.deSerializeScalar(outArg(t_),iStream);
113 stateSerializer.deSerializeScalar(outArg(t_old_),iStream);
114 stateSerializer.deSerializeScalar(outArg(dt_),iStream);
115 stateSerializer.deSerializeInt(outArg(numSteps_),iStream);
116 stateSerializer.deSerializeBool(outArg(isInitialized_),iStream);
117 stateSerializer.deSerializeBool(outArg(haveInitialCondition_),iStream);
118 if (is_null(parameterList_)) {
119 parameterList_ = Teuchos::parameterList();
120 }
121 stateSerializer.deSerializeParameterList(outArg(*parameterList_),iStream);
122}
123
124template<class Scalar>
125RCP<MomentoBase<Scalar> > ForwardEulerStepperMomento<Scalar>::clone() const
126{
127 RCP<ForwardEulerStepperMomento<Scalar> > m = rcp(new ForwardEulerStepperMomento<Scalar>());
128 m->set_solution_vector(solution_vector_);
129 m->set_residual_vector(residual_vector_);
130 m->set_solution_vector_old(solution_vector_old_);
131 m->set_t(t_);
132 m->set_t_old(t_old_);
133 m->set_dt(dt_);
134 m->set_numSteps(numSteps_);
135 m->set_isInitialized(isInitialized_);
136 m->set_haveInitialCondition(haveInitialCondition_);
137 m->set_parameterList(parameterList_);
138 if (!Teuchos::is_null(this->getMyParamList())) {
139 m->setParameterList(Teuchos::parameterList(*(this->getMyParamList())));
140 }
141 m->set_model(model_);
142 m->set_basePoint(basePoint_);
143 // How do I copy the VerboseObject data?
144 // 07/10/09 tscoffe: Its not set up in Teuchos to do this yet
145 return m;
146}
147
148template<class Scalar>
149void ForwardEulerStepperMomento<Scalar>::set_solution_vector(const RCP<const VectorBase<Scalar> >& solution_vector )
150{
151 solution_vector_ = Teuchos::null;
152 if (!Teuchos::is_null(solution_vector)) {
153 solution_vector_ = solution_vector->clone_v();
154 }
155}
156
157template<class Scalar>
158RCP<VectorBase<Scalar> > ForwardEulerStepperMomento<Scalar>::get_solution_vector() const
159{
160 return solution_vector_;
161}
162
163template<class Scalar>
164void ForwardEulerStepperMomento<Scalar>::set_residual_vector(const RCP<const VectorBase<Scalar> >& residual_vector)
165{
166 residual_vector_ = Teuchos::null;
167 if (!Teuchos::is_null(residual_vector)) {
168 residual_vector_ = residual_vector->clone_v();
169 }
170}
171
172template<class Scalar>
173RCP<VectorBase<Scalar> > ForwardEulerStepperMomento<Scalar>::get_residual_vector() const
174{
175 return residual_vector_;
176}
177
178template<class Scalar>
179void ForwardEulerStepperMomento<Scalar>::set_solution_vector_old(const RCP<const VectorBase<Scalar> >& solution_vector_old )
180{
181 solution_vector_old_ = Teuchos::null;
182 if (!Teuchos::is_null(solution_vector_old)) {
183 solution_vector_old_ = solution_vector_old->clone_v();
184 }
185}
186
187template<class Scalar>
188RCP<VectorBase<Scalar> > ForwardEulerStepperMomento<Scalar>::get_solution_vector_old() const
189{
190 return solution_vector_old_;
191}
192
193template<class Scalar>
194void ForwardEulerStepperMomento<Scalar>::set_t(const Scalar & t)
195{
196 t_ = t;
197}
198
199template<class Scalar>
200Scalar ForwardEulerStepperMomento<Scalar>::get_t() const
201{
202 return t_;
203}
204
205template<class Scalar>
206void ForwardEulerStepperMomento<Scalar>::set_t_old(const Scalar & t_old)
207{
208 t_old_ = t_old;
209}
210
211template<class Scalar>
212Scalar ForwardEulerStepperMomento<Scalar>::get_t_old() const
213{
214 return t_old_;
215}
216
217template<class Scalar>
218void ForwardEulerStepperMomento<Scalar>::set_dt(const Scalar & dt)
219{
220 dt_ = dt;
221}
222
223template<class Scalar>
224Scalar ForwardEulerStepperMomento<Scalar>::get_dt() const
225{
226 return dt_;
227}
228
229template<class Scalar>
230void ForwardEulerStepperMomento<Scalar>::set_numSteps(const int & numSteps)
231{
232 numSteps_ = numSteps;
233}
234
235template<class Scalar>
236int ForwardEulerStepperMomento<Scalar>::get_numSteps() const
237{
238 return numSteps_;
239}
240
241template<class Scalar>
242void ForwardEulerStepperMomento<Scalar>::set_isInitialized(const bool & isInitialized)
243{
244 isInitialized_ = isInitialized;
245}
246
247template<class Scalar>
248bool ForwardEulerStepperMomento<Scalar>::get_isInitialized() const
249{
250 return isInitialized_;
251}
252
253template<class Scalar>
254void ForwardEulerStepperMomento<Scalar>::set_haveInitialCondition(const bool & haveInitialCondition)
255{
256 haveInitialCondition_ = haveInitialCondition;
257}
258
259template<class Scalar>
260bool ForwardEulerStepperMomento<Scalar>::get_haveInitialCondition() const
261{
262 return haveInitialCondition_;
263}
264
265template<class Scalar>
266void ForwardEulerStepperMomento<Scalar>::set_parameterList(const RCP<const ParameterList>& pl)
267{
268 parameterList_ = Teuchos::null;
269 if (!Teuchos::is_null(pl)) {
270 parameterList_ = Teuchos::parameterList(*pl);
271 }
272}
273
274template<class Scalar>
275RCP<ParameterList> ForwardEulerStepperMomento<Scalar>::get_parameterList() const
276{
277 return parameterList_;
278}
279
280template<class Scalar>
281void ForwardEulerStepperMomento<Scalar>::setParameterList(const RCP<ParameterList>& paramList)
282{
283 this->setMyParamList(paramList);
284}
285
286template<class Scalar>
287RCP<const ParameterList> ForwardEulerStepperMomento<Scalar>::getValidParameters() const
288{
289 return Teuchos::null;
290}
291
292template<class Scalar>
293void ForwardEulerStepperMomento<Scalar>::set_model(const RCP<const Thyra::ModelEvaluator<Scalar> >& model)
294{
295 model_ = model;
296}
297
298template<class Scalar>
299RCP<const Thyra::ModelEvaluator<Scalar> > ForwardEulerStepperMomento<Scalar>::get_model() const
300{
301 return model_;
302}
303
304template<class Scalar>
305void ForwardEulerStepperMomento<Scalar>::set_basePoint(const RCP<const Thyra::ModelEvaluatorBase::InArgs<Scalar> >& basePoint)
306{
307 basePoint_ = basePoint;
308}
309
310template<class Scalar>
311RCP<const Thyra::ModelEvaluatorBase::InArgs<Scalar> > ForwardEulerStepperMomento<Scalar>::get_basePoint() const
312{
313 return basePoint_;
314}
315
316// --------------------------------------------------------------------
317// ForwardEulerStepper definitions:
318// --------------------------------------------------------------------
319
320// Nonmember constructor
321template<class Scalar>
322RCP<ForwardEulerStepper<Scalar> > forwardEulerStepper()
323{
324 RCP<ForwardEulerStepper<Scalar> > stepper = rcp(new ForwardEulerStepper<Scalar>());
325 return stepper;
326}
327
328// Nonmember constructor
329template<class Scalar>
330RCP<ForwardEulerStepper<Scalar> > forwardEulerStepper(
331 const RCP<Thyra::ModelEvaluator<Scalar> >& model
332 )
333{
334 RCP<ForwardEulerStepper<Scalar> > stepper = forwardEulerStepper<Scalar>();
335 {
336 RCP<StepperBase<Scalar> > stepperBase =
337 Teuchos::rcp_dynamic_cast<StepperBase<Scalar> >(stepper,true);
338 setStepperModel(Teuchos::inOutArg(*stepperBase),model);
339 }
340 return stepper;
341}
342
343template<class Scalar>
345{
346 this->defaultInitializAll_();
347 dt_ = Teuchos::ScalarTraits<Scalar>::zero();
348 numSteps_ = 0;
349}
350
351template<class Scalar>
353{
354 model_ = Teuchos::null;
355 solution_vector_ = Teuchos::null;
356 residual_vector_ = Teuchos::null;
357 t_ = ST::nan();
358 dt_ = ST::nan();
359 t_old_ = ST::nan();
360 solution_vector_old_ = Teuchos::null;
361 //basePoint_;
362 numSteps_ = -1;
363 haveInitialCondition_ = false;
364 parameterList_ = Teuchos::null;
365 isInitialized_ = false;
366}
367
368template<class Scalar>
369void ForwardEulerStepper<Scalar>::initialize_()
370{
371 if (!isInitialized_) {
372 TEUCHOS_TEST_FOR_EXCEPTION( is_null(model_), std::logic_error,
373 "Error! Please set a model on the stepper.\n"
374 );
375 residual_vector_ = Thyra::createMember(model_->get_f_space());
376 isInitialized_ = true;
377 }
378}
379
380template<class Scalar>
384
385template<class Scalar>
386RCP<const Thyra::VectorSpaceBase<Scalar> > ForwardEulerStepper<Scalar>::get_x_space() const
387{
388 TEUCHOS_TEST_FOR_EXCEPTION(!haveInitialCondition_,std::logic_error,"Error, attempting to call get_x_space before setting an initial condition!\n");
389 return(solution_vector_->space());
390}
391
392template<class Scalar>
393Scalar ForwardEulerStepper<Scalar>::takeStep(Scalar dt, StepSizeType flag)
394{
395 TEUCHOS_TEST_FOR_EXCEPTION( !haveInitialCondition_, std::logic_error,
396 "Error! Attempting to call takeStep before setting an initial condition!\n"
397 );
398 this->initialize_();
399 if (flag == STEP_TYPE_VARIABLE) {
400 // print something out about this method not supporting automatic variable step-size
401 return(-ST::one());
402 }
403 //Thyra::eval_f<Scalar>(*model_,*solution_vector_,t_+dt,&*residual_vector_);
404 eval_model_explicit<Scalar>(*model_,basePoint_,*solution_vector_,t_+dt,Teuchos::outArg(*residual_vector_));
405
406 // solution_vector_old_ = solution_vector_
407 Thyra::V_V(Teuchos::outArg(*solution_vector_old_),*solution_vector_);
408 // solution_vector = solution_vector + dt*residual_vector
409 Thyra::Vp_StV(solution_vector_.ptr(),dt,*residual_vector_);
410 t_old_ = t_;
411 t_ += dt;
412 dt_ = dt;
413 numSteps_++;
414
415 return(dt);
416}
417
418template<class Scalar>
419const StepStatus<Scalar> ForwardEulerStepper<Scalar>::getStepStatus() const
420{
421 StepStatus<Scalar> stepStatus;
422
423 if (!haveInitialCondition_) {
424 stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
425 }
426 else if (numSteps_ == 0) {
427 stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
428 stepStatus.order = this->getOrder();
429 stepStatus.time = t_;
430 stepStatus.solution = solution_vector_;
431 }
432 else {
433 stepStatus.stepStatus = STEP_STATUS_CONVERGED;
434 stepStatus.stepSize = dt_;
435 stepStatus.order = this->getOrder();
436 stepStatus.time = t_;
437 stepStatus.stepLETValue = Scalar(-ST::one());
438 stepStatus.solution = solution_vector_;
439 stepStatus.residual = residual_vector_;
440 }
441
442 return(stepStatus);
443}
444
445template<class Scalar>
447{
448 std::string name = "Rythmos::ForwardEulerStepper";
449 return(name);
450}
451
452template<class Scalar>
454 Teuchos::FancyOStream &out,
455 const Teuchos::EVerbosityLevel verbLevel
456 ) const
457{
458 if ( (static_cast<int>(verbLevel) == static_cast<int>(Teuchos::VERB_DEFAULT) ) ||
459 (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW) )
460 ) {
461 out << description() << "::describe" << std::endl;
462 out << "model = " << model_->description() << std::endl;
463 } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)) {
464 } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM)) {
465 } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_HIGH)) {
466 out << "model = " << std::endl;
467 model_->describe(out,verbLevel);
468 out << "solution_vector = " << std::endl;
469 solution_vector_->describe(out,verbLevel);
470 out << "residual_vector = " << std::endl;
471 residual_vector_->describe(out,verbLevel);
472 }
473}
474
475template<class Scalar>
477 const Array<Scalar>& /* time_vec */
478 ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& /* x_vec */
479 ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& /* xdot_vec */
480 )
481{
482 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, addPoints is not implemented for ForwardEulerStepper.\n");
483}
484
485template<class Scalar>
487 const Array<Scalar>& time_vec
488 ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec
489 ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
490 ,Array<ScalarMag>* accuracy_vec) const
491{
492 TEUCHOS_ASSERT( haveInitialCondition_ );
493 using Teuchos::constOptInArg;
494 using Teuchos::null;
495 defaultGetPoints<Scalar>(
496 t_old_,
497 constOptInArg(*solution_vector_old_),
498 Ptr<const VectorBase<Scalar> >(null),
499 t_,
500 constOptInArg(*solution_vector_),
501 Ptr<const VectorBase<Scalar> >(null),
502 time_vec,
503 ptr(x_vec),
504 ptr(xdot_vec),
505 ptr(accuracy_vec),
506 Ptr<InterpolatorBase<Scalar> >(null)
507 );
508}
509
510template<class Scalar>
512{
513 if (!haveInitialCondition_) {
514 return(invalidTimeRange<Scalar>());
515 } else {
516 return(TimeRange<Scalar>(t_old_,t_));
517 }
518}
519
520template<class Scalar>
521void ForwardEulerStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
522{
523 TEUCHOS_ASSERT( time_vec != NULL );
524 time_vec->clear();
525 if (!haveInitialCondition_) {
526 return;
527 } else {
528 time_vec->push_back(t_old_);
529 }
530 if (numSteps_ > 0) {
531 time_vec->push_back(t_);
532 }
533}
534
535template<class Scalar>
536void ForwardEulerStepper<Scalar>::removeNodes(Array<Scalar>& /* time_vec */)
537{
538 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for ForwardEulerStepper.\n");
539}
540
541template<class Scalar>
543{
544 return(1);
545}
546
547template <class Scalar>
548void ForwardEulerStepper<Scalar>::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
549{
550 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
551 paramList->validateParametersAndSetDefaults(*this->getValidParameters());
552 parameterList_ = paramList;
553 Teuchos::readVerboseObjectSublist(&*parameterList_,this);
554}
555
556template <class Scalar>
557Teuchos::RCP<Teuchos::ParameterList> ForwardEulerStepper<Scalar>::getNonconstParameterList()
558{
559 return(parameterList_);
560}
561
562template <class Scalar>
563Teuchos::RCP<Teuchos::ParameterList> ForwardEulerStepper<Scalar>::unsetParameterList()
564{
565 Teuchos::RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
566 parameterList_ = Teuchos::null;
567 return(temp_param_list);
568}
569
570
571template<class Scalar>
572RCP<const Teuchos::ParameterList>
574{
575 using Teuchos::ParameterList;
576 static RCP<const ParameterList> validPL;
577 if (is_null(validPL)) {
578 RCP<ParameterList> pl = Teuchos::parameterList();
579 Teuchos::setupVerboseObjectSublist(&*pl);
580 validPL = pl;
581 }
582 return validPL;
583}
584
585
586template<class Scalar>
587void ForwardEulerStepper<Scalar>::setModel(const RCP<const Thyra::ModelEvaluator<Scalar> >& model)
588{
589 TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
590 assertValidModel( *this, *model );
591 model_ = model;
592}
593
594
595template<class Scalar>
596void ForwardEulerStepper<Scalar>::setNonconstModel(const RCP<Thyra::ModelEvaluator<Scalar> >& model)
597{
598 this->setModel(model); // TODO 09/09/09 tscoffe: use ConstNonconstObjectContainer!
599}
600
601
602template<class Scalar>
603Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
605{
606 return model_;
607}
608
609template<class Scalar>
610Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >
612{
613 return Teuchos::null;
614}
615
616template<class Scalar>
618 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
619 )
620{
621 basePoint_ = initialCondition;
622 t_ = initialCondition.get_t();
623 t_old_ = t_;
624 solution_vector_ = initialCondition.get_x()->clone_v();
625 solution_vector_old_ = solution_vector_->clone_v();
626 haveInitialCondition_ = true;
627}
628
629
630template<class Scalar>
631Thyra::ModelEvaluatorBase::InArgs<Scalar>
633{
634 return basePoint_;
635}
636
637
638template<class Scalar>
640{
641 return true;
642}
643
644template<class Scalar>
645RCP<StepperBase<Scalar> >
647{
648
649 // Just use the interface to clone the algorithm in a basically
650 // uninitialized state
651
652 RCP<StepperBase<Scalar> >
653 stepper = Teuchos::rcp(new ForwardEulerStepper<Scalar>());
654
655 if (!is_null(model_)) {
656 setStepperModel(Teuchos::inOutArg(*stepper),model_); // Shallow copy is okay!
657 }
658
659 if (!is_null(parameterList_)) {
660 stepper->setParameterList(Teuchos::parameterList(*parameterList_));
661 }
662
663 return stepper;
664
665}
666
667template<class Scalar>
668RCP<const MomentoBase<Scalar> >
670{
671 RCP<ForwardEulerStepperMomento<Scalar> > momento = Teuchos::rcp(new ForwardEulerStepperMomento<Scalar>());
672 momento->set_solution_vector(solution_vector_);
673 momento->set_solution_vector_old(solution_vector_old_);
674 momento->set_residual_vector(residual_vector_);
675 momento->set_t(t_);
676 momento->set_t_old(t_old_);
677 momento->set_dt(dt_);
678 momento->set_numSteps(numSteps_);
679 momento->set_isInitialized(isInitialized_);
680 momento->set_haveInitialCondition(haveInitialCondition_);
681 momento->set_parameterList(parameterList_);
682 momento->set_model(model_);
683 RCP<Thyra::ModelEvaluatorBase::InArgs<Scalar> > bp = Teuchos::rcp(new Thyra::ModelEvaluatorBase::InArgs<Scalar>(basePoint_));
684 momento->set_basePoint(bp);
685 return momento;
686}
687
688template<class Scalar>
689void ForwardEulerStepper<Scalar>::setMomento( const Ptr<const MomentoBase<Scalar> >& momentoPtr )
690{
691 Ptr<const ForwardEulerStepperMomento<Scalar> > feMomentoPtr =
692 Teuchos::ptr_dynamic_cast<const ForwardEulerStepperMomento<Scalar> >(momentoPtr,true);
693 TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::is_null(feMomentoPtr->get_model()), std::logic_error,
694 "Error! Rythmos::ForwardEulerStepper::setMomento: The momento must have a valid model through set_model(...) prior to calling ForwardEulerStepper::setMomento(...)."
695 );
696 TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::is_null(feMomentoPtr->get_basePoint()), std::logic_error,
697 "Error! Rythmos::ForwardEulerStepper::setMomento: The momento must have a valid base point through set_basePoint(...) prior to calling ForwardEulerStepper::setMomento(...)."
698 );
699 model_ = feMomentoPtr->get_model();
700 basePoint_ = *(feMomentoPtr->get_basePoint());
701 const ForwardEulerStepperMomento<Scalar>& feMomento = *feMomentoPtr;
702 solution_vector_ = feMomento.get_solution_vector();
703 solution_vector_old_ = feMomento.get_solution_vector_old();
704 residual_vector_ = feMomento.get_residual_vector();
705 t_ = feMomento.get_t();
706 t_old_ = feMomento.get_t_old();
707 dt_ = feMomento.get_dt();
708 numSteps_ = feMomento.get_numSteps();
709 isInitialized_ = feMomento.get_isInitialized();
710 haveInitialCondition_ = feMomento.get_haveInitialCondition();
711 parameterList_ = feMomento.get_parameterList();
712 this->checkConsistentState_();
713}
714
715template<class Scalar>
717{
718 if (isInitialized_) {
719 TEUCHOS_ASSERT( !Teuchos::is_null(model_) );
720 TEUCHOS_ASSERT( !Teuchos::is_null(residual_vector_) );
721 }
722 if (haveInitialCondition_) {
723 TEUCHOS_ASSERT( !ST::isnaninf(t_) );
724 TEUCHOS_ASSERT( !ST::isnaninf(t_old_) );
725 TEUCHOS_ASSERT( !Teuchos::is_null(solution_vector_) );
726 TEUCHOS_ASSERT( !Teuchos::is_null(solution_vector_old_) );
727 TEUCHOS_ASSERT( t_ >= basePoint_.get_t() );
728 TEUCHOS_ASSERT( t_old_ >= basePoint_.get_t() );
729 }
730 if (numSteps_ > 0) {
731 TEUCHOS_ASSERT(isInitialized_);
732 TEUCHOS_ASSERT(haveInitialCondition_);
733 }
734}
735
736//
737// Explicit Instantiation macro
738//
739// Must be expanded from within the Rythmos namespace!
740//
741
742#define RYTHMOS_FORWARD_EULER_STEPPER_INSTANT(SCALAR) \
743 \
744 template class ForwardEulerStepperMomento< SCALAR >; \
745 template class ForwardEulerStepper< SCALAR >; \
746 \
747 template RCP<ForwardEulerStepper< SCALAR > > forwardEulerStepper(); \
748 template RCP<ForwardEulerStepper< SCALAR > > forwardEulerStepper( \
749 const RCP<Thyra::ModelEvaluator< SCALAR > >& model \
750 );
751
752
753} // namespace Rythmos
754
755#endif //Rythmos_FORWARDEULER_STEPPER_DEF_H
RCP< const Thyra::VectorBase< Scalar > > get_x(const InterpolationBufferBase< Scalar > &interpBuffer, const Scalar &t)
Get a single point x(t) from an interpolation buffer.