From ce7308670920e4500ef32b73d6d6da5d62679fc7 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Mon, 11 Jul 2016 23:52:03 +0200
Subject: [PATCH] Tests: Use proper obstacles for ::ProjectedBlockGSStep
---
dune/solvers/test/gssteptest.cc | 51 +++++++++++++++++++++++++++++----
1 file changed, 46 insertions(+), 5 deletions(-)
diff --git a/dune/solvers/test/gssteptest.cc b/dune/solvers/test/gssteptest.cc
index 763e11d0..b9cae7fb 100644
--- a/dune/solvers/test/gssteptest.cc
+++ b/dune/solvers/test/gssteptest.cc
@@ -46,7 +46,7 @@ struct GSTestSuite {
std::vector<BoxConstraint<double, BlockSize>> obstacles(size);
for (size_t i = 0; i < size; ++i)
for (size_t j = 0; j < BlockSize; ++j)
- obstacles[i][j] = { -1, +1 };
+ obstacles[i][j] = { p.u_ex[i][j] + 2.0, p.u_ex[i][j] + 3.0 };
using LoopSolver = ::LoopSolver<Vector>;
using Step = LinearIterationStep<Matrix, Vector>;
@@ -182,6 +182,8 @@ struct GSTestSuite {
// test projected block GS
if (trivialDirichletOnly) // TODO: missing feature in ProjectedBlockGS
{
+ // Test 1: Without obstacles, we should get the same solution as with an
+ // unconstrained solver
size_t size = p.u.size();
using GSStep = ProjectedBlockGSStep<Matrix, Vector>;
typename GSStep::HasObstacle hasObstacle(size, false);
@@ -191,11 +193,50 @@ struct GSTestSuite {
test(&gsStep, "ProjectedBlockGS free");
}
hasObstacle.setAll();
+ // Test2: Different initial iterates should lead to the same solution,
+ // which should, furthermore, not intersect the obstacles
{
- ProjectedBlockGSStep<Matrix, Vector> gsStep;
- gsStep.hasObstacle_ = &hasObstacle;
- gsStep.obstacles_ = &obstacles;
- test(&gsStep, "ProjectedBlockGS obstacle");
+ ProjectedBlockGSStep<Matrix, Vector> pGSStep1;
+ pGSStep1.hasObstacle_ = &hasObstacle;
+ pGSStep1.obstacles_ = &obstacles;
+ ProjectedBlockGSStep<Matrix, Vector> pGSStep2;
+ pGSStep2.hasObstacle_ = &hasObstacle;
+ pGSStep2.obstacles_ = &obstacles;
+
+ Vector u_copy1 = p.u;
+ for (size_t i = 0; i < p.u.size(); ++i)
+ for (size_t j = 0; j < blocksize; ++j)
+ if (p.ignore[i][j])
+ u_copy1[i][j] += 2.25;
+ Vector u_copy2 = p.u;
+ for (size_t i = 0; i < p.u.size(); ++i)
+ for (size_t j = 0; j < blocksize; ++j)
+ if (p.ignore[i][j])
+ u_copy2[i][j] += 2.75;
+
+ auto result1 = solve(&pGSStep1, u_copy1, 1e-12, 2000);
+ auto result2 = solve(&pGSStep2, u_copy2, 1e-12, 2000);
+
+ for (size_t i = 0; i < p.u.size(); ++i)
+ for (size_t j = 0; j < blocksize; ++j)
+ if (std::abs(result1[i][j] -
+ obstacles[i][j].projectIn(result1[i][j])) > 1e-12 or
+ std::abs(result2[i][j] -
+ obstacles[i][j].projectIn(result2[i][j])) > 1e-12) {
+ std::cerr << "### error: obstacles not respected!" << std::endl;
+ passed = false;
+ break;
+ }
+ auto normDiff = p.energyNorm.diff(result1, result2);
+ // We cannot expect too much here (hence the low tolerance):
+ // We are solving a potentially large system using nothing but
+ // a smoother. Convergence will be very slow.
+ if (normDiff > 1e-6) {
+ std::cerr << "### error: Different initial iterates lead to "
+ "different solutions! (diff = "
+ << normDiff << ")" << std::endl;
+ passed = false;
+ }
}
}
--
GitLab