Merge pull request #142 from memmett/feature/rel_residual

Compute relative residuals as well as absolute residuals.

Author: memmett

examples/advection_diffusion/vanilla_sdc.cpp

@@ -23,7 +23,7 @@ namespace pfasst
   {
     namespace advection_diffusion
     {
-      error_map run_vanilla_sdc(double abs_residual_tol)
+      error_map run_vanilla_sdc(double abs_residual_tol, double rel_residual_tol=0.0)
       {
         SDC<> sdc;
 
@@ -38,7 +38,7 @@ namespace pfasst
 
         sweeper->set_quadrature(quad);
         sweeper->set_factory(factory);
-        sweeper->set_residual_tolerances(abs_residual_tol, 0.0);
+        sweeper->set_residual_tolerances(abs_residual_tol, rel_residual_tol);
 
         sdc.add_level(sweeper);
         sdc.set_duration(0.0, 4*0.01, 0.01, 4);

include/pfasst/encap/encap_sweeper.hpp

@@ -250,21 +250,21 @@ namespace pfasst
         {
           if (this->abs_residual_tol > 0.0 || this->rel_residual_tol > 0.0) {
             if (this->residuals.size() == 0) {
-              for (auto x: this->get_nodes()) {
-                UNUSED(x);
+              for (size_t m = 0; m < this->get_nodes().size(); m++) {
                 this->residuals.push_back(this->get_factory()->create(pfasst::encap::solution));
               }
             }
             this->residual(this->get_controller()->get_time_step(), this->residuals);
-            vector<time> rnorms;
-            for (auto r: this->residuals) {
-              rnorms.push_back(r->norm0());
+            vector<time> anorms, rnorms;
+            for (size_t m = 0; m < this->get_nodes().size(); m++) {
+              anorms.push_back(this->residuals[m]->norm0());
+              rnorms.push_back(anorms.back() / this->get_state(m)->norm0());
             }
+            auto amax = *std::max_element(anorms.begin(), anorms.end());
             auto rmax = *std::max_element(rnorms.begin(), rnorms.end());
-            if (rmax < this->abs_residual_tol) {
+            if (amax < this->abs_residual_tol || rmax < this->rel_residual_tol) {
               return true;
             }
-            // XXX: check rel norms too
           }
           return false;
         }

tests/examples/advection_diffusion/test_advection_diffusion.cpp

@@ -74,6 +74,29 @@ TEST(AdaptiveErrorTest, VanillaSDC)
   ASSERT_EQ(max_iter, (size_t) 2);
 }
 
+TEST(RelativeAdaptiveErrorTest, VanillaSDC)
+{
+  typedef error_map::value_type vtype;
+
+  auto get_iter  = [](const vtype x) { return get<1>(get<0>(x)); };
+  auto get_error = [](const vtype x) { return get<1>(x); };
+
+  auto errors = run_vanilla_sdc(0.0, 1.e-6);
+  auto max_iter = get_iter(*std::max_element(errors.begin(), errors.end(),
+                                             [get_iter](const vtype p1, const vtype p2) { return get_iter(p1) < get_iter(p2); }));
+
+  vector<double> tol = { 5e-8, 5e-8, 5e-8, 5e-8 };
+  vector<double> err;
+  for (auto& x: errors) {
+    if (get_iter(x) == max_iter) {
+      err.push_back(get_error(x));
+    }
+  }
+
+  EXPECT_THAT(err, testing::Pointwise(DoubleLess(), tol));
+  ASSERT_EQ(max_iter, (size_t) 2);
+}
+
 TEST(ErrorTest, SerialMLSDC)
 {
   typedef error_map::value_type vtype;