]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Convert step-31 to a simplex test 11503/head
authorPeter Munch <peterrmuench@gmail.com>
Mon, 11 Jan 2021 13:28:04 +0000 (14:28 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Tue, 2 Mar 2021 15:43:05 +0000 (16:43 +0100)
tests/simplex/step-31.cc [new file with mode: 0644]
tests/simplex/step-31.with_simplex_support=on.with_trilinos=true.output [new file with mode: 0644]

diff --git a/tests/simplex/step-31.cc b/tests/simplex/step-31.cc
new file mode 100644 (file)
index 0000000..5fa28f8
--- /dev/null
@@ -0,0 +1,1128 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Step-31 on a simplex mesh. Following incompatible modifications had to be
+// made:
+//  - number of maximum quadrature points had to be limited
+//  - no local refinement
+
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/block_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+#include <limits>
+#include <memory>
+
+#include "../tests.h"
+
+namespace Step31
+{
+  using namespace dealii;
+  namespace EquationData
+  {
+    constexpr double eta     = 1;
+    constexpr double kappa   = 1e-6;
+    constexpr double beta    = 10;
+    constexpr double density = 1;
+    template <int dim>
+    class TemperatureInitialValues : public Function<dim>
+    {
+    public:
+      TemperatureInitialValues()
+        : Function<dim>(1)
+      {}
+      virtual double
+      value(const Point<dim> & /*p*/,
+            const unsigned int /*component*/ = 0) const override
+      {
+        return 0;
+      }
+      virtual void
+      vector_value(const Point<dim> &p, Vector<double> &value) const override
+      {
+        for (unsigned int c = 0; c < this->n_components; ++c)
+          value(c) = TemperatureInitialValues<dim>::value(p, c);
+      }
+    };
+    template <int dim>
+    class TemperatureRightHandSide : public Function<dim>
+    {
+    public:
+      TemperatureRightHandSide()
+        : Function<dim>(1)
+      {}
+      virtual double
+      value(const Point<dim> & p,
+            const unsigned int component = 0) const override
+      {
+        (void)component;
+        Assert(component == 0,
+               ExcMessage("Invalid operation for a scalar function."));
+        Assert((dim == 2) || (dim == 3), ExcNotImplemented());
+        static const Point<dim> source_centers[3] = {
+          (dim == 2 ? Point<dim>(.3, .1) : Point<dim>(.3, .5, .1)),
+          (dim == 2 ? Point<dim>(.45, .1) : Point<dim>(.45, .5, .1)),
+          (dim == 2 ? Point<dim>(.75, .1) : Point<dim>(.75, .5, .1))};
+        static const double source_radius = (dim == 2 ? 1. / 32 : 1. / 8);
+        return ((source_centers[0].distance(p) < source_radius) ||
+                    (source_centers[1].distance(p) < source_radius) ||
+                    (source_centers[2].distance(p) < source_radius) ?
+                  1 :
+                  0);
+      }
+      virtual void
+      vector_value(const Point<dim> &p, Vector<double> &value) const override
+      {
+        for (unsigned int c = 0; c < this->n_components; ++c)
+          value(c) = TemperatureRightHandSide<dim>::value(p, c);
+      }
+    };
+  } // namespace EquationData
+  namespace LinearSolvers
+  {
+    template <class MatrixType, class PreconditionerType>
+    class InverseMatrix : public Subscriptor
+    {
+    public:
+      InverseMatrix(const MatrixType &        m,
+                    const PreconditionerType &preconditioner);
+      template <typename VectorType>
+      void
+      vmult(VectorType &dst, const VectorType &src) const;
+
+    private:
+      const SmartPointer<const MatrixType> matrix;
+      const PreconditionerType &           preconditioner;
+    };
+    template <class MatrixType, class PreconditionerType>
+    InverseMatrix<MatrixType, PreconditionerType>::InverseMatrix(
+      const MatrixType &        m,
+      const PreconditionerType &preconditioner)
+      : matrix(&m)
+      , preconditioner(preconditioner)
+    {}
+    template <class MatrixType, class PreconditionerType>
+    template <typename VectorType>
+    void
+    InverseMatrix<MatrixType, PreconditionerType>::vmult(
+      VectorType &      dst,
+      const VectorType &src) const
+    {
+      SolverControl        solver_control(src.size(), 1e-7 * src.l2_norm());
+      SolverCG<VectorType> cg(solver_control);
+      dst = 0;
+      try
+        {
+          cg.solve(*matrix, dst, src, preconditioner);
+        }
+      catch (std::exception &e)
+        {
+          Assert(false, ExcMessage(e.what()));
+        }
+    }
+    template <class PreconditionerTypeA, class PreconditionerTypeMp>
+    class BlockSchurPreconditioner : public Subscriptor
+    {
+    public:
+      BlockSchurPreconditioner(
+        const TrilinosWrappers::BlockSparseMatrix &S,
+        const InverseMatrix<TrilinosWrappers::SparseMatrix,
+                            PreconditionerTypeMp> &Mpinv,
+        const PreconditionerTypeA &                Apreconditioner);
+      void
+      vmult(TrilinosWrappers::MPI::BlockVector &      dst,
+            const TrilinosWrappers::MPI::BlockVector &src) const;
+
+    private:
+      const SmartPointer<const TrilinosWrappers::BlockSparseMatrix>
+        stokes_matrix;
+      const SmartPointer<const InverseMatrix<TrilinosWrappers::SparseMatrix,
+                                             PreconditionerTypeMp>>
+                                            m_inverse;
+      const PreconditionerTypeA &           a_preconditioner;
+      mutable TrilinosWrappers::MPI::Vector tmp;
+    };
+    template <class PreconditionerTypeA, class PreconditionerTypeMp>
+    BlockSchurPreconditioner<PreconditionerTypeA, PreconditionerTypeMp>::
+      BlockSchurPreconditioner(
+        const TrilinosWrappers::BlockSparseMatrix &S,
+        const InverseMatrix<TrilinosWrappers::SparseMatrix,
+                            PreconditionerTypeMp> &Mpinv,
+        const PreconditionerTypeA &                Apreconditioner)
+      : stokes_matrix(&S)
+      , m_inverse(&Mpinv)
+      , a_preconditioner(Apreconditioner)
+      , tmp(complete_index_set(stokes_matrix->block(1, 1).m()))
+    {}
+    template <class PreconditionerTypeA, class PreconditionerTypeMp>
+    void
+    BlockSchurPreconditioner<PreconditionerTypeA, PreconditionerTypeMp>::vmult(
+      TrilinosWrappers::MPI::BlockVector &      dst,
+      const TrilinosWrappers::MPI::BlockVector &src) const
+    {
+      a_preconditioner.vmult(dst.block(0), src.block(0));
+      stokes_matrix->block(1, 0).residual(tmp, dst.block(0), src.block(1));
+      tmp *= -1;
+      m_inverse->vmult(dst.block(1), tmp);
+    }
+  } // namespace LinearSolvers
+  template <int dim>
+  class BoussinesqFlowProblem
+  {
+  public:
+    BoussinesqFlowProblem();
+    void
+    run();
+
+  private:
+    void
+    setup_dofs();
+    void
+    assemble_stokes_preconditioner();
+    void
+    build_stokes_preconditioner();
+    void
+    assemble_stokes_system();
+    void
+    assemble_temperature_system(const double maximal_velocity);
+    void
+    assemble_temperature_matrix();
+    double
+    get_maximal_velocity() const;
+    std::pair<double, double>
+    get_extrapolated_temperature_range() const;
+    void
+    solve();
+    void
+    output_results() const;
+    void
+    refine_mesh(const unsigned int max_grid_level);
+    double
+    compute_viscosity(
+      const std::vector<double> &        old_temperature,
+      const std::vector<double> &        old_old_temperature,
+      const std::vector<Tensor<1, dim>> &old_temperature_grads,
+      const std::vector<Tensor<1, dim>> &old_old_temperature_grads,
+      const std::vector<double> &        old_temperature_laplacians,
+      const std::vector<double> &        old_old_temperature_laplacians,
+      const std::vector<Tensor<1, dim>> &old_velocity_values,
+      const std::vector<Tensor<1, dim>> &old_old_velocity_values,
+      const std::vector<double> &        gamma_values,
+      const double                       global_u_infty,
+      const double                       global_T_variation,
+      const double                       cell_diameter) const;
+    Triangulation<dim>                  triangulation;
+    double                              global_Omega_diameter;
+    MappingFE<dim>                      mapping;
+    const unsigned int                  stokes_degree;
+    FESystem<dim>                       stokes_fe;
+    DoFHandler<dim>                     stokes_dof_handler;
+    AffineConstraints<double>           stokes_constraints;
+    std::vector<IndexSet>               stokes_partitioning;
+    TrilinosWrappers::BlockSparseMatrix stokes_matrix;
+    TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix;
+    TrilinosWrappers::MPI::BlockVector  stokes_solution;
+    TrilinosWrappers::MPI::BlockVector  old_stokes_solution;
+    TrilinosWrappers::MPI::BlockVector  stokes_rhs;
+    const unsigned int                  temperature_degree;
+    FE_SimplexP<dim>                    temperature_fe;
+    DoFHandler<dim>                     temperature_dof_handler;
+
+    QGaussSimplex<dim> quadrature_temperature;
+    QGaussSimplex<dim> quadrature_stokes;
+
+    AffineConstraints<double>      temperature_constraints;
+    TrilinosWrappers::SparseMatrix temperature_mass_matrix;
+    TrilinosWrappers::SparseMatrix temperature_stiffness_matrix;
+    TrilinosWrappers::SparseMatrix temperature_matrix;
+    TrilinosWrappers::MPI::Vector  temperature_solution;
+    TrilinosWrappers::MPI::Vector  old_temperature_solution;
+    TrilinosWrappers::MPI::Vector  old_old_temperature_solution;
+    TrilinosWrappers::MPI::Vector  temperature_rhs;
+    double                         time_step;
+    double                         old_time_step;
+    unsigned int                   timestep_number;
+    std::shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
+    std::shared_ptr<TrilinosWrappers::PreconditionIC>  Mp_preconditioner;
+    bool                                               rebuild_stokes_matrix;
+    bool rebuild_temperature_matrices;
+    bool rebuild_stokes_preconditioner;
+  };
+  template <int dim>
+  BoussinesqFlowProblem<dim>::BoussinesqFlowProblem()
+    : triangulation(/*Triangulation<dim>::maximum_smoothing*/)
+    , global_Omega_diameter(std::numeric_limits<double>::quiet_NaN())
+    , mapping(FE_SimplexP<dim>(1))
+    , stokes_degree(1)
+    , stokes_fe(FE_SimplexP<dim>(stokes_degree + 1),
+                dim,
+                FE_SimplexP<dim>(stokes_degree),
+                1)
+    , stokes_dof_handler(triangulation)
+    , temperature_degree(2)
+    , temperature_fe(temperature_degree)
+    , temperature_dof_handler(triangulation)
+    , quadrature_temperature(temperature_degree + 1)
+    , quadrature_stokes(stokes_degree + 2)
+    , time_step(0)
+    , old_time_step(0)
+    , timestep_number(0)
+    , rebuild_stokes_matrix(true)
+    , rebuild_temperature_matrices(true)
+    , rebuild_stokes_preconditioner(true)
+  {}
+  template <int dim>
+  double
+  BoussinesqFlowProblem<dim>::get_maximal_velocity() const
+  {
+    const auto &                     quadrature_formula = quadrature_stokes;
+    const unsigned int               n_q_points = quadrature_formula.size();
+    FEValues<dim>                    fe_values(mapping,
+                            stokes_fe,
+                            quadrature_formula,
+                            update_values);
+    std::vector<Tensor<1, dim>>      velocity_values(n_q_points);
+    double                           max_velocity = 0;
+    const FEValuesExtractors::Vector velocities(0);
+    for (const auto &cell : stokes_dof_handler.active_cell_iterators())
+      {
+        fe_values.reinit(cell);
+        fe_values[velocities].get_function_values(stokes_solution,
+                                                  velocity_values);
+        for (unsigned int q = 0; q < n_q_points; ++q)
+          max_velocity = std::max(max_velocity, velocity_values[q].norm());
+      }
+    return max_velocity;
+  }
+  template <int dim>
+  std::pair<double, double>
+  BoussinesqFlowProblem<dim>::get_extrapolated_temperature_range() const
+  {
+    const auto &        quadrature_formula = quadrature_temperature;
+    const unsigned int  n_q_points         = quadrature_formula.size();
+    FEValues<dim>       fe_values(mapping,
+                            temperature_fe,
+                            quadrature_formula,
+                            update_values);
+    std::vector<double> old_temperature_values(n_q_points);
+    std::vector<double> old_old_temperature_values(n_q_points);
+    if (timestep_number != 0)
+      {
+        double min_temperature = std::numeric_limits<double>::max(),
+               max_temperature = -std::numeric_limits<double>::max();
+        for (const auto &cell : temperature_dof_handler.active_cell_iterators())
+          {
+            fe_values.reinit(cell);
+            fe_values.get_function_values(old_temperature_solution,
+                                          old_temperature_values);
+            fe_values.get_function_values(old_old_temperature_solution,
+                                          old_old_temperature_values);
+            for (unsigned int q = 0; q < n_q_points; ++q)
+              {
+                const double temperature =
+                  (1. + time_step / old_time_step) * old_temperature_values[q] -
+                  time_step / old_time_step * old_old_temperature_values[q];
+                min_temperature = std::min(min_temperature, temperature);
+                max_temperature = std::max(max_temperature, temperature);
+              }
+          }
+        return std::make_pair(min_temperature, max_temperature);
+      }
+    else
+      {
+        double min_temperature = std::numeric_limits<double>::max(),
+               max_temperature = -std::numeric_limits<double>::max();
+        for (const auto &cell : temperature_dof_handler.active_cell_iterators())
+          {
+            fe_values.reinit(cell);
+            fe_values.get_function_values(old_temperature_solution,
+                                          old_temperature_values);
+            for (unsigned int q = 0; q < n_q_points; ++q)
+              {
+                const double temperature = old_temperature_values[q];
+                min_temperature = std::min(min_temperature, temperature);
+                max_temperature = std::max(max_temperature, temperature);
+              }
+          }
+        return std::make_pair(min_temperature, max_temperature);
+      }
+  }
+  template <int dim>
+  double
+  BoussinesqFlowProblem<dim>::compute_viscosity(
+    const std::vector<double> &        old_temperature,
+    const std::vector<double> &        old_old_temperature,
+    const std::vector<Tensor<1, dim>> &old_temperature_grads,
+    const std::vector<Tensor<1, dim>> &old_old_temperature_grads,
+    const std::vector<double> &        old_temperature_laplacians,
+    const std::vector<double> &        old_old_temperature_laplacians,
+    const std::vector<Tensor<1, dim>> &old_velocity_values,
+    const std::vector<Tensor<1, dim>> &old_old_velocity_values,
+    const std::vector<double> &        gamma_values,
+    const double                       global_u_infty,
+    const double                       global_T_variation,
+    const double                       cell_diameter) const
+  {
+    constexpr double beta  = 0.017 * dim;
+    constexpr double alpha = 1.0;
+    if (global_u_infty == 0)
+      return 5e-3 * cell_diameter;
+    const unsigned int n_q_points   = old_temperature.size();
+    double             max_residual = 0;
+    double             max_velocity = 0;
+    for (unsigned int q = 0; q < n_q_points; ++q)
+      {
+        const Tensor<1, dim> u =
+          (old_velocity_values[q] + old_old_velocity_values[q]) / 2;
+        const double dT_dt =
+          (old_temperature[q] - old_old_temperature[q]) / old_time_step;
+        const double u_grad_T =
+          u * (old_temperature_grads[q] + old_old_temperature_grads[q]) / 2;
+        const double kappa_Delta_T =
+          EquationData::kappa *
+          (old_temperature_laplacians[q] + old_old_temperature_laplacians[q]) /
+          2;
+        const double residual =
+          std::abs((dT_dt + u_grad_T - kappa_Delta_T - gamma_values[q]) *
+                   std::pow((old_temperature[q] + old_old_temperature[q]) / 2,
+                            alpha - 1.));
+        max_residual = std::max(residual, max_residual);
+        max_velocity = std::max(std::sqrt(u * u), max_velocity);
+      }
+    const double c_R            = std::pow(2., (4. - 2 * alpha) / dim);
+    const double global_scaling = c_R * global_u_infty * global_T_variation *
+                                  std::pow(global_Omega_diameter, alpha - 2.);
+    return (
+      beta * max_velocity *
+      std::min(cell_diameter,
+               std::pow(cell_diameter, alpha) * max_residual / global_scaling));
+  }
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::setup_dofs()
+  {
+    std::vector<unsigned int> stokes_sub_blocks(dim + 1, 0);
+    stokes_sub_blocks[dim] = 1;
+    {
+      stokes_dof_handler.distribute_dofs(stokes_fe);
+      DoFRenumbering::component_wise(stokes_dof_handler, stokes_sub_blocks);
+      stokes_constraints.clear();
+      DoFTools::make_hanging_node_constraints(stokes_dof_handler,
+                                              stokes_constraints);
+      std::set<types::boundary_id> no_normal_flux_boundaries;
+      no_normal_flux_boundaries.insert(0);
+      VectorTools::compute_no_normal_flux_constraints(stokes_dof_handler,
+                                                      0,
+                                                      no_normal_flux_boundaries,
+                                                      stokes_constraints,
+                                                      mapping);
+      stokes_constraints.close();
+    }
+    {
+      temperature_dof_handler.distribute_dofs(temperature_fe);
+      temperature_constraints.clear();
+      DoFTools::make_hanging_node_constraints(temperature_dof_handler,
+                                              temperature_constraints);
+      temperature_constraints.close();
+    }
+    const std::vector<types::global_dof_index> stokes_dofs_per_block =
+      DoFTools::count_dofs_per_fe_block(stokes_dof_handler, stokes_sub_blocks);
+    const unsigned int n_u = stokes_dofs_per_block[0],
+                       n_p = stokes_dofs_per_block[1],
+                       n_T = temperature_dof_handler.n_dofs();
+    deallog << "Number of active cells: " << triangulation.n_active_cells()
+            << " (on " << triangulation.n_levels() << " levels)" << std::endl
+            << "Number of degrees of freedom: " << n_u + n_p + n_T << " ("
+            << n_u << '+' << n_p << '+' << n_T << ')' << std::endl
+            << std::endl;
+    stokes_partitioning.resize(2);
+    stokes_partitioning[0] = complete_index_set(n_u);
+    stokes_partitioning[1] = complete_index_set(n_p);
+    {
+      stokes_matrix.clear();
+      BlockDynamicSparsityPattern dsp(2, 2);
+      dsp.block(0, 0).reinit(n_u, n_u);
+      dsp.block(0, 1).reinit(n_u, n_p);
+      dsp.block(1, 0).reinit(n_p, n_u);
+      dsp.block(1, 1).reinit(n_p, n_p);
+      dsp.collect_sizes();
+      Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
+      for (unsigned int c = 0; c < dim + 1; ++c)
+        for (unsigned int d = 0; d < dim + 1; ++d)
+          if (!((c == dim) && (d == dim)))
+            coupling[c][d] = DoFTools::always;
+          else
+            coupling[c][d] = DoFTools::none;
+      DoFTools::make_sparsity_pattern(
+        stokes_dof_handler, coupling, dsp, stokes_constraints, false);
+      stokes_matrix.reinit(dsp);
+    }
+    {
+      Amg_preconditioner.reset();
+      Mp_preconditioner.reset();
+      stokes_preconditioner_matrix.clear();
+      BlockDynamicSparsityPattern dsp(2, 2);
+      dsp.block(0, 0).reinit(n_u, n_u);
+      dsp.block(0, 1).reinit(n_u, n_p);
+      dsp.block(1, 0).reinit(n_p, n_u);
+      dsp.block(1, 1).reinit(n_p, n_p);
+      dsp.collect_sizes();
+      Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
+      for (unsigned int c = 0; c < dim + 1; ++c)
+        for (unsigned int d = 0; d < dim + 1; ++d)
+          if (c == d)
+            coupling[c][d] = DoFTools::always;
+          else
+            coupling[c][d] = DoFTools::none;
+      DoFTools::make_sparsity_pattern(
+        stokes_dof_handler, coupling, dsp, stokes_constraints, false);
+      stokes_preconditioner_matrix.reinit(dsp);
+    }
+    {
+      temperature_mass_matrix.clear();
+      temperature_stiffness_matrix.clear();
+      temperature_matrix.clear();
+      DynamicSparsityPattern dsp(n_T, n_T);
+      DoFTools::make_sparsity_pattern(temperature_dof_handler,
+                                      dsp,
+                                      temperature_constraints,
+                                      false);
+      temperature_matrix.reinit(dsp);
+      temperature_mass_matrix.reinit(temperature_matrix);
+      temperature_stiffness_matrix.reinit(temperature_matrix);
+    }
+    IndexSet temperature_partitioning = complete_index_set(n_T);
+    stokes_solution.reinit(stokes_partitioning, MPI_COMM_WORLD);
+    old_stokes_solution.reinit(stokes_partitioning, MPI_COMM_WORLD);
+    stokes_rhs.reinit(stokes_partitioning, MPI_COMM_WORLD);
+    temperature_solution.reinit(temperature_partitioning, MPI_COMM_WORLD);
+    old_temperature_solution.reinit(temperature_partitioning, MPI_COMM_WORLD);
+    old_old_temperature_solution.reinit(temperature_partitioning,
+                                        MPI_COMM_WORLD);
+    temperature_rhs.reinit(temperature_partitioning, MPI_COMM_WORLD);
+  }
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::assemble_stokes_preconditioner()
+  {
+    stokes_preconditioner_matrix          = 0;
+    const auto &       quadrature_formula = quadrature_stokes;
+    FEValues<dim>      stokes_fe_values(mapping,
+                                   stokes_fe,
+                                   quadrature_formula,
+                                   update_JxW_values | update_values |
+                                     update_gradients);
+    const unsigned int dofs_per_cell = stokes_fe.n_dofs_per_cell();
+    const unsigned int n_q_points    = quadrature_formula.size();
+    FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+    std::vector<Tensor<2, dim>>          grad_phi_u(dofs_per_cell);
+    std::vector<double>                  phi_p(dofs_per_cell);
+    const FEValuesExtractors::Vector     velocities(0);
+    const FEValuesExtractors::Scalar     pressure(dim);
+    for (const auto &cell : stokes_dof_handler.active_cell_iterators())
+      {
+        stokes_fe_values.reinit(cell);
+        local_matrix = 0;
+        for (unsigned int q = 0; q < n_q_points; ++q)
+          {
+            for (unsigned int k = 0; k < dofs_per_cell; ++k)
+              {
+                grad_phi_u[k] = stokes_fe_values[velocities].gradient(k, q);
+                phi_p[k]      = stokes_fe_values[pressure].value(k, q);
+              }
+            for (unsigned int i = 0; i < dofs_per_cell; ++i)
+              for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                local_matrix(i, j) +=
+                  (EquationData::eta *
+                     scalar_product(grad_phi_u[i], grad_phi_u[j]) +
+                   (1. / EquationData::eta) * phi_p[i] * phi_p[j]) *
+                  stokes_fe_values.JxW(q);
+          }
+        cell->get_dof_indices(local_dof_indices);
+        stokes_constraints.distribute_local_to_global(
+          local_matrix, local_dof_indices, stokes_preconditioner_matrix);
+      }
+  }
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::build_stokes_preconditioner()
+  {
+    if (rebuild_stokes_preconditioner == false)
+      return;
+    deallog << "   Rebuilding Stokes preconditioner..." << std::flush;
+    assemble_stokes_preconditioner();
+    Amg_preconditioner = std::make_shared<TrilinosWrappers::PreconditionAMG>();
+    std::vector<std::vector<bool>> constant_modes;
+    FEValuesExtractors::Vector     velocity_components(0);
+    DoFTools::extract_constant_modes(stokes_dof_handler,
+                                     stokes_fe.component_mask(
+                                       velocity_components),
+                                     constant_modes);
+    TrilinosWrappers::PreconditionAMG::AdditionalData amg_data;
+    amg_data.constant_modes        = constant_modes;
+    amg_data.elliptic              = true;
+    amg_data.higher_order_elements = true;
+    amg_data.smoother_sweeps       = 2;
+    amg_data.aggregation_threshold = 0.02;
+    Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0, 0),
+                                   amg_data);
+    Mp_preconditioner = std::make_shared<TrilinosWrappers::PreconditionIC>();
+    Mp_preconditioner->initialize(stokes_preconditioner_matrix.block(1, 1));
+    deallog << std::endl;
+    rebuild_stokes_preconditioner = false;
+  }
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::assemble_stokes_system()
+  {
+    deallog << "   Assembling..." << std::flush;
+    if (rebuild_stokes_matrix == true)
+      stokes_matrix = 0;
+    stokes_rhs                       = 0;
+    const auto &  quadrature_formula = quadrature_stokes;
+    FEValues<dim> stokes_fe_values(
+      mapping,
+      stokes_fe,
+      quadrature_formula,
+      update_values | update_quadrature_points | update_JxW_values |
+        (rebuild_stokes_matrix == true ? update_gradients : UpdateFlags(0)));
+    FEValues<dim>      temperature_fe_values(mapping,
+                                        temperature_fe,
+                                        quadrature_formula,
+                                        update_values);
+    const unsigned int dofs_per_cell = stokes_fe.n_dofs_per_cell();
+    const unsigned int n_q_points    = quadrature_formula.size();
+    FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
+    Vector<double>     local_rhs(dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+    std::vector<double>                  old_temperature_values(n_q_points);
+    std::vector<Tensor<1, dim>>          phi_u(dofs_per_cell);
+    std::vector<SymmetricTensor<2, dim>> grads_phi_u(dofs_per_cell);
+    std::vector<double>                  div_phi_u(dofs_per_cell);
+    std::vector<double>                  phi_p(dofs_per_cell);
+    const FEValuesExtractors::Vector     velocities(0);
+    const FEValuesExtractors::Scalar     pressure(dim);
+    auto       cell             = stokes_dof_handler.begin_active();
+    const auto endc             = stokes_dof_handler.end();
+    auto       temperature_cell = temperature_dof_handler.begin_active();
+    for (; cell != endc; ++cell, ++temperature_cell)
+      {
+        stokes_fe_values.reinit(cell);
+        temperature_fe_values.reinit(temperature_cell);
+        local_matrix = 0;
+        local_rhs    = 0;
+        temperature_fe_values.get_function_values(old_temperature_solution,
+                                                  old_temperature_values);
+        for (unsigned int q = 0; q < n_q_points; ++q)
+          {
+            const double old_temperature = old_temperature_values[q];
+            for (unsigned int k = 0; k < dofs_per_cell; ++k)
+              {
+                phi_u[k] = stokes_fe_values[velocities].value(k, q);
+                if (rebuild_stokes_matrix)
+                  {
+                    grads_phi_u[k] =
+                      stokes_fe_values[velocities].symmetric_gradient(k, q);
+                    div_phi_u[k] =
+                      stokes_fe_values[velocities].divergence(k, q);
+                    phi_p[k] = stokes_fe_values[pressure].value(k, q);
+                  }
+              }
+            if (rebuild_stokes_matrix)
+              for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                  local_matrix(i, j) +=
+                    (EquationData::eta * 2 * (grads_phi_u[i] * grads_phi_u[j]) -
+                     div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) *
+                    stokes_fe_values.JxW(q);
+            const Point<dim> gravity =
+              -((dim == 2) ? (Point<dim>(0, 1)) : (Point<dim>(0, 0, 1)));
+            for (unsigned int i = 0; i < dofs_per_cell; ++i)
+              local_rhs(i) += (-EquationData::density * EquationData::beta *
+                               gravity * phi_u[i] * old_temperature) *
+                              stokes_fe_values.JxW(q);
+          }
+        cell->get_dof_indices(local_dof_indices);
+        if (rebuild_stokes_matrix == true)
+          stokes_constraints.distribute_local_to_global(local_matrix,
+                                                        local_rhs,
+                                                        local_dof_indices,
+                                                        stokes_matrix,
+                                                        stokes_rhs);
+        else
+          stokes_constraints.distribute_local_to_global(local_rhs,
+                                                        local_dof_indices,
+                                                        stokes_rhs);
+      }
+    rebuild_stokes_matrix = false;
+    deallog << std::endl;
+  }
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::assemble_temperature_matrix()
+  {
+    if (rebuild_temperature_matrices == false)
+      return;
+    temperature_mass_matrix               = 0;
+    temperature_stiffness_matrix          = 0;
+    const auto &       quadrature_formula = quadrature_temperature;
+    FEValues<dim>      temperature_fe_values(mapping,
+                                        temperature_fe,
+                                        quadrature_formula,
+                                        update_values | update_gradients |
+                                          update_JxW_values);
+    const unsigned int dofs_per_cell = temperature_fe.n_dofs_per_cell();
+    const unsigned int n_q_points    = quadrature_formula.size();
+    FullMatrix<double> local_mass_matrix(dofs_per_cell, dofs_per_cell);
+    FullMatrix<double> local_stiffness_matrix(dofs_per_cell, dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+    std::vector<double>                  phi_T(dofs_per_cell);
+    std::vector<Tensor<1, dim>>          grad_phi_T(dofs_per_cell);
+    for (const auto &cell : temperature_dof_handler.active_cell_iterators())
+      {
+        local_mass_matrix      = 0;
+        local_stiffness_matrix = 0;
+        temperature_fe_values.reinit(cell);
+        for (unsigned int q = 0; q < n_q_points; ++q)
+          {
+            for (unsigned int k = 0; k < dofs_per_cell; ++k)
+              {
+                grad_phi_T[k] = temperature_fe_values.shape_grad(k, q);
+                phi_T[k]      = temperature_fe_values.shape_value(k, q);
+              }
+            for (unsigned int i = 0; i < dofs_per_cell; ++i)
+              for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                {
+                  local_mass_matrix(i, j) +=
+                    (phi_T[i] * phi_T[j] * temperature_fe_values.JxW(q));
+                  local_stiffness_matrix(i, j) +=
+                    (EquationData::kappa * grad_phi_T[i] * grad_phi_T[j] *
+                     temperature_fe_values.JxW(q));
+                }
+          }
+        cell->get_dof_indices(local_dof_indices);
+        temperature_constraints.distribute_local_to_global(
+          local_mass_matrix, local_dof_indices, temperature_mass_matrix);
+        temperature_constraints.distribute_local_to_global(
+          local_stiffness_matrix,
+          local_dof_indices,
+          temperature_stiffness_matrix);
+      }
+    rebuild_temperature_matrices = false;
+  }
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::assemble_temperature_system(
+    const double maximal_velocity)
+  {
+    const bool use_bdf2_scheme = (timestep_number != 0);
+    if (use_bdf2_scheme == true)
+      {
+        temperature_matrix.copy_from(temperature_mass_matrix);
+        temperature_matrix *=
+          (2 * time_step + old_time_step) / (time_step + old_time_step);
+        temperature_matrix.add(time_step, temperature_stiffness_matrix);
+      }
+    else
+      {
+        temperature_matrix.copy_from(temperature_mass_matrix);
+        temperature_matrix.add(time_step, temperature_stiffness_matrix);
+      }
+    temperature_rhs                       = 0;
+    const auto &       quadrature_formula = quadrature_temperature;
+    FEValues<dim>      temperature_fe_values(mapping,
+                                        temperature_fe,
+                                        quadrature_formula,
+                                        update_values | update_gradients |
+                                          update_hessians |
+                                          update_quadrature_points |
+                                          update_JxW_values);
+    FEValues<dim>      stokes_fe_values(stokes_fe,
+                                   quadrature_formula,
+                                   update_values);
+    const unsigned int dofs_per_cell = temperature_fe.n_dofs_per_cell();
+    const unsigned int n_q_points    = quadrature_formula.size();
+    Vector<double>     local_rhs(dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+    std::vector<Tensor<1, dim>>          old_velocity_values(n_q_points);
+    std::vector<Tensor<1, dim>>          old_old_velocity_values(n_q_points);
+    std::vector<double>                  old_temperature_values(n_q_points);
+    std::vector<double>                  old_old_temperature_values(n_q_points);
+    std::vector<Tensor<1, dim>>          old_temperature_grads(n_q_points);
+    std::vector<Tensor<1, dim>>          old_old_temperature_grads(n_q_points);
+    std::vector<double>                  old_temperature_laplacians(n_q_points);
+    std::vector<double> old_old_temperature_laplacians(n_q_points);
+    EquationData::TemperatureRightHandSide<dim> temperature_right_hand_side;
+    std::vector<double>                         gamma_values(n_q_points);
+    std::vector<double>                         phi_T(dofs_per_cell);
+    std::vector<Tensor<1, dim>>                 grad_phi_T(dofs_per_cell);
+    const std::pair<double, double>             global_T_range =
+      get_extrapolated_temperature_range();
+    const FEValuesExtractors::Vector velocities(0);
+    auto       cell        = temperature_dof_handler.begin_active();
+    const auto endc        = temperature_dof_handler.end();
+    auto       stokes_cell = stokes_dof_handler.begin_active();
+    for (; cell != endc; ++cell, ++stokes_cell)
+      {
+        local_rhs = 0;
+        temperature_fe_values.reinit(cell);
+        stokes_fe_values.reinit(stokes_cell);
+        temperature_fe_values.get_function_values(old_temperature_solution,
+                                                  old_temperature_values);
+        temperature_fe_values.get_function_values(old_old_temperature_solution,
+                                                  old_old_temperature_values);
+        temperature_fe_values.get_function_gradients(old_temperature_solution,
+                                                     old_temperature_grads);
+        temperature_fe_values.get_function_gradients(
+          old_old_temperature_solution, old_old_temperature_grads);
+        temperature_fe_values.get_function_laplacians(
+          old_temperature_solution, old_temperature_laplacians);
+        temperature_fe_values.get_function_laplacians(
+          old_old_temperature_solution, old_old_temperature_laplacians);
+        temperature_right_hand_side.value_list(
+          temperature_fe_values.get_quadrature_points(), gamma_values);
+        stokes_fe_values[velocities].get_function_values(stokes_solution,
+                                                         old_velocity_values);
+        stokes_fe_values[velocities].get_function_values(
+          old_stokes_solution, old_old_velocity_values);
+        const double nu =
+          compute_viscosity(old_temperature_values,
+                            old_old_temperature_values,
+                            old_temperature_grads,
+                            old_old_temperature_grads,
+                            old_temperature_laplacians,
+                            old_old_temperature_laplacians,
+                            old_velocity_values,
+                            old_old_velocity_values,
+                            gamma_values,
+                            maximal_velocity,
+                            global_T_range.second - global_T_range.first,
+                            cell->diameter());
+        for (unsigned int q = 0; q < n_q_points; ++q)
+          {
+            for (unsigned int k = 0; k < dofs_per_cell; ++k)
+              {
+                grad_phi_T[k] = temperature_fe_values.shape_grad(k, q);
+                phi_T[k]      = temperature_fe_values.shape_value(k, q);
+              }
+            const double T_term_for_rhs =
+              (use_bdf2_scheme ?
+                 (old_temperature_values[q] * (1 + time_step / old_time_step) -
+                  old_old_temperature_values[q] * (time_step * time_step) /
+                    (old_time_step * (time_step + old_time_step))) :
+                 old_temperature_values[q]);
+            const Tensor<1, dim> ext_grad_T =
+              (use_bdf2_scheme ?
+                 (old_temperature_grads[q] * (1 + time_step / old_time_step) -
+                  old_old_temperature_grads[q] * time_step / old_time_step) :
+                 old_temperature_grads[q]);
+            const Tensor<1, dim> extrapolated_u =
+              (use_bdf2_scheme ?
+                 (old_velocity_values[q] * (1 + time_step / old_time_step) -
+                  old_old_velocity_values[q] * time_step / old_time_step) :
+                 old_velocity_values[q]);
+            for (unsigned int i = 0; i < dofs_per_cell; ++i)
+              local_rhs(i) +=
+                (T_term_for_rhs * phi_T[i] -
+                 time_step * extrapolated_u * ext_grad_T * phi_T[i] -
+                 time_step * nu * ext_grad_T * grad_phi_T[i] +
+                 time_step * gamma_values[q] * phi_T[i]) *
+                temperature_fe_values.JxW(q);
+          }
+        cell->get_dof_indices(local_dof_indices);
+        temperature_constraints.distribute_local_to_global(local_rhs,
+                                                           local_dof_indices,
+                                                           temperature_rhs);
+      }
+  }
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::solve()
+  {
+    deallog << "   Solving..." << std::endl;
+    {
+      const LinearSolvers::InverseMatrix<TrilinosWrappers::SparseMatrix,
+                                         TrilinosWrappers::PreconditionIC>
+        mp_inverse(stokes_preconditioner_matrix.block(1, 1),
+                   *Mp_preconditioner);
+      const LinearSolvers::BlockSchurPreconditioner<
+        TrilinosWrappers::PreconditionAMG,
+        TrilinosWrappers::PreconditionIC>
+                    preconditioner(stokes_matrix, mp_inverse, *Amg_preconditioner);
+      SolverControl solver_control(stokes_matrix.m(),
+                                   1e-6 * stokes_rhs.l2_norm());
+      SolverGMRES<TrilinosWrappers::MPI::BlockVector> gmres(
+        solver_control,
+        SolverGMRES<TrilinosWrappers::MPI::BlockVector>::AdditionalData(100));
+      for (unsigned int i = 0; i < stokes_solution.size(); ++i)
+        if (stokes_constraints.is_constrained(i))
+          stokes_solution(i) = 0;
+      gmres.solve(stokes_matrix, stokes_solution, stokes_rhs, preconditioner);
+      stokes_constraints.distribute(stokes_solution);
+      deallog << "   " << solver_control.last_step()
+              << " GMRES iterations for Stokes subsystem." << std::endl;
+    }
+    old_time_step                 = time_step;
+    const double maximal_velocity = get_maximal_velocity();
+    if (maximal_velocity >= 0.01)
+      time_step = 1. / (1.7 * dim * std::sqrt(1. * dim)) / temperature_degree *
+                  GridTools::minimal_cell_diameter(triangulation) /
+                  maximal_velocity;
+    else
+      time_step = 1. / (1.7 * dim * std::sqrt(1. * dim)) / temperature_degree *
+                  GridTools::minimal_cell_diameter(triangulation) / .01;
+    deallog << "   "
+            << "Time step: " << time_step << std::endl;
+    temperature_solution = old_temperature_solution;
+    assemble_temperature_system(maximal_velocity);
+    {
+      SolverControl solver_control(temperature_matrix.m(),
+                                   1e-8 * temperature_rhs.l2_norm());
+      SolverCG<TrilinosWrappers::MPI::Vector> cg(solver_control);
+      TrilinosWrappers::PreconditionIC        preconditioner;
+      preconditioner.initialize(temperature_matrix);
+      cg.solve(temperature_matrix,
+               temperature_solution,
+               temperature_rhs,
+               preconditioner);
+      temperature_constraints.distribute(temperature_solution);
+      deallog << "   " << solver_control.last_step()
+              << " CG iterations for temperature." << std::endl;
+      double min_temperature = temperature_solution(0),
+             max_temperature = temperature_solution(0);
+      for (unsigned int i = 0; i < temperature_solution.size(); ++i)
+        {
+          min_temperature =
+            std::min<double>(min_temperature, temperature_solution(i));
+          max_temperature =
+            std::max<double>(max_temperature, temperature_solution(i));
+        }
+      deallog << "   Temperature range: " << min_temperature << ' '
+              << max_temperature << std::endl;
+    }
+  }
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::output_results() const
+  {
+    if (timestep_number % 10 != 0)
+      return;
+    std::vector<std::string> stokes_names(dim, "velocity");
+    stokes_names.emplace_back("p");
+    std::vector<DataComponentInterpretation::DataComponentInterpretation>
+      stokes_component_interpretation(
+        dim + 1, DataComponentInterpretation::component_is_scalar);
+    for (unsigned int i = 0; i < dim; ++i)
+      stokes_component_interpretation[i] =
+        DataComponentInterpretation::component_is_part_of_vector;
+    DataOut<dim> data_out;
+    data_out.add_data_vector(stokes_dof_handler,
+                             stokes_solution,
+                             stokes_names,
+                             stokes_component_interpretation);
+    data_out.add_data_vector(temperature_dof_handler,
+                             temperature_solution,
+                             "T");
+    data_out.build_patches(mapping,
+                           std::max(stokes_degree, temperature_degree));
+    std::ofstream output("solution-" +
+                         Utilities::int_to_string(timestep_number, 4) + ".vtk");
+    data_out.write_vtk(output);
+  }
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::refine_mesh(const unsigned int max_grid_level)
+  {
+    return;
+
+    Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
+    KellyErrorEstimator<dim>::estimate(temperature_dof_handler,
+                                       QGauss<dim - 1>(temperature_degree + 1),
+                                       {},
+                                       temperature_solution,
+                                       estimated_error_per_cell);
+    GridRefinement::refine_and_coarsen_fixed_fraction(triangulation,
+                                                      estimated_error_per_cell,
+                                                      0.8,
+                                                      0.1);
+    if (triangulation.n_levels() > max_grid_level)
+      for (auto &cell :
+           triangulation.active_cell_iterators_on_level(max_grid_level))
+        cell->clear_refine_flag();
+    std::vector<TrilinosWrappers::MPI::Vector> x_temperature(2);
+    x_temperature[0]                            = temperature_solution;
+    x_temperature[1]                            = old_temperature_solution;
+    TrilinosWrappers::MPI::BlockVector x_stokes = stokes_solution;
+    SolutionTransfer<dim, TrilinosWrappers::MPI::Vector> temperature_trans(
+      temperature_dof_handler);
+    SolutionTransfer<dim, TrilinosWrappers::MPI::BlockVector> stokes_trans(
+      stokes_dof_handler);
+    triangulation.prepare_coarsening_and_refinement();
+    temperature_trans.prepare_for_coarsening_and_refinement(x_temperature);
+    stokes_trans.prepare_for_coarsening_and_refinement(x_stokes);
+    triangulation.execute_coarsening_and_refinement();
+    setup_dofs();
+    std::vector<TrilinosWrappers::MPI::Vector> tmp(2);
+    tmp[0].reinit(temperature_solution);
+    tmp[1].reinit(temperature_solution);
+    temperature_trans.interpolate(x_temperature, tmp);
+    temperature_solution     = tmp[0];
+    old_temperature_solution = tmp[1];
+    temperature_constraints.distribute(temperature_solution);
+    temperature_constraints.distribute(old_temperature_solution);
+    stokes_trans.interpolate(x_stokes, stokes_solution);
+    stokes_constraints.distribute(stokes_solution);
+    rebuild_stokes_matrix         = true;
+    rebuild_temperature_matrices  = true;
+    rebuild_stokes_preconditioner = true;
+  }
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::run()
+  {
+    const unsigned int initial_refinement     = (dim == 2 ? 4 : 2);
+    const unsigned int n_pre_refinement_steps = (dim == 2 ? 0 : 0);
+
+    Triangulation<dim> tria_temp;
+    GridGenerator::hyper_cube(tria_temp);
+
+    GridGenerator::convert_hypercube_to_simplex_mesh(tria_temp, triangulation);
+
+    global_Omega_diameter = GridTools::diameter(triangulation);
+    triangulation.refine_global(initial_refinement);
+    setup_dofs();
+    unsigned int pre_refinement_step = 0;
+  start_time_iteration:
+    VectorTools::project(mapping,
+                         temperature_dof_handler,
+                         temperature_constraints,
+                         quadrature_temperature,
+                         EquationData::TemperatureInitialValues<dim>(),
+                         old_temperature_solution);
+    timestep_number = 0;
+    time_step = old_time_step = 0;
+    double time               = 0;
+    do
+      {
+        deallog << "Timestep " << timestep_number << ":  t=" << time
+                << std::endl;
+        assemble_stokes_system();
+        build_stokes_preconditioner();
+        assemble_temperature_matrix();
+        solve();
+        output_results();
+        deallog << std::endl;
+        if ((timestep_number == 0) &&
+            (pre_refinement_step < n_pre_refinement_steps))
+          {
+            refine_mesh(initial_refinement + n_pre_refinement_steps);
+            ++pre_refinement_step;
+            goto start_time_iteration;
+          }
+        else if ((timestep_number > 0) && (timestep_number % 5 == 0))
+          refine_mesh(initial_refinement + n_pre_refinement_steps);
+        time += time_step;
+        ++timestep_number;
+        old_stokes_solution          = stokes_solution;
+        old_old_temperature_solution = old_temperature_solution;
+        old_temperature_solution     = temperature_solution;
+      }
+    while (time <= 2);
+  }
+} // namespace Step31
+int
+main(int argc, char *argv[])
+{
+  try
+    {
+      using namespace dealii;
+      using namespace Step31;
+      Utilities::MPI::MPI_InitFinalize mpi_initialization(
+        argc, argv, numbers::invalid_unsigned_int);
+      initlog();
+      deallog.depth_file(1);
+
+      AssertThrow(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1,
+                  ExcMessage(
+                    "This program can only be run in serial, use ./step-31"));
+      BoussinesqFlowProblem<2> flow_problem;
+      flow_problem.run();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  return 0;
+}
diff --git a/tests/simplex/step-31.with_simplex_support=on.with_trilinos=true.output b/tests/simplex/step-31.with_simplex_support=on.with_trilinos=true.output
new file mode 100644 (file)
index 0000000..07eb3ff
--- /dev/null
@@ -0,0 +1,45 @@
+
+DEAL::Number of active cells: 2048 (on 5 levels)
+DEAL::Number of degrees of freedom: 13764 (8450+1089+4225)
+DEAL::
+DEAL::Timestep 0:  t=0.00000
+DEAL::   Assembling...
+DEAL::   Rebuilding Stokes preconditioner...
+DEAL::   Solving...
+DEAL::   0 GMRES iterations for Stokes subsystem.
+DEAL::   Time step: 0.459559
+DEAL::   9 CG iterations for temperature.
+DEAL::   Temperature range: -0.141015 0.586975
+DEAL::
+DEAL::Timestep 1:  t=0.459559
+DEAL::   Assembling...
+DEAL::   Solving...
+DEAL::   24 GMRES iterations for Stokes subsystem.
+DEAL::   Time step: 0.459559
+DEAL::   9 CG iterations for temperature.
+DEAL::   Temperature range: -0.301941 1.17152
+DEAL::
+DEAL::Timestep 2:  t=0.919118
+DEAL::   Assembling...
+DEAL::   Solving...
+DEAL::   23 GMRES iterations for Stokes subsystem.
+DEAL::   Time step: 0.459559
+DEAL::   9 CG iterations for temperature.
+DEAL::   Temperature range: -0.467381 1.71946
+DEAL::
+DEAL::Timestep 3:  t=1.37868
+DEAL::   Assembling...
+DEAL::   Solving...
+DEAL::   23 GMRES iterations for Stokes subsystem.
+DEAL::   Time step: 0.459559
+DEAL::   9 CG iterations for temperature.
+DEAL::   Temperature range: -0.612772 2.18114
+DEAL::
+DEAL::Timestep 4:  t=1.83824
+DEAL::   Assembling...
+DEAL::   Solving...
+DEAL::   22 GMRES iterations for Stokes subsystem.
+DEAL::   Time step: 0.459559
+DEAL::   9 CG iterations for temperature.
+DEAL::   Temperature range: -0.705842 2.76592
+DEAL::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.