FEI Package Browser (Single Doxygen Collection)  Version of the Day
beam.cpp
Go to the documentation of this file.
1 /*
2 // @HEADER
3 // ************************************************************************
4 // FEI: Finite Element Interface to Linear Solvers
5 // Copyright (2005) Sandia Corporation.
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
8 // U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Alan Williams (william@sandia.gov)
38 //
39 // ************************************************************************
40 // @HEADER
41 */
42 
43 
44 //
45 // This is a simple program to exercise FEI classes for the
46 // purposes of testing, code tuning and scaling studies.
47 //
48 
49 #include "fei_iostream.hpp"
50 
51 //Including the header fei_base.hpp pulls in the declaration for
52 //various classes in the 'fei' namespace.
53 
54 #include "fei_base.hpp"
55 
56 
57 //Make provision for using any one of several solver libraries. This is
58 //handled by the code in LibraryFactory.{hpp,cpp}.
59 
61 
62 
63 //we need to include some 'test_utils' headers which are used for
64 //setting up the data for the hex-beam example problem.
65 
66 #include "test_utils/HexBeam.hpp"
67 #include "test_utils/HexBeamCR.hpp"
68 
69 
70 //fei_test_utils.hpp declares things like the intialize_mpi function.
71 
73 
74 //
75 //Include definitions of macros like 'CHK_ERR' to call functions and check
76 //the return code.
77 //
78 #undef fei_file
79 #define fei_file "beam.cpp"
80 #include "fei_ErrMacros.hpp"
81 
82 //==============================================================================
83 //Here's the main...
84 //==============================================================================
85 int main(int argc, char** argv)
86 {
87  MPI_Comm comm;
88  int numProcs, localProc;
90  comm = MPI_COMM_WORLD;
91 
92  double start_time = fei::utils::cpu_time();
93 
94  //read input parameters from a file specified on the command-line with
95  // '-i file'
96  std::vector<std::string> stdstrings;
98  comm, localProc,
99  stdstrings) );
100  fei::ParameterSet params;
101  //parse the strings from the input file into a fei::ParameterSet object.
102  fei::utils::parse_strings(stdstrings, " ", params);
103 
104  std::string solverName;
105  std::string datasource;
106  std::string constraintform;
107  int W = 0;
108  int D = 0;
109  int DofPerNode = 0;
110 
111  int errcode = 0;
112  errcode += params.getStringParamValue("SOLVER_LIBRARY", solverName);
113  errcode += params.getStringParamValue("DATA_SOURCE", datasource);
114  errcode += params.getIntParamValue("W", W);
115  errcode += params.getIntParamValue("D", D);
116  errcode += params.getIntParamValue("DofPerNode", DofPerNode);
117 
118  if (errcode != 0) {
119  fei::console_out() << "Failed to find one or more required parameters in input-file."
120  << std::endl << "Required parameters:"<<std::endl
121  << "SOLVER_LIBRARY" << std::endl
122  << "DATA_SOURCE" << std::endl
123  << "CONSTRAINT_FORM" << std::endl
124  << "W" << std::endl << "D" << std::endl << "DofPerNode" << std::endl;
125 
126 #ifndef FEI_SER
127  MPI_Finalize();
128 #endif
129 
130  return(-1);
131  }
132 
133  params.getStringParamValue("CONSTRAINT_FORM", constraintform);
134 
135  bool slave_constraints = (constraintform == "slaves");
136 
137  HexBeam* hexcubeptr = NULL;
138  if (datasource == "HexBeam") {
139  hexcubeptr = new HexBeam(W, D, DofPerNode,
141  }
142  else {
143  hexcubeptr = new HexBeamCR(W, D, DofPerNode,
145  }
146 
147  HexBeam& hexcube = *hexcubeptr;
148 
149  if (localProc == 0) {
150  int numCRs = (W+1)*(W+1)* ((numProcs*2)-1);
151  if (hexcube.getNumCRs() < 1) numCRs = 0;
152  //macros std::cout and std::endl are aliases for std::cout and std::endl,
153  //defined in fei_iostream.hpp.
154  std::cout << std::endl;
155  std::cout << "========================================================\n";
156  std::cout << "FEI version: " << fei::utils::version() << std::endl;
157  std::cout << "--------------------------------------------------------\n";
158  std::cout << "Size W: " << W << " (num-elements-along-side-of-cube)\n";
159  std::cout << "Size D: " << D << " (num-elements-along-depth-of-cube)\n";
160  std::cout << "DOF per node: " << DofPerNode <<"\n";
161  std::cout << "Num local elements: " << hexcube.localNumElems_ << "\n";
162  std::cout << "Num global elements: " << hexcube.totalNumElems_ << "\n";
163  std::cout << "Num local DOF: " << hexcube.numLocalDOF_ << "\n";
164  std::cout << "Num global DOF: " << hexcube.numGlobalDOF_ << "\n";
165  std::cout << "Num global CRs: " << numCRs << "\n";
166  std::cout << "========================================================"
167  << std::endl;
168  }
169 
170  double start_init_time = fei::utils::cpu_time();
171 
173  try {
174  factory = fei::create_fei_Factory(comm, solverName.c_str());
175  }
176  catch (...) {
177  std::cout << "library " << solverName << " not available."<<std::endl;
178 
179 #ifndef FEI_SER
180  MPI_Finalize();
181 #endif
182 
183  return(-1);
184  }
185  if (factory.get() == NULL) {
186  std::cout << "fei::Factory allocation failed." << std::endl;
187 
188 #ifndef FEI_SER
189  MPI_Finalize();
190 #endif
191 
192  return(-1);
193  }
194 
195  factory->parameters(params);
196 
198  factory->createVectorSpace(comm, "beam");
199 
202  factory->createMatrixGraph(nodeSpace, dummy, "beam");
203 
204  //load some solver-control parameters.
205  nodeSpace->setParameters(params);
206  matrixGraph->setParameters(params);
207 
208  int fieldID = 0;
209  int fieldSize = hexcube.numDofPerNode();
210  int nodeIDType = 0;
211  int crIDType = 1;
212 
213  nodeSpace->defineFields( 1, &fieldID, &fieldSize );
214  nodeSpace->defineIDTypes(1, &nodeIDType );
215  nodeSpace->defineIDTypes(1, &crIDType );
216 
217  //HexBeam_Functions::init_elem_connectivities demonstrates the process
218  //of initializing element-to-node connectivities on the fei::MatrixGraph.
220  init_elem_connectivities( matrixGraph.get(), hexcube ) );
221 
222  //HexBeam_Functions::init_shared_nodes demonstrates the process of
223  //initializing shared nodes on the fei::MatrixGraph.
225  init_shared_nodes( matrixGraph.get(), hexcube ) );
226 
227  int firstLocalCRID = 0;
228  if (slave_constraints) {
230  init_slave_constraints( matrixGraph.get(), hexcube) );
231  }
232  else {
234  init_constraints( matrixGraph.get(), hexcube, localProc, firstLocalCRID) );
235  }
236 
237  CHK_ERR( matrixGraph->initComplete() );
238 
239  double fei_init_time = fei::utils::cpu_time() - start_init_time;
240 
241  if (localProc == 0) {
242  std::cout.setf(IOS_FIXED, IOS_FLOATFIELD);
243  std::cout << "Initialization cpu time: " << fei_init_time << std::endl;
244  }
245 
246  //Now the initialization phase is complete. Next we'll do the load phase,
247  //which for this problem just consists of loading the element data
248  //(element-wise stiffness arrays and load vectors) and the boundary
249  //condition data.
250 
251  double fei_creatematrix_start_time = fei::utils::cpu_time();
252 
253  fei::SharedPtr<fei::Matrix> mat = factory->createMatrix(matrixGraph);
254 
255  double fei_creatematrix_time = fei::utils::cpu_time() - fei_creatematrix_start_time;
256  if (localProc == 0) {
257  std::cout.setf(IOS_FIXED, IOS_FLOATFIELD);
258  std::cout << "Create-Matrix cpu time: " << fei_creatematrix_time << std::endl;
259  }
260 
261  double start_load_time = fei::utils::cpu_time();
262 
263  fei::SharedPtr<fei::Vector> solnVec = factory->createVector(matrixGraph, true);
264  fei::SharedPtr<fei::Vector> rhsVec = factory->createVector(matrixGraph);
266  factory->createLinearSystem(matrixGraph);
267 
268  CHK_ERR( linSys->parameters(params));
269 
270  linSys->setMatrix(mat);
271  linSys->setSolutionVector(solnVec);
272  linSys->setRHS(rhsVec);
273 
275  load_elem_data(matrixGraph.get(), mat.get(), rhsVec.get(), hexcube) );
276 
277  //temporarily disable boundary-conditions
278  // CHK_ERR( load_BC_data(linSys, hexcube) );
279 
280  if (!slave_constraints) {
282  load_constraints(linSys.get(), hexcube, firstLocalCRID) );
283  }
284 
285  CHK_ERR( linSys->loadComplete() );
286 
287  double fei_load_time = fei::utils::cpu_time() - start_load_time;
288 
289  if (localProc == 0) {
290  //IOS macros are defined in fei_stdinc.h
291  std::cout.setf(IOS_FIXED, IOS_FLOATFIELD);
292  std::cout << "Coef. loading cpu time: " << fei_load_time << std::endl;
293  std::cout << "Total assembly wall time: "
294  << fei_init_time + fei_creatematrix_time + fei_load_time << std::endl;
295  }
296 
297  //
298  //now the load phase is complete, so we're ready to launch the underlying
299  //solver and solve Ax=b
300  //
301 
302  //First, check whether the params (set from the input .i file) specify
303  //a "Trilinos_Solver" string (which may be Amesos...). If not, it won't
304  //matter but if so, the factory may switch based on this value, when
305  //creating the solver object instance.
306 
307  std::string solver_name_value;
308  params.getStringParamValue("Trilinos_Solver", solver_name_value);
309 
310  const char* charptr_solvername =
311  solver_name_value.empty() ? NULL : solver_name_value.c_str();
312 
313  fei::SharedPtr<fei::Solver> solver = factory->createSolver(charptr_solvername);
314 
315  int status;
316  int itersTaken = 0;
317 
318  if (localProc==0) std::cout << "solve..." << std::endl;
319  double start_solve_time = fei::utils::cpu_time();
320 
321  int err = solver->solve(linSys.get(),
322  NULL, //preconditioningMatrix
323  params,
324  itersTaken,
325  status);
326 
327  double solve_time = fei::utils::cpu_time()-start_solve_time;
328 
329  if (err!=0 && localProc==0) {
330  std::cout << "solve returned err: " << err <<", status: "
331  << status << std::endl;
332  }
333 
334  if (localProc==0) {
335  std::cout << " cpu-time in solve: " << solve_time << std::endl;
336  }
337 
338  CHK_ERR( solnVec->scatterToOverlap() );
339 
340  delete hexcubeptr;
341 
342  double elapsed_cpu_time = fei::utils::cpu_time() - start_time;
343 
344  //The following IOS_... macros are defined in base/fei_iostream.hpp
345  std::cout.setf(IOS_FIXED, IOS_FLOATFIELD);
346  if (localProc==0) {
347  std::cout << "Proc0 cpu times (seconds):" << std::endl
348  << " FEI initialize: " << fei_init_time << std::endl
349  << " FEI create-matrix: " << fei_creatematrix_time << std::endl
350  << " FEI load: " << fei_load_time << std::endl
351  << " solve: " << solve_time << std::endl
352  << "Total program time: " << elapsed_cpu_time << std::endl;
353 
354  }
355 
356 
357 #ifndef FEI_SER
358  MPI_Finalize();
359 #endif
360 
361  return(0);
362 }
363 
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
int getStringParamValue(const char *name, std::string &paramValue) const
#define MPI_COMM_WORLD
Definition: fei_mpi.h:58
fei::SharedPtr< fei::Factory > create_fei_Factory(MPI_Comm comm, const char *libraryName)
T * get() const
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
#define IOS_FLOATFIELD
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 initialize_mpi(int argc, char **argv, int &localProc, int &numProcs)
virtual fei::SharedPtr< fei::Vector > createVector(fei::SharedPtr< fei::VectorSpace > vecSpace, int numVectors=1)=0
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 getIntParamValue(const char *name, int &paramValue) const
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
int main(int argc, char **argv)
Definition: beam.cpp:85
int numProcs(MPI_Comm comm)