• Main Page
  • Related Pages
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

Main_Solve.cpp

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 
00011 #include "Albany_Utils.hpp"
00012 #include "Albany_SolverFactory.hpp"
00013 
00014 #include "Piro_PerformSolve.hpp"
00015 #include "Teuchos_ParameterList.hpp"
00016 
00017 #include "Teuchos_GlobalMPISession.hpp"
00018 #include "Teuchos_TimeMonitor.hpp"
00019 #include "Teuchos_VerboseObject.hpp"
00020 #include "Teuchos_StandardCatchMacros.hpp"
00021 #include "Epetra_Map.h"  //Needed for serial, somehow
00022 
00023 //This header is for debug output -- writing of solution (xfinal) to MatrixMarket file
00024 #include "EpetraExt_MultiVectorOut.h"
00025 #include "EpetraExt_BlockMapOut.h"
00026 
00027 // Uncomment for run time nan checking
00028 // This is set in the toplevel CMakeLists.txt file
00029 
00030 #ifdef ENABLE_CHECK_FPE
00031 #include <math.h>
00032 //#include <Accelerate/Accelerate.h>
00033 #include <xmmintrin.h>
00034 #endif
00035 
00036 #include "Thyra_EpetraThyraWrappers.hpp"
00037 
00038 Teuchos::RCP<const Epetra_Vector>
00039 epetraVectorFromThyra(
00040   const Teuchos::RCP<const Epetra_Comm> &comm,
00041   const Teuchos::RCP<const Thyra::VectorBase<double> > &thyra)
00042 {
00043   Teuchos::RCP<const Epetra_Vector> result;
00044   if (Teuchos::nonnull(thyra)) {
00045     const Teuchos::RCP<const Epetra_Map> epetra_map = Thyra::get_Epetra_Map(*thyra->space(), comm);
00046     result = Thyra::get_Epetra_Vector(*epetra_map, thyra);
00047   }
00048   return result;
00049 }
00050 
00051 Teuchos::RCP<const Epetra_MultiVector>
00052 epetraMultiVectorFromThyra(
00053   const Teuchos::RCP<const Epetra_Comm> &comm,
00054   const Teuchos::RCP<const Thyra::MultiVectorBase<double> > &thyra)
00055 {
00056   Teuchos::RCP<const Epetra_MultiVector> result;
00057   if (Teuchos::nonnull(thyra)) {
00058     const Teuchos::RCP<const Epetra_Map> epetra_map = Thyra::get_Epetra_Map(*thyra->range(), comm);
00059     result = Thyra::get_Epetra_MultiVector(*epetra_map, thyra);
00060   }
00061   return result;
00062 }
00063 
00064 void epetraFromThyra(
00065   const Teuchos::RCP<const Epetra_Comm> &comm,
00066   const Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<double> > > &thyraResponses,
00067   const Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<double> > > > &thyraSensitivities,
00068   Teuchos::Array<Teuchos::RCP<const Epetra_Vector> > &responses,
00069   Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Epetra_MultiVector> > > &sensitivities)
00070 {
00071   responses.clear();
00072   responses.reserve(thyraResponses.size());
00073   typedef Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<double> > > ThyraResponseArray;
00074   for (ThyraResponseArray::const_iterator it_begin = thyraResponses.begin(),
00075       it_end = thyraResponses.end(),
00076       it = it_begin;
00077       it != it_end;
00078       ++it) {
00079     responses.push_back(epetraVectorFromThyra(comm, *it));
00080   }
00081 
00082   sensitivities.clear();
00083   sensitivities.reserve(thyraSensitivities.size());
00084   typedef Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<double> > > > ThyraSensitivityArray;
00085   for (ThyraSensitivityArray::const_iterator it_begin = thyraSensitivities.begin(),
00086       it_end = thyraSensitivities.end(),
00087       it = it_begin;
00088       it != it_end;
00089       ++it) {
00090     ThyraSensitivityArray::const_reference sens_thyra = *it;
00091     Teuchos::Array<Teuchos::RCP<const Epetra_MultiVector> > sens;
00092     sens.reserve(sens_thyra.size());
00093     for (ThyraSensitivityArray::value_type::const_iterator jt = sens_thyra.begin(),
00094         jt_end = sens_thyra.end();
00095         jt != jt_end;
00096         ++jt) {
00097         sens.push_back(epetraMultiVectorFromThyra(comm, *jt));
00098     }
00099     sensitivities.push_back(sens);
00100   }
00101 }
00102 
00103 int main(int argc, char *argv[]) {
00104 
00105   int status=0; // 0 = pass, failures are incremented
00106   bool success = true;
00107 
00108 #ifdef ALBANY_DEBUG
00109   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00110 #else // bypass printing process startup info
00111   Teuchos::GlobalMPISession mpiSession(&argc, &argv, NULL);
00112 #endif
00113 
00114 #ifdef ENABLE_CHECK_FPE
00115    // Catch FPEs
00116    _mm_setcsr(_MM_MASK_MASK &~
00117     (_MM_MASK_OVERFLOW | _MM_MASK_INVALID | _MM_MASK_DIV_ZERO) );
00118 #endif
00119 
00120   using Teuchos::RCP;
00121   using Teuchos::rcp;
00122 
00123   RCP<Teuchos::FancyOStream> out(Teuchos::VerboseObjectBase::getDefaultOStream());
00124 
00125   // Command-line argument for input file
00126   std::string xmlfilename;
00127   if(argc > 1){
00128 
00129     if(!strcmp(argv[1],"--help")){
00130       printf("albany [inputfile.xml]\n");
00131       exit(1);
00132     }
00133     else
00134       xmlfilename = argv[1];
00135 
00136   }
00137   else
00138     xmlfilename = "input.xml";
00139 
00140   try {
00141     RCP<Teuchos::Time> totalTime =
00142       Teuchos::TimeMonitor::getNewTimer("Albany: ***Total Time***");
00143 
00144     RCP<Teuchos::Time> setupTime =
00145       Teuchos::TimeMonitor::getNewTimer("Albany: Setup Time");
00146     Teuchos::TimeMonitor totalTimer(*totalTime); //start timer
00147     Teuchos::TimeMonitor setupTimer(*setupTime); //start timer
00148 
00149     Albany::SolverFactory slvrfctry(xmlfilename, Albany_MPI_COMM_WORLD);
00150     RCP<Epetra_Comm> appComm = Albany::createEpetraCommFromMpiComm(Albany_MPI_COMM_WORLD);
00151     RCP<Albany::Application> app;
00152     const RCP<Thyra::ModelEvaluator<double> > solver =
00153       slvrfctry.createThyraSolverAndGetAlbanyApp(app, appComm, appComm);
00154 
00155     setupTimer.~TimeMonitor();
00156 
00157     Teuchos::ParameterList &solveParams =
00158       slvrfctry.getAnalysisParameters().sublist("Solve", /*mustAlreadyExist =*/ false);
00159     // By default, request the sensitivities if not explicitly disabled
00160     solveParams.get("Compute Sensitivities", true);
00161 
00162     Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<double> > > thyraResponses;
00163     Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<double> > > > thyraSensitivities;
00164 
00165        // The PoissonSchrodinger_SchroPo and PoissonSchroMosCap1D tests seg fault as albanyApp is null -
00166        // For now, do not resize the response vectors. FIXME sort out this issue.
00167     if(Teuchos::nonnull(app))
00168       Piro::PerformSolveBase(*solver, solveParams, thyraResponses, thyraSensitivities, app->getAdaptSolMgr()->getSolObserver());
00169     else
00170       Piro::PerformSolveBase(*solver, solveParams, thyraResponses, thyraSensitivities);
00171 
00172     Teuchos::Array<Teuchos::RCP<const Epetra_Vector> > responses;
00173     Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Epetra_MultiVector> > > sensitivities;
00174     epetraFromThyra(appComm, thyraResponses, thyraSensitivities, responses, sensitivities);
00175 
00176     const int num_p = solver->Np(); // Number of *vectors* of parameters
00177     const int num_g = solver->Ng(); // Number of *vectors* of responses
00178 
00179     *out << "Finished eval of first model: Params, Responses "
00180       << std::setprecision(12) << std::endl;
00181 
00182     const Thyra::ModelEvaluatorBase::InArgs<double> nominal = solver->getNominalValues();
00183     for (int i=0; i<num_p; i++) {
00184       const Teuchos::RCP<const Epetra_Vector> p_init = epetraVectorFromThyra(appComm, nominal.get_p(i));
00185       p_init->Print(*out << "\nParameter vector " << i << ":\n");
00186     }
00187 
00188     for (int i=0; i<num_g-1; i++) {
00189       const RCP<const Epetra_Vector> g = responses[i];
00190       bool is_scalar = true;
00191 
00192       if (app != Teuchos::null)
00193         is_scalar = app->getResponse(i)->isScalarResponse();
00194 
00195       if (is_scalar) {
00196         g->Print(*out << "\nResponse vector " << i << ":\n");
00197 
00198         if (num_p == 0) {
00199           // Just calculate regression data
00200           status += slvrfctry.checkSolveTestResults(i, 0, g.get(), NULL);
00201         } else {
00202           for (int j=0; j<num_p; j++) {
00203             const RCP<const Epetra_MultiVector> dgdp = sensitivities[i][j];
00204             if (Teuchos::nonnull(dgdp)) {
00205               dgdp->Print(*out << "\nSensitivities (" << i << "," << j << "):!\n");
00206             }
00207             status += slvrfctry.checkSolveTestResults(i, j, g.get(), dgdp.get());
00208           }
00209         }
00210       }
00211     }
00212 
00213     const RCP<const Epetra_Vector> xfinal = responses.back();
00214     double mnv; xfinal->MeanValue(&mnv);
00215 
00216     // Create debug output object
00217     Teuchos::ParameterList &debugParams =
00218       slvrfctry.getParameters().sublist("Debug Output", true);
00219     bool writeToMatrixMarketSoln = debugParams.get("Write Solution to MatrixMarket", false);
00220     bool writeToCoutSoln = debugParams.get("Write Solution to Standard Output", false);
00221     if (writeToMatrixMarketSoln == true) { 
00222 
00223       //create serial map that puts the whole solution on processor 0
00224       int numMyElements = (xfinal->Comm().MyPID() == 0) ? app->getDiscretization()->getMap()->NumGlobalElements() : 0;
00225       const Epetra_Map serial_map(-1, numMyElements, 0, xfinal->Comm());
00226 
00227       //create importer from parallel map to serial map and populate serial solution xfinal_serial
00228       Epetra_Import importOperator(serial_map, *app->getDiscretization()->getMap());
00229       Epetra_Vector xfinal_serial(serial_map);
00230       xfinal_serial.Import(*app->getDiscretization()->getSolutionField(), importOperator, Insert);
00231 
00232       //writing to MatrixMarket file
00233       EpetraExt::MultiVectorToMatrixMarketFile("xfinal.mm", xfinal_serial);
00234     }
00235     if (writeToCoutSoln == true) 
00236        std::cout << "xfinal: " << *xfinal << std::endl;
00237 
00238     *out << "Main_Solve: MeanValue of final solution " << mnv << std::endl;
00239     *out << "\nNumber of Failed Comparisons: " << status << std::endl;
00240   }
00241   TEUCHOS_STANDARD_CATCH_STATEMENTS(true, std::cerr, success);
00242   if (!success) status+=10000;
00243 
00244   Teuchos::TimeMonitor::summarize(*out,false,true,false/*zero timers*/);
00245   return status;
00246 }

Generated on Wed Mar 26 2014 18:36:39 for Albany: a Trilinos-based PDE code by  doxygen 1.7.1