FEI Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
poisson_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 the FEI, for the purposes of
11 // testing, code tuning and scaling studies.
12 //
13 // This program assembles a linear system from a 2D square Poisson
14 // problem, using 4-node square elements. There is only 1 degree-of-
15 // freedom per node.
16 //
17 // This problem was coded up with the help of Ray Tuminaro.
18 // Alan Williams 12-20-2000
19 //
20 // The input file for this program should provide the following:
21 // WHICH_FEI <OLDFEI|fei::FEI_Impl>
22 // SOLVER_LIBRARY <library-name> -- e.g., Aztec
23 // L <int> -- the global length (num-elements) on a side of the 2D square
24 //
25 //
26 
27 #include <fei_iostream.hpp>
28 #include <cmath>
29 #include <fei_base.hpp>
30 #include <FEI_Implementation.hpp>
31 
32 //Now make provision for using any one of several solver libraries. This is
33 //handled by the code in LibraryFactory.{hpp,cpp}.
34 
36 #include <fei_LibraryWrapper.hpp>
37 
38 //And, we need to include some headers for utility classes which are simply
39 //used for setting up the data for the example problem.
40 
43 
44 #include <test_utils/ElemBlock.hpp>
45 #include <test_utils/CRSet.hpp>
48 
49 #include <test_utils/SolnCheck.hpp>
50 
52 #include <snl_fei_Utils.hpp>
53 //
54 //Include definitions of macros to call functions and check the return code.
55 //
56 #undef fei_file
57 #define fei_file "poisson_main.cpp"
58 #include <fei_ErrMacros.hpp>
59 
60 //============================================================================
61 //Here's the main...
62 //============================================================================
63 int poisson_main(int argc, char** argv,
64  MPI_Comm comm, int numProcs, int localProc){
65  int outputLevel;
66  int masterProc = 0, err = 0;
67 
68  double start_time = fei::utils::cpu_time();
69 
70  std::vector<std::string> stdstrings;
72  comm, localProc,
73  stdstrings) );
74  const char** params = NULL;
75  int numParams = 0;
76  fei::utils::strings_to_char_ptrs(stdstrings, numParams, params);
77 
78  fei::ParameterSet paramset;
79  fei::utils::parse_strings(stdstrings, " ", paramset);
80 
81  std::string which_fei;
82  std::string solverName;
83  int L = 0;
84 
85  int errcode = 0;
86  errcode += paramset.getStringParamValue("SOLVER_LIBRARY", solverName);
87  errcode += paramset.getStringParamValue("WHICH_FEI", which_fei);
88  errcode += paramset.getIntParamValue("L", L);
89 
90  if (errcode != 0) {
91  fei::console_out() << "Failed to find one or more required parameters in input-file."
92  << FEI_ENDL << "Required parameters:"<<FEI_ENDL
93  << "SOLVER_LIBRARY" << FEI_ENDL
94  << "WHICH_FEI" << FEI_ENDL
95  << "L" << FEI_ENDL;
96  return(-1);
97  }
98 
99  if (localProc == 0) {
100  int nodes = (L+1)*(L+1);
101  int eqns = nodes;
102  FEI_COUT << FEI_ENDL;
103  FEI_COUT << "========================================================"
104  << FEI_ENDL;
105  FEI_COUT << "Square size L: " << L << " elements." << FEI_ENDL;
106  FEI_COUT << "Global number of elements: " << L*L << FEI_ENDL;
107  FEI_COUT << "Global number of nodes: " << nodes << FEI_ENDL;
108  FEI_COUT << "Global number of equations: " << eqns <<FEI_ENDL;
109  FEI_COUT << "========================================================"
110  << FEI_ENDL;
111  }
112 
113  outputLevel = fei_test_utils::whichArg(numParams, params, "outputLevel 1");
114  if (outputLevel >= 0) outputLevel = 1;
115  if (outputLevel < 0) outputLevel = 0;
116 
117  if ((masterProc == localProc)&&(outputLevel>0)) {
118  fei_test_utils::print_args(argc, argv);
119  }
120 
121  if (outputLevel == 1) {
122  if (localProc != 0) outputLevel = 0;
123  }
124 
125  //PoissonData is the object that will be in charge of generating the
126  //data to pump into the FEI object.
127 
128  PoissonData poissonData(L, numProcs, localProc, outputLevel);
129 
130  double start_init_time = fei::utils::cpu_time();
131 
135 
136  if (which_fei == "OLDFEI") {
137  try {
138  wrapper = fei::create_LibraryWrapper(comm, solverName.c_str());
139  }
140  catch (std::runtime_error& exc) {
141  fei::console_out() << exc.what()<<FEI_ENDL;
142  ERReturn(-1);
143  }
144  fei.reset(new FEI_Implementation(wrapper, comm));
145  }
146  else if (which_fei == "fei::FEI_Impl") {
147  try {
148  factory = fei::create_fei_Factory(comm, solverName.c_str());
149  }
150  catch (std::runtime_error& exc) {
151  fei::console_out() << exc.what()<<FEI_ENDL;
152  ERReturn(-1);
153  }
154  fei = factory->createFEI(comm);
155  }
156  else {
157  fei::console_out() << "poisson_main ERROR, value of 'WHICH_FEI' must be 'OLDFEI' or 'fei::FEI_Impl'"<< FEI_ENDL;
158  ERReturn(-1);
159  }
160 
161  const char* feiVersionString;
162  CHK_ERR( fei->version(feiVersionString) );
163 
164  if (localProc==0) FEI_COUT << feiVersionString << FEI_ENDL;
165 
166  //load some parameters.
167  CHK_ERR( fei->parameters( numParams, params ) );
168 
169  delete [] params;
170 
171  if (outputLevel>0 && localProc==0) FEI_COUT << "setSolveType" << FEI_ENDL;
173 
174 
175  int numFields = poissonData.getNumFields();
176  int* fieldSizes = poissonData.getFieldSizes();
177  int* fieldIDs = poissonData.getFieldIDs();
178 
179  if (outputLevel>0 && localProc==0) FEI_COUT << "initFields" << FEI_ENDL;
180  CHK_ERR( fei->initFields( numFields, fieldSizes, fieldIDs ) );
181 
182 
183  CHK_ERR( init_elem_connectivities(fei.get(), poissonData) );
184 
185  CHK_ERR( set_shared_nodes(fei.get(), poissonData) );
186 
187  //The FEI_COUT and IOS_... macros are defined in base/fei_iostream.hpp
189 
190  if (outputLevel>0 && localProc==0) FEI_COUT << "initComplete" << FEI_ENDL;
191  CHK_ERR( fei->initComplete() );
192 
193  double fei_init_time = fei::utils::cpu_time() - start_init_time;
194 
195  //Now the initialization phase is complete. Next we'll do the load phase,
196  //which for this problem just consists of loading the element data
197  //(element-wise stiffness arrays and load vectors) and the boundary
198  //condition data.
199  //This simple problem doesn't have an constraint relations, etc.
200 
201  double start_load_time = fei::utils::cpu_time();
202 
203  CHK_ERR( load_elem_data(fei.get(), poissonData) );
204 
205  CHK_ERR( load_BC_data(fei.get(), poissonData) );
206 
207  double fei_load_time = fei::utils::cpu_time() - start_load_time;
208 
209  //
210  //now the load phase is complete, so we're ready to launch the underlying
211  //solver and solve Ax=b
212  //
213  int status;
214  if (outputLevel>0 && localProc==0) FEI_COUT << "solve..." << FEI_ENDL;
215  double start_solve_time = fei::utils::cpu_time();
216  err = fei->solve(status);
217  if (err) {
218  if (localProc==0) FEI_COUT << "solve returned err: " << err << FEI_ENDL;
219  }
220 
221  double iTime, lTime, sTime, rTime;
222  CHK_ERR(fei->cumulative_cpu_times(iTime, lTime, sTime, rTime) );
223 
224  double solve_time = fei::utils::cpu_time() - start_solve_time;
225 
226  if (localProc == 0) {
227  FEI_COUT << "FEI cpu-times:" << FEI_ENDL
228  << " init. phase: " << iTime << FEI_ENDL
229  << " load phase: " << lTime << FEI_ENDL
230  << " solve time: " << sTime << FEI_ENDL;
231  }
232 
233  double norm = 0.0;
235  CHK_ERR( fei->residualNorm(1, 1, &fieldIDs[0], &norm) );
236  if (localProc == 0) {
237  FEI_COUT << "returned residual norm: " << norm << FEI_ENDL;
238  }
239 
240  int itersTaken = 0;
241  CHK_ERR( fei->iterations(itersTaken) );
242 
243  //
244  //We oughtta make sure the solution we just computed is correct...
245  //
246 
247  int numNodes = 0;
248  CHK_ERR( fei->getNumLocalNodes(numNodes) );
249 
250  double maxErr = 0.0;
251  if (numNodes > 0) {
252  int lenNodeIDs = numNodes;
253  GlobalID* nodeIDs = new GlobalID[lenNodeIDs];
254  double* soln = new double[lenNodeIDs];
255  if (nodeIDs != NULL && soln != NULL) {
256  CHK_ERR( fei->getLocalNodeIDList(numNodes, nodeIDs, lenNodeIDs) );
257 
258  int fieldID = 1;
259  CHK_ERR( fei->getNodalFieldSolution(fieldID, numNodes, nodeIDs, soln));
260 
261  for(int i=0; i<numNodes; i++) {
262  int nID = (int)nodeIDs[i];
263  double x = (1.0* ((nID-1)%(L+1)))/L;
264  double y = (1.0* ((nID-1)/(L+1)))/L;
265 
266  double exactSoln = x*x + y*y;
267  double error = std::abs(exactSoln - soln[i]);
268  if (maxErr < error) maxErr = error;
269  }
270 
271  delete [] nodeIDs;
272  delete [] soln;
273  }
274  else {
275  fei::console_out() << "allocation of nodeIDs or soln failed." << FEI_ENDL;
276  }
277 
278  }
279 
280 #ifndef FEI_SER
281  double globalMaxErr = 0.0;
282  MPI_Allreduce(&maxErr, &globalMaxErr, 1, MPI_DOUBLE, MPI_MAX, comm);
283  maxErr = globalMaxErr;
284 #endif
285  bool testPassed = true;
286  if (maxErr > 1.e-6) testPassed = false;
287 
288  if (testPassed && localProc == 0) {
289  FEI_COUT << "poisson: TEST PASSED, maxErr = " << maxErr << ", iterations: "
290  << itersTaken << FEI_ENDL;
291  //This is something the SIERRA runtest tool looks for in test output...
292  FEI_COUT << "SIERRA execution successful" << FEI_ENDL;
293  }
294  if (testPassed == false && localProc == 0) {
295  FEI_COUT << "maxErr = " << maxErr << ", TEST FAILED" << FEI_ENDL;
296  FEI_COUT << "(Test is deemed to have passed if the maximum difference"
297  << " between the exact and computed solutions is 1.e-6 or less.)"
298  << FEI_ENDL;
299  }
300 
301  double elapsed_cpu_time = fei::utils::cpu_time() - start_time;
302 
303  //The following IOS_... macros are defined in base/fei_macros.h
305  if (localProc==0) {
306  FEI_COUT << "Proc0 cpu times (seconds):" << FEI_ENDL
307  << " FEI initialize: " << fei_init_time << FEI_ENDL
308  << " FEI load: " << fei_load_time << FEI_ENDL
309  << " solve: " << solve_time << FEI_ENDL
310  << "Total program time: " << elapsed_cpu_time << FEI_ENDL;
311  }
312 
313  wrapper.reset();
314  fei.reset();
315  factory.reset();
316  //If Prometheus is being used, we need to make sure that the
317  //LibraryWrapper inside factory is deleted before MPI_Finalize() is called.
318  //(That's why we call the 'reset' methods on these shared-pointers rather
319  //than just letting them destroy themselves when they go out of scope.)
320 
321  return(0);
322 }
323 
void strings_to_char_ptrs(std::vector< std::string > &stdstrings, int &numStrings, const char **&charPtrs)
Definition: fei_utils.cpp:178
#define IOS_SCIENTIFIC
#define FEI_COUT
virtual int parameters(int numParams, const char *const *paramStrings)=0
virtual int initComplete()=0
fei::SharedPtr< fei::Factory > create_fei_Factory(MPI_Comm comm, const char *libraryName)
int GlobalID
Definition: fei_defs.h:60
virtual int solve(int &status)=0
int poisson_main(int argc, char **argv, MPI_Comm comm, int numProcs, int localProc)
virtual int getNodalFieldSolution(int fieldID, int numNodes, const GlobalID *nodeIDs, double *results)=0
#define FEI_Implementation
Definition: fei_version.h:66
int whichArg(int argc, const char *const *argv, const char *findarg)
void reset(T *p=0)
virtual int getNumLocalNodes(int &numNodes)=0
#define MPI_Comm
Definition: fei_mpi.h:56
virtual int version(const char *&versionString)=0
int init_elem_connectivities(FEI *fei, HexBeam &hexcube)
Definition: HexBeam.cpp:280
#define FEI_SINGLE_SYSTEM
Definition: fei_defs.h:65
int load_elem_data(FEI *fei, HexBeam &hexcube)
Definition: HexBeam.cpp:440
int getStringParamValue(const char *name, std::string &paramValue) const
#define IOS_FLOATFIELD
int get_filename_and_read_input(int argc, char **argv, MPI_Comm comm, int localProc, std::vector< std::string > &stdstrings)
int set_shared_nodes(FEI *fei, PoissonData &poissonData)
virtual int initFields(int numFields, const int *fieldSizes, const int *fieldIDs, const int *fieldTypes=NULL)=0
#define ERReturn(a)
T * get() const
virtual int setSolveType(int solveType)=0
std::string error
#define FEI_ENDL
void print_args(int argc, char **argv)
std::ostream & console_out()
void parse_strings(std::vector< std::string > &stdstrings, const char *separator_string, fei::ParameterSet &paramset)
Definition: fei_utils.cpp:191
int load_BC_data(FEI *fei, HexBeam &hexcube)
Definition: HexBeam.cpp:484
double cpu_time()
Definition: fei_utils.cpp:46
int localProc(MPI_Comm comm)
int * getFieldIDs()
Definition: PoissonData.hpp:45
int getNumFields()
Definition: PoissonData.hpp:43
virtual int residualNorm(int whichNorm, int numFields, int *fieldIDs, double *norms)=0
virtual int iterations(int &itersTaken) const =0
int * getFieldSizes()
Definition: PoissonData.hpp:44
virtual fei::SharedPtr< FEI > createFEI(fei::SharedPtr< LibraryWrapper > wrapper, MPI_Comm comm)
Definition: fei_Factory.cpp:65
#define CHK_ERR(a)
virtual int getLocalNodeIDList(int &numNodes, GlobalID *nodeIDs, int lenNodeIDs)=0
#define IOS_FIXED
fei::SharedPtr< LibraryWrapper > create_LibraryWrapper(MPI_Comm comm, const char *libraryName)
virtual int cumulative_cpu_times(double &initPhase, double &loadPhase, double &solve, double &solnReturn)=0
int numProcs(MPI_Comm comm)
int getIntParamValue(const char *name, int &paramValue) const