Merge pull request #64 from danielru/feature/scalar_tests

Very nice!  High fives!

Author: memmett

examples/scalar/scalar_sdc.cpp

@@ -1,3 +1,10 @@
+/*
+ * Scalar test equation example using an encapsulated IMEX sweeper.
+ * Solves Dahlquist's test equation
+ * y' = a y + i b y
+ * treating the real part implicitly and the imaginary part explicitly.
+ */
+ 
 #include<complex>
 
 #include<pfasst.hpp>
@@ -6,32 +13,45 @@
 
 #include "scalar_sweeper.hpp"
 
-int main(int, char**)
+/*
+ * A simple run routine with preset parameters
+ */
+
+double run_scalar_sdc(const size_t nsteps, const double dt, const size_t nnodes,
+  const size_t niters, const complex<double> lambda)
 {
+
   pfasst::SDC<> sdc;
 
-  const size_t nsteps = 2;
-  const double dt     = 1.0;
-  const size_t nnodes = 4;
-  const size_t niters = 6;
-  const complex<double> lambda = complex<double>(-1.0, 1.0);
   const complex<double> y0     = complex<double>(1.0, 0.0);
-
-  // so far, only lobatto nodes appear to be working
-  auto nodes   = pfasst::compute_nodes(nnodes, "gauss-lobatto");
+  auto nodes = pfasst::compute_nodes(nnodes, "gauss-lobatto");
   auto factory = make_shared<pfasst::encap::VectorFactory<complex<double>>>(1);
-  auto sweeper = make_shared<ScalarSweeper<>>(lambda, y0);
-
+  auto sweeper = make_shared<ScalarSweeper<>>(lambda,y0);
+  
   sweeper->set_nodes(nodes);
   sweeper->set_factory(factory);
-
+  
   sdc.add_level(sweeper);
-  sdc.set_duration(0.0, dt * nsteps, dt, niters);
+  sdc.set_duration(0.0, dt*nsteps, dt, niters);
   sdc.setup();
-
+  
   auto q0 = sweeper->get_state(0);
   sweeper->exact(q0, 0.0);
-
+  
   sdc.run();
+  
+  return sweeper->get_errors();
+}
 
+#ifndef PFASST_UNIT_TESTING
+int main(int /*argc*/, char** /*argv*/)
+{
+  const size_t nsteps = 2;
+  const double dt     = 1.0;
+  const size_t nnodes = 4;
+  const size_t niters = 6;
+  const complex<double> lambda = complex<double>(-1.0, 1.0);
+  
+  run_scalar_sdc(nsteps, dt, nnodes, niters, lambda);
 }
+#endif

examples/scalar/scalar_sweeper.hpp

@@ -32,6 +32,7 @@ class ScalarSweeper : public pfasst::encap::IMEXSweeper<time>
 
     virtual ~ScalarSweeper()
     {
+      cout << "Final error:               " << scientific << this->_error << endl;
       cout << "Number of calls to f1eval: " << this->_nf1eval << endl;
       cout << "Number of calls to f2eval: " << this->_nf2eval << endl;
       cout << "Number of calls to f2comp: " << this->_nf2comp << endl;
@@ -46,6 +47,12 @@ class ScalarSweeper : public pfasst::encap::IMEXSweeper<time>
       this->exact(qex, t);
       double max_err = abs(qend[0] - qex[0]) / abs(qex[0]);
       cout << "err: " << scientific << max_err << endl;
+      this->_error = max_err;
+    }
+    
+    double get_errors()
+    {
+      return this->_error;    
     }
 
     void predict(bool initial) override
@@ -121,6 +128,7 @@ class ScalarSweeper : public pfasst::encap::IMEXSweeper<time>
     complex<double> _lambda, _y0;
     int _nf1eval, _nf2eval, _nf2comp;
     const complex<double> i_complex = complex<double>(0, 1);
+    double _error;
 
 };
 #endif

tests/examples/CMakeLists.txt

@@ -1,2 +1,3 @@
 
 add_subdirectory(advection_diffusion)
+add_subdirectory(scalar)

tests/examples/scalar/CMakeLists.txt

@@ -0,0 +1,25 @@
+# Building and Running Tests
+include_directories(
+    ${3rdparty_INCLUDES}
+    ${pfasst_INCLUDES}
+)
+
+set(TESTS
+    test_scalar
+)
+
+foreach(test ${TESTS})
+    message(STATUS "  ${test}")
+    add_executable(${test} ${test}.cpp)
+    add_dependencies(${test} googlemock)
+    target_link_libraries(${test}
+        ${3rdparty_DEPENDEND_LIBS}
+        ${pfasst_DEPENDED_LIBS}
+    )
+    set_target_properties(${test}
+        PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/tests/examples/scalar
+    )
+    add_test(NAME ${test}
+        COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/tests/examples/scalar/${test} --gtest_output=xml:${test}_out.xml
+    )
+endforeach(test)

tests/examples/scalar/test_scalar.cpp

@@ -0,0 +1,71 @@
+/*
+ * Tests for the scalar example solving the test equation
+ */ 
+
+#include <cmath>
+#include <gtest/gtest.h>
+#include <gmock/gmock.h>
+
+#define PFASST_UNIT_TESTING
+#include "../examples/scalar/scalar_sdc.cpp"
+#undef PFASST_UNIT_TESTING
+
+MATCHER (DoubleMore, "")
+{
+  return get<0>(arg) > get<1>(arg);
+}
+
+
+/*
+ * For Lobatto nodes, the resulting method should of order 2*M-2 with M=number of nodes
+ * The test below verifies that the code approximately (up to a safety factor) reproduces
+ * the theoretically expected rate of convergence 
+ */
+TEST(ConvergenceTest, ScalarSDC)
+{
+    
+  const complex<double> lambda = complex<double>(-1.0, 1.0);
+  const double Tend = 4.0;
+  const vector<size_t> nsteps = { 2, 5, 10, 15, 20 };
+  size_t nsteps_l = nsteps.size();
+
+  vector<double> err(nsteps.size());
+  vector<double> convrate(nsteps.size()-1);
+    
+  double dt;
+    
+  // Test converge up to 6 nodes: For more nodes, errors are of the order of machine 
+  // precision and convergence can no longer be monitored  
+  for ( size_t nnodes = 2; nnodes<=6; ++nnodes)
+  {
+    // Expect convergence rate of 2*nodes-2 from collocation formula,
+    // doing an identical number of iteration should suffice to reach this as each
+    // iteration should increase order by one
+    size_t niters = 2*nnodes-2;
+  
+    for ( size_t i = 0; i<=nsteps_l-1; ++i)
+    {
+      dt = Tend/double(nsteps[i]);
+      err[i] = run_scalar_sdc(nsteps[i], dt, nnodes, niters, lambda);
+    }
+  
+    for ( size_t i = 0; i<=nsteps_l-2; ++i)
+    {
+      convrate[i] = log10(err[i+1]/err[i])/log10(double(nsteps[i])/double(nsteps[i+1]));
+      
+      EXPECT_THAT(convrate[i], testing::DoubleNear(double(2*nnodes-2), 0.99)) << "Convergence rate at node " << i << " not within expected range";
+    }
+  
+    // NOTE: There is probably a much more elegant way to test this, because
+    // expected_cr contains the same value in all entries. But I could not so far figure
+    // out how to build a more clever MATCHER here so far....
+    // EXPECT_THAT(convrate, testing::Pointwise(DoubleMore(), expected_cr ));
+  
+  }
+}
+
+int main(int argc, char** argv)
+{
+  testing::InitGoogleTest(&argc, argv);
+  return RUN_ALL_TESTS();
+}