Albany: a Trilinos-based PDE code: Main

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 
00007 #include <iostream>
00008 #include <string>
00009 
00010 #include "Albany_Utils.hpp"
00011 #include "Albany_SolverFactory.hpp"
00012 #include "Piro_Epetra_StokhosSolver.hpp"
00013 #include "Stokhos_EpetraVectorOrthogPoly.hpp"
00014 #include "Teuchos_VerboseObject.hpp"
00015 #include "Teuchos_StandardCatchMacros.hpp"
00016 #include "Stokhos.hpp"
00017 #include "Stokhos_Epetra.hpp"
00018 
00019 int main(int argc, char *argv[]) {
00020 
00021   int status=0; // 0 = pass, failures are incremented
00022   bool success = true;
00023   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00024   Teuchos::RCP<Teuchos::FancyOStream> out(Teuchos::VerboseObjectBase::getDefaultOStream());
00025 
00026   // Command-line argument for input file
00027   char * xmlfilename=0;
00028   char * sg_xmlfilename=0;
00029   char * adjsg_xmlfilename=0;
00030   char defaultfile[10]={"input.xml"};
00031   char sg_defaultfile[12]={"inputSG.xml"};
00032   char adjsg_defaultfile[20]={"inputSG_adjoint.xml"};
00033   bool do_initial_guess;
00034   if(argc>1){
00035     if(!strcmp(argv[1],"--help")){
00036       printf("albanySG [inputfile.xml inputfileSG.xml inputfileSGadjoint.xml]\n");
00037       exit(1);
00038     }
00039     else {
00040       if (argc == 2) {
00041   sg_xmlfilename = argv[1];
00042   do_initial_guess = false;
00043       }
00044       else {
00045   xmlfilename=argv[1];
00046   sg_xmlfilename = argv[2];
00047   adjsg_xmlfilename = argv[3];
00048   do_initial_guess = true;
00049       }
00050     }
00051   }
00052   else {
00053     xmlfilename=defaultfile;
00054     sg_xmlfilename=sg_defaultfile;
00055     adjsg_xmlfilename=adjsg_defaultfile;
00056     do_initial_guess = true;
00057   }
00058        
00059   
00060   try {
00061 
00062     Teuchos::RCP<Teuchos::Time> totalTime =
00063       Teuchos::TimeMonitor::getNewTimer("AlbanySGAdjoint: ***Total Time***");
00064     Teuchos::TimeMonitor totalTimer(*totalTime); //start timer
00065     
00066     // Setup communication objects
00067     Teuchos::RCP<Epetra_Comm> globalComm = 
00068       Albany::createEpetraCommFromMpiComm(Albany_MPI_COMM_WORLD);
00069 
00070     Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_forward_solution;
00071 
00072     //
00073     // Solve forward problem
00074     //
00075     {
00076       
00077     Teuchos::RCP<Teuchos::Time> forwardTime =
00078       Teuchos::TimeMonitor::getNewTimer("AlbanySGAdjoint: ***Forward Solver Time***");
00079     Teuchos::TimeMonitor forwardTimer(*forwardTime); //start timer
00080 
00081     // Parse parameters
00082     Albany::SolverFactory sg_slvrfctry(sg_xmlfilename, Albany_MPI_COMM_WORLD);
00083     Teuchos::ParameterList& albanyParams = sg_slvrfctry.getParameters();
00084     Teuchos::RCP< Teuchos::ParameterList> piroParams = 
00085       Teuchos::rcp(&(albanyParams.sublist("Piro")),false);
00086     
00087     // Create stochastic Galerkin solver
00088     Teuchos::RCP<Piro::Epetra::StokhosSolver> sg_solver =
00089       Teuchos::rcp(new Piro::Epetra::StokhosSolver(piroParams, globalComm));
00090 
00091     // Get comm for spatial problem
00092     Teuchos::RCP<const Epetra_Comm> app_comm = sg_solver->getSpatialComm();
00093 
00094     // Compute initial guess if requested
00095     Teuchos::RCP<Epetra_Vector> ig;
00096     if (do_initial_guess) {
00097 
00098       // Create solver
00099       Albany::SolverFactory slvrfctry(
00100   xmlfilename,
00101   Albany::getMpiCommFromEpetraComm(*app_comm));
00102       Teuchos::RCP<EpetraExt::ModelEvaluator> solver = 
00103   slvrfctry.create(app_comm, app_comm);
00104 
00105       // Setup in/out args
00106       EpetraExt::ModelEvaluator::InArgs params_in = solver->createInArgs();
00107       EpetraExt::ModelEvaluator::OutArgs responses_out = 
00108   solver->createOutArgs();
00109       int np = params_in.Np();
00110       for (int i=0; i<np; i++) {
00111   Teuchos::RCP<const Epetra_Vector> p = solver->get_p_init(i);
00112   params_in.set_p(i, p);
00113       }
00114       int ng = responses_out.Ng();
00115       for (int i=0; i<ng; i++) {
00116   Teuchos::RCP<Epetra_Vector> g = 
00117     Teuchos::rcp(new Epetra_Vector(*(solver->get_g_map(i))));
00118   responses_out.set_g(i, g);
00119       }
00120 
00121       // Evaluate model
00122       solver->evalModel(params_in, responses_out);
00123 
00124       // Print responses (not last one since that is x)
00125       *out << std::endl;
00126       out->precision(8);
00127       for (int i=0; i<ng-1; i++) {
00128   if (responses_out.get_g(i) != Teuchos::null)
00129     *out << "Response " << i << " = " << std::endl 
00130          << *(responses_out.get_g(i)) << std::endl;
00131       }
00132 
00133     }
00134 
00135     // Create SG solver
00136     Teuchos::RCP<Albany::Application> app;
00137     Teuchos::RCP<EpetraExt::ModelEvaluator> model = 
00138       sg_slvrfctry.createAlbanyAppAndModel(app, app_comm, ig);
00139     sg_solver->setup(model);
00140 
00141     // Evaluate SG responses at SG parameters
00142     EpetraExt::ModelEvaluator::InArgs sg_inArgs = sg_solver->createInArgs();
00143     EpetraExt::ModelEvaluator::OutArgs sg_outArgs = 
00144       sg_solver->createOutArgs();
00145     int np = sg_inArgs.Np();
00146     for (int i=0; i<np; i++) {
00147       if (sg_inArgs.supports(EpetraExt::ModelEvaluator::IN_ARG_p_sg, i)) {
00148   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> p_sg = 
00149     sg_solver->get_p_sg_init(i);
00150   sg_inArgs.set_p_sg(i, p_sg);
00151       }
00152     }
00153 
00154     int ng = sg_outArgs.Ng();
00155     for (int i=0; i<ng; i++) {
00156       if (sg_outArgs.supports(EpetraExt::ModelEvaluator::OUT_ARG_g_sg, i)) {
00157   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> g_sg =
00158     sg_solver->create_g_sg(i);
00159   sg_outArgs.set_g_sg(i, g_sg);
00160       }
00161     }
00162 
00163     sg_solver->evalModel(sg_inArgs, sg_outArgs);
00164 
00165     for (int i=0; i<ng-1; i++) {
00166       // Don't loop over last g which is x, since it is a long vector
00167       // to print out.
00168       if (sg_outArgs.supports(EpetraExt::ModelEvaluator::OUT_ARG_g_sg, i)) {
00169 
00170   // Print mean and standard deviation      
00171   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> g_sg = 
00172     sg_outArgs.get_g_sg(i);
00173   Epetra_Vector g_mean(*(sg_solver->get_g_map(i)));
00174   Epetra_Vector g_std_dev(*(sg_solver->get_g_map(i)));
00175   g_sg->computeMean(g_mean);
00176   g_sg->computeStandardDeviation(g_std_dev);
00177   out->precision(12);
00178   *out << "Response " << i << " Mean =      " << std::endl 
00179        << g_mean << std::endl;
00180   *out << "Response " << i << " Std. Dev. = " << std::endl 
00181        << g_std_dev << std::endl;
00182 
00183   status += sg_slvrfctry.checkSGTestResults(0, g_sg);
00184       }
00185     }
00186     *out << "\nNumber of Failed Comparisons: " << status << std::endl;
00187   
00188     sg_forward_solution = sg_outArgs.get_g_sg(ng-1);
00189 
00190     }
00191 
00192 
00193     /* Space reserved for the projection of the forward solution onto
00194        the higher order basis for the adjoint solution.  
00195        In general, this will require projecting is physical space
00196        as well as in stochastic space.
00197     */
00198 
00199 
00200 
00201     // 
00202     // Solve adjoint problem
00203     //
00204     {
00205 
00206     Teuchos::RCP<Teuchos::Time> adjointTime =
00207       Teuchos::TimeMonitor::getNewTimer("AlbanySG: ***Adjoint Solver Time***");
00208     Teuchos::TimeMonitor adjtotalTimer(*adjointTime); //start timer
00209 
00210     // Parse parameters
00211     Albany::SolverFactory sg_slvrfctry(adjsg_xmlfilename, 
00212                Albany_MPI_COMM_WORLD);
00213     Teuchos::ParameterList& albanyParams = sg_slvrfctry.getParameters();
00214     Teuchos::RCP< Teuchos::ParameterList> piroParams = 
00215       Teuchos::rcp(&(albanyParams.sublist("Piro")),false);
00216     
00217     // Create stochastic Galerkin solver
00218     Teuchos::RCP<Piro::Epetra::StokhosSolver> sg_solver =
00219       Teuchos::rcp(new Piro::Epetra::StokhosSolver(piroParams, globalComm));
00220 
00221     // Get comm for spatial problem
00222     Teuchos::RCP<const Epetra_Comm> app_comm = sg_solver->getSpatialComm();
00223 
00224     // Create SG solver
00225     Teuchos::RCP<Albany::Application> app;
00226     Teuchos::RCP<EpetraExt::ModelEvaluator> model = 
00227       sg_slvrfctry.createAlbanyAppAndModel(app, app_comm);
00228     sg_solver->setup(model);
00229 
00230     // Set projected forward solution as the initial guess
00231     sg_solver->set_x_sg_init(*sg_forward_solution);
00232 
00233     // Evaluate SG responses at SG parameters
00234     EpetraExt::ModelEvaluator::InArgs sg_inArgs = sg_solver->createInArgs();
00235     EpetraExt::ModelEvaluator::OutArgs sg_outArgs = 
00236       sg_solver->createOutArgs();
00237     int np = sg_inArgs.Np();
00238     for (int i=0; i<np; i++) {
00239       if (sg_inArgs.supports(EpetraExt::ModelEvaluator::IN_ARG_p_sg, i)) {
00240   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> p_sg = 
00241     sg_solver->get_p_sg_init(i);
00242   sg_inArgs.set_p_sg(i, p_sg);
00243       }
00244     }
00245 
00246     int ng = sg_outArgs.Ng();
00247     for (int i=0; i<ng; i++) {
00248       if (sg_outArgs.supports(EpetraExt::ModelEvaluator::OUT_ARG_g_sg, i)) {
00249   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> g_sg =
00250     sg_solver->create_g_sg(i);
00251   sg_outArgs.set_g_sg(i, g_sg);
00252       }
00253     }
00254 
00255     sg_solver->evalModel(sg_inArgs, sg_outArgs);
00256 
00257     for (int i=0; i<ng-1; i++) {
00258       // Don't loop over last g which is x, since it is a long vector
00259       // to print out.
00260       if (sg_outArgs.supports(EpetraExt::ModelEvaluator::OUT_ARG_g_sg, i)) {
00261 
00262   // Print mean and standard deviation      
00263   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> g_sg = 
00264     sg_outArgs.get_g_sg(i);
00265   Epetra_Vector g_mean(*(sg_solver->get_g_map(i)));
00266   Epetra_Vector g_std_dev(*(sg_solver->get_g_map(i)));
00267   g_sg->computeMean(g_mean);
00268   g_sg->computeStandardDeviation(g_std_dev);
00269   out->precision(12);
00270   *out << "Response " << i << " Mean =      " << std::endl 
00271        << g_mean << std::endl;
00272   *out << "Response " << i << " Std. Dev. = " << std::endl 
00273        << g_std_dev << std::endl;
00274 
00275   status += sg_slvrfctry.checkSGTestResults(0, g_sg);
00276       }
00277     }
00278     *out << "\nNumber of Failed Comparisons: " << status << std::endl;
00279 
00280     /* Space reserved for computing the error representation which involves
00281        integrating over both physical and stochastic space and may require
00282        a number of projections.
00283     */
00284 
00285     }
00286 
00287     totalTimer.~TimeMonitor();
00288     Teuchos::TimeMonitor::summarize(std::cout,false,true,false);
00289 
00290   }
00291   TEUCHOS_STANDARD_CATCH_STATEMENTS(true, std::cerr, success);
00292   if (!success) status+=10000;
00293 
00294   return status;
00295 }
Main_SGAdjoint.cpp
