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