FEI Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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//==============================================================================
95int main(int argc, char** argv)
96{
97 MPI_Comm comm;
98 int numProcs = 1, localProc = 0;
99 CHK_ERR( fei_test_utils::initialize_mpi(argc, argv, localProc, numProcs) );
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
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()
Definition: PoissonData.hpp:45
int getNumFields()
Definition: PoissonData.hpp:43
int * getFieldSizes()
Definition: PoissonData.hpp:44
int getIntParamValue(const char *name, int &paramValue) const
int getStringParamValue(const char *name, std::string &paramValue) const
T * get() const
#define CHK_ERR(a)
int GlobalID
Definition: fei_defs.h:60
#define FEI_ENDL
#define FEI_COUT
#define IOS_FIXED
#define IOS_FLOATFIELD
#define IOS_SCIENTIFIC
#define MPI_COMM_WORLD
Definition: fei_mpi.h:58
#define MPI_Comm
Definition: fei_mpi.h:56
void parse_strings(std::vector< std::string > &stdstrings, const char *separator_string, fei::ParameterSet &paramset)
Definition: fei_utils.cpp:191
const char * version()
Definition: fei_utils.hpp:53
double cpu_time()
Definition: fei_utils.cpp:46
void print_args(int argc, char **argv)
int initialize_mpi(int argc, char **argv, int &localProc, int &numProcs)
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()
int main(int argc, char **argv)
Definition: poisson.cpp:95