MultiPeriodNLPThyraEpetraAdvDiffReactOptMain.cpp

Go to the documentation of this file.
00001 #include "GLpApp_AdvDiffReactOptModelCreator.hpp"
00002 #include "MoochoPack_MoochoThyraSolver.hpp"
00003 #include "Thyra_EpetraModelEvaluator.hpp"
00004 #include "Thyra_EpetraLinearOp.hpp"
00005 #include "Thyra_DefaultClusteredSpmdProductVectorSpace.hpp"
00006 #include "Thyra_DefaultMultiPeriodModelEvaluator.hpp"
00007 #include "Thyra_VectorSpaceTester.hpp"
00008 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
00009 #include "Thyra_DefaultInverseModelEvaluator.hpp"
00010 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00011 #include "Thyra_DefaultMultipliedLinearOp.hpp"
00012 #include "Thyra_TestingTools.hpp"
00013 #include "Teuchos_OpaqueWrapper.hpp"
00014 #include "Teuchos_GlobalMPISession.hpp"
00015 #include "Teuchos_CommandLineProcessor.hpp"
00016 #include "Teuchos_StandardCatchMacros.hpp"
00017 #include "Teuchos_VerboseObject.hpp"
00018 #include "Teuchos_Tuple.hpp"
00019 #include "Teuchos_Utils.hpp"
00020 #include "Teuchos_DefaultComm.hpp"
00021 #ifdef HAVE_MPI
00022 #  include "Teuchos_DefaultMpiComm.hpp"
00023 #  include "Epetra_MpiComm.h"
00024 #else
00025 #  include "Teuchos_DefaultSerialComm.hpp"
00026 #  include "Epetra_SerialComm.h"
00027 #endif
00028 
00029 namespace {
00030 
00031 typedef AbstractLinAlgPack::value_type  Scalar;
00032 
00033 } // namespace
00034 
00035 int main( int argc, char* argv[] )
00036 {
00037 
00038   using Teuchos::rcp;
00039   using Teuchos::rcp_dynamic_cast;
00040   using Teuchos::rcp_implicit_cast;
00041   using Teuchos::null;
00042   using Teuchos::RCP;
00043   using Teuchos::outArg;
00044   using Teuchos::Array;
00045   using Teuchos::tuple;
00046   using Teuchos::ParameterList;
00047   using Teuchos::OpaqueWrapper;
00048   using Teuchos::OSTab;
00049   using Teuchos::CommandLineProcessor;
00050   using Teuchos::toString;
00051   using Thyra::VectorBase;
00052   using Thyra::ProductVectorBase;
00053   typedef Thyra::ModelEvaluatorBase MEB;
00054   typedef Thyra::Index Index;
00055   using Thyra::ModelEvaluator;
00056   using MoochoPack::MoochoSolver;
00057   using MoochoPack::MoochoThyraSolver;
00058 
00059   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00060 
00061   const int procRank = mpiSession.getRank();
00062   const int numProcs = mpiSession.getNProc();
00063 
00064   Teuchos::Time timer("");
00065   
00066   bool result, success = true;
00067 
00068   RCP<Teuchos::FancyOStream>
00069     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00070 
00071   try {
00072   
00073     Stratimikos::DefaultLinearSolverBuilder   lowsfCreator;
00074     GLpApp::AdvDiffReactOptModelCreator     epetraModelCreator;
00075 
00076     // Create the solver object
00077     MoochoThyraSolver solver;
00078 
00079     //
00080     // Get options from the command line
00081     //
00082 
00083     CommandLineProcessor  clp;
00084     clp.throwExceptions(false);
00085     clp.addOutputSetupOptions(true);
00086 
00087     epetraModelCreator.setupCLP(&clp);
00088     lowsfCreator.setupCLP(&clp);
00089     solver.setupCLP(&clp);
00090 
00091     int numProcsPerCluster = -1;
00092     clp.setOption( "num-procs-per-cluster", &numProcsPerCluster,
00093       "Number of processes in a cluster (<=0 means only one cluster)." );
00094     int numPeriodsPerCluster = 1;
00095     clp.setOption( "num-periods-per-cluster", &numPeriodsPerCluster,
00096       "Number of periods in a cluster." );
00097     bool dumpAll = false;
00098     clp.setOption( "dump-all", "no-dump-all", &dumpAll,
00099       "Set to true, then a bunch of debugging output will be created for the clustered vector tests." );
00100     bool skipSolve = false;
00101     clp.setOption( "skip-solve", "no-skip-solve", &skipSolve,
00102       "Temporary flag for skip solve for testing." );
00103     double perturbedParamScaling = 1.0;
00104     clp.setOption( "p-perturb-scaling", &perturbedParamScaling,
00105       "Scaling for perturbed paramters from the initial forward solve." );
00106     bool doMultiPeriod = true;
00107     clp.setOption( "multi-period", "no-multi-period", &doMultiPeriod,
00108       "Do a mulit-period solve or not." );
00109     bool useOuterInverse = true;
00110     clp.setOption( "use-outer-inverse", "use-inner-inverse", &useOuterInverse,
00111       "Determines if the outer inverse model will be used or the inner inverse." );
00112     double periodParamScale = 1.0;
00113     clp.setOption( "period-param-scale", &periodParamScale,
00114       "Sets the scaling factor to scale z[i] from one period to the next." );
00115     bool initialSolveContinuation = false;
00116     clp.setOption( "init-solve-continuation", "init-solve-all-at-once", &initialSolveContinuation,
00117       "Determines if the inital solve is done using continuation or all at once." );
00118     bool useStatelessPeriodModel = false;
00119     clp.setOption( "use-stateless-period-model", "use-statefull-period-model", &useStatelessPeriodModel,
00120       "Determines if a stateless or a statefull period model should be used or not." );
00121     double stateInvError = 1e-8;
00122     clp.setOption( "state-inv-error", &stateInvError,
00123       "The error in the l2 norm of the state inverse solution error." );
00124     double paramInvError = 1e-8;
00125     clp.setOption( "param-inv-error", &paramInvError,
00126       "The error in the l2 norm of the parameter inverse solution error." );
00127 
00128     CommandLineProcessor::EParseCommandLineReturn
00129       parse_return = clp.parse(argc,argv,&std::cerr);
00130 
00131     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
00132       return parse_return;
00133 
00134     lowsfCreator.readParameters(out.get());
00135     solver.readParameters(out.get());
00136 
00137     *out
00138       << "\n***"
00139       << "\n*** NLPThyraEpetraAdvDiffReactOptMain, Global numProcs = "<<numProcs
00140       << "\n***\n";
00141 
00142     int clusterRank = -1;
00143     int numClusters = -1;
00144 
00145 #ifdef HAVE_MPI
00146 
00147     RCP<OpaqueWrapper<MPI_Comm> >
00148       intraClusterMpiComm = Teuchos::opaqueWrapper<MPI_Comm>(MPI_COMM_WORLD),
00149       interClusterMpiComm = Teuchos::null;
00150     
00151     {
00152       if ( numProcsPerCluster <= 0 ) {
00153         *out
00154           << "\nnumProcsPerCluster = " << numProcsPerCluster
00155           << " <= 0: Setting to " << numProcs << "...\n";
00156         numProcsPerCluster = numProcs;
00157       }
00158       *out << "\nCreating communicator for local cluster of "<<numProcsPerCluster<<" processes ...\n";
00159       numClusters = numProcs/numProcsPerCluster;
00160       const int remainingProcs = numProcs%numProcsPerCluster;
00161       TEST_FOR_EXCEPTION(
00162         remainingProcs!=0,std::logic_error
00163         ,"Error, The number of processes per cluster numProcsPerCluster="<<numProcsPerCluster
00164         << " does not divide into the global number of processes numProcs="<<numProcs
00165         << " and instead has remainder="<<remainingProcs<<"!"
00166         );
00167       // Determine which cluster this process is part of and what the global
00168       // process ranges are.
00169       clusterRank = procRank / numProcsPerCluster; // Integer division!
00170       *out << "\nclusterRank = " << clusterRank << "\n";
00171       const int firstClusterProcRank = clusterRank * numProcsPerCluster;
00172       const int lastClusterProcRank = firstClusterProcRank + numProcsPerCluster - 1;
00173       *out << "\nclusterProcRange = ["<<firstClusterProcRank<<","<<lastClusterProcRank<<"]\n";
00174       // Create the communicator for this cluster of processes
00175       *out << "\nCreating intraClusterMpiComm ...";
00176       MPI_Comm rawIntraClusterMpiComm = MPI_COMM_NULL;
00177       MPI_Comm_split(
00178         MPI_COMM_WORLD        // comm
00179         ,clusterRank          // color (will all be put in the same output comm)
00180         ,0                    // key (not important here)
00181         ,&rawIntraClusterMpiComm // newcomm
00182         );
00183       intraClusterMpiComm = Teuchos::opaqueWrapper(rawIntraClusterMpiComm,MPI_Comm_free);
00184       {
00185         *out << "\nintraClusterMpiComm:";
00186         Teuchos::OSTab tab(out);
00187         int rank, size;
00188         MPI_Comm_size(*intraClusterMpiComm,&size);
00189         MPI_Comm_rank(*intraClusterMpiComm,&rank);
00190         *out << "\nsize="<<size;
00191         *out << "\nrank="<<rank;
00192         *out << "\n";
00193       }
00194       // Create the communicator for just the root process in each cluster
00195       *out << "\nCreating interClusterMpiComm ...";
00196       MPI_Comm rawInterClusterMpiComm = MPI_COMM_NULL;
00197       MPI_Comm_split(
00198         MPI_COMM_WORLD                                  // comm
00199         ,procRank==firstClusterProcRank?0:MPI_UNDEFINED // color
00200         ,0                                              // key
00201         ,&rawInterClusterMpiComm                           // newcomm
00202         );
00203       if(rawInterClusterMpiComm!=MPI_COMM_NULL)
00204         interClusterMpiComm = Teuchos::opaqueWrapper(rawInterClusterMpiComm,MPI_Comm_free);
00205       else
00206         interClusterMpiComm = Teuchos::opaqueWrapper(rawInterClusterMpiComm);
00207       {
00208         *out << "\ninterClusterMpiComm:";
00209         Teuchos::OSTab tab(out);
00210         if(*interClusterMpiComm==MPI_COMM_NULL) {
00211           *out << " NULL\n";
00212         }
00213         else {
00214           int rank, size;
00215           MPI_Comm_size(*interClusterMpiComm,&size);
00216           MPI_Comm_rank(*interClusterMpiComm,&rank);
00217           *out << "\nsize="<<size;
00218           *out << "\nrank="<<rank;
00219           *out << "\n";
00220         }
00221       }
00222     }
00223 
00224 #endif
00225 
00226     RCP<Epetra_Comm> comm = Teuchos::null;
00227 #ifdef HAVE_MPI
00228     comm = Teuchos::rcp(new Epetra_MpiComm(*intraClusterMpiComm));
00229     Teuchos::set_extra_data(intraClusterMpiComm, "mpiComm", outArg(comm));
00230 #else
00231     comm = Teuchos::rcp(new Epetra_SerialComm());
00232 #endif
00233     
00234     //
00235     // Create the Thyra::ModelEvaluator object
00236     //
00237     
00238     *out << "\nCreate the GLpApp::AdvDiffReactOptModel wrapper object ...\n";
00239     
00240     RCP<GLpApp::AdvDiffReactOptModel>
00241       epetraModel = epetraModelCreator.createModel(comm);
00242 
00243     *out << "\nCreate the Thyra::LinearOpWithSolveFactory object ...\n";
00244 
00245     RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
00246       lowsFactory = lowsfCreator.createLinearSolveStrategy("");
00247     // ToDo: Set the output stream before calling above!
00248     
00249     *out << "\nCreate the Thyra::EpetraModelEvaluator wrapper object ...\n";
00250     
00251     RCP<Thyra::EpetraModelEvaluator>
00252       epetraThyraModel = rcp(new Thyra::EpetraModelEvaluator()); // Sets default options!
00253     epetraThyraModel->setOStream(out);
00254     epetraThyraModel->initialize(epetraModel,lowsFactory);
00255 
00256     *out
00257       << "\nnx = " << epetraThyraModel->get_x_space()->dim()
00258       << "\nnp = " << epetraThyraModel->get_p_space(0)->dim() << "\n";
00259 
00260     //
00261     // Create the parallel product spaces for x and f
00262     //
00263 
00264     RCP<const Thyra::ProductVectorSpaceBase<Scalar> >
00265       x_bar_space, f_bar_space;
00266     
00267 #ifdef HAVE_MPI
00268 
00269     // For now just build and test these vector spaces if we are not doing a
00270     // solve!  We have a lot more work to do to the "Clustered" support
00271     // software before this will work for a solve.
00272     
00273     if (skipSolve) {
00274       
00275       *out << "\nCreate block parallel vector spaces for multi-period model.x and model.f ...\n";
00276       RCP<const Teuchos::Comm<Index> >
00277         intraClusterComm = rcp(new Teuchos::MpiComm<Index>(intraClusterMpiComm)),
00278         interClusterComm = Teuchos::createMpiComm<Index>(interClusterMpiComm);
00279       x_bar_space = Teuchos::rcp(
00280         new Thyra::DefaultClusteredSpmdProductVectorSpace<Scalar>(
00281           intraClusterComm
00282           ,0 // clusterRootRank
00283           ,interClusterComm
00284           ,1 // numBlocks
00285           ,tuple<RCP<const Thyra::VectorSpaceBase<Scalar> > >(
00286             epetraThyraModel->get_x_space()
00287             ).getRawPtr()
00288           )
00289         );
00290       f_bar_space = Teuchos::rcp(
00291         new Thyra::DefaultClusteredSpmdProductVectorSpace<Scalar>(
00292           intraClusterComm
00293           ,0 // clusterRootRank
00294           ,interClusterComm
00295           ,1 // numBlocks
00296           ,tuple<RCP<const Thyra::VectorSpaceBase<Scalar> > >(
00297             epetraThyraModel->get_f_space()
00298             ).getRawPtr()
00299           )
00300         );
00301       
00302       Thyra::VectorSpaceTester<Scalar> vectorSpaceTester;
00303       vectorSpaceTester.show_all_tests(true);
00304       vectorSpaceTester.dump_all(dumpAll);
00305 
00306 #ifdef RTOPPACK_SPMD_APPLY_OP_DUMP
00307       RTOpPack::show_mpi_apply_op_dump = dumpAll;
00308 #endif
00309 #ifdef THYRA_SPMD_VECTOR_BASE_DUMP
00310       Thyra::SpmdVectorBase<Scalar>::show_dump = dumpAll;
00311 #endif
00312 
00313       *out << "\nTesting the vector space x_bar_space ...\n";
00314       result = vectorSpaceTester.check(*x_bar_space,OSTab(out).get());
00315       if(!result) success = false;
00316 
00317       *out << "\nTesting the vector space f_bar_space ...\n";
00318       result = vectorSpaceTester.check(*f_bar_space,OSTab(out).get());
00319       if(!result) success = false;
00320       
00321       RCP<const VectorBase<Scalar> >
00322         x0 = epetraThyraModel->getNominalValues().get_x();
00323       double nrm_x0;
00324       
00325 /* 2008/02/21: rabartl: I have commented this out since it is causing an MPI error?
00326 
00327       *out << "\nTiming a global reduction across just this cluster: ||x0||_1 = ";
00328       timer.start(true);
00329       nrm_x0 = Thyra::norm_1(*x0);
00330       *out << nrm_x0 << "\n";
00331       timer.stop();
00332       *out << "\n    time = " << timer.totalElapsedTime() << " seconds\n";
00333       
00334       *out << "\nTiming a global reduction across the entire set of processes: ||x0||_1 = ";
00335       timer.start(true);
00336       RTOpPack::ROpNorm1<Scalar> norm_1_op;
00337       RCP<RTOpPack::ReductTarget> norm_1_targ = norm_1_op.reduct_obj_create();
00338       const Teuchos::RCP<const Teuchos::Comm<Teuchos_Index> > comm
00339         = Teuchos::DefaultComm<Index>::getComm();
00340       const Teuchos::ArrayView<const Teuchos::Ptr<const VectorBase<Scalar> > >
00341         vecs = Teuchos::tuple(Teuchos::ptrInArg(*x0));
00342       Teuchos::dyn_cast<const Thyra::SpmdVectorBase<Scalar> >(*x0).applyOpImplWithComm(
00343         comm.ptr(),
00344         norm_1_op, vecs, Teuchos::null, norm_1_targ.ptr(),
00345         0, -1, 0
00346         );
00347       nrm_x0 = norm_1_op(*norm_1_targ);
00348       *out << nrm_x0 << "\n";
00349       timer.stop();
00350       *out << "\n    time = " << timer.totalElapsedTime() << " seconds\n";
00351 
00352 */
00353 
00354 #ifdef RTOPPACK_SPMD_APPLY_OP_DUMP
00355       RTOpPack::show_mpi_apply_op_dump = false;
00356 #endif
00357 #ifdef THYRA_SPMD_VECTOR_BASE_DUMP
00358       Thyra::SpmdVectorBase<Scalar>::show_dump = false;
00359 #endif
00360 
00361     }
00362 
00363 #endif // HAVE_MPI
00364     
00365     if(skipSolve) {
00366 
00367       if(success)
00368         *out << "\nEnd Result: TEST PASSED" << endl;
00369       else
00370         *out << "\nEnd Result: TEST FAILED" << endl;
00371 
00372       return ( success ? 0 : 1 );
00373 
00374     }
00375 
00376     const int N = numPeriodsPerCluster;
00377 
00378     Array<RCP<Thyra::ModelEvaluator<Scalar> > >
00379       inverseThyraModels(N);
00380     if (useOuterInverse) {
00381       *out << "\nUsing Thyra::DefaultInverseModelEvaluator for the objective function ...\n";
00382       if ( useStatelessPeriodModel ) {
00383         *out << "\nBuilding a single Thyra::DefaultInverseModelEvaluator object where the matching vector will be maintained externally ...\n";
00384       }
00385       else {
00386         *out << "\nBuilding multiple Thyra::DefaultInverseModelEvaluator objects where the matching vector is held internally ...\n";
00387       }
00388       RCP<ParameterList> invMEPL = Teuchos::parameterList();
00389       invMEPL->set( "Observation Multiplier", 1.0 );
00390       invMEPL->set( "Parameter Multiplier", epetraModel->getDataPool()->getbeta() );
00391       if ( useStatelessPeriodModel )
00392         invMEPL->set( "Observation Target as Parameter", true );
00393       RCP<const Thyra::EpetraLinearOp>
00394         H = Thyra::epetraLinearOp(epetraModel->getDataPool()->getH(),"H"),
00395         R = Thyra::epetraLinearOp(epetraModel->getDataPool()->getR(),"R");
00396       RCP<const Thyra::MultiVectorBase<Scalar> >
00397         B_bar = Thyra::create_MultiVector(
00398           epetraModel->get_B_bar(),
00399           R->spmdRange()
00400           );
00401       RCP<const Thyra::LinearOpBase<Scalar> >
00402         R_bar = Thyra::multiply<Scalar>(Thyra::adjoint<Scalar>(B_bar),R,B_bar);
00403       for ( int i = 0; i < N; ++i ) {
00404         if ( ( useStatelessPeriodModel && i==0 ) || !useStatelessPeriodModel ) {
00405           RCP<Thyra::DefaultInverseModelEvaluator<Scalar> >
00406             _inverseThyraModel = Thyra::defaultInverseModelEvaluator<Scalar>(
00407               epetraThyraModel );
00408           _inverseThyraModel->setParameterList(invMEPL);
00409           _inverseThyraModel->set_observationMatchWeightingOp(H);
00410           _inverseThyraModel->set_parameterRegularizationWeightingOp(R_bar);
00411           inverseThyraModels[i] = _inverseThyraModel;
00412         }
00413         else {
00414 #ifdef TEUCHOS_DEBUG
00415           TEST_FOR_EXCEPT( ! ( useStatelessPeriodModel && i > 0 ) );
00416 #endif
00417           inverseThyraModels[i] = inverseThyraModels[0];
00418         }
00419       }
00420     }
00421     else {
00422       *out << "\nUsing built-in inverse objective function ...\n";
00423       TEST_FOR_EXCEPTION(
00424         N != 1, std::logic_error,
00425         "Error, you can't have N = "<<N<<" > 1\n"
00426         "and be using an internal inverse objective!" );
00427       inverseThyraModels[0] = epetraThyraModel;
00428     }
00429 
00430     const int p_index = 0;
00431     const int z_index = 1;
00432     const int z_p_index = 1; // Index of the reaction rate parameter parameter subvector
00433     const int z_x_index = 2; // Index of the state matching subvector parameter
00434     Array<int> z_indexes;
00435     if (useStatelessPeriodModel)
00436       z_indexes = tuple<int>(z_p_index, z_x_index);
00437     else
00438       z_indexes = tuple<int>(z_p_index);
00439     Array<Array<RCP<const VectorBase<Scalar> > > > z;
00440     const int g_index = ( useOuterInverse ? 1 : 0 );
00441     Array<Scalar> weights;
00442     RCP<VectorBase<Scalar> >
00443       z_base = inverseThyraModels[0]->getNominalValues().get_p(z_index)->clone_v();
00444     *out << "\nz_base =\n" << Teuchos::describe(*z_base,Teuchos::VERB_EXTREME);
00445     Scalar scale_z_i = 1.0;
00446     for ( int i = 0; i < N; ++i ) {
00447       weights.push_back(1.0);
00448       RCP<VectorBase<Scalar> > z_i = z_base->clone_v();
00449       Vt_S( &*z_i, scale_z_i );
00450       *out << "\nz["<<i<<"] =\n" << Teuchos::describe(*z_i,Teuchos::VERB_EXTREME);
00451       if ( useStatelessPeriodModel ) {
00452         z.push_back(
00453           tuple<RCP<const VectorBase<Scalar> > >(
00454             z_i,
00455             null // We will set this again later!
00456             )
00457           );
00458       }
00459       else {
00460         z.push_back(
00461           tuple<RCP<const VectorBase<Scalar> > >(z_i)
00462           );
00463       }
00464       scale_z_i *= periodParamScale;
00465     }
00466 
00467     RCP<Thyra::ModelEvaluator<Scalar> >
00468       thyraModel = inverseThyraModels[0];
00469 
00470     if (doMultiPeriod) {
00471       thyraModel =
00472         rcp(
00473           new Thyra::DefaultMultiPeriodModelEvaluator<Scalar>(
00474             N, inverseThyraModels, z_indexes, z, g_index, weights,
00475             x_bar_space, f_bar_space
00476             )
00477           );
00478     }
00479     
00480     MoochoSolver::ESolutionStatus solution_status;
00481 
00482     //
00483     *out << "\n***\n*** Solving the initial forward problem\n***\n";
00484     //
00485 
00486     // Set the solve mode to solve the forward problem
00487     solver.setSolveMode(MoochoThyraSolver::SOLVE_MODE_FORWARD);
00488 
00489     // Save the solution for model.x and model.p to be used later
00490     RCP<const VectorBase<Scalar> >
00491       x_opt, // Will be set below
00492       x_init, // Will be set below
00493       p_opt = inverseThyraModels[0]->getNominalValues().get_p(0)->clone_v();
00494 
00495     *out << "\np_opt =\n" << Teuchos::describe(*p_opt,Teuchos::VERB_EXTREME);
00496     
00497     if ( initialSolveContinuation ) {
00498       
00499       *out << "\nSolving individual period problems one at time using continuation ...\n";
00500 
00501       RCP<ProductVectorBase<Scalar> >
00502         x_opt_prod = rcp_dynamic_cast<ProductVectorBase<Scalar> >(
00503          createMember( thyraModel->get_x_space() ), true );
00504 
00505       RCP<const VectorBase<Scalar> > period_x;
00506 
00507       for ( int i = 0; i < N; ++i ) {
00508 
00509         *out << "\nSolving period i = " << i << " using guess from last period ...\n";
00510       
00511         // Set the deliminator for the output files!
00512         solver.getSolver().set_output_context("fwd-init-"+toString(i));
00513 
00514         // Set the period model
00515         solver.setModel(inverseThyraModels[i]);
00516         
00517         // Set the initial guess and the parameter values
00518         MEB::InArgs<Scalar> initialGuess = inverseThyraModels[i]->createInArgs();
00519         initialGuess.set_p(z_index,z[i][0]->clone_v());
00520         initialGuess.set_p(p_index,p_opt->clone_v());
00521         if ( i == 0 ) {
00522           // For the first period just use whatever initial guess is built
00523           // into the model
00524         }
00525         else {
00526           // Set the final solution for x from the last period!
00527           initialGuess.set_x(period_x);
00528         }
00529         solver.setInitialGuess(initialGuess);
00530 
00531         // Solve the period model
00532         solution_status = solver.solve();
00533         TEST_FOR_EXCEPT( solution_status != MoochoSolver::SOLVE_RETURN_SOLVED );
00534 
00535         // Save the final solution for the next period!
00536         period_x = solver.getFinalPoint().get_x()->clone_v();
00537         assign( &*x_opt_prod->getNonconstVectorBlock(i), *period_x );
00538         if ( useStatelessPeriodModel )
00539           z[i][1] = period_x->clone_v(); // This is our matching vector!
00540         
00541       }
00542 
00543       x_opt = x_opt_prod;
00544       x_init = x_opt->clone_v();
00545 
00546       if ( useStatelessPeriodModel ) {
00547         rcp_dynamic_cast<Thyra::DefaultMultiPeriodModelEvaluator<Scalar> >(
00548           thyraModel
00549           )->reset_z(z);
00550       }
00551 
00552     }
00553     else {
00554 
00555       *out << "\nSolving all periods simultaniously ...\n";
00556       
00557       // Set the deliminator for the output files!
00558       solver.getSolver().set_output_context("fwd-init");
00559       
00560       // Set the model
00561       solver.setModel(thyraModel);
00562       
00563       // Set the initial guess from files (if specified on commandline)
00564       solver.readInitialGuess(out.get());
00565       
00566       // Solve the initial forward problem
00567       solution_status = solver.solve();
00568       TEST_FOR_EXCEPT( solution_status != MoochoSolver::SOLVE_RETURN_SOLVED );
00569 
00570       // Save the solution for model.x and model.p to be used later
00571       x_opt = solver.getFinalPoint().get_x()->clone_v();
00572       x_init = solver.getFinalPoint().get_x()->clone_v();
00573       
00574     }
00575     
00576     //
00577     *out << "\n***\n*** Solving the perturbed forward problem\n***\n";
00578     //
00579     
00580     // Set the deliminator for the output files!
00581     solver.getSolver().set_output_context("fwd");
00582     
00583     // Set the solve mode to solve the forward problem
00584     solver.setSolveMode(MoochoThyraSolver::SOLVE_MODE_FORWARD);
00585     
00586     // Set the model
00587     solver.setModel(thyraModel);
00588     
00589     // Set the initial guess and the perturbed parameters
00590     RCP<VectorBase<Scalar> >
00591       p_init = p_opt->clone_v();
00592     {
00593       MEB::InArgs<Scalar> initialGuess = thyraModel->createInArgs();
00594       initialGuess.setArgs(thyraModel->getNominalValues());
00595       initialGuess.set_x(x_init);
00596       Thyra::Vt_S(&*p_init,perturbedParamScaling);
00597       initialGuess.set_p(0,p_init);
00598       //*out << "\nInitial Guess:\n" << Teuchos::describe(initialGuess,Teuchos::VERB_EXTREME);
00599       solver.setInitialGuess(initialGuess);
00600     }
00601 
00602     // Solve the perturbed forward problem
00603     solution_status = solver.solve();
00604     
00605     // Save the solution for model.x and model.p to be used later
00606     x_init = solver.getFinalPoint().get_x()->clone_v();
00607     p_init = solver.getFinalPoint().get_p(0)->clone_v();
00608 
00609     *out
00610       << "\nrelVectorErr(x_perturb,x_opt) = " << Thyra::relVectorErr(*x_init,*x_opt)
00611       << "\nrelVectorErr(p_perturb,p_opt) = " << Thyra::relVectorErr(*p_init,*p_opt)
00612       << "\n";
00613     
00614     //
00615     *out << "\n***\n*** Solving the perturbed inverse problem\n***\n";
00616     //
00617 
00618     // Set the deliminator for the output files!
00619     solver.getSolver().set_output_context("inv");
00620 
00621     //TEST_FOR_EXCEPT("ToDo: We need to use the DefaultInverseModelEvaluator to set the matching vector correctly!");
00622 
00623     // Set the matching vector
00624     if ( N > 1 ) {
00625       TEST_FOR_EXCEPTION(
00626         !useOuterInverse, std::logic_error,
00627         "Error, if N > 1, you have to use the outer inverse objective function\n"
00628         "since each target vector will be different!" );
00629       RCP<const ProductVectorBase<Scalar> >
00630         x_opt_prod = rcp_dynamic_cast<const ProductVectorBase<Scalar> >(
00631           rcp_implicit_cast<const VectorBase<Scalar> >(x_opt), true
00632           ); // This cast can *not* fail!
00633       for ( int i = 0; i < N; ++i ) {
00634         rcp_dynamic_cast<Thyra::DefaultInverseModelEvaluator<Scalar> >(
00635           inverseThyraModels[i], true
00636           )->set_observationTarget(x_opt_prod->getVectorBlock(i));
00637       }
00638     }
00639     else if ( 1 == N ) {
00640       RCP<const ProductVectorBase<Scalar> >
00641         x_opt_prod = rcp_dynamic_cast<const ProductVectorBase<Scalar> >(
00642           rcp_implicit_cast<const VectorBase<Scalar> >(x_opt)
00643           ); // This cast can fail!
00644       RCP<const VectorBase<Scalar> >
00645         x_opt_i =
00646         ( !is_null(x_opt_prod)
00647           ? x_opt_prod->getVectorBlock(0)
00648           : rcp_implicit_cast<const VectorBase<Scalar> >(x_opt)
00649           );
00650       if (useOuterInverse) {
00651         rcp_dynamic_cast<Thyra::DefaultInverseModelEvaluator<Scalar> >(
00652           inverseThyraModels[0], true
00653           )->set_observationTarget(x_opt_i);
00654       }
00655       else {
00656         epetraModel->set_q(
00657           Thyra::get_Epetra_Vector(*epetraModel->get_x_map(), x_opt_i)
00658           );
00659       }
00660     }
00661     else {
00662       TEST_FOR_EXCEPT("Error, should not get here!");
00663     }
00664     
00665     // Set the solve mode to solve the inverse problem
00666     solver.setSolveMode(MoochoThyraSolver::SOLVE_MODE_OPTIMIZE);
00667    
00668     // Set the model
00669     solver.setModel(thyraModel);
00670 
00671     // Set the initial guess for model.x and model.p
00672     {
00673       MEB::InArgs<Scalar> initialGuess = thyraModel->createInArgs();
00674       initialGuess.setArgs(thyraModel->getNominalValues());
00675       initialGuess.set_x(x_init);
00676       initialGuess.set_p(0,p_init);
00677       //*out << "\nInitial Guess:\n" << Teuchos::describe(initialGuess,Teuchos::VERB_EXTREME);
00678       solver.setInitialGuess(initialGuess);
00679     }
00680     
00681     // Solve the inverse problem
00682     solution_status = solver.solve();
00683     TEST_FOR_EXCEPT( solution_status != MoochoSolver::SOLVE_RETURN_SOLVED );
00684 
00685     //
00686     *out << "\n***\n*** Testing the error in the inversion\n***\n";
00687     //
00688     
00689     // Get the inverted for solution and compare it to the optimal solution
00690 
00691     RCP<const VectorBase<Scalar> >
00692       x_inv = solver.getFinalPoint().get_x(),
00693       p_inv = solver.getFinalPoint().get_p(0);
00694 
00695     *out << "\np_opt =\n" << Teuchos::describe(*p_opt,Teuchos::VERB_EXTREME);
00696     *out << "\np_inv =\n" << Teuchos::describe(*p_inv,Teuchos::VERB_EXTREME);
00697 
00698     const Scalar
00699       x_err = Thyra::relVectorErr( *x_inv, *x_opt ),
00700       p_err = Thyra::relVectorErr( *p_inv, *p_opt );
00701 
00702     const bool
00703       x_test_passed = ( x_err <= stateInvError ),
00704       p_test_passed = ( p_err <= paramInvError );
00705 
00706     *out
00707       << "\nrelVectorErr(x_inv,x_opt) = " << x_err << " <= " << stateInvError
00708       << " : " << Thyra::passfail(x_test_passed)
00709       << "\nrelVectorErr(p_inv,p_opt) = " << p_err << " <= " << paramInvError
00710       << " : " << Thyra::passfail(p_test_passed)
00711       << "\n";
00712     
00713     // Write the final solution
00714     solver.writeFinalSolution(out.get());
00715     
00716     // Write the final parameters to file
00717     lowsfCreator.writeParamsFile(*lowsFactory);
00718     solver.writeParamsFile();
00719     
00720     TEST_FOR_EXCEPT(
00721       solution_status != MoochoSolver::SOLVE_RETURN_SOLVED
00722       );
00723     
00724   }
00725   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
00726 
00727   if(success)
00728     *out << "\nEnd Result: TEST PASSED" << endl;
00729   else
00730     *out << "\nEnd Result: TEST FAILED" << endl;
00731   
00732   return ( success ? 0 : 1 );
00733   
00734 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:11:00 2011 for MOOCHO (Single Doxygen Collection) by  doxygen 1.6.3