FEI Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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>
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
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
45#include <test_utils/CRSet.hpp>
48
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//============================================================================
63int 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;
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;
172 CHK_ERR( fei->setSolveType(FEI_SINGLE_SYSTEM) );
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
int init_elem_connectivities(FEI *fei, PoissonData &poissonData)
int load_elem_data(FEI *fei, PoissonData &poissonData)
int set_shared_nodes(FEI *fei, PoissonData &poissonData)
int load_BC_data(FEI *fei, PoissonData &poissonData)
int * getFieldIDs()
int getNumFields()
int * getFieldSizes()
int getIntParamValue(const char *name, int &paramValue) const
int getStringParamValue(const char *name, std::string &paramValue) const
void reset(T *p=0)
#define ERReturn(a)
#define CHK_ERR(a)
int GlobalID
Definition fei_defs.h:60
#define FEI_SINGLE_SYSTEM
Definition fei_defs.h:65
#define FEI_ENDL
#define FEI_COUT
#define IOS_FIXED
#define IOS_FLOATFIELD
#define IOS_SCIENTIFIC
#define MPI_Comm
Definition fei_mpi.h:56
#define FEI_Implementation
Definition fei_version.h:66
void parse_strings(std::vector< std::string > &stdstrings, const char *separator_string, fei::ParameterSet &paramset)
double cpu_time()
Definition fei_utils.cpp:46
void strings_to_char_ptrs(std::vector< std::string > &stdstrings, int &numStrings, const char **&charPtrs)
void print_args(int argc, char **argv)
int whichArg(int argc, const char *const *argv, const char *findarg)
int get_filename_and_read_input(int argc, char **argv, MPI_Comm comm, int localProc, std::vector< std::string > &stdstrings)
fei::SharedPtr< fei::Factory > create_fei_Factory(MPI_Comm comm, const char *libraryName)
std::ostream & console_out()
fei::SharedPtr< LibraryWrapper > create_LibraryWrapper(MPI_Comm comm, const char *libraryName)
int poisson_main(int argc, char **argv, MPI_Comm comm, int numProcs, int localProc)