Merge branch 'pr/81/head' into release/v0.2.0

* pr/81/head:
  removed const keywords in scalar_sweeper; they should be there, no idea why it worked even with them
  more documentation for scalar example; added mathjax directives at multiple points
  improved docu for scalar example; added ** at several points
  extended documentation for scalar example

Signed-off-by: Torbjörn Klatt <[email protected]>

Author: torbjoernk

examples/scalar/scalar_sdc.cpp

@@ -1,8 +1,13 @@
-/*
+/**
  * 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.
+ *
+ * Solves Dahlquist's test equation with single level SDC
+ *
+ * \\( u' = \\lambda*u \\quad\\text{ , } u(0) = u_0 \\)
+ *
+ * with complex lambda, treating the real part implicitly and the imaginary part 
+ * explicitly.
+ *
  */
 
 #include<complex>
@@ -13,26 +18,32 @@
 
 #include "scalar_sweeper.hpp"
 
-/*
- * 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, 
                       const pfasst::QuadratureType nodetype)
 {
   pfasst::SDC<> sdc;
-
+  
+  /*
+   * For test equation, set initial value to \\( 1+0i )\\
+   */
   const complex<double> y0 = complex<double>(1.0, 0.0);
 
   auto nodes = pfasst::compute_nodes(nnodes, nodetype);
+  
+  /*
+   * This is a scalar example, so we use the encap::VectorFactory with fixed
+   * length of 1 and complex type.
+   */
   auto factory = make_shared<pfasst::encap::VectorFactory<complex<double>>>(1);
   auto sweeper = make_shared<ScalarSweeper<>>(lambda, y0);
 
   sweeper->set_nodes(nodes);
   sweeper->set_factory(factory);
 
   sdc.add_level(sweeper);
+  
+  // Final time Tend = dt*nsteps
   sdc.set_duration(0.0, dt*nsteps, dt, niters);
   sdc.setup();
 
@@ -45,6 +56,9 @@ double run_scalar_sdc(const size_t nsteps, const double dt, const size_t nnodes,
 }
 
 #ifndef PFASST_UNIT_TESTING
+/**
+ * Main routine running the scalar example with a preset parameters
+ */
 int main(int /*argc*/, char** /*argv*/)
 {
   const size_t nsteps = 2;

examples/scalar/scalar_sweeper.hpp

@@ -1,7 +1,3 @@
-/*
- * Sweeper for scalar test equation
- */
-
 #ifndef _EXAMPLES__SCALAR__SCALAR_SWEEPER_HPP_
 #define _EXAMPLES__SCALAR__SCALAR_SWEEPER_HPP_
 
@@ -14,20 +10,45 @@
 using namespace std;
 
 template<typename time = pfasst::time_precision>
+
+/**
+ * Sweeper for scalar test equation
+ *
+ * \\( u' = \\lambda*u \\quad\\text{ , } u(0) = u_0 \\)
+ *
+ * with complex lambda using an IMEX scheme. Derived from the generic imex_sweeper.
+ *
+ */
 class ScalarSweeper
   : public pfasst::encap::IMEXSweeper<time>
 {
+
   private:
     typedef pfasst::encap::Encapsulation<time> encap_type;
+    
+    //! Define a type for a complex PFASST vector encapsulation.
     typedef pfasst::encap::VectorEncapsulation<complex<double>> complex_vector_type;
-
+    
+     //! Parameter lambda and initial value u0 
     complex<double> lambda, u0;
-    size_t n_f_expl_eval, n_f_impl_eval, n_impl_solve;
-    const complex<double> i_complex = complex<double>(0, 1);
+    
+    //! The complex unit i = sqrt(-1)
+    complex<double> i_complex = complex<double>(0, 1);
+    
+     //! Error at the final time. For the scalar example, an analytical solution is known.
     double error;
+    
+     //! Counters for how often f_expl_eval, f_impl_eval and impl_solve are called.
+    size_t n_f_expl_eval, n_f_impl_eval, n_impl_solve;
 
   public:
-    ScalarSweeper(complex<double> lambda, complex<double> u0)
+  
+    /**
+     * Generic constructor; initialize all function call counters with zero.
+     * @param[in] lambda coefficient in test equation
+     * @param[in] u0initial value at \\( t=0 \\)
+     */ 
+    ScalarSweeper(const complex<double> lambda, const complex<double> u0)
       :   lambda(lambda)
         , u0(u0)
         , n_f_expl_eval(0)
@@ -36,6 +57,9 @@ class ScalarSweeper
         , error(0.0)
     {}
 
+    /**
+     * Upon destruction, report final error and number of function calls
+     */
     virtual ~ScalarSweeper()
     {
       cout << "Final error:                    " << scientific << this->error << endl;
@@ -44,6 +68,10 @@ class ScalarSweeper
       cout << "Number of implicit solves:      " << this->n_impl_solve << endl;
     }
 
+    /**
+     * Compute error between last state and exact solution at time tand print it to cout
+     * @param[in] Time t
+     */
     void echo_error(time t)
     {
       auto& qend = pfasst::encap::as_vector<complex<double>, time>(this->get_end_state());
@@ -56,11 +84,17 @@ class ScalarSweeper
       this->error = max_err;
     }
 
+    /**
+     * Returns error, but does not update it!
+     */
     double get_errors()
     {
       return this->error;
     }
 
+    /**
+     * Prediction step and update of error. Uses predictor as provided by IMEXSweeper.
+     */
     void predict(bool initial) override
     {
       pfasst::encap::IMEXSweeper<time>::predict(initial);
@@ -69,6 +103,9 @@ class ScalarSweeper
       this->echo_error(t + dt);
     }
 
+    /**
+     * Perform a sweep and update error. Uses sweep as provided by IMEXSweeper.
+     */
     void sweep() override
     {
       pfasst::encap::IMEXSweeper<time>::sweep();
@@ -77,6 +114,11 @@ class ScalarSweeper
       this->echo_error(t + dt);
     }
 
+    /**
+     * Computes the exact solution \\( u_0 \\exp \\left( \\lambda*t \\right) \\) 
+     * at a given time t.
+     * @param[in] Time t
+     */
     void exact(complex_vector_type& q, time t)
     {
       q[0] = this->u0 * exp(this->lambda * t);
@@ -88,6 +130,10 @@ class ScalarSweeper
       this->exact(q, t);
     }
 
+    /**
+     * Evaluate the explicit part of the right hand side: Multiply with 
+     * \\( \\text{imag}(\\lambda) \\)
+     */
     void f_expl_eval(shared_ptr<encap_type> f_encap,
                      shared_ptr<encap_type> q_encap, time t) override
     {
@@ -101,6 +147,10 @@ class ScalarSweeper
       this->n_f_expl_eval++;
     }
 
+    /**
+     * Evaluate the implicit part of the right hand side: Multiply with 
+     * \\( \\text{real}(\\lambda) \\)
+     */
     void f_impl_eval(shared_ptr<encap_type> f_encap,
                      shared_ptr<encap_type> q_encap, time t) override
     {
@@ -114,6 +164,11 @@ class ScalarSweeper
       this->n_f_impl_eval++;
     }
 
+    /**
+     * For given \\( b \\), solve 
+     * \\( \\left( \\mathbb{I}_d - \\Delta t \\text{real}(\\lambda) \\right) u = b \\) 
+     * for \\( u \\) and set f_encap to \\( \\text{real}(\\lambda) u \\)
+     */
     void impl_solve(shared_ptr<encap_type> f_encap,
                     shared_ptr<encap_type> q_encap, time t, time dt,
                     shared_ptr<encap_type> rhs_encap) override
@@ -126,6 +181,8 @@ class ScalarSweeper
       // invert f_impl = multiply with inverse of real part of lambda
       double inv = 1.0 / (1.0 - double(dt) * real(this->lambda));
       q[0] = inv * rhs[0];
+      
+      // compute f_impl_eval of q[0], i.e. multiply with real(lambda)
       f[0] = real(this->lambda) * q[0];
 
       this->n_impl_solve++;