FEI Version of the Day
poisson3_main.cpp
00001 /*
00002 // @HEADER
00003 // ************************************************************************
00004 //             FEI: Finite Element Interface to Linear Solvers
00005 //                  Copyright (2005) Sandia Corporation.
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
00008 // U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Alan Williams (william@sandia.gov) 
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 */
00042 
00043 
00044 //
00045 // This is a simple program to exercise FEI classes,
00046 // for the purposes of testing, code tuning and scaling studies.
00047 //
00048 // This program assembles a linear system from a 2D square Poisson
00049 // problem, using 4-node square elements. There is only 1 degree-of-
00050 // freedom per node, and only one element-block per processor.
00051 //
00052 // This problem was coded up with Ray Tuminaro.
00053 //
00054 // The input file for this program should provide the following:
00055 //   SOLVER_LIBRARY <library-name> -- e.g., Aztec
00056 //   L <int> -- the global length (num-elements) on a side of the 2D square 
00057 //
00058 // Alan Williams 03-13-2002
00059 //
00060 
00061 #include <fei_iostream.hpp>
00062 #include <cmath>
00063 
00064 //Including the header fei_base.hpp gets us the declaration for
00065 //the various FEI classes.
00066 
00067 #include <fei_base.hpp>
00068 
00069 //Now make provision for using any one of several solver libraries. This is
00070 //handled by the code in LibraryFactory.[hC].
00071 
00072 #include <test_utils/LibraryFactory.hpp>
00073 
00074 //And, we need to include some headers for utility classes which are simply
00075 //used for setting up the data for the example problem.
00076 
00077 #include <test_utils/Poisson_Elem.hpp>
00078 #include <test_utils/PoissonData.hpp>
00079 
00080 #include <test_utils/ElemBlock.hpp>
00081 #include <test_utils/CRSet.hpp>
00082 #include <test_utils/CommNodeSet.hpp>
00083 #include <test_utils/DataReader.hpp>
00084 
00085 #include <test_utils/SolnCheck.hpp>
00086 #include <test_utils/fei_test_utils.hpp>
00087 #include <snl_fei_Utils.hpp>
00088 //
00089 //Include definitions of macros to call functions and check the return code.
00090 //
00091 #include <fei_ErrMacros.hpp>
00092 
00093 
00094 //==============================================================================
00095 //Here's the main...
00096 //==============================================================================
00097 int poisson3_main(int argc, char** argv,
00098                   MPI_Comm comm, int numProcs, int localProc){
00099   int masterProc = 0;
00100 
00101   double start_time = fei::utils::cpu_time();
00102 
00103   std::vector<std::string> stdstrings;
00104   CHK_ERR( fei_test_utils::get_filename_and_read_input(argc, argv,
00105             comm, localProc,
00106             stdstrings) );
00107   const char** params = NULL;
00108   int numParams = 0;
00109   fei::utils::strings_to_char_ptrs(stdstrings, numParams, params);
00110   
00111   fei::ParameterSet paramset;
00112   fei::utils::parse_strings(stdstrings, " ", paramset);
00113 
00114   std::string solverName;
00115   int L = 0;
00116   int outputLevel = 0;
00117 
00118   int errcode = 0;
00119   errcode += paramset.getStringParamValue("SOLVER_LIBRARY", solverName);
00120   errcode += paramset.getIntParamValue("L", L);
00121   paramset.getIntParamValue("outputLevel", outputLevel);
00122 
00123   if (errcode != 0) {
00124     fei::console_out() << "Failed to find one or more required parameters in input-file."
00125        << FEI_ENDL << "Required parameters:"<<FEI_ENDL
00126        << "SOLVER_LIBRARY" << FEI_ENDL
00127        << "L" << FEI_ENDL;
00128     return(-1);
00129   }
00130 
00131   if (localProc == 0) {
00132     int nodes = (L+1)*(L+1);
00133     int eqns = nodes;
00134     FEI_COUT << FEI_ENDL;
00135     FEI_COUT << "========================================================" 
00136    << FEI_ENDL;
00137     FEI_COUT << "FEI version: " << fei::utils::version() << FEI_ENDL;
00138     FEI_COUT << "Square size     L: " << L << " elements." << FEI_ENDL;
00139     FEI_COUT << "Global number of elements: " << L*L << FEI_ENDL;
00140     FEI_COUT << "Global number of nodes: " << nodes << FEI_ENDL;
00141     FEI_COUT << "Global number of equations: " << eqns <<FEI_ENDL;
00142     FEI_COUT << "========================================================" 
00143    << FEI_ENDL;
00144   }
00145 
00146   if ((masterProc == localProc)&&(outputLevel>0)) {
00147     fei_test_utils::print_args(argc, argv);
00148   }
00149 
00150   if (outputLevel == 1) {
00151     if (localProc != 0) outputLevel = 0;
00152   }
00153 
00154   //PoissonData is the object that will be in charge of generating the
00155   //data to pump into the FEI objects.
00156 
00157   PoissonData poissonData(L, numProcs, localProc, outputLevel);
00158 
00159   double start_init_time = fei::utils::cpu_time();
00160 
00161   fei::SharedPtr<fei::Factory> factory;
00162   try {
00163     factory = fei::create_fei_Factory(comm, solverName.c_str());
00164   }
00165   catch (...) {
00166     FEI_COUT << "library " << solverName << " not available."<<FEI_ENDL;
00167     return(-1);
00168   }
00169 
00170   if (factory.get() == NULL) {
00171     FEI_COUT << "snl_fei::Factory creation failed." << FEI_ENDL;
00172     return(-1);
00173   }
00174 
00175   factory->parameters(paramset);
00176 
00177   fei::SharedPtr<fei::VectorSpace> nodeSpace =
00178     factory->createVectorSpace(comm, "poisson3");
00179 
00180   fei::SharedPtr<fei::VectorSpace> dummy;
00181   fei::SharedPtr<fei::MatrixGraph> matrixGraph =
00182     factory->createMatrixGraph(nodeSpace, dummy, "poisson3");
00183 
00184   //load some control parameters.
00185   matrixGraph->setParameters(paramset);
00186 
00187 
00188   int numFields = poissonData.getNumFields();
00189   int* fieldSizes = poissonData.getFieldSizes();
00190   int* fieldIDs = poissonData.getFieldIDs();
00191   int nodeIDType = 0;
00192 
00193   if (outputLevel>0 && localProc==0) FEI_COUT << "defineFields" << FEI_ENDL;
00194   nodeSpace->defineFields( numFields, fieldIDs, fieldSizes );
00195 
00196   if (outputLevel>0 && localProc==0) FEI_COUT << "defineIDTypes" << FEI_ENDL;
00197   nodeSpace->defineIDTypes( 1, &nodeIDType );
00198 
00199   CHK_ERR( init_elem_connectivities(matrixGraph.get(), poissonData) );
00200 
00201   CHK_ERR( set_shared_nodes(nodeSpace.get(), poissonData) );
00202 
00203 
00204   //The following IOS_... macros are defined in base/fei_macros.h
00205   FEI_COUT.setf(IOS_FIXED, IOS_FLOATFIELD);
00206   if (outputLevel>0 && localProc==0) FEI_COUT << "initComplete" << FEI_ENDL;
00207   CHK_ERR( matrixGraph->initComplete() );
00208 
00209   double fei_init_time = fei::utils::cpu_time() - start_init_time;
00210 
00211   //Now the initialization phase is complete. Next we'll do the load phase,
00212   //which for this problem just consists of loading the element data
00213   //(element-wise stiffness arrays and load vectors) and the boundary
00214   //condition data.
00215   //This simple problem doesn't have any constraint relations, etc.
00216 
00217   double start_load_time = fei::utils::cpu_time();
00218 
00219 
00220   fei::SharedPtr<fei::Matrix> mat = factory->createMatrix(matrixGraph);
00221   fei::SharedPtr<fei::Vector> solnVec = factory->createVector(nodeSpace, true);
00222   fei::SharedPtr<fei::Vector> rhsVec  = factory->createVector(nodeSpace);
00223   fei::SharedPtr<fei::LinearSystem> linSys= factory->createLinearSystem(matrixGraph);
00224 
00225   linSys->setMatrix(mat);
00226   linSys->setSolutionVector(solnVec);
00227   linSys->setRHS(rhsVec);
00228 
00229   CHK_ERR( linSys->parameters(paramset));
00230   CHK_ERR( load_elem_data(matrixGraph.get(), mat.get(), rhsVec.get(), poissonData) );
00231 
00232   CHK_ERR( load_BC_data(linSys.get(), poissonData) );
00233 
00234   CHK_ERR( linSys->loadComplete() );
00235 
00236   double fei_load_time = fei::utils::cpu_time() - start_load_time;
00237 
00238   //
00239   //now the load phase is complete, so we're ready to launch the underlying
00240   //solver and solve Ax=b
00241   //
00242 
00243   fei::SharedPtr<fei::Solver> solver = factory->createSolver();
00244 
00245   int status;
00246   int itersTaken = 0;
00247 
00248   if (outputLevel>0 && localProc==0) FEI_COUT << "solve..." << FEI_ENDL;
00249   double start_solve_time = fei::utils::cpu_time();
00250 
00251   int err = solver->solve(linSys.get(),
00252         NULL, //preconditioningMatrix
00253         paramset, itersTaken, status);
00254 
00255   double solve_time = fei::utils::cpu_time() - start_solve_time;
00256 
00257   if (err!=0) {
00258     if (localProc==0) FEI_COUT << "solve returned err: " << err <<", status: "
00259          << status << FEI_ENDL;
00260   }
00261 
00262   CHK_ERR( solnVec->scatterToOverlap() );
00263 
00264   //
00265   //We oughtta make sure the solution we just computed is correct...
00266   //
00267 
00268   int numNodes = nodeSpace->getNumOwnedAndSharedIDs(nodeIDType);
00269 
00270   double maxErr = 0.0;
00271   if (numNodes > 0) {
00272     int lenNodeIDs = numNodes;
00273     GlobalID* nodeIDs = new GlobalID[lenNodeIDs];
00274     double* soln = new double[lenNodeIDs];
00275     if (nodeIDs != NULL && soln != NULL) {
00276       CHK_ERR( nodeSpace->getOwnedAndSharedIDs(nodeIDType, numNodes,
00277               nodeIDs, lenNodeIDs) );
00278 
00279       int fieldID = 1;
00280       CHK_ERR( solnVec->copyOutFieldData(fieldID, nodeIDType,
00281               numNodes, nodeIDs, soln));
00282 
00283       for(int i=0; i<numNodes; i++) {
00284   int nID = (int)nodeIDs[i];
00285   double x = (1.0* ((nID-1)%(L+1)))/L;
00286   double y = (1.0* ((nID-1)/(L+1)))/L;
00287 
00288   double exactSoln = x*x + y*y;
00289   double error = std::abs(exactSoln - soln[i]);
00290   if (maxErr < error) maxErr = error;
00291       }
00292 
00293       delete [] nodeIDs;
00294       delete [] soln;
00295     }
00296     else {
00297       fei::console_out() << "allocation of nodeIDs or soln failed." << FEI_ENDL; 
00298     }
00299 
00300   }
00301 
00302 #ifndef FEI_SER
00303   double globalMaxErr = 0.0;
00304   MPI_Allreduce(&maxErr, &globalMaxErr, 1, MPI_DOUBLE, MPI_MAX, comm);
00305   maxErr = globalMaxErr;
00306 #endif
00307   bool testPassed = true;
00308   if (maxErr > 1.e-6) testPassed = false;
00309 
00310   double elapsed_cpu_time = fei::utils::cpu_time() - start_time;
00311   int returnValue = 0;
00312 
00313   //The following IOS_... macros are defined in base/fei_macros.h
00314   FEI_COUT.setf(IOS_FIXED, IOS_FLOATFIELD);
00315   if (localProc==0) {
00316     FEI_COUT << "Proc0 cpu times (seconds):" << FEI_ENDL
00317    << "   FEI initialize:  " << fei_init_time << FEI_ENDL
00318          << "   FEI load:        " << fei_load_time << FEI_ENDL
00319          << "      solve:        " << solve_time << FEI_ENDL
00320          << "Total program time: " << elapsed_cpu_time << FEI_ENDL;
00321 
00322 #if defined(FEI_PLATFORM) && defined(FEI_OPT_LEVEL)
00323     double benchmark = fei_init_time+fei_load_time;
00324     FEI_OSTRINGSTREAM testname;
00325     testname << "poisson3_"<<L<<"_"<<solverName<<"_np"<<numProcs<<"_"
00326        <<FEI_PLATFORM<<"_"<<FEI_OPT_LEVEL;
00327  
00328     double file_value = 0.0;
00329     bool file_value_available = true;
00330     try {
00331       file_value = fei_test_utils::get_file_benchmark("./poisson3_timings.txt",
00332                  testname.str().c_str());
00333     }
00334     catch(std::runtime_error& exc) {
00335       file_value_available = false;
00336     }
00337 
00338     if (file_value_available) {
00339       bool test_passed = fei_test_utils::check_and_cout_test_result(testname.str(),
00340                    benchmark,
00341                    file_value, 10);
00342       returnValue = test_passed ? 0 : 1;
00343     }
00344 #endif
00345 
00346   }
00347 
00348   if (testPassed && returnValue==0 && localProc == 0) {
00349     FEI_COUT.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
00350     FEI_COUT << "poisson: TEST PASSED, maxErr = " << maxErr << ", iterations: "
00351    << itersTaken << FEI_ENDL;
00352     //This is something the SIERRA runtest tool looks for in test output...
00353     FEI_COUT << "SIERRA execution successful" << FEI_ENDL;
00354   }
00355   if ((testPassed == false || returnValue != 0) && localProc == 0) {
00356     if (testPassed==true) {
00357       FEI_COUT << "maxErr = " << maxErr << ", but time-taken outside margin. TEST FAILED" << FEI_ENDL;
00358     }
00359     else {
00360       FEI_COUT << "maxErr = " << maxErr << ", TEST FAILED"<<FEI_ENDL;
00361     }
00362     FEI_COUT << "(Test is deemed to have passed if the maximum difference"
00363    << " between the exact and computed solutions is 1.e-6 or less, *AND*"
00364    << " time-taken matches file-benchmark if available.)"
00365    << FEI_ENDL;
00366   }
00367 
00368   return(returnValue);
00369 }
00370 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends