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

Generated on Wed May 12 21:52:31 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7