Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_StepperValidator.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_STEPPER_VALIDATOR_H
30#define Rythmos_STEPPER_VALIDATOR_H
31
32#include "Teuchos_ParameterList.hpp"
33#include "Teuchos_ParameterListAcceptor.hpp"
34#include "Teuchos_Describable.hpp"
35#include "Teuchos_VerboseObject.hpp"
36#include "Rythmos_IntegratorBuilder.hpp"
37#include "Thyra_StateFuncModelEvaluatorBase.hpp"
38#include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
39#include "Thyra_ModelEvaluator.hpp"
40
41#include "Rythmos_StepperBase.hpp"
42#include "Rythmos_Types.hpp"
43#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
44#include "Thyra_DetachedVectorView.hpp"
45#include "Thyra_DefaultSpmdVectorSpace.hpp"
46#include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp"
47#include "Rythmos_TimeStepNonlinearSolver.hpp"
48
49#include "Teuchos_StandardCatchMacros.hpp"
50
51
52namespace Rythmos {
53
54
62template<class Scalar>
64 : virtual public Teuchos::Describable
65 , virtual public Teuchos::ParameterListAcceptor
66 , virtual public Teuchos::VerboseObject<StepperValidator<Scalar> >
67{
68public:
69
71
72 virtual ~StepperValidator();
73
76 const RCP<IntegratorBuilder<Scalar> > &integratorBuilder
77 );
78
80 void validateStepper() const;
81
84
86 void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
87
89 RCP<Teuchos::ParameterList> getNonconstParameterList();
90
92 RCP<Teuchos::ParameterList> unsetParameterList();
93
95 RCP<const Teuchos::ParameterList> getValidParameters() const;
96
98
101
103 std::string description() const;
104
106 void describe(
107 Teuchos::FancyOStream &out,
108 const Teuchos::EVerbosityLevel verbLevel
109 ) const;
110
112
113private:
114
115 // Private member functions:
116 void defaultInitializeAll_();
117 RCP<StepperBase<Scalar> > getStepper_(
118 const RCP<Thyra::ModelEvaluator<Scalar> >& model,
119 const Thyra::ModelEvaluatorBase::InArgs<Scalar>& initialCondition,
120 const RCP<Thyra::NonlinearSolverBase<Scalar> >& nlSolver
121 ) const;
122 bool isImplicitStepper_() const;
123 Thyra::ModelEvaluatorBase::InArgs<Scalar> getSomeIC_(
124 const Thyra::ModelEvaluator<Scalar>& model) const;
125
126 // Validate that parameters set through setInitialCondition actually get set
127 // on the model when evalModel is called.
128 void validateIC_() const;
129 // Validate that getTimeRange and getStepStatus return the correct states of
130 // initialization.
131 void validateStates_() const;
132 // Validate that we can get the initial condition through getPoints after
133 // setInitialCondition has been set and after the first step.
134 void validateGetIC_() const;
135 // Validate that we can get the initial condition through
136 // getInitialCondition
137 void validateGetIC2_() const;
138 // Validate that the stepper supports getNodes, which is used by the
139 // Trailing Interpolation Buffer feature of the Integrator.
140 void validateGetNodes_() const;
141
142 // Private member data:
143 RCP<IntegratorBuilder<Scalar> > integratorBuilder_;
144 RCP<ParameterList> paramList_;
145 mutable std::string stepperName_;
146
147};
148
149//
150// StepperValidator nonmember constructor:
151//
152template<class Scalar>
153RCP<StepperValidator<Scalar> > stepperValidator()
154{
155 return rcp(new StepperValidator<Scalar>() );
156}
157
158//
159// Mock Model class for validating a stepper
160//
161template<class Scalar>
162class StepperValidatorMockModel
163 : public Thyra::StateFuncModelEvaluatorBase<Scalar>,
164 public Teuchos::ParameterListAcceptorDefaultBase
165{
166public:
167
168 StepperValidatorMockModel();
169
170 virtual ~StepperValidatorMockModel();
171
172 RCP<const std::vector<Thyra::ModelEvaluatorBase::InArgs<Scalar> > > getPassedInArgs() const;
173
176
178 RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
180 RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
182 Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
184 RCP<Thyra::LinearOpWithSolveBase<Scalar> > create_W() const;
186 RCP<Thyra::LinearOpBase<Scalar> > create_W_op() const;
188 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > get_W_factory() const;
190 Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
191
193 RCP<const Thyra::VectorSpaceBase<Scalar> > get_p_space(int l) const;
195 RCP<const Teuchos::Array<std::string> > get_p_names(int l) const;
197 RCP<const Thyra::VectorSpaceBase<Scalar> > get_g_space(int j) const;
198
200
203
205 void setParameterList(RCP<ParameterList> const& paramList);
206
208 RCP<const ParameterList> getValidParameters() const;
209
211
212private:
213
216
218 Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
219
221 void evalModelImpl(
222 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs_bar,
223 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs_bar
224 ) const;
225
227
228 void defaultInitializeAll_();
229 void initialize_();
230
231 mutable RCP<std::vector<Thyra::ModelEvaluatorBase::InArgs<Scalar> > > passedInArgs_;
232 bool isImplicit_;
233 bool isInitialized_;
234 Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs_;
235 Thyra::ModelEvaluatorBase::OutArgs<Scalar> outArgs_;
236 Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
237 RCP<const Thyra::VectorSpaceBase<Scalar> > x_space_;
238 RCP<const Thyra::VectorSpaceBase<Scalar> > f_space_;
239 RCP<const Thyra::VectorSpaceBase<Scalar> > p_space_;
240 int Np_; // number of parameter vectors (1)
241 int Ng_; // number of observation functions (0)
242 int np_; // length of parameter vector (1)
243 int dim_; // length of solution vector
244
245};
246
247//
248// StepperValidatorMockModel nonmember constructors:
249//
250template<class Scalar>
251RCP<StepperValidatorMockModel<Scalar> > stepperValidatorMockModel()
252{
253 return(rcp(new StepperValidatorMockModel<Scalar>()));
254}
255
256template<class Scalar>
257RCP<StepperValidatorMockModel<Scalar> > stepperValidatorMockModel(
258 bool isImplicit
259 )
260{
261 RCP<StepperValidatorMockModel<Scalar> > model = rcp(new StepperValidatorMockModel<Scalar>());
262 RCP<ParameterList> pl = Teuchos::parameterList();
263 pl->set("Is Implicit",isImplicit);
264 model->setParameterList(pl);
265 return model;
266}
267
268//
269// StepperValidatorMockModel Definitions:
270//
271template<class Scalar>
272StepperValidatorMockModel<Scalar>::StepperValidatorMockModel()
273{
274 this->defaultInitializeAll_();
275 Np_ = 1;
276 Ng_ = 0;
277 np_ = 1;
278 dim_ = 1;
279 // Create x_space and f_space
280 x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
281 f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
282 // Create p_space
283 p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(np_);
284 passedInArgs_ = rcp(new std::vector<Thyra::ModelEvaluatorBase::InArgs<Scalar> >);
285 this->initialize_();
286}
287
288template<class Scalar>
289void StepperValidatorMockModel<Scalar>::defaultInitializeAll_()
290{
291 passedInArgs_ = Teuchos::null;
292 isImplicit_ = false;
293 isInitialized_ = false;
294 //inArgs_;
295 //outArgs_;
296 //nominalValues_;
297 x_space_ = Teuchos::null;
298 f_space_ = Teuchos::null;
299 p_space_ = Teuchos::null;
300 Np_ = -1;
301 Ng_ = -1;
302 np_ = -1;
303 dim_ = -1;
304}
305
306template<class Scalar>
307StepperValidatorMockModel<Scalar>::~StepperValidatorMockModel()
308{ }
309
310template<class Scalar>
311void StepperValidatorMockModel<Scalar>::initialize_()
312{
313 if (!isInitialized_) {
314 // Set up prototypical InArgs
315 Thyra::ModelEvaluatorBase::InArgsSetup<Scalar> inArgs;
316 inArgs.setModelEvalDescription(this->description());
317 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_t );
318 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x );
319 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_beta );
320#ifdef HAVE_THYRA_ME_POLYNOMIAL
321 // For ExplicitTaylorPolynomialStepper
322 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x_poly );
323#endif // HAVE_THYRA_ME_POLYNOMIAL
324 inArgs.set_Np(1);
325 if (isImplicit_) {
326 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x_dot );
327 inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_alpha );
328 }
329 inArgs_ = inArgs;
330 // Set up prototypical OutArgs
331 Thyra::ModelEvaluatorBase::OutArgsSetup<Scalar> outArgs;
332 outArgs.setModelEvalDescription(this->description());
333 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_f );
334 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_W_op );
335#ifdef HAVE_THYRA_ME_POLYNOMIAL
336 // For ExplicitTaylorPolynomialStepper
337 outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_f_poly );
338#endif // HAVE_THYRA_ME_POLYNOMIAL
339 outArgs.set_Np_Ng(Np_,Ng_);
340 outArgs_ = outArgs;
341 // Set up nominal values
342 nominalValues_ = inArgs_;
343 nominalValues_.set_t(Scalar(1.0));
344 const RCP<VectorBase<Scalar> > x_ic = Thyra::createMember(x_space_);
345 Thyra::V_S(Teuchos::outArg(*x_ic),Scalar(2.0));
346 nominalValues_.set_x(x_ic);
347 const RCP<VectorBase<Scalar> > p_ic = Thyra::createMember(p_space_);
348 Thyra::V_S(Teuchos::outArg(*p_ic),Scalar(3.0));
349 nominalValues_.set_p(0,p_ic);
350 isInitialized_ = true;
351 }
352}
353
354
355template<class Scalar>
356RCP<const std::vector<Thyra::ModelEvaluatorBase::InArgs<Scalar> > > StepperValidatorMockModel<Scalar>::getPassedInArgs() const
357{
358 return passedInArgs_;
359}
360
361
362template<class Scalar>
363RCP<const Thyra::VectorSpaceBase<Scalar> > StepperValidatorMockModel<Scalar>::get_x_space() const
364{
365 TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
366 return x_space_;
367}
368
369
370template<class Scalar>
371RCP<const Thyra::VectorSpaceBase<Scalar> > StepperValidatorMockModel<Scalar>::get_f_space() const
372{
373 TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
374 return f_space_;
375}
376
377
378template<class Scalar>
379Thyra::ModelEvaluatorBase::InArgs<Scalar> StepperValidatorMockModel<Scalar>::getNominalValues() const
380{
381 TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
382 return nominalValues_;
383}
384
385
386template<class Scalar>
387RCP<Thyra::LinearOpWithSolveBase<Scalar> > StepperValidatorMockModel<Scalar>::create_W() const
388{
389 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory = this->get_W_factory();
390 RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
391 {
392 // 01/20/09 tscoffe: This is a total hack to provide a full rank matrix to
393 // linearOpWithSolve because it ends up factoring the matrix during
394 // initialization, which it really shouldn't do, or I'm doing something
395 // wrong here. The net effect is that I get exceptions thrown in
396 // optimized mode due to the matrix being rank deficient unless I do this.
397 RCP<Thyra::MultiVectorBase<Scalar> > multivec = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,true);
398 {
399 RCP<Thyra::VectorBase<Scalar> > vec = Thyra::createMember(x_space_);
400 {
401 Thyra::DetachedVectorView<Scalar> vec_view( *vec );
402 vec_view[0] = 1.0;
403 }
404 V_V(multivec->col(0).ptr(),*vec);
405 }
406 }
407 RCP<Thyra::LinearOpWithSolveBase<Scalar> > W =
408 Thyra::linearOpWithSolve<Scalar>(
409 *W_factory,
410 matrix
411 );
412 return W;
413}
414
415
416template<class Scalar>
417RCP<Thyra::LinearOpBase<Scalar> > StepperValidatorMockModel<Scalar>::create_W_op() const
418{
419 RCP<Thyra::MultiVectorBase<Scalar> > matrix = Thyra::createMembers(x_space_, dim_);
420 return(matrix);
421}
422
423
424template<class Scalar>
425RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > StepperValidatorMockModel<Scalar>::get_W_factory() const
426{
427 RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory =
428 Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
429 return W_factory;
430}
431
432
433template<class Scalar>
434Thyra::ModelEvaluatorBase::InArgs<Scalar> StepperValidatorMockModel<Scalar>::createInArgs() const
435{
436 TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
437 return inArgs_;
438}
439
440
441template<class Scalar>
442RCP<const Thyra::VectorSpaceBase<Scalar> > StepperValidatorMockModel<Scalar>::get_p_space(int l) const
443{
444 TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
445 return p_space_;
446}
447
448
449template<class Scalar>
450RCP<const Teuchos::Array<std::string> > StepperValidatorMockModel<Scalar>::get_p_names(int l) const
451{
452 RCP<Teuchos::Array<std::string> > p_strings =
453 Teuchos::rcp(new Teuchos::Array<std::string>());
454 p_strings->push_back("Model Coefficient");
455 return p_strings;
456}
457
458
459template<class Scalar>
460RCP<const Thyra::VectorSpaceBase<Scalar> > StepperValidatorMockModel<Scalar>::get_g_space(int j) const
461{
462 return Teuchos::null;
463}
464
465
466template<class Scalar>
467void StepperValidatorMockModel<Scalar>::setParameterList(RCP<ParameterList> const& paramList)
468{
469 TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
470 paramList->validateParameters(*this->getValidParameters());
471 Teuchos::readVerboseObjectSublist(&*paramList,this);
472 bool test_isImplicit = paramList->get("Is Implicit",false);
473 if (isImplicit_ != test_isImplicit) {
474 isImplicit_ = test_isImplicit;
475 isInitialized_ = false;
476 this->initialize_();
477 }
478 setMyParamList(paramList);
479}
480template<class Scalar>
481RCP<const ParameterList> StepperValidatorMockModel<Scalar>::getValidParameters() const
482{
483 static RCP<const ParameterList> validPL;
484 if (is_null(validPL)) {
485 RCP<ParameterList> pl = Teuchos::parameterList();
486 pl->set("Is Implicit",false);
487 Teuchos::setupVerboseObjectSublist(&*pl);
488 validPL = pl;
489 }
490 return validPL;
491}
492template<class Scalar>
493Thyra::ModelEvaluatorBase::OutArgs<Scalar> StepperValidatorMockModel<Scalar>::createOutArgsImpl() const
494{
495 TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
496 return outArgs_;
497}
498
499template<class Scalar>
500void StepperValidatorMockModel<Scalar>::evalModelImpl(
501 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
502 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
503 ) const
504{
505 typedef Teuchos::ScalarTraits<Scalar> ST;
506 passedInArgs_->push_back(inArgs);
507 // Fill f with zeros.
508 RCP<VectorBase<Scalar> > f_out = outArgs.get_f();
509 if (!is_null(f_out)) {
510 Thyra::V_S(Teuchos::outArg(*f_out),ST::zero());
511 }
512#ifdef HAVE_THYRA_ME_POLYNOMIAL
513 if (outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_f_poly)) {
514 RCP<Teuchos::Polynomial<VectorBase<Scalar> > > f_poly_out = outArgs.get_f_poly();
515 if (!is_null(f_poly_out)) {
516 //Thyra::V_S(Teuchos::outArg(*f_poly_out),ST::zero());
517 }
518 }
519#endif // HAVE_THYRA_ME_POLYNOMIAL
520
521}
522
523
524//
525// StepperValidator Definitions:
526//
527
528template<class Scalar>
529StepperValidator<Scalar>::StepperValidator()
530{
531 this->defaultInitializeAll_();
532}
533
534template<class Scalar>
535void StepperValidator<Scalar>::defaultInitializeAll_()
536{
537 integratorBuilder_ = Teuchos::null;
538 paramList_ = Teuchos::null;
539 stepperName_ = "";
540}
541
542template<class Scalar>
543 StepperValidator<Scalar>::~StepperValidator()
544{ }
545
546template<class Scalar>
548 const RCP<IntegratorBuilder<Scalar> > &integratorBuilder
549 )
550{
551 TEUCHOS_ASSERT( !is_null(integratorBuilder) );
552 integratorBuilder_ = integratorBuilder;
553}
554
555template<class Scalar>
557{
558 // Extract the name of the stepper for later
559 {
560 RCP<const ParameterList> pl = integratorBuilder_->getParameterList();
561 if (is_null(pl)) {
562 pl = integratorBuilder_->getValidParameters();
563 }
564 const Teuchos::ParameterList& stepperSelectionPL = pl->sublist("Stepper Settings").sublist("Stepper Selection");
565 stepperName_ = stepperSelectionPL.get<std::string>("Stepper Type");
566 }
567 bool verbose = true;
568 Array<bool> success_array;
569 bool local_success = true;
570
571 // Verify that the stepper passes parameters to the model in evalModel:
572 local_success = true;
573 try {
574 this->validateIC_();
575 }
576 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,local_success);
577 success_array.push_back(local_success);
578
579 // Verify that the stepper states are correct
580 // uninitialized => getTimeRange == invalidTimeRange
581 // initialized, but no step => getTimeRange.length() == 0, [t_ic,t_ic]
582 // initialized, step taken => getTimeRange.length() > 0
583 local_success = true;
584 try {
585 this->validateStates_();
586 }
587 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,local_success);
588 success_array.push_back(local_success);
589
590 // Verify that getPoints(t_ic) returns the IC after initialization
591 // and after the first step
592 local_success = true;
593 try {
594 this->validateGetIC_();
595 }
596 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,local_success);
597 success_array.push_back(local_success);
598
599 // Verify that getInitialCondition() returns the IC after
600 // setInitialCondition(...)
601 local_success = true;
602 try {
603 this->validateGetIC2_();
604 }
605 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,local_success);
606 success_array.push_back(local_success);
607
608 // Validate that the stepper supports getNodes, which is used by
609 // the Trailing Interpolation Buffer feature of the Integrator.
610 local_success = true;
611 try {
612 this->validateGetNodes_();
613 }
614 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,local_success);
615 success_array.push_back(local_success);
616
617 // Verify that getPoints(t) returns the same vectors as getStepStatus
618 // TODO
619
620 bool global_success = true;
621 for (int i=0 ; i < Teuchos::as<int>(success_array.size()) ; ++i) {
622 global_success = global_success && success_array[i];
623 }
624
625 TEUCHOS_TEST_FOR_EXCEPTION( !global_success, std::logic_error,
626 "Error! StepperValidator: The stepper " << stepperName_
627 << " did not pass stepper validation.");
628}
629
630template<class Scalar>
631 void StepperValidator<Scalar>::setParameterList(RCP<Teuchos::ParameterList> const& paramList)
632{
633 TEUCHOS_TEST_FOR_EXCEPT(true);
634}
635
636template<class Scalar>
638{
639 TEUCHOS_TEST_FOR_EXCEPT(true);
640 return Teuchos::null;
641}
642
643template<class Scalar>
645{
646 TEUCHOS_TEST_FOR_EXCEPT(true);
647 return Teuchos::null;
648}
649
650template<class Scalar>
651 RCP<const Teuchos::ParameterList> StepperValidator<Scalar>::getValidParameters() const
652{
653 TEUCHOS_TEST_FOR_EXCEPT(true);
654 return Teuchos::null;
655}
656
657template<class Scalar>
659{
660 TEUCHOS_TEST_FOR_EXCEPT(true);
661 return("");
662}
663
664template<class Scalar>
666 Teuchos::FancyOStream &out,
667 const Teuchos::EVerbosityLevel verbLevel
668 ) const
669{
670 TEUCHOS_TEST_FOR_EXCEPT(true);
671}
672
673template<class Scalar>
674RCP<StepperBase<Scalar> > StepperValidator<Scalar>::getStepper_(
675 const RCP<Thyra::ModelEvaluator<Scalar> >& model,
676 const Thyra::ModelEvaluatorBase::InArgs<Scalar>& initialCondition,
677 const RCP<Thyra::NonlinearSolverBase<Scalar> >& nlSolver
678 ) const
679{
680 RCP<IntegratorBase<Scalar> > integrator = integratorBuilder_->create(model,initialCondition,nlSolver);
681 RCP<StepperBase<Scalar> > stepper = integrator->getNonconstStepper();
682 return stepper;
683}
684
685template<class Scalar>
686bool StepperValidator<Scalar>::isImplicitStepper_() const
687{
688 RCP<const StepperBuilder<Scalar> > sb = integratorBuilder_->getStepperBuilder();
689 RCP<StepperBase<Scalar> > stepper = sb->create(stepperName_);
690 return stepper->isImplicit();
691}
692
693template<class Scalar>
694Thyra::ModelEvaluatorBase::InArgs<Scalar> StepperValidator<Scalar>::getSomeIC_(
695 const Thyra::ModelEvaluator<Scalar>& model
696 ) const
697{
698 // Set up some initial condition:
699 Thyra::ModelEvaluatorBase::InArgs<Scalar> model_ic = model.createInArgs();
700 if (model_ic.supports(Thyra::ModelEvaluatorBase::IN_ARG_t)) {
701 Scalar time = Scalar(0.125);
702 model_ic.set_t(time);
703 }
704 if (model_ic.supports(Thyra::ModelEvaluatorBase::IN_ARG_x)) {
705 RCP<VectorBase<Scalar> > x_ic = Thyra::createMember(*(model.get_x_space()));
706 {
707 Thyra::DetachedVectorView<Scalar> x_ic_view( *x_ic );
708 x_ic_view[0] = Scalar(10.0);
709 }
710 model_ic.set_x(x_ic);
711 }
712 for (int i=0 ; i<model_ic.Np() ; ++i) {
713 RCP<VectorBase<Scalar> > p_ic = Thyra::createMember(*(model.get_p_space(i)));
714 {
715 Thyra::DetachedVectorView<Scalar> p_ic_view( *p_ic );
716 p_ic_view[i] = Scalar(11.0+i);
717 }
718 model_ic.set_p(i,p_ic);
719 }
720 if (model_ic.supports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot)) {
721 RCP<VectorBase<Scalar> > xdot_ic = Thyra::createMember(*(model.get_x_space()));
722 {
723 Thyra::DetachedVectorView<Scalar> xdot_ic_view( *xdot_ic );
724 xdot_ic_view[0] = Scalar(12.0);
725 }
726 model_ic.set_x_dot(xdot_ic);
727 }
728 return model_ic;
729}
730
731
732template<class Scalar>
733void StepperValidator<Scalar>::validateIC_() const
734{
735 // typedef Teuchos::ScalarTraits<Scalar> ST; // unused
736 // Determine if the stepper is implicit or not:
737 bool isImplicit = this->isImplicitStepper_();
738 RCP<StepperValidatorMockModel<Scalar> > model =
739 stepperValidatorMockModel<Scalar>(isImplicit);
740 // Set up some initial condition:
741 Thyra::ModelEvaluatorBase::InArgs<Scalar> stepper_ic = this->getSomeIC_(*model);
742 // Create nonlinear solver (if needed)
743 RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver;
744 if (isImplicit) {
745 nlSolver = Rythmos::timeStepNonlinearSolver<Scalar>();
746 }
747 RCP<StepperBase<Scalar> > stepper = this->getStepper_(model,stepper_ic,nlSolver);
748 Scalar dt = Scalar(0.1);
749 stepper->takeStep(dt,STEP_TYPE_FIXED);
750 // Verify that the parameter got set on the model by asking the model for the
751 // inArgs passed to it through evalModel
752 RCP<const std::vector<Thyra::ModelEvaluatorBase::InArgs<Scalar> > >
753 passedInArgs_ptr = model->getPassedInArgs();
754 const std::vector<Thyra::ModelEvaluatorBase::InArgs<Scalar> >&
755 passedInArgs = *passedInArgs_ptr;
756 bool valid_t = false;
757 // Technically this will fail for any Runge Kutta Butcher Tableau where:
758 // c(0) != 0 and c(s) != 1, where s = # of stages.
759 if ( (passedInArgs[0].get_t() == stepper_ic.get_t() ) ||
760 (passedInArgs[0].get_t() == stepper_ic.get_t()+dt) )
761 {
762 valid_t = true;
763 }
764 TEUCHOS_TEST_FOR_EXCEPTION( !valid_t, std::logic_error,
765 "Error! StepperValidator::validateIC: Time did not get correctly set on"
766 " the model through StepperBase::setInitialCondition!"
767 );
768 // If a stepper uses a predictor, then the x passed into the model will not
769 // be the same as the IC, so we can't check it.
770 RCP<const VectorBase<Scalar> > p_out = passedInArgs[0].get_p(0);
771 TEUCHOS_TEST_FOR_EXCEPTION( is_null(p_out), std::logic_error,
772 "Error! StepperValidator::validateIC: Parameter 0 did not get set on the"
773 " model through StepperBase::setInitialCondition!"
774 );
775 {
776 Thyra::ConstDetachedVectorView<Scalar> p_out_view( *p_out );
777 TEUCHOS_TEST_FOR_EXCEPTION( p_out_view[0] != Scalar(11.0), std::logic_error,
778 "Error! StepperValidator::validateIC: Parameter 0 did not get set correctly on the model through StepperBase::setInitialCondition!"
779 );
780 }
781 if (isImplicit) {
782 // Each time integration method will approximate xdot according to its own
783 // algorithm, so we can't really test it here.
784 }
785}
786
787template<class Scalar>
788void StepperValidator<Scalar>::validateStates_() const
789{
790 // typedef Teuchos::ScalarTraits<Scalar> ST; // unused
791 RCP<const StepperBuilder<Scalar> > sb =
792 integratorBuilder_->getStepperBuilder();
793 RCP<StepperBase<Scalar> > stepper = sb->create(stepperName_);
794 {
795 // Uninitialized Stepper:
796 TimeRange<Scalar> tr = stepper->getTimeRange();
797 TEUCHOS_TEST_FOR_EXCEPTION( tr.isValid(), std::logic_error,
798 "Error! StepperValidator::validateStates: Uninitialized "
799 "stepper returned a valid time range!"
800 );
801 const StepStatus<Scalar> ss = stepper->getStepStatus();
802 TEUCHOS_TEST_FOR_EXCEPTION( ss.stepStatus != STEP_STATUS_UNINITIALIZED,
803 std::logic_error,
804 "Error! StepperValidator::validateStates: Uninitialized "
805 "stepper returned a valid step status!"
806 );
807 }
808 bool isImplicit = stepper->isImplicit();
809 RCP<StepperValidatorMockModel<Scalar> > model =
810 stepperValidatorMockModel<Scalar>(isImplicit);
811 Thyra::ModelEvaluatorBase::InArgs<Scalar> stepper_ic =
812 this->getSomeIC_(*model);
813 stepper->setInitialCondition(stepper_ic);
814 {
815 // Has initial condition:
816 TimeRange<Scalar> tr = stepper->getTimeRange();
817 TEUCHOS_TEST_FOR_EXCEPTION( compareTimeValues(tr.lower(),tr.upper()) != 0,
818 std::logic_error,
819 "Error! StepperValidator::validateStates: Stepper with "
820 "initial condition returned a non zero time range!"
821 );
822// const StepStatus<Scalar> ss = stepper->getStepStatus();
823// TEUCHOS_TEST_FOR_EXCEPTION( ss.stepStatus != STEP_STATUS_UNKNOWN,
824// std::logic_error,
825// "Error! StepperValidator::validateStates: Stepper with "
826// "initial condition did not return STEP_STATUS_UNKNOWN!"
827// );
828 }
829 // 04/16/09 tscoffe: Now we use the integratorBuilder to create a fully
830 // initialized stepper which we can use for taking a step. We can't just
831 // continue setting them up because we don't know what they all require.
832 RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver;
833 if (isImplicit) {
834 nlSolver = timeStepNonlinearSolver<Scalar>();
835 }
836 stepper = this->getStepper_(model,stepper_ic,nlSolver);
837 {
838 // Still has initial condition:
839 TimeRange<Scalar> tr = stepper->getTimeRange();
840 TEUCHOS_TEST_FOR_EXCEPTION( compareTimeValues(tr.lower(),tr.upper()) != 0,
841 std::logic_error,
842 "Error! StepperValidator::validateStates: Fully initialized "
843 "stepper returned a non zero time range!"
844 );
845// const StepStatus<Scalar> ss = stepper->getStepStatus();
846// TEUCHOS_TEST_FOR_EXCEPTION( ss.stepStatus != STEP_STATUS_UNKNOWN,
847// std::logic_error,
848// "Error! StepperValidator::validateStates: Fully initialized "
849// "stepper, prior to taking a step, did not return STEP_STATUS_UNKNOWN!"
850// );
851 }
852 Scalar dt = Scalar(0.1);
853 stepper->takeStep(dt,STEP_TYPE_FIXED);
854 {
855 // Taken Step:
856 TimeRange<Scalar> tr = stepper->getTimeRange();
857 TEUCHOS_TEST_FOR_EXCEPTION( compareTimeValues(tr.lower(),tr.upper()) >= 0,
858 std::logic_error,
859 "Error! StepperValidator::validateStates: Stepper returned "
860 "a zero (or invalid) time range after taking a step!"
861 );
862 const StepStatus<Scalar> ss = stepper->getStepStatus();
863 TEUCHOS_TEST_FOR_EXCEPTION( ss.stepStatus != STEP_STATUS_CONVERGED,
864 std::logic_error,
865 "Error! StepperValidator::validateStates: Stepper did not "
866 "return converged step status after taking a step!"
867 );
868 }
869}
870
871template<class Scalar>
872void StepperValidator<Scalar>::validateGetIC_() const
873{
874 // typedef Teuchos::ScalarTraits<Scalar> ST; // unused
875 typedef typename ScalarTraits<Scalar>::magnitudeType ScalarMag;
876 // Determine if the stepper is implicit or not:
877 bool isImplicit = this->isImplicitStepper_();
878 RCP<StepperValidatorMockModel<Scalar> > model =
879 stepperValidatorMockModel<Scalar>(isImplicit);
880 // Set up some initial condition:
881 Thyra::ModelEvaluatorBase::InArgs<Scalar> stepper_ic = this->getSomeIC_(*model);
882 // Create nonlinear solver (if needed)
883 RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver;
884 if (isImplicit) {
885 nlSolver = Rythmos::timeStepNonlinearSolver<Scalar>();
886 }
887 RCP<StepperBase<Scalar> > stepper = this->getStepper_(model,stepper_ic,nlSolver);
888 // Verify we can get the IC through getPoints prior to taking a step:
889 {
890 Array<Scalar> time_vec;
891 Array<RCP<const VectorBase<Scalar> > > x_vec;
892 Array<RCP<const VectorBase<Scalar> > > xdot_vec;
893 Array<ScalarMag> accuracy_vec;
894 time_vec.push_back(stepper_ic.get_t());
895 stepper->getPoints(time_vec,&x_vec,&xdot_vec,&accuracy_vec);
896 {
897 Thyra::ConstDetachedVectorView<Scalar> x_view( *x_vec[0] );
898 TEUCHOS_TEST_FOR_EXCEPTION( compareTimeValues<Scalar>(x_view[0],Scalar(10.0)) != 0,
899 std::logic_error,
900 "Error! StepperValidator::validateGetIC: Stepper did not return the initial"
901 " condition for X through getPoints prior to taking a step!"
902 );
903 }
904 if (isImplicit && !is_null(xdot_vec[0])) {
905 Thyra::ConstDetachedVectorView<Scalar> xdot_view( *xdot_vec[0] );
906 TEUCHOS_TEST_FOR_EXCEPTION( compareTimeValues<Scalar>(xdot_view[0],Scalar(12.0)) != 0,
907 std::logic_error,
908 "Error! StepperValidator::validateGetIC: Stepper did not return the initial"
909 " condition for XDOT through getPoints prior to taking a step!"
910 );
911 }
912 }
913 // Verify we can get the IC through getPoints after taking a step:
914 {
915 Scalar dt = Scalar(0.1);
916 stepper->takeStep(dt,STEP_TYPE_FIXED);
917 Array<Scalar> time_vec;
918 Array<RCP<const VectorBase<Scalar> > > x_vec;
919 Array<RCP<const VectorBase<Scalar> > > xdot_vec;
920 Array<ScalarMag> accuracy_vec;
921 time_vec.push_back(stepper_ic.get_t());
922 stepper->getPoints(time_vec,&x_vec,&xdot_vec,&accuracy_vec);
923 {
924 Thyra::ConstDetachedVectorView<Scalar> x_view( *x_vec[0] );
925 TEUCHOS_TEST_FOR_EXCEPTION( compareTimeValues<Scalar>(x_view[0],Scalar(10.0)) != 0,
926 std::logic_error,
927 "Error! StepperValidator::validateGetIC: Stepper did not return the initial"
928 " condition for X through getPoints after taking a step!"
929 );
930 }
931 if (isImplicit && !is_null(xdot_vec[0])) {
932 Thyra::ConstDetachedVectorView<Scalar> xdot_view( *xdot_vec[0] );
933 TEUCHOS_TEST_FOR_EXCEPTION( compareTimeValues<Scalar>(xdot_view[0],Scalar(12.0)) != 0,
934 std::logic_error,
935 "Error! StepperValidator::validateGetIC: Stepper did not return the initial"
936 " condition for XDOT through getPoints after taking a step!"
937 );
938 }
939 }
940}
941
942
943template<class Scalar>
944void StepperValidator<Scalar>::validateGetIC2_() const
945{
946 // typedef Teuchos::ScalarTraits<Scalar> ST; // unused
947 // typedef typename ScalarTraits<Scalar>::magnitudeType ScalarMag; // unused
948 // Determine if the stepper is implicit or not:
949 bool isImplicit = this->isImplicitStepper_();
950 RCP<StepperValidatorMockModel<Scalar> > model =
951 stepperValidatorMockModel<Scalar>(isImplicit);
952 // Set up some initial condition:
953 Thyra::ModelEvaluatorBase::InArgs<Scalar> stepper_ic = this->getSomeIC_(*model);
954 // Create nonlinear solver (if needed)
955 RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver;
956 if (isImplicit) {
957 nlSolver = Rythmos::timeStepNonlinearSolver<Scalar>();
958 }
959 RCP<StepperBase<Scalar> > stepper = this->getStepper_(model,stepper_ic,nlSolver);
960 // Verify we can get the IC back through getInitialCondition:
961 Thyra::ModelEvaluatorBase::InArgs<Scalar> new_ic = stepper->getInitialCondition();
962 //TEUCHOS_ASSERT( new_ic == stepper_ic );
963 // Verify new_ic == stepper_ic
964 {
965 TEUCHOS_ASSERT( new_ic.get_t() == stepper_ic.get_t() );
966 TEUCHOS_ASSERT( new_ic.get_x() == stepper_ic.get_x() );
967 for (int i=0 ; i<stepper_ic.Np() ; ++i) {
968 TEUCHOS_ASSERT( new_ic.get_p(i) == stepper_ic.get_p(i) );
969 }
970 if (isImplicit) {
971 TEUCHOS_ASSERT( new_ic.get_x_dot() == stepper_ic.get_x_dot() );
972 }
973 }
974}
975
976
977template<class Scalar>
978void StepperValidator<Scalar>::validateGetNodes_() const
979{
980 // typedef Teuchos::ScalarTraits<Scalar> ST; // unused
981 // Create uninitialized stepper and verify we get no nodes back
982 {
983 RCP<const StepperBuilder<Scalar> > sb = integratorBuilder_->getStepperBuilder();
984 RCP<StepperBase<Scalar> > stepper = sb->create(stepperName_);
985 Array<Scalar> nodes;
986 stepper->getNodes(&nodes);
987 TEUCHOS_TEST_FOR_EXCEPTION( nodes.size() != 0, std::logic_error,
988 "Error! StepperValidator::validateGetNodes: Uninitialized stepper returned non-empty node list!"
989 );
990 }
991 // Create fully initialize stepper and verify we get back one node for IC
992 bool isImplicit = this->isImplicitStepper_();
993 RCP<StepperValidatorMockModel<Scalar> > model =
994 stepperValidatorMockModel<Scalar>(isImplicit);
995 Thyra::ModelEvaluatorBase::InArgs<Scalar> stepper_ic = this->getSomeIC_(*model);
996 RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver;
997 if (isImplicit) {
998 nlSolver = Rythmos::timeStepNonlinearSolver<Scalar>();
999 }
1000 RCP<StepperBase<Scalar> > stepper = this->getStepper_(model,stepper_ic,nlSolver);
1001 {
1002 Array<Scalar> nodes;
1003 stepper->getNodes(&nodes);
1004 TEUCHOS_TEST_FOR_EXCEPTION( nodes.size() == 0, std::logic_error,
1005 "Error! StepperValidator::validateGetNodes: Initialized stepper returned empty node list!"
1006 );
1007 TEUCHOS_TEST_FOR_EXCEPTION( nodes.size() > 1, std::logic_error,
1008 "Error! StepperValidator::validateGetNodes: Initialized stepper returned node list with more than one node!"
1009 );
1010 }
1011 // Take a step with the stepper and verify we get back two nodes
1012 Scalar dt = Scalar(0.1);
1013 stepper->takeStep(dt,STEP_TYPE_FIXED);
1014 {
1015 Array<Scalar> nodes;
1016 stepper->getNodes(&nodes);
1017 TEUCHOS_TEST_FOR_EXCEPTION( nodes.size() == 0, std::logic_error,
1018 "Error! StepperValidator::validateGetNodes: After taking a step, stepper returned empty node list!"
1019 );
1020 TEUCHOS_TEST_FOR_EXCEPTION( nodes.size() == 1, std::logic_error,
1021 "Error! StepperValidator::validateGetNodes: After taking a step, stepper returned node list with only one node!"
1022 );
1023 TEUCHOS_TEST_FOR_EXCEPTION( nodes.size() > 2, std::logic_error,
1024 "Error! StepperValidator::validateGetNodes: After taking a step, stepper returned node list with more than two nodes!"
1025 );
1026 }
1027}
1028
1029} // namespace Rythmos
1030
1031#endif //Rythmos_STEPPER_VALIDATOR_H
Class for validating steppers.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
void setParameterList(RCP< Teuchos::ParameterList > const &paramList)
RCP< const Teuchos::ParameterList > getValidParameters() const
RCP< Teuchos::ParameterList > unsetParameterList()
void validateStepper() const
Validate the stepper by the StepperBuilder.
void setIntegratorBuilder(const RCP< IntegratorBuilder< Scalar > > &integratorBuilder)
Set a Integrator Builder object.
RCP< Teuchos::ParameterList > getNonconstParameterList()