Implement explicit mass lumping in MassOperator.

diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h
index 3608a37ffed1d1365502f83073eeb3f0ab4cf5e0..017af8daa43c65a66ccfd3aa3ef715f9a412e5c5 100644 (file)
--- a/include/deal.II/matrix_free/operators.h
+++ b/include/deal.II/matrix_free/operators.h
@@ -763,12 +763,41 @@ namespace MatrixFreeOperators
     MassOperator();
 
     /**
-     * For preconditioning, we store a lumped mass matrix at the diagonal
-     * entries.
+     * Same as the base class.
      */
     virtual void
     compute_diagonal() override;
 
+    /**
+     * Compute the lumped mass matrix. This is equal to the mass matrix times a
+     * vector of all ones and is equivalent to approximating the mass matrix
+     * with a nodal quadrature rule.
+     *
+     * The lumped mass matrix is an excellent preconditioner for mass matrices
+     * corresponding to FE_Q elements on axis-aligned cells. However, some
+     * elements (like FE_SimplexP with degrees higher than 1) have basis
+     * functions whose integrals are zero or negative (and therefore their
+     * lumped mass matrix entries are zero or negative). For such elements a
+     * lumped mass matrix is a very poor approximation of the operator - the
+     * diagonal should be used instead. If you are interested in using mass
+     * lumping with simplices then use FE_SimplexP_Bubbles instead of
+     * FE_SimplexP.
+     */
+    void
+    compute_lumped_diagonal();
+
+    /**
+     * Get read access to the lumped diagonal of this operator.
+     */
+    const std::shared_ptr<DiagonalMatrix<VectorType>> &
+    get_matrix_lumped_diagonal() const;
+
+    /**
+     * Get read access to the inverse lumped diagonal of this operator.
+     */
+    const std::shared_ptr<DiagonalMatrix<VectorType>> &
+    get_matrix_lumped_diagonal_inverse() const;
+
   private:
     /**
      * Applies the mass matrix operation on an input vector. It is
@@ -787,6 +816,18 @@ namespace MatrixFreeOperators
       VectorType &                                            dst,
       const VectorType &                                      src,
       const std::pair<unsigned int, unsigned int> &           cell_range) const;
+
+    /**
+     * A shared pointer to a diagonal matrix that stores the
+     * lumped diagonal elements as a vector.
+     */
+    std::shared_ptr<DiagonalMatrix<VectorType>> lumped_diagonal_entries;
+
+    /**
+     * A shared pointer to a diagonal matrix that stores the inverse of
+     * lumped diagonal elements as a vector.
+     */
+    std::shared_ptr<DiagonalMatrix<VectorType>> inverse_lumped_diagonal_entries;
   };
 
 
@@ -1861,6 +1902,107 @@ namespace MatrixFreeOperators
 
 
 
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename VectorType,
+            typename VectorizedArrayType>
+  void
+  MassOperator<dim,
+               fe_degree,
+               n_q_points_1d,
+               n_components,
+               VectorType,
+               VectorizedArrayType>::compute_lumped_diagonal()
+  {
+    using Number =
+      typename Base<dim, VectorType, VectorizedArrayType>::value_type;
+    Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
+           ExcNotInitialized());
+    Assert(this->selected_rows == this->selected_columns,
+           ExcMessage("This function is only implemented for square (not "
+                      "rectangular) operators."));
+
+    inverse_lumped_diagonal_entries =
+      std::make_shared<DiagonalMatrix<VectorType>>();
+    lumped_diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
+    VectorType &inverse_lumped_diagonal_vector =
+      inverse_lumped_diagonal_entries->get_vector();
+    VectorType &lumped_diagonal_vector = lumped_diagonal_entries->get_vector();
+    this->initialize_dof_vector(inverse_lumped_diagonal_vector);
+    this->initialize_dof_vector(lumped_diagonal_vector);
+
+    // Re-use the inverse_lumped_diagonal_vector as the vector of 1s
+    inverse_lumped_diagonal_vector = Number(1.);
+    apply_add(lumped_diagonal_vector, inverse_lumped_diagonal_vector);
+    this->set_constrained_entries_to_one(lumped_diagonal_vector);
+
+    const size_type locally_owned_size =
+      inverse_lumped_diagonal_vector.locally_owned_size();
+    // A caller may request a lumped diagonal matrix when it doesn't make sense
+    // (e.g., an element with negative-mean basis functions). Avoid division by
+    // zero so we don't cause a floating point exception but permit negative
+    // entries here.
+    for (size_type i = 0; i < locally_owned_size; ++i)
+      {
+        if (inverse_lumped_diagonal_vector.local_element(i) == Number(0.))
+          inverse_lumped_diagonal_vector.local_element(i) = Number(1.);
+        else
+          inverse_lumped_diagonal_vector.local_element(i) =
+            Number(1.) / lumped_diagonal_vector.local_element(i);
+      }
+
+    inverse_lumped_diagonal_vector.update_ghost_values();
+    lumped_diagonal_vector.update_ghost_values();
+  }
+
+
+
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename VectorType,
+            typename VectorizedArrayType>
+  const std::shared_ptr<DiagonalMatrix<VectorType>> &
+  MassOperator<dim,
+               fe_degree,
+               n_q_points_1d,
+               n_components,
+               VectorType,
+               VectorizedArrayType>::get_matrix_lumped_diagonal_inverse() const
+  {
+    Assert(inverse_lumped_diagonal_entries.get() != nullptr &&
+             inverse_lumped_diagonal_entries->m() > 0,
+           ExcNotInitialized());
+    return inverse_lumped_diagonal_entries;
+  }
+
+
+
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename VectorType,
+            typename VectorizedArrayType>
+  const std::shared_ptr<DiagonalMatrix<VectorType>> &
+  MassOperator<dim,
+               fe_degree,
+               n_q_points_1d,
+               n_components,
+               VectorType,
+               VectorizedArrayType>::get_matrix_lumped_diagonal() const
+  {
+    Assert(lumped_diagonal_entries.get() != nullptr &&
+             lumped_diagonal_entries->m() > 0,
+           ExcNotInitialized());
+    return lumped_diagonal_entries;
+  }
+
+
+
   template <int dim,
             int fe_degree,
             int n_q_points_1d,

diff --git a/tests/matrix_free/mass_operator_02.cc b/tests/matrix_free/mass_operator_02.cc
index 80e2c920388aae4baad4495dbb1aec415ca6a400..dc5100081fb5406b6b42e1923ed6d609ba987463 100644 (file)
--- a/tests/matrix_free/mass_operator_02.cc
+++ b/tests/matrix_free/mass_operator_02.cc
@@ -168,6 +168,7 @@ test()
   }
   sparse_matrix.compress(VectorOperation::add);
 
+  // Check the diagonal:
   for (unsigned int i = 0; i < ref.local_size(); ++i)
     {
       const auto glob_index = owned_set.nth_index_in_set(i);
@@ -179,13 +180,32 @@ test()
   ref.compress(VectorOperation::insert);
 
   out -= ref;
-  const double diff_norm = out.linfty_norm();
 
-  deallog << "Norm of difference: " << diff_norm << std::endl;
-  deallog << "l2_norm: " << inverse_diagonal.l2_norm() << std::endl;
-  deallog << "l1_norm: " << inverse_diagonal.l1_norm() << std::endl;
-  deallog << "linfty_norm: " << inverse_diagonal.linfty_norm() << std::endl
-          << std::endl;
+  deallog << "Norm of difference: " << out.linfty_norm() << std::endl;
+  deallog << "l2_norm: " << ref.l2_norm() << std::endl;
+  deallog << "l1_norm: " << ref.l1_norm() << std::endl;
+  deallog << "linfty_norm: " << ref.linfty_norm() << std::endl << std::endl;
+
+  // Check the lumped diagonal:
+  mf.compute_lumped_diagonal();
+  out = mf.get_matrix_lumped_diagonal_inverse()->get_vector();
+  sparse_matrix.vmult(ref, in);
+  for (unsigned int i = 0; i < ref.local_size(); ++i)
+    {
+      const auto glob_index = owned_set.nth_index_in_set(i);
+      if (constraints.is_constrained(glob_index))
+        ref.local_element(i) = 1.;
+      else
+        ref.local_element(i) = 1. / ref.local_element(i);
+    }
+  ref.compress(VectorOperation::insert);
+
+  out -= ref;
+
+  deallog << "Norm of difference: " << out.linfty_norm() << std::endl;
+  deallog << "l2_norm: " << ref.l2_norm() << std::endl;
+  deallog << "l1_norm: " << ref.l1_norm() << std::endl;
+  deallog << "linfty_norm: " << ref.linfty_norm() << std::endl << std::endl;
 }
 
 

diff --git a/tests/matrix_free/mass_operator_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output b/tests/matrix_free/mass_operator_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output
index 3a595b1c9b95d050b6a0e997f7c9abd535a58202..f6ca581abe6bf06c58ca5385416d49f591892e8e 100644 (file)
--- a/tests/matrix_free/mass_operator_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output
+++ b/tests/matrix_free/mass_operator_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output
@@ -5,21 +5,41 @@ DEAL:2d::l2_norm: 1008.02
 DEAL:2d::l1_norm: 7088.00
 DEAL:2d::linfty_norm: 144.000
 DEAL:2d::
+DEAL:2d::Norm of difference: 2.84217e-14
+DEAL:2d::l2_norm: 504.352
+DEAL:2d::l1_norm: 3536.64
+DEAL:2d::linfty_norm: 92.1600
+DEAL:2d::
 DEAL:2d::Testing FE_Q<2>(2)
 DEAL:2d::Norm of difference: 3.41061e-13
 DEAL:2d::l2_norm: 8100.00
 DEAL:2d::l1_norm: 108964.
 DEAL:2d::linfty_norm: 900.000
 DEAL:2d::
+DEAL:2d::Norm of difference: 2.84217e-13
+DEAL:2d::l2_norm: 5051.63
+DEAL:2d::l1_norm: 68866.9
+DEAL:2d::linfty_norm: 576.000
+DEAL:2d::
 DEAL:3d::Testing FE_Q<3>(1)
 DEAL:3d::Norm of difference: 2.27374e-13
 DEAL:3d::l2_norm: 32003.0
 DEAL:3d::l1_norm: 593090.
 DEAL:3d::linfty_norm: 1728.00
 DEAL:3d::
+DEAL:3d::Norm of difference: 2.27374e-13
+DEAL:3d::l2_norm: 11325.6
+DEAL:3d::l1_norm: 207861.
+DEAL:3d::linfty_norm: 884.736
+DEAL:3d::
 DEAL:3d::Testing FE_Q<3>(2)
 DEAL:3d::Norm of difference: 2.54659e-11
 DEAL:3d::l2_norm: 729000.
 DEAL:3d::l1_norm: 3.59385e+07
 DEAL:3d::linfty_norm: 27000.0
 DEAL:3d::
+DEAL:3d::Norm of difference: 1.27329e-11
+DEAL:3d::l2_norm: 359043.
+DEAL:3d::l1_norm: 1.80487e+07
+DEAL:3d::linfty_norm: 13824.0
+DEAL:3d::

diff --git a/tests/matrix_free/mass_operator_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=4.output b/tests/matrix_free/mass_operator_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=4.output
index 602077b3a4ab3c023d4b389cab6bac08e9e6d2b7..d4f80e1caef5f814e11061ac0ad0c5c7f2b65bdd 100644 (file)
--- a/tests/matrix_free/mass_operator_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=4.output
+++ b/tests/matrix_free/mass_operator_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=4.output
@@ -5,21 +5,41 @@ DEAL:2d::l2_norm: 1008.02
 DEAL:2d::l1_norm: 7088.00
 DEAL:2d::linfty_norm: 144.000
 DEAL:2d::
+DEAL:2d::Norm of difference: 7.10543e-14
+DEAL:2d::l2_norm: 504.352
+DEAL:2d::l1_norm: 3536.64
+DEAL:2d::linfty_norm: 92.1600
+DEAL:2d::
 DEAL:2d::Testing FE_Q<2>(2)
 DEAL:2d::Norm of difference: 7.95808e-13
 DEAL:2d::l2_norm: 8100.00
 DEAL:2d::l1_norm: 108964.
 DEAL:2d::linfty_norm: 900.000
 DEAL:2d::
+DEAL:2d::Norm of difference: 7.95808e-13
+DEAL:2d::l2_norm: 5051.63
+DEAL:2d::l1_norm: 68866.9
+DEAL:2d::linfty_norm: 576.000
+DEAL:2d::
 DEAL:3d::Testing FE_Q<3>(1)
 DEAL:3d::Norm of difference: 6.82121e-13
 DEAL:3d::l2_norm: 32003.0
 DEAL:3d::l1_norm: 593090.
 DEAL:3d::linfty_norm: 1728.00
 DEAL:3d::
+DEAL:3d::Norm of difference: 5.68434e-13
+DEAL:3d::l2_norm: 11325.6
+DEAL:3d::l1_norm: 207861.
+DEAL:3d::linfty_norm: 884.736
+DEAL:3d::
 DEAL:3d::Testing FE_Q<3>(2)
 DEAL:3d::Norm of difference: 4.00178e-11
 DEAL:3d::l2_norm: 729000.
 DEAL:3d::l1_norm: 3.59385e+07
 DEAL:3d::linfty_norm: 27000.0
 DEAL:3d::
+DEAL:3d::Norm of difference: 3.81988e-11
+DEAL:3d::l2_norm: 359043.
+DEAL:3d::l1_norm: 1.80487e+07
+DEAL:3d::linfty_norm: 13824.0
+DEAL:3d::

diff --git a/tests/matrix_free/mass_operator_03.cc b/tests/matrix_free/mass_operator_03.cc
index da50f183fcb064281f43448f663f65ad585f2772..21e766516e04a59aad43ca7a328fa6adf6d5f0b1 100644 (file)
--- a/tests/matrix_free/mass_operator_03.cc
+++ b/tests/matrix_free/mass_operator_03.cc
@@ -93,7 +93,7 @@ test()
     mf;
   mf.initialize(mf_data);
   mf.compute_diagonal();
-  const LinearAlgebra::distributed::Vector<double> &diagonal =
+  const LinearAlgebra::distributed::Vector<double> &inverse_diagonal =
     mf.get_matrix_diagonal_inverse()->get_vector();
 
   LinearAlgebra::distributed::Vector<number> in, out, ref;
@@ -110,7 +110,7 @@ test()
     }
 
   in.update_ghost_values();
-  out = diagonal;
+  out = inverse_diagonal;
 
   // assemble trilinos sparse matrix with
   // (v, u) for reference
@@ -165,6 +165,7 @@ test()
   }
   sparse_matrix.compress(VectorOperation::add);
 
+  // Check the diagonal:
   for (unsigned int i = 0; i < ref.local_size(); ++i)
     {
       const auto glob_index = owned_set.nth_index_in_set(i);
@@ -176,13 +177,32 @@ test()
   ref.compress(VectorOperation::insert);
 
   out -= ref;
-  const double diff_norm = out.linfty_norm();
 
-  deallog << "Norm of difference: " << diff_norm << std::endl;
-  deallog << "l2_norm: " << diagonal.l2_norm() << std::endl;
-  deallog << "l1_norm: " << diagonal.l1_norm() << std::endl;
-  deallog << "linfty_norm: " << diagonal.linfty_norm() << std::endl
-          << std::endl;
+  deallog << "Norm of difference: " << out.linfty_norm() << std::endl;
+  deallog << "l2_norm: " << ref.l2_norm() << std::endl;
+  deallog << "l1_norm: " << ref.l1_norm() << std::endl;
+  deallog << "linfty_norm: " << ref.linfty_norm() << std::endl << std::endl;
+
+  // Check the lumped diagonal:
+  mf.compute_lumped_diagonal();
+  out = mf.get_matrix_lumped_diagonal_inverse()->get_vector();
+  sparse_matrix.vmult(ref, in);
+  for (unsigned int i = 0; i < ref.local_size(); ++i)
+    {
+      const auto glob_index = owned_set.nth_index_in_set(i);
+      if (constraints.is_constrained(glob_index))
+        ref.local_element(i) = 1.;
+      else
+        ref.local_element(i) = 1. / ref.local_element(i);
+    }
+  ref.compress(VectorOperation::insert);
+
+  out -= ref;
+
+  deallog << "Norm of difference: " << out.linfty_norm() << std::endl;
+  deallog << "l2_norm: " << ref.l2_norm() << std::endl;
+  deallog << "l1_norm: " << ref.l1_norm() << std::endl;
+  deallog << "linfty_norm: " << ref.linfty_norm() << std::endl << std::endl;
 }
 
 

diff --git a/tests/matrix_free/mass_operator_03.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output b/tests/matrix_free/mass_operator_03.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output
index 5e1cd4dedd1f98e1c38a95817e5e8a1e873bd5d1..eef201a60191f7f84ad3ac06a60da769f959b428 100644 (file)
--- a/tests/matrix_free/mass_operator_03.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output
+++ b/tests/matrix_free/mass_operator_03.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output
@@ -5,21 +5,41 @@ DEAL:2d::l2_norm: 2160.00
 DEAL:2d::l1_norm: 17424.0
 DEAL:2d::linfty_norm: 576.000
 DEAL:2d::
+DEAL:2d::Norm of difference: 2.84217e-14
+DEAL:2d::l2_norm: 960.000
+DEAL:2d::l1_norm: 7744.00
+DEAL:2d::linfty_norm: 256.000
+DEAL:2d::
 DEAL:2d::Testing FE_Q_Hierarchical<2>(2)
 DEAL:2d::Norm of difference: 2.27374e-13
 DEAL:2d::l2_norm: 3960.00
 DEAL:2d::l1_norm: 63504.0
 DEAL:2d::linfty_norm: 576.000
 DEAL:2d::
+DEAL:2d::Norm of difference: 4.27463e-11
+DEAL:2d::l2_norm: 37440.0
+DEAL:2d::l1_norm: 553536.
+DEAL:2d::linfty_norm: 3600.00
+DEAL:2d::
 DEAL:3d::Testing FE_Q_Hierarchical<3>(1)
 DEAL:3d::Norm of difference: 3.63798e-12
 DEAL:3d::l2_norm: 100388.
 DEAL:3d::l1_norm: 2.29997e+06
 DEAL:3d::linfty_norm: 13824.0
 DEAL:3d::
+DEAL:3d::Norm of difference: 9.09495e-13
+DEAL:3d::l2_norm: 29744.5
+DEAL:3d::l1_norm: 681472.
+DEAL:3d::linfty_norm: 4096.00
+DEAL:3d::
 DEAL:3d::Testing FE_Q_Hierarchical<3>(2)
 DEAL:3d::Norm of difference: 1.63709e-11
 DEAL:3d::l2_norm: 249197.
 DEAL:3d::l1_norm: 1.60030e+07
 DEAL:3d::linfty_norm: 13824.0
 DEAL:3d::
+DEAL:3d::Norm of difference: 1.40572e-08
+DEAL:3d::l2_norm: 7.24442e+06
+DEAL:3d::l1_norm: 4.11831e+08
+DEAL:3d::linfty_norm: 216000.
+DEAL:3d::

diff --git a/tests/matrix_free/mass_operator_03.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=4.output b/tests/matrix_free/mass_operator_03.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=4.output
index c199e6b0c571cfa5de0fcf28facb232be2258d10..2f99b53c2a804d9d65ccc95c5a4011eac8429cd0 100644 (file)
--- a/tests/matrix_free/mass_operator_03.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=4.output
+++ b/tests/matrix_free/mass_operator_03.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=4.output
@@ -5,21 +5,41 @@ DEAL:2d::l2_norm: 2160.00
 DEAL:2d::l1_norm: 17424.0
 DEAL:2d::linfty_norm: 576.000
 DEAL:2d::
+DEAL:2d::Norm of difference: 1.13687e-13
+DEAL:2d::l2_norm: 960.000
+DEAL:2d::l1_norm: 7744.00
+DEAL:2d::linfty_norm: 256.000
+DEAL:2d::
 DEAL:2d::Testing FE_Q_Hierarchical<2>(2)
 DEAL:2d::Norm of difference: 5.68434e-13
 DEAL:2d::l2_norm: 3960.00
 DEAL:2d::l1_norm: 63504.0
 DEAL:2d::linfty_norm: 576.000
 DEAL:2d::
+DEAL:2d::Norm of difference: 4.27463e-11
+DEAL:2d::l2_norm: 37440.0
+DEAL:2d::l1_norm: 553536.
+DEAL:2d::linfty_norm: 3600.00
+DEAL:2d::
 DEAL:3d::Testing FE_Q_Hierarchical<3>(1)
 DEAL:3d::Norm of difference: 5.45697e-12
 DEAL:3d::l2_norm: 100388.
 DEAL:3d::l1_norm: 2.29997e+06
 DEAL:3d::linfty_norm: 13824.0
 DEAL:3d::
+DEAL:3d::Norm of difference: 1.36424e-12
+DEAL:3d::l2_norm: 29744.5
+DEAL:3d::l1_norm: 681472.
+DEAL:3d::linfty_norm: 4096.00
+DEAL:3d::
 DEAL:3d::Testing FE_Q_Hierarchical<3>(2)
 DEAL:3d::Norm of difference: 3.63798e-11
 DEAL:3d::l2_norm: 249197.
 DEAL:3d::l1_norm: 1.60030e+07
 DEAL:3d::linfty_norm: 13824.0
 DEAL:3d::
+DEAL:3d::Norm of difference: 1.40572e-08
+DEAL:3d::l2_norm: 7.24442e+06
+DEAL:3d::l1_norm: 4.11831e+08
+DEAL:3d::linfty_norm: 216000.
+DEAL:3d::