FEI Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
beam_main.cpp
Go to the documentation of this file.
1 /*--------------------------------------------------------------------*/
2 /* Copyright 2005 Sandia Corporation. */
3 /* Under the terms of Contract DE-AC04-94AL85000, there is a */
4 /* non-exclusive license for use of this work by or on behalf */
5 /* of the U.S. Government. Export of this program may require */
6 /* a license from the United States Government. */
7 /*--------------------------------------------------------------------*/
8 
9 //
10 // This is a simple program to exercise FEI classes,
11 // for the purposes of testing, code tuning and scaling studies.
12 //
13 
14 #include <fei_macros.hpp> //includes std stuff like iostream, etc.
15 
16 //Including the header fei_base.hpp gets us the declaration for
17 //the various fei:: and snl_fei:: classes.
18 
19 #include <fei_base.hpp>
20 
21 //Now make provision for using any one of several solver libraries. This is
22 //handled by the code in LibraryFactory.{hpp,cpp}.
23 
25 
26 //And finally, we need to include some 'utils' headers which are used for
27 //setting up the data for the example problem.
28 
29 #include <test_utils/HexBeam.hpp>
30 #include <test_utils/HexBeamCR.hpp>
31 #include <snl_fei_Utils.hpp>
33 //
34 //Include definitions of macros to call functions and check the return code.
35 //
36 #undef fei_file
37 #define fei_file "cube3.cpp"
38 #include <fei_ErrMacros.hpp>
39 
40 //==============================================================================
41 //Here's the main...
42 //==============================================================================
43 int beam_main(int argc, char** argv,
44  MPI_Comm comm, int numProcs, int localProc){
45 
46  double start_time = fei::utils::cpu_time();
47 
48  std::vector<std::string> stdstrings;
50  comm, localProc,
51  stdstrings) );
52  fei::ParameterSet params;
53  fei::utils::parse_strings(stdstrings, " ", params);
54 
55  std::string solverName;
56  std::string datasource;
57  std::string constraintform;
58  int W = 0;
59  int D = 0;
60  int DofPerNode = 0;
61 
62  int errcode = 0;
63  errcode += params.getStringParamValue("SOLVER_LIBRARY", solverName);
64  errcode += params.getStringParamValue("DATA_SOURCE", datasource);
65  errcode += params.getIntParamValue("W", W);
66  errcode += params.getIntParamValue("D", D);
67  errcode += params.getIntParamValue("DofPerNode", DofPerNode);
68 
69  if (errcode != 0) {
70  fei::console_out() << "Failed to find one or more required parameters in input-file."
71  << FEI_ENDL << "Required parameters:"<<FEI_ENDL
72  << "SOLVER_LIBRARY" << FEI_ENDL
73  << "DATA_SOURCE" << FEI_ENDL
74  << "CONSTRAINT_FORM" << FEI_ENDL
75  << "W" << FEI_ENDL << "D" << FEI_ENDL << "DofPerNode" << FEI_ENDL;
76  return(-1);
77  }
78 
79  params.getStringParamValue("CONSTRAINT_FORM", constraintform);
80 
81  bool slave_constraints = false;
82  if ("slaves" == constraintform) {
83  slave_constraints = true;
84  }
85 
86  HexBeam* hexcubeptr = NULL;
87  if (datasource == "HexBeam") {
88  hexcubeptr = new HexBeam(W, D, DofPerNode,
89  HexBeam::OneD, numProcs, localProc);
90  }
91  else {
92  hexcubeptr = new HexBeamCR(W, D, DofPerNode,
93  HexBeam::OneD, numProcs, localProc);
94  }
95 
96  HexBeam& hexcube = *hexcubeptr;
97 
98  if (localProc == 0) {
99  int numCRs = (W+1)*(W+1)* ((numProcs*2)-1);
100  if (hexcube.getNumCRs() < 1) numCRs = 0;
101  FEI_COUT << FEI_ENDL;
102  FEI_COUT << "========================================================"
103  << FEI_ENDL;
104  FEI_COUT << "FEI version: " << fei::utils::version() << FEI_ENDL;
105  FEI_COUT << "--------------------------------------------------------"<<FEI_ENDL;
106  FEI_COUT << "Size W: " << W << " (num-elements-along-side-of-cube)"<<FEI_ENDL;
107  FEI_COUT << "Size D: " << D << " (num-elements-along-depth-of-cube)"<<FEI_ENDL;
108  FEI_COUT << "DOF per node: " << DofPerNode <<FEI_ENDL;
109  FEI_COUT << "Num local elements: " << hexcube.localNumElems_ << FEI_ENDL;
110  FEI_COUT << "Num global elements: " << hexcube.totalNumElems_ << FEI_ENDL;
111  FEI_COUT << "Num local DOF: " << hexcube.numLocalDOF_ << FEI_ENDL;
112  FEI_COUT << "Num global DOF: " << hexcube.numGlobalDOF_ << FEI_ENDL;
113  FEI_COUT << "Num global CRs: " << numCRs << FEI_ENDL;
114  FEI_COUT << "========================================================"
115  << FEI_ENDL;
116  }
117 
118  double start_init_time = fei::utils::cpu_time();
119 
121  try {
122  factory = fei::create_fei_Factory(comm, solverName.c_str());
123  }
124  catch (...) {
125  FEI_COUT << "library " << solverName << " not available."<<FEI_ENDL;
126  return(-1);
127  }
128  if (factory.get() == NULL) {
129  FEI_COUT << "snl_fei::Factory allocation failed." << FEI_ENDL;
130  return(-1);
131  }
132 
133  factory->parameters(params);
134 
136  factory->createVectorSpace(comm, "cube3");
137 
140  factory->createMatrixGraph(nodeSpace, dummy, "cube3");
141 
142  //load some solver-control parameters.
143  nodeSpace->setParameters(params);
144  matrixGraph->setParameters(params);
145 
146  int fieldID = 0;
147  int fieldSize = hexcube.numDofPerNode();
148  int nodeIDType = 0;
149  int crIDType = 1;
150 
151  nodeSpace->defineFields( 1, &fieldID, &fieldSize );
152  nodeSpace->defineIDTypes(1, &nodeIDType );
153  nodeSpace->defineIDTypes(1, &crIDType );
154 
155  CHK_ERR( HexBeam_Functions::
156  init_elem_connectivities( matrixGraph.get(), hexcube ) );
157 
158  CHK_ERR( HexBeam_Functions::
159  init_shared_nodes( matrixGraph.get(), hexcube ) );
160 
161  int firstLocalCRID = 0;
162  if (slave_constraints) {
163  CHK_ERR( HexBeam_Functions::
164  init_slave_constraints( matrixGraph.get(), hexcube) );
165  }
166  else {
167  CHK_ERR( HexBeam_Functions::
168  init_constraints( matrixGraph.get(), hexcube, localProc, firstLocalCRID) );
169  }
170 
171  CHK_ERR( matrixGraph->initComplete() );
172 
173  double fei_init_time = fei::utils::cpu_time() - start_init_time;
174 
175  if (localProc == 0) {
177  FEI_COUT << "Initialization cpu time: " << fei_init_time << FEI_ENDL;
178  }
179 
180  //Now the initialization phase is complete. Next we'll do the load phase,
181  //which for this problem just consists of loading the element data
182  //(element-wise stiffness arrays and load vectors) and the boundary
183  //condition data.
184 
185  double fei_creatematrix_start_time = fei::utils::cpu_time();
186 
187  fei::SharedPtr<fei::Matrix> mat = factory->createMatrix(matrixGraph);
188 
189  double fei_creatematrix_time = fei::utils::cpu_time() - fei_creatematrix_start_time;
190  if (localProc == 0) {
192  FEI_COUT << "Create-Matrix cpu time: " << fei_creatematrix_time << FEI_ENDL;
193  }
194 
195  double start_load_time = fei::utils::cpu_time();
196 
197  fei::SharedPtr<fei::Vector> solnVec = factory->createVector(matrixGraph, true);
198  fei::SharedPtr<fei::Vector> rhsVec = factory->createVector(matrixGraph);
200  factory->createLinearSystem(matrixGraph);
201 
202  CHK_ERR( linSys->parameters(params));
203 
204  linSys->setMatrix(mat);
205  linSys->setSolutionVector(solnVec);
206  linSys->setRHS(rhsVec);
207 
208  CHK_ERR( HexBeam_Functions::
209  load_elem_data(matrixGraph.get(), mat.get(), rhsVec.get(), hexcube) );
210 
211  //temporarily disable boundary-conditions
212  // CHK_ERR( load_BC_data(linSys, hexcube) );
213 
214  if (!slave_constraints) {
215  CHK_ERR( HexBeam_Functions::
216  load_constraints(linSys.get(), hexcube, firstLocalCRID) );
217  }
218 
219  CHK_ERR( linSys->loadComplete() );
220 
221  double fei_load_time = fei::utils::cpu_time() - start_load_time;
222 
223  if (localProc == 0) {
224  //IOS macros are defined in fei_macros.h
226  FEI_COUT << "Coef. loading cpu time: " << fei_load_time << FEI_ENDL;
227  FEI_COUT << "Total assembly wall time: "
228  << fei_init_time + fei_creatematrix_time + fei_load_time << FEI_ENDL;
229  }
230 
231  //
232  //now the load phase is complete, so we're ready to launch the underlying
233  //solver and solve Ax=b
234  //
235 
236  //First, check whether the params (set from the input .i file) specify
237  //a "Trilinos_Solver" string (which may be Amesos...). If not, it won't
238  //matter but if so, the factory may switch based on this value, when
239  //creating the solver object instance.
240 
241  std::string solver_name_value;
242  params.getStringParamValue("Trilinos_Solver", solver_name_value);
243 
244  const char* charptr_solvername =
245  solver_name_value.empty() ? 0 : solver_name_value.c_str();
246 
247  fei::SharedPtr<fei::Solver> solver = factory->createSolver(charptr_solvername);
248 
249  int status;
250  int itersTaken = 0;
251 
252  if (localProc==0) FEI_COUT << "solve..." << FEI_ENDL;
253  double start_solve_time = fei::utils::cpu_time();
254 
255  int err = solver->solve(linSys.get(),
256  NULL, //preconditioningMatrix
257  params,
258  itersTaken,
259  status);
260 
261  double solve_time = fei::utils::cpu_time()-start_solve_time;
262 
263  if (err!=0) {
264  if (localProc==0) FEI_COUT << "solve returned err: " << err <<", status: "
265  << status << FEI_ENDL;
266  }
267 
268  if (localProc==0) {
269  FEI_COUT << " cpu-time in solve: " << solve_time << FEI_ENDL;
270  }
271 
272  CHK_ERR( solnVec->scatterToOverlap() );
273 
274  delete hexcubeptr;
275 
276  double elapsed_cpu_time = fei::utils::cpu_time() - start_time;
277  int returnValue = 0;
278 
279  //The following IOS_... macros are defined in base/fei_macros.h
281  if (localProc==0) {
282  FEI_COUT << "Proc0 cpu times (seconds):" << FEI_ENDL
283  << " FEI initialize: " << fei_init_time << FEI_ENDL
284  << " FEI create-matrix: " << fei_creatematrix_time << FEI_ENDL
285  << " FEI load: " << fei_load_time << FEI_ENDL
286  << " solve: " << solve_time << FEI_ENDL
287  << "Total program time: " << elapsed_cpu_time << FEI_ENDL;
288 
289 #if defined(FEI_PLATFORM) && defined(FEI_OPT_LEVEL)
290  double benchmark = fei_init_time;
291 
292  std::string slavestr;
293  if (!constraintform.empty()) {
294  slavestr = constraintform;
295  }
296  if (slavestr.size() > 0) slavestr += "_";
297 
298  FEI_OSTRINGSTREAM testname_init;
299  testname_init << "cube3_init_"<<slavestr<<W<<"_"<<D<<"_"<<DofPerNode<<"_"
300  <<solverName<<"_np"<<numProcs<<"_"
301  <<FEI_PLATFORM<<"_"<<FEI_OPT_LEVEL;
302  FEI_OSTRINGSTREAM testname_create;
303  testname_create << "cube3_creatematrix_"<<slavestr<<W<<"_"<<D<<"_"<<DofPerNode
304  <<"_"<<solverName<<"_np"<<numProcs<<"_"
305  <<FEI_PLATFORM<<"_"<<FEI_OPT_LEVEL;
306  FEI_OSTRINGSTREAM testname_load;
307  testname_load << "cube3_load_"<<slavestr<<W<<"_"<<D<<"_"<<DofPerNode<<"_"
308  <<solverName<<"_np"<<numProcs<<"_"
309  <<FEI_PLATFORM<<"_"<<FEI_OPT_LEVEL;
310 
311  double file_init, file_create, file_load;
312  bool file_benchmarks_available = true;
313  try {
314  file_init = fei_test_utils::get_file_benchmark("./cube3_timings.txt",
315  testname_init.str().c_str());
316  file_create = fei_test_utils::get_file_benchmark("./cube3_timings.txt",
317  testname_create.str().c_str());
318  file_load = fei_test_utils::get_file_benchmark("./cube3_timings.txt",
319  testname_load.str().c_str());
320  }
321  catch (std::runtime_error& exc) {
322  file_benchmarks_available = false;
323  }
324 
325  if (file_benchmarks_available) {
326 
327  bool init_test_passed =
328  fei_test_utils::check_and_cout_test_result(testname_init.str(), benchmark,
329  file_init, 10);
330 
331  benchmark = fei_creatematrix_time;
332  bool create_test_passed =
333  fei_test_utils::check_and_cout_test_result(testname_create.str(), benchmark,
334  file_create, 10);
335 
336  benchmark = fei_load_time;
337  bool load_test_passed =
338  fei_test_utils::check_and_cout_test_result(testname_load.str(), benchmark,
339  file_load, 10);
340 
341  returnValue = init_test_passed&&create_test_passed&&load_test_passed ? 0 : 1;
342  }
343 #endif
344 
345  }
346 
347  bool testPassed = returnValue==0;
348  if (testPassed && localProc == 0) {
349  //A string for the SIERRA runtest tool to search for in test output...
350  FEI_COUT << "FEI test successful" << FEI_ENDL;
351  }
352 
353  return(returnValue);
354 }
int load_constraints(FEI *fei, HexBeam &hexcube, int firstLocalCRID)
Definition: HexBeam.cpp:397
virtual void setMatrix(fei::SharedPtr< fei::Matrix > &matrix)
virtual fei::SharedPtr< fei::Solver > createSolver(const char *name=0)=0
virtual void setParameters(const fei::ParameterSet &params)=0
#define FEI_COUT
fei::SharedPtr< fei::Factory > create_fei_Factory(MPI_Comm comm, const char *libraryName)
virtual int parameters(int numParams, const char *const *paramStrings)=0
virtual void parameters(const fei::ParameterSet &paramset)
Definition: fei_Factory.cpp:38
int init_slave_constraints(fei::MatrixGraph *matrixGraph, HexBeam &hexcube)
Definition: HexBeam.cpp:636
virtual void setSolutionVector(fei::SharedPtr< fei::Vector > &soln)
virtual fei::SharedPtr< fei::LinearSystem > createLinearSystem(fei::SharedPtr< fei::MatrixGraph > &matrixGraph)
void defineFields(int numFields, const int *fieldIDs, const int *fieldSizes, const int *fieldTypes=NULL)
int numGlobalDOF_
Definition: HexBeam.hpp:92
virtual void setRHS(fei::SharedPtr< fei::Vector > &rhs)
virtual int scatterToOverlap()=0
int init_constraints(FEI *fei, HexBeam &hexcube, int &firstLocalCRID)
Definition: HexBeam.cpp:355
int totalNumElems_
Definition: HexBeam.hpp:75
virtual int loadComplete(bool applyBCs=true, bool globalAssemble=true)=0
#define MPI_Comm
Definition: fei_mpi.h:56
void setParameters(const fei::ParameterSet &paramset)
virtual fei::SharedPtr< VectorSpace > createVectorSpace(MPI_Comm, const char *name)
int init_elem_connectivities(FEI *fei, HexBeam &hexcube)
Definition: HexBeam.cpp:280
int load_elem_data(FEI *fei, HexBeam &hexcube)
Definition: HexBeam.cpp:440
int getStringParamValue(const char *name, std::string &paramValue) const
#define IOS_FLOATFIELD
double get_file_benchmark(const char *filename, const char *testname)
int get_filename_and_read_input(int argc, char **argv, MPI_Comm comm, int localProc, std::vector< std::string > &stdstrings)
virtual int initComplete()=0
virtual int numDofPerNode()
Definition: HexBeam.hpp:36
int numLocalDOF_
Definition: HexBeam.hpp:91
int init_shared_nodes(FEI *fei, HexBeam &hexcube)
Definition: HexBeam.cpp:328
void defineIDTypes(int numIDTypes, const int *idTypes)
int beam_main(int argc, char **argv, MPI_Comm comm, int numProcs, int localProc)
Definition: beam_main.cpp:43
T * get() const
virtual fei::SharedPtr< fei::Vector > createVector(fei::SharedPtr< fei::VectorSpace > vecSpace, int numVectors=1)=0
#define FEI_ENDL
std::ostream & console_out()
void parse_strings(std::vector< std::string > &stdstrings, const char *separator_string, fei::ParameterSet &paramset)
Definition: fei_utils.cpp:191
double cpu_time()
Definition: fei_utils.cpp:46
int localProc(MPI_Comm comm)
const char * version()
Definition: fei_utils.hpp:53
virtual int getNumCRs()
Definition: HexBeam.hpp:63
virtual int solve(fei::LinearSystem *linearSystem, fei::Matrix *preconditioningMatrix, const fei::ParameterSet &parameterSet, int &iterationsTaken, int &status)
Definition: fei_Solver.cpp:65
virtual fei::SharedPtr< fei::Matrix > createMatrix(fei::SharedPtr< fei::MatrixGraph > matrixGraph)=0
#define CHK_ERR(a)
virtual fei::SharedPtr< fei::MatrixGraph > createMatrixGraph(fei::SharedPtr< fei::VectorSpace > rowSpace, fei::SharedPtr< fei::VectorSpace > columnSpace, const char *name)=0
int localNumElems_
Definition: HexBeam.hpp:77
#define IOS_FIXED
#define FEI_OSTRINGSTREAM
Definition: fei_sstream.hpp:32
int numProcs(MPI_Comm comm)
int getIntParamValue(const char *name, int &paramValue) const
bool check_and_cout_test_result(std::string testname, double value, double file_value, unsigned margin)