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 "Albany_NOXObserver.hpp"
00013 
00014 #include "Piro_Epetra_StokhosSolver.hpp"
00015 #include "Stokhos_EpetraVectorOrthogPoly.hpp"
00016 #include "Teuchos_VerboseObject.hpp"
00017 #include "Teuchos_StandardCatchMacros.hpp"
00018 #include "Stokhos.hpp"
00019 #include "Stokhos_Epetra.hpp"
00020 
00021 int main(int argc, char *argv[]) {
00022 
00023   int status=0; // 0 = pass, failures are incremented
00024   bool success = true;
00025   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00026   Teuchos::RCP<Teuchos::FancyOStream> out(Teuchos::VerboseObjectBase::getDefaultOStream());
00027 
00028   // Command-line argument for input file
00029   char * xmlfilename=0;
00030   char * sg_xmlfilename=0;
00031   char defaultfile[10]={"input.xml"};
00032   char sg_defaultfile[12]={"inputSG.xml"};
00033   bool do_initial_guess;
00034   if(argc>1){
00035     if(!strcmp(argv[1],"--help")){
00036       printf("albanySG [inputfile.xml inputfileSG.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   do_initial_guess = true;
00048       }
00049     }
00050   }
00051   else {
00052     xmlfilename=defaultfile;
00053     sg_xmlfilename=sg_defaultfile;
00054     do_initial_guess = true;
00055   }
00056        
00057   
00058   try {
00059 
00060     Teuchos::RCP<Teuchos::Time> totalTime =
00061       Teuchos::TimeMonitor::getNewTimer("AlbanySG: ***Total Time***");
00062     Teuchos::TimeMonitor totalTimer(*totalTime); //start timer
00063 
00064     // Setup communication objects
00065     Teuchos::RCP<Epetra_Comm> globalComm = 
00066       Albany::createEpetraCommFromMpiComm(Albany_MPI_COMM_WORLD);
00067 
00068     // Parse parameters
00069     Albany::SolverFactory sg_slvrfctry(sg_xmlfilename, Albany_MPI_COMM_WORLD);
00070     Teuchos::ParameterList& albanyParams = sg_slvrfctry.getParameters();
00071     Teuchos::RCP< Teuchos::ParameterList> piroParams = 
00072       Teuchos::rcp(&(albanyParams.sublist("Piro")),false);
00073     
00074     // Create stochastic Galerkin solver
00075     Teuchos::RCP<Piro::Epetra::StokhosSolver> sg_solver =
00076       Teuchos::rcp(new Piro::Epetra::StokhosSolver(piroParams, globalComm));
00077 
00078     // Get comm for spatial problem
00079     Teuchos::RCP<const Epetra_Comm> app_comm = sg_solver->getSpatialComm();
00080 
00081     // Compute initial guess if requested
00082     Teuchos::RCP<Epetra_Vector> ig;
00083     if (do_initial_guess) {
00084 
00085       // Create solver
00086       Albany::SolverFactory slvrfctry(
00087   xmlfilename,
00088   Albany::getMpiCommFromEpetraComm(*app_comm));
00089       Teuchos::RCP<EpetraExt::ModelEvaluator> solver = 
00090   slvrfctry.create(app_comm, app_comm);
00091 
00092       // Setup in/out args
00093       EpetraExt::ModelEvaluator::InArgs params_in = solver->createInArgs();
00094       EpetraExt::ModelEvaluator::OutArgs responses_out = 
00095   solver->createOutArgs();
00096       int np = params_in.Np();
00097       for (int i=0; i<np; i++) {
00098   Teuchos::RCP<const Epetra_Vector> p = solver->get_p_init(i);
00099   params_in.set_p(i, p);
00100       }
00101       int ng = responses_out.Ng();
00102       for (int i=0; i<ng; i++) {
00103   Teuchos::RCP<Epetra_Vector> g = 
00104     Teuchos::rcp(new Epetra_Vector(*(solver->get_g_map(i))));
00105   responses_out.set_g(i, g);
00106       }
00107 
00108       // Evaluate model
00109       solver->evalModel(params_in, responses_out);
00110 
00111       // Print responses (not last one since that is x)
00112       *out << std::endl;
00113       out->precision(8);
00114       for (int i=0; i<ng-1; i++) {
00115   if (responses_out.get_g(i) != Teuchos::null)
00116     *out << "Response " << i << " = " << std::endl 
00117          << *(responses_out.get_g(i)) << std::endl;
00118       }
00119 
00120       // Get final solution as initial guess
00121       ig = responses_out.get_g(ng-1);
00122 
00123       Teuchos::TimeMonitor::summarize(std::cout,false,true,false);
00124     }
00125 
00126     // Create SG solver
00127     Teuchos::RCP<Albany::Application> app;
00128     Teuchos::RCP<EpetraExt::ModelEvaluator> model = 
00129       sg_slvrfctry.createAlbanyAppAndModel(app, app_comm, ig);
00130 
00131     // Hack in rigid body modes for ML
00132     {
00133       Teuchos::RCP<Teuchos::ParameterList> sg_solver_params =
00134         Teuchos::sublist(Teuchos::sublist(piroParams, "Stochastic Galerkin"), "SG Solver Parameters");
00135       Teuchos::RCP<Teuchos::ParameterList> sg_prec_params =
00136         Teuchos::sublist(sg_solver_params, "SG Preconditioner");
00137 
00138       if (sg_prec_params->isParameter("Mean Preconditioner Type")) {
00139         if (sg_prec_params->get<std::string>("Mean Preconditioner Type") == "ML") {
00140 
00141           Teuchos::RCP<Teuchos::ParameterList> ml_params =
00142             Teuchos::sublist(sg_prec_params, "Mean Preconditioner Parameters");
00143           app->getProblem()->getNullSpace()->updateMLPL(ml_params);
00144           sg_solver->resetSolverParameters(*sg_solver_params);
00145         }
00146       }
00147     }
00148 
00149     // Setup SG solver
00150     {
00151       const Teuchos::RCP<NOX::Epetra::Observer > NOX_observer =
00152         Teuchos::rcp(new Albany_NOXObserver(app));
00153       sg_solver->setup(model, NOX_observer);
00154     }
00155 
00156     // Evaluate SG responses at SG parameters
00157     EpetraExt::ModelEvaluator::InArgs sg_inArgs = sg_solver->createInArgs();
00158     EpetraExt::ModelEvaluator::OutArgs sg_outArgs = 
00159       sg_solver->createOutArgs();
00160     int np = sg_inArgs.Np();
00161     for (int i=0; i<np; i++) {
00162       if (sg_inArgs.supports(EpetraExt::ModelEvaluator::IN_ARG_p_sg, i)) {
00163   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> p_sg = 
00164     sg_solver->get_p_sg_init(i);
00165   sg_inArgs.set_p_sg(i, p_sg);
00166       }
00167     }
00168 
00169     // By default, request the sensitivities if not explicitly disabled
00170     const bool computeSensitivities =
00171       sg_slvrfctry.getAnalysisParameters().sublist("Solve").get("Compute Sensitivities", true);
00172     int ng = sg_outArgs.Ng();
00173     for (int i=0; i<ng; i++) {
00174       if (sg_outArgs.supports(EpetraExt::ModelEvaluator::OUT_ARG_g_sg, i)) {
00175   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> g_sg =
00176     sg_solver->create_g_sg(i);
00177   sg_outArgs.set_g_sg(i, g_sg);
00178       }
00179 
00180       for (int j=0; j<np; j++) {
00181   EpetraExt::ModelEvaluator::DerivativeSupport ds =
00182     sg_outArgs.supports(EpetraExt::ModelEvaluator::OUT_ARG_DgDp_sg,i,j);
00183   if (computeSensitivities && 
00184       ds.supports(EpetraExt::ModelEvaluator::DERIV_MV_BY_COL)) {
00185     int ncol = sg_solver->get_p_map(j)->NumMyElements();
00186     Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdp_sg =
00187       sg_solver->create_g_mv_sg(i,ncol);
00188     sg_outArgs.set_DgDp_sg(i, j, dgdp_sg);
00189   }
00190       }
00191     }
00192 
00193     sg_solver->evalModel(sg_inArgs, sg_outArgs);
00194 
00195     bool printResponse = 
00196       albanyParams.sublist("Problem").get("Print Response Expansion", true);
00197     for (int i=0; i<ng-1; i++) {
00198       // Don't loop over last g which is x, since it is a long vector
00199       // to print out.
00200       if (sg_outArgs.supports(EpetraExt::ModelEvaluator::OUT_ARG_g_sg, i)) {
00201 
00202   // Print mean and standard deviation      
00203   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> g_sg = 
00204     sg_outArgs.get_g_sg(i);
00205   if (g_sg != Teuchos::null && app->getResponse(i)->isScalarResponse()) {
00206     Epetra_Vector g_mean(*(sg_solver->get_g_map(i)));
00207     Epetra_Vector g_std_dev(*(sg_solver->get_g_map(i)));
00208     g_sg->computeMean(g_mean);
00209     g_sg->computeStandardDeviation(g_std_dev);
00210     out->precision(12);
00211     out->setf(std::ios::scientific);
00212     *out << "Response " << i << " Mean =      " << std::endl 
00213          << g_mean << std::endl;
00214     *out << "Response " << i << " Std. Dev. = " << std::endl 
00215          << g_std_dev << std::endl;
00216     if (printResponse) {
00217       *out << "Response " << i << "           = " << std::endl 
00218      << *g_sg << std::endl;
00219       for (int j=0; j<np; j++) {
00220         EpetraExt::ModelEvaluator::DerivativeSupport ds =
00221     sg_outArgs.supports(EpetraExt::ModelEvaluator::OUT_ARG_DgDp_sg,i,j);
00222         if (!ds.none()) {
00223     Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdp_sg =
00224       sg_outArgs.get_DgDp_sg(i,j).getMultiVector();
00225     if (dgdp_sg != Teuchos::null)
00226       *out << "Response " << i << " Derivative " << j << " = " 
00227            << std::endl << *dgdp_sg << std::endl;
00228         }
00229       }
00230     }
00231 
00232     status += sg_slvrfctry.checkSGTestResults(i, g_sg, &g_mean, &g_std_dev);
00233   }
00234       }
00235     }
00236     *out << "\nNumber of Failed Comparisons: " << status << std::endl;
00237 
00238     totalTimer.~TimeMonitor();
00239     Teuchos::TimeMonitor::summarize(std::cout,false,true,false);
00240     Teuchos::TimeMonitor::zeroOutTimers();
00241 
00242   }
00243   TEUCHOS_STANDARD_CATCH_STATEMENTS(true, std::cerr, success);
00244   if (!success) status+=10000;
00245   
00246   return status;
00247 }
Main_SGSolve.cpp
