poisson_main.cpp

00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 //
00010 // This is a simple program to exercise the FEI, for the purposes of
00011 // testing, code tuning and scaling studies.
00012 //
00013 // This program assembles a linear system from a 2D square Poisson
00014 // problem, using 4-node square elements. There is only 1 degree-of-
00015 // freedom per node.
00016 //
00017 // This problem was coded up with the help of Ray Tuminaro.
00018 // Alan Williams 12-20-2000
00019 //
00020 // The input file for this program should provide the following:
00021 //   WHICH_FEI <OLDFEI|fei::FEI_Impl>
00022 //   SOLVER_LIBRARY <library-name> -- e.g., Aztec
00023 //   L <int> -- the global length (num-elements) on a side of the 2D square 
00024 //
00025 //
00026 
00027 #include <fei_iostream.hpp>
00028 #include <cmath>
00029 #include <fei_base.hpp>
00030 #include <FEI_Implementation.hpp>
00031 
00032 //Now make provision for using any one of several solver libraries. This is
00033 //handled by the code in LibraryFactory.{hpp,cpp}.
00034 
00035 #include <test_utils/LibraryFactory.hpp>
00036 #include <fei_LibraryWrapper.hpp>
00037 
00038 //And, we need to include some headers for utility classes which are simply
00039 //used for setting up the data for the example problem.
00040 
00041 #include <test_utils/Poisson_Elem.hpp>
00042 #include <test_utils/PoissonData.hpp>
00043 
00044 #include <test_utils/ElemBlock.hpp>
00045 #include <test_utils/CRSet.hpp>
00046 #include <test_utils/CommNodeSet.hpp>
00047 #include <test_utils/DataReader.hpp>
00048 
00049 #include <test_utils/SolnCheck.hpp>
00050 
00051 #include <test_utils/fei_test_utils.hpp>
00052 #include <snl_fei_Utils.hpp>
00053 //
00054 //Include definitions of macros to call functions and check the return code.
00055 //
00056 #undef fei_file
00057 #define fei_file "poisson_main.cpp"
00058 #include <fei_ErrMacros.hpp>
00059 
00060 //============================================================================
00061 //Here's the main...
00062 //============================================================================
00063 int poisson_main(int argc, char** argv,
00064                  MPI_Comm comm, int numProcs, int localProc){
00065   int outputLevel;
00066   int masterProc = 0, err = 0;
00067 
00068   double start_time = fei::utils::cpu_time();
00069 
00070   std::vector<std::string> stdstrings;
00071   CHK_ERR( fei_test_utils::get_filename_and_read_input(argc, argv,
00072             comm, localProc,
00073             stdstrings) );
00074   const char** params = NULL;
00075   int numParams = 0;
00076   fei::utils::strings_to_char_ptrs(stdstrings, numParams, params);
00077   
00078   fei::ParameterSet paramset;
00079   fei::utils::parse_strings(stdstrings, " ", paramset);
00080 
00081   std::string which_fei;
00082   std::string solverName;
00083   int L = 0;
00084 
00085   int errcode = 0;
00086   errcode += paramset.getStringParamValue("SOLVER_LIBRARY", solverName);
00087   errcode += paramset.getStringParamValue("WHICH_FEI", which_fei);
00088   errcode += paramset.getIntParamValue("L", L);
00089 
00090   if (errcode != 0) {
00091     FEI_CERR << "Failed to find one or more required parameters in input-file."
00092        << FEI_ENDL << "Required parameters:"<<FEI_ENDL
00093        << "SOLVER_LIBRARY" << FEI_ENDL
00094        << "WHICH_FEI" << FEI_ENDL
00095        << "L" << FEI_ENDL;
00096     return(-1);
00097   }
00098 
00099   if (localProc == 0) {
00100     int nodes = (L+1)*(L+1);
00101     int eqns = nodes;
00102     FEI_COUT << FEI_ENDL;
00103     FEI_COUT << "========================================================" 
00104    << FEI_ENDL;
00105     FEI_COUT << "Square size     L: " << L << " elements." << FEI_ENDL;
00106     FEI_COUT << "Global number of elements: " << L*L << FEI_ENDL;
00107     FEI_COUT << "Global number of nodes: " << nodes << FEI_ENDL;
00108     FEI_COUT << "Global number of equations: " << eqns <<FEI_ENDL;
00109     FEI_COUT << "========================================================" 
00110    << FEI_ENDL;
00111   }
00112 
00113   outputLevel = fei_test_utils::whichArg(numParams, params, "outputLevel 1");
00114   if (outputLevel >= 0) outputLevel = 1;
00115   if (outputLevel < 0) outputLevel = 0;
00116 
00117   if ((masterProc == localProc)&&(outputLevel>0)) {
00118     fei_test_utils::print_args(argc, argv);
00119   }
00120 
00121   if (outputLevel == 1) {
00122     if (localProc != 0) outputLevel = 0;
00123   }
00124 
00125   //PoissonData is the object that will be in charge of generating the
00126   //data to pump into the FEI object.
00127 
00128   PoissonData poissonData(L, numProcs, localProc, outputLevel);
00129 
00130   double start_init_time = fei::utils::cpu_time();
00131 
00132   fei::SharedPtr<fei::Factory> factory;
00133   fei::SharedPtr<LibraryWrapper> wrapper;
00134   fei::SharedPtr<FEI> fei;
00135 
00136   if (which_fei == "OLDFEI") {
00137     try {
00138       wrapper = fei::create_LibraryWrapper(comm, solverName.c_str());
00139     }
00140     catch (std::runtime_error& exc) {
00141       FEI_CERR << exc.what()<<FEI_ENDL;
00142       ERReturn(-1);
00143     }
00144     fei.reset(new FEI_Implementation(wrapper, comm));
00145   }
00146   else if (which_fei == "fei::FEI_Impl") {
00147     try {
00148       factory = fei::create_fei_Factory(comm, solverName.c_str());
00149     }
00150     catch (std::runtime_error& exc) {
00151       FEI_CERR << exc.what()<<FEI_ENDL;
00152       ERReturn(-1);
00153     }
00154     fei = factory->createFEI(comm);
00155   }
00156   else {
00157     FEI_CERR << "poisson_main ERROR, value of 'WHICH_FEI' must be 'OLDFEI' or 'fei::FEI_Impl'"<< FEI_ENDL;
00158     ERReturn(-1);
00159   }
00160 
00161   const char* feiVersionString;
00162   CHK_ERR( fei->version(feiVersionString) );
00163 
00164   if (localProc==0) FEI_COUT << feiVersionString << FEI_ENDL;
00165 
00166   //load some parameters.
00167   CHK_ERR( fei->parameters( numParams, params ) );
00168 
00169   delete [] params;
00170 
00171   if (outputLevel>0 && localProc==0) FEI_COUT << "setSolveType" << FEI_ENDL;
00172   CHK_ERR( fei->setSolveType(FEI_SINGLE_SYSTEM) );
00173 
00174 
00175   int numFields = poissonData.getNumFields();
00176   int* fieldSizes = poissonData.getFieldSizes();
00177   int* fieldIDs = poissonData.getFieldIDs();
00178 
00179   if (outputLevel>0 && localProc==0) FEI_COUT << "initFields" << FEI_ENDL;
00180   CHK_ERR( fei->initFields( numFields, fieldSizes, fieldIDs ) );
00181 
00182 
00183   CHK_ERR( init_elem_connectivities(fei.get(), poissonData) );
00184 
00185   CHK_ERR( set_shared_nodes(fei.get(), poissonData) );
00186 
00187   //The FEI_COUT and IOS_... macros are defined in base/fei_iostream.hpp
00188   FEI_COUT.setf(IOS_FIXED, IOS_FLOATFIELD);
00189 
00190   if (outputLevel>0 && localProc==0) FEI_COUT << "initComplete" << FEI_ENDL;
00191   CHK_ERR( fei->initComplete() );
00192 
00193   double fei_init_time = fei::utils::cpu_time() - start_init_time;
00194 
00195   //Now the initialization phase is complete. Next we'll do the load phase,
00196   //which for this problem just consists of loading the element data
00197   //(element-wise stiffness arrays and load vectors) and the boundary
00198   //condition data.
00199   //This simple problem doesn't have an constraint relations, etc.
00200 
00201   double start_load_time = fei::utils::cpu_time();
00202 
00203   CHK_ERR( load_elem_data(fei.get(), poissonData) );
00204 
00205   CHK_ERR( load_BC_data(fei.get(), poissonData) );
00206 
00207   double fei_load_time = fei::utils::cpu_time() - start_load_time;
00208 
00209   //
00210   //now the load phase is complete, so we're ready to launch the underlying
00211   //solver and solve Ax=b
00212   //
00213   int status;
00214   if (outputLevel>0 && localProc==0) FEI_COUT << "solve..." << FEI_ENDL;
00215   double start_solve_time = fei::utils::cpu_time();
00216   err = fei->solve(status);
00217   if (err) {
00218     if (localProc==0) FEI_COUT << "solve returned err: " << err << FEI_ENDL;
00219   }
00220 
00221   double iTime, lTime, sTime, rTime;
00222   CHK_ERR(fei->cumulative_cpu_times(iTime, lTime, sTime, rTime) );
00223 
00224   double solve_time = fei::utils::cpu_time() - start_solve_time;
00225 
00226   if (localProc == 0) {
00227     FEI_COUT << "FEI cpu-times:" << FEI_ENDL
00228    << "    init. phase: " << iTime << FEI_ENDL
00229    << "    load  phase: " << lTime << FEI_ENDL
00230    << "    solve  time: " << sTime << FEI_ENDL;
00231   }
00232 
00233   double norm = 0.0;
00234   FEI_COUT.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
00235   CHK_ERR( fei->residualNorm(1, 1, &fieldIDs[0], &norm) );
00236   if (localProc == 0) {
00237     FEI_COUT << "returned residual norm: " << norm << FEI_ENDL;
00238   }
00239 
00240   int itersTaken = 0;
00241   CHK_ERR( fei->iterations(itersTaken) );
00242 
00243   //
00244   //We oughtta make sure the solution we just computed is correct...
00245   //
00246 
00247   int numNodes = 0;
00248   CHK_ERR( fei->getNumLocalNodes(numNodes) );
00249 
00250   double maxErr = 0.0;
00251   if (numNodes > 0) {
00252     int lenNodeIDs = numNodes;
00253     GlobalID* nodeIDs = new GlobalID[lenNodeIDs];
00254     double* soln = new double[lenNodeIDs];
00255     if (nodeIDs != NULL && soln != NULL) {
00256       CHK_ERR( fei->getLocalNodeIDList(numNodes, nodeIDs, lenNodeIDs) );
00257 
00258       int fieldID = 1;
00259       CHK_ERR( fei->getNodalFieldSolution(fieldID, numNodes, nodeIDs, soln));
00260 
00261       for(int i=0; i<numNodes; i++) {
00262   int nID = (int)nodeIDs[i];
00263   double x = (1.0* ((nID-1)%(L+1)))/L;
00264   double y = (1.0* ((nID-1)/(L+1)))/L;
00265 
00266   double exactSoln = x*x + y*y;
00267   double error = std::abs(exactSoln - soln[i]);
00268   if (maxErr < error) maxErr = error;
00269       }
00270 
00271       delete [] nodeIDs;
00272       delete [] soln;
00273     }
00274     else {
00275       FEI_CERR << "allocation of nodeIDs or soln failed." << FEI_ENDL; 
00276     }
00277 
00278   }
00279 
00280 #ifndef FEI_SER
00281   double globalMaxErr = 0.0;
00282   MPI_Allreduce(&maxErr, &globalMaxErr, 1, MPI_DOUBLE, MPI_MAX, comm);
00283   maxErr = globalMaxErr;
00284 #endif
00285   bool testPassed = true;
00286   if (maxErr > 1.e-6) testPassed = false;
00287 
00288   if (testPassed && localProc == 0) {
00289     FEI_COUT << "poisson: TEST PASSED, maxErr = " << maxErr << ", iterations: "
00290    << itersTaken << FEI_ENDL;
00291     //This is something the SIERRA runtest tool looks for in test output...
00292     FEI_COUT << "SIERRA execution successful" << FEI_ENDL;
00293   }
00294   if (testPassed == false && localProc == 0) {
00295     FEI_COUT << "maxErr = " << maxErr << ", TEST FAILED" << FEI_ENDL;
00296     FEI_COUT << "(Test is deemed to have passed if the maximum difference"
00297    << " between the exact and computed solutions is 1.e-6 or less.)"
00298    << FEI_ENDL;
00299   }
00300 
00301   double elapsed_cpu_time = fei::utils::cpu_time() - start_time;
00302 
00303   //The following IOS_... macros are defined in base/fei_macros.h
00304   FEI_COUT.setf(IOS_FIXED, IOS_FLOATFIELD);
00305   if (localProc==0) {
00306     FEI_COUT << "Proc0 cpu times (seconds):" << FEI_ENDL
00307    << "   FEI initialize:  " << fei_init_time << FEI_ENDL
00308          << "   FEI load:        " << fei_load_time << FEI_ENDL
00309          << "      solve:        " << solve_time << FEI_ENDL
00310          << "Total program time: " << elapsed_cpu_time << FEI_ENDL;
00311   }
00312 
00313   wrapper.reset();
00314   fei.reset();
00315   factory.reset();
00316   //If Prometheus is being used, we need to make sure that the
00317   //LibraryWrapper inside factory is deleted before MPI_Finalize() is called.
00318   //(That's why we call the 'reset' methods on these shared-pointers rather
00319   //than just letting them destroy themselves when they go out of scope.)
00320 
00321   return(0);
00322 }
00323 

Generated on Tue Jul 13 09:27:46 2010 for FEI by  doxygen 1.4.7