FEI Version of the Day
beam_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 
00049 #include <fei_macros.hpp> //includes std stuff like iostream, etc.
00050 
00051 //Including the header fei_base.hpp gets us the declaration for
00052 //the various fei:: and snl_fei:: classes.
00053 
00054 #include <fei_base.hpp>
00055 
00056 //Now make provision for using any one of several solver libraries. This is
00057 //handled by the code in LibraryFactory.{hpp,cpp}.
00058 
00059 #include <test_utils/LibraryFactory.hpp>
00060 
00061 //And finally, we need to include some 'utils' headers which are used for
00062 //setting up the data for the example problem.
00063 
00064 #include <test_utils/HexBeam.hpp>
00065 #include <test_utils/HexBeamCR.hpp>
00066 #include <snl_fei_Utils.hpp>
00067 #include <test_utils/fei_test_utils.hpp>
00068 //
00069 //Include definitions of macros to call functions and check the return code.
00070 //
00071 #undef fei_file
00072 #define fei_file "cube3.cpp"
00073 #include <fei_ErrMacros.hpp>
00074 
00075 //==============================================================================
00076 //Here's the main...
00077 //==============================================================================
00078 int beam_main(int argc, char** argv,
00079          MPI_Comm comm, int numProcs, int localProc){
00080 
00081   double start_time = fei::utils::cpu_time();
00082 
00083   std::vector<std::string> stdstrings;
00084   CHK_ERR( fei_test_utils::get_filename_and_read_input(argc, argv,
00085             comm, localProc,
00086             stdstrings) );
00087   fei::ParameterSet params;
00088   fei::utils::parse_strings(stdstrings, " ", params);
00089 
00090   std::string solverName;
00091   std::string datasource;
00092   std::string constraintform;
00093   int W = 0;
00094   int D = 0;
00095   int DofPerNode = 0;
00096 
00097   int errcode = 0;
00098   errcode += params.getStringParamValue("SOLVER_LIBRARY", solverName);
00099   errcode += params.getStringParamValue("DATA_SOURCE", datasource);
00100   errcode += params.getIntParamValue("W", W);
00101   errcode += params.getIntParamValue("D", D);
00102   errcode += params.getIntParamValue("DofPerNode", DofPerNode);
00103 
00104   if (errcode != 0) {
00105     fei::console_out() << "Failed to find one or more required parameters in input-file."
00106        << FEI_ENDL << "Required parameters:"<<FEI_ENDL
00107        << "SOLVER_LIBRARY" << FEI_ENDL
00108        << "DATA_SOURCE" << FEI_ENDL
00109        << "CONSTRAINT_FORM" << FEI_ENDL
00110        << "W" << FEI_ENDL << "D" << FEI_ENDL << "DofPerNode" << FEI_ENDL;
00111     return(-1);
00112   }
00113 
00114   params.getStringParamValue("CONSTRAINT_FORM", constraintform);
00115 
00116   bool slave_constraints = false;
00117   if ("slaves" == constraintform) {
00118     slave_constraints = true;
00119   }
00120 
00121   HexBeam* hexcubeptr = NULL;
00122   if (datasource == "HexBeam") {
00123     hexcubeptr = new HexBeam(W, D, DofPerNode,
00124            HexBeam::OneD, numProcs, localProc);
00125   }
00126   else {
00127     hexcubeptr = new HexBeamCR(W, D, DofPerNode,
00128              HexBeam::OneD, numProcs, localProc);
00129   }
00130 
00131   HexBeam& hexcube = *hexcubeptr;
00132 
00133   if (localProc == 0) {
00134     int numCRs = (W+1)*(W+1)* ((numProcs*2)-1);
00135     if (hexcube.getNumCRs() < 1) numCRs = 0;
00136     FEI_COUT << FEI_ENDL;
00137     FEI_COUT << "========================================================" 
00138    << FEI_ENDL;
00139     FEI_COUT << "FEI version: " << fei::utils::version() << FEI_ENDL;
00140     FEI_COUT << "--------------------------------------------------------"<<FEI_ENDL;
00141     FEI_COUT << "Size W: " << W << " (num-elements-along-side-of-cube)"<<FEI_ENDL;
00142     FEI_COUT << "Size D: " << D << " (num-elements-along-depth-of-cube)"<<FEI_ENDL;
00143     FEI_COUT << "DOF per node: " << DofPerNode <<FEI_ENDL;
00144     FEI_COUT << "Num local  elements: " << hexcube.localNumElems_ << FEI_ENDL;
00145     FEI_COUT << "Num global elements: " << hexcube.totalNumElems_ << FEI_ENDL;
00146     FEI_COUT << "Num local  DOF: " << hexcube.numLocalDOF_ << FEI_ENDL;
00147     FEI_COUT << "Num global DOF: " << hexcube.numGlobalDOF_ << FEI_ENDL;
00148     FEI_COUT << "Num global CRs: " << numCRs << FEI_ENDL;
00149     FEI_COUT << "========================================================" 
00150    << FEI_ENDL;
00151   }
00152 
00153   double start_init_time = fei::utils::cpu_time();
00154 
00155   fei::SharedPtr<fei::Factory> factory;
00156   try {
00157     factory = fei::create_fei_Factory(comm, solverName.c_str());
00158   }
00159   catch (...) {
00160     FEI_COUT << "library " << solverName << " not available."<<FEI_ENDL;
00161     return(-1);
00162   }
00163   if (factory.get() == NULL) {
00164     FEI_COUT << "snl_fei::Factory allocation failed." << FEI_ENDL;
00165     return(-1);
00166   }
00167 
00168   factory->parameters(params);
00169 
00170   fei::SharedPtr<fei::VectorSpace> nodeSpace =
00171     factory->createVectorSpace(comm, "cube3");
00172 
00173   fei::SharedPtr<fei::VectorSpace> dummy;
00174   fei::SharedPtr<fei::MatrixGraph> matrixGraph =
00175     factory->createMatrixGraph(nodeSpace, dummy, "cube3");
00176 
00177   //load some solver-control parameters.
00178   nodeSpace->setParameters(params);
00179   matrixGraph->setParameters(params);
00180 
00181   int fieldID = 0;
00182   int fieldSize = hexcube.numDofPerNode();
00183   int nodeIDType = 0;
00184   int crIDType = 1;
00185 
00186   nodeSpace->defineFields( 1, &fieldID, &fieldSize );
00187   nodeSpace->defineIDTypes(1, &nodeIDType );
00188   nodeSpace->defineIDTypes(1, &crIDType );
00189 
00190   CHK_ERR( HexBeam_Functions::
00191      init_elem_connectivities( matrixGraph.get(), hexcube ) );
00192 
00193   CHK_ERR( HexBeam_Functions::
00194      init_shared_nodes( matrixGraph.get(), hexcube ) );
00195 
00196   int firstLocalCRID = 0;
00197   if (slave_constraints) {
00198     CHK_ERR( HexBeam_Functions::
00199        init_slave_constraints( matrixGraph.get(), hexcube) );
00200   }
00201   else {
00202     CHK_ERR( HexBeam_Functions::
00203        init_constraints( matrixGraph.get(), hexcube, localProc, firstLocalCRID) );
00204   }
00205 
00206   CHK_ERR( matrixGraph->initComplete() );
00207 
00208   double fei_init_time = fei::utils::cpu_time() - start_init_time;
00209 
00210   if (localProc == 0) {
00211     FEI_COUT.setf(IOS_FIXED, IOS_FLOATFIELD);
00212     FEI_COUT << "Initialization cpu time:   " << fei_init_time << FEI_ENDL;
00213   }
00214 
00215   //Now the initialization phase is complete. Next we'll do the load phase,
00216   //which for this problem just consists of loading the element data
00217   //(element-wise stiffness arrays and load vectors) and the boundary
00218   //condition data.
00219 
00220   double fei_creatematrix_start_time = fei::utils::cpu_time();
00221 
00222   fei::SharedPtr<fei::Matrix> mat = factory->createMatrix(matrixGraph);
00223 
00224   double fei_creatematrix_time = fei::utils::cpu_time() - fei_creatematrix_start_time;
00225   if (localProc == 0) {
00226     FEI_COUT.setf(IOS_FIXED, IOS_FLOATFIELD);
00227     FEI_COUT << "Create-Matrix cpu time:   " << fei_creatematrix_time << FEI_ENDL;
00228   }
00229 
00230   double start_load_time = fei::utils::cpu_time();
00231 
00232   fei::SharedPtr<fei::Vector> solnVec = factory->createVector(matrixGraph, true);
00233   fei::SharedPtr<fei::Vector> rhsVec  = factory->createVector(matrixGraph);
00234   fei::SharedPtr<fei::LinearSystem> linSys =
00235     factory->createLinearSystem(matrixGraph);
00236 
00237   CHK_ERR( linSys->parameters(params));
00238 
00239   linSys->setMatrix(mat);
00240   linSys->setSolutionVector(solnVec);
00241   linSys->setRHS(rhsVec);
00242 
00243   CHK_ERR( HexBeam_Functions::
00244      load_elem_data(matrixGraph.get(), mat.get(), rhsVec.get(), hexcube) );
00245 
00246   //temporarily disable boundary-conditions
00247   //  CHK_ERR( load_BC_data(linSys, hexcube) );
00248 
00249   if (!slave_constraints) {
00250     CHK_ERR( HexBeam_Functions::
00251        load_constraints(linSys.get(), hexcube, firstLocalCRID) );
00252   }
00253 
00254   CHK_ERR( linSys->loadComplete() );
00255 
00256   double fei_load_time = fei::utils::cpu_time() - start_load_time;
00257 
00258   if (localProc == 0) {
00259     //IOS macros are defined in fei_macros.h
00260     FEI_COUT.setf(IOS_FIXED, IOS_FLOATFIELD);
00261     FEI_COUT << "Coef. loading cpu time:    " << fei_load_time << FEI_ENDL;
00262     FEI_COUT << "Total assembly wall time:   "
00263    << fei_init_time + fei_creatematrix_time + fei_load_time << FEI_ENDL;
00264   }
00265 
00266   //
00267   //now the load phase is complete, so we're ready to launch the underlying
00268   //solver and solve Ax=b
00269   //
00270 
00271   //First, check whether the params (set from the input .i file) specify
00272   //a "Trilinos_Solver" string (which may be Amesos...). If not, it won't
00273   //matter but if so, the factory may switch based on this value, when
00274   //creating the solver object instance.
00275 
00276   std::string solver_name_value;
00277   params.getStringParamValue("Trilinos_Solver", solver_name_value);
00278 
00279   const char* charptr_solvername =
00280     solver_name_value.empty() ? 0 : solver_name_value.c_str();
00281 
00282   fei::SharedPtr<fei::Solver> solver = factory->createSolver(charptr_solvername);
00283 
00284   int status;
00285   int itersTaken = 0;
00286 
00287   if (localProc==0) FEI_COUT << "solve..." << FEI_ENDL;
00288   double start_solve_time = fei::utils::cpu_time();
00289 
00290   int err = solver->solve(linSys.get(),
00291         NULL, //preconditioningMatrix
00292         params,
00293         itersTaken,
00294         status);
00295 
00296   double solve_time = fei::utils::cpu_time()-start_solve_time;
00297 
00298   if (err!=0) {
00299     if (localProc==0) FEI_COUT << "solve returned err: " << err <<", status: "
00300          << status << FEI_ENDL;
00301   }
00302 
00303   if (localProc==0) {
00304     FEI_COUT << " cpu-time in solve: " << solve_time << FEI_ENDL;
00305   }
00306 
00307   CHK_ERR( solnVec->scatterToOverlap() );
00308 
00309   delete hexcubeptr;
00310 
00311   double elapsed_cpu_time = fei::utils::cpu_time() - start_time;
00312   int returnValue = 0;
00313 
00314   //The following IOS_... macros are defined in base/fei_macros.h
00315   FEI_COUT.setf(IOS_FIXED, IOS_FLOATFIELD);
00316   if (localProc==0) {
00317     FEI_COUT << "Proc0 cpu times (seconds):" << FEI_ENDL
00318    << "   FEI initialize:    " << fei_init_time << FEI_ENDL
00319          << "   FEI create-matrix: " << fei_creatematrix_time << FEI_ENDL
00320          << "   FEI load:          " << fei_load_time << FEI_ENDL
00321          << "      solve:          " << solve_time << FEI_ENDL
00322          << "Total program time:   " << elapsed_cpu_time << FEI_ENDL;
00323 
00324 #if defined(FEI_PLATFORM) && defined(FEI_OPT_LEVEL)
00325     double benchmark = fei_init_time;
00326 
00327     std::string slavestr;
00328     if (!constraintform.empty()) {
00329       slavestr = constraintform;
00330     }
00331     if (slavestr.size() > 0) slavestr += "_";
00332 
00333     FEI_OSTRINGSTREAM testname_init;
00334     testname_init << "cube3_init_"<<slavestr<<W<<"_"<<D<<"_"<<DofPerNode<<"_"
00335       <<solverName<<"_np"<<numProcs<<"_"
00336       <<FEI_PLATFORM<<"_"<<FEI_OPT_LEVEL;
00337     FEI_OSTRINGSTREAM testname_create;
00338     testname_create << "cube3_creatematrix_"<<slavestr<<W<<"_"<<D<<"_"<<DofPerNode
00339         <<"_"<<solverName<<"_np"<<numProcs<<"_"
00340         <<FEI_PLATFORM<<"_"<<FEI_OPT_LEVEL;
00341     FEI_OSTRINGSTREAM testname_load;
00342     testname_load << "cube3_load_"<<slavestr<<W<<"_"<<D<<"_"<<DofPerNode<<"_"
00343       <<solverName<<"_np"<<numProcs<<"_"
00344       <<FEI_PLATFORM<<"_"<<FEI_OPT_LEVEL;
00345 
00346     double file_init, file_create, file_load;
00347     bool file_benchmarks_available = true;
00348     try {
00349       file_init = fei_test_utils::get_file_benchmark("./cube3_timings.txt",
00350                 testname_init.str().c_str());
00351       file_create = fei_test_utils::get_file_benchmark("./cube3_timings.txt",
00352             testname_create.str().c_str());
00353       file_load = fei_test_utils::get_file_benchmark("./cube3_timings.txt",
00354                 testname_load.str().c_str());
00355     }
00356     catch (std::runtime_error& exc) {
00357       file_benchmarks_available = false;
00358     }
00359 
00360     if (file_benchmarks_available) {
00361 
00362       bool init_test_passed =
00363   fei_test_utils::check_and_cout_test_result(testname_init.str(), benchmark,
00364               file_init, 10);
00365 
00366       benchmark = fei_creatematrix_time;
00367       bool create_test_passed =
00368   fei_test_utils::check_and_cout_test_result(testname_create.str(), benchmark,
00369               file_create, 10);
00370 
00371       benchmark = fei_load_time;
00372       bool load_test_passed =
00373   fei_test_utils::check_and_cout_test_result(testname_load.str(), benchmark,
00374               file_load, 10);
00375 
00376       returnValue = init_test_passed&&create_test_passed&&load_test_passed ? 0 : 1;
00377     }
00378 #endif
00379 
00380   }
00381 
00382   bool testPassed = returnValue==0;
00383   if (testPassed && localProc == 0) {
00384     //A string for the SIERRA runtest tool to search for in test output...
00385     FEI_COUT << "FEI test successful" << FEI_ENDL;
00386   }
00387 
00388   return(returnValue);
00389 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends