extended documentation for scalar example

Author: Daniel Ruprecht

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
+ *
+ *    y' = lambda y = a y + i b y, y(0) = 1.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+i0
+   */
   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,5 +1,10 @@
 /*
  * Sweeper for scalar test equation
+ *
+ * u' = lambda*u, u(0) = u0
+ *
+ * with complex lambda using an IMEX scheme. Derived from the generic imex_sweeper
+ *
  */
 
 #ifndef _EXAMPLES__SCALAR__SCALAR_SWEEPER_HPP_
@@ -14,19 +19,47 @@
 using namespace std;
 
 template<typename time = pfasst::time_precision>
+
+/*
+ * ScalarSweeper is derived from the generic IMEXSweeper
+ */
 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;
+    
+    /*
+     * The complex unit i = sqrt(-1)
+     */
     const 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:
+  
+    /*
+     * Generic constructor; initialize all function call counters with zero.
+     */ 
     ScalarSweeper(complex<double> lambda, complex<double> u0)
       :   lambda(lambda)
         , u0(u0)
@@ -36,6 +69,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 +80,9 @@ class ScalarSweeper
       cout << "Number of implicit solves:      " << this->n_impl_solve << endl;
     }
 
+    /*
+     * Compute error and print it to cout
+     */
     void echo_error(time t)
     {
       auto& qend = pfasst::encap::as_vector<complex<double>, time>(this->get_end_state());
@@ -56,11 +95,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 +114,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 +125,9 @@ class ScalarSweeper
       this->echo_error(t + dt);
     }
 
+    /*
+     * Computes the exact solution u0*exp(lambda*t) at a given time t.
+     */
     void exact(complex_vector_type& q, time t)
     {
       q[0] = this->u0 * exp(this->lambda * t);
@@ -88,6 +139,9 @@ class ScalarSweeper
       this->exact(q, t);
     }
 
+    /*
+     * Evaluate the explicit part of the right hand side: Multiply with imag(lambda)
+     */
     void f_expl_eval(shared_ptr<encap_type> f_encap,
                      shared_ptr<encap_type> q_encap, time t) override
     {
@@ -101,6 +155,9 @@ class ScalarSweeper
       this->n_f_expl_eval++;
     }
 
+    /*
+     * Evaluate the implicit part of the right hand side: Multiply with real(lambda)
+     */
     void f_impl_eval(shared_ptr<encap_type> f_encap,
                      shared_ptr<encap_type> q_encap, time t) override
     {
@@ -114,6 +171,9 @@ class ScalarSweeper
       this->n_f_impl_eval++;
     }
 
+    /*
+     * For given b, solve (Id - dt*real(lambda))*u = b for 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 +186,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++;