poisson3_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 FEI classes,
00011 // for the purposes of 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, and only one element-block per processor.
00016 //
00017 // This problem was coded up with Ray Tuminaro.
00018 //
00019 // The input file for this program should provide the following:
00020 //   SOLVER_LIBRARY <library-name> -- e.g., Aztec
00021 //   L <int> -- the global length (num-elements) on a side of the 2D square 
00022 //
00023 // Alan Williams 03-13-2002
00024 //
00025 
00026 #include <fei_iostream.hpp>
00027 #include <cmath>
00028 
00029 //Including the header fei_base.hpp gets us the declaration for
00030 //the various FEI classes.
00031 
00032 #include <fei_base.hpp>
00033 
00034 //Now make provision for using any one of several solver libraries. This is
00035 //handled by the code in LibraryFactory.[hC].
00036 
00037 #include <test_utils/LibraryFactory.hpp>
00038 
00039 //And, we need to include some headers for utility classes which are simply
00040 //used for setting up the data for the example problem.
00041 
00042 #include <test_utils/Poisson_Elem.hpp>
00043 #include <test_utils/PoissonData.hpp>
00044 
00045 #include <test_utils/ElemBlock.hpp>
00046 #include <test_utils/CRSet.hpp>
00047 #include <test_utils/CommNodeSet.hpp>
00048 #include <test_utils/DataReader.hpp>
00049 
00050 #include <test_utils/SolnCheck.hpp>
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 #include <fei_ErrMacros.hpp>
00057 
00058 
00059 //==============================================================================
00060 //Here's the main...
00061 //==============================================================================
00062 int poisson3_main(int argc, char** argv,
00063                   MPI_Comm comm, int numProcs, int localProc){
00064   int masterProc = 0;
00065 
00066   double start_time = fei::utils::cpu_time();
00067 
00068   std::vector<std::string> stdstrings;
00069   CHK_ERR( fei_test_utils::get_filename_and_read_input(argc, argv,
00070             comm, localProc,
00071             stdstrings) );
00072   const char** params = NULL;
00073   int numParams = 0;
00074   fei::utils::strings_to_char_ptrs(stdstrings, numParams, params);
00075   
00076   fei::ParameterSet paramset;
00077   fei::utils::parse_strings(stdstrings, " ", paramset);
00078 
00079   std::string solverName;
00080   int L = 0;
00081   int outputLevel = 0;
00082 
00083   int errcode = 0;
00084   errcode += paramset.getStringParamValue("SOLVER_LIBRARY", solverName);
00085   errcode += paramset.getIntParamValue("L", L);
00086   paramset.getIntParamValue("outputLevel", outputLevel);
00087 
00088   if (errcode != 0) {
00089     FEI_CERR << "Failed to find one or more required parameters in input-file."
00090        << FEI_ENDL << "Required parameters:"<<FEI_ENDL
00091        << "SOLVER_LIBRARY" << FEI_ENDL
00092        << "L" << FEI_ENDL;
00093     return(-1);
00094   }
00095 
00096   if (localProc == 0) {
00097     int nodes = (L+1)*(L+1);
00098     int eqns = nodes;
00099     FEI_COUT << FEI_ENDL;
00100     FEI_COUT << "========================================================" 
00101    << FEI_ENDL;
00102     FEI_COUT << "FEI version: " << fei::utils::version() << FEI_ENDL;
00103     FEI_COUT << "Square size     L: " << L << " elements." << FEI_ENDL;
00104     FEI_COUT << "Global number of elements: " << L*L << FEI_ENDL;
00105     FEI_COUT << "Global number of nodes: " << nodes << FEI_ENDL;
00106     FEI_COUT << "Global number of equations: " << eqns <<FEI_ENDL;
00107     FEI_COUT << "========================================================" 
00108    << FEI_ENDL;
00109   }
00110 
00111   if ((masterProc == localProc)&&(outputLevel>0)) {
00112     fei_test_utils::print_args(argc, argv);
00113   }
00114 
00115   if (outputLevel == 1) {
00116     if (localProc != 0) outputLevel = 0;
00117   }
00118 
00119   //PoissonData is the object that will be in charge of generating the
00120   //data to pump into the FEI objects.
00121 
00122   PoissonData poissonData(L, numProcs, localProc, outputLevel);
00123 
00124   double start_init_time = fei::utils::cpu_time();
00125 
00126   fei::SharedPtr<fei::Factory> factory;
00127   try {
00128     factory = fei::create_fei_Factory(comm, solverName.c_str());
00129   }
00130   catch (...) {
00131     FEI_COUT << "library " << solverName << " not available."<<FEI_ENDL;
00132     return(-1);
00133   }
00134 
00135   if (factory.get() == NULL) {
00136     FEI_COUT << "snl_fei::Factory creation failed." << FEI_ENDL;
00137     return(-1);
00138   }
00139 
00140   factory->parameters(paramset);
00141 
00142   fei::SharedPtr<fei::VectorSpace> nodeSpace =
00143     factory->createVectorSpace(comm, "poisson3");
00144 
00145   fei::SharedPtr<fei::VectorSpace> dummy;
00146   fei::SharedPtr<fei::MatrixGraph> matrixGraph =
00147     factory->createMatrixGraph(nodeSpace, dummy, "poisson3");
00148 
00149   //load some control parameters.
00150   matrixGraph->setParameters(paramset);
00151 
00152 
00153   int numFields = poissonData.getNumFields();
00154   int* fieldSizes = poissonData.getFieldSizes();
00155   int* fieldIDs = poissonData.getFieldIDs();
00156   int nodeIDType = 0;
00157 
00158   if (outputLevel>0 && localProc==0) FEI_COUT << "defineFields" << FEI_ENDL;
00159   nodeSpace->defineFields( numFields, fieldIDs, fieldSizes );
00160 
00161   if (outputLevel>0 && localProc==0) FEI_COUT << "defineIDTypes" << FEI_ENDL;
00162   nodeSpace->defineIDTypes( 1, &nodeIDType );
00163 
00164   CHK_ERR( init_elem_connectivities(matrixGraph.get(), poissonData) );
00165 
00166   CHK_ERR( set_shared_nodes(nodeSpace.get(), poissonData) );
00167 
00168 
00169   //The following IOS_... macros are defined in base/fei_macros.h
00170   FEI_COUT.setf(IOS_FIXED, IOS_FLOATFIELD);
00171   if (outputLevel>0 && localProc==0) FEI_COUT << "initComplete" << FEI_ENDL;
00172   CHK_ERR( matrixGraph->initComplete() );
00173 
00174   double fei_init_time = fei::utils::cpu_time() - start_init_time;
00175 
00176   //Now the initialization phase is complete. Next we'll do the load phase,
00177   //which for this problem just consists of loading the element data
00178   //(element-wise stiffness arrays and load vectors) and the boundary
00179   //condition data.
00180   //This simple problem doesn't have any constraint relations, etc.
00181 
00182   double start_load_time = fei::utils::cpu_time();
00183 
00184 
00185   fei::SharedPtr<fei::Matrix> mat = factory->createMatrix(matrixGraph);
00186   fei::SharedPtr<fei::Vector> solnVec = factory->createVector(nodeSpace, true);
00187   fei::SharedPtr<fei::Vector> rhsVec  = factory->createVector(nodeSpace);
00188   fei::SharedPtr<fei::LinearSystem> linSys= factory->createLinearSystem(matrixGraph);
00189 
00190   linSys->setMatrix(mat);
00191   linSys->setSolutionVector(solnVec);
00192   linSys->setRHS(rhsVec);
00193 
00194   CHK_ERR( linSys->parameters(paramset));
00195   CHK_ERR( load_elem_data(matrixGraph.get(), mat.get(), rhsVec.get(), poissonData) );
00196 
00197   CHK_ERR( load_BC_data(linSys.get(), poissonData) );
00198 
00199   CHK_ERR( linSys->loadComplete() );
00200 
00201   double fei_load_time = fei::utils::cpu_time() - start_load_time;
00202 
00203   //
00204   //now the load phase is complete, so we're ready to launch the underlying
00205   //solver and solve Ax=b
00206   //
00207 
00208   fei::SharedPtr<fei::Solver> solver = factory->createSolver();
00209 
00210   int status;
00211   int itersTaken = 0;
00212 
00213   if (outputLevel>0 && localProc==0) FEI_COUT << "solve..." << FEI_ENDL;
00214   double start_solve_time = fei::utils::cpu_time();
00215 
00216   int err = solver->solve(linSys.get(),
00217         NULL, //preconditioningMatrix
00218         paramset, itersTaken, status);
00219 
00220   double solve_time = fei::utils::cpu_time() - start_solve_time;
00221 
00222   if (err!=0) {
00223     if (localProc==0) FEI_COUT << "solve returned err: " << err <<", status: "
00224          << status << FEI_ENDL;
00225   }
00226 
00227   CHK_ERR( solnVec->scatterToOverlap() );
00228 
00229   //
00230   //We oughtta make sure the solution we just computed is correct...
00231   //
00232 
00233   int numNodes = nodeSpace->getNumOwnedAndSharedIDs(nodeIDType);
00234 
00235   double maxErr = 0.0;
00236   if (numNodes > 0) {
00237     int lenNodeIDs = numNodes;
00238     GlobalID* nodeIDs = new GlobalID[lenNodeIDs];
00239     double* soln = new double[lenNodeIDs];
00240     if (nodeIDs != NULL && soln != NULL) {
00241       CHK_ERR( nodeSpace->getOwnedAndSharedIDs(nodeIDType, numNodes,
00242               nodeIDs, lenNodeIDs) );
00243 
00244       int fieldID = 1;
00245       CHK_ERR( solnVec->copyOutFieldData(fieldID, nodeIDType,
00246               numNodes, nodeIDs, soln));
00247 
00248       for(int i=0; i<numNodes; i++) {
00249   int nID = (int)nodeIDs[i];
00250   double x = (1.0* ((nID-1)%(L+1)))/L;
00251   double y = (1.0* ((nID-1)/(L+1)))/L;
00252 
00253   double exactSoln = x*x + y*y;
00254   double error = std::abs(exactSoln - soln[i]);
00255   if (maxErr < error) maxErr = error;
00256       }
00257 
00258       delete [] nodeIDs;
00259       delete [] soln;
00260     }
00261     else {
00262       FEI_CERR << "allocation of nodeIDs or soln failed." << FEI_ENDL; 
00263     }
00264 
00265   }
00266 
00267 #ifndef FEI_SER
00268   double globalMaxErr = 0.0;
00269   MPI_Allreduce(&maxErr, &globalMaxErr, 1, MPI_DOUBLE, MPI_MAX, comm);
00270   maxErr = globalMaxErr;
00271 #endif
00272   bool testPassed = true;
00273   if (maxErr > 1.e-6) testPassed = false;
00274 
00275   double elapsed_cpu_time = fei::utils::cpu_time() - start_time;
00276   int returnValue = 0;
00277 
00278   //The following IOS_... macros are defined in base/fei_macros.h
00279   FEI_COUT.setf(IOS_FIXED, IOS_FLOATFIELD);
00280   if (localProc==0) {
00281     FEI_COUT << "Proc0 cpu times (seconds):" << FEI_ENDL
00282    << "   FEI initialize:  " << fei_init_time << FEI_ENDL
00283          << "   FEI load:        " << fei_load_time << FEI_ENDL
00284          << "      solve:        " << solve_time << FEI_ENDL
00285          << "Total program time: " << elapsed_cpu_time << FEI_ENDL;
00286 
00287 #if defined(FEI_PLATFORM) && defined(FEI_OPT_LEVEL)
00288     double benchmark = fei_init_time+fei_load_time;
00289     FEI_OSTRINGSTREAM testname;
00290     testname << "poisson3_"<<L<<"_"<<solverName<<"_np"<<numProcs<<"_"
00291        <<FEI_PLATFORM<<"_"<<FEI_OPT_LEVEL;
00292  
00293     double file_value = 0.0;
00294     bool file_value_available = true;
00295     try {
00296       file_value = fei_test_utils::get_file_benchmark("./poisson3_timings.txt",
00297                  testname.str().c_str());
00298     }
00299     catch(std::runtime_error& exc) {
00300       file_value_available = false;
00301     }
00302 
00303     if (file_value_available) {
00304       bool test_passed = fei_test_utils::check_and_cout_test_result(testname.str(),
00305                    benchmark,
00306                    file_value, 10);
00307       returnValue = test_passed ? 0 : 1;
00308     }
00309 #endif
00310 
00311   }
00312 
00313   if (testPassed && returnValue==0 && localProc == 0) {
00314     FEI_COUT.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
00315     FEI_COUT << "poisson: TEST PASSED, maxErr = " << maxErr << ", iterations: "
00316    << itersTaken << FEI_ENDL;
00317     //This is something the SIERRA runtest tool looks for in test output...
00318     FEI_COUT << "SIERRA execution successful" << FEI_ENDL;
00319   }
00320   if ((testPassed == false || returnValue != 0) && localProc == 0) {
00321     if (testPassed==true) {
00322       FEI_COUT << "maxErr = " << maxErr << ", but time-taken outside margin. TEST FAILED" << FEI_ENDL;
00323     }
00324     else {
00325       FEI_COUT << "maxErr = " << maxErr << ", TEST FAILED"<<FEI_ENDL;
00326     }
00327     FEI_COUT << "(Test is deemed to have passed if the maximum difference"
00328    << " between the exact and computed solutions is 1.e-6 or less, *AND*"
00329    << " time-taken matches file-benchmark if available.)"
00330    << FEI_ENDL;
00331   }
00332 
00333   return(returnValue);
00334 }
00335 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends
Generated on Wed Apr 13 10:08:24 2011 for FEI by  doxygen 1.6.3