Rythmos - Transient Integration for Differential Equations  Version of the Day
 All Classes Functions Variables Typedefs Pages
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 
52 namespace Rythmos {
53 
54 
62 template<class Scalar>
63 class StepperValidator
64  : virtual public Teuchos::Describable
65  , virtual public Teuchos::ParameterListAcceptor
66  , virtual public Teuchos::VerboseObject<StepperValidator<Scalar> >
67 {
68 public:
69 
70  StepperValidator();
71 
72  virtual ~StepperValidator();
73 
75  void setIntegratorBuilder(
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 
113 private:
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 //
152 template<class Scalar>
153 RCP<StepperValidator<Scalar> > stepperValidator()
154 {
155  return rcp(new StepperValidator<Scalar>() );
156 }
157 
158 //
159 // Mock Model class for validating a stepper
160 //
161 template<class Scalar>
162 class StepperValidatorMockModel
163  : public Thyra::StateFuncModelEvaluatorBase<Scalar>,
164  public Teuchos::ParameterListAcceptorDefaultBase
165 {
166 public:
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 
212 private:
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 //
250 template<class Scalar>
251 RCP<StepperValidatorMockModel<Scalar> > stepperValidatorMockModel()
252 {
253  return(rcp(new StepperValidatorMockModel<Scalar>()));
254 }
255 
256 template<class Scalar>
257 RCP<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 //
271 template<class Scalar>
272 StepperValidatorMockModel<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 
288 template<class Scalar>
289 void 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 
306 template<class Scalar>
307 StepperValidatorMockModel<Scalar>::~StepperValidatorMockModel()
308 { }
309 
310 template<class Scalar>
311 void 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 
355 template<class Scalar>
356 RCP<const std::vector<Thyra::ModelEvaluatorBase::InArgs<Scalar> > > StepperValidatorMockModel<Scalar>::getPassedInArgs() const
357 {
358  return passedInArgs_;
359 }
360 
361 
362 template<class Scalar>
363 RCP<const Thyra::VectorSpaceBase<Scalar> > StepperValidatorMockModel<Scalar>::get_x_space() const
364 {
365  TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
366  return x_space_;
367 }
368 
369 
370 template<class Scalar>
371 RCP<const Thyra::VectorSpaceBase<Scalar> > StepperValidatorMockModel<Scalar>::get_f_space() const
372 {
373  TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
374  return f_space_;
375 }
376 
377 
378 template<class Scalar>
379 Thyra::ModelEvaluatorBase::InArgs<Scalar> StepperValidatorMockModel<Scalar>::getNominalValues() const
380 {
381  TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
382  return nominalValues_;
383 }
384 
385 
386 template<class Scalar>
387 RCP<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 
416 template<class Scalar>
417 RCP<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 
424 template<class Scalar>
425 RCP<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 
433 template<class Scalar>
434 Thyra::ModelEvaluatorBase::InArgs<Scalar> StepperValidatorMockModel<Scalar>::createInArgs() const
435 {
436  TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
437  return inArgs_;
438 }
439 
440 
441 template<class Scalar>
442 RCP<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 
449 template<class Scalar>
450 RCP<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 
459 template<class Scalar>
460 RCP<const Thyra::VectorSpaceBase<Scalar> > StepperValidatorMockModel<Scalar>::get_g_space(int j) const
461 {
462  return Teuchos::null;
463 }
464 
465 
466 template<class Scalar>
467 void 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 }
480 template<class Scalar>
481 RCP<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 }
492 template<class Scalar>
493 Thyra::ModelEvaluatorBase::OutArgs<Scalar> StepperValidatorMockModel<Scalar>::createOutArgsImpl() const
494 {
495  TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
496  return outArgs_;
497 }
498 
499 template<class Scalar>
500 void 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 
528 template<class Scalar>
529 StepperValidator<Scalar>::StepperValidator()
530 {
531  this->defaultInitializeAll_();
532 }
533 
534 template<class Scalar>
535 void StepperValidator<Scalar>::defaultInitializeAll_()
536 {
537  integratorBuilder_ = Teuchos::null;
538  paramList_ = Teuchos::null;
539  stepperName_ = "";
540 }
541 
542 template<class Scalar>
543  StepperValidator<Scalar>::~StepperValidator()
544 { }
545 
546 template<class Scalar>
547  void StepperValidator<Scalar>::setIntegratorBuilder(
548  const RCP<IntegratorBuilder<Scalar> > &integratorBuilder
549  )
550 {
551  TEUCHOS_ASSERT( !is_null(integratorBuilder) );
552  integratorBuilder_ = integratorBuilder;
553 }
554 
555 template<class Scalar>
556  void StepperValidator<Scalar>::validateStepper() const
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 
630 template<class Scalar>
631  void StepperValidator<Scalar>::setParameterList(RCP<Teuchos::ParameterList> const& paramList)
632 {
633  TEUCHOS_TEST_FOR_EXCEPT(true);
634 }
635 
636 template<class Scalar>
637  RCP<Teuchos::ParameterList> StepperValidator<Scalar>::getNonconstParameterList()
638 {
639  TEUCHOS_TEST_FOR_EXCEPT(true);
640  return Teuchos::null;
641 }
642 
643 template<class Scalar>
644  RCP<Teuchos::ParameterList> StepperValidator<Scalar>::unsetParameterList()
645 {
646  TEUCHOS_TEST_FOR_EXCEPT(true);
647  return Teuchos::null;
648 }
649 
650 template<class Scalar>
651  RCP<const Teuchos::ParameterList> StepperValidator<Scalar>::getValidParameters() const
652 {
653  TEUCHOS_TEST_FOR_EXCEPT(true);
654  return Teuchos::null;
655 }
656 
657 template<class Scalar>
658  std::string StepperValidator<Scalar>::description() const
659 {
660  TEUCHOS_TEST_FOR_EXCEPT(true);
661  return("");
662 }
663 
664 template<class Scalar>
665  void StepperValidator<Scalar>::describe(
666  Teuchos::FancyOStream &out,
667  const Teuchos::EVerbosityLevel verbLevel
668  ) const
669 {
670  TEUCHOS_TEST_FOR_EXCEPT(true);
671 }
672 
673 template<class Scalar>
674 RCP<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 
685 template<class Scalar>
686 bool 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 
693 template<class Scalar>
694 Thyra::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 
732 template<class Scalar>
733 void 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 
787 template<class Scalar>
788 void 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 
871 template<class Scalar>
872 void 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 
943 template<class Scalar>
944 void 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 
977 template<class Scalar>
978 void 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
Rythmos::IntegratorBuilder
Concrete integrator builder class.
Definition: Rythmos_IntegratorBuilder_decl.hpp:63
Rythmos::StepperValidator::unsetParameterList
RCP< Teuchos::ParameterList > unsetParameterList()
Definition: Rythmos_StepperValidator.hpp:644
Rythmos::StepperValidator
Class for validating steppers.
Definition: Rythmos_StepperValidator.hpp:63
Rythmos::StepperValidator::getValidParameters
RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: Rythmos_StepperValidator.hpp:651
Rythmos::StepperValidator::getNonconstParameterList
RCP< Teuchos::ParameterList > getNonconstParameterList()
Definition: Rythmos_StepperValidator.hpp:637
Rythmos::StepperValidator::setParameterList
void setParameterList(RCP< Teuchos::ParameterList > const &paramList)
Definition: Rythmos_StepperValidator.hpp:631
Rythmos::StepperValidator::validateStepper
void validateStepper() const
Validate the stepper by the StepperBuilder.
Definition: Rythmos_StepperValidator.hpp:556
Rythmos::StepperValidator::setIntegratorBuilder
void setIntegratorBuilder(const RCP< IntegratorBuilder< Scalar > > &integratorBuilder)
Set a Integrator Builder object.
Definition: Rythmos_StepperValidator.hpp:547
Rythmos::StepperValidator::describe
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Definition: Rythmos_StepperValidator.hpp:665
Generated on Fri Sep 1 2023 07:55:50 for Rythmos - Transient Integration for Differential Equations by   doxygen 1.8.5