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 "Teuchos_GlobalMPISession.hpp"
00013 #include "Teuchos_TimeMonitor.hpp"
00014 #include "Teuchos_VerboseObject.hpp"
00015 #include "Teuchos_StandardCatchMacros.hpp"
00016 #include "Epetra_Map.h"  //Needed for serial, somehow
00017 
00018 int main(int argc, char *argv[]) {
00019 
00020   int status=0; // 0 = pass, failures are incremented
00021   bool success = true;
00022   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00023 
00024   using Teuchos::RCP;
00025   using Teuchos::rcp;
00026 
00027   RCP<Teuchos::FancyOStream> out(Teuchos::VerboseObjectBase::getDefaultOStream());
00028 
00029   // Command-line argument for input file
00030   char * xmlfilename=0;
00031   char defaultfile[10]={"input.xml"};
00032   if(argc>1){
00033     if(!strcmp(argv[1],"--help")){
00034       printf("albany [inputfile.xml]\n");
00035       exit(1);
00036     }
00037     else
00038       xmlfilename=argv[1];
00039   }
00040   else
00041     xmlfilename=defaultfile;
00042   
00043   char * xmladjfilename=0;
00044   char defaultadjfile[18]={"input_adjoint.xml"};
00045   if(argc>1){
00046     if(!strcmp(argv[1],"--help")){
00047       printf("albany [inputadjfile.xml]\n");
00048       exit(1);
00049     }
00050     else
00051       xmladjfilename=argv[2];
00052   }
00053   else
00054     xmladjfilename=defaultadjfile;
00055 
00056 
00057 
00058   try {
00059 
00060     RCP<Teuchos::Time> totalTime = 
00061       Teuchos::TimeMonitor::getNewTimer("Albany: ***Total Time***");
00062     RCP<Teuchos::Time> setupTime = 
00063       Teuchos::TimeMonitor::getNewTimer("Albany: Setup Time");
00064     Teuchos::TimeMonitor totalTimer(*totalTime); //start timer
00065     Teuchos::TimeMonitor setupTimer(*setupTime); //start timer
00066 
00067     Albany::SolverFactory slvrfctry(xmlfilename, Albany_MPI_COMM_WORLD);
00068     RCP<Epetra_Comm> appComm = Albany::createEpetraCommFromMpiComm(Albany_MPI_COMM_WORLD);
00069     RCP<EpetraExt::ModelEvaluator> App = slvrfctry.create(appComm, appComm);
00070 
00071     EpetraExt::ModelEvaluator::InArgs params_in = App->createInArgs();
00072     EpetraExt::ModelEvaluator::OutArgs responses_out = App->createOutArgs();
00073     int num_p = params_in.Np();     // Number of *vectors* of parameters
00074     int num_g = responses_out.Ng(); // Number of *vectors* of responses
00075     RCP<Epetra_Vector> p1;
00076     RCP<Epetra_Vector> g1;
00077 
00078     if (num_p > 0)
00079       p1 = rcp(new Epetra_Vector(*(App->get_p_init(0))));
00080     if (num_g > 1)
00081       g1 = rcp(new Epetra_Vector(*(App->get_g_map(0))));
00082     RCP<Epetra_Vector> xfinal =
00083       rcp(new Epetra_Vector(*(App->get_g_map(num_g-1)),true) );
00084 
00085     // Sensitivities
00086     RCP<Epetra_MultiVector> dgdp;
00087     if (num_p>0 && num_g>1) {
00088       // By default, request the sensitivities if not explicitly disabled
00089       const bool requestedSensitivities =
00090         slvrfctry.getAnalysisParameters().sublist("Solve").get("Compute Sensitivities", true);
00091 
00092       if (requestedSensitivities) {
00093         *out << "Main: DgDp sensititivies requested" << std::endl;
00094         *out << " Num Responses: " << g1->GlobalLength()
00095           << ",   Num Parameters: " << p1->GlobalLength() << std::endl;
00096 
00097         const bool supportsSensitivities =
00098           !responses_out.supports(EpetraExt::ModelEvaluator::OUT_ARG_DgDp, 0, 0).none();
00099 
00100         if (supportsSensitivities && p1->GlobalLength() > 0) {
00101           *out << "Main: Model supports requested DgDp sensitivities" << std::endl;
00102           dgdp = rcp(new Epetra_MultiVector(g1->Map(), p1->GlobalLength()));
00103         } else {
00104           *out << "Main: Model does not supports requested DgDp sensitivities" << std::endl;
00105         }
00106       }
00107     }
00108 
00109     // Evaluate first model
00110     if (num_p > 0)  params_in.set_p(0,p1);
00111     if (num_g > 1)  responses_out.set_g(0,g1);
00112     responses_out.set_g(num_g-1,xfinal);
00113 
00114     if (Teuchos::nonnull(dgdp)) responses_out.set_DgDp(0,0,dgdp);
00115     setupTimer.~TimeMonitor();
00116     App->evalModel(params_in, responses_out);
00117 
00118     *out << "Finished eval of first model: Params, Responses " 
00119          << std::setprecision(12) << std::endl;
00120     if (num_p>0) p1->Print(*out << "\nParameters!\n");
00121     if (num_g>1) g1->Print(*out << "\nResponses!\n");
00122     if (Teuchos::nonnull(dgdp))
00123       dgdp->Print(*out << "\nSensitivities!\n");
00124     double mnv; xfinal->MeanValue(&mnv);
00125     *out << "Main_Solve: MeanValue of final solution " << mnv << std::endl;
00126 
00127     //cout << "Final Solution \n" << *xfinal << std::endl;
00128 
00129     status += slvrfctry.checkSolveTestResults(0, 0, g1.get(), dgdp.get());
00130     *out << "\nNumber of Failed Comparisons: " << status << std::endl;
00131 
00132     // write out the current mesh to an exodus file
00133 
00134     // promote the mesh to higher order
00135 
00136 
00137     // Start adjoint solve
00138 
00139     RCP<Epetra_Vector> xinit =
00140       rcp(new Epetra_Vector(*(App->get_g_map(num_g-1)),true) );
00141     
00142     RCP<Teuchos::Time> totalAdjTime = 
00143       Teuchos::TimeMonitor::getNewTimer("Albany: ***Total Time***");
00144     RCP<Teuchos::Time> setupAdjTime = 
00145       Teuchos::TimeMonitor::getNewTimer("Albany: Setup Time");
00146     Teuchos::TimeMonitor totalAdjTimer(*totalAdjTime); //start timer
00147     Teuchos::TimeMonitor setupAdjTimer(*setupAdjTime); //start timer
00148 
00149     Albany::SolverFactory adjslvrfctry(xmladjfilename, Albany_MPI_COMM_WORLD);
00150     RCP<Epetra_Comm> adjappComm = Albany::createEpetraCommFromMpiComm(Albany_MPI_COMM_WORLD);
00151     RCP<EpetraExt::ModelEvaluator> AdjApp = adjslvrfctry.create(adjappComm, adjappComm, xinit);
00152 
00153     EpetraExt::ModelEvaluator::InArgs adj_params_in = AdjApp->createInArgs();
00154     EpetraExt::ModelEvaluator::OutArgs adj_responses_out = AdjApp->createOutArgs();
00155     int adj_num_p = adj_params_in.Np();     // Number of *vectors* of parameters
00156     int adj_num_g = adj_responses_out.Ng(); // Number of *vectors* of responses
00157     RCP<Epetra_Vector> adj_p1;
00158     RCP<Epetra_Vector> adj_g1;
00159 
00160     if (adj_num_p > 0)
00161       adj_p1 = rcp(new Epetra_Vector(*(AdjApp->get_p_init(0))));
00162     if (adj_num_g > 1)
00163       adj_g1 = rcp(new Epetra_Vector(*(AdjApp->get_g_map(0))));
00164     RCP<Epetra_Vector> adj_xfinal =
00165       rcp(new Epetra_Vector(*(AdjApp->get_g_map(adj_num_g-1)),true) );
00166 
00167     // Adjoint sensitivities
00168     RCP<Epetra_MultiVector> adj_dgdp;
00169     if (adj_num_p>0 && adj_num_g>1) {
00170       // By default, request the sensitivities if not explicitly disabled
00171       const bool requestedSensitivities =
00172         adjslvrfctry.getAnalysisParameters().sublist("Solve").get("Compute Sensitivities", true);
00173 
00174       if (requestedSensitivities) {
00175         *out << "Main: Adjoint DgDp sensititivies requested" << std::endl;
00176         *out << " Num Responses: " << g1->GlobalLength()
00177           << ",   Num Parameters: " << p1->GlobalLength() << std::endl;
00178 
00179         const bool adj_supportsSensitivities =
00180           !adj_responses_out.supports(EpetraExt::ModelEvaluator::OUT_ARG_DgDp, 0, 0).none();
00181 
00182         if (adj_supportsSensitivities && adj_p1->GlobalLength() > 0) {
00183           *out << "Main: Model supports requested adjoint DgDp sensitivities" << std::endl;
00184           adj_dgdp = rcp(new Epetra_MultiVector(adj_g1->Map(), adj_p1->GlobalLength() ));
00185         } else {
00186           *out << "Main: Model does not supports requested adjoint DgDp sensitivities" << std::endl;
00187         }
00188       }
00189     }
00190 
00191 
00192     /* Set the initial guess for the adjoint to be the forward solution
00193        projected onto the adjoint degrees of freedom.  The adjoint problem
00194        is linear, so this is only performed so that the derivative of the
00195        response function (which forms the RHS) is correct.
00196     */
00197 
00198     // Teuchos::RCP<ENAT::SGNOXSolver> AdjApp =
00199     //  Teuchos::rcp_dynamic_cast<ENAT::SGNOXSolver>(adjslvrfctry.create(adjapp_comm, adjapp_comm, xfinal));
00200 
00201 
00202 
00203     // Evaluate adjoint model
00204     if (adj_num_p > 0)  adj_params_in.set_p(0,adj_p1);
00205     if (adj_num_g > 1)  adj_responses_out.set_g(0,adj_g1);
00206     adj_responses_out.set_g(adj_num_g-1,adj_xfinal);
00207 
00208     if (Teuchos::nonnull(adj_dgdp)) adj_responses_out.set_DgDp(0,0,adj_dgdp);
00209     setupAdjTimer.~TimeMonitor();
00210     AdjApp->evalModel(adj_params_in, adj_responses_out);
00211 
00212     *out << "Finished eval of adjoint model: Params, Responses " 
00213          << std::setprecision(12) << std::endl;
00214     if (adj_num_p>0) adj_p1->Print(*out << "\nParameters!\n");
00215     if (adj_num_g>1) adj_g1->Print(*out << "\nResponses!\n");
00216     if (Teuchos::nonnull(adj_dgdp))
00217       adj_dgdp->Print(*out << "\nSensitivities!\n");
00218     double adj_mnv; adj_xfinal->MeanValue(&adj_mnv);
00219     *out << "Main_Solve: MeanValue of final solution " << adj_mnv << std::endl;
00220 
00221     //cout << "Final Solution \n" << *xfinal << std::endl;
00222 
00223     status += adjslvrfctry.checkSolveTestResults(0, 0, adj_g1.get(), adj_dgdp.get());
00224     *out << "\nNumber of Failed Comparisons: " << status << std::endl;
00225   }
00226 
00227   TEUCHOS_STANDARD_CATCH_STATEMENTS(true, std::cerr, success);
00228   if (!success) status+=10000;
00229 
00230   Teuchos::TimeMonitor::summarize(*out,false,true,false/*zero timers*/);
00231   return status;
00232 }
Main_Adjoint.cpp
