]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Commit bug in Trilinos precondition
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 30 Jun 2014 10:23:34 +0000 (10:23 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 30 Jun 2014 10:23:34 +0000 (10:23 +0000)
git-svn-id: https://svn.dealii.org/trunk@33099 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/lac/trilinos_precondition.h
deal.II/source/lac/trilinos_precondition.cc
tests/trilinos/precondition_amg_dgp.cc [new file with mode: 0644]
tests/trilinos/precondition_amg_dgp.with_64bit_indices=off.output [new file with mode: 0644]

index ba57ec5880d9cb2d750f2fefe8d414b91fe9d871..95df70da1b7e22d28ca64ed4ebade26dd296130c 100644 (file)
@@ -148,6 +148,13 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Fixed: TrilinosWrappers::PreconditionAMG did not read user-provided
+  constant modes (aka null space) when the null space dimension is one but not
+  just the trivial one vector. This is now fixed.
+  <br>
+  (Martin Kronbichler, 2014/06/30)
+  </li>
+
   <li> Simplified interfaces for FEEvaluation: Previously, the user had to
   select the appropriate kernel (FEEvaluation, FEEvaluationGeneral,
   FEEvaluationDGP, FEEvaluationGL) for the matrix-free evaluation
index ba1d6ce23ec67dac146352e672d3fcbd6af09abb..9ab30fb98f92292488c6fc5661febcb3f4ca03c7 100644 (file)
@@ -1301,7 +1301,7 @@ namespace TrilinosWrappers
                       const unsigned int                     n_cycles = 1,
                       const bool                             w_cyle = false,
                       const double                           aggregation_threshold = 1e-4,
-                      const std::vector<std::vector<bool> > &constant_modes = std::vector<std::vector<bool> > (1),
+                      const std::vector<std::vector<bool> > &constant_modes = std::vector<std::vector<bool> > (0),
                       const unsigned int                     smoother_sweeps = 2,
                       const unsigned int                     smoother_overlap = 0,
                       const bool                             output_details = false,
index a526b310d76fe851b2fb523f3f5e1b571854f288..4c42b28e7875b6fd3399654481eeb7ed8c617d7c 100644 (file)
@@ -611,7 +611,7 @@ namespace TrilinosWrappers
                                                    constant_modes_dimension);
     std::vector<double> dummy (constant_modes_dimension);
 
-    if (constant_modes_dimension > 1)
+    if (constant_modes_dimension > 0)
       {
         const bool constant_modes_are_global =
           additional_data.constant_modes[0].size() == n_rows;
diff --git a/tests/trilinos/precondition_amg_dgp.cc b/tests/trilinos/precondition_amg_dgp.cc
new file mode 100644 (file)
index 0000000..dacb279
--- /dev/null
@@ -0,0 +1,254 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2013 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// solves a 2D Poisson equation for linear FE_DGP elements (SIP
+// discretization) with AMG preconditioner
+
+#include "../tests.h"
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/base/function.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/meshworker/dof_info.h>
+#include <deal.II/meshworker/integration_info.h>
+#include <deal.II/meshworker/assembler.h>
+#include <deal.II/meshworker/loop.h>
+#include <deal.II/integrators/laplace.h>
+
+
+#include <fstream>
+#include <iomanip>
+
+
+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_Q<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 (7);
+}
+
+
+
+template <int dim>
+void Step4<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (fe);
+
+  CompressedSparsityPattern 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());
+
+  MappingQ1<dim> mapping;
+  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 ()
+{
+
+  // variant 1: solve with AMG
+  deallog.push(Utilities::int_to_string(dof_handler.n_dofs(),5));
+  TrilinosWrappers::PreconditionAMG preconditioner;
+  TrilinosWrappers::PreconditionAMG::AdditionalData data;
+  DoFTools::extract_constant_modes(dof_handler, std::vector<bool>(1,true),
+                                   data.constant_modes);
+  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)
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+  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.with_64bit_indices=off.output b/tests/trilinos/precondition_amg_dgp.with_64bit_indices=off.output
new file mode 100644 (file)
index 0000000..45b1d6c
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL:16641:cg::Starting value 5039.73
+DEAL:16641:cg::Convergence step 15 value 0
+DEAL:66049:cg::Starting value 407234.
+DEAL:66049:cg::Convergence step 17 value 0

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.