From: David Wells Date: Thu, 29 Jul 2021 21:40:18 +0000 (-0400) Subject: Make MassOperator's diagonals true diagonals, not lumped ones. X-Git-Tag: v9.4.0-rc1~1071^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7af998f7d482c5f6994def49951900c6306cf0b8;p=dealii.git Make MassOperator's diagonals true diagonals, not lumped ones. The previous implementation used row sums (i.e., mass lumping) instead of computing the diagonal. This should be changed for two reasons: 1. It's not really the diagonal (and its inconsistent with LaplaceOperator). 2. Lumped mass matrices are only positive for either low-order elements (where the basis functions are all positive) or very well-behaved elements (like FE_Q). In particular, the values on the diagonal are either zero or negative for TET10 if we use lumping, which isn't going to work. --- diff --git a/doc/news/changes/incompatibilities/20210729DavidWells b/doc/news/changes/incompatibilities/20210729DavidWells new file mode 100644 index 0000000000..e4ddc17dfe --- /dev/null +++ b/doc/news/changes/incompatibilities/20210729DavidWells @@ -0,0 +1,4 @@ +Changed: MatrixFreeOperators::MassOperator now computes its true diagonal (and +not a lumped mass matrix) as its diagonal. +
+(David Wells, 2021/07/29) diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index 8465ea3415..3608a37ffe 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -29,6 +29,7 @@ #include #include +#include #include @@ -40,7 +41,7 @@ namespace MatrixFreeOperators { namespace BlockHelper { - // workaroud for unifying non-block vector and block vector implementations + // workaround for unifying non-block vector and block vector implementations // a non-block vector has one block and the only subblock is the vector // itself template @@ -357,6 +358,11 @@ namespace MatrixFreeOperators * A derived class needs to implement this function and resize and fill * the protected member inverse_diagonal_entries and/or diagonal_entries * accordingly. + * + * @note Since the diagonal is frequently used as a smoother or + * preconditioner, entries corresponding to constrained DoFs are set to 1 + * (instead of the correct value of 0) so that the diagonal matrix is + * invertible. */ virtual void compute_diagonal() = 0; @@ -1793,6 +1799,9 @@ namespace MatrixFreeOperators 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.")); this->inverse_diagonal_entries = std::make_shared>(); @@ -1802,17 +1811,49 @@ namespace MatrixFreeOperators VectorType &diagonal_vector = this->diagonal_entries->get_vector(); this->initialize_dof_vector(inverse_diagonal_vector); this->initialize_dof_vector(diagonal_vector); - inverse_diagonal_vector = Number(1.); - apply_add(diagonal_vector, inverse_diagonal_vector); + // Set up the action of the mass matrix in a way that's compatible with + // MatrixFreeTools::compute_diagonal: + auto diagonal_evaluation = [](auto &integrator) { + integrator.evaluate(EvaluationFlags::values); + for (unsigned int q = 0; q < integrator.n_q_points; ++q) + integrator.submit_value(integrator.get_value(q), q); + integrator.integrate(EvaluationFlags::values); + }; + + std::function::value_type, + VectorizedArrayType> &)> + diagonal_evaluation_f(diagonal_evaluation); + + Assert(this->selected_rows.size() > 0, ExcInternalError()); + for (unsigned int block_n = 0; block_n < this->selected_rows.size(); + ++block_n) + MatrixFreeTools::compute_diagonal(*this->data, + BlockHelper::subblock(diagonal_vector, + block_n), + diagonal_evaluation_f, + this->selected_rows[block_n]); + + // Constrained entries will create zeros on the main diagonal, which we + // don't want this->set_constrained_entries_to_one(diagonal_vector); + inverse_diagonal_vector = diagonal_vector; - const unsigned int locally_owned_size = - inverse_diagonal_vector.locally_owned_size(); - for (unsigned int i = 0; i < locally_owned_size; ++i) - inverse_diagonal_vector.local_element(i) = - Number(1.) / inverse_diagonal_vector.local_element(i); + for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size(); + ++i) + { + Assert(diagonal_vector.local_element(i) > Number(0), + ExcInternalError()); + inverse_diagonal_vector.local_element(i) = + 1. / inverse_diagonal_vector.local_element(i); + } inverse_diagonal_vector.update_ghost_values(); diagonal_vector.update_ghost_values();