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