Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_ForwardSensitivityStepperTester_def.hpp
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //                           Rythmos Package
00005 //                 Copyright (2006) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #ifndef Rythmos_FORWARD_SENSITIVITY_STEPPER_TESTER_DEF_H
00030 #define Rythmos_FORWARD_SENSITIVITY_STEPPER_TESTER_DEF_H
00031 
00032 
00033 #include "Rythmos_ForwardSensitivityStepperTester_decl.hpp"
00034 #include "Rythmos_ForwardSensitivityStepper.hpp"
00035 #include "Rythmos_StepperAsModelEvaluator.hpp"
00036 #include "Thyra_DirectionalFiniteDiffCalculator.hpp"
00037 #include "Thyra_TestingTools.hpp"
00038 
00039 
00040 template<class Scalar>
00041 Teuchos::RCP<Rythmos::ForwardSensitivityStepperTester<Scalar> >
00042 Rythmos::forwardSensitivityStepperTester()
00043 {
00044   return Teuchos::rcp(new ForwardSensitivityStepperTester<Scalar>);
00045 }
00046 
00047 
00048 template<class Scalar>
00049 Teuchos::RCP<Rythmos::ForwardSensitivityStepperTester<Scalar> >
00050 Rythmos::forwardSensitivityStepperTester(const RCP<ParameterList> &paramList)
00051 {
00052   const RCP<Rythmos::ForwardSensitivityStepperTester<Scalar> > fsst =
00053     forwardSensitivityStepperTester<Scalar>();
00054   fsst->setParameterList(paramList);
00055   return fsst;
00056 }
00057 
00058 
00059 namespace Rythmos {
00060 
00061 
00062 // Overridden from ParameterListAcceptor
00063 
00064 
00065 template<class Scalar>
00066 void ForwardSensitivityStepperTester<Scalar>::setParameterList(
00067   RCP<ParameterList> const& paramList
00068   )
00069 {
00070   namespace FSSTU = ForwardSensitivityStepperTesterUtils;
00071   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00072   this->setMyParamList(paramList);
00073   errorTol_ = Teuchos::getParameter<double>(*paramList, FSSTU::ErrorTol_name);
00074 }
00075 
00076 
00077 template<class Scalar>
00078 RCP<const ParameterList>
00079 ForwardSensitivityStepperTester<Scalar>::getValidParameters() const
00080 {
00081   namespace FSSTU = ForwardSensitivityStepperTesterUtils;
00082   static RCP<const ParameterList> validPL;
00083   if (is_null(validPL)) {
00084     const RCP<ParameterList> pl = Teuchos::parameterList();
00085     pl->set(FSSTU::ErrorTol_name, FSSTU::ErrorTol_default);
00086     pl->sublist(FSSTU::FdCalc_name).disableRecursiveValidation().setParameters(
00087       *Thyra::DirectionalFiniteDiffCalculator<Scalar>().getValidParameters()
00088       );
00089     validPL = pl;
00090   }
00091   return validPL;
00092 }
00093 
00094 
00095 template<class Scalar>
00096 bool ForwardSensitivityStepperTester<Scalar>::testForwardSens(
00097   const Ptr<IntegratorBase<Scalar> > &fwdSensIntegrator
00098   )
00099 {
00100 
00101   using Teuchos::outArg;
00102   using Teuchos::rcp_dynamic_cast;
00103   using Teuchos::OSTab;
00104   using Teuchos::sublist;
00105   typedef Thyra::ModelEvaluatorBase MEB;
00106   namespace FSSTU = ForwardSensitivityStepperTesterUtils;
00107 
00108   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00109   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00110 
00111   const RCP<ParameterList> paramList = this->getMyNonconstParamList();
00112 
00113   bool success = true;
00114 
00115   // A) Extract out and dynamic cast the forward sens stepper
00116   RCP<ForwardSensitivityStepper<Scalar> > fwdSensStepper =
00117     rcp_dynamic_cast<ForwardSensitivityStepper<Scalar> >(
00118       fwdSensIntegrator->getNonconstStepper()
00119       );
00120 
00121   // B) Extract out the fwd state stepper
00122   RCP<StepperBase<Scalar> > stateStepper = fwdSensStepper->getNonconstStateStepper();
00123 
00124   // C) Get the initial condition for the state
00125 
00126   MEB::InArgs<Scalar> state_and_sens_ic = fwdSensStepper->getInitialCondition();
00127   MEB::InArgs<Scalar> state_ic = 
00128     extractStateInitialCondition(*fwdSensStepper, state_and_sens_ic);
00129 
00130   // D) Create a StepperAsModelEvalutor for the fwd state calculation
00131 
00132   RCP<StepperAsModelEvaluator<Scalar> >
00133     stateIntegratorAsModel = stepperAsModelEvaluator(
00134       stateStepper, fwdSensIntegrator->cloneIntegrator(), state_ic
00135       );
00136   //stateIntegratorAsModel->setVerbLevel(verbLevel);
00137 
00138   // E) Compute discrete forward sensitivities using fwdSensIntegrator
00139 
00140   const Scalar finalTime = fwdSensIntegrator->getFwdTimeRange().upper();
00141 
00142   *out << "\nCompute discrete forward sensitivities using fwdSensIntegrator ...\n";
00143 
00144   RCP<const VectorBase<Scalar> > x_bar_final, x_bar_dot_final;
00145   {
00146     OSTab tab(out);
00147     get_fwd_x_and_x_dot( *fwdSensIntegrator, finalTime,
00148       outArg(x_bar_final), outArg(x_bar_dot_final) );
00149     *out
00150       << "\nx_bar_final = x_bar(p,finalTime) evaluated using stateAndSensIntegratorAsModel:\n"
00151       << describe(*x_bar_final, verbLevel);
00152   }
00153   
00154   // F) Compute discrete forward sensitivities by finite differences
00155 
00156   *out << "\nCompute discrete forward sensitivities by finite differences ...\n";
00157 
00158   RCP<Thyra::MultiVectorBase<Scalar> > DxDp_fd_final;
00159   {
00160     OSTab tab(out);
00161 
00162     Thyra::DirectionalFiniteDiffCalculator<Scalar> fdCalc;
00163     if (nonnull(paramList))
00164       fdCalc.setParameterList(sublist(paramList, FSSTU::FdCalc_name));
00165     fdCalc.setOStream(out);
00166     fdCalc.setVerbLevel(Teuchos::VERB_NONE);
00167     
00168     MEB::InArgs<Scalar>
00169       fdBasePoint = stateIntegratorAsModel->createInArgs();
00170     
00171     fdBasePoint.setArgs(state_ic, true); // Set the parameters
00172     fdBasePoint.set_t(finalTime);
00173     
00174     const int g_index = 0;
00175     const int p_index = fwdSensStepper->getFwdSensModel()->get_p_index();
00176     
00177     DxDp_fd_final =
00178       createMembers(
00179         stateIntegratorAsModel->get_g_space(g_index),
00180         stateIntegratorAsModel->get_p_space(p_index)->dim()
00181         );
00182     
00183     typedef Thyra::DirectionalFiniteDiffCalculatorTypes::SelectedDerivatives
00184       SelectedDerivatives;
00185     
00186     MEB::OutArgs<Scalar> fdOutArgs =
00187       fdCalc.createOutArgs(
00188         *stateIntegratorAsModel,
00189         SelectedDerivatives().supports(MEB::OUT_ARG_DgDp, g_index, p_index)
00190         );
00191     fdOutArgs.set_DgDp(g_index, p_index, DxDp_fd_final);
00192     
00193     // Silence the model evaluators that are called.  The fdCal object
00194     // will show all of the inputs and outputs for each call.
00195     stateStepper->setVerbLevel(Teuchos::VERB_NONE);
00196     stateIntegratorAsModel->setVerbLevel(Teuchos::VERB_NONE);
00197     
00198     fdCalc.calcDerivatives(
00199       *stateIntegratorAsModel, fdBasePoint,
00200       stateIntegratorAsModel->createOutArgs(), // Don't bother with function value
00201       fdOutArgs
00202       );
00203     
00204     *out
00205       << "\nFinite difference DxDp_fd_final = DxDp(p,finalTime): "
00206       << describe(*DxDp_fd_final, verbLevel);
00207 
00208   }
00209 
00210   // G) Compare the integrated and finite difference sensitivities
00211 
00212   *out << "\nChecking that integrated DxDp(p,finalTime) and finite-diff DxDp(p,finalTime) are similar ...\n";
00213   
00214   {
00215     
00216     Teuchos::OSTab tab(out);
00217     
00218     RCP<const Thyra::VectorBase<Scalar> >
00219       DxDp_vec_final = Thyra::productVectorBase<Scalar>(x_bar_final)->getVectorBlock(1);
00220     
00221     RCP<const Thyra::VectorBase<Scalar> >
00222       DxDp_fd_vec_final = Thyra::multiVectorProductVector(
00223         rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> >(
00224           DxDp_vec_final->range()
00225           ),
00226         DxDp_fd_final
00227         );
00228 
00229     const bool result = Thyra::testRelNormDiffErr(
00230       "DxDp_vec_final", *DxDp_vec_final,
00231       "DxDp_fd_vec_final", *DxDp_fd_vec_final,
00232       "maxSensError", errorTol_,
00233       "warningTol", 1.0, // Don't warn
00234       &*out, verbLevel
00235       );
00236     if (!result)
00237       success = false;
00238     
00239   }
00240 
00241   return success;
00242 
00243 }
00244 
00245 
00246 template<class Scalar>
00247 ForwardSensitivityStepperTester<Scalar>::ForwardSensitivityStepperTester()
00248   : errorTol_(ForwardSensitivityStepperTesterUtils::ErrorTol_default)
00249 {}
00250 
00251 
00252 } // namespace Rythmos
00253 
00254 
00255 // 
00256 // Explicit Instantiation macro
00257 //
00258 // Must be expanded from within the Rythmos namespace!
00259 //
00260 
00261 #define FORWARD_SENSITIVITY_STEPPER_INSTANT(SCALAR) \
00262   \
00263   template class ForwardSensitivityStepperTester<SCALAR >; \
00264   \
00265   template RCP<ForwardSensitivityStepperTester<SCALAR > > \
00266   forwardSensitivityStepperTester(); \
00267   \
00268   template RCP<ForwardSensitivityStepperTester<SCALAR> > \
00269   forwardSensitivityStepperTester(const RCP<ParameterList> &paramList);
00270 
00271 
00272 #endif //Rythmos_FORWARD_SENSITIVITY_STEPPER_TESTER_DEF_H
 All Classes Functions Variables Typedefs Friends