FEI Package Browser (Single Doxygen Collection)  Version of the Day
poisson.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,
46 // for the purposes of testing, code tuning and scaling studies.
47 //
48 // This program assembles a linear system from a 2D square Poisson
49 // problem, using 4-node square elements. There is only 1 degree-of-
50 // freedom per node, and only one element-block per processor.
51 //
52 // This problem was coded up with Ray Tuminaro.
53 //
54 // The input file for this program should provide the following:
55 // SOLVER_LIBRARY <library-name> -- e.g., Aztec
56 // L <int> -- the global length (num-elements) on a side of the 2D square
57 //
58 //
59 // Alan Williams 03-13-2002
60 //
61 
62 #include "fei_iostream.hpp"
63 #include <cmath>
64 
65 //Including the header fei_base.hpp gets us the declaration for
66 //various classes and functions in the 'fei' namespace.
67 
68 #include "fei_base.hpp"
69 
70 
71 //Now make provision for using any one of several solver libraries. This is
72 //handled by the code in LibraryFactory.{hpp,cpp}.
73 
75 
76 
77 //And, we need to include some headers for utility classes which are simply
78 //used for setting up the data for the example problem.
79 
82 
84 
85 //
86 //Include definitions of macros like 'CHK_ERR' to call functions and check
87 //the return code.
88 //
89 #include "fei_ErrMacros.hpp"
90 
91 
92 //==============================================================================
93 //Here's the main...
94 //==============================================================================
95 int main(int argc, char** argv)
96 {
97  MPI_Comm comm;
98  int numProcs = 1, localProc = 0;
100  comm = MPI_COMM_WORLD;
101 
102  double start_time = fei::utils::cpu_time();
103 
104  //read input parameters from a file specified on the command-line with
105  // '-i file'
106  std::vector<std::string> stdstrings;
108  comm, localProc,
109  stdstrings) );
110 
111  //parse the strings from the input file into a fei::ParameterSet object.
112  fei::ParameterSet paramset;
113  fei::utils::parse_strings(stdstrings, " ", paramset);
114 
115  std::string solverName;
116  int L = 0;
117  int outputLevel = 0;
118 
119  int errcode = 0;
120  errcode += paramset.getStringParamValue("SOLVER_LIBRARY", solverName);
121  errcode += paramset.getIntParamValue("L", L);
122  paramset.getIntParamValue("outputLevel", outputLevel);
123 
124  if (errcode != 0) {
125  fei::console_out() << "Failed to find one or more required parameters in input-file."
126  << FEI_ENDL << "Required parameters:"<<FEI_ENDL
127  << "SOLVER_LIBRARY" << FEI_ENDL
128  << "L" << FEI_ENDL;
129 #ifndef FEI_SER
130  MPI_Finalize();
131 #endif
132  return(-1);
133  }
134 
135  if (localProc == 0) {
136  int nodes = (L+1)*(L+1);
137  int eqns = nodes;
138  //macros FEI_COUT and FEI_ENDL are aliases for std::cout and std::endl,
139  //defined in fei_iostream.hpp.
140  FEI_COUT << "\n========================================================\n";
141  FEI_COUT << "FEI version: " << fei::utils::version() << "\n";
142  FEI_COUT << "Square size L: " << L << " elements.\n";
143  FEI_COUT << "Global number of elements: " << L*L << "\n";
144  FEI_COUT << "Global number of nodes: " << nodes << "\n";
145  FEI_COUT << "Global number of equations: " << eqns <<"\n";
146  FEI_COUT << "========================================================"
147  << FEI_ENDL;
148  }
149 
150  if (outputLevel == 1) {
151  if (localProc != 0) outputLevel = 0;
152  }
153 
154  if (outputLevel>0) {
155  fei_test_utils::print_args(argc, argv);
156  }
157 
158  //PoissonData is the object that will be in charge of generating the
159  //data to pump into the FEI objects.
160 
161  PoissonData poissonData(L, numProcs, localProc, outputLevel);
162 
163  double start_init_time = fei::utils::cpu_time();
164 
166  try {
167  factory = fei::create_fei_Factory(comm, solverName.c_str());
168  }
169  catch (std::runtime_error& exc) {
170  FEI_COUT << "library " << solverName << " not available."<<FEI_ENDL;
171 #ifndef FEI_SER
172  MPI_Finalize();
173 #endif
174  return(-1);
175  }
176 
177  if (factory.get() == NULL) {
178  FEI_COUT << "fei::Factory creation failed." << FEI_ENDL;
179 #ifndef FEI_SER
180  MPI_Finalize();
181 #endif
182  return(-1);
183  }
184 
185  factory->parameters(paramset);
186 
188  factory->createVectorSpace(comm, "poisson3");
189 
192  factory->createMatrixGraph(nodeSpace, dummy, "poisson3");
193 
194  //load some control parameters.
195  matrixGraph->setParameters(paramset);
196 
197 
198  int numFields = poissonData.getNumFields();
199  int* fieldSizes = poissonData.getFieldSizes();
200  int* fieldIDs = poissonData.getFieldIDs();
201  int nodeIDType = 0;
202 
203  if (outputLevel>0 && localProc==0) FEI_COUT << "defineFields" << FEI_ENDL;
204  nodeSpace->defineFields( numFields, fieldIDs, fieldSizes );
205 
206  if (outputLevel>0 && localProc==0) FEI_COUT << "defineIDTypes" << FEI_ENDL;
207  nodeSpace->defineIDTypes( 1, &nodeIDType );
208 
209  CHK_ERR( init_elem_connectivities(matrixGraph.get(), poissonData) );
210 
211  CHK_ERR( set_shared_nodes(nodeSpace.get(), poissonData) );
212 
213 
214  //The following IOS_... macros are defined in base/fei_macros.h
216  if (outputLevel>0 && localProc==0) FEI_COUT << "initComplete" << FEI_ENDL;
217  CHK_ERR( matrixGraph->initComplete() );
218 
219  double fei_init_time = fei::utils::cpu_time() - start_init_time;
220 
221  //Now the initialization phase is complete. Next we'll do the load phase,
222  //which for this problem just consists of loading the element data
223  //(element-wise stiffness arrays and load vectors) and the boundary
224  //condition data.
225  //This simple problem doesn't have any constraint relations, etc.
226 
227  double start_load_time = fei::utils::cpu_time();
228 
229 
230  fei::SharedPtr<fei::Matrix> mat = factory->createMatrix(matrixGraph);
231  fei::SharedPtr<fei::Vector> solnVec = factory->createVector(nodeSpace, true);
232  fei::SharedPtr<fei::Vector> rhsVec = factory->createVector(nodeSpace);
233  fei::SharedPtr<fei::LinearSystem> linSys= factory->createLinearSystem(matrixGraph);
234 
235  linSys->setMatrix(mat);
236  linSys->setSolutionVector(solnVec);
237  linSys->setRHS(rhsVec);
238 
239  CHK_ERR( linSys->parameters(paramset));
240  CHK_ERR( load_elem_data(matrixGraph.get(), mat.get(), rhsVec.get(), poissonData) );
241 
242  CHK_ERR( load_BC_data(linSys.get(), poissonData) );
243 
244  CHK_ERR( linSys->loadComplete() );
245 
246  double fei_load_time = fei::utils::cpu_time() - start_load_time;
247 
248  //
249  //now the load phase is complete, so we're ready to launch the underlying
250  //solver and solve Ax=b
251  //
252 
253  fei::SharedPtr<fei::Solver> solver = factory->createSolver();
254 
255  int status;
256  int itersTaken = 0;
257 
258  if (outputLevel>0 && localProc==0) FEI_COUT << "solve..." << FEI_ENDL;
259  double start_solve_time = fei::utils::cpu_time();
260 
261  int err = solver->solve(linSys.get(),
262  NULL, //preconditioningMatrix
263  paramset, itersTaken, status);
264 
265  double solve_time = fei::utils::cpu_time() - start_solve_time;
266 
267  if (err!=0) {
268  if (localProc==0) FEI_COUT << "solve returned err: " << err <<", status: "
269  << status << FEI_ENDL;
270  }
271 
272  CHK_ERR( solnVec->scatterToOverlap() );
273 
274  //
275  //We should make sure the solution we just computed is correct...
276  //
277 
278  int numNodes = nodeSpace->getNumOwnedAndSharedIDs(nodeIDType);
279 
280  double maxErr = 0.0;
281  if (numNodes > 0) {
282  int lenNodeIDs = numNodes;
283  GlobalID* nodeIDs = new GlobalID[lenNodeIDs];
284  double* soln = new double[lenNodeIDs];
285  if (nodeIDs != NULL && soln != NULL) {
286  CHK_ERR( nodeSpace->getOwnedAndSharedIDs(nodeIDType, numNodes,
287  nodeIDs, lenNodeIDs) );
288 
289  int fieldID = 1;
290  CHK_ERR( solnVec->copyOutFieldData(fieldID, nodeIDType,
291  numNodes, nodeIDs, soln));
292 
293  for(int i=0; i<numNodes; i++) {
294  int nID = (int)nodeIDs[i];
295  double x = (1.0* ((nID-1)%(L+1)))/L;
296  double y = (1.0* ((nID-1)/(L+1)))/L;
297 
298  double exactSoln = x*x + y*y;
299  double error = std::abs(exactSoln - soln[i]);
300  if (maxErr < error) maxErr = error;
301  }
302 
303  delete [] nodeIDs;
304  delete [] soln;
305  }
306  else {
307  fei::console_out() << "allocation of nodeIDs or soln failed." << FEI_ENDL;
308  }
309 
310  }
311 
312 #ifndef FEI_SER
313  double globalMaxErr = 0.0;
314  MPI_Allreduce(&maxErr, &globalMaxErr, 1, MPI_DOUBLE, MPI_MAX, comm);
315  maxErr = globalMaxErr;
316 #endif
317  bool testPassed = true;
318  if (maxErr > 1.e-6) testPassed = false;
319 
320  double elapsed_cpu_time = fei::utils::cpu_time() - start_time;
321  int returnValue = 0;
322 
323  //The following IOS_... macros are defined in base/fei_macros.h
325  if (localProc==0) {
326  FEI_COUT << "Proc0 cpu times (seconds):" << FEI_ENDL
327  << " FEI initialize: " << fei_init_time << FEI_ENDL
328  << " FEI load: " << fei_load_time << FEI_ENDL
329  << " solve: " << solve_time << FEI_ENDL
330  << "Total program time: " << elapsed_cpu_time << FEI_ENDL;
331  }
332 
333  if (testPassed && returnValue==0 && localProc == 0) {
335  FEI_COUT << "poisson: TEST PASSED, maxErr = " << maxErr << ", iterations: "
336  << itersTaken << FEI_ENDL;
337  FEI_COUT << "Poisson test successful" << FEI_ENDL;
338  }
339  if ((testPassed == false || returnValue != 0) && localProc == 0) {
340  FEI_COUT << "maxErr = " << maxErr << ", TEST FAILED\n";
341  FEI_COUT << "(Test is deemed to have passed if the maximum difference"
342  << " between the exact and computed solutions is 1.e-6 or less, *AND*"
343  << " time-taken matches file-benchmark if available.)"
344  << FEI_ENDL;
345  }
346 
347 #ifndef FEI_SER
348  MPI_Finalize();
349 #endif
350  return(returnValue);
351 }
352 
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 IOS_SCIENTIFIC
int getStringParamValue(const char *name, std::string &paramValue) const
#define FEI_COUT
#define MPI_COMM_WORLD
Definition: fei_mpi.h:58
fei::SharedPtr< fei::Factory > create_fei_Factory(MPI_Comm comm, const char *libraryName)
int GlobalID
Definition: fei_defs.h:60
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
virtual int copyOutFieldData(int fieldID, int idType, int numIDs, const int *IDs, double *data, int vectorIndex=0)=0
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)
virtual void setRHS(fei::SharedPtr< fei::Vector > &rhs)
virtual int scatterToOverlap()=0
virtual int loadComplete(bool applyBCs=true, bool globalAssemble=true)=0
int main(int argc, char **argv)
Definition: poisson.cpp:95
#define MPI_Comm
Definition: fei_mpi.h:56
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
int set_shared_nodes(FEI *fei, PoissonData &poissonData)
int getOwnedAndSharedIDs(int idtype, int lenList, int *IDs, int &numOwnedAndSharedIDs)
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
#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 getNumOwnedAndSharedIDs(int idType)
int getIntParamValue(const char *name, int &paramValue) const
int localProc(MPI_Comm comm)
int * getFieldIDs()
Definition: PoissonData.hpp:45
int getNumFields()
Definition: PoissonData.hpp:43
const char * version()
Definition: fei_utils.hpp:53
virtual int solve(fei::LinearSystem *linearSystem, fei::Matrix *preconditioningMatrix, const fei::ParameterSet &parameterSet, int &iterationsTaken, int &status)
Definition: fei_Solver.cpp:65
int * getFieldSizes()
Definition: PoissonData.hpp:44
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
#define IOS_FIXED
int numProcs(MPI_Comm comm)