diagonalTransientMain.cpp

00001 //@HEADER
00002 
00003 // ***********************************************************************
00004 //
00005 //                     Rythmos Package
00006 //                 Copyright (2006) Sandia Corporation
00007 //
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 //
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00026 //
00027 // ***********************************************************************
00028 //@HEADER 
00029 
00030 #include "EpetraExt_DiagonalTransientModel.hpp"
00031 #include "Rythmos_BackwardEulerStepper.hpp"
00032 #include "Rythmos_ImplicitBDFStepper.hpp"
00033 #include "Rythmos_ForwardSensitivityStepper.hpp"
00034 #include "Rythmos_TimeStepNonlinearSolver.hpp"
00035 #include "Rythmos_SimpleIntegrator.hpp"
00036 #include "Rythmos_SimpleIntegrationControlStrategy.hpp"
00037 #include "Rythmos_StepperAsModelEvaluator.hpp"
00038 #include "Thyra_DefaultRealLinearSolverBuilder.hpp"
00039 #include "Thyra_EpetraModelEvaluator.hpp"
00040 #include "Thyra_DirectionalFiniteDiffCalculator.hpp"
00041 #include "Thyra_TestingTools.hpp"
00042 #include "Teuchos_StandardCatchMacros.hpp"
00043 #include "Teuchos_GlobalMPISession.hpp"
00044 #include "Teuchos_DefaultComm.hpp"
00045 #include "Teuchos_VerboseObject.hpp"
00046 #include "Teuchos_XMLParameterListHelpers.hpp"
00047 #include "Teuchos_CommandLineProcessor.hpp"
00048 #include "Teuchos_VerbosityLevelCommandLineProcessorHelpers.hpp"
00049 #include "Teuchos_as.hpp"
00050 
00051 #ifdef HAVE_MPI
00052 #  include "Epetra_MpiComm.h"
00053 #else
00054 #  include "Epetra_SerialComm.h"
00055 #endif // HAVE_MPI
00056 
00057 namespace {
00058 
00059 
00060 const std::string TimeStepNonlinearSolver_name = "TimeStepNonlinearSolver";
00061 
00062 const std::string Stratimikos_name = "Stratimikos";
00063 
00064 const std::string DiagonalTransientModel_name = "DiagonalTransientModel";
00065 
00066 const std::string RythmosStepper_name = "Rythmos Stepper";
00067 
00068 const std::string RythmosIntegrator_name = "Rythmos Integrator";
00069 
00070 const std::string FdCalc_name = "FD Calc";
00071 
00072 
00073 Teuchos::RCP<const Teuchos::ParameterList>
00074 getValidParameters()
00075 {
00076   using Teuchos::RCP; using Teuchos::ParameterList;
00077   static RCP<const ParameterList> validPL;
00078   if (is_null(validPL)) {
00079     RCP<ParameterList> pl = Teuchos::parameterList();
00080     pl->sublist(TimeStepNonlinearSolver_name);
00081     pl->sublist(Stratimikos_name);
00082     pl->sublist(DiagonalTransientModel_name);
00083     pl->sublist(RythmosStepper_name);
00084     pl->sublist(RythmosIntegrator_name);
00085     pl->sublist(FdCalc_name);
00086     validPL = pl;
00087   }
00088   return validPL;
00089 }
00090 
00091 
00092 } // namespace
00093 
00094 int main(int argc, char *argv[])
00095 {
00096 
00097   using std::endl;
00098   typedef double Scalar;
00099   typedef double ScalarMag;
00100   using Teuchos::describe;
00101   using Teuchos::RCP;
00102   using Teuchos::rcp;
00103   using Teuchos::rcp_implicit_cast;
00104   using Teuchos::rcp_dynamic_cast;
00105   using Teuchos::as;
00106   using Teuchos::ParameterList;
00107   using Teuchos::CommandLineProcessor;
00108   typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
00109   typedef Thyra::ModelEvaluatorBase MEB;
00110   typedef Thyra::DefaultMultiVectorProductVectorSpace<Scalar> DMVPVS;
00111   using Thyra::productVectorBase;
00112   
00113   bool result, success = true;
00114 
00115   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00116 
00117   RCP<Epetra_Comm> epetra_comm;
00118 #ifdef HAVE_MPI
00119   epetra_comm = rcp( new Epetra_MpiComm(MPI_COMM_WORLD) );
00120 #else
00121   epetra_comm = rcp( new Epetra_SerialComm );
00122 #endif // HAVE_MPI
00123 
00124   RCP<Teuchos::FancyOStream>
00125     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00126 
00127   try {
00128 
00129     //
00130     // Read commandline options
00131     //
00132 
00133     CommandLineProcessor clp;
00134     clp.throwExceptions(false);
00135     clp.addOutputSetupOptions(true);
00136 
00137     std::string linearSolverParamsFile = "";
00138     clp.setOption( "linear-solver-params-file", &linearSolverParamsFile,
00139       "File name for XML linear solver parameters for Stratimikos" );
00140 
00141     std::string linearSolverExtraParams = "";
00142     clp.setOption( "linear-solver-extra-params", &linearSolverExtraParams,
00143       "Extra XML parameter list string for linear solver parameters for Stratimikos" );
00144 
00145     double maxStateError = 1e-6;
00146     clp.setOption( "max-state-error", &maxStateError,
00147       "The maximum allowed error in the integrated state in relation to the exact state solution" );
00148 
00149     double finalTime = 1e-3;
00150     clp.setOption( "final-time", &finalTime,
00151       "Final integration time (initial time is 0.0)" );
00152 
00153     int numTimeSteps = 10;
00154     clp.setOption( "num-time-steps", &numTimeSteps,
00155       "Number of (fixed) time steps.  If <= 0.0, then variable time steps are taken" );
00156 
00157     bool useBDF = false;
00158     clp.setOption( "use-BDF", "use-BE", &useBDF,
00159       "Use BDF or Backward Euler (BE)" );
00160 
00161     bool doFwdSensSolve = false;
00162     clp.setOption( "fwd-sens-solve", "state-solve", &doFwdSensSolve,
00163       "Do the forward sensitivity solve or just the state solve" );
00164 
00165     double maxRestateError = 0.0;
00166     clp.setOption( "max-restate-error", &maxRestateError,
00167       "The maximum allowed error between the state integrated by itself verses integrated along with DxDp" );
00168 
00169     double maxSensError = 1e-4;
00170     clp.setOption( "max-sens-error", &maxSensError,
00171       "The maximum allowed error in the integrated sensitivity in relation to"
00172       " the finite-difference sensitivity" );
00173 
00174     Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_DEFAULT;
00175     setVerbosityLevelOption( "verb-level", &verbLevel,
00176       "Top-level verbosity level.  By default, this gets deincremented as you go deeper into numerical objects.",
00177       &clp );
00178 
00179     bool testExactSensitivity = false;
00180     clp.setOption( "test-exact-sens", "no-test-exact-sens", &testExactSensitivity,
00181       "Test the exact sensitivity with finite differences or not." );
00182 
00183     bool dumpFinalSolutions = false;
00184     clp.setOption(
00185       "dump-final-solutions", "no-dump-final-solutions", &dumpFinalSolutions,
00186       "Determine if the final solutions are dumpped or not." );
00187     
00188     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00189     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00190     
00191     if ( Teuchos::VERB_DEFAULT == verbLevel )
00192       verbLevel = Teuchos::VERB_LOW;
00193 
00194     const Teuchos::EVerbosityLevel
00195       solnVerbLevel = ( dumpFinalSolutions ? Teuchos::VERB_EXTREME : verbLevel );
00196 
00197     //
00198     // Get the base parameter list that all other parameter lists will be read
00199     // from.
00200     //
00201     
00202     RCP<ParameterList>
00203       paramList = Teuchos::parameterList();
00204     if (linearSolverParamsFile.length())
00205       updateParametersFromXmlFile( linearSolverParamsFile, &*paramList );
00206     if (linearSolverExtraParams.length())
00207       updateParametersFromXmlString( linearSolverExtraParams, &*paramList );
00208 
00209     if (testExactSensitivity) {
00210       paramList->sublist(DiagonalTransientModel_name).set("Exact Solution as Response",true);
00211     }
00212 
00213     paramList->validateParameters(*getValidParameters(),0); // Only validate top level lists!
00214 
00215     //
00216     // Create the Stratimikos linear solver factory.
00217     //
00218     // This is the linear solve strategy that will be used to solve for the
00219     // linear system with the W.
00220     //
00221 
00222     Thyra::DefaultRealLinearSolverBuilder linearSolverBuilder;
00223     linearSolverBuilder.setParameterList(sublist(paramList,Stratimikos_name));
00224     RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
00225       W_factory = createLinearSolveStrategy(linearSolverBuilder);
00226     
00227     //
00228     // Create the underlying EpetraExt::ModelEvaluator
00229     //
00230 
00231     RCP<EpetraExt::DiagonalTransientModel>
00232       epetraStateModel = EpetraExt::diagonalTransientModel(
00233         epetra_comm,
00234         sublist(paramList,DiagonalTransientModel_name)
00235         );
00236 
00237     *out <<"\nepetraStateModel valid options:\n";
00238     epetraStateModel->getValidParameters()->print(
00239       *out, PLPrintOptions().indent(2).showTypes(true).showDoc(true)
00240       );
00241 
00242     //
00243     // Create the Thyra-wrapped ModelEvaluator
00244     //
00245     
00246     RCP<Thyra::ModelEvaluator<double> >
00247       stateModel = epetraModelEvaluator(epetraStateModel,W_factory);
00248 
00249     *out << "\nParameter names = " << *stateModel->get_p_names(0) << "\n";
00250     
00251     //
00252     // Create the Rythmos stateStepper
00253     //
00254 
00255     RCP<Rythmos::TimeStepNonlinearSolver<double> >
00256       nonlinearSolver = Teuchos::rcp(new Rythmos::TimeStepNonlinearSolver<double>());
00257     RCP<ParameterList>
00258       nonlinearSolverPL = sublist(paramList,TimeStepNonlinearSolver_name);
00259     nonlinearSolverPL->get("Default Tol",1e-3*maxStateError); // Set default if not set
00260     nonlinearSolver->setParameterList(nonlinearSolverPL);
00261 
00262     RCP<Rythmos::StepperBase<Scalar> > stateStepper;
00263 
00264     if (useBDF) {
00265       stateStepper = rcp(
00266         new Rythmos::ImplicitBDFStepper<double>(
00267           stateModel, nonlinearSolver
00268           )
00269         );
00270     }
00271     else {
00272       stateStepper = rcp(
00273         new Rythmos::BackwardEulerStepper<double>(
00274           stateModel, nonlinearSolver
00275           )
00276         );
00277     }
00278 
00279     *out <<"\nstateStepper:\n" << describe(*stateStepper,verbLevel);
00280     *out <<"\nstateStepper valid options:\n";
00281     stateStepper->getValidParameters()->print(
00282       *out, PLPrintOptions().indent(2).showTypes(true).showDoc(true)
00283       );
00284 
00285     stateStepper->setParameterList(sublist(paramList,RythmosStepper_name));
00286 
00287     //
00288     // Setup finite difference objects that will be used for tests
00289     //
00290 
00291     Thyra::DirectionalFiniteDiffCalculator<Scalar> fdCalc;
00292     fdCalc.setParameterList(sublist(paramList,FdCalc_name));
00293     fdCalc.setOStream(out);
00294     fdCalc.setVerbLevel(verbLevel);
00295 
00296     //
00297     // Use a StepperAsModelEvaluator to integrate the state
00298     //
00299 
00300     const MEB::InArgs<Scalar>
00301       state_ic = stateModel->getNominalValues();
00302     *out << "\nstate_ic:\n" << describe(state_ic,verbLevel);
00303 
00304     RCP<Rythmos::IntegratorBase<Scalar> > integrator;
00305     {
00306       RCP<ParameterList>
00307         integratorPL = sublist(paramList,RythmosIntegrator_name);
00308       integratorPL->set( "Take Variable Steps", as<bool>(numTimeSteps < 0) );
00309       integratorPL->set( "Fixed dt", as<double>((finalTime - state_ic.get_t())/numTimeSteps) );
00310       RCP<Rythmos::IntegratorBase<Scalar> >
00311         simpleIntegrator = Rythmos::controlledSimpleIntegrator<Scalar>(
00312           Rythmos::simpleIntegrationControlStrategy<Scalar>(integratorPL)
00313           );
00314       integrator = simpleIntegrator;
00315     }
00316 
00317     RCP<Rythmos::StepperAsModelEvaluator<Scalar> >
00318       stateIntegratorAsModel = Rythmos::stepperAsModelEvaluator(
00319         stateStepper, integrator, state_ic
00320         );
00321     stateIntegratorAsModel->setVerbLevel(verbLevel);
00322     
00323     *out << "\nUse the StepperAsModelEvaluator to integrate state x(p,finalTime) ... \n";
00324     
00325     RCP<Thyra::VectorBase<Scalar> > x_final;
00326 
00327     {
00328       
00329       Teuchos::OSTab tab(out);
00330 
00331       x_final = createMember(stateIntegratorAsModel->get_g_space(0));
00332       
00333       eval_g(
00334         *stateIntegratorAsModel,
00335         0, *state_ic.get_p(0),
00336         finalTime,
00337         0, &*x_final
00338         );
00339 
00340       *out
00341         << "\nx_final = x(p,finalTime) evaluated using stateIntegratorAsModel:\n"
00342         << describe(*x_final,solnVerbLevel);
00343 
00344     }
00345 
00346     //
00347     // Test the integrated state against the exact analytical state solution
00348     //
00349 
00350     RCP<const Thyra::VectorBase<Scalar> >
00351       exact_x_final = create_Vector(
00352         epetraStateModel->getExactSolution(finalTime),
00353         stateModel->get_x_space()
00354         );
00355     
00356     result = Thyra::testRelNormDiffErr(
00357       "exact_x_final", *exact_x_final, "x_final", *x_final,
00358       "maxStateError", maxStateError, "warningTol", 1.0, // Don't warn
00359       &*out, solnVerbLevel
00360       );
00361     if (!result) success = false;
00362 
00363     //
00364     // Solve and test the forward sensitivity computation
00365     //
00366       
00367     if (doFwdSensSolve) {
00368 
00369       //
00370       // Create the forward sensitivity stepper
00371       //
00372       
00373       RCP<Rythmos::ForwardSensitivityStepper<Scalar> >
00374         stateAndSensStepper = Rythmos::forwardSensitivityStepper<Scalar>(
00375           stateModel, 0, stateModel->getNominalValues(),
00376           stateStepper, nonlinearSolver
00377           );
00378       // The above call will result in stateStepper and nonlinearSolver being
00379       // cloned.  This helps to ensure consistency between the state and
00380       // sensitivity computations!
00381 
00382       //
00383       // Set the initial condition for the state and forward sensitivities
00384       //
00385 
00386       RCP<Thyra::VectorBase<Scalar> > s_bar_init
00387         = createMember(stateAndSensStepper->getFwdSensModel()->get_x_space());
00388       assign( &*s_bar_init, 0.0 );
00389       RCP<Thyra::VectorBase<Scalar> > s_bar_dot_init
00390         = createMember(stateAndSensStepper->getFwdSensModel()->get_x_space());
00391       assign( &*s_bar_dot_init, 0.0 );
00392       // Above, I believe that these are the correct initial conditions for
00393       // s_bar and s_bar_dot given how the EpetraExt::DiagonalTransientModel
00394       // is currently implemented!
00395 
00396       RCP<const Rythmos::StateAndForwardSensitivityModelEvaluator<Scalar> >
00397         stateAndSensModel = stateAndSensStepper->getStateAndFwdSensModel();
00398 
00399       MEB::InArgs<Scalar>
00400         state_and_sens_ic = stateAndSensStepper->getModel()->createInArgs();
00401 
00402       // Copy time, parameters etc.
00403       state_and_sens_ic.setArgs(state_ic);
00404       // Set initial condition for x_bar = [ x; s_bar ]
00405       state_and_sens_ic.set_x(
00406         stateAndSensModel->create_x_bar_vec(state_ic.get_x(),s_bar_init)
00407         );
00408       // Set initial condition for x_bar_dot = [ x_dot; s_bar_dot ]
00409       state_and_sens_ic.set_x_dot(
00410         stateAndSensModel->create_x_bar_vec(state_ic.get_x_dot(),s_bar_dot_init)
00411         );
00412 
00413       *out << "\nstate_and_sens_ic:\n" << describe(state_and_sens_ic,verbLevel);
00414  
00415       stateAndSensStepper->setInitialCondition(state_and_sens_ic);
00416 
00417       //
00418       // Use a StepperAsModelEvaluator to integrate the state+sens
00419       //
00420     
00421       RCP<Rythmos::StepperAsModelEvaluator<Scalar> >
00422         stateAndSensIntegratorAsModel = Rythmos::stepperAsModelEvaluator(
00423           rcp_implicit_cast<Rythmos::StepperBase<Scalar> >(stateAndSensStepper),
00424           integrator, state_and_sens_ic
00425           );
00426       stateAndSensIntegratorAsModel->setVerbLevel(verbLevel);
00427     
00428       *out << "\nUse the StepperAsModelEvaluator to integrate state + sens x_bar(p,finalTime) ... \n";
00429     
00430       RCP<Thyra::VectorBase<Scalar> > x_bar_final;
00431 
00432       {
00433       
00434         Teuchos::OSTab tab(out);
00435 
00436         x_bar_final = createMember(stateAndSensIntegratorAsModel->get_g_space(0));
00437       
00438         eval_g(
00439           *stateAndSensIntegratorAsModel,
00440           0, *state_ic.get_p(0),
00441           finalTime,
00442           0, &*x_bar_final
00443           );
00444 
00445         *out
00446           << "\nx_bar_final = x_bar(p,finalTime) evaluated using stateAndSensIntegratorAsModel:\n"
00447           << describe(*x_bar_final,solnVerbLevel);
00448 
00449       }
00450 
00451       //
00452       // Test that the state computed above is same as computed initially!
00453       //
00454 
00455       *out << "\nChecking that x(p,finalTime) computed as part of x_bar above is the same ...\n";
00456 
00457       {
00458 
00459         Teuchos::OSTab tab(out);
00460 
00461         RCP<const Thyra::VectorBase<Scalar> >
00462           x_in_x_bar_final = productVectorBase<Scalar>(x_bar_final)->getVectorBlock(0);
00463 
00464         result = Thyra::testRelNormDiffErr<Scalar>(
00465           "x_final", *x_final,
00466           "x_in_x_bar_final", *x_in_x_bar_final,
00467           "maxRestateError", maxRestateError,
00468           "warningTol", 1.0, // Don't warn
00469           &*out, solnVerbLevel
00470           );
00471         if (!result) success = false;
00472 
00473       }
00474 
00475       //
00476       // Compute DxDp using finite differences
00477       //
00478 
00479       *out << "\nApproximating DxDp(p,t) using directional finite differences of integrator for x(p,t) ...\n";
00480     
00481       RCP<Thyra::MultiVectorBase<Scalar> > DxDp_fd_final;
00482 
00483       {
00484       
00485         Teuchos::OSTab tab(out);
00486       
00487       
00488         MEB::InArgs<Scalar>
00489           fdBasePoint = stateIntegratorAsModel->createInArgs();
00490       
00491         fdBasePoint.set_t(finalTime);
00492         fdBasePoint.set_p(0,stateModel->getNominalValues().get_p(0));
00493       
00494         DxDp_fd_final = createMembers(
00495           stateIntegratorAsModel->get_g_space(0),
00496           stateIntegratorAsModel->get_p_space(0)->dim()
00497           );
00498       
00499         typedef Thyra::DirectionalFiniteDiffCalculatorTypes::SelectedDerivatives
00500           SelectedDerivatives;
00501       
00502         MEB::OutArgs<Scalar> fdOutArgs =
00503           fdCalc.createOutArgs(
00504             *stateIntegratorAsModel,
00505             SelectedDerivatives().supports(MEB::OUT_ARG_DgDp,0,0)
00506             );
00507         fdOutArgs.set_DgDp(0,0,DxDp_fd_final);
00508       
00509         // Silence the model evaluators that are called.  The fdCal object
00510         // will show all of the inputs and outputs for each call.
00511         stateStepper->setVerbLevel(Teuchos::VERB_NONE);
00512         stateIntegratorAsModel->setVerbLevel(Teuchos::VERB_NONE);
00513       
00514         fdCalc.calcDerivatives(
00515           *stateIntegratorAsModel, fdBasePoint,
00516           stateIntegratorAsModel->createOutArgs(), // Don't bother with function value
00517           fdOutArgs
00518           );
00519         
00520         *out
00521           << "\nFinite difference DxDp_fd_final = DxDp(p,finalTime): "
00522           << describe(*DxDp_fd_final,solnVerbLevel);
00523 
00524       }
00525 
00526       //
00527       // Test that the integrated sens and the F.D. sens are similar
00528       //
00529 
00530       *out << "\nChecking that integrated DxDp(p,finalTime) and finite-diff DxDp(p,finalTime) are similar ...\n";
00531 
00532       {
00533 
00534         Teuchos::OSTab tab(out);
00535 
00536         RCP<const Thyra::VectorBase<Scalar> >
00537           DxDp_vec_final = Thyra::productVectorBase<Scalar>(x_bar_final)->getVectorBlock(1);
00538       
00539         RCP<const Thyra::VectorBase<Scalar> >
00540           DxDp_fd_vec_final = Thyra::multiVectorProductVector(
00541             rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> >(
00542               DxDp_vec_final->range()
00543               ),
00544             DxDp_fd_final
00545             );
00546 
00547         result = Thyra::testRelNormDiffErr(
00548           "DxDp_vec_final", *DxDp_vec_final,
00549           "DxDp_fd_vec_final", *DxDp_fd_vec_final,
00550           "maxSensError", maxSensError,
00551           "warningTol", 1.0, // Don't warn
00552           &*out, solnVerbLevel
00553           );
00554         if (!result) success = false;
00555 
00556       }
00557 
00558     }
00559     
00560   }
00561   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success);
00562 
00563   if(success)
00564     *out << "\nEnd Result: TEST PASSED" << endl;
00565   else
00566     *out << "\nEnd Result: TEST FAILED" << endl;
00567   
00568   return ( success ? 0 : 1 );
00569 
00570 } // end main() [Doxygen looks for this!]
00571 

Generated on Tue Oct 20 12:46:07 2009 for Rythmos - Transient Integration for Differential Equations by doxygen 1.4.7