Merge branch 'development' into feature/scalar2

Author: Matthew Emmett

examples/advection_diffusion/advection_diffusion_sweeper.hpp

@@ -5,12 +5,15 @@
 #ifndef _ADVECTION_DIFFUSION_SWEEPER_HPP_
 #define _ADVECTION_DIFFUSION_SWEEPER_HPP_
 
+#include <cstdlib>
 #include <cassert>
 #include <complex>
 #include <map>
 #include <ostream>
 #include <vector>
+#include <utility>
 
+#include <pfasst/globals.hpp>
 #include <pfasst/encap/imex_sweeper.hpp>
 
 #include "fft.hpp"
@@ -24,42 +27,56 @@ using pfasst::encap::Encapsulation;
 using pfasst::encap::as_vector;
 
 
+/**
+ * errors at different iterations and time nodes
+ * 
+ * Mapping a pair of step/iteration indices onto the error of the solution.
+ */
 typedef map<pair<size_t, size_t>, double> error_map;
 
 template<typename time = pfasst::time_precision>
 class AdvectionDiffusionSweeper
   : public pfasst::encap::IMEXSweeper<time>
 {
+    //! @{
     FFT fft;
-
     vector<complex<double>> ddx, lap;
+    //! @}
+
+    //! @{
     error_map errors;
+    //! @}
 
+    //! @{
     double v  = 1.0;
     time   t0 = 1.0;
     double nu = 0.02;
     size_t nf1evals = 0;
+    //! @}
 
   public:
+    //! @{
     AdvectionDiffusionSweeper(size_t nvars)
     {
-      ddx.resize(nvars);
-      lap.resize(nvars);
+      this->ddx.resize(nvars);
+      this->lap.resize(nvars);
       for (size_t i = 0; i < nvars; i++) {
         double kx = 2 * PI * ((i <= nvars / 2) ? int(i) : int(i) - int(nvars));
-        ddx[i] = complex<double>(0.0, 1.0) * kx;
-        lap[i] = (kx * kx < 1e-13) ? 0.0 : -kx * kx;
+        this->ddx[i] = complex<double>(0.0, 1.0) * kx;
+        this->lap[i] = (kx * kx < 1e-13) ? 0.0 : -kx * kx;
       }
     }
 
-    ~AdvectionDiffusionSweeper()
+    virtual ~AdvectionDiffusionSweeper()
     {
-      cout << "number of f1 evals: " << nf1evals << endl;
+      cout << "number of f1 evals: " << this->nf1evals << endl;
     }
+    //! @}
 
+    //! @{
     void exact(shared_ptr<Encapsulation<time>> q, time t)
     {
-      this->exact(as_vector<double,time>(q), t);
+      this->exact(as_vector<double, time>(q), t);
     }
 
     void exact(DVectorT& q, time t)
@@ -81,10 +98,10 @@ class AdvectionDiffusionSweeper
 
     void echo_error(time t, bool predict = false)
     {
-      auto& qend = as_vector<double,time>(this->get_end_state());
+      auto& qend = as_vector<double, time>(this->get_end_state());
       DVectorT qex(qend.size());
 
-      exact(qex, t);
+      this->exact(qex, t);
 
       double max = 0.0;
       for (size_t i = 0; i < qend.size(); i++) {
@@ -98,80 +115,111 @@ class AdvectionDiffusionSweeper
            << " (" << qend.size() << ", " << predict << ")"
            << endl;
 
-      errors.insert(pair<pair<size_t, size_t>, double>
-		    (pair<size_t, size_t>(n, k), max));
+      this->errors.insert(pair<pair<size_t, size_t>, double>(pair<size_t, size_t>(n, k), max));
     }
 
+    /**
+     * retrieve errors at iterations and time nodes
+     */
     error_map get_errors()
     {
-      return errors;
+      return this->errors;
     }
+    //! @}
 
-    void predict(bool initial)
+    //! @{
+    /**
+     * @copybrief pfasst::encap::IMEXSweeper::predict()
+     */
+    void predict(bool initial) override
     {
       pfasst::encap::IMEXSweeper<time>::predict(initial);
       time t  = this->get_controller()->get_time();
       time dt = this->get_controller()->get_time_step();
-      echo_error(t + dt, true);
+      this->echo_error(t + dt, true);
     }
 
-    void sweep()
+    /**
+     * @copybrief pfasst::encap::IMEXSweeper::sweep()
+     */
+    void sweep() override
     {
       pfasst::encap::IMEXSweeper<time>::sweep();
       time t  = this->get_controller()->get_time();
       time dt = this->get_controller()->get_time_step();
-      echo_error(t + dt);
+      this->echo_error(t + dt);
     }
-
-    void f1eval(shared_ptr<Encapsulation<time>> _f, shared_ptr<Encapsulation<time>> _q, time /*t*/)
+    //! @}
+
+    //! @{
+    /**
+     * @copybrief pfasst::encap::IMEXSweeper::f_expl_eval()
+     */
+    void f_expl_eval(shared_ptr<Encapsulation<time>> f_expl_encap,
+                     shared_ptr<Encapsulation<time>> u_encap,
+                     time t) override
     {
-      auto& q = as_vector<double,time>(_q);
-      auto& f = as_vector<double,time>(_f);
+      UNUSED(t);
+      auto& u = as_vector<double, time>(u_encap);
+      auto& f_expl = as_vector<double, time>(f_expl_encap);
 
-      double c = -v / double(q.size());
+      double c = -v / double(u.size());
 
-      auto* z = fft.forward(q);
-      for (size_t i = 0; i < q.size(); i++) {
-        z[i] *= c * ddx[i];
+      auto* z = this->fft.forward(u);
+      for (size_t i = 0; i < u.size(); i++) {
+        z[i] *= c * this->ddx[i];
       }
-      fft.backward(f);
+      this->fft.backward(f_expl);
 
-      nf1evals++;
+      this->nf1evals++;
     }
 
-    void f2eval(shared_ptr<Encapsulation<time>> _f, shared_ptr<Encapsulation<time>> _q, time /*t*/)
+    /**
+     * @copybrief pfasst::encap::IMEXSweeper::f_impl_eval()
+     */
+    void f_impl_eval(shared_ptr<Encapsulation<time>> f_impl_encap,
+                     shared_ptr<Encapsulation<time>> u_encap,
+                     time t) override
     {
-      auto& q = as_vector<double,time>(_q);
-      auto& f = as_vector<double,time>(_f);
+      UNUSED(t);
+      auto& u = as_vector<double, time>(u_encap);
+      auto& f_impl = as_vector<double, time>(f_impl_encap);
 
-      double c = nu / double(q.size());
+      double c = nu / double(u.size());
 
-      auto* z = fft.forward(q);
-      for (size_t i = 0; i < q.size(); i++) {
-        z[i] *= c * lap[i];
+      auto* z = this->fft.forward(u);
+      for (size_t i = 0; i < u.size(); i++) {
+        z[i] *= c * this->lap[i];
       }
-      fft.backward(f);
+      this->fft.backward(f_impl);
     }
 
-    void f2comp(shared_ptr<Encapsulation<time>> _f, shared_ptr<Encapsulation<time>> _q, time /*t*/, time dt,
-                shared_ptr<Encapsulation<time>> _rhs)
+    /**
+     * @copybrief pfasst::encap::IMEXSweeper::impl_solve()
+     */
+    void impl_solve(shared_ptr<Encapsulation<time>> f_impl_encap,
+                    shared_ptr<Encapsulation<time>> u_encap,
+                    time t, time dt,
+                    shared_ptr<Encapsulation<time>> rhs_encap) override
     {
-      auto& q = as_vector<double,time>(_q);
-      auto& f = as_vector<double,time>(_f);
-      auto& rhs = as_vector<double,time>(_rhs);
+      UNUSED(t);
+      auto& u = as_vector<double, time>(u_encap);
+      auto& f_impl = as_vector<double, time>(f_impl_encap);
+      auto& rhs = as_vector<double, time>(rhs_encap);
 
-      auto* z = fft.forward(rhs);
-      for (size_t i = 0; i < q.size(); i++) {
-        z[i] /= (1.0 - nu * double(dt) * lap[i]) * double(q.size());
-      }
-      fft.backward(q);
+      double c = nu * double(dt);
 
-      for (size_t i = 0; i < q.size(); i++) {
-        f[i] = (q[i] - rhs[i]) / double(dt);
+      auto* z = this->fft.forward(rhs);
+      for (size_t i = 0; i < u.size(); i++) {
+        z[i] /= (1.0 - c * this->lap[i]) * double(u.size());
       }
+      this->fft.backward(u);
 
+      for (size_t i = 0; i < u.size(); i++) {
+        f_impl[i] = (u[i] - rhs[i]) / double(dt);
+      }
     }
-
+    //! @}
 };
 
 #endif

examples/advection_diffusion/fft.hpp

@@ -72,7 +72,6 @@ class FFT
         x[i] = real(wk->z[i]);
       }
     }
-
 };
 
 #endif

examples/advection_diffusion/mpi_pfasst.cpp

@@ -24,10 +24,8 @@ using namespace pfasst;
 using namespace pfasst::encap;
 using namespace pfasst::mpi;
 
-int main(int argc, char** argv)
+error_map run_mpi_pfasst()
 {
-  MPI_Init(&argc, &argv);
-
   const size_t nsteps = 4;
   const double dt     = 0.01;
   const size_t niters = 4;
@@ -64,7 +62,16 @@ int main(int argc, char** argv)
   pf.set_nsweeps({2, 1});
   pf.run();
 
-  fftw_cleanup();
+  auto fine = pf.get_finest<AdvectionDiffusionSweeper<>>();
+  return fine->get_errors();
+}
 
+#ifndef PFASST_UNIT_TESTING
+int main(int argc, char** argv)
+{
+  MPI_Init(&argc, &argv);
+  run_mpi_pfasst();
+  fftw_cleanup();
   MPI_Finalize();
 }
+#endif

examples/advection_diffusion/serial_mlsdc_autobuild.cpp

@@ -7,8 +7,11 @@
  * controller.
  */
 
-#include <memory>
+#include <cstdlib>
 #include <cassert>
+#include <memory>
+#include <vector>
+#include <string>
 
 #include <fftw3.h>
 
@@ -24,6 +27,7 @@ using namespace std;
 using namespace pfasst;
 using namespace pfasst::encap;
 
+
 int main(int /*argc*/, char** /*argv*/)
 {
   MLSDC<> mlsdc;

examples/advection_diffusion/spectral_transfer_1d.hpp

@@ -6,13 +6,17 @@
 #define _SPECTRAL_TRANSFER_1D_HPP_
 
 #include <cassert>
+#include <cstdlib>
 #include <memory>
 
 #include <pfasst/encap/vector.hpp>
 #include <pfasst/encap/poly_interp.hpp>
 
 #include "fft.hpp"
 
+using namespace std;
+
+
 template<typename time = pfasst::time_precision>
 class SpectralTransfer1D
   : public pfasst::encap::PolyInterpMixin<time>
@@ -22,14 +26,13 @@ class SpectralTransfer1D
     FFT fft;
 
   public:
-
-    void interpolate(shared_ptr<Encapsulation> dst, shared_ptr<const Encapsulation> src)
+    void interpolate(shared_ptr<Encapsulation> dst, shared_ptr<const Encapsulation> src) override
     {
-      auto& fine = as_vector<double,time>(dst);
-      auto& crse = as_vector<double,time>(src);
+      auto& fine = pfasst::encap::as_vector<double, time>(dst);
+      auto& crse = pfasst::encap::as_vector<double, time>(src);
 
-      auto* crse_z = fft.forward(crse);
-      auto* fine_z = fft.get_workspace(fine.size())->z;
+      auto* crse_z = this->fft.forward(crse);
+      auto* fine_z = this->fft.get_workspace(fine.size())->z;
 
       for (size_t i = 0; i < fine.size(); i++) {
         fine_z[i] = 0.0;
@@ -45,21 +48,20 @@ class SpectralTransfer1D
         fine_z[fine.size() - crse.size() / 2 + i] = c * crse_z[crse.size() / 2 + i];
       }
 
-      fft.backward(fine);
+      this->fft.backward(fine);
     }
 
-    void restrict(shared_ptr<Encapsulation> dst, shared_ptr<const Encapsulation> src)
+    void restrict(shared_ptr<Encapsulation> dst, shared_ptr<const Encapsulation> src) override
     {
-      auto& fine = as_vector<double,time>(src);
-      auto& crse = as_vector<double,time>(dst);
+      auto& fine = pfasst::encap::as_vector<double, time>(src);
+      auto& crse = pfasst::encap::as_vector<double, time>(dst);
 
       size_t xrat = fine.size() / crse.size();
 
       for (size_t i = 0; i < crse.size(); i++) {
         crse[i] = fine[xrat*i];
       }
     }
-
 };
 
 #endif

examples/advection_diffusion/vanilla_sdc.cpp

@@ -4,6 +4,7 @@
  * This example uses a vanilla SDC sweeper.
  */
 
+#include <cstdlib>
 #include <memory>
 
 #include <fftw3.h>