00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "SundanceProblemTesting.hpp"
00032
00033 #include "SundanceMaximalCellFilter.hpp"
00034 #include "SundanceDimensionalCellFilter.hpp"
00035 #include "SundancePositionalCellPredicate.hpp"
00036 #include "SundanceMesh.hpp"
00037 #include "SundanceMeshSource.hpp"
00038 #include "SundanceMeshType.hpp"
00039 #include "SundanceBasicSimplicialMeshType.hpp"
00040 #include "SundancePartitionedLineMesher.hpp"
00041 #include "SundancePartitionedRectangleMesher.hpp"
00042 #include "PlayaEpetraVectorType.hpp"
00043 #include "PlayaLinearSolverBuilder.hpp"
00044 #include "SundanceGaussianQuadrature.hpp"
00045 #include "SundanceCoordExpr.hpp"
00046 #include "SundanceCellDiameterExpr.hpp"
00047 #include "SundanceLinearProblem.hpp"
00048 #include "SundanceIntegral.hpp"
00049
00050 namespace Sundance
00051 {
00052
00053 bool checkErrorNorms(
00054 const Mesh& mesh,
00055 const CellFilter& filter,
00056 const Expr& numSoln,
00057 const Expr& exactSoln,
00058 const QuadratureFamily& quad,
00059 double L2Tol,
00060 double H1SemiTol,
00061 double H1Tol)
00062 {
00063 Tabs tab;
00064 Tabs tab1;
00065 Expr df = numSoln - exactSoln;
00066
00067 double L2Err = L2Norm(mesh, filter, df, quad);
00068 Out::root() << tab << "L2 norm check:" << std::endl
00069 << tab1 << setw(16) << setprecision(6) << L2Err
00070 << " tol=" << setw(12) << L2Tol ;
00071 if (fabs(L2Err) > L2Tol) Out::root() << " <==== FAIL!" << std::endl;
00072 else Out::root() << std::endl;
00073
00074 double H1SemiErr = H1Seminorm(mesh, filter, df, quad);
00075 Out::root() << tab << "H1 seminorm check:" << std::endl
00076 << tab1 << setw(16) << setprecision(6) << H1SemiErr
00077 << " tol=" << setw(12) << H1SemiTol ;
00078 if (fabs(H1SemiErr) > H1SemiTol) Out::root() << " <==== FAIL!" << std::endl;
00079 else Out::root() << std::endl;
00080
00081
00082 double H1Err = H1Norm(mesh, filter, df, quad);
00083 Out::root() << tab << "H1 norm check:" << std::endl
00084 << tab1 << setw(16) << setprecision(6) << H1Err
00085 << " tol=" << setw(12) << H1Tol;
00086 if (fabs(H1Err) > H1Tol) Out::root() << " <==== FAIL!" << std::endl;
00087 else Out::root() << std::endl;
00088
00089 return (fabs(L2Err) <= L2Tol)
00090 && (fabs(H1SemiErr) <= H1SemiTol)
00091 && (fabs(H1Err) <= H1Tol) ;
00092 }
00093
00094 double fitPower(const Array<double>& h, const Array<double>& err)
00095 {
00096 Array<double> x(h.size());
00097 Array<double> y(h.size());
00098 double xBar = 0.0;
00099 double yBar = 0.0;
00100 for (int i=0; i<h.size(); i++)
00101 {
00102 x[i] = log(h[i]);
00103 y[i] = log(err[i]);
00104 xBar += x[i];
00105 yBar += y[i];
00106 }
00107
00108 xBar /= h.size();
00109 yBar /= h.size();
00110
00111 double u = 0.0;
00112 double v = 0.0;
00113 for (int i=0; i<h.size(); i++)
00114 {
00115 u += (x[i]-xBar)*(y[i]-yBar);
00116 v += pow(x[i]-xBar,2.0);
00117 }
00118
00119 return u/v;
00120 }
00121
00122
00123
00124 LineDomain::LineDomain(const Array<int>& nx)
00125 : a_(0.0), b_(1.0), nx_(nx), interior_(new MaximalCellFilter()),
00126 left_(), right_(), mesh_()
00127 {init();}
00128
00129
00130 LineDomain::LineDomain(double a, double b, const Array<int>& nx)
00131 : a_(a), b_(b), nx_(nx), interior_(new MaximalCellFilter()),
00132 left_(), right_(), mesh_()
00133 {init();}
00134
00135 void LineDomain::init()
00136 {
00137 int np = MPIComm::world().getNProc();
00138 MeshType meshType = new BasicSimplicialMeshType();
00139 mesh_.resize(nx_.size());
00140 for (int i=0; i<nx_.size(); i++)
00141 {
00142 MeshSource mesher = new PartitionedLineMesher(a_, b_, np*nx_[i], meshType);
00143 mesh_[i] = mesher.getMesh();
00144 }
00145
00146 CellFilter points = new DimensionalCellFilter(0);
00147 left_ = points.subset(new CoordinateValueCellPredicate(0,a_));
00148 right_ = points.subset(new CoordinateValueCellPredicate(0,b_));
00149 }
00150
00151
00152
00153 RectangleDomain::RectangleDomain(const Array<int>& n)
00154 : ax_(0.0), bx_(1.0), nx_(n),
00155 ay_(0.0), by_(1.0), ny_(n),
00156 interior_(new MaximalCellFilter()),
00157 north_(), south_(), east_(), west_(),
00158 mesh_()
00159 {init();}
00160
00161 RectangleDomain::RectangleDomain(
00162 double ax, double bx, const Array<int>& nx,
00163 double ay, double by, const Array<int>& ny
00164 )
00165 : ax_(ax), bx_(bx), nx_(nx),
00166 ay_(ay), by_(by), ny_(ny),
00167 interior_(new MaximalCellFilter()),
00168 north_(), south_(), east_(), west_(),
00169 mesh_()
00170 {init();}
00171
00172
00173 void RectangleDomain::init()
00174 {
00175 TEUCHOS_TEST_FOR_EXCEPT(nx_.size() != ny_.size());
00176 int np = MPIComm::world().getNProc();
00177 int npx = -1;
00178 int npy = -1;
00179
00180 PartitionedRectangleMesher::balanceXY(np, &npx, &npy);
00181 TEUCHOS_TEST_FOR_EXCEPT(npx < 1);
00182 TEUCHOS_TEST_FOR_EXCEPT(npy < 1);
00183 TEUCHOS_TEST_FOR_EXCEPT(npx * npy != np);
00184
00185 MeshType meshType = new BasicSimplicialMeshType();
00186 mesh_.resize(nx_.size());
00187
00188 for (int i=0; i<nx_.size(); i++)
00189 {
00190 MeshSource mesher
00191 = new PartitionedRectangleMesher(
00192 ax_, bx_, nx_[i], npx,
00193 ay_, by_, ny_[i], npy,
00194 meshType);
00195 mesh_[i] = mesher.getMesh();
00196 }
00197
00198 CellFilter edges = new DimensionalCellFilter(1);
00199 north_ = edges.subset(new CoordinateValueCellPredicate(1,by_));
00200 south_ = edges.subset(new CoordinateValueCellPredicate(1,ay_));
00201 west_ = edges.subset(new CoordinateValueCellPredicate(0,ax_));
00202 east_ = edges.subset(new CoordinateValueCellPredicate(0,bx_));
00203 }
00204
00205
00206
00207
00208
00209
00210
00211 Array<double> L2NormCalculator::computeNorms(
00212 const ForwardProblemTestBase* prob,
00213 int meshIndex,
00214 const Expr& numSoln, const Expr& exactSoln) const
00215 {
00216 Expr errFunc = (numSoln - exactSoln).flatten();
00217
00218 Mesh mesh = prob->getMesh(meshIndex);
00219 CellFilter interior = prob->interior();
00220
00221 Array<int> p = prob->pExpected();
00222 TEUCHOS_TEST_FOR_EXCEPTION(p.size() != errFunc.size(),
00223 std::runtime_error,
00224 "size mismatch between array of expected orders (p=" << p << ") and "
00225 "array of solutions: " << errFunc.size());
00226
00227 Array<double> rtn(p.size());
00228 for (int i=0; i<errFunc.size(); i++)
00229 {
00230 QuadratureFamily quad = new GaussianQuadrature(2*p[i]);
00231 double L2Err = L2Norm(mesh, interior, errFunc[i], quad);
00232 rtn[i] = L2Err;
00233 }
00234
00235 return rtn;
00236 }
00237
00238
00239 Array<LPTestSpec> LPTestBase::specs() const
00240 {
00241 return tuple(
00242 LPTestSpec("amesos.xml", 1.0e-10, makeSet<int>(1)),
00243 LPTestSpec("aztec-ifpack.xml", 1.0e-10),
00244 LPTestSpec("aztec-ml.xml", 1.0e-10),
00245 LPTestSpec("belos-ifpack.xml", 1.0e-8),
00246 LPTestSpec("belos-ml.xml", 1.0e-10)
00247 );
00248 }
00249
00250
00251 std::ostream& operator<<(std::ostream& os, const LPTestSpec& spec)
00252 {
00253 os << "LPTestSpec(tol=" << spec.tol() << ", solver=" << spec.solverFile()
00254 << std::endl;
00255 return os;
00256 }
00257
00258
00259
00260
00261 bool LPTestSuite::run() const
00262 {
00263 int np = MPIComm::world().getNProc();
00264
00265 bool allOK = true;
00266
00267 for (int i=0; i<tests_.size(); i++)
00268 {
00269 Tabs tab(0);
00270 Array<LPTestSpec> specs = tests_[i]->specs();
00271 for (int j=0; j<specs.size(); j++)
00272 {
00273 if (specs[j].nProcIsAllowed(np))
00274 {
00275 Out::root() << std::endl;
00276 Out::root() << std::endl;
00277 Out::root() << std::endl;
00278 Out::root() << tab
00279 << "-------------------------------------"
00280 "-------------------------------------"
00281 << std::endl;
00282 Out::root() << tab << "running test " << tests_[i]->name()
00283 << " with spec " << specs[j] << std::endl;
00284 Out::root() << tab
00285 << "-------------------------------------"
00286 "-------------------------------------"
00287 << std::endl;
00288
00289 std::string solverFile = specs[j].solverFile();
00290 double tol = specs[j].tol();
00291 bool pass = tests_[i]->run(solverFile, tol);
00292 allOK = pass && allOK;
00293 }
00294 else
00295 {
00296 Out::root() << tab << "skipping test " << tests_[i]->name()
00297 << " with spec=" << specs[j] << std::endl;
00298 }
00299 }
00300 Out::root() << std::endl;
00301 Out::root() << std::endl;
00302 }
00303 return allOK;
00304 }
00305
00306
00307 LPTestSuite::LPTestSuite()
00308 : tests_(), testSpecs_() {}
00309
00310
00311 void LPTestSuite::registerTest(const RCP<LPTestBase>& test)
00312 {
00313 tests_.append(test);
00314 }
00315
00316 VectorType<double> ForwardProblemTestBase::vecType() const
00317 {
00318 return new EpetraVectorType();
00319 }
00320
00321 Expr ForwardProblemTestBase::coord(int d) const
00322 {
00323 TEUCHOS_TEST_FOR_EXCEPT(d<0 || d>2);
00324 return new CoordExpr(d);
00325 }
00326
00327 double ForwardProblemTestBase::cellSize(int i) const
00328 {
00329 Expr h = new CellDiameterExpr();
00330 Expr hExpr = Integral(interior(), h, new GaussianQuadrature(1));
00331 Expr AExpr = Integral(interior(), 1.0, new GaussianQuadrature(1));
00332
00333 double area = evaluateIntegral(getMesh(i), AExpr);
00334 double hMean = evaluateIntegral(getMesh(i), hExpr)/area;
00335
00336 return hMean;
00337 }
00338
00339 RCP<ErrNormCalculatorBase> ForwardProblemTestBase::normCalculator() const
00340 {
00341 return rcp(new L2NormCalculator());
00342
00343 }
00344 bool ForwardProblemTestBase::run(const std::string& solverFile,
00345 double tol) const
00346 {
00347 if (numMeshes()==1) return runSingleTest(solverFile, tol);
00348 else return runTestSequence(solverFile, tol);
00349 }
00350
00351
00352 bool ForwardProblemTestBase::runTestSequence(const std::string& solverFile,
00353 const double& pTol) const
00354 {
00355 Tabs tab(0);
00356 Out::root() << tab << "Running test sequence " << name() << std::endl;
00357
00358 LinearSolver<double> solver
00359 = LinearSolverBuilder::createSolver(solverFile);
00360
00361 Expr exact = exactSoln();
00362
00363 Array<double> hList(numMeshes());
00364 Array<Array<double> > errList(numMeshes());
00365
00366 for (int i=0; i<numMeshes(); i++)
00367 {
00368 Tabs tab1;
00369 Out::root() << tab1 << "running on mesh #" << i << std::endl;
00370 Expr soln;
00371 bool solveOK = solve(getMesh(i), solver, soln);
00372 if (!solveOK) return false;
00373
00374 postRunCallback(i, getMesh(i), solverFile, soln);
00375
00376 double h = cellSize(i);
00377 Array<double> err = normCalculator()->computeNorms(
00378 this, i, soln, exact);
00379 hList[i] = h;
00380 errList[i] = err;
00381 }
00382
00383 bool allPass = true;
00384 Out::root() << tab << "Error results: " << std::endl;
00385 for (int f=0; f < pExpected().size(); f++)
00386 {
00387 Tabs tab1;
00388 Out::root() << tab1 << "Variable #" << f << std::endl;
00389 Out::root() << std::endl << tab1 << "Error norm versus h" << std::endl;
00390 Array<double> err;
00391 for (int i=0; i<numMeshes(); i++)
00392 {
00393 Tabs tab2;
00394 Out::root() << tab2 << setw(20) << setprecision(5) << hList[i]
00395 << " " << setw(20) << setprecision(5) << errList[i][f]
00396 << std::endl;
00397 err.append(errList[i][f]);
00398 }
00399
00400 double p = fitPower(hList, err);
00401 Out::root() << tab << "Measured exponent: " << p << std::endl;
00402 Out::root() << tab << "Expected exponent: " << pExpected()[f] << std::endl;
00403 double pErr = fabs(p - pExpected()[f]) ;
00404 Out::root() << tab << "Difference: " << setw(12) << pErr
00405 << " Tolerance: " << setw(12) << pTol;
00406 bool pass = pErr <= pTol;
00407 if (!pass) Out::root() << " <==== FAIL!";
00408 Out::root() << std::endl;
00409 allPass = allPass && pass;
00410 }
00411
00412 return allPass;
00413 }
00414
00415
00416
00417 bool ForwardProblemTestBase::runSingleTest(const std::string& solverFile,
00418 const double& tol) const
00419 {
00420 Tabs tab(0);
00421
00422 LinearSolver<double> solver
00423 = LinearSolverBuilder::createSolver(solverFile);
00424
00425 Expr exact = exactSoln();
00426
00427 Expr soln;
00428 bool solveOK = solve(getMesh(0), solver, soln);
00429 if (!solveOK) return false;
00430
00431 Array<double> err = normCalculator()->computeNorms(this, 0,
00432 soln, exact);
00433
00434 for (int i=0; i<err.size(); i++)
00435 {
00436 Out::root() << tab << setw(20) << setprecision(5) << err[i]
00437 << " " << setw(20) << setprecision(5) << tol
00438 << std::endl;
00439 }
00440 return err[0] <= tol;
00441 }
00442
00443 bool LPTestBase::solve(const Mesh& mesh,
00444 const LinearSolver<double>& solver,
00445 Expr& soln) const
00446 {
00447 Tabs tab(0);
00448 LinearProblem::solveFailureIsFatal() = false;
00449 SolverState<double> state = prob(mesh).solve(solver, soln);
00450
00451 bool converged = (state.finalState() == SolveConverged) ;
00452 if (!converged)
00453 {
00454 Out::root() << tab << "solve failed to converge!" << std::endl;
00455 Tabs tab1;
00456 Out::root() << tab1 << "state = " << state << std::endl;
00457 }
00458 return converged;
00459 }
00460
00461 LP1DTestBase::LP1DTestBase(const Array<int>& nx)
00462 : domain_(nx) {}
00463
00464 LP1DTestBase::LP1DTestBase(double a, double b, const Array<int>& nx)
00465 : domain_(a, b, nx) {}
00466
00467 LPRectTestBase::LPRectTestBase(const Array<int>& n)
00468 : domain_(n) {}
00469 }
