FEI Version of the Day
Loading...
Searching...
No Matches
poisson3_main.cpp
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// 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, and only one element-block per processor.
16//
17// This problem was coded up with Ray Tuminaro.
18//
19// The input file for this program should provide the following:
20// SOLVER_LIBRARY <library-name> -- e.g., Aztec
21// L <int> -- the global length (num-elements) on a side of the 2D square
22//
23// Alan Williams 03-13-2002
24//
25
26#include <fei_iostream.hpp>
27#include <cmath>
28
29//Including the header fei_base.hpp gets us the declaration for
30//the various FEI classes.
31
32#include <fei_base.hpp>
33
34//Now make provision for using any one of several solver libraries. This is
35//handled by the code in LibraryFactory.[hC].
36
37#include <test_utils/LibraryFactory.hpp>
38
39//And, we need to include some headers for utility classes which are simply
40//used for setting up the data for the example problem.
41
42#include <test_utils/Poisson_Elem.hpp>
43#include <test_utils/PoissonData.hpp>
44
45#include <test_utils/ElemBlock.hpp>
46#include <test_utils/CRSet.hpp>
47#include <test_utils/CommNodeSet.hpp>
48#include <test_utils/DataReader.hpp>
49
50#include <test_utils/SolnCheck.hpp>
51#include <test_utils/fei_test_utils.hpp>
52#include <snl_fei_Utils.hpp>
53//
54//Include definitions of macros to call functions and check the return code.
55//
56#include <fei_ErrMacros.hpp>
57
58
59//==============================================================================
60//Here's the main...
61//==============================================================================
62int poisson3_main(int argc, char** argv,
63 MPI_Comm comm, int numProcs, int localProc){
64 int masterProc = 0;
65
66 double start_time = fei::utils::cpu_time();
67
68 std::vector<std::string> stdstrings;
70 comm, localProc,
71 stdstrings) );
72 const char** params = NULL;
73 int numParams = 0;
74 fei::utils::strings_to_char_ptrs(stdstrings, numParams, params);
75
76 fei::ParameterSet paramset;
77 fei::utils::parse_strings(stdstrings, " ", paramset);
78
79 std::string solverName;
80 int L = 0;
81 int outputLevel = 0;
82
83 int errcode = 0;
84 errcode += paramset.getStringParamValue("SOLVER_LIBRARY", solverName);
85 errcode += paramset.getIntParamValue("L", L);
86 paramset.getIntParamValue("outputLevel", outputLevel);
87
88 if (errcode != 0) {
89 fei::console_out() << "Failed to find one or more required parameters in input-file."
90 << FEI_ENDL << "Required parameters:"<<FEI_ENDL
91 << "SOLVER_LIBRARY" << FEI_ENDL
92 << "L" << FEI_ENDL;
93 return(-1);
94 }
95
96 if (localProc == 0) {
97 int nodes = (L+1)*(L+1);
98 int eqns = nodes;
99 FEI_COUT << FEI_ENDL;
100 FEI_COUT << "========================================================"
101 << FEI_ENDL;
102 FEI_COUT << "FEI version: " << fei::utils::version() << FEI_ENDL;
103 FEI_COUT << "Square size L: " << L << " elements." << FEI_ENDL;
104 FEI_COUT << "Global number of elements: " << L*L << FEI_ENDL;
105 FEI_COUT << "Global number of nodes: " << nodes << FEI_ENDL;
106 FEI_COUT << "Global number of equations: " << eqns <<FEI_ENDL;
107 FEI_COUT << "========================================================"
108 << FEI_ENDL;
109 }
110
111 if ((masterProc == localProc)&&(outputLevel>0)) {
112 fei_test_utils::print_args(argc, argv);
113 }
114
115 if (outputLevel == 1) {
116 if (localProc != 0) outputLevel = 0;
117 }
118
119 //PoissonData is the object that will be in charge of generating the
120 //data to pump into the FEI objects.
121
122 PoissonData poissonData(L, numProcs, localProc, outputLevel);
123
124 double start_init_time = fei::utils::cpu_time();
125
127 try {
128 factory = fei::create_fei_Factory(comm, solverName.c_str());
129 }
130 catch (...) {
131 FEI_COUT << "library " << solverName << " not available."<<FEI_ENDL;
132 return(-1);
133 }
134
135 if (factory.get() == NULL) {
136 FEI_COUT << "snl_fei::Factory creation failed." << FEI_ENDL;
137 return(-1);
138 }
139
140 factory->parameters(paramset);
141
143 factory->createVectorSpace(comm, "poisson3");
144
147 factory->createMatrixGraph(nodeSpace, dummy, "poisson3");
148
149 //load some control parameters.
150 matrixGraph->setParameters(paramset);
151
152
153 int numFields = poissonData.getNumFields();
154 int* fieldSizes = poissonData.getFieldSizes();
155 int* fieldIDs = poissonData.getFieldIDs();
156 int nodeIDType = 0;
157
158 if (outputLevel>0 && localProc==0) FEI_COUT << "defineFields" << FEI_ENDL;
159 nodeSpace->defineFields( numFields, fieldIDs, fieldSizes );
160
161 if (outputLevel>0 && localProc==0) FEI_COUT << "defineIDTypes" << FEI_ENDL;
162 nodeSpace->defineIDTypes( 1, &nodeIDType );
163
164 CHK_ERR( init_elem_connectivities(matrixGraph.get(), poissonData) );
165
166 CHK_ERR( set_shared_nodes(nodeSpace.get(), poissonData) );
167
168
169 //The following IOS_... macros are defined in base/fei_macros.h
170 FEI_COUT.setf(IOS_FIXED, IOS_FLOATFIELD);
171 if (outputLevel>0 && localProc==0) FEI_COUT << "initComplete" << FEI_ENDL;
172 CHK_ERR( matrixGraph->initComplete() );
173
174 double fei_init_time = fei::utils::cpu_time() - start_init_time;
175
176 //Now the initialization phase is complete. Next we'll do the load phase,
177 //which for this problem just consists of loading the element data
178 //(element-wise stiffness arrays and load vectors) and the boundary
179 //condition data.
180 //This simple problem doesn't have any constraint relations, etc.
181
182 double start_load_time = fei::utils::cpu_time();
183
184
185 fei::SharedPtr<fei::Matrix> mat = factory->createMatrix(matrixGraph);
186 fei::SharedPtr<fei::Vector> solnVec = factory->createVector(nodeSpace, true);
187 fei::SharedPtr<fei::Vector> rhsVec = factory->createVector(nodeSpace);
188 fei::SharedPtr<fei::LinearSystem> linSys= factory->createLinearSystem(matrixGraph);
189
190 linSys->setMatrix(mat);
191 linSys->setSolutionVector(solnVec);
192 linSys->setRHS(rhsVec);
193
194 CHK_ERR( linSys->parameters(paramset));
195 CHK_ERR( load_elem_data(matrixGraph.get(), mat.get(), rhsVec.get(), poissonData) );
196
197 CHK_ERR( load_BC_data(linSys.get(), poissonData) );
198
199 CHK_ERR( linSys->loadComplete() );
200
201 double fei_load_time = fei::utils::cpu_time() - start_load_time;
202
203 //
204 //now the load phase is complete, so we're ready to launch the underlying
205 //solver and solve Ax=b
206 //
207
208 fei::SharedPtr<fei::Solver> solver = factory->createSolver();
209
210 int status;
211 int itersTaken = 0;
212
213 if (outputLevel>0 && localProc==0) FEI_COUT << "solve..." << FEI_ENDL;
214 double start_solve_time = fei::utils::cpu_time();
215
216 int err = solver->solve(linSys.get(),
217 NULL, //preconditioningMatrix
218 paramset, itersTaken, status);
219
220 double solve_time = fei::utils::cpu_time() - start_solve_time;
221
222 if (err!=0) {
223 if (localProc==0) FEI_COUT << "solve returned err: " << err <<", status: "
224 << status << FEI_ENDL;
225 }
226
227 CHK_ERR( solnVec->scatterToOverlap() );
228
229 //
230 //We oughtta make sure the solution we just computed is correct...
231 //
232
233 int numNodes = nodeSpace->getNumOwnedAndSharedIDs(nodeIDType);
234
235 double maxErr = 0.0;
236 if (numNodes > 0) {
237 int lenNodeIDs = numNodes;
238 GlobalID* nodeIDs = new GlobalID[lenNodeIDs];
239 double* soln = new double[lenNodeIDs];
240 if (nodeIDs != NULL && soln != NULL) {
241 CHK_ERR( nodeSpace->getOwnedAndSharedIDs(nodeIDType, numNodes,
242 nodeIDs, lenNodeIDs) );
243
244 int fieldID = 1;
245 CHK_ERR( solnVec->copyOutFieldData(fieldID, nodeIDType,
246 numNodes, nodeIDs, soln));
247
248 for(int i=0; i<numNodes; i++) {
249 int nID = (int)nodeIDs[i];
250 double x = (1.0* ((nID-1)%(L+1)))/L;
251 double y = (1.0* ((nID-1)/(L+1)))/L;
252
253 double exactSoln = x*x + y*y;
254 double error = std::abs(exactSoln - soln[i]);
255 if (maxErr < error) maxErr = error;
256 }
257
258 delete [] nodeIDs;
259 delete [] soln;
260 }
261 else {
262 fei::console_out() << "allocation of nodeIDs or soln failed." << FEI_ENDL;
263 }
264
265 }
266
267#ifndef FEI_SER
268 double globalMaxErr = 0.0;
269 MPI_Allreduce(&maxErr, &globalMaxErr, 1, MPI_DOUBLE, MPI_MAX, comm);
270 maxErr = globalMaxErr;
271#endif
272 bool testPassed = true;
273 if (maxErr > 1.e-6) testPassed = false;
274
275 double elapsed_cpu_time = fei::utils::cpu_time() - start_time;
276 int returnValue = 0;
277
278 //The following IOS_... macros are defined in base/fei_macros.h
279 FEI_COUT.setf(IOS_FIXED, IOS_FLOATFIELD);
280 if (localProc==0) {
281 FEI_COUT << "Proc0 cpu times (seconds):" << FEI_ENDL
282 << " FEI initialize: " << fei_init_time << FEI_ENDL
283 << " FEI load: " << fei_load_time << FEI_ENDL
284 << " solve: " << solve_time << FEI_ENDL
285 << "Total program time: " << elapsed_cpu_time << FEI_ENDL;
286
287#if defined(FEI_PLATFORM) && defined(FEI_OPT_LEVEL)
288 double benchmark = fei_init_time+fei_load_time;
289 FEI_OSTRINGSTREAM testname;
290 testname << "poisson3_"<<L<<"_"<<solverName<<"_np"<<numProcs<<"_"
291 <<FEI_PLATFORM<<"_"<<FEI_OPT_LEVEL;
292
293 double file_value = 0.0;
294 bool file_value_available = true;
295 try {
296 file_value = fei_test_utils::get_file_benchmark("./poisson3_timings.txt",
297 testname.str().c_str());
298 }
299 catch(std::runtime_error& exc) {
300 file_value_available = false;
301 }
302
303 if (file_value_available) {
304 bool test_passed = fei_test_utils::check_and_cout_test_result(testname.str(),
305 benchmark,
306 file_value, 10);
307 returnValue = test_passed ? 0 : 1;
308 }
309#endif
310
311 }
312
313 if (testPassed && returnValue==0 && localProc == 0) {
314 FEI_COUT.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
315 FEI_COUT << "poisson: TEST PASSED, maxErr = " << maxErr << ", iterations: "
316 << itersTaken << FEI_ENDL;
317 //This is something the SIERRA runtest tool looks for in test output...
318 FEI_COUT << "SIERRA execution successful" << FEI_ENDL;
319 }
320 if ((testPassed == false || returnValue != 0) && localProc == 0) {
321 if (testPassed==true) {
322 FEI_COUT << "maxErr = " << maxErr << ", but time-taken outside margin. TEST FAILED" << FEI_ENDL;
323 }
324 else {
325 FEI_COUT << "maxErr = " << maxErr << ", TEST FAILED"<<FEI_ENDL;
326 }
327 FEI_COUT << "(Test is deemed to have passed if the maximum difference"
328 << " between the exact and computed solutions is 1.e-6 or less, *AND*"
329 << " time-taken matches file-benchmark if available.)"
330 << FEI_ENDL;
331 }
332
333 return(returnValue);
334}
335
int getIntParamValue(const char *name, int &paramValue) const
int getStringParamValue(const char *name, std::string &paramValue) const
void parse_strings(std::vector< std::string > &stdstrings, const char *separator_string, fei::ParameterSet &paramset)
const char * version()
Definition fei_utils.hpp:53
double cpu_time()
Definition fei_utils.cpp:46
void strings_to_char_ptrs(std::vector< std::string > &stdstrings, int &numStrings, const char **&charPtrs)
int get_filename_and_read_input(int argc, char **argv, MPI_Comm comm, int localProc, std::vector< std::string > &stdstrings)
double get_file_benchmark(const char *filename, const char *testname)
fei::SharedPtr< fei::Factory > create_fei_Factory(MPI_Comm comm, const char *libraryName)
std::ostream & console_out()
int numProcs(MPI_Comm comm)