Merge PR #205 into development

torbjoernk · torbjoernk · commit 5f94593ce8fa · 2015-05-03T10:40:45.000+02:00
* upstream/pr/205/head:
  implicit: Split decl/impl.
  imex: Tidy includes.

Signed-off-by: Torbjörn Klatt &lt;t.klatt@fz-juelich.de&gt;
diff --git a/include/pfasst/encap/imex_sweeper.hpp b/include/pfasst/encap/imex_sweeper.hpp
@@ -1,14 +1,13 @@
 #ifndef _PFASST_ENCAP_IMEX_SWEEPER_HPP_
 #define _PFASST_ENCAP_IMEX_SWEEPER_HPP_
 
-#include <vector>
 #include <memory>
-using namespace std;
+#include <vector>
 
 #include "pfasst/encap/encapsulation.hpp"
 #include "pfasst/encap/encap_sweeper.hpp"
-#include "pfasst/encap/vector.hpp"
 
+using namespace std;
 
 namespace pfasst
 {
@@ -21,12 +20,12 @@ namespace pfasst
      *
      * This IMEX sweeper is for ODEs of the form
      * \\( \\dot{U} = F_{\\rm expl}(t,U) + F_{\\rm impl}(t, U) \\).
-     * To reduce complexity and computational effort the non-stiff part is treated explicitly and 
+     * To reduce complexity and computational effort the non-stiff part is treated explicitly and
      * the stiff part implicitly.
      *
-     * This sweeper requires three interfaces to be implemented: two routines to evaluate the 
-     * explicit \\( F_{\\rm expl} \\) and implicit \\( F_{\\rm impl} \\) pieces for a given state, 
-     * and one routine that solves (perhaps with an external solver) the backward-Euler equation 
+     * This sweeper requires three interfaces to be implemented: two routines to evaluate the
+     * explicit \\( F_{\\rm expl} \\) and implicit \\( F_{\\rm impl} \\) pieces for a given state,
+     * and one routine that solves (perhaps with an external solver) the backward-Euler equation
      * \\( U^{n+1} - \\Delta t F_{\\rm impl}(U^{n+1}) = RHS \\) for \\( U^{n+1} \\).
      *
      * @tparam time precision type of the time dimension
diff --git a/include/pfasst/encap/implicit_sweeper.hpp b/include/pfasst/encap/implicit_sweeper.hpp
@@ -1,66 +1,16 @@
-
 #ifndef _PFASST_ENCAP_IMPLICIT_SWEEPER_HPP_
 #define _PFASST_ENCAP_IMPLICIT_SWEEPER_HPP_
 
-#include <cstdlib>
-#include <cassert>
-#include <vector>
 #include <memory>
+#include <vector>
 
-#include "../globals.hpp"
-#include "../quadrature.hpp"
-#include "encapsulation.hpp"
-#include "encap_sweeper.hpp"
-#include "vector.hpp"
+#include "pfasst/encap/encapsulation.hpp"
+#include "pfasst/encap/encap_sweeper.hpp"
 
 using namespace std;
 
 namespace pfasst
 {
-
-  template<typename scalar>
-  using lu_pair = pair< Matrix<scalar>, Matrix<scalar> >;
-
-  template<typename scalar>
-  static lu_pair<scalar> lu_decomposition(const Matrix<scalar>& A)
-  {
-    assert(A.rows() == A.cols());
-
-    auto n = A.rows();
-
-    Matrix<scalar> L = Matrix<scalar>::Zero(n, n);
-    Matrix<scalar> U = Matrix<scalar>::Zero(n, n);
-
-    if (A.rows() == 1) {
-
-      L(0, 0) = 1.0;
-      U(0, 0) = A(0,0);
-
-    } else {
-
-      // first row of U is first row of A
-      auto U12 = A.block(0, 1, 1, n-1);
-
-      // first column of L is first column of A / a11
-      auto L21 = A.block(1, 0, n-1, 1) / A(0, 0);
-
-      // remove first row and column and recurse
-      auto A22  = A.block(1, 1, n-1, n-1);
-      Matrix<scalar> tmp = A22 - L21 * U12;
-      auto LU22 = lu_decomposition(tmp);
-
-      L(0, 0) = 1.0;
-      U(0, 0) = A(0, 0);
-      L.block(1, 0, n-1, 1) = L21;
-      U.block(0, 1, 1, n-1) = U12;
-      L.block(1, 1, n-1, n-1) = get<0>(LU22);
-      U.block(1, 1, n-1, n-1) = get<1>(LU22);
-
-    }
-
-    return lu_pair<scalar>(L, U);
-  }
-
   namespace encap
   {
     using pfasst::encap::Encapsulation;
@@ -90,32 +40,7 @@ namespace pfasst
 
         Matrix<time> q_tilde;
 
-        /**
-        * Set end state to \\( U_0 + \\int F_{expl} + F_{expl} \\).
-        */
-        void set_end_state()
-        {
-          if (this->quadrature->right_is_node()) {
-            this->end_state->copy(this->state.back());
-          } else {
-            vector<shared_ptr<Encapsulation<time>>> dst = { this->end_state };
-            dst[0]->copy(this->start_state);
-            dst[0]->mat_apply(dst, this->get_controller()->get_time_step(), this->quadrature->get_b_mat(), this->fs_impl, false);
-          }
-        }
-
-        /**
-        * Augment nodes: nodes <- [t0] + dt * nodes
-        */
-        vector<time> augment(time t0, time dt, vector<time> const & nodes)
-        {
-          vector<time> t(1 + nodes.size());
-          t[0] = t0;
-          for (size_t m = 0; m < nodes.size(); m++) {
-            t[m+1] = t0 + dt * nodes[m];
-          }
-          return t;
-        }
+        void set_end_state();
 
       public:
         //! @{
@@ -126,37 +51,8 @@ namespace pfasst
         //! @{
         /**
          * @copydoc ISweeper::setup(bool)
-         *
          */
-        virtual void setup(bool coarse) override
-        {
-          pfasst::encap::EncapSweeper<time>::setup(coarse);
-
-          auto const nodes = this->quadrature->get_nodes();
-          auto const num_nodes = this->quadrature->get_num_nodes();
-
-          if (this->quadrature->left_is_node()) {
-            CLOG(INFO, "Sweeper") << "implicit sweeper shouldn't include left endpoint";
-            throw ValueError("implicit sweeper shouldn't include left endpoint");
-          }
-
-          for (size_t m = 0; m < num_nodes; m++) {
-            this->s_integrals.push_back(this->get_factory()->create(pfasst::encap::solution));
-            this->fs_impl.push_back(this->get_factory()->create(pfasst::encap::function));
-          }
-
-          Matrix<time> QT = this->quadrature->get_q_mat().transpose();
-          auto lu = lu_decomposition(QT);
-          auto L = get<0>(lu);
-          auto U = get<1>(lu);
-          this->q_tilde = U.transpose();
-
-          CLOG(DEBUG, "Sweeper") << "Q':" << endl << QT;
-          CLOG(DEBUG, "Sweeper") << "L:" << endl << L;
-          CLOG(DEBUG, "Sweeper") << "U:" << endl << U;
-          CLOG(DEBUG, "Sweeper") << "LU:" << endl << L * U;
-          CLOG(DEBUG, "Sweeper") << "q_tilde:" << endl << this->q_tilde;
-        }
+        virtual void setup(bool coarse) override;
 
         /**
          * Compute low-order provisional solution.
@@ -166,103 +62,33 @@ namespace pfasst
          * @param[in] initial if `true` the explicit and implicit part of the right hand side of the
          *     ODE get evaluated with the initial value
          */
-        virtual void predict(bool initial) override
-        {
-          UNUSED(initial);
-
-          auto const dt = this->get_controller()->get_time_step();
-          auto const t  = this->get_controller()->get_time();
-
-          CLOG(INFO, "Sweeper") << "predicting step " << this->get_controller()->get_step() + 1
-                                << " (t=" << t << ", dt=" << dt << ")";
-
-          auto const anodes = augment(t, dt, this->quadrature->get_nodes());
-          for (size_t m = 0; m < anodes.size() - 1; ++m) {
-            this->impl_solve(this->fs_impl[m], this->state[m], anodes[m], anodes[m+1] - anodes[m],
-                             m == 0 ? this->get_start_state() : this->state[m-1]);
-          }
-
-          this->set_end_state();
-        }
+        virtual void predict(bool initial) override;
 
         /**
          * Perform one SDC sweep/iteration.
          *
          * This computes a high-order solution from the previous iteration's function values and
          * corrects it using forward/backward Euler steps across the nodes.
          */
-        virtual void sweep() override
-        {
-          auto const dt = this->get_controller()->get_time_step();
-          auto const t  = this->get_controller()->get_time();
-
-          CLOG(INFO, "Sweeper") << "sweeping on step " << this->get_controller()->get_step() + 1
-                                << " in iteration " << this->get_controller()->get_iteration()
-                                << " (dt=" << dt << ")";
-
-          this->s_integrals[0]->mat_apply(this->s_integrals, dt, this->quadrature->get_s_mat(), this->fs_impl, true);
-          if (this->fas_corrections.size() > 0) {
-            for (size_t m = 0; m < this->s_integrals.size(); m++) {
-              this->s_integrals[m]->saxpy(1.0, this->fas_corrections[m]);
-            }
-          }
-
-          for (size_t m = 0; m < this->s_integrals.size(); m++) {
-            for (size_t n = 0; n < m; n++) {
-              this->s_integrals[m]->saxpy(-dt*this->q_tilde(m, n), this->fs_impl[n]);
-            }
-          }
-
-          shared_ptr<Encapsulation<time>> rhs = this->get_factory()->create(pfasst::encap::solution);
-
-          auto const anodes = augment(t, dt, this->quadrature->get_nodes());
-          for (size_t m = 0; m < anodes.size() - 1; ++m) {
-            auto const ds = anodes[m+1] - anodes[m];
-            rhs->copy(m == 0 ? this->get_start_state() : this->state[m-1]);
-            rhs->saxpy(1.0, this->s_integrals[m]);
-            rhs->saxpy(-ds, this->fs_impl[m]);
-            for (size_t n = 0; n < m; n++) {
-              rhs->saxpy(dt*this->q_tilde(m, n), this->fs_impl[n]);
-            }
-            this->impl_solve(this->fs_impl[m], this->state[m], anodes[m], ds, rhs);
-          }
-          this->set_end_state();
-        }
+        virtual void sweep() override;
 
         /**
          * Advance the end solution to start solution.
          */
-        virtual void advance() override
-        {
-          this->start_state->copy(this->end_state);
-        }
+        virtual void advance() override;
 
         /**
          * @copybrief EncapSweeper::evaluate()
          */
-        virtual void reevaluate(bool initial_only) override
-        {
-          if (initial_only) {
-            return;
-          }
-          auto const dt = this->get_controller()->get_time_step();
-          auto const t0 = this->get_controller()->get_time();
-          auto const nodes = this->quadrature->get_nodes();
-          for (size_t m = 0; m < nodes.size(); m++) {
-            this->f_impl_eval(this->fs_impl[m], this->state[m], t0 + dt * nodes[m]);
-          }
-        }
+        virtual void reevaluate(bool initial_only) override;
 
         /**
          * @copybrief EncapSweeper::integrate()
          *
          * @param[in] dt width of time interval to integrate over
          * @param[in,out] dst integrated values; will get zeroed out beforehand
          */
-        virtual void integrate(time dt, vector<shared_ptr<Encapsulation<time>>> dst) const override
-        {
-          dst[0]->mat_apply(dst, dt, this->quadrature->get_q_mat(), this->fs_impl, true);
-        }
+        virtual void integrate(time dt, vector<shared_ptr<Encapsulation<time>>> dst) const override;
         //! @}
 
         //! @{
@@ -316,4 +142,6 @@ namespace pfasst
   }  // ::pfasst::encap
 }  // ::pfasst
 
+#include "pfasst/encap/implicit_sweeper_impl.hpp"
+
 #endif
diff --git a/src/pfasst/encap/implicit_sweeper_impl.hpp b/src/pfasst/encap/implicit_sweeper_impl.hpp
