Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_BasicDiscreteAdjointStepperTester_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., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #ifndef Rythmos_BASIC_DISCRETE_ADJOINT_STEPPER_TESTER_DEF_H
00030 #define Rythmos_BASIC_DISCRETE_ADJOINT_STEPPER_TESTER_DEF_H
00031 
00032 
00033 #include "Rythmos_BasicDiscreteAdjointStepperTester_decl.hpp"
00034 #include "Rythmos_ForwardSensitivityStepper.hpp"
00035 #include "Rythmos_AdjointModelEvaluator.hpp"
00036 #include "Rythmos_DefaultIntegrator.hpp" // ToDo: Remove when we can!
00037 #include "Thyra_LinearNonlinearSolver.hpp"
00038 #include "Thyra_VectorBase.hpp"
00039 #include "Thyra_VectorStdOps.hpp"
00040 #include "Thyra_TestingTools.hpp"
00041 
00042 
00043 template<class Scalar>
00044 Teuchos::RCP<Rythmos::BasicDiscreteAdjointStepperTester<Scalar> >
00045 Rythmos::basicDiscreteAdjointStepperTester()
00046 {
00047   return Teuchos::rcp(new BasicDiscreteAdjointStepperTester<Scalar>);
00048 }
00049 
00050 
00051 template<class Scalar>
00052 Teuchos::RCP<Rythmos::BasicDiscreteAdjointStepperTester<Scalar> >
00053 Rythmos::basicDiscreteAdjointStepperTester(const RCP<ParameterList> &paramList)
00054 {
00055   const RCP<Rythmos::BasicDiscreteAdjointStepperTester<Scalar> > bdast =
00056     basicDiscreteAdjointStepperTester<Scalar>();
00057   bdast->setParameterList(paramList);
00058   return bdast;
00059 }
00060 
00061 
00062 namespace Rythmos {
00063 
00064 
00065 // Overridden from ParameterListAcceptor
00066 
00067 
00068 template<class Scalar>
00069 void BasicDiscreteAdjointStepperTester<Scalar>::setParameterList(
00070   RCP<ParameterList> const& paramList
00071   )
00072 {
00073   namespace BDASTU = BasicDiscreteAdjointStepperTesterUtils;
00074   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00075   this->setMyParamList(paramList);
00076   errorTol_ = Teuchos::getParameter<double>(*paramList, BDASTU::ErrorTol_name);
00077 }
00078 
00079 
00080 template<class Scalar>
00081 RCP<const ParameterList>
00082 BasicDiscreteAdjointStepperTester<Scalar>::getValidParameters() const
00083 {
00084   namespace BDASTU = BasicDiscreteAdjointStepperTesterUtils;
00085   static RCP<const ParameterList> validPL;
00086   if (is_null(validPL)) {
00087     const RCP<ParameterList> pl = Teuchos::parameterList();
00088     pl->set(BDASTU::ErrorTol_name, BDASTU::ErrorTol_default);
00089     validPL = pl;
00090   }
00091   return validPL;
00092 }
00093 
00094 
00095 template<class Scalar>
00096 bool BasicDiscreteAdjointStepperTester<Scalar>::testAdjointStepper(
00097   Thyra::ModelEvaluator<Scalar>& adjointModel,
00098   const Ptr<IntegratorBase<Scalar> >& forwardIntegrator
00099   )
00100 {
00101 
00102   using Teuchos::describe;
00103   using Teuchos::outArg;
00104   using Teuchos::rcp_dynamic_cast;
00105   using Teuchos::dyn_cast;
00106   using Teuchos::dyn_cast;
00107   using Teuchos::OSTab;
00108   using Teuchos::sublist;
00109   typedef Thyra::ModelEvaluatorBase MEB;
00110   using Thyra::createMember;
00111   namespace BDASTU = BasicDiscreteAdjointStepperTesterUtils;
00112 
00113   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00114   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00115 
00116   const RCP<ParameterList> paramList = this->getMyNonconstParamList();
00117 
00118   bool success = true;
00119 
00120   OSTab tab(out);
00121 
00122   //
00123   *out << "\n*** Entering BasicDiscreteAdjointStepperTester<Scalar>::testAdjointStepper(...) ...\n";
00124   //
00125 
00126   const TimeRange<Scalar> fwdTimeRange = forwardIntegrator->getFwdTimeRange();
00127   const RCP<Rythmos::StepperBase<Scalar> > fwdStepper =
00128     forwardIntegrator->getNonconstStepper();
00129   const RCP<const Thyra::ModelEvaluator<Scalar> > fwdModel =
00130     fwdStepper->getModel();
00131   const RCP<const Thyra::VectorSpaceBase<Scalar> > f_space = fwdModel->get_f_space();
00132 
00133   //
00134   *out << "\nA) Construct the IC basis B ...\n";
00135   //
00136 
00137   const RCP<Thyra::MultiVectorBase<Scalar> > B =
00138     createMembers(fwdModel->get_x_space(), 1, "B");
00139   Thyra::seed_randomize<Scalar>(0);
00140   Thyra::randomize<Scalar>( -1.0, +1.0, B.ptr() );
00141   
00142   *out << "\nB: " << describe(*B, verbLevel);
00143 
00144   //
00145   *out << "\nB) Construct the forward sensitivity integrator ...\n";
00146   //
00147 
00148   const RCP<ForwardSensitivityStepper<Scalar> > fwdSensStepper =
00149     forwardSensitivityStepper<Scalar>();
00150   fwdSensStepper->initializeSyncedSteppersInitCondOnly(
00151     fwdModel,
00152     B->domain(), // p_space
00153     fwdStepper->getInitialCondition(),
00154     fwdStepper,
00155     dyn_cast<SolverAcceptingStepperBase<Scalar> >(*fwdStepper).getNonconstSolver()
00156     );
00157   *out << "\nfwdSensStepper: " << describe(*fwdSensStepper, verbLevel);
00158 
00159   const MEB::InArgs<Scalar> fwdIC = 
00160     forwardIntegrator->getStepper()->getInitialCondition();
00161 
00162   MEB::InArgs<Scalar> fwdSensIC = 
00163     createStateAndSensInitialCondition<Scalar>(*fwdSensStepper, fwdIC, B);
00164   fwdSensStepper->setInitialCondition(fwdSensIC);
00165 
00166   const RCP<IntegratorBase<Scalar> > fwdSensIntegrator =
00167     forwardIntegrator->cloneIntegrator();
00168   fwdSensIntegrator->setStepper(fwdSensStepper,
00169     forwardIntegrator->getFwdTimeRange().upper());
00170   *out << "\nfwdSensIntegrator: " << describe(*fwdSensIntegrator, verbLevel);
00171 
00172   //
00173   *out << "\nC) Construct the adjoint sensitivity integrator ...\n";
00174   //
00175 
00176   RCP<AdjointModelEvaluator<Scalar> > adjModel =
00177     adjointModelEvaluator<Scalar>(fwdModel, fwdTimeRange);
00178   adjModel->setFwdStateSolutionBuffer(
00179     dyn_cast<DefaultIntegrator<Scalar> >(*forwardIntegrator).getTrailingInterpolationBuffer()
00180     );
00186   *out << "\nadjModel: " << describe(*adjModel, verbLevel);
00187   
00188   // lambda(t_final) = 0.0 (for now)
00189   const RCP<Thyra::VectorBase<Scalar> > lambda_ic = createMember(f_space);
00190   V_S( lambda_ic.ptr(), 0.0 );
00191   
00192   // lambda_dot(t_final,i) = 0.0
00193   const RCP<Thyra::VectorBase<Scalar> > lambda_dot_ic = createMember(f_space);
00194   Thyra::V_S( lambda_dot_ic.ptr(), 0.0 );
00195   
00196   MEB::InArgs<Scalar> adj_ic = adjModel->getNominalValues();
00197   adj_ic.set_x(lambda_ic);
00198   adj_ic.set_x_dot(lambda_dot_ic);
00199   *out << "\nadj_ic: " << describe(adj_ic, Teuchos::VERB_EXTREME);
00200   
00201   const RCP<Thyra::LinearNonlinearSolver<Scalar> > adjTimeStepSolver =
00202     Thyra::linearNonlinearSolver<Scalar>();
00203   const RCP<Rythmos::StepperBase<Scalar> > adjStepper =
00204     forwardIntegrator->getStepper()->cloneStepperAlgorithm();
00205   *out << "\nadjStepper: " << describe(*adjStepper, verbLevel);
00206  
00207   adjStepper->setInitialCondition(adj_ic);
00208  
00209   const RCP<IntegratorBase<Scalar> > adjIntegrator = forwardIntegrator->cloneIntegrator();
00210   adjIntegrator->setStepper(adjStepper, forwardIntegrator->getFwdTimeRange().length());
00211 
00212   //
00213   *out << "\nD) Solve the forawrd problem storing the state along the way ...\n";
00214   //
00215 
00216   const double t_final = fwdTimeRange.upper();
00217   RCP<const Thyra::VectorBase<Scalar> > x_final, x_dot_final;
00218   get_fwd_x_and_x_dot( *forwardIntegrator, t_final, outArg(x_final), outArg(x_dot_final) );
00219 
00220   *out << "\nt_final = " << t_final << "\n";
00221   *out << "\nx_final: " << *x_final;
00222   *out << "\nx_dot_final: " << *x_dot_final;
00223 
00224   //
00225   *out << "\nE) Solve the forawrd sensitivity equations and compute fwd reduced sens ...\n";
00226   //
00227 
00228   // Set the initial condition
00229   TEUCHOS_TEST_FOR_EXCEPT(true);
00230 
00231   // Solve the fwd sens equations
00232   TEUCHOS_TEST_FOR_EXCEPT(true);
00233 
00234   // Compute the reduced gradient at t_f
00235   TEUCHOS_TEST_FOR_EXCEPT(true);
00236 
00237   //
00238   *out << "\nF) Solve the adjoint equations and compute adj reduced sens ...\n";
00239   //
00240 
00241   // Compute and set the adjoint initial condition
00242   TEUCHOS_TEST_FOR_EXCEPT(true);
00243 
00244   // Solve the adjoint equations
00245   TEUCHOS_TEST_FOR_EXCEPT(true);
00246 
00247   // Compute the reduced gradient at t_0
00248   TEUCHOS_TEST_FOR_EXCEPT(true);
00249 
00250   //
00251   *out << "\nG) Compare forward and adjoint reduced sens ...\n";
00252   //
00253 
00254   TEUCHOS_TEST_FOR_EXCEPT(true);
00255 
00256   //
00257   *out << "\n*** Leaving BasicDiscreteAdjointStepperTester<Scalar>::testAdjointStepper(...) ...\n";
00258   //
00259 
00260   return success;
00261 
00262 }
00263 
00264 
00265 // private:
00266 
00267 
00268 template<class Scalar>
00269 BasicDiscreteAdjointStepperTester<Scalar>::BasicDiscreteAdjointStepperTester()
00270   : errorTol_(BasicDiscreteAdjointStepperTesterUtils::ErrorTol_default)
00271 {}
00272 
00273 
00274 } // namespace Rythmos
00275 
00276 
00277 // 
00278 // Explicit Instantiation macro
00279 //
00280 // Must be expanded from within the Rythmos namespace!
00281 //
00282 
00283 #define FORWARD_SENSITIVITY_STEPPER_INSTANT(SCALAR) \
00284   \
00285   template class BasicDiscreteAdjointStepperTester<SCALAR >; \
00286   \
00287   template RCP<BasicDiscreteAdjointStepperTester<SCALAR > > \
00288   basicDiscreteAdjointStepperTester(); \
00289   \
00290   template RCP<BasicDiscreteAdjointStepperTester<SCALAR> > \
00291   basicDiscreteAdjointStepperTester(const RCP<ParameterList> &paramList);
00292 
00293 
00294 #endif //Rythmos_BASIC_DISCRETE_ADJOINT_STEPPER_TESTER_DEF_H
 All Classes Functions Variables Typedefs Friends