]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add tests for inverse_operator() using Trilinos payload
authorJean-Paul Pelteret <jppelteret@gmail.com>
Sun, 12 Sep 2021 18:01:31 +0000 (20:01 +0200)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Tue, 14 Sep 2021 19:32:34 +0000 (21:32 +0200)
12 files changed:
tests/lac/step-40-linear_operator_01.cc [new file with mode: 0644]
tests/lac/step-40-linear_operator_01.with_trilinos=true.mpirun=3.with_p4est=true.output [new file with mode: 0644]
tests/lac/step-40-linear_operator_02.cc [new file with mode: 0644]
tests/lac/step-40-linear_operator_02.with_trilinos=true.mpirun=3.with_p4est=true.output [new file with mode: 0644]
tests/lac/step-40-linear_operator_03.cc [new file with mode: 0644]
tests/lac/step-40-linear_operator_03.with_trilinos=true.mpirun=3.with_p4est=true.output [new file with mode: 0644]
tests/lac/step-40-linear_operator_04.cc [new file with mode: 0644]
tests/lac/step-40-linear_operator_04.with_trilinos=true.mpirun=3.with_p4est=true.output [new file with mode: 0644]
tests/lac/step-40-linear_operator_05.cc [new file with mode: 0644]
tests/lac/step-40-linear_operator_05.with_trilinos=true.mpirun=3.with_p4est=true.output [new file with mode: 0644]
tests/lac/step-40-linear_operator_06.cc [new file with mode: 0644]
tests/lac/step-40-linear_operator_06.with_trilinos=true.mpirun=3.with_p4est=true.output [new file with mode: 0644]

diff --git a/tests/lac/step-40-linear_operator_01.cc b/tests/lac/step-40-linear_operator_01.cc
new file mode 100644 (file)
index 0000000..7897131
--- /dev/null
@@ -0,0 +1,373 @@
+
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2009 - 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.
+//
+// ---------------------------------------------------------------------
+
+// A lightly adapted version of the step-40 tutorial program.
+// Tests different ways of settings up a Trilinos inverse_operator(),
+// specifically targeting the correct configuration of the
+// TrilinosPayload::inverse_payload() when using different exemplar matrices.
+
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.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/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/linear_operator_tools.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+namespace Step40
+{
+  using namespace dealii;
+
+  template <int dim>
+  class LaplaceProblem
+  {
+  public:
+    LaplaceProblem();
+
+    void
+    run();
+
+  private:
+    void
+    setup_system();
+    void
+    assemble_system();
+    void
+    solve();
+    void
+    refine_grid();
+
+    MPI_Comm mpi_communicator;
+
+    parallel::distributed::Triangulation<dim> triangulation;
+
+    FE_Q<dim>       fe;
+    DoFHandler<dim> dof_handler;
+
+    IndexSet locally_owned_dofs;
+    IndexSet locally_relevant_dofs;
+
+    AffineConstraints<double> constraints;
+
+    TrilinosWrappers::SparseMatrix system_matrix;
+    TrilinosWrappers::MPI::Vector  locally_relevant_solution;
+    TrilinosWrappers::MPI::Vector  system_rhs;
+
+    ConditionalOStream pcout;
+  };
+
+  template <int dim>
+  LaplaceProblem<dim>::LaplaceProblem()
+    : mpi_communicator(MPI_COMM_WORLD)
+    , triangulation(mpi_communicator,
+                    typename Triangulation<dim>::MeshSmoothing(
+                      Triangulation<dim>::smoothing_on_refinement |
+                      Triangulation<dim>::smoothing_on_coarsening))
+    , fe(2)
+    , dof_handler(triangulation)
+    , pcout(Utilities::MPI::this_mpi_process(mpi_communicator) == 0 ?
+              deallog.get_file_stream() :
+              std::cout,
+            (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
+  {}
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+
+    locally_owned_dofs = dof_handler.locally_owned_dofs();
+    DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+    locally_relevant_solution.reinit(locally_owned_dofs,
+                                     locally_relevant_dofs,
+                                     mpi_communicator);
+    system_rhs.reinit(locally_owned_dofs, mpi_communicator);
+
+    constraints.clear();
+    constraints.reinit(locally_relevant_dofs);
+    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             0,
+                                             Functions::ZeroFunction<dim>(),
+                                             constraints);
+    constraints.close();
+
+    DynamicSparsityPattern dsp(locally_relevant_dofs);
+
+    DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
+    SparsityTools::distribute_sparsity_pattern(dsp,
+                                               dof_handler.locally_owned_dofs(),
+                                               mpi_communicator,
+                                               locally_relevant_dofs);
+
+    system_matrix.reinit(locally_owned_dofs,
+                         locally_owned_dofs,
+                         dsp,
+                         mpi_communicator);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::assemble_system()
+  {
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
+
+    FEValues<dim> fe_values(fe,
+                            quadrature_formula,
+                            update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values);
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int n_q_points    = quadrature_formula.size();
+
+    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+    Vector<double>     cell_rhs(dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      if (cell->is_locally_owned())
+        {
+          cell_matrix = 0.;
+          cell_rhs    = 0.;
+
+          fe_values.reinit(cell);
+
+          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+            {
+              const double rhs_value =
+                (fe_values.quadrature_point(q_point)[1] >
+                     0.5 +
+                       0.25 * std::sin(4.0 * numbers::PI *
+                                       fe_values.quadrature_point(q_point)[0]) ?
+                   1. :
+                   -1.);
+
+              for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                {
+                  for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                    cell_matrix(i, j) += fe_values.shape_grad(i, q_point) *
+                                         fe_values.shape_grad(j, q_point) *
+                                         fe_values.JxW(q_point);
+
+                  cell_rhs(i) += rhs_value *                         //
+                                 fe_values.shape_value(i, q_point) * //
+                                 fe_values.JxW(q_point);
+                }
+            }
+
+          cell->get_dof_indices(local_dof_indices);
+          constraints.distribute_local_to_global(cell_matrix,
+                                                 cell_rhs,
+                                                 local_dof_indices,
+                                                 system_matrix,
+                                                 system_rhs);
+        }
+
+    system_matrix.compress(VectorOperation::add);
+    system_rhs.compress(VectorOperation::add);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::solve()
+  {
+    TrilinosWrappers::MPI::Vector completely_distributed_solution(
+      locally_owned_dofs, mpi_communicator);
+
+    SolverControl solver_control(dof_handler.n_dofs(), 1e-12, false, false);
+
+    TrilinosWrappers::SolverCG solver(solver_control);
+
+    TrilinosWrappers::PreconditionAMG preconditioner;
+
+    TrilinosWrappers::PreconditionAMG::AdditionalData data;
+
+    preconditioner.initialize(system_matrix, data);
+
+    auto K     = linear_operator<TrilinosWrappers::MPI::Vector>(system_matrix);
+    auto K_inv = inverse_operator(K, solver, preconditioner);
+
+    check_solver_within_range(completely_distributed_solution =
+                                K_inv * system_rhs,
+                              solver_control.last_step(),
+                              20,
+                              25);
+
+    pcout << "   Solved in " << solver_control.last_step() << " iterations."
+          << std::endl;
+
+    constraints.distribute(completely_distributed_solution);
+
+    locally_relevant_solution = completely_distributed_solution;
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::refine_grid()
+  {
+    triangulation.refine_global(1);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::run()
+  {
+    const unsigned int n_cycles = 2;
+    for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
+      {
+        pcout << "Cycle " << cycle << ':' << std::endl;
+
+        if (cycle == 0)
+          {
+            GridGenerator::hyper_cube(triangulation);
+            triangulation.refine_global(5);
+          }
+        else
+          refine_grid();
+
+        setup_system();
+
+        pcout << "   Number of active cells:       "
+              << triangulation.n_global_active_cells() << std::endl
+              << "      ";
+        const auto n_locally_owned_active_cells_per_processor =
+          Utilities::MPI::all_gather(
+            triangulation.get_communicator(),
+            triangulation.n_locally_owned_active_cells());
+        for (unsigned int i = 0;
+             i < Utilities::MPI::n_mpi_processes(mpi_communicator);
+             ++i)
+          pcout << n_locally_owned_active_cells_per_processor[i] << '+';
+        pcout << std::endl;
+
+        pcout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
+              << std::endl
+              << "      ";
+        for (unsigned int i = 0;
+             i < Utilities::MPI::n_mpi_processes(mpi_communicator);
+             ++i)
+          pcout << Utilities::MPI::all_gather(
+                     MPI_COMM_WORLD, dof_handler.n_locally_owned_dofs())[i]
+                << '+';
+        pcout << std::endl;
+
+        assemble_system();
+        solve();
+
+        pcout << std::endl;
+      }
+  }
+} // namespace Step40
+
+
+int
+test_mpi()
+{
+  try
+    {
+      using namespace Step40;
+
+
+      {
+        LaplaceProblem<2> laplace_problem_2d;
+        laplace_problem_2d.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;
+}
+
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    {
+      initlog();
+
+      deallog.push("mpi");
+      test_mpi();
+      deallog.pop();
+    }
+  else
+    test_mpi();
+}
diff --git a/tests/lac/step-40-linear_operator_01.with_trilinos=true.mpirun=3.with_p4est=true.output b/tests/lac/step-40-linear_operator_01.with_trilinos=true.mpirun=3.with_p4est=true.output
new file mode 100644 (file)
index 0000000..f1fc100
--- /dev/null
@@ -0,0 +1,17 @@
+
+Cycle 0:
+   Number of active cells:       1024
+      340+344+340+
+   Number of degrees of freedom: 4225
+      1453+1412+1360+
+DEAL:mpi::Solver stopped within 20 - 25 iterations
+   Solved in 21 iterations.
+
+Cycle 1:
+   Number of active cells:       4096
+      1364+1368+1364+
+   Number of degrees of freedom: 16641
+      5645+5540+5456+
+DEAL:mpi::Solver stopped within 20 - 25 iterations
+   Solved in 22 iterations.
+
diff --git a/tests/lac/step-40-linear_operator_02.cc b/tests/lac/step-40-linear_operator_02.cc
new file mode 100644 (file)
index 0000000..b8bec22
--- /dev/null
@@ -0,0 +1,376 @@
+
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2009 - 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.
+//
+// ---------------------------------------------------------------------
+
+// A lightly adapted version of the step-40 tutorial program.
+// Tests different ways of settings up a Trilinos inverse_operator(),
+// specifically targeting the correct configuration of the
+// TrilinosPayload::inverse_payload() when using different exemplar matrices.
+
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.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/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/linear_operator_tools.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+namespace Step40
+{
+  using namespace dealii;
+
+  template <int dim>
+  class LaplaceProblem
+  {
+  public:
+    LaplaceProblem();
+
+    void
+    run();
+
+  private:
+    void
+    setup_system();
+    void
+    assemble_system();
+    void
+    solve();
+    void
+    refine_grid();
+
+    MPI_Comm mpi_communicator;
+
+    parallel::distributed::Triangulation<dim> triangulation;
+
+    FE_Q<dim>       fe;
+    DoFHandler<dim> dof_handler;
+
+    IndexSet locally_owned_dofs;
+    IndexSet locally_relevant_dofs;
+
+    AffineConstraints<double> constraints;
+
+    TrilinosWrappers::SparseMatrix system_matrix;
+    TrilinosWrappers::MPI::Vector  locally_relevant_solution;
+    TrilinosWrappers::MPI::Vector  system_rhs;
+
+    ConditionalOStream pcout;
+  };
+
+  template <int dim>
+  LaplaceProblem<dim>::LaplaceProblem()
+    : mpi_communicator(MPI_COMM_WORLD)
+    , triangulation(mpi_communicator,
+                    typename Triangulation<dim>::MeshSmoothing(
+                      Triangulation<dim>::smoothing_on_refinement |
+                      Triangulation<dim>::smoothing_on_coarsening))
+    , fe(2)
+    , dof_handler(triangulation)
+    , pcout(Utilities::MPI::this_mpi_process(mpi_communicator) == 0 ?
+              deallog.get_file_stream() :
+              std::cout,
+            (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
+  {}
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+
+    locally_owned_dofs = dof_handler.locally_owned_dofs();
+    DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+    locally_relevant_solution.reinit(locally_owned_dofs,
+                                     locally_relevant_dofs,
+                                     mpi_communicator);
+    system_rhs.reinit(locally_owned_dofs, mpi_communicator);
+
+    constraints.clear();
+    constraints.reinit(locally_relevant_dofs);
+    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             0,
+                                             Functions::ZeroFunction<dim>(),
+                                             constraints);
+    constraints.close();
+
+    DynamicSparsityPattern dsp(locally_relevant_dofs);
+
+    DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
+    SparsityTools::distribute_sparsity_pattern(dsp,
+                                               dof_handler.locally_owned_dofs(),
+                                               mpi_communicator,
+                                               locally_relevant_dofs);
+
+    system_matrix.reinit(locally_owned_dofs,
+                         locally_owned_dofs,
+                         dsp,
+                         mpi_communicator);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::assemble_system()
+  {
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
+
+    FEValues<dim> fe_values(fe,
+                            quadrature_formula,
+                            update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values);
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int n_q_points    = quadrature_formula.size();
+
+    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+    Vector<double>     cell_rhs(dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      if (cell->is_locally_owned())
+        {
+          cell_matrix = 0.;
+          cell_rhs    = 0.;
+
+          fe_values.reinit(cell);
+
+          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+            {
+              const double rhs_value =
+                (fe_values.quadrature_point(q_point)[1] >
+                     0.5 +
+                       0.25 * std::sin(4.0 * numbers::PI *
+                                       fe_values.quadrature_point(q_point)[0]) ?
+                   1. :
+                   -1.);
+
+              for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                {
+                  for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                    cell_matrix(i, j) += fe_values.shape_grad(i, q_point) *
+                                         fe_values.shape_grad(j, q_point) *
+                                         fe_values.JxW(q_point);
+
+                  cell_rhs(i) += rhs_value *                         //
+                                 fe_values.shape_value(i, q_point) * //
+                                 fe_values.JxW(q_point);
+                }
+            }
+
+          cell->get_dof_indices(local_dof_indices);
+          constraints.distribute_local_to_global(cell_matrix,
+                                                 cell_rhs,
+                                                 local_dof_indices,
+                                                 system_matrix,
+                                                 system_rhs);
+        }
+
+    system_matrix.compress(VectorOperation::add);
+    system_rhs.compress(VectorOperation::add);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::solve()
+  {
+    TrilinosWrappers::MPI::Vector completely_distributed_solution(
+      locally_owned_dofs, mpi_communicator);
+
+    SolverControl solver_control(dof_handler.n_dofs(), 1e-12, false, false);
+
+    TrilinosWrappers::SolverCG solver(solver_control);
+
+    TrilinosWrappers::PreconditionAMG preconditioner;
+
+    TrilinosWrappers::PreconditionAMG::AdditionalData data;
+
+    preconditioner.initialize(system_matrix, data);
+
+    auto K = linear_operator<TrilinosWrappers::MPI::Vector>(system_matrix);
+    auto prec_K =
+      linear_operator<TrilinosWrappers::MPI::Vector>(system_matrix,
+                                                     preconditioner);
+    auto K_inv = inverse_operator(K, solver, prec_K);
+
+    check_solver_within_range(completely_distributed_solution =
+                                K_inv * system_rhs,
+                              solver_control.last_step(),
+                              20,
+                              25);
+
+    pcout << "   Solved in " << solver_control.last_step() << " iterations."
+          << std::endl;
+
+    constraints.distribute(completely_distributed_solution);
+
+    locally_relevant_solution = completely_distributed_solution;
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::refine_grid()
+  {
+    triangulation.refine_global(1);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::run()
+  {
+    const unsigned int n_cycles = 2;
+    for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
+      {
+        pcout << "Cycle " << cycle << ':' << std::endl;
+
+        if (cycle == 0)
+          {
+            GridGenerator::hyper_cube(triangulation);
+            triangulation.refine_global(5);
+          }
+        else
+          refine_grid();
+
+        setup_system();
+
+        pcout << "   Number of active cells:       "
+              << triangulation.n_global_active_cells() << std::endl
+              << "      ";
+        const auto n_locally_owned_active_cells_per_processor =
+          Utilities::MPI::all_gather(
+            triangulation.get_communicator(),
+            triangulation.n_locally_owned_active_cells());
+        for (unsigned int i = 0;
+             i < Utilities::MPI::n_mpi_processes(mpi_communicator);
+             ++i)
+          pcout << n_locally_owned_active_cells_per_processor[i] << '+';
+        pcout << std::endl;
+
+        pcout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
+              << std::endl
+              << "      ";
+        for (unsigned int i = 0;
+             i < Utilities::MPI::n_mpi_processes(mpi_communicator);
+             ++i)
+          pcout << Utilities::MPI::all_gather(
+                     MPI_COMM_WORLD, dof_handler.n_locally_owned_dofs())[i]
+                << '+';
+        pcout << std::endl;
+
+        assemble_system();
+        solve();
+
+        pcout << std::endl;
+      }
+  }
+} // namespace Step40
+
+
+int
+test_mpi()
+{
+  try
+    {
+      using namespace Step40;
+
+
+      {
+        LaplaceProblem<2> laplace_problem_2d;
+        laplace_problem_2d.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;
+}
+
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    {
+      initlog();
+
+      deallog.push("mpi");
+      test_mpi();
+      deallog.pop();
+    }
+  else
+    test_mpi();
+}
diff --git a/tests/lac/step-40-linear_operator_02.with_trilinos=true.mpirun=3.with_p4est=true.output b/tests/lac/step-40-linear_operator_02.with_trilinos=true.mpirun=3.with_p4est=true.output
new file mode 100644 (file)
index 0000000..f1fc100
--- /dev/null
@@ -0,0 +1,17 @@
+
+Cycle 0:
+   Number of active cells:       1024
+      340+344+340+
+   Number of degrees of freedom: 4225
+      1453+1412+1360+
+DEAL:mpi::Solver stopped within 20 - 25 iterations
+   Solved in 21 iterations.
+
+Cycle 1:
+   Number of active cells:       4096
+      1364+1368+1364+
+   Number of degrees of freedom: 16641
+      5645+5540+5456+
+DEAL:mpi::Solver stopped within 20 - 25 iterations
+   Solved in 22 iterations.
+
diff --git a/tests/lac/step-40-linear_operator_03.cc b/tests/lac/step-40-linear_operator_03.cc
new file mode 100644 (file)
index 0000000..d8438bb
--- /dev/null
@@ -0,0 +1,377 @@
+
+
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2009 - 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.
+//
+// ---------------------------------------------------------------------
+
+// A lightly adapted version of the step-40 tutorial program.
+// Tests different ways of settings up a Trilinos inverse_operator(),
+// specifically targeting the correct configuration of the
+// TrilinosPayload::inverse_payload() when using different exemplar matrices.
+
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.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/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/linear_operator_tools.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+namespace Step40
+{
+  using namespace dealii;
+
+  template <int dim>
+  class LaplaceProblem
+  {
+  public:
+    LaplaceProblem();
+
+    void
+    run();
+
+  private:
+    void
+    setup_system();
+    void
+    assemble_system();
+    void
+    solve();
+    void
+    refine_grid();
+
+    MPI_Comm mpi_communicator;
+
+    parallel::distributed::Triangulation<dim> triangulation;
+
+    FE_Q<dim>       fe;
+    DoFHandler<dim> dof_handler;
+
+    IndexSet locally_owned_dofs;
+    IndexSet locally_relevant_dofs;
+
+    AffineConstraints<double> constraints;
+
+    TrilinosWrappers::SparseMatrix system_matrix;
+    TrilinosWrappers::MPI::Vector  locally_relevant_solution;
+    TrilinosWrappers::MPI::Vector  system_rhs;
+
+    ConditionalOStream pcout;
+  };
+
+  template <int dim>
+  LaplaceProblem<dim>::LaplaceProblem()
+    : mpi_communicator(MPI_COMM_WORLD)
+    , triangulation(mpi_communicator,
+                    typename Triangulation<dim>::MeshSmoothing(
+                      Triangulation<dim>::smoothing_on_refinement |
+                      Triangulation<dim>::smoothing_on_coarsening))
+    , fe(2)
+    , dof_handler(triangulation)
+    , pcout(Utilities::MPI::this_mpi_process(mpi_communicator) == 0 ?
+              deallog.get_file_stream() :
+              std::cout,
+            (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
+  {}
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+
+    locally_owned_dofs = dof_handler.locally_owned_dofs();
+    DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+    locally_relevant_solution.reinit(locally_owned_dofs,
+                                     locally_relevant_dofs,
+                                     mpi_communicator);
+    system_rhs.reinit(locally_owned_dofs, mpi_communicator);
+
+    constraints.clear();
+    constraints.reinit(locally_relevant_dofs);
+    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             0,
+                                             Functions::ZeroFunction<dim>(),
+                                             constraints);
+    constraints.close();
+
+    DynamicSparsityPattern dsp(locally_relevant_dofs);
+
+    DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
+    SparsityTools::distribute_sparsity_pattern(dsp,
+                                               dof_handler.locally_owned_dofs(),
+                                               mpi_communicator,
+                                               locally_relevant_dofs);
+
+    system_matrix.reinit(locally_owned_dofs,
+                         locally_owned_dofs,
+                         dsp,
+                         mpi_communicator);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::assemble_system()
+  {
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
+
+    FEValues<dim> fe_values(fe,
+                            quadrature_formula,
+                            update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values);
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int n_q_points    = quadrature_formula.size();
+
+    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+    Vector<double>     cell_rhs(dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      if (cell->is_locally_owned())
+        {
+          cell_matrix = 0.;
+          cell_rhs    = 0.;
+
+          fe_values.reinit(cell);
+
+          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+            {
+              const double rhs_value =
+                (fe_values.quadrature_point(q_point)[1] >
+                     0.5 +
+                       0.25 * std::sin(4.0 * numbers::PI *
+                                       fe_values.quadrature_point(q_point)[0]) ?
+                   1. :
+                   -1.);
+
+              for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                {
+                  for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                    cell_matrix(i, j) += fe_values.shape_grad(i, q_point) *
+                                         fe_values.shape_grad(j, q_point) *
+                                         fe_values.JxW(q_point);
+
+                  cell_rhs(i) += rhs_value *                         //
+                                 fe_values.shape_value(i, q_point) * //
+                                 fe_values.JxW(q_point);
+                }
+            }
+
+          cell->get_dof_indices(local_dof_indices);
+          constraints.distribute_local_to_global(cell_matrix,
+                                                 cell_rhs,
+                                                 local_dof_indices,
+                                                 system_matrix,
+                                                 system_rhs);
+        }
+
+    system_matrix.compress(VectorOperation::add);
+    system_rhs.compress(VectorOperation::add);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::solve()
+  {
+    TrilinosWrappers::MPI::Vector completely_distributed_solution(
+      locally_owned_dofs, mpi_communicator);
+
+    SolverControl solver_control(dof_handler.n_dofs(), 1e-12, false, false);
+
+    TrilinosWrappers::SolverCG solver(solver_control);
+
+    TrilinosWrappers::PreconditionAMG preconditioner;
+
+    TrilinosWrappers::PreconditionAMG::AdditionalData data;
+
+    preconditioner.initialize(system_matrix, data);
+
+    auto K = linear_operator<TrilinosWrappers::MPI::Vector>(system_matrix);
+    auto prec_K =
+      TrilinosWrappers::linear_operator<TrilinosWrappers::MPI::Vector>(
+        system_matrix, preconditioner);
+    auto K_inv = inverse_operator(K, solver, prec_K);
+
+    check_solver_within_range(completely_distributed_solution =
+                                K_inv * system_rhs,
+                              solver_control.last_step(),
+                              20,
+                              25);
+
+    pcout << "   Solved in " << solver_control.last_step() << " iterations."
+          << std::endl;
+
+    constraints.distribute(completely_distributed_solution);
+
+    locally_relevant_solution = completely_distributed_solution;
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::refine_grid()
+  {
+    triangulation.refine_global(1);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::run()
+  {
+    const unsigned int n_cycles = 2;
+    for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
+      {
+        pcout << "Cycle " << cycle << ':' << std::endl;
+
+        if (cycle == 0)
+          {
+            GridGenerator::hyper_cube(triangulation);
+            triangulation.refine_global(5);
+          }
+        else
+          refine_grid();
+
+        setup_system();
+
+        pcout << "   Number of active cells:       "
+              << triangulation.n_global_active_cells() << std::endl
+              << "      ";
+        const auto n_locally_owned_active_cells_per_processor =
+          Utilities::MPI::all_gather(
+            triangulation.get_communicator(),
+            triangulation.n_locally_owned_active_cells());
+        for (unsigned int i = 0;
+             i < Utilities::MPI::n_mpi_processes(mpi_communicator);
+             ++i)
+          pcout << n_locally_owned_active_cells_per_processor[i] << '+';
+        pcout << std::endl;
+
+        pcout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
+              << std::endl
+              << "      ";
+        for (unsigned int i = 0;
+             i < Utilities::MPI::n_mpi_processes(mpi_communicator);
+             ++i)
+          pcout << Utilities::MPI::all_gather(
+                     MPI_COMM_WORLD, dof_handler.n_locally_owned_dofs())[i]
+                << '+';
+        pcout << std::endl;
+
+        assemble_system();
+        solve();
+
+        pcout << std::endl;
+      }
+  }
+} // namespace Step40
+
+
+int
+test_mpi()
+{
+  try
+    {
+      using namespace Step40;
+
+
+      {
+        LaplaceProblem<2> laplace_problem_2d;
+        laplace_problem_2d.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;
+}
+
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    {
+      initlog();
+
+      deallog.push("mpi");
+      test_mpi();
+      deallog.pop();
+    }
+  else
+    test_mpi();
+}
diff --git a/tests/lac/step-40-linear_operator_03.with_trilinos=true.mpirun=3.with_p4est=true.output b/tests/lac/step-40-linear_operator_03.with_trilinos=true.mpirun=3.with_p4est=true.output
new file mode 100644 (file)
index 0000000..f1fc100
--- /dev/null
@@ -0,0 +1,17 @@
+
+Cycle 0:
+   Number of active cells:       1024
+      340+344+340+
+   Number of degrees of freedom: 4225
+      1453+1412+1360+
+DEAL:mpi::Solver stopped within 20 - 25 iterations
+   Solved in 21 iterations.
+
+Cycle 1:
+   Number of active cells:       4096
+      1364+1368+1364+
+   Number of degrees of freedom: 16641
+      5645+5540+5456+
+DEAL:mpi::Solver stopped within 20 - 25 iterations
+   Solved in 22 iterations.
+
diff --git a/tests/lac/step-40-linear_operator_04.cc b/tests/lac/step-40-linear_operator_04.cc
new file mode 100644 (file)
index 0000000..d7fe515
--- /dev/null
@@ -0,0 +1,377 @@
+
+
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2009 - 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.
+//
+// ---------------------------------------------------------------------
+
+// A lightly adapted version of the step-40 tutorial program.
+// Tests different ways of settings up a Trilinos inverse_operator(),
+// specifically targeting the correct configuration of the
+// TrilinosPayload::inverse_payload() when using different exemplar matrices.
+
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.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/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/linear_operator_tools.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+namespace Step40
+{
+  using namespace dealii;
+
+  template <int dim>
+  class LaplaceProblem
+  {
+  public:
+    LaplaceProblem();
+
+    void
+    run();
+
+  private:
+    void
+    setup_system();
+    void
+    assemble_system();
+    void
+    solve();
+    void
+    refine_grid();
+
+    MPI_Comm mpi_communicator;
+
+    parallel::distributed::Triangulation<dim> triangulation;
+
+    FE_Q<dim>       fe;
+    DoFHandler<dim> dof_handler;
+
+    IndexSet locally_owned_dofs;
+    IndexSet locally_relevant_dofs;
+
+    AffineConstraints<double> constraints;
+
+    TrilinosWrappers::SparseMatrix system_matrix;
+    TrilinosWrappers::MPI::Vector  locally_relevant_solution;
+    TrilinosWrappers::MPI::Vector  system_rhs;
+
+    ConditionalOStream pcout;
+  };
+
+  template <int dim>
+  LaplaceProblem<dim>::LaplaceProblem()
+    : mpi_communicator(MPI_COMM_WORLD)
+    , triangulation(mpi_communicator,
+                    typename Triangulation<dim>::MeshSmoothing(
+                      Triangulation<dim>::smoothing_on_refinement |
+                      Triangulation<dim>::smoothing_on_coarsening))
+    , fe(2)
+    , dof_handler(triangulation)
+    , pcout(Utilities::MPI::this_mpi_process(mpi_communicator) == 0 ?
+              deallog.get_file_stream() :
+              std::cout,
+            (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
+  {}
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+
+    locally_owned_dofs = dof_handler.locally_owned_dofs();
+    DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+    locally_relevant_solution.reinit(locally_owned_dofs,
+                                     locally_relevant_dofs,
+                                     mpi_communicator);
+    system_rhs.reinit(locally_owned_dofs, mpi_communicator);
+
+    constraints.clear();
+    constraints.reinit(locally_relevant_dofs);
+    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             0,
+                                             Functions::ZeroFunction<dim>(),
+                                             constraints);
+    constraints.close();
+
+    DynamicSparsityPattern dsp(locally_relevant_dofs);
+
+    DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
+    SparsityTools::distribute_sparsity_pattern(dsp,
+                                               dof_handler.locally_owned_dofs(),
+                                               mpi_communicator,
+                                               locally_relevant_dofs);
+
+    system_matrix.reinit(locally_owned_dofs,
+                         locally_owned_dofs,
+                         dsp,
+                         mpi_communicator);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::assemble_system()
+  {
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
+
+    FEValues<dim> fe_values(fe,
+                            quadrature_formula,
+                            update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values);
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int n_q_points    = quadrature_formula.size();
+
+    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+    Vector<double>     cell_rhs(dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      if (cell->is_locally_owned())
+        {
+          cell_matrix = 0.;
+          cell_rhs    = 0.;
+
+          fe_values.reinit(cell);
+
+          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+            {
+              const double rhs_value =
+                (fe_values.quadrature_point(q_point)[1] >
+                     0.5 +
+                       0.25 * std::sin(4.0 * numbers::PI *
+                                       fe_values.quadrature_point(q_point)[0]) ?
+                   1. :
+                   -1.);
+
+              for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                {
+                  for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                    cell_matrix(i, j) += fe_values.shape_grad(i, q_point) *
+                                         fe_values.shape_grad(j, q_point) *
+                                         fe_values.JxW(q_point);
+
+                  cell_rhs(i) += rhs_value *                         //
+                                 fe_values.shape_value(i, q_point) * //
+                                 fe_values.JxW(q_point);
+                }
+            }
+
+          cell->get_dof_indices(local_dof_indices);
+          constraints.distribute_local_to_global(cell_matrix,
+                                                 cell_rhs,
+                                                 local_dof_indices,
+                                                 system_matrix,
+                                                 system_rhs);
+        }
+
+    system_matrix.compress(VectorOperation::add);
+    system_rhs.compress(VectorOperation::add);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::solve()
+  {
+    TrilinosWrappers::MPI::Vector completely_distributed_solution(
+      locally_owned_dofs, mpi_communicator);
+
+    SolverControl solver_control(dof_handler.n_dofs(), 1e-12, false, false);
+
+    TrilinosWrappers::SolverCG solver(solver_control);
+
+    TrilinosWrappers::PreconditionAMG preconditioner;
+
+    TrilinosWrappers::PreconditionAMG::AdditionalData data;
+
+    preconditioner.initialize(system_matrix, data);
+
+    auto K = linear_operator<TrilinosWrappers::MPI::Vector>(system_matrix);
+    auto prec_K =
+      TrilinosWrappers::linear_operator<TrilinosWrappers::MPI::Vector>(
+        K, preconditioner);
+    auto K_inv = inverse_operator(K, solver, prec_K);
+
+    check_solver_within_range(completely_distributed_solution =
+                                K_inv * system_rhs,
+                              solver_control.last_step(),
+                              20,
+                              25);
+
+    pcout << "   Solved in " << solver_control.last_step() << " iterations."
+          << std::endl;
+
+    constraints.distribute(completely_distributed_solution);
+
+    locally_relevant_solution = completely_distributed_solution;
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::refine_grid()
+  {
+    triangulation.refine_global(1);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::run()
+  {
+    const unsigned int n_cycles = 2;
+    for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
+      {
+        pcout << "Cycle " << cycle << ':' << std::endl;
+
+        if (cycle == 0)
+          {
+            GridGenerator::hyper_cube(triangulation);
+            triangulation.refine_global(5);
+          }
+        else
+          refine_grid();
+
+        setup_system();
+
+        pcout << "   Number of active cells:       "
+              << triangulation.n_global_active_cells() << std::endl
+              << "      ";
+        const auto n_locally_owned_active_cells_per_processor =
+          Utilities::MPI::all_gather(
+            triangulation.get_communicator(),
+            triangulation.n_locally_owned_active_cells());
+        for (unsigned int i = 0;
+             i < Utilities::MPI::n_mpi_processes(mpi_communicator);
+             ++i)
+          pcout << n_locally_owned_active_cells_per_processor[i] << '+';
+        pcout << std::endl;
+
+        pcout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
+              << std::endl
+              << "      ";
+        for (unsigned int i = 0;
+             i < Utilities::MPI::n_mpi_processes(mpi_communicator);
+             ++i)
+          pcout << Utilities::MPI::all_gather(
+                     MPI_COMM_WORLD, dof_handler.n_locally_owned_dofs())[i]
+                << '+';
+        pcout << std::endl;
+
+        assemble_system();
+        solve();
+
+        pcout << std::endl;
+      }
+  }
+} // namespace Step40
+
+
+int
+test_mpi()
+{
+  try
+    {
+      using namespace Step40;
+
+
+      {
+        LaplaceProblem<2> laplace_problem_2d;
+        laplace_problem_2d.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;
+}
+
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    {
+      initlog();
+
+      deallog.push("mpi");
+      test_mpi();
+      deallog.pop();
+    }
+  else
+    test_mpi();
+}
diff --git a/tests/lac/step-40-linear_operator_04.with_trilinos=true.mpirun=3.with_p4est=true.output b/tests/lac/step-40-linear_operator_04.with_trilinos=true.mpirun=3.with_p4est=true.output
new file mode 100644 (file)
index 0000000..f1fc100
--- /dev/null
@@ -0,0 +1,17 @@
+
+Cycle 0:
+   Number of active cells:       1024
+      340+344+340+
+   Number of degrees of freedom: 4225
+      1453+1412+1360+
+DEAL:mpi::Solver stopped within 20 - 25 iterations
+   Solved in 21 iterations.
+
+Cycle 1:
+   Number of active cells:       4096
+      1364+1368+1364+
+   Number of degrees of freedom: 16641
+      5645+5540+5456+
+DEAL:mpi::Solver stopped within 20 - 25 iterations
+   Solved in 22 iterations.
+
diff --git a/tests/lac/step-40-linear_operator_05.cc b/tests/lac/step-40-linear_operator_05.cc
new file mode 100644 (file)
index 0000000..d752d26
--- /dev/null
@@ -0,0 +1,375 @@
+
+
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2009 - 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.
+//
+// ---------------------------------------------------------------------
+
+// A lightly adapted version of the step-40 tutorial program.
+// Tests different ways of settings up a Trilinos inverse_operator(),
+// specifically targeting the correct configuration of the
+// TrilinosPayload::inverse_payload() when using different exemplar matrices.
+
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.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/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/linear_operator_tools.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+namespace Step40
+{
+  using namespace dealii;
+
+  template <int dim>
+  class LaplaceProblem
+  {
+  public:
+    LaplaceProblem();
+
+    void
+    run();
+
+  private:
+    void
+    setup_system();
+    void
+    assemble_system();
+    void
+    solve();
+    void
+    refine_grid();
+
+    MPI_Comm mpi_communicator;
+
+    parallel::distributed::Triangulation<dim> triangulation;
+
+    FE_Q<dim>       fe;
+    DoFHandler<dim> dof_handler;
+
+    IndexSet locally_owned_dofs;
+    IndexSet locally_relevant_dofs;
+
+    AffineConstraints<double> constraints;
+
+    TrilinosWrappers::SparseMatrix system_matrix;
+    TrilinosWrappers::MPI::Vector  locally_relevant_solution;
+    TrilinosWrappers::MPI::Vector  system_rhs;
+
+    ConditionalOStream pcout;
+  };
+
+  template <int dim>
+  LaplaceProblem<dim>::LaplaceProblem()
+    : mpi_communicator(MPI_COMM_WORLD)
+    , triangulation(mpi_communicator,
+                    typename Triangulation<dim>::MeshSmoothing(
+                      Triangulation<dim>::smoothing_on_refinement |
+                      Triangulation<dim>::smoothing_on_coarsening))
+    , fe(2)
+    , dof_handler(triangulation)
+    , pcout(Utilities::MPI::this_mpi_process(mpi_communicator) == 0 ?
+              deallog.get_file_stream() :
+              std::cout,
+            (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
+  {}
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+
+    locally_owned_dofs = dof_handler.locally_owned_dofs();
+    DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+    locally_relevant_solution.reinit(locally_owned_dofs,
+                                     locally_relevant_dofs,
+                                     mpi_communicator);
+    system_rhs.reinit(locally_owned_dofs, mpi_communicator);
+
+    constraints.clear();
+    constraints.reinit(locally_relevant_dofs);
+    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             0,
+                                             Functions::ZeroFunction<dim>(),
+                                             constraints);
+    constraints.close();
+
+    DynamicSparsityPattern dsp(locally_relevant_dofs);
+
+    DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
+    SparsityTools::distribute_sparsity_pattern(dsp,
+                                               dof_handler.locally_owned_dofs(),
+                                               mpi_communicator,
+                                               locally_relevant_dofs);
+
+    system_matrix.reinit(locally_owned_dofs,
+                         locally_owned_dofs,
+                         dsp,
+                         mpi_communicator);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::assemble_system()
+  {
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
+
+    FEValues<dim> fe_values(fe,
+                            quadrature_formula,
+                            update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values);
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int n_q_points    = quadrature_formula.size();
+
+    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+    Vector<double>     cell_rhs(dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      if (cell->is_locally_owned())
+        {
+          cell_matrix = 0.;
+          cell_rhs    = 0.;
+
+          fe_values.reinit(cell);
+
+          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+            {
+              const double rhs_value =
+                (fe_values.quadrature_point(q_point)[1] >
+                     0.5 +
+                       0.25 * std::sin(4.0 * numbers::PI *
+                                       fe_values.quadrature_point(q_point)[0]) ?
+                   1. :
+                   -1.);
+
+              for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                {
+                  for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                    cell_matrix(i, j) += fe_values.shape_grad(i, q_point) *
+                                         fe_values.shape_grad(j, q_point) *
+                                         fe_values.JxW(q_point);
+
+                  cell_rhs(i) += rhs_value *                         //
+                                 fe_values.shape_value(i, q_point) * //
+                                 fe_values.JxW(q_point);
+                }
+            }
+
+          cell->get_dof_indices(local_dof_indices);
+          constraints.distribute_local_to_global(cell_matrix,
+                                                 cell_rhs,
+                                                 local_dof_indices,
+                                                 system_matrix,
+                                                 system_rhs);
+        }
+
+    system_matrix.compress(VectorOperation::add);
+    system_rhs.compress(VectorOperation::add);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::solve()
+  {
+    TrilinosWrappers::MPI::Vector completely_distributed_solution(
+      locally_owned_dofs, mpi_communicator);
+
+    SolverControl solver_control(dof_handler.n_dofs(), 1e-12, false, false);
+
+    TrilinosWrappers::SolverCG solver(solver_control);
+
+    TrilinosWrappers::PreconditionAMG preconditioner;
+
+    TrilinosWrappers::PreconditionAMG::AdditionalData data;
+
+    preconditioner.initialize(system_matrix, data);
+
+    auto K      = linear_operator<TrilinosWrappers::MPI::Vector>(system_matrix);
+    auto prec_K = TrilinosWrappers::linear_operator(K, preconditioner);
+    auto K_inv  = inverse_operator(K, solver, prec_K);
+
+    check_solver_within_range(completely_distributed_solution =
+                                K_inv * system_rhs,
+                              solver_control.last_step(),
+                              20,
+                              25);
+
+    pcout << "   Solved in " << solver_control.last_step() << " iterations."
+          << std::endl;
+
+    constraints.distribute(completely_distributed_solution);
+
+    locally_relevant_solution = completely_distributed_solution;
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::refine_grid()
+  {
+    triangulation.refine_global(1);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::run()
+  {
+    const unsigned int n_cycles = 2;
+    for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
+      {
+        pcout << "Cycle " << cycle << ':' << std::endl;
+
+        if (cycle == 0)
+          {
+            GridGenerator::hyper_cube(triangulation);
+            triangulation.refine_global(5);
+          }
+        else
+          refine_grid();
+
+        setup_system();
+
+        pcout << "   Number of active cells:       "
+              << triangulation.n_global_active_cells() << std::endl
+              << "      ";
+        const auto n_locally_owned_active_cells_per_processor =
+          Utilities::MPI::all_gather(
+            triangulation.get_communicator(),
+            triangulation.n_locally_owned_active_cells());
+        for (unsigned int i = 0;
+             i < Utilities::MPI::n_mpi_processes(mpi_communicator);
+             ++i)
+          pcout << n_locally_owned_active_cells_per_processor[i] << '+';
+        pcout << std::endl;
+
+        pcout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
+              << std::endl
+              << "      ";
+        for (unsigned int i = 0;
+             i < Utilities::MPI::n_mpi_processes(mpi_communicator);
+             ++i)
+          pcout << Utilities::MPI::all_gather(
+                     MPI_COMM_WORLD, dof_handler.n_locally_owned_dofs())[i]
+                << '+';
+        pcout << std::endl;
+
+        assemble_system();
+        solve();
+
+        pcout << std::endl;
+      }
+  }
+} // namespace Step40
+
+
+int
+test_mpi()
+{
+  try
+    {
+      using namespace Step40;
+
+
+      {
+        LaplaceProblem<2> laplace_problem_2d;
+        laplace_problem_2d.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;
+}
+
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    {
+      initlog();
+
+      deallog.push("mpi");
+      test_mpi();
+      deallog.pop();
+    }
+  else
+    test_mpi();
+}
diff --git a/tests/lac/step-40-linear_operator_05.with_trilinos=true.mpirun=3.with_p4est=true.output b/tests/lac/step-40-linear_operator_05.with_trilinos=true.mpirun=3.with_p4est=true.output
new file mode 100644 (file)
index 0000000..f1fc100
--- /dev/null
@@ -0,0 +1,17 @@
+
+Cycle 0:
+   Number of active cells:       1024
+      340+344+340+
+   Number of degrees of freedom: 4225
+      1453+1412+1360+
+DEAL:mpi::Solver stopped within 20 - 25 iterations
+   Solved in 21 iterations.
+
+Cycle 1:
+   Number of active cells:       4096
+      1364+1368+1364+
+   Number of degrees of freedom: 16641
+      5645+5540+5456+
+DEAL:mpi::Solver stopped within 20 - 25 iterations
+   Solved in 22 iterations.
+
diff --git a/tests/lac/step-40-linear_operator_06.cc b/tests/lac/step-40-linear_operator_06.cc
new file mode 100644 (file)
index 0000000..91f3f25
--- /dev/null
@@ -0,0 +1,375 @@
+
+
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2009 - 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.
+//
+// ---------------------------------------------------------------------
+
+// A lightly adapted version of the step-40 tutorial program.
+// Tests different ways of settings up a Trilinos inverse_operator(),
+// specifically targeting the correct configuration of the
+// TrilinosPayload::inverse_payload() when using different exemplar matrices.
+
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.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/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/linear_operator_tools.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+namespace Step40
+{
+  using namespace dealii;
+
+  template <int dim>
+  class LaplaceProblem
+  {
+  public:
+    LaplaceProblem();
+
+    void
+    run();
+
+  private:
+    void
+    setup_system();
+    void
+    assemble_system();
+    void
+    solve();
+    void
+    refine_grid();
+
+    MPI_Comm mpi_communicator;
+
+    parallel::distributed::Triangulation<dim> triangulation;
+
+    FE_Q<dim>       fe;
+    DoFHandler<dim> dof_handler;
+
+    IndexSet locally_owned_dofs;
+    IndexSet locally_relevant_dofs;
+
+    AffineConstraints<double> constraints;
+
+    TrilinosWrappers::SparseMatrix system_matrix;
+    TrilinosWrappers::MPI::Vector  locally_relevant_solution;
+    TrilinosWrappers::MPI::Vector  system_rhs;
+
+    ConditionalOStream pcout;
+  };
+
+  template <int dim>
+  LaplaceProblem<dim>::LaplaceProblem()
+    : mpi_communicator(MPI_COMM_WORLD)
+    , triangulation(mpi_communicator,
+                    typename Triangulation<dim>::MeshSmoothing(
+                      Triangulation<dim>::smoothing_on_refinement |
+                      Triangulation<dim>::smoothing_on_coarsening))
+    , fe(2)
+    , dof_handler(triangulation)
+    , pcout(Utilities::MPI::this_mpi_process(mpi_communicator) == 0 ?
+              deallog.get_file_stream() :
+              std::cout,
+            (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
+  {}
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+
+    locally_owned_dofs = dof_handler.locally_owned_dofs();
+    DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+    locally_relevant_solution.reinit(locally_owned_dofs,
+                                     locally_relevant_dofs,
+                                     mpi_communicator);
+    system_rhs.reinit(locally_owned_dofs, mpi_communicator);
+
+    constraints.clear();
+    constraints.reinit(locally_relevant_dofs);
+    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             0,
+                                             Functions::ZeroFunction<dim>(),
+                                             constraints);
+    constraints.close();
+
+    DynamicSparsityPattern dsp(locally_relevant_dofs);
+
+    DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
+    SparsityTools::distribute_sparsity_pattern(dsp,
+                                               dof_handler.locally_owned_dofs(),
+                                               mpi_communicator,
+                                               locally_relevant_dofs);
+
+    system_matrix.reinit(locally_owned_dofs,
+                         locally_owned_dofs,
+                         dsp,
+                         mpi_communicator);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::assemble_system()
+  {
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
+
+    FEValues<dim> fe_values(fe,
+                            quadrature_formula,
+                            update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values);
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int n_q_points    = quadrature_formula.size();
+
+    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+    Vector<double>     cell_rhs(dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      if (cell->is_locally_owned())
+        {
+          cell_matrix = 0.;
+          cell_rhs    = 0.;
+
+          fe_values.reinit(cell);
+
+          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+            {
+              const double rhs_value =
+                (fe_values.quadrature_point(q_point)[1] >
+                     0.5 +
+                       0.25 * std::sin(4.0 * numbers::PI *
+                                       fe_values.quadrature_point(q_point)[0]) ?
+                   1. :
+                   -1.);
+
+              for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                {
+                  for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                    cell_matrix(i, j) += fe_values.shape_grad(i, q_point) *
+                                         fe_values.shape_grad(j, q_point) *
+                                         fe_values.JxW(q_point);
+
+                  cell_rhs(i) += rhs_value *                         //
+                                 fe_values.shape_value(i, q_point) * //
+                                 fe_values.JxW(q_point);
+                }
+            }
+
+          cell->get_dof_indices(local_dof_indices);
+          constraints.distribute_local_to_global(cell_matrix,
+                                                 cell_rhs,
+                                                 local_dof_indices,
+                                                 system_matrix,
+                                                 system_rhs);
+        }
+
+    system_matrix.compress(VectorOperation::add);
+    system_rhs.compress(VectorOperation::add);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::solve()
+  {
+    TrilinosWrappers::MPI::Vector completely_distributed_solution(
+      locally_owned_dofs, mpi_communicator);
+
+    SolverControl solver_control(dof_handler.n_dofs(), 1e-12, false, false);
+
+    TrilinosWrappers::SolverCG solver(solver_control);
+
+    TrilinosWrappers::PreconditionAMG preconditioner;
+
+    TrilinosWrappers::PreconditionAMG::AdditionalData data;
+
+    preconditioner.initialize(system_matrix, data);
+
+    auto K      = linear_operator<TrilinosWrappers::MPI::Vector>(system_matrix);
+    auto prec_K = linear_operator(K, preconditioner);
+    auto K_inv  = inverse_operator(K, solver, prec_K);
+
+    check_solver_within_range(completely_distributed_solution =
+                                K_inv * system_rhs,
+                              solver_control.last_step(),
+                              20,
+                              25);
+
+    pcout << "   Solved in " << solver_control.last_step() << " iterations."
+          << std::endl;
+
+    constraints.distribute(completely_distributed_solution);
+
+    locally_relevant_solution = completely_distributed_solution;
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::refine_grid()
+  {
+    triangulation.refine_global(1);
+  }
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::run()
+  {
+    const unsigned int n_cycles = 2;
+    for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
+      {
+        pcout << "Cycle " << cycle << ':' << std::endl;
+
+        if (cycle == 0)
+          {
+            GridGenerator::hyper_cube(triangulation);
+            triangulation.refine_global(5);
+          }
+        else
+          refine_grid();
+
+        setup_system();
+
+        pcout << "   Number of active cells:       "
+              << triangulation.n_global_active_cells() << std::endl
+              << "      ";
+        const auto n_locally_owned_active_cells_per_processor =
+          Utilities::MPI::all_gather(
+            triangulation.get_communicator(),
+            triangulation.n_locally_owned_active_cells());
+        for (unsigned int i = 0;
+             i < Utilities::MPI::n_mpi_processes(mpi_communicator);
+             ++i)
+          pcout << n_locally_owned_active_cells_per_processor[i] << '+';
+        pcout << std::endl;
+
+        pcout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
+              << std::endl
+              << "      ";
+        for (unsigned int i = 0;
+             i < Utilities::MPI::n_mpi_processes(mpi_communicator);
+             ++i)
+          pcout << Utilities::MPI::all_gather(
+                     MPI_COMM_WORLD, dof_handler.n_locally_owned_dofs())[i]
+                << '+';
+        pcout << std::endl;
+
+        assemble_system();
+        solve();
+
+        pcout << std::endl;
+      }
+  }
+} // namespace Step40
+
+
+int
+test_mpi()
+{
+  try
+    {
+      using namespace Step40;
+
+
+      {
+        LaplaceProblem<2> laplace_problem_2d;
+        laplace_problem_2d.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;
+}
+
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    {
+      initlog();
+
+      deallog.push("mpi");
+      test_mpi();
+      deallog.pop();
+    }
+  else
+    test_mpi();
+}
diff --git a/tests/lac/step-40-linear_operator_06.with_trilinos=true.mpirun=3.with_p4est=true.output b/tests/lac/step-40-linear_operator_06.with_trilinos=true.mpirun=3.with_p4est=true.output
new file mode 100644 (file)
index 0000000..f1fc100
--- /dev/null
@@ -0,0 +1,17 @@
+
+Cycle 0:
+   Number of active cells:       1024
+      340+344+340+
+   Number of degrees of freedom: 4225
+      1453+1412+1360+
+DEAL:mpi::Solver stopped within 20 - 25 iterations
+   Solved in 21 iterations.
+
+Cycle 1:
+   Number of active cells:       4096
+      1364+1368+1364+
+   Number of degrees of freedom: 16641
+      5645+5540+5456+
+DEAL:mpi::Solver stopped within 20 - 25 iterations
+   Solved in 22 iterations.
+

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.