MultiPeriodNLPThyraEpetraAdvDiffReactOptMain.cpp

Go to the documentation of this file.
00001 #include "GLpApp_AdvDiffReactOptModelCreator.hpp"
00002 #include "MoochoPack_ThyraModelEvaluatorSolver.hpp"
00003 #include "Thyra_EpetraModelEvaluator.hpp"
00004 #include "Thyra_SpmdMultiVectorFileIO.hpp"
00005 #include "Thyra_DefaultClusteredSpmdProductVectorSpace.hpp"
00006 #include "Thyra_DefaultMultiPeriodModelEvaluator.hpp"
00007 #include "Thyra_VectorSpaceTester.hpp"
00008 #include "Thyra_DefaultRealLinearSolverBuilder.hpp"
00009 #include "RTOpPack_MPI_apply_op_decl.hpp"
00010 #include "Teuchos_OpaqueWrapper.hpp"
00011 #include "Teuchos_GlobalMPISession.hpp"
00012 #include "Teuchos_CommandLineProcessor.hpp"
00013 #include "Teuchos_StandardCatchMacros.hpp"
00014 #include "Teuchos_VerboseObject.hpp"
00015 #include "Teuchos_arrayArg.hpp"
00016 #include "Teuchos_DefaultComm.hpp"
00017 #ifdef HAVE_MPI
00018 #  include "Teuchos_DefaultMpiComm.hpp"
00019 #  include "Epetra_MpiComm.h"
00020 #else
00021 #  include "Teuchos_DefaultSerialComm.hpp"
00022 #  include "Epetra_SerialComm.h"
00023 #endif
00024 
00025 namespace {
00026 
00027 typedef AbstractLinAlgPack::value_type  Scalar;
00028 
00029 } // namespace
00030 
00031 int main( int argc, char* argv[] )
00032 {
00033   using Teuchos::rcp;
00034   using Teuchos::null;
00035   using Teuchos::RefCountPtr;
00036   using Teuchos::OpaqueWrapper;
00037   using Teuchos::OSTab;
00038   using MoochoPack::MoochoSolver;
00039   using MoochoPack::ThyraModelEvaluatorSolver;
00040   using Teuchos::CommandLineProcessor;
00041   typedef Thyra::ModelEvaluatorBase MEB;
00042   typedef Thyra::Index Index;
00043 
00044   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00045 
00046   const int procRank = mpiSession.getRank();
00047   const int numProcs = mpiSession.getNProc();
00048 
00049   Teuchos::Time timer("");
00050   
00051   bool result, success = true;
00052 
00053   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00054     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00055 
00056   try {
00057   
00058     Thyra::DefaultRealLinearSolverBuilder   lowsfCreator;
00059     GLpApp::AdvDiffReactOptModelCreator     epetraModelCreator;
00060 
00061     // Create the solver object
00062     ThyraModelEvaluatorSolver solver;
00063 
00064     //
00065     // Get options from the command line
00066     //
00067 
00068     int                 numProcsPerCluster     = -1;
00069     double              perturbedParamScaling  = 1.0;
00070     bool                dumpAll                = false;
00071     bool                skipSolve              = false;
00072 
00073     CommandLineProcessor  clp;
00074     clp.throwExceptions(false);
00075     clp.addOutputSetupOptions(true);
00076 
00077     epetraModelCreator.setupCLP(&clp);
00078     lowsfCreator.setupCLP(&clp);
00079     solver.setupCLP(&clp);
00080 
00081     clp.setOption( "num-procs-per-cluster", &numProcsPerCluster, "Number of processes in a cluster (<=0 means only one cluster)." );
00082     clp.setOption( "p-perturb-scaling", &perturbedParamScaling, "Scaling for perturbed paramters from the initial forward solve." );
00083     clp.setOption( "dump-all", "no-dump-all", &dumpAll, "Set to true, then a bunch of debugging output will be created." );
00084     clp.setOption( "skip-solve", "no-skip-solve", &skipSolve, "Temporary flag for skip solve for testing." );
00085 
00086     CommandLineProcessor::EParseCommandLineReturn
00087       parse_return = clp.parse(argc,argv,&std::cerr);
00088 
00089     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
00090       return parse_return;
00091 
00092     lowsfCreator.readParameters(out.get());
00093 
00094     *out
00095       << "\n***"
00096       << "\n*** NLPThyraEpetraAdvDiffReactOptMain, Global numProcs = "<<numProcs
00097       << "\n***\n";
00098 
00099     int clusterRank = -1;
00100     int numClusters = -1;
00101 #ifdef HAVE_MPI
00102     RefCountPtr<OpaqueWrapper<MPI_Comm> >
00103       intraClusterMpiComm = Teuchos::opaqueWrapper<MPI_Comm>(MPI_COMM_WORLD),
00104       interClusterMpiComm = Teuchos::null;
00105     //if( numProcsPerCluster > 0 ) {
00106     if(1) {
00107       *out << "\nCreating communicator for local cluster of "<<numProcsPerCluster<<" processes ...\n";
00108       numClusters = numProcs/numProcsPerCluster;
00109       const int remainingProcs = numProcs%numProcsPerCluster;
00110       TEST_FOR_EXCEPTION(
00111         remainingProcs!=0,std::logic_error
00112         ,"Error, The number of processes per cluster numProcsPerCluster="<<numProcsPerCluster
00113         << " does not divide into the global number of processes numProcs="<<numProcs
00114         << " and instead has remainder="<<remainingProcs<<"!"
00115         );
00116       // Determine which cluster this process is part of and what the global
00117       // process ranges are.
00118       clusterRank = procRank / numProcsPerCluster; // Integer division!
00119       *out << "\nclusterRank = " << clusterRank << "\n";
00120       const int firstClusterProcRank = clusterRank * numProcsPerCluster;
00121       const int lastClusterProcRank = firstClusterProcRank + numProcsPerCluster - 1;
00122       *out << "\nclusterProcRange = ["<<firstClusterProcRank<<","<<lastClusterProcRank<<"]\n";
00123       // Create the communicator for this cluster of processes
00124       *out << "\nCreating intraClusterMpiComm ...";
00125       MPI_Comm rawIntraClusterMpiComm = MPI_COMM_NULL;
00126       MPI_Comm_split(
00127         MPI_COMM_WORLD        // comm
00128         ,clusterRank          // color (will all be put in the same output comm)
00129         ,0                    // key (not important here)
00130         ,&rawIntraClusterMpiComm // newcomm
00131         );
00132       intraClusterMpiComm = Teuchos::opaqueWrapper(rawIntraClusterMpiComm,MPI_Comm_free);
00133       if(1) {
00134         *out << "\nintraClusterMpiComm:";
00135         Teuchos::OSTab tab(out);
00136         int rank, size;
00137         MPI_Comm_size(*intraClusterMpiComm,&size);
00138         MPI_Comm_rank(*intraClusterMpiComm,&rank);
00139         *out << "\nsize="<<size;
00140         *out << "\nrank="<<rank;
00141         *out << "\n";
00142       }
00143       // Create the communicator for just the root process in each cluster
00144       *out << "\nCreating interClusterMpiComm ...";
00145       MPI_Comm rawInterClusterMpiComm = MPI_COMM_NULL;
00146       MPI_Comm_split(
00147         MPI_COMM_WORLD                                  // comm
00148         ,procRank==firstClusterProcRank?0:MPI_UNDEFINED // color
00149         ,0                                              // key
00150         ,&rawInterClusterMpiComm                           // newcomm
00151         );
00152       if(rawInterClusterMpiComm!=MPI_COMM_NULL)
00153         interClusterMpiComm = Teuchos::opaqueWrapper(rawInterClusterMpiComm,MPI_Comm_free);
00154       else
00155         interClusterMpiComm = Teuchos::opaqueWrapper(rawInterClusterMpiComm);
00156       if(1) {
00157         *out << "\ninterClusterMpiComm:";
00158         Teuchos::OSTab tab(out);
00159         if(*interClusterMpiComm==MPI_COMM_NULL) {
00160           *out << " NULL\n";
00161         }
00162         else {
00163           int rank, size;
00164           MPI_Comm_size(*interClusterMpiComm,&size);
00165           MPI_Comm_rank(*interClusterMpiComm,&rank);
00166           *out << "\nsize="<<size;
00167           *out << "\nrank="<<rank;
00168           *out << "\n";
00169         }
00170       }
00171     }
00172 #endif
00173 
00174     Teuchos::RefCountPtr<Epetra_Comm> comm = Teuchos::null;
00175 #ifdef HAVE_MPI
00176     comm = Teuchos::rcp(new Epetra_MpiComm(*intraClusterMpiComm));
00177     Teuchos::set_extra_data(intraClusterMpiComm,"mpiComm",&comm);
00178 #else
00179     comm = Teuchos::rcp(new Epetra_SerialComm());
00180 #endif
00181     
00182     //
00183     // Create the Thyra::ModelEvaluator object
00184     //
00185     
00186     *out << "\nCreate the GLpApp::AdvDiffReactOptModel wrapper object ...\n";
00187     
00188     Teuchos::RefCountPtr<GLpApp::AdvDiffReactOptModel>
00189       epetraModel = epetraModelCreator.createModel(comm);
00190 
00191     *out << "\nCreate the Thyra::LinearOpWithSolveFactory object ...\n";
00192 
00193     Teuchos::RefCountPtr<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
00194       lowsFactory = lowsfCreator.createLinearSolveStrategy("");
00195     // ToDo: Set the output stream before calling above!
00196     
00197     *out << "\nCreate the Thyra::EpetraModelEvaluator wrapper object ...\n";
00198     
00199     Teuchos::RefCountPtr<Thyra::EpetraModelEvaluator>
00200       epetraThyraModel = rcp(new Thyra::EpetraModelEvaluator()); // Sets default options!
00201     epetraThyraModel->setOStream(out);
00202     epetraThyraModel->initialize(epetraModel,lowsFactory);
00203 
00204     *out
00205       << "\nnx = " << epetraThyraModel->get_x_space()->dim()
00206       << "\nnp = " << epetraThyraModel->get_p_space(0)->dim() << "\n";
00207 
00208     Teuchos::RefCountPtr<Thyra::ModelEvaluator<Scalar> >
00209       thyraModel = epetraThyraModel;
00210     
00211 #ifdef HAVE_MPI
00212 
00213     //if( numClusters > 0 ) {
00214     if(1) {
00215       
00216       *out << "\nCreate block parallel vector spaces for multi-period model.x and model.f ...\n";
00217       Teuchos::RefCountPtr<const Teuchos::Comm<Index> >
00218         intraClusterComm = rcp(new Teuchos::MpiComm<Index>(intraClusterMpiComm)),
00219         interClusterComm = Teuchos::createMpiComm<Index>(interClusterMpiComm);
00220       Teuchos::RefCountPtr<Thyra::DefaultClusteredSpmdProductVectorSpace<Scalar> >
00221         x_bar_space = Teuchos::rcp(
00222           new Thyra::DefaultClusteredSpmdProductVectorSpace<Scalar>(
00223             intraClusterComm
00224             ,0 // clusterRootRank
00225             ,interClusterComm
00226             ,1 // numBlocks
00227             ,Teuchos::arrayArg<Teuchos::RefCountPtr<const Thyra::VectorSpaceBase<Scalar> > >(
00228               epetraThyraModel->get_x_space()
00229               )()
00230             )
00231           ),
00232         f_bar_space = Teuchos::rcp(
00233           new Thyra::DefaultClusteredSpmdProductVectorSpace<Scalar>(
00234             intraClusterComm
00235             ,0 // clusterRootRank
00236             ,interClusterComm
00237             ,1 // numBlocks
00238             ,Teuchos::arrayArg<Teuchos::RefCountPtr<const Thyra::VectorSpaceBase<Scalar> > >(
00239               epetraThyraModel->get_f_space()
00240               )()
00241             )
00242           );
00243 
00244       Thyra::VectorSpaceTester<Scalar> vectorSpaceTester;
00245       vectorSpaceTester.show_all_tests(true);
00246       vectorSpaceTester.dump_all(dumpAll);
00247 
00248 #ifdef RTOPPACK_SPMD_APPLY_OP_DUMP
00249       RTOpPack::show_mpi_apply_op_dump = dumpAll;
00250 #endif
00251 #ifdef THYRA_SPMD_VECTOR_BASE_DUMP
00252       Thyra::SpmdVectorBase<Scalar>::show_dump = dumpAll;
00253 #endif
00254 
00255       *out << "\nTesting the vector space x_bar_space ...\n";
00256       result = vectorSpaceTester.check(*x_bar_space,&*OSTab(out).getOStream());
00257       if(!result) success = false;
00258 
00259       *out << "\nTesting the vector space f_bar_space ...\n";
00260       result = vectorSpaceTester.check(*f_bar_space,&*OSTab(out).getOStream());
00261       if(!result) success = false;
00262       
00263       Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> >
00264         x0 = epetraThyraModel->getNominalValues().get_x();
00265       double nrm_x0;
00266       
00267       *out << "\nTiming a global reduction across just this cluster: ||x0||_1 = ";
00268       timer.start(true);
00269       nrm_x0 = Thyra::norm_1(*x0);
00270       *out << nrm_x0 << "\n";
00271       timer.stop();
00272       *out << "\n    time = " << timer.totalElapsedTime() << " seconds\n";
00273       
00274       *out << "\nTiming a global reduction across the entire set of processes: ||x0||_1 = ";
00275       timer.start(true);
00276       RTOpPack::ROpNorm1<Scalar> norm_1_op;
00277       Teuchos::RefCountPtr<RTOpPack::ReductTarget> norm_1_targ = norm_1_op.reduct_obj_create();
00278       const Thyra::VectorBase<Scalar>* vecs[] = { &*x0 };
00279       Teuchos::dyn_cast<const Thyra::SpmdVectorBase<Scalar> >(*x0).applyOp(
00280         &*Teuchos::DefaultComm<Index>::getComm()
00281         ,norm_1_op,1,vecs,0,static_cast<Thyra::VectorBase<Scalar>**>(NULL),&*norm_1_targ
00282         ,0,-1,0
00283         );
00284       nrm_x0 = norm_1_op(*norm_1_targ);
00285       *out << nrm_x0 << "\n";
00286       timer.stop();
00287       *out << "\n    time = " << timer.totalElapsedTime() << " seconds\n";
00288 
00289 #ifdef RTOPPACK_SPMD_APPLY_OP_DUMP
00290       RTOpPack::show_mpi_apply_op_dump = false;
00291 #endif
00292 #ifdef THYRA_SPMD_VECTOR_BASE_DUMP
00293       Thyra::SpmdVectorBase<Scalar>::show_dump = false;
00294 #endif
00295 
00296       const int N = 1;
00297       const int z_index = 1;
00298       Teuchos::Array<Teuchos::RefCountPtr<Thyra::ModelEvaluator<Scalar> > >
00299         models(N);
00300       Teuchos::Array<Scalar>
00301         weights(N);
00302       Teuchos::Array<Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > >
00303         z(N);
00304       for( int i = 0; i < N; ++i ) {
00305         models[i] = epetraThyraModel;
00306         weights[i] = 1.0;
00307         z[i] = epetraThyraModel->getNominalValues().get_p(z_index)->clone_v();
00308       }
00309       thyraModel =
00310         rcp(
00311           new Thyra::DefaultMultiPeriodModelEvaluator<Scalar>(
00312             1,&models[0],&weights[0],&z[0],z_index
00313             ,x_bar_space,f_bar_space
00314             )
00315           );
00316 
00317     }
00318 
00319 #endif // HAVE_MPI
00320 
00321     if(skipSolve)
00322       return ( success ? 0 : 1 );
00323 
00324     MoochoSolver::ESolutionStatus solution_status;
00325 
00326     //
00327     *out << "\n***\n*** Solving the initial forward problem\n***\n";
00328     //
00329 
00330     // Set the deliminator for the output files!
00331     solver.getSolver().set_output_context("fwd-init");
00332 
00333     // Set the solve mode to solve the forward problem
00334     solver.setDoSim(true);
00335     
00336     // Set the model
00337     solver.setModel(thyraModel);
00338     
00339     // Set the initial guess from files (if specified on commandline)
00340     solver.readInitialGuess(out.get());
00341     
00342     // Solve the initial forward problem
00343     solution_status = solver.solve();
00344 
00345     // Save the solution for model.x and model.p to be used later
00346     RefCountPtr<Thyra::VectorBase<Scalar> >
00347       x_opt = solver.getFinalPoint().get_x()->clone_v(),
00348       x_init = solver.getFinalPoint().get_x()->clone_v(),
00349       p_init = solver.getFinalPoint().get_p(0)->clone_v();
00350     
00351     //
00352     *out << "\n***\n*** Solving the perturbed forward problem\n***\n";
00353     //
00354 
00355     // Set the deliminator for the output files!
00356     solver.getSolver().set_output_context("fwd");
00357 
00358     // Set the solve mode to solve the forward problem
00359     solver.setDoSim(true);
00360    
00361     // Set the model
00362     solver.setModel(thyraModel);
00363 
00364     // Set the initial guess and the perturbed parameters
00365     if(1) {
00366       MEB::InArgs<Scalar> initialGuess = thyraModel->createInArgs();
00367       initialGuess.setArgs(thyraModel->getNominalValues());
00368       initialGuess.set_x(x_init);
00369       Thyra::Vt_S(&*p_init,perturbedParamScaling);
00370       initialGuess.set_p(0,p_init);
00371       //*out << "\nInitial Guess:\n" << Teuchos::describe(initialGuess,Teuchos::VERB_EXTREME);
00372       solver.setInitialGuess(initialGuess);
00373     }
00374 
00375     // Solve the perturbed forward problem
00376     solution_status = solver.solve();
00377 
00378     // Save the solution for model.x and model.p to be used later
00379     x_init = solver.getFinalPoint().get_x()->clone_v();
00380     p_init = solver.getFinalPoint().get_p(0)->clone_v();
00381     
00382     //
00383     *out << "\n***\n*** Solving the perturbed inverse problem\n***\n";
00384     //
00385 
00386     // Set the deliminator for the output files!
00387     solver.getSolver().set_output_context("inv");
00388 
00389     // Set the matching vector
00390     epetraModel->set_q(Thyra::get_Epetra_Vector(*epetraModel->get_x_map(),x_opt));
00391 
00392     // Set the solve mode to solve the inverse problem
00393     solver.setDoSim(false);
00394    
00395     // Set the model
00396     solver.setModel(thyraModel);
00397 
00398     // Set the initial guess for model.x and model.p
00399     if(1) {
00400       MEB::InArgs<Scalar> initialGuess = thyraModel->createInArgs();
00401       initialGuess.setArgs(thyraModel->getNominalValues());
00402       initialGuess.set_x(x_init);
00403       initialGuess.set_p(0,p_init);
00404       //*out << "\nInitial Guess:\n" << Teuchos::describe(initialGuess,Teuchos::VERB_EXTREME);
00405       solver.setInitialGuess(initialGuess);
00406     }
00407     
00408     // Solve the inverse problem
00409     solution_status = solver.solve();
00410     
00411     // Write the final solution
00412     solver.writeFinalSolution(out.get());
00413     
00414     // Write the final parameters to file
00415     lowsfCreator.writeParamsFile(*lowsFactory);
00416     
00417     //
00418     // Return the solution status (0 if successful)
00419     //
00420     
00421     return solution_status;
00422     
00423   }
00424   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
00425 
00426   return MoochoSolver::SOLVE_RETURN_EXCEPTION;
00427   
00428 }

Generated on Thu Sep 18 12:35:18 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1