From: tcclevenger <tcleven@clemson.edu>
Date: Tue, 19 Feb 2019 23:15:20 +0000 (-0700)
Subject: add step-63.cc (no doc)
X-Git-Tag: v9.1.0-rc1~150^2~6
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=88d865167aaceba18ffeee838085c780eee40502;p=dealii.git

add step-63.cc (no doc)
---

diff --git a/examples/step-63/block_jacobi.prm b/examples/step-63/block_jacobi.prm
new file mode 100644
index 0000000000..a8dc7f5566
--- /dev/null
+++ b/examples/step-63/block_jacobi.prm
@@ -0,0 +1,14 @@
+# Finite Element degree
+set fe degree       =  1 
+
+# Smoother Type: sor|jacobi|block sor|block jacobi
+set smoother type   = block jacobi
+
+# Dof renumbering: no|random|downstream|upstream
+set dof renumbering = downstream
+
+# With streamline diffusion: true|false
+set with sd         = true
+
+# Output: true|false
+set output          = false
diff --git a/examples/step-63/block_sor.prm b/examples/step-63/block_sor.prm
new file mode 100644
index 0000000000..f5991834c4
--- /dev/null
+++ b/examples/step-63/block_sor.prm
@@ -0,0 +1,14 @@
+# Finite Element degree
+set fe degree       =  1 
+
+# Smoother Type: sor|jacobi|block sor|block jacobi
+set smoother type   = block sor
+
+# Dof renumbering: none|random|downstream|upstream
+set dof renumbering = downstream 
+
+# With streamline diffusion: true|false
+set with sd         = true
+
+# Output: true|false
+set output          = false
diff --git a/examples/step-63/jacobi.prm b/examples/step-63/jacobi.prm
new file mode 100644
index 0000000000..d6a55ff53b
--- /dev/null
+++ b/examples/step-63/jacobi.prm
@@ -0,0 +1,14 @@
+# Finite Element degree
+set fe degree       =  1
+
+# Smoother Type: sor|jacobi|block sor|block jacobi
+set smoother type   = jacobi
+
+# Dof renumbering: no|random|downstream|upstream
+set dof renumbering = downstream
+
+# With streamline diffusion: true|false
+set with sd         = true
+
+# Output: true|false
+set output          = false
diff --git a/examples/step-63/sor.prm b/examples/step-63/sor.prm
new file mode 100644
index 0000000000..47ab13c3d5
--- /dev/null
+++ b/examples/step-63/sor.prm
@@ -0,0 +1,15 @@
+# Finite Element degree
+set fe degree       = 1
+
+# Smoother Type: sor|jacobi|block sor|block jacobi
+set smoother type   = sor
+
+# Dof renumbering: no|random|downstream|upstream
+set dof renumbering = downstream  
+
+
+# With streamline diffusion: true|false
+set with sd         = true
+
+# Output: true|false
+set output          = false 
diff --git a/examples/step-63/step-63.cc b/examples/step-63/step-63.cc
index 2534596dc6..dcfcd1fa54 100644
--- a/examples/step-63/step-63.cc
+++ b/examples/step-63/step-63.cc
@@ -19,16 +19,1091 @@
 
 // @note: This is work in progress and will be an example for block smoothers
 // in geometric multigrid.
+#include <deal.II/base/tensor_function.h>
+#include <deal.II/base/work_stream.h>
+#include <deal.II/base/std_cxx14/memory.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/base/path_search.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/relaxation_block.h>
 
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/meshworker/mesh_loop.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+#include <deal.II/multigrid/multigrid.h>
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_smoother.h>
+#include <deal.II/multigrid/mg_matrix.h>
+
+
+
+#include <fstream>
 #include <iostream>
-#include <deal.II/base/exceptions.h>
 
 
-int main()
+
+#include <boost/random.hpp>
+#include <boost/random/uniform_int_distribution.hpp>
+
+
+
+namespace Step100
+{
+  using namespace dealii;
+
+
+  template <int dim>
+  struct ScratchData
+  {
+    ScratchData(const FiniteElement<dim> &fe,
+                const unsigned int        quadrature_degree)
+      : fe_values(fe,
+                  QGauss<dim>(quadrature_degree),
+                  update_values | update_gradients | update_hessians |
+                    update_quadrature_points | update_JxW_values)
+    {}
+
+    ScratchData(const ScratchData<dim> &scratch_data)
+      : fe_values(scratch_data.fe_values.get_fe(),
+                  scratch_data.fe_values.get_quadrature(),
+                  update_values | update_gradients | update_hessians |
+                    update_quadrature_points | update_JxW_values)
+    {}
+
+    FEValues<dim> fe_values;
+  };
+  struct CopyData
+  {
+    unsigned int level;
+    unsigned int dofs_per_cell;
+
+    FullMatrix<double>                   cell_matrix;
+    Vector<double>                       cell_rhs;
+    std::vector<types::global_dof_index> local_dof_indices;
+  };
+
+
+  struct Settings
+  {
+    bool try_parse(const std::string &prm_filename);
+
+    unsigned int fe_degree;
+    std::string  smoother_type;
+    std::string  dof_renum;
+    bool         with_sd;
+    bool         output;
+  };
+
+
+  bool Settings::try_parse(const std::string &prm_filename)
+  {
+    ParameterHandler prm;
+
+    prm.declare_entry("fe degree",
+                      "1",
+                      Patterns::Integer(0),
+                      "Finite Element degree");
+    prm.declare_entry("smoother type",
+                      "block sor",
+                      Patterns::Selection("sor|jacobi|block sor|block jacobi"),
+                      "Smoother Type: sor|jacobi|block sor|block jacobi");
+    prm.declare_entry("dof renumbering",
+                      "downstream",
+                      Patterns::Selection("none|random|downstream|upstream"),
+                      "Dof renumbering: none|random|downstream|upstream");
+    prm.declare_entry("with sd",
+                      "true",
+                      Patterns::Bool(),
+                      "With streamline diffusion: true|false");
+    prm.declare_entry("output",
+                      "true",
+                      Patterns::Bool(),
+                      "Generate graphical output: true|false");
+
+    try
+      {
+        prm.parse_input(prm_filename);
+      }
+    catch (const dealii::PathSearch::ExcFileNotFound &)
+      {
+        if (prm_filename.size() > 0)
+          std::cerr << "ERRROR: could not open the .prm file '" << prm_filename
+                    << "'" << std::endl;
+        else
+          std::cerr << "Usage: please pass a .prm file as the first argument"
+                    << std::endl;
+
+        prm.print_parameters(std::cout, ParameterHandler::Text);
+        return false;
+      }
+    this->fe_degree     = prm.get_integer("fe degree");
+    this->smoother_type = prm.get("smoother type");
+    this->dof_renum     = prm.get("dof renumbering");
+    this->with_sd       = prm.get_bool("with sd");
+    this->output        = prm.get_bool("output");
+
+    return true;
+  }
+
+
+
+  namespace
+  {
+    template <class Iterator, int dim>
+    struct CompareDownstream
+    {
+      /**
+       * Constructor.
+       */
+      CompareDownstream(const Tensor<1, dim> &dir)
+        : dir(dir)
+      {}
+      /**
+       * Return true if c1 less c2.
+       */
+      bool operator()(const Iterator &c1, const Iterator &c2) const
+      {
+        const Tensor<1, dim> diff = c2->center() - c1->center();
+        return (diff * dir > 0);
+      }
+
+    private:
+      /**
+       * Flow direction.
+       */
+      const Tensor<1, dim> dir;
+    };
+
+    // Functions for creating permutation of cells for output and Block
+    // smoothers
+    template <int dim>
+    std::vector<unsigned int>
+    create_downstream_order(const DoFHandler<dim> &dof,
+                            const Tensor<1, dim>   direction,
+                            const unsigned int     level)
+    {
+      std::vector<typename DoFHandler<dim>::level_cell_iterator> ordered_cells;
+      ordered_cells.reserve(dof.get_triangulation().n_cells(level));
+      const CompareDownstream<typename DoFHandler<dim>::level_cell_iterator,
+                              dim>
+        comparator(direction);
+
+      typename DoFHandler<dim>::level_cell_iterator cell = dof.begin(level);
+      typename DoFHandler<dim>::level_cell_iterator endc = dof.end(level);
+      for (; cell != endc; ++cell)
+        ordered_cells.push_back(cell);
+
+      std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
+
+      std::vector<unsigned> ordered_indices;
+      ordered_indices.reserve(dof.get_triangulation().n_cells(level));
+
+      for (unsigned int i = 0; i < ordered_cells.size(); ++i)
+        ordered_indices.push_back(ordered_cells[i]->index());
+
+      return ordered_indices;
+    }
+
+    template <int dim>
+    std::vector<unsigned int>
+    create_downstream_order(const DoFHandler<dim> &dof,
+                            const Tensor<1, dim>   direction)
+    {
+      std::vector<typename DoFHandler<dim>::active_cell_iterator> ordered_cells;
+      ordered_cells.reserve(dof.get_triangulation().n_active_cells());
+      const CompareDownstream<typename DoFHandler<dim>::active_cell_iterator,
+                              dim>
+        comparator(direction);
+
+      typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
+      typename DoFHandler<dim>::active_cell_iterator endc = dof.end();
+      for (; cell != endc; ++cell)
+        ordered_cells.push_back(cell);
+
+      std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
+
+      std::vector<unsigned int> ordered_indices;
+      ordered_indices.reserve(dof.get_triangulation().n_active_cells());
+
+      for (unsigned int i = 0; i < ordered_cells.size(); ++i)
+        ordered_indices.push_back(ordered_cells[i]->index());
+
+      return ordered_indices;
+    }
+
+
+
+    template <int dim>
+    std::vector<unsigned int> create_random_order(const DoFHandler<dim> &dof,
+                                                  const unsigned int     level)
+    {
+      const unsigned int n_cells = dof.get_triangulation().n_cells(level);
+
+      std::vector<unsigned int> ordered_cells;
+      ordered_cells.reserve(n_cells);
+
+      typename DoFHandler<dim>::cell_iterator cell = dof.begin(level);
+      typename DoFHandler<dim>::cell_iterator endc = dof.end(level);
+      for (; cell != endc; ++cell)
+        ordered_cells.push_back(cell->index());
+
+      // shuffle the elements; the following is essentially std::shuffle (which
+      // is new in C++11) but with a boost URNG
+      ::boost::mt19937 random_number_generator;
+      for (unsigned int i = 1; i < n_cells; ++i)
+        {
+          // get a random number between 0 and i (inclusive)
+          const unsigned int j =
+            ::boost::random::uniform_int_distribution<>(0, i)(
+              random_number_generator);
+
+          // if possible, swap the elements
+          if (i != j)
+            std::swap(ordered_cells[i], ordered_cells[j]);
+        }
+
+      return ordered_cells;
+    }
+
+    template <int dim>
+    std::vector<unsigned int> create_random_order(const DoFHandler<dim> &dof)
+    {
+      const unsigned int n_cells = dof.get_triangulation().n_active_cells();
+
+      std::vector<unsigned int> ordered_cells;
+      ordered_cells.reserve(n_cells);
+
+      typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
+      typename DoFHandler<dim>::active_cell_iterator endc = dof.end();
+      for (; cell != endc; ++cell)
+        ordered_cells.push_back(cell->index());
+
+      // shuffle the elements; the following is essentially std::shuffle (which
+      // is new in C++11) but with a boost URNG
+      ::boost::mt19937 random_number_generator;
+      for (unsigned int i = 1; i < n_cells; ++i)
+        {
+          // get a random number between 0 and i (inclusive)
+          const unsigned int j =
+            ::boost::random::uniform_int_distribution<>(0, i)(
+              random_number_generator);
+
+          // if possible, swap the elements
+          if (i != j)
+            std::swap(ordered_cells[i], ordered_cells[j]);
+        }
+
+      return ordered_cells;
+    }
+
+  } // namespace
+
+
+
+  // RHS and boundary is an adaptation from one in
+  // Finite Elements and Fast Iterative Solvers: with Applications
+  // in Incompressible Fluid Dynamics
+  template <int dim>
+  class RightHandSide : public Function<dim>
+  {
+  public:
+    RightHandSide()
+      : Function<dim>()
+    {}
+
+    virtual double value(const Point<dim> & p,
+                         const unsigned int component = 0) const;
+
+    virtual void value_list(const std::vector<Point<dim>> &points,
+                            std::vector<double> &          values,
+                            const unsigned int             component = 0) const;
+
+  private:
+    static const Point<dim> center_point;
+  };
+
+  template <int dim>
+  double RightHandSide<dim>::value(const Point<dim> & p,
+                                   const unsigned int component) const
+  {
+    Assert(component == 0, ExcIndexRange(component, 0, 1));
+    (void)component;
+    (void)p;
+
+    return 0.0;
+  }
+
+
+
+  template <int dim>
+  void RightHandSide<dim>::value_list(const std::vector<Point<dim>> &points,
+                                      std::vector<double> &          values,
+                                      const unsigned int component) const
+  {
+    Assert(values.size() == points.size(),
+           ExcDimensionMismatch(values.size(), points.size()));
+
+    for (unsigned int i = 0; i < points.size(); ++i)
+      values[i] = RightHandSide<dim>::value(points[i], component);
+  }
+
+
+
+  template <int dim>
+  class BoundaryValues : public Function<dim>
+  {
+  public:
+    BoundaryValues()
+      : Function<dim>()
+    {}
+
+    virtual double value(const Point<dim> & p,
+                         const unsigned int component = 0) const;
+
+    virtual void value_list(const std::vector<Point<dim>> &points,
+                            std::vector<double> &          values,
+                            const unsigned int             component = 0) const;
+  };
+
+
+  template <int dim>
+  double BoundaryValues<dim>::value(const Point<dim> & p,
+                                    const unsigned int component) const
+  {
+    Assert(component == 0, ExcIndexRange(component, 0, 1));
+    (void)component;
+
+    if (std::fabs(p[0] * p[0] + p[1] * p[1] - 0.3 * 0.3) <
+          1e-8                                        // around cylinder
+        || std::fabs(p[0] + 1) < 1e-8                 // x == -1
+        || std::fabs(p[1] - 1) < 1e-8                 // y == 1
+        || (std::fabs(p[1] + 1) < 1e-8 && p[0] < 0.5) // y == -1, x <= 0.5
+    )
+      {
+        return 0.0;
+      }
+    else if (std::fabs(p[0] - 1) < 1e-8                     // x = 1
+             || (std::fabs(p[1] + 1) < 1e-8 && p[0] >= 0.5) // y == -1, x > 0.5
+    )
+      {
+        return 1.0;
+      }
+    else
+      {
+        return 0.0;
+      }
+  }
+
+
+
+  template <int dim>
+  void BoundaryValues<dim>::value_list(const std::vector<Point<dim>> &points,
+                                       std::vector<double> &          values,
+                                       const unsigned int component) const
+  {
+    Assert(values.size() == points.size(),
+           ExcDimensionMismatch(values.size(), points.size()));
+
+    for (unsigned int i = 0; i < points.size(); ++i)
+      values[i] = BoundaryValues<dim>::value(points[i], component);
+  }
+
+  template <int dim>
+  double delta_value(const double         hk,
+                     const double         eps,
+                     const Tensor<1, dim> dir,
+                     const double         pk)
+  {
+    // Value defined in 'On discontinuity–capturing methods for
+    // convection–diffusion equations'
+    double Peclet = dir.norm() * hk / (2.0 * eps * pk);
+    double coth =
+      (1.0 + std::exp(-2.0 * Peclet)) / (1.0 - std::exp(-2.0 * Peclet));
+
+    return hk / (2.0 * dir.norm() * pk) * (coth - 1.0 / Peclet);
+  }
+
+
+  template <int dim>
+  class AdvectionProblem
+  {
+  public:
+    AdvectionProblem(Settings settings);
+    void run();
+
+  private:
+    void setup_system();
+
+    template <class IteratorType>
+    void assemble_cell(const IteratorType &cell,
+                       ScratchData<dim> &  scratch_data,
+                       CopyData &          copy_data);
+    void assemble_system_and_multigrid();
+
+    std::unique_ptr<MGSmoother<Vector<double>>> create_smoother();
+
+    void solve();
+    void refine_grid();
+    void output_results(const unsigned int cycle) const;
+
+    Triangulation<dim> triangulation;
+    DoFHandler<dim>    dof_handler;
+
+    FE_Q<dim>     fe;
+    MappingQ<dim> mapping;
+    unsigned int  quad_degree;
+
+    AffineConstraints<double> constraints;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    Vector<double> solution;
+    Vector<double> system_rhs;
+
+    MGLevelObject<SparsityPattern> mg_sparsity_patterns;
+    MGLevelObject<SparsityPattern> mg_interface_sparsity_patterns;
+
+    MGLevelObject<SparseMatrix<double>> mg_matrices;
+    MGLevelObject<SparseMatrix<double>> mg_interface_in;
+    MGLevelObject<SparseMatrix<double>> mg_interface_out;
+
+    MGConstrainedDoFs mg_constrained_dofs;
+
+    Settings       settings;
+    const double   epsilon;
+    Tensor<1, dim> advection_direction;
+  };
+
+
+
+  template <int dim>
+  AdvectionProblem<dim>::AdvectionProblem(Settings settings)
+    : triangulation(Triangulation<dim>::limit_level_difference_at_vertices)
+    , dof_handler(triangulation)
+    , fe(settings.fe_degree)
+    , mapping(settings.fe_degree)
+    , quad_degree(2 * fe.degree + 2)
+    , settings(settings)
+    , epsilon(0.005)
+  {
+    // Set Advection direction (problem is an adaptation from one in
+    // Finite Elements and Fast Iterative Solvers: with Applications
+    // in Incompressible Fluid Dynamics)
+    advection_direction[0] = -std::sin(numbers::PI / 6.0);
+    if (dim > 1)
+      advection_direction[1] = std::cos(numbers::PI / 6.0);
+    if (dim > 2)
+      advection_direction[2] = std::sin(numbers::PI / 6.0);
+  }
+
+
+
+  template <int dim>
+  void AdvectionProblem<dim>::setup_system()
+  {
+    const unsigned int n_levels = triangulation.n_levels();
+
+    dof_handler.distribute_dofs(fe);
+
+    // We could renumber the active dofs with DoFRenumbering::downstream()
+    // here, but the smoothers only act on multigrid levels and as such, this
+    // wouldn't matter. Instead, we will renumber the DoFs on each multigrid
+    // level below.
+
+    solution.reinit(dof_handler.n_dofs());
+    system_rhs.reinit(dof_handler.n_dofs());
+
+    constraints.clear();
+    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+
+    VectorTools::interpolate_boundary_values(
+      mapping, dof_handler, 0, BoundaryValues<dim>(), constraints);
+    VectorTools::interpolate_boundary_values(
+      mapping, dof_handler, 1, BoundaryValues<dim>(), constraints);
+    constraints.close();
+
+    DynamicSparsityPattern dsp(dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler,
+                                    dsp,
+                                    constraints,
+                                    /*keep_constrained_dofs = */ false);
+
+    sparsity_pattern.copy_from(dsp);
+
+    system_matrix.reinit(sparsity_pattern);
+
+
+    // Setup GMG DoFs
+    dof_handler.distribute_mg_dofs();
+
+    // Renumber DoFs on each level in downstream or upstream direction if
+    // needed. This is only necessary for point smoothers (SOR and Jacobi) as
+    // the block smoothers operate on cells (see create_smoother()):
+    if (settings.smoother_type == "sor" || settings.smoother_type == "jacobi")
+      {
+        if (settings.dof_renum == "downstream" ||
+            settings.dof_renum == "upstream")
+          {
+            const Tensor<1, dim> direction =
+              (settings.dof_renum == "upstream" ? -1.0 : 1.0) *
+              advection_direction;
+
+            for (unsigned int level = 0; level < n_levels; ++level)
+              DoFRenumbering::downstream(dof_handler,
+                                         level,
+                                         direction,
+                                         /*dof_wise_renumbering = */ true);
+          }
+        else if (settings.dof_renum == "random")
+          {
+            for (unsigned int level = 0; level < n_levels; ++level)
+              DoFRenumbering::random(dof_handler, level);
+          }
+        else
+          Assert(false, ExcNotImplemented());
+      }
+
+    mg_constrained_dofs.clear();
+    mg_constrained_dofs.initialize(dof_handler);
+
+    std::set<types::boundary_id> dirichlet_boundary_ids = {0, 1};
+    mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,
+                                                       dirichlet_boundary_ids);
+
+    mg_matrices.resize(0, n_levels - 1);
+    mg_matrices.clear_elements();
+    mg_interface_in.resize(0, n_levels - 1);
+    mg_interface_in.clear_elements();
+    mg_interface_out.resize(0, n_levels - 1);
+    mg_interface_out.clear_elements();
+    mg_sparsity_patterns.resize(0, n_levels - 1);
+    mg_interface_sparsity_patterns.resize(0, n_levels - 1);
+
+    for (unsigned int level = 0; level < n_levels; ++level)
+      {
+        {
+          DynamicSparsityPattern dsp(dof_handler.n_dofs(level),
+                                     dof_handler.n_dofs(level));
+          MGTools::make_sparsity_pattern(dof_handler, dsp, level);
+          mg_sparsity_patterns[level].copy_from(dsp);
+          mg_matrices[level].reinit(mg_sparsity_patterns[level]);
+        }
+        {
+          DynamicSparsityPattern dsp(dof_handler.n_dofs(level),
+                                     dof_handler.n_dofs(level));
+          MGTools::make_interface_sparsity_pattern(dof_handler,
+                                                   mg_constrained_dofs,
+                                                   dsp,
+                                                   level);
+          mg_interface_sparsity_patterns[level].copy_from(dsp);
+
+          // We need both interface in and out matrices since our problem is not
+          // symmetric
+          mg_interface_in[level].reinit(mg_interface_sparsity_patterns[level]);
+          mg_interface_out[level].reinit(mg_interface_sparsity_patterns[level]);
+        }
+      }
+  }
+
+
+  template <int dim>
+  template <class IteratorType>
+  void AdvectionProblem<dim>::assemble_cell(const IteratorType &cell,
+                                            ScratchData<dim> &  scratch_data,
+                                            CopyData &          copy_data)
+  {
+    const unsigned int level = cell->level();
+    copy_data.level          = level;
+
+    const unsigned int dofs_per_cell =
+      scratch_data.fe_values.get_fe().dofs_per_cell;
+    copy_data.dofs_per_cell = dofs_per_cell;
+
+    const unsigned int n_q_points =
+      scratch_data.fe_values.get_quadrature().size();
+    copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+
+    if (!cell->is_level_cell())
+      copy_data.cell_rhs.reinit(dofs_per_cell);
+
+    copy_data.local_dof_indices.resize(dofs_per_cell);
+    cell->get_active_or_mg_dof_indices(copy_data.local_dof_indices);
+
+    scratch_data.fe_values.reinit(cell);
+
+    const RightHandSide<dim> right_hand_side;
+    std::vector<double>      rhs_values(n_q_points);
+
+    right_hand_side.value_list(scratch_data.fe_values.get_quadrature_points(),
+                               rhs_values);
+
+    const double delta = settings.with_sd ? delta_value(cell->diameter(),
+                                                        epsilon,
+                                                        advection_direction,
+                                                        settings.fe_degree) :
+                                            0.0;
+
+    for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        {
+          for (unsigned int j = 0; j < dofs_per_cell; ++j)
+            {
+              copy_data.cell_matrix(i, j) +=
+                (epsilon * scratch_data.fe_values.shape_grad(j, q_point) *
+                 scratch_data.fe_values.shape_grad(i, q_point) *
+                 scratch_data.fe_values.JxW(q_point)) +
+                ((advection_direction *
+                  scratch_data.fe_values.shape_grad(j, q_point)) *
+                 scratch_data.fe_values.shape_value(i, q_point)) *
+                  scratch_data.fe_values.JxW(q_point);
+
+              if (settings.with_sd)
+                copy_data.cell_matrix(i, j) +=
+                  delta *
+                    (advection_direction *
+                     scratch_data.fe_values.shape_grad(j, q_point)) *
+                    (advection_direction *
+                     scratch_data.fe_values.shape_grad(i, q_point)) *
+                    scratch_data.fe_values.JxW(q_point) -
+                  delta * epsilon *
+                    trace(scratch_data.fe_values.shape_hessian(j, q_point)) *
+                    (advection_direction *
+                     scratch_data.fe_values.shape_grad(i, q_point)) *
+                    scratch_data.fe_values.JxW(q_point);
+            }
+          if (!cell->is_level_cell())
+            {
+              copy_data.cell_rhs(i) +=
+                scratch_data.fe_values.shape_value(i, q_point) *
+                rhs_values[q_point] * scratch_data.fe_values.JxW(q_point);
+              if (settings.with_sd)
+                copy_data.cell_rhs(i) +=
+                  delta * rhs_values[q_point] * advection_direction *
+                  scratch_data.fe_values.shape_grad(i, q_point) *
+                  scratch_data.fe_values.JxW(q_point);
+            }
+        }
+  }
+
+
+  template <int dim>
+  void AdvectionProblem<dim>::assemble_system_and_multigrid()
+  {
+    auto cell_worker_active =
+      [&](const decltype(dof_handler.begin_active()) &cell,
+          ScratchData<dim> &                          scratch_data,
+          CopyData &                                  copy_data) {
+        this->assemble_cell(cell, scratch_data, copy_data);
+      };
+
+
+
+    auto copier_active = [&](const CopyData &copy_data) {
+      constraints.distribute_local_to_global(copy_data.cell_matrix,
+                                             copy_data.cell_rhs,
+                                             copy_data.local_dof_indices,
+                                             system_matrix,
+                                             system_rhs);
+    };
+
+
+    MeshWorker::mesh_loop(dof_handler.begin_active(),
+                          dof_handler.end(),
+                          cell_worker_active,
+                          copier_active,
+                          ScratchData<dim>(fe, quad_degree),
+                          CopyData(),
+                          MeshWorker::assemble_own_cells);
+
+
+    // Assemble GMG
+    std::vector<AffineConstraints<double>> boundary_constraints(
+      triangulation.n_global_levels());
+    for (unsigned int level = 0; level < triangulation.n_global_levels();
+         ++level)
+      {
+        IndexSet dofset;
+        DoFTools::extract_locally_relevant_level_dofs(dof_handler,
+                                                      level,
+                                                      dofset);
+        boundary_constraints[level].reinit(dofset);
+        boundary_constraints[level].add_lines(
+          mg_constrained_dofs.get_refinement_edge_indices(level));
+        boundary_constraints[level].add_lines(
+          mg_constrained_dofs.get_boundary_indices(level));
+        boundary_constraints[level].close();
+      }
+
+    auto cell_worker_mg = [&](const decltype(dof_handler.begin_mg()) &cell,
+                              ScratchData<dim> &scratch_data,
+                              CopyData &        copy_data) {
+      this->assemble_cell(cell, scratch_data, copy_data);
+    };
+
+    auto copier_mg = [&](const CopyData &copy_data) {
+      boundary_constraints[copy_data.level].distribute_local_to_global(
+        copy_data.cell_matrix,
+        copy_data.local_dof_indices,
+        mg_matrices[copy_data.level]);
+
+      // If (i,j) is an interface_out dof pair, then (j,i) is an interface_in
+      // dof pair. Note: for interface_in, we load the transpose of the
+      // interface entries, i.e., the entry for dof pair (j,i) is stored in
+      // interface_in(i,j).
+      for (unsigned int i = 0; i < copy_data.dofs_per_cell; ++i)
+        for (unsigned int j = 0; j < copy_data.dofs_per_cell; ++j)
+          if (mg_constrained_dofs.is_interface_matrix_entry(
+                copy_data.level,
+                copy_data.local_dof_indices[i],
+                copy_data.local_dof_indices[j]))
+            {
+              mg_interface_out[copy_data.level].add(
+                copy_data.local_dof_indices[i],
+                copy_data.local_dof_indices[j],
+                copy_data.cell_matrix(i, j));
+              mg_interface_in[copy_data.level].add(
+                copy_data.local_dof_indices[i],
+                copy_data.local_dof_indices[j],
+                copy_data.cell_matrix(j, i));
+            }
+    };
+
+    MeshWorker::mesh_loop(dof_handler.begin_mg(),
+                          dof_handler.end_mg(),
+                          cell_worker_mg,
+                          copier_mg,
+                          ScratchData<dim>(fe, quad_degree),
+                          CopyData(),
+                          MeshWorker::assemble_own_cells);
+  }
+
+
+  template <int dim>
+  std::unique_ptr<MGSmoother<Vector<double>>>
+  AdvectionProblem<dim>::create_smoother()
+  {
+    if (settings.smoother_type == "sor")
+      {
+        typedef PreconditionSOR<SparseMatrix<double>> Smoother;
+
+        auto smoother =
+          std_cxx14::make_unique<MGSmootherPrecondition<SparseMatrix<double>,
+                                                        Smoother,
+                                                        Vector<double>>>();
+        smoother->initialize(mg_matrices,
+                             Smoother::AdditionalData(fe.degree == 1 ? 1.0 :
+                                                                       0.62));
+        smoother->set_steps(2);
+        return smoother;
+      }
+    else if (settings.smoother_type == "jacobi")
+      {
+        typedef PreconditionJacobi<SparseMatrix<double>> Smoother;
+        auto                                             smoother =
+          std_cxx14::make_unique<MGSmootherPrecondition<SparseMatrix<double>,
+                                                        Smoother,
+                                                        Vector<double>>>();
+        smoother->initialize(mg_matrices,
+                             Smoother::AdditionalData(fe.degree == 1 ? 0.6667 :
+                                                                       0.47));
+        smoother->set_steps(4);
+        return smoother;
+      }
+    else if (settings.smoother_type == "block sor")
+      {
+        typedef RelaxationBlockSOR<SparseMatrix<double>, double, Vector<double>>
+          Smoother;
+
+        static MGLevelObject<typename Smoother::AdditionalData> smoother_data;
+        smoother_data.resize(0, triangulation.n_levels() - 1);
+
+        for (unsigned int level = 0; level < triangulation.n_levels(); ++level)
+          {
+            DoFTools::make_cell_patches(smoother_data[level].block_list,
+                                        dof_handler,
+                                        level);
+
+            smoother_data[level].relaxation = 1.0;
+            smoother_data[level].inversion = PreconditionBlockBase<double>::svd;
+
+            std::vector<unsigned int> ordered_indices;
+            if (settings.dof_renum == "downstream")
+              ordered_indices = create_downstream_order(dof_handler,
+                                                        advection_direction,
+                                                        level);
+            else if (settings.dof_renum == "upstream")
+              ordered_indices =
+                create_downstream_order(dof_handler,
+                                        -1.0 * advection_direction,
+                                        level);
+            else if (settings.dof_renum == "random")
+              ordered_indices = create_random_order(dof_handler, level);
+            else if (settings.dof_renum == "none")
+              {
+                // Do nothing
+              }
+            else
+              AssertThrow(false, ExcNotImplemented());
+
+            smoother_data[level].order =
+              std::vector<std::vector<unsigned int>>(1, ordered_indices);
+          }
+
+        auto smoother =
+          std_cxx14::make_unique<MGSmootherPrecondition<SparseMatrix<double>,
+                                                        Smoother,
+                                                        Vector<double>>>();
+        smoother->initialize(mg_matrices, smoother_data);
+        smoother->set_steps(1);
+        return smoother;
+      }
+    else if (settings.smoother_type == "block jacobi")
+      {
+        typedef RelaxationBlockJacobi<SparseMatrix<double>,
+                                      double,
+                                      Vector<double>>
+          Smoother;
+
+        static MGLevelObject<typename Smoother::AdditionalData> smoother_data;
+        smoother_data.resize(0, triangulation.n_levels() - 1);
+
+        for (unsigned int level = 0; level < triangulation.n_levels(); ++level)
+          {
+            DoFTools::make_cell_patches(smoother_data[level].block_list,
+                                        dof_handler,
+                                        level);
+
+            smoother_data[level].relaxation = 0.25;
+            smoother_data[level].inversion = PreconditionBlockBase<double>::svd;
+
+            std::vector<unsigned int> ordered_indices;
+            if (settings.dof_renum == "downstream")
+              ordered_indices = create_downstream_order(dof_handler,
+                                                        advection_direction,
+                                                        level);
+            else if (settings.dof_renum == "upstream")
+              ordered_indices =
+                create_downstream_order(dof_handler,
+                                        -1.0 * advection_direction,
+                                        level);
+            else if (settings.dof_renum == "random")
+              ordered_indices = create_random_order(dof_handler, level);
+            else if (settings.dof_renum == "none")
+              {
+                // Do nothing
+              }
+            else
+              AssertThrow(false, ExcNotImplemented());
+
+            smoother_data[level].order =
+              std::vector<std::vector<unsigned int>>(1, ordered_indices);
+          }
+
+        auto smoother =
+          std_cxx14::make_unique<MGSmootherPrecondition<SparseMatrix<double>,
+                                                        Smoother,
+                                                        Vector<double>>>();
+        smoother->initialize(mg_matrices, smoother_data);
+        smoother->set_steps(2);
+        return smoother;
+      }
+    else
+      AssertThrow(false, ExcNotImplemented());
+  }
+
+
+  template <int dim>
+  void AdvectionProblem<dim>::solve()
+  {
+    Timer time;
+
+    const double       solve_tol = 1e-8 * system_rhs.l2_norm();
+    const unsigned int max_iters = 200;
+    SolverControl      solver_control(max_iters, solve_tol, true, true);
+    solver_control.enable_history_data();
+
+    typedef MGTransferPrebuilt<Vector<double>> Transfer;
+    Transfer                                   mg_transfer(mg_constrained_dofs);
+    mg_transfer.build_matrices(dof_handler);
+
+    FullMatrix<double> coarse_matrix;
+    coarse_matrix.copy_from(mg_matrices[0]);
+    MGCoarseGridHouseholder<double, Vector<double>> coarse_grid_solver;
+    coarse_grid_solver.initialize(coarse_matrix);
+
+    std::unique_ptr<MGSmoother<Vector<double>>> mg_smoother = create_smoother();
+
+    mg::Matrix<Vector<double>> mg_matrix(mg_matrices);
+    mg::Matrix<Vector<double>> mg_interface_matrix_in(mg_interface_in);
+    mg::Matrix<Vector<double>> mg_interface_matrix_out(mg_interface_out);
+
+    Multigrid<Vector<double>> mg(
+      mg_matrix, coarse_grid_solver, mg_transfer, *mg_smoother, *mg_smoother);
+    mg.set_edge_matrices(mg_interface_matrix_out, mg_interface_matrix_in);
+
+    PreconditionMG<dim, Vector<double>, Transfer> preconditioner(dof_handler,
+                                                                 mg,
+                                                                 mg_transfer);
+
+
+    std::cout << "     Solving with GMRES to tol " << solve_tol << "..."
+              << std::endl;
+    SolverGMRES<> solver(solver_control);
+
+    time.restart();
+    solver.solve(system_matrix, solution, system_rhs, preconditioner);
+    time.stop();
+
+    std::cout << "          converged in " << solver_control.last_step()
+              << " iterations"
+              << " in " << time.last_wall_time() << " seconds " << std::endl;
+
+    constraints.distribute(solution);
+  }
+
+
+
+  template <int dim>
+  void AdvectionProblem<dim>::output_results(const unsigned int cycle) const
+  {
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler(dof_handler);
+    data_out.add_data_vector(solution, "solution");
+
+    // Here we generate an index for each cell to visualize the ordering used
+    // by the smoothers. Note that we do this only for the active cells
+    // instead of the levels, where the smoothers are actually used. For the
+    // point smoothers we renumber DoFs instead of cells, so this is only an
+    // approximation of what happens in reality. Finally, the random ordering
+    // is not the random ordering we actually use (see create_smoother() for
+    // that).
+    const unsigned int n_active_cells = triangulation.n_active_cells();
+    Vector<double>     cell_indices(n_active_cells);
+
+    {
+      // First generate a permutation vector for the cell indices:
+      std::vector<unsigned int> ordered_indices;
+      if (settings.dof_renum == "downstream")
+        {
+          ordered_indices =
+            create_downstream_order(dof_handler, advection_direction);
+        }
+      else if (settings.dof_renum == "upstream")
+        {
+          ordered_indices =
+            create_downstream_order(dof_handler, -1.0 * advection_direction);
+        }
+      else if (settings.dof_renum == "random")
+        {
+          ordered_indices = create_random_order(dof_handler);
+        }
+      else if (settings.dof_renum == "none")
+        {
+          ordered_indices.resize(n_active_cells);
+          for (unsigned int i = 0; i < n_active_cells; ++i)
+            ordered_indices[i] = i;
+        }
+      else
+        AssertThrow(false, ExcNotImplemented());
+
+      // Then copy the permutation in ordered_indices into an output vector:
+      for (unsigned int i = 0; i < n_active_cells; ++i)
+        cell_indices(ordered_indices[i]) = static_cast<double>(i);
+    }
+
+    data_out.add_data_vector(cell_indices, "cell_index");
+
+    data_out.build_patches();
+
+    std::string filename =
+      "solution-" + Utilities::int_to_string(cycle) + ".vtu";
+    std::ofstream output(filename.c_str());
+    data_out.write_vtu(output);
+  }
+
+
+  template <int dim>
+  void AdvectionProblem<dim>::run()
+  {
+    for (unsigned int cycle = 0; cycle < (settings.fe_degree == 1 ? 7 : 5);
+         ++cycle)
+      {
+        std::cout << "  Cycle " << cycle << ':' << std::endl;
+
+        if (cycle == 0)
+          {
+            GridGenerator::hyper_cube_with_cylindrical_hole(
+              triangulation, 0.3, 1.0, 0.5, 1, false);
+            static const SphericalManifold<dim> manifold_description(
+              Point<dim>(0, 0));
+            triangulation.set_manifold(1, manifold_description);
+          }
+
+        triangulation.refine_global();
+
+        setup_system();
+
+        std::cout << "     Number of active cells:       "
+                  << triangulation.n_active_cells() << " ("
+                  << triangulation.n_levels() << " levels)" << std::endl;
+        std::cout << "     Number of degrees of freedom: "
+                  << dof_handler.n_dofs() << std::endl;
+
+        assemble_system_and_multigrid();
+
+        solve();
+
+        if (settings.output)
+          output_results(cycle);
+
+        std::cout << std::endl;
+      }
+  }
+} // namespace Step100
+
+
+int main(int argc, char *argv[])
 {
   try
     {
-      // do nothing.
+      Step100::Settings settings;
+      if (!settings.try_parse((argc > 1) ? (argv[1]) : ""))
+        return 0;
+
+      Step100::AdvectionProblem<2> advection_problem_2d(settings);
+      advection_problem_2d.run();
     }
   catch (std::exception &exc)
     {
@@ -41,7 +1116,6 @@ int main()
                 << "Aborting!" << std::endl
                 << "----------------------------------------------------"
                 << std::endl;
-
       return 1;
     }
   catch (...)