From: Jean-Paul Pelteret <jppelteret@gmail.com>
Date: Fri, 15 Mar 2019 08:40:03 +0000 (+0100)
Subject: Remove dummy vector when initializing AMG constant modes
X-Git-Tag: v9.1.0-rc1~252^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F7793%2Fhead;p=dealii.git

Remove dummy vector when initializing AMG constant modes
---

diff --git a/source/lac/trilinos_precondition_ml.cc b/source/lac/trilinos_precondition_ml.cc
index 848315efaa..1e65a430e9 100644
--- a/source/lac/trilinos_precondition_ml.cc
+++ b/source/lac/trilinos_precondition_ml.cc
@@ -178,18 +178,8 @@ namespace TrilinosWrappers
         parameter_list.set("null space: type", "pre-computed");
         parameter_list.set("null space: dimension",
                            distributed_constant_modes.NumVectors());
-        if (my_size > 0)
-          parameter_list.set("null space: vectors",
-                             distributed_constant_modes.Values());
-        else
-          {
-            // We need to set a valid pointer to data even if there is no data
-            // on the current processor. Therefore, pass a dummy in that case
-            static std::vector<double> dummy;
-            if (dummy.size() != constant_modes_dimension)
-              dummy.resize(constant_modes_dimension);
-            parameter_list.set("null space: vectors", dummy.data());
-          }
+        parameter_list.set("null space: vectors",
+                           distributed_constant_modes.Values());
       }
   }
 
diff --git a/tests/trilinos/precondition_amg_dgp_03.cc b/tests/trilinos/precondition_amg_dgp_03.cc
new file mode 100644
index 0000000000..1fcf871403
--- /dev/null
+++ b/tests/trilinos/precondition_amg_dgp_03.cc
@@ -0,0 +1,273 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2013 - 2019 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// solves a 2D Poisson equation for linear FE_DGP elements (SIP
+// discretization) with AMG preconditioner.
+//
+// This variant of the test trilinos/precondition_amg_dgp_01 does not initialize
+// the constant modes vector passes to the AMG preconditioner.
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgp.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/integrators/laplace.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/meshworker/assembler.h>
+#include <deal.II/meshworker/dof_info.h>
+#include <deal.II/meshworker/integration_info.h>
+#include <deal.II/meshworker/loop.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+class MatrixIntegrator : public MeshWorker::LocalIntegrator<dim>
+{
+public:
+  void
+  cell(MeshWorker::DoFInfo<dim> &                 dinfo,
+       typename MeshWorker::IntegrationInfo<dim> &info) const;
+  void
+  boundary(MeshWorker::DoFInfo<dim> &                 dinfo,
+           typename MeshWorker::IntegrationInfo<dim> &info) const;
+  void
+  face(MeshWorker::DoFInfo<dim> &                 dinfo1,
+       MeshWorker::DoFInfo<dim> &                 dinfo2,
+       typename MeshWorker::IntegrationInfo<dim> &info1,
+       typename MeshWorker::IntegrationInfo<dim> &info2) const;
+};
+
+template <int dim>
+void
+MatrixIntegrator<dim>::cell(
+  MeshWorker::DoFInfo<dim> &                 dinfo,
+  typename MeshWorker::IntegrationInfo<dim> &info) const
+{
+  LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0, false).matrix,
+                                         info.fe_values());
+}
+
+template <int dim>
+void
+MatrixIntegrator<dim>::boundary(
+  MeshWorker::DoFInfo<dim> &                 dinfo,
+  typename MeshWorker::IntegrationInfo<dim> &info) const
+{
+  const unsigned int deg = info.fe_values(0).get_fe().degree;
+  LocalIntegrators::Laplace ::nitsche_matrix(
+    dinfo.matrix(0, false).matrix,
+    info.fe_values(0),
+    LocalIntegrators::Laplace::compute_penalty(dinfo, dinfo, deg, deg));
+}
+
+template <int dim>
+void
+MatrixIntegrator<dim>::face(
+  MeshWorker::DoFInfo<dim> &                 dinfo1,
+  MeshWorker::DoFInfo<dim> &                 dinfo2,
+  typename MeshWorker::IntegrationInfo<dim> &info1,
+  typename MeshWorker::IntegrationInfo<dim> &info2) const
+{
+  const unsigned int deg = info1.fe_values(0).get_fe().degree;
+  LocalIntegrators::Laplace ::ip_matrix(
+    dinfo1.matrix(0, false).matrix,
+    dinfo1.matrix(0, true).matrix,
+    dinfo2.matrix(0, true).matrix,
+    dinfo2.matrix(0, false).matrix,
+    info1.fe_values(0),
+    info2.fe_values(0),
+    LocalIntegrators::Laplace::compute_penalty(dinfo1, dinfo2, deg, deg));
+}
+
+
+template <int dim>
+class Step4
+{
+public:
+  Step4();
+  void
+  run();
+
+private:
+  void
+  make_grid();
+  void
+  setup_system();
+  void
+  solve();
+
+  Triangulation<dim> triangulation;
+  FE_DGP<dim>        fe;
+  DoFHandler<dim>    dof_handler;
+
+  TrilinosWrappers::SparseMatrix system_matrix;
+
+  Vector<double> solution;
+  Vector<double> system_rhs;
+};
+
+
+
+template <int dim>
+Step4<dim>::Step4()
+  : fe(1)
+  , dof_handler(triangulation)
+{}
+
+
+template <int dim>
+void
+Step4<dim>::make_grid()
+{
+  GridGenerator::hyper_cube(triangulation, -1, 1);
+  triangulation.refine_global(6);
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::setup_system()
+{
+  dof_handler.distribute_dofs(fe);
+
+  DynamicSparsityPattern c_sparsity(dof_handler.n_dofs());
+  DoFTools::make_flux_sparsity_pattern(dof_handler, c_sparsity);
+  system_matrix.reinit(c_sparsity);
+
+  solution.reinit(dof_handler.n_dofs());
+  system_rhs.reinit(dof_handler.n_dofs());
+
+  MappingQGeneric<dim>                mapping(1);
+  MeshWorker::IntegrationInfoBox<dim> info_box;
+  UpdateFlags update_flags = update_values | update_gradients;
+  info_box.add_update_flags_all(update_flags);
+  info_box.initialize(fe, mapping);
+
+  MeshWorker::DoFInfo<dim> dof_info(dof_handler);
+  MeshWorker::Assembler::MatrixSimple<TrilinosWrappers::SparseMatrix> assembler;
+  assembler.initialize(system_matrix);
+  MatrixIntegrator<dim> integrator;
+  MeshWorker::integration_loop<dim, dim>(dof_handler.begin_active(),
+                                         dof_handler.end(),
+                                         dof_info,
+                                         info_box,
+                                         integrator,
+                                         assembler);
+
+  system_matrix.compress(VectorOperation::add);
+
+  for (unsigned int i = 0; i < system_rhs.size(); ++i)
+    system_rhs(i) = 0.01 * i - 0.000001 * i * i;
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::solve()
+{
+  deallog.push(Utilities::int_to_string(dof_handler.n_dofs(), 5));
+  TrilinosWrappers::PreconditionAMG                 preconditioner;
+  TrilinosWrappers::PreconditionAMG::AdditionalData data;
+  data.smoother_sweeps = 2;
+  {
+    solution = 0;
+    SolverControl solver_control(1000, 1e-10);
+    SolverCG<>    solver(solver_control);
+    preconditioner.initialize(system_matrix, data);
+    solver.solve(system_matrix, solution, system_rhs, preconditioner);
+  }
+  deallog.pop();
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::run()
+{
+  for (unsigned int cycle = 0; cycle < 2; ++cycle)
+    {
+      if (cycle == 0)
+        make_grid();
+      else
+        triangulation.refine_global(1);
+
+      setup_system();
+      solve();
+    }
+}
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, testing_max_num_threads());
+
+  try
+    {
+      Step4<2> test;
+      test.run();
+    }
+  catch (std::exception &exc)
+    {
+      deallog << std::endl
+              << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Exception on processing: " << std::endl
+              << exc.what() << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      deallog << std::endl
+              << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Unknown exception!" << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/trilinos/precondition_amg_dgp_03.with_64bit_indices=off.output b/tests/trilinos/precondition_amg_dgp_03.with_64bit_indices=off.output
new file mode 100644
index 0000000000..5399b9277c
--- /dev/null
+++ b/tests/trilinos/precondition_amg_dgp_03.with_64bit_indices=off.output
@@ -0,0 +1,5 @@
+
+DEAL:12288:cg::Starting value 1970.22
+DEAL:12288:cg::Convergence step 121 value 0
+DEAL:49152:cg::Starting value 179304.
+DEAL:49152:cg::Convergence step 271 value 0