From 1e3ce0cce493fe6c29445f838c215b4492f1fd5e Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 31 Jul 2021 12:14:17 -0400 Subject: [PATCH] Implement explicit mass lumping in MassOperator. --- include/deal.II/matrix_free/operators.h | 146 +++++++++++++++++- tests/matrix_free/mass_operator_02.cc | 32 +++- ...h_mpi=true.with_p4est=true.mpirun=1.output | 20 +++ ...h_mpi=true.with_p4est=true.mpirun=4.output | 20 +++ tests/matrix_free/mass_operator_03.cc | 36 ++++- ...h_mpi=true.with_p4est=true.mpirun=1.output | 20 +++ ...h_mpi=true.with_p4est=true.mpirun=4.output | 20 +++ 7 files changed, 278 insertions(+), 16 deletions(-) diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index 3608a37ffe..017af8daa4 100644 --- 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> & + get_matrix_lumped_diagonal() const; + + /** + * Get read access to the inverse lumped diagonal of this operator. + */ + const std::shared_ptr> & + 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 & cell_range) const; + + /** + * A shared pointer to a diagonal matrix that stores the + * lumped diagonal elements as a vector. + */ + std::shared_ptr> lumped_diagonal_entries; + + /** + * A shared pointer to a diagonal matrix that stores the inverse of + * lumped diagonal elements as a vector. + */ + std::shared_ptr> inverse_lumped_diagonal_entries; }; @@ -1861,6 +1902,107 @@ namespace MatrixFreeOperators + template + void + MassOperator::compute_lumped_diagonal() + { + using Number = + typename Base::value_type; + Assert((Base::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>(); + lumped_diagonal_entries = std::make_shared>(); + 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 + const std::shared_ptr> & + MassOperator::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 + const std::shared_ptr> & + MassOperator::get_matrix_lumped_diagonal() const + { + Assert(lumped_diagonal_entries.get() != nullptr && + lumped_diagonal_entries->m() > 0, + ExcNotInitialized()); + return lumped_diagonal_entries; + } + + + template 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 3a595b1c9b..f6ca581abe 100644 --- 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 602077b3a4..d4f80e1cae 100644 --- 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 da50f183fc..21e766516e 100644 --- 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 &diagonal = + const LinearAlgebra::distributed::Vector &inverse_diagonal = mf.get_matrix_diagonal_inverse()->get_vector(); LinearAlgebra::distributed::Vector 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 5e1cd4dedd..eef201a601 100644 --- 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 c199e6b0c5..2f99b53c2a 100644 --- 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:: -- 2.39.5