MOOCHO Version of the Day
DiagonalTransientInversionMain.cpp
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00006 //                  Copyright (2003) 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 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00039 // 
00040 // ***********************************************************************
00041 // @HEADER
00042 */
00043 
00044 
00045 #include "EpetraExt_DiagonalTransientModel.hpp"
00046 #include "MoochoPack_MoochoThyraSolver.hpp"
00047 #include "Rythmos_BackwardEulerStepper.hpp"
00048 #include "Rythmos_ImplicitBDFStepper.hpp"
00049 #include "Rythmos_ForwardSensitivityStepper.hpp"
00050 #include "Rythmos_TimeStepNonlinearSolver.hpp"
00051 #include "Rythmos_DefaultIntegrator.hpp"
00052 #include "Rythmos_SimpleIntegrationControlStrategy.hpp"
00053 #include "Rythmos_StepperAsModelEvaluator.hpp"
00054 #include "Rythmos_ForwardSensitivityIntegratorAsModelEvaluator.hpp"
00055 #include "Thyra_EpetraModelEvaluator.hpp"
00056 #include "Thyra_DefaultSpmdMultiVectorFileIO.hpp"
00057 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
00058 #include "Thyra_DefaultInverseModelEvaluator.hpp"
00059 #include "Thyra_ModelEvaluatorHelpers.hpp"
00060 #include "Thyra_DefaultFiniteDifferenceModelEvaluator.hpp"
00061 #include "Thyra_TestingTools.hpp"
00062 #include "Teuchos_GlobalMPISession.hpp"
00063 #include "Teuchos_CommandLineProcessor.hpp"
00064 #include "Teuchos_VerbosityLevelCommandLineProcessorHelpers.hpp"
00065 #include "Teuchos_StandardCatchMacros.hpp"
00066 #include "Teuchos_VerboseObject.hpp"
00067 #include "Teuchos_XMLParameterListHelpers.hpp"
00068 #ifdef HAVE_MPI
00069 #  include "Epetra_MpiComm.h"
00070 #else
00071 #  include "Epetra_SerialComm.h"
00072 #endif
00073 
00074 
00075 namespace {
00076 
00077 
00078 typedef AbstractLinAlgPack::value_type  Scalar;
00079 
00080 
00081 } // namespace
00082 
00083 
00084 int main( int argc, char* argv[] )
00085 {
00086 
00087   using std::endl;
00088   using Teuchos::rcp;
00089   using Teuchos::RCP;
00090   using Teuchos::rcp_dynamic_cast;
00091   using Teuchos::OSTab;
00092   using Teuchos::Array;
00093   using Teuchos::includesVerbLevel;
00094   using Teuchos::ParameterList;
00095   using Teuchos::CommandLineProcessor;
00096   typedef Teuchos::ParameterList::PrintOptions PLPO;
00097   typedef Thyra::ModelEvaluatorBase MEB;
00098   typedef Thyra::DefaultMultiVectorProductVectorSpace<Scalar> DMVPVS;
00099   using MoochoPack::MoochoSolver;
00100   using MoochoPack::MoochoThyraSolver;
00101 
00102   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00103 
00104   RCP<Epetra_Comm> epetra_comm;
00105 #ifdef HAVE_MPI
00106   epetra_comm = rcp( new Epetra_MpiComm(MPI_COMM_WORLD) );
00107 #else
00108   epetra_comm = rcp( new Epetra_SerialComm );
00109 #endif // HAVE_MPI
00110   
00111   bool success = true, result;
00112 
00113   Teuchos::RCP<Teuchos::FancyOStream>
00114     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00115   
00116   // Rythmos uses long names so enlarge the size of this!
00117   out->setMaxLenLinePrefix(20);
00118 
00119   try {
00120   
00121     // Create the solver objects
00122     Stratimikos::DefaultLinearSolverBuilder lowsfCreator;
00123     MoochoThyraSolver moochoThyraSolver;
00124 
00125     //
00126     // A) Get options from the command line and get the parameter list
00127     //
00128 
00129     CommandLineProcessor  clp;
00130     clp.throwExceptions(false);
00131     clp.addOutputSetupOptions(true);
00132 
00133     moochoThyraSolver.setupCLP(&clp);
00134 
00135     std::string paramsFileName = "";
00136     clp.setOption( "params-file", &paramsFileName,
00137       "File name for varaious parameter lists in xml format.",
00138       true );
00139 
00140     bool useBDF = false;
00141     clp.setOption( "use-BDF", "use-BE", &useBDF,
00142       "Use BDF or Backward Euler (BE)" );
00143 
00144     double finalTime = 1e-3;
00145     clp.setOption( "final-time", &finalTime,
00146       "Final integration time (initial time is 0.0)" );
00147     
00148     int numObservationPoints= 10;
00149     clp.setOption( "num-observations", &numObservationPoints,
00150       "Number of state observations to take as basis for inversion" );
00151 
00152     double scaleParamsBy = 1e-2;
00153     clp.setOption( "scale-params-by", &scaleParamsBy,
00154       "Amount to scale the parameters by to perturb them" );
00155 
00156     Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_DEFAULT;
00157     setVerbosityLevelOption( "verb-level", &verbLevel,
00158       "Top-level verbosity level.  By default, this gets deincremented"
00159       " as you go deeper into numerical objects.",
00160       &clp );
00161 
00162     bool printValidOptions = false;
00163     clp.setOption( "print-valid-options", "no-print-valid-options",
00164       &printValidOptions, "Print the valid options or not" );
00165 
00166     bool testTransientModelEvaluator = true;
00167     clp.setOption( "test-transient-model", "no-test-transient-model", &testTransientModelEvaluator,
00168       "Test the Rythmos::ForwardSensitivityIntegratorAsModelEvaluator object or not." );
00169 
00170     double maxSensError = 1e-4;
00171     clp.setOption( "max-sens-error", &maxSensError,
00172       "The maximum allowed error in the integrated sensitivity in relation to"
00173       " the finite-difference sensitivity" );
00174 
00175     bool doInverseProblem = true;
00176     clp.setOption( "do-inv-prob", "no-do-inv-prob", &doInverseProblem,
00177       "Run and test the inverse problem or not" );
00178 
00179     double p_inv_err_tol = 1e-8;
00180     clp.setOption( "max-p-inv-error", &p_inv_err_tol,
00181       "The maximum allowed error in the inverted for parameters" );
00182     
00183     CommandLineProcessor::EParseCommandLineReturn
00184       parse_return = clp.parse(argc,argv,&std::cerr);
00185 
00186     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
00187       return parse_return;
00188 
00189     moochoThyraSolver.readParameters(&*out);
00190 
00191     RCP<ParameterList>
00192       paramList = Teuchos::parameterList();
00193     if (paramsFileName.length())
00194       updateParametersFromXmlFile( paramsFileName, paramList.ptr() );
00195     *out << "\nRead in parameters list:\n";
00196     paramList->print(*out, PLPO().indent(2).showTypes(true).showDoc(true));
00197     
00198     //
00199     // B) Create the Thyra::ModelEvaluator object for the ODE
00200     //
00201     
00202     *out << "\nCreate the EpetraExt::DiagonalTransientModel wrapper object ...\n";
00203 
00204     RCP<EpetraExt::DiagonalTransientModel>
00205       epetraStateModel = EpetraExt::diagonalTransientModel(
00206         epetra_comm,
00207         sublist(paramList,"DiagonalTransientModel",true)
00208         );
00209 
00210     if (printValidOptions) {
00211       *out <<"\nepetraStateModel valid options:\n";
00212       epetraStateModel->getValidParameters()->print(
00213         *out, PLPO().indent(2).showTypes(true).showDoc(true)
00214         );
00215     }
00216 
00217     Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
00218     linearSolverBuilder.setParameterList(sublist(paramList,"Stratimikos",true));
00219     RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
00220       W_factory = createLinearSolveStrategy(linearSolverBuilder);
00221 
00222     RCP<Thyra::ModelEvaluator<Scalar> >
00223       stateModel = epetraModelEvaluator(epetraStateModel,W_factory);
00224     
00225     //
00226     // C) Create and initialize the state stepper and state integrator objects
00227     //
00228 
00229     RCP<ParameterList>
00230       nonlinearSolverPL = sublist(paramList,"TimeStepNonlinearSolver",true);
00231     //nonlinearSolverPL->get("Default Tol",1e-3*maxStateError); // Set default if not set
00232     RCP<Rythmos::TimeStepNonlinearSolver<Scalar> >
00233       nonlinearSolver = rcp(new Rythmos::TimeStepNonlinearSolver<Scalar>());
00234     nonlinearSolver->setParameterList(nonlinearSolverPL);
00235 
00236     RCP<Rythmos::StepperBase<Scalar> > stateStepper;
00237     if (useBDF) {
00238       stateStepper = rcp(
00239         new Rythmos::ImplicitBDFStepper<Scalar>(
00240           stateModel, nonlinearSolver
00241           )
00242         );
00243     }
00244     else {
00245       stateStepper = rcp(
00246         new Rythmos::BackwardEulerStepper<Scalar>(
00247           stateModel, nonlinearSolver
00248           )
00249         );
00250     }
00251     stateStepper->setParameterList(sublist(paramList,"Rythmos Stepper",true));
00252 
00253     *out <<"\nstateStepper: " << describe(*stateStepper,verbLevel);
00254     if (printValidOptions) {
00255       *out <<"\nstateStepper valid options:\n";
00256       stateStepper->getValidParameters()->print(
00257         *out, PLPO().indent(2).showTypes(true).showDoc(true)
00258         );
00259     }
00260       
00261     RCP<Rythmos::IntegratorBase<Scalar> > integrator;
00262     {
00263       integrator = Rythmos::controlledDefaultIntegrator<Scalar>(
00264         Rythmos::simpleIntegrationControlStrategy<Scalar>(
00265           sublist(paramList,"Rythmos Integration Control",true)
00266           )
00267         );
00268     }
00269     
00270     const MEB::InArgs<Scalar>
00271       state_ic = stateModel->getNominalValues();
00272     *out << "\nstate_ic: " << describe(state_ic,verbLevel);
00273 
00274     RCP<const Thyra::VectorBase<Scalar> >
00275       p = state_ic.get_p(0);
00276     
00277     stateStepper->setInitialCondition(state_ic);
00278     integrator->setStepper(stateStepper,finalTime);
00279       
00280     if (printValidOptions) {
00281       *out <<"\nintegrator valid options:\n";
00282       integrator->getValidParameters()->print(
00283         *out, PLPO().indent(2).showTypes(true).showDoc(true)
00284         );
00285     }
00286 
00287     //
00288     // D) Generate the matching vectors by running the forward problem through
00289     //    the state integrator
00290     //
00291 
00292     Array<Scalar> observationTimes;
00293     Array<RCP<const Thyra::VectorBase<Scalar> > > stateObservations;
00294     {
00295       const double obs_dt = (finalTime - state_ic.get_t())/numObservationPoints;
00296       *out << "\nCollecting observations at intervals of obs_dt = "
00297            << obs_dt << " ...\n";
00298       OSTab tab(out);
00299       int obs_idx;
00300       double curr_t;
00301       for (
00302         obs_idx = 0, curr_t = state_ic.get_t() + obs_dt; 
00303         obs_idx < numObservationPoints;
00304         ++obs_idx, curr_t += obs_dt
00305         )
00306       {
00307         // Make sure we don't step over final time
00308         curr_t = std::min(curr_t,finalTime);
00309 
00310         *out << "\nCollecting state observation at curr_t = " << curr_t << " ...\n";
00311         observationTimes.push_back(curr_t);
00312         stateObservations.push_back(get_fwd_x(*integrator,curr_t)->clone_v());
00313         if (includesVerbLevel(verbLevel,Teuchos::VERB_EXTREME)) {
00314           *out << "\nx("<<curr_t<<") = "
00315                << Teuchos::describe(*stateObservations.back(),Teuchos::VERB_EXTREME);
00316         }
00317       }
00318     }
00319     
00320     //
00321     // F) Create the forward sensitivity stepper and reduced sensitivity model
00322     //    evaluator
00323     //
00324     
00325     //
00326     // F.1) Create the forward sensitivity stepper
00327     //
00328     
00329     RCP<Rythmos::ForwardSensitivityStepper<Scalar> >
00330       stateAndSensStepper = Rythmos::forwardSensitivityStepper<Scalar>(
00331         stateModel, 0, stateModel->getNominalValues(),
00332         stateStepper, nonlinearSolver
00333         );
00334     // The above call will result in stateStepper and nonlinearSolver being
00335     // cloned.  This helps to ensure consistency between the state and
00336     // sensitivity computations!
00337 
00338     const RCP<const DMVPVS>
00339       s_bar_space = rcp_dynamic_cast<const DMVPVS>(
00340         stateAndSensStepper->getFwdSensModel()->get_x_space()
00341         );
00342 
00343     //
00344     // F.2) Set the initial condition for the state and forward sensitivities
00345     //
00346 
00347     RCP<Thyra::VectorBase<Scalar> > s_bar_init
00348       = createMember(stateAndSensStepper->getFwdSensModel()->get_x_space());
00349     assign( s_bar_init.ptr(), 0.0 );
00350     RCP<Thyra::VectorBase<Scalar> > s_bar_dot_init
00351       = createMember(stateAndSensStepper->getFwdSensModel()->get_x_space());
00352     assign( s_bar_dot_init.ptr(), 0.0 );
00353     // Above, I believe that these are the correct initial conditions for
00354     // s_bar and s_bar_dot given how the EpetraExt::DiagonalTransientModel
00355     // is currently implemented!
00356 
00357     RCP<const Rythmos::StateAndForwardSensitivityModelEvaluator<Scalar> >
00358       stateAndSensModel = stateAndSensStepper->getStateAndFwdSensModel();
00359 
00360     MEB::InArgs<Scalar>
00361       state_and_sens_ic = stateAndSensStepper->getModel()->createInArgs();
00362 
00363     // Copy time, parameters etc.
00364     state_and_sens_ic.setArgs(state_ic);
00365     // Set initial condition for x_bar = [ x; s_bar ]
00366     state_and_sens_ic.set_x(
00367       stateAndSensModel->create_x_bar_vec(state_ic.get_x(),s_bar_init)
00368       );
00369     // Set initial condition for x_bar_dot = [ x_dot; s_bar_dot ]
00370     state_and_sens_ic.set_x_dot(
00371       stateAndSensModel->create_x_bar_vec(state_ic.get_x_dot(),s_bar_dot_init)
00372       );
00373     
00374     *out << "\nstate_and_sens_ic:\n" << describe(state_and_sens_ic,verbLevel);
00375     
00376     stateAndSensStepper->setInitialCondition(state_and_sens_ic);
00377 
00378     //
00379     // F.3) Setup the state-matching response functions
00380     //
00381 
00382     Array<RCP<const Thyra::ModelEvaluator<Scalar> > > matchingFuncs;
00383     Array<Thyra::ModelEvaluatorBase::InArgs<Scalar> > matchingFuncBasePoints;
00384 
00385     for (int obs_idx = 0; obs_idx < numObservationPoints; ++obs_idx ) {
00386       RCP<Thyra::DefaultInverseModelEvaluator<Scalar> >
00387         matchingFunc = Thyra::defaultInverseModelEvaluator(stateModel);
00388       matchingFunc->setParameterList(
00389         sublist(paramList,"DefaultInverseModelEvalautor",true));
00390       matchingFunc->set_observationTarget(stateObservations[obs_idx]);
00391       matchingFuncs.push_back(matchingFunc);
00392       matchingFuncBasePoints.push_back(state_ic);
00393     }
00394 
00395     //
00396     // F.4) Setup the final integrator that will evaluate the reduced response
00397     //      function and it's derivatives
00398     
00399     namespace FSIAMET = Rythmos::ForwardSensitivityIntegratorAsModelEvaluatorTypes;
00400     const RCP<Rythmos::IntegratorBase<Scalar> > nullIntegrator;
00401     RCP<Rythmos::ForwardSensitivityIntegratorAsModelEvaluator<Scalar> >
00402       sensIntegatorAsModelEvaluator =
00403       Rythmos::forwardSensitivityIntegratorAsModelEvaluator<Scalar>(
00404         stateStepper, integrator, stateAndSensStepper, nullIntegrator,
00405         state_and_sens_ic, observationTimes, matchingFuncs,
00406         matchingFuncBasePoints, FSIAMET::RESPONSE_TYPE_SUM
00407         );
00408     sensIntegatorAsModelEvaluator->setParameterList(
00409       sublist(paramList,"ForwardSensitivityIntegratorAsModelEvaluator",true) );
00410     
00411     //
00412     // G) Test the ForwardSensitivityIntegratorAsModelEvaluator object
00413     //
00414 
00415     if (testTransientModelEvaluator) {
00416 
00417       //
00418       // G.1) Test that the base evaluation computes a zero response function.
00419       //
00420 
00421       // Here, we check that if we evaluation p -> g(p) for the base p, then we
00422       // should get zero since the state should match the read observations
00423       // exactly.  We also check that the computed state is exactly the same
00424       // also.  We can cheat and get the final state from stateStepper which is
00425       // used internally.
00426       //
00427 
00428       *out << "\nTest that sensIntegatorAsModelEvaluator computes the right base point ... \n";
00429       {
00430       
00431         RCP<Thyra::VectorBase<Scalar> >
00432           d_hat = createMember(sensIntegatorAsModelEvaluator->get_g_space(0));
00433 
00434         eval_g(
00435           *sensIntegatorAsModelEvaluator,
00436           0, *state_ic.get_p(0),
00437           0, d_hat.ptr()
00438           );
00439 
00440         *out
00441           << "\nd_hat(p) evaluated using sensIntegatorAsModelEvaluator:\n"
00442           << *d_hat;
00443 
00444         // Test that d_hat == 0.0 to all digits since 
00445 
00446         const Scalar d_hat_nrm = norm(*d_hat);
00447 
00448         const bool d_hat_norm_is_zero = (d_hat_nrm == 0.0);
00449 
00450         *out << "\nCheck: ||d_hat|| = " << d_hat_nrm << " == 0 : "
00451              << Thyra::passfail(d_hat_norm_is_zero) << endl;
00452 
00453         if (!d_hat_norm_is_zero)
00454           success = false;
00455 
00456         // Test that:
00457         //
00458         //  get_x(stateStepper,finalTime) == stateObservations[numObservationPoints-1]
00459         //
00460         // to the given precision.
00461 
00462         *out
00463           << "\nChecking that get_x(stateStepper,finalTime) == stateObservations[numObservationPoints-1] ... \n";
00464 
00465         RCP<const Thyra::VectorBase<Scalar> >
00466           x_final = get_x(*stateStepper,observationTimes[numObservationPoints-1]);
00467       
00468         *out << endl;
00469         result = Thyra::testRelNormDiffErr(
00470           "get_x(stateStepper,finalTime)", *x_final,
00471           "stateObservations[numObservationPoints-1]", *stateObservations[numObservationPoints-1],
00472           "errorTol", 0.0, "warningTol", 1e-8,
00473           &OSTab(out).o(), verbLevel
00474           );
00475         if (!result) success = false;
00476 
00477       }
00478 
00479       *out << "\nBase parameters p = " << *p;
00480 
00481       RCP<Thyra::VectorBase<Scalar> >
00482         p_perturb = createMember(p->space());
00483 
00484       V_StV( p_perturb.ptr(), scaleParamsBy, *p );
00485 
00486       *out << "\nPerturbed parameters p_perturb = " << *p_perturb;
00487 
00488       //
00489       // G.2) Test the sensitivity D(d_hat)/d(p) by finite differences
00490       //
00491 
00492       //
00493       *out << "\nCompute analytic D(d_hat)/d(p) at p_perturb ...\n";
00494       //
00495 
00496       MEB::DerivativeMultiVector<Scalar> D_d_hat_D_p_exact;
00497       {
00498 
00499         D_d_hat_D_p_exact = Thyra::create_DgDp_mv(
00500           *sensIntegatorAsModelEvaluator, 0, 0, MEB::DERIV_TRANS_MV_BY_ROW );
00501 
00502         MEB::InArgs<Scalar>
00503           inArgs = sensIntegatorAsModelEvaluator->createInArgs();
00504         inArgs.set_p(0,p_perturb);
00505 
00506         MEB::OutArgs<Scalar>
00507           outArgs = sensIntegatorAsModelEvaluator->createOutArgs();
00508         outArgs.set_DgDp(0,0,D_d_hat_D_p_exact);
00509 
00510         sensIntegatorAsModelEvaluator->evalModel(inArgs,outArgs);
00511       
00512         *out
00513           << "\nAnalytic D(d_hat)/d(p)^T = "
00514           << describe(*D_d_hat_D_p_exact.getMultiVector(),verbLevel);
00515 
00516       }
00517  
00518       //
00519       *out << "\nCompute finite-diff D(d_hat)/d(p) at p_perturb ...\n";
00520       //
00521 
00522       MEB::DerivativeMultiVector<Scalar> D_d_hat_D_p_fd;
00523       {
00524 
00525         RCP<Thyra::DirectionalFiniteDiffCalculator<Scalar> > fdCalc =
00526           Thyra::directionalFiniteDiffCalculator<Scalar>(
00527             sublist(paramList,"DirectionalFiniteDiffCalculator",true)
00528             );
00529 
00530         RCP<Thyra::DefaultFiniteDifferenceModelEvaluator<Scalar> > fdModelEval =
00531           Thyra::defaultFiniteDifferenceModelEvaluator<Scalar>(
00532             sensIntegatorAsModelEvaluator,
00533             fdCalc
00534             );
00535 
00536         MEB::InArgs<Scalar>
00537           inArgs = fdModelEval->createInArgs();
00538         inArgs.set_p(0,p_perturb);
00539       
00540         D_d_hat_D_p_fd = Thyra::create_DgDp_mv(
00541           *fdModelEval, 0, 0, MEB::DERIV_TRANS_MV_BY_ROW
00542           );
00543       
00544         MEB::OutArgs<Scalar>
00545           outArgs = fdModelEval->createOutArgs();
00546         outArgs.set_DgDp( 0, 0, D_d_hat_D_p_fd );
00547       
00548         // Silence the model evaluators that are called.  The fdCal object
00549         // will show all of the inputs and outputs for each call.
00550         stateStepper->setVerbLevel(Teuchos::VERB_NONE);
00551         sensIntegatorAsModelEvaluator->setVerbLevel(Teuchos::VERB_NONE);
00552 
00553         fdModelEval->evalModel(inArgs,outArgs);
00554       
00555         *out
00556           << "\nFinite difference D(d_hat)/d(p)^T = "
00557           << describe(*D_d_hat_D_p_fd.getMultiVector(),verbLevel);
00558       
00559       }
00560     
00561       //
00562       *out << "\nCompare analytic and finite-diff D(d_hat)/d(p)^T ...\n";
00563       //
00564 
00565       RCP<const Thyra::VectorBase<Scalar> >
00566         D_d_hat_D_p_exact_vec = D_d_hat_D_p_exact.getMultiVector()->col(0),
00567         D_d_hat_D_p_fd_vec = D_d_hat_D_p_fd.getMultiVector()->col(0);
00568     
00569       result = Thyra::testRelNormDiffErr(
00570         "D_d_hat_D_p_exact_vec", *D_d_hat_D_p_exact_vec,
00571         "D_d_hat_D_p_fd_vec", *D_d_hat_D_p_fd_vec,
00572         "maxSensError", maxSensError,
00573         "warningTol", 1.0, // Don't warn
00574         &*out, verbLevel
00575         );
00576       if (!result) success = false;
00577 
00578     }
00579 
00580     if (doInverseProblem) {
00581     
00582       //
00583       // H) Setup the MoochoThyraSolver object
00584       //
00585       
00586       *out << "\nSetup the MoochoThyraSolver object ...\n";
00587       
00588       moochoThyraSolver.setParameterList(
00589         sublist(paramList,"MoochoThyraSolver",true) );
00590       
00591       moochoThyraSolver.setModel(sensIntegatorAsModelEvaluator);
00592       
00593       //
00594       // I) Solve the transient inverse problem using MOOCHO.
00595       //
00596       
00597       *out << "\nSetup the transient inverse problem ...\n";
00598       
00599       moochoThyraSolver.readInitialGuess(&*out);
00600       
00601       const MoochoSolver::ESolutionStatus
00602         solution_status = moochoThyraSolver.solve();
00603       
00604       if ( solution_status != MoochoSolver::SOLVE_RETURN_SOLVED )
00605         success = false;
00606       
00607       moochoThyraSolver.writeFinalSolution(&*out);
00608       
00609       //
00610       // J) Check the inverse solution against the true solution
00611       //
00612       
00613       RCP<const Thyra::VectorBase<Scalar> >
00614         p_inv = moochoThyraSolver.getFinalPoint().get_p(0);
00615       
00616       result = Thyra::testRelNormDiffErr(
00617         "p_inv", *p_inv,
00618         "p", *p,
00619         "errorTol", p_inv_err_tol,
00620         "warningTol", p_inv_err_tol * 1e+2,
00621         &OSTab(out).o(), Teuchos::VERB_EXTREME
00622         );
00623       if (!result) success = false;
00624 
00625     }
00626       
00627   }
00628   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success);
00629 
00630   if (success)
00631     *out << "\nEnd Result: TEST PASSED\n";
00632   else
00633     *out << "\nEnd Result: TEST FAILED\n";
00634   
00635   return ( success ? 0 : 1 );
00636 
00637 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends