From: Peter Munch Date: Fri, 25 Feb 2022 18:18:04 +0000 (+0100) Subject: Improve asserts in parallel::distributed::SolutionTransfer X-Git-Tag: v9.4.0-rc1~238^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0ba4fca4e6acf74dcb72478f4533b0e03e46909a;p=dealii.git Improve asserts in parallel::distributed::SolutionTransfer --- diff --git a/doc/news/changes/minor/20220225Munch b/doc/news/changes/minor/20220225Munch new file mode 100644 index 0000000000..1e461f06ac --- /dev/null +++ b/doc/news/changes/minor/20220225Munch @@ -0,0 +1,9 @@ +Improved: The asserts in parallel::distributed::SolutionTransfer have +been improved. Before this change, an assert has been only called +in the parallel case for distributed vectors for cells neighboring +interprocess boundaries. Now, we check for all DoFs that the +contributions from all cells have the same values. This condition +might be not given if the solution has high gradients at hanging +nodes. +
+(Peter Munch, Magdalena Schreter, 2022/02/25) diff --git a/include/deal.II/base/partitioner.templates.h b/include/deal.II/base/partitioner.templates.h index 221b3e5542..34fe115dac 100644 --- a/include/deal.II/base/partitioner.templates.h +++ b/include/deal.II/base/partitioner.templates.h @@ -620,7 +620,10 @@ namespace Utilities // p::d::SolutionTransfer. The rationale is that during // interpolation on two elements sharing the face, values on // this face obtained from each side might be different due to - // additions being done in different order. + // additions being done in different order. If the local + // value is zero, it indicates that the local process has not + // set the value during the cell loop and its value can be + // safely overridden. Assert(*read_position == Number() || internal::get_abs(locally_owned_array[j] - *read_position) <= diff --git a/include/deal.II/dofs/dof_accessor.h b/include/deal.II/dofs/dof_accessor.h index 39ee53ef31..2f561c44ed 100644 --- a/include/deal.II/dofs/dof_accessor.h +++ b/include/deal.II/dofs/dof_accessor.h @@ -1743,6 +1743,15 @@ public: * of the respective finite element class for a description of what the * prolongation matrices represent in this case. * + * @note Cells set the values of DoFs independently and might overwrite + * previously set values in the global vector, for example when calling the + * same function earlier from a different cell. By setting @p perform_check, + * you can enable a check that the previous value and the one to be set here + * are at least roughly the same. In practice, they might be slightly + * different because they are computed in a way that theoretically ensures + * that they are the same, but in practice they are only equal up to + * round-off. + * * @note Unlike the get_dof_values() function, this function is only * available on cells, rather than on lines, quads, and hexes, since * interpolation is presently only provided for cells by the finite element @@ -1754,7 +1763,8 @@ public: const Vector &local_values, OutputVector & values, const unsigned int fe_index = - DoFHandler::invalid_fe_index) const; + DoFHandler::invalid_fe_index, + const bool perform_check = false) const; /** * Distribute a local (cell based) vector to a global one by mapping the diff --git a/source/distributed/solution_transfer.cc b/source/distributed/solution_transfer.cc index 99dcf55b42..889658aad0 100644 --- a/source/distributed/solution_transfer.cc +++ b/source/distributed/solution_transfer.cc @@ -424,7 +424,8 @@ namespace parallel for (; it_input != dof_values.cend(); ++it_input, ++it_output) cell->set_dof_values_by_interpolation(*it_input, *(*it_output), - fe_index); + fe_index, + true); } } // namespace distributed } // namespace parallel diff --git a/source/dofs/dof_accessor_set.cc b/source/dofs/dof_accessor_set.cc index ca5dc3cf96..bd2f50e5c6 100644 --- a/source/dofs/dof_accessor_set.cc +++ b/source/dofs/dof_accessor_set.cc @@ -41,6 +41,229 @@ DEAL_II_NAMESPACE_OPEN +template +DeclException2(ExcNonMatchingElementsSetDofValuesByInterpolation, + Number, + Number, + << "Called set_dof_values_by_interpolation(), but" + << " the element to be set, value " << std::setprecision(16) + << arg1 << ", does not match with the non-zero value " + << std::setprecision(16) << arg2 << " already set before."); + +namespace internal +{ + /** + * In the set_dof_values(), we need to invoke abs() also on unsigned data + * types, which is ill-formed on newer C++ standards. To avoid this, we use + * std::abs on default types, but simply return the number on unsigned types. + */ + template + std::enable_if_t::value, + typename numbers::NumberTraits::real_type> + get_abs(const Number a) + { + return std::abs(a); + } + + template + std::enable_if_t::value, Number> + get_abs(const Number a) + { + return a; + } + + /** + * Check if a vector is a deal.II vector. + */ + template + constexpr bool is_dealii_vector = + std::is_same>::value || + std::is_same>::value || + std::is_same< + VectorType, + dealii::LinearAlgebra::Vector>::value || + std::is_same>::value || + std::is_same>::value; + + /** + * Helper functions that call set_ghost_state() if the vector supports this + * operation. + */ + template + using set_ghost_state_t = + decltype(std::declval().set_ghost_state(std::declval())); + + template + constexpr bool has_set_ghost_state = + is_supported_operation; + + template , + VectorType>::type * = nullptr> + void + set_ghost_state(VectorType &vector, const bool ghosted) + { + vector.set_ghost_state(ghosted); + } + + template , + VectorType>::type * = nullptr> + void + set_ghost_state(VectorType &, const bool) + { + // serial vector: nothing to do + } + + /** + * Helper function that sets the values on a cell, but also checks if the + * new values are similar to the old values. + */ + template + void + set_dof_values(const DoFCellAccessor &cell, + const Vector & local_values, + OutputVector & values, + const bool perform_check) + { + (void)perform_check; + +#ifdef DEBUG + if (perform_check && is_dealii_vector) + { + const bool old_ghost_state = values.has_ghost_elements(); + set_ghost_state(values, true); + + Vector local_values_old(cell.get_fe().n_dofs_per_cell()); + cell.get_dof_values(values, local_values_old); + + for (unsigned int i = 0; i < cell.get_fe().n_dofs_per_cell(); ++i) + { + // a check consistent with the one in + // Utilities::MPI::Partitioner::import_from_ghosted_array_finish() + Assert(local_values_old[i] == number() || + get_abs(local_values_old[i] - local_values[i]) <= + get_abs(local_values_old[i] + local_values[i]) * + 100000. * + std::numeric_limits::real_type>::epsilon(), + ExcNonMatchingElementsSetDofValuesByInterpolation( + local_values[i], local_values_old[i])); + } + + set_ghost_state(values, old_ghost_state); + } +#endif + + cell.set_dof_values(local_values, values); + } + + + template + void + process_by_interpolation( + const DoFCellAccessor & cell, + const Vector & local_values, + OutputVector & values, + const unsigned int fe_index_, + const std::function &cell, + const Vector &local_values, + OutputVector & values)> &processor) + { + const unsigned int fe_index = + (cell.get_dof_handler().has_hp_capabilities() == false && + fe_index_ == DoFHandler::invalid_fe_index) ? + DoFHandler::default_fe_index : + fe_index_; + + if (cell.is_active() && !cell.is_artificial()) + { + if ((cell.get_dof_handler().has_hp_capabilities() == false) || + // for hp-DoFHandlers, we need to require that on + // active cells, you either don't specify an fe_index, + // or that you specify the correct one + (fe_index == cell.active_fe_index()) || + (fe_index == DoFHandler::invalid_fe_index)) + // simply set the values on this cell + processor(cell, local_values, values); + else + { + Assert(local_values.size() == + cell.get_dof_handler().get_fe(fe_index).n_dofs_per_cell(), + ExcMessage("Incorrect size of local_values vector.")); + + FullMatrix interpolation( + cell.get_fe().n_dofs_per_cell(), + cell.get_dof_handler().get_fe(fe_index).n_dofs_per_cell()); + + cell.get_fe().get_interpolation_matrix( + cell.get_dof_handler().get_fe(fe_index), interpolation); + + // do the interpolation to the target space. for historical + // reasons, matrices are set to size 0x0 internally even if + // we reinit as 4x0, so we have to treat this case specially + Vector tmp(cell.get_fe().n_dofs_per_cell()); + if ((tmp.size() > 0) && (local_values.size() > 0)) + interpolation.vmult(tmp, local_values); + + // now set the dof values in the global vector + processor(cell, tmp, values); + } + } + else + // otherwise distribute them to the children + { + Assert((cell.get_dof_handler().has_hp_capabilities() == false) || + (fe_index != DoFHandler::invalid_fe_index), + ExcMessage( + "You cannot call this function on non-active cells " + "of DoFHandler objects unless you provide an explicit " + "finite element index because they do not have naturally " + "associated finite element spaces associated: degrees " + "of freedom are only distributed on active cells for which " + "the active FE index has been set.")); + + const FiniteElement &fe = + cell.get_dof_handler().get_fe(fe_index); + const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); + + Assert(local_values.size() == dofs_per_cell, + (typename DoFCellAccessor::BaseClass:: + ExcVectorDoesNotMatch())); + Assert(values.size() == cell.get_dof_handler().n_dofs(), + (typename DoFCellAccessor::BaseClass:: + ExcVectorDoesNotMatch())); + + Vector tmp(dofs_per_cell); + + for (unsigned int child = 0; child < cell.n_children(); ++child) + { + if (tmp.size() > 0) + fe.get_prolongation_matrix(child, cell.refinement_case()) + .vmult(tmp, local_values); + process_by_interpolation( + *cell.child(child), tmp, values, fe_index, processor); + } + } + } + +} // namespace internal + + template template @@ -48,84 +271,19 @@ void DoFCellAccessor::set_dof_values_by_interpolation( const Vector &local_values, OutputVector & values, - const unsigned int fe_index_) const + const unsigned int fe_index_, + const bool perform_check) const { - const unsigned int fe_index = - (this->dof_handler->hp_capability_enabled == false && - fe_index_ == DoFHandler::invalid_fe_index) ? - DoFHandler::default_fe_index : - fe_index_; - - if (this->is_active() && !this->is_artificial()) - { - if ((this->dof_handler->hp_capability_enabled == false) || - // for hp-DoFHandlers, we need to require that on - // active cells, you either don't specify an fe_index, - // or that you specify the correct one - (fe_index == this->active_fe_index()) || - (fe_index == DoFHandler::invalid_fe_index)) - // simply set the values on this cell - this->set_dof_values(local_values, values); - else - { - Assert(local_values.size() == - this->dof_handler->get_fe(fe_index).n_dofs_per_cell(), - ExcMessage("Incorrect size of local_values vector.")); - - FullMatrix interpolation( - this->get_fe().n_dofs_per_cell(), - this->dof_handler->get_fe(fe_index).n_dofs_per_cell()); - - this->get_fe().get_interpolation_matrix( - this->dof_handler->get_fe(fe_index), interpolation); - - // do the interpolation to the target space. for historical - // reasons, matrices are set to size 0x0 internally even - // we reinit as 4x0, so we have to treat this case specially - Vector tmp(this->get_fe().n_dofs_per_cell()); - if ((tmp.size() > 0) && (local_values.size() > 0)) - interpolation.vmult(tmp, local_values); - - // now set the dof values in the global vector - this->set_dof_values(tmp, values); - } - } - else - // otherwise distribute them to the children - { - Assert((this->dof_handler->hp_capability_enabled == false) || - (fe_index != DoFHandler::invalid_fe_index), - ExcMessage( - "You cannot call this function on non-active cells " - "of DoFHandler objects unless you provide an explicit " - "finite element index because they do not have naturally " - "associated finite element spaces associated: degrees " - "of freedom are only distributed on active cells for which " - "the active FE index has been set.")); - - const FiniteElement &fe = - this->get_dof_handler().get_fe(fe_index); - const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); - - Assert(this->dof_handler != nullptr, - typename BaseClass::ExcInvalidObject()); - Assert(local_values.size() == dofs_per_cell, - typename BaseClass::ExcVectorDoesNotMatch()); - Assert(values.size() == this->dof_handler->n_dofs(), - typename BaseClass::ExcVectorDoesNotMatch()); - - Vector tmp(dofs_per_cell); - - for (unsigned int child = 0; child < this->n_children(); ++child) - { - if (tmp.size() > 0) - fe.get_prolongation_matrix(child, this->refinement_case()) - .vmult(tmp, local_values); - this->child(child)->set_dof_values_by_interpolation(tmp, - values, - fe_index); - } - } + internal::process_by_interpolation( + *this, + local_values, + values, + fe_index_, + [perform_check](const DoFCellAccessor &cell, + const Vector & local_values, + OutputVector & values) { + internal::set_dof_values(cell, local_values, values, perform_check); + }); } diff --git a/source/dofs/dof_accessor_set.inst.in b/source/dofs/dof_accessor_set.inst.in index cc02b5b864..6bdc91315f 100644 --- a/source/dofs/dof_accessor_set.inst.in +++ b/source/dofs/dof_accessor_set.inst.in @@ -20,7 +20,8 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; lda : BOOL) template void DoFCellAccessor:: set_dof_values_by_interpolation(const Vector &, VEC &, - const unsigned int fe_index) const; + const unsigned int fe_index, + const bool) const; #if deal_II_dimension != 3 @@ -28,15 +29,18 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; lda : BOOL) DoFCellAccessor:: set_dof_values_by_interpolation(const Vector &, VEC &, - const unsigned int fe_index) const; + const unsigned int fe_index, + const bool) const; #endif #if deal_II_dimension == 3 template void DoFCellAccessor<1, 3, lda>::set_dof_values_by_interpolation( - const Vector &, VEC &, const unsigned int fe_index) - const; + const Vector &, + VEC &, + const unsigned int fe_index, + const bool) const; #endif } diff --git a/source/numerics/solution_transfer.cc b/source/numerics/solution_transfer.cc index a7175306cf..685c486763 100644 --- a/source/numerics/solution_transfer.cc +++ b/source/numerics/solution_transfer.cc @@ -192,7 +192,8 @@ SolutionTransfer::refine_interpolate( in, (*pointerstruct->second.indices_ptr)[i]); cell->set_dof_values_by_interpolation(local_values, out, - this_fe_index); + this_fe_index, + true); } } } @@ -520,7 +521,8 @@ SolutionTransfer::interpolate( cell->set_dof_values_by_interpolation(tmp, all_out[j], - old_fe_index); + old_fe_index, + true); } } else if (valuesptr) diff --git a/tests/codim_one/solution_transfer_01.cc b/tests/codim_one/solution_transfer_01.cc index 9dc5e2f8ae..01e95e7471 100644 --- a/tests/codim_one/solution_transfer_01.cc +++ b/tests/codim_one/solution_transfer_01.cc @@ -88,7 +88,6 @@ main() // get the interpolated solution // back Vector tmp(dh.n_dofs()); - tmp = 2; soltrans.interpolate(solution, tmp); deallog << "New values:" << std::endl;