From a2dc8737ba26b74a51f89e4ebadcdbf65befaa1b Mon Sep 17 00:00:00 2001 From: Sebastian Kinnewig Date: Fri, 26 Jan 2024 11:18:59 +0100 Subject: [PATCH] Rework import_elements() in read_write_vector for TpetraWrapper::Vector --- .../deal.II/lac/read_write_vector.templates.h | 41 ++++++++++++------- 1 file changed, 26 insertions(+), 15 deletions(-) diff --git a/include/deal.II/lac/read_write_vector.templates.h b/include/deal.II/lac/read_write_vector.templates.h index 652990854b..3991ae6af6 100644 --- a/include/deal.II/lac/read_write_vector.templates.h +++ b/include/deal.II/lac/read_write_vector.templates.h @@ -559,8 +559,7 @@ namespace LinearAlgebra -#ifdef DEAL_II_WITH_TRILINOS -# ifdef DEAL_II_TRILINOS_WITH_TPETRA +#ifdef DEAL_II_TRILINOS_WITH_TPETRA template template std::enable_if_t && @@ -614,26 +613,37 @@ namespace LinearAlgebra Tpetra::Vector target_vector( tpetra_export.getSourceMap()); - target_vector.doImport(vector, tpetra_export, Tpetra::REPLACE); - const auto *new_values = target_vector.getData().get(); - const auto size = target_vector.getLocalLength(); + // Communicate the vector to the correct map. + // Remark: We use here doImport on an Export object since we have to use + // the communication plan stored in the tpetra_comm_patern backward. + target_vector.doImport(vector, tpetra_export, Tpetra::INSERT); + +# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0) + auto vector_2d = target_vector.template getLocalView( + Tpetra::Access::ReadOnly); +# else + target_vector.template sync(); + auto vector_2d = target_vector.template getLocalView(); +# endif + auto new_values = Kokkos::subview(vector_2d, Kokkos::ALL(), 0); + auto size = target_vector.getLocalLength(); using size_type = std::decay_t; - Assert(size == 0 || values != nullptr, ExcInternalError("Export failed.")); + Assert(size == 0 || values != nullptr, ExcInternalError("Import failed.")); AssertDimension(size, stored_elements.n_elements()); switch (operation) { case VectorOperation::insert: for (size_type i = 0; i < size; ++i) - values[i] = new_values[i]; + values[i] = new_values(i); break; case VectorOperation::add: for (size_type i = 0; i < size; ++i) - values[i] += new_values[i]; + values[i] += new_values(i); break; case VectorOperation::min: @@ -644,15 +654,15 @@ namespace LinearAlgebra for (size_type i = 0; i < size; ++i) { Assert( - std::imag(new_values[i]) == 0., + std::imag(new_values(i)) == 0., ExcMessage( "VectorOperation::min is not defined if there is an imaginary part!)")); Assert( std::imag(values[i]) == 0., ExcMessage( "VectorOperation::min is not defined if there is an imaginary part!)")); - if (std::real(new_values[i]) - std::real(values[i]) < 0.0) - values[i] = new_values[i]; + if (std::real(new_values(i)) - std::real(values[i]) < 0.0) + values[i] = new_values(i); } break; @@ -660,15 +670,15 @@ namespace LinearAlgebra for (size_type i = 0; i < size; ++i) { Assert( - std::imag(new_values[i]) == 0., + std::imag(new_values(i)) == 0., ExcMessage( "VectorOperation::max is not defined if there is an imaginary part!)")); Assert( std::imag(values[i]) == 0., ExcMessage( "VectorOperation::max is not defined if there is an imaginary part!)")); - if (std::real(new_values[i]) - std::real(values[i]) > 0.0) - values[i] = new_values[i]; + if (std::real(new_values(i)) - std::real(values[i]) > 0.0) + values[i] = new_values(i); } break; @@ -676,10 +686,11 @@ namespace LinearAlgebra AssertThrow(false, ExcNotImplemented()); } } -# endif +#endif +#ifdef DEAL_II_WITH_TRILINOS template void ReadWriteVector::import_elements( -- 2.39.5