diff --git a/.clang-tidy b/.clang-tidy index 447cf868cead..3e520f69fea2 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -13,6 +13,10 @@ Checks: '*, # complains about Catch macros having infinite loops, -bugprone-infinite-loop, -bugprone-macro-parentheses, +# useless as it triggers on calls to value() which already do a check, +# and is too aggressive as it requires the check to be immediately before the +# access. +-bugprone-unchecked-optional-access, # Official justification is not public, unofficial justification has nothing, # specific to iterators and is basically "never create mutable temporaries",, # which is against standard practice everywhere else., diff --git a/src/Domain/Structure/OrientationMapHelpers.cpp b/src/Domain/Structure/OrientationMapHelpers.cpp index 73755d1d0570..c380c4a2b43f 100644 --- a/src/Domain/Structure/OrientationMapHelpers.cpp +++ b/src/Domain/Structure/OrientationMapHelpers.cpp @@ -315,6 +315,30 @@ void orient_each_component( } } // namespace +template +Mesh orient_mesh_on_slice( + const Mesh& mesh_on_slice, const size_t sliced_dim, + const OrientationMap& orientation_of_neighbor) { + if constexpr (VolumeDim < 3) { + return mesh_on_slice; + } else { + if (orientation_of_neighbor.is_aligned()) { + return mesh_on_slice; + } + const size_t first_dim_of_slice = sliced_dim == 0 ? 1 : 0; + const size_t second_dim_of_slice = sliced_dim == 2 ? 1 : 2; + if (orientation_of_neighbor(first_dim_of_slice) > + orientation_of_neighbor(second_dim_of_slice)) { + return Mesh<2>{ + {{mesh_on_slice.extents(1), mesh_on_slice.extents(0)}}, + {{mesh_on_slice.basis(1), mesh_on_slice.basis(0)}}, + {{mesh_on_slice.quadrature(1), mesh_on_slice.quadrature(0)}}}; + } else { + return mesh_on_slice; + } + } +} + template void orient_variables( const gsl::not_null result, const DataVector& variables, @@ -428,6 +452,9 @@ DataVector orient_variables_on_slice( #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) #define INSTANTIATION(r, data) \ + template Mesh orient_mesh_on_slice( \ + const Mesh& mesh_on_slice, const size_t sliced_dim, \ + const OrientationMap& orientation_of_neighbor); \ template void orient_variables( \ const gsl::not_null result, const DataVector& variables, \ const Index& extents, \ diff --git a/src/Domain/Structure/OrientationMapHelpers.hpp b/src/Domain/Structure/OrientationMapHelpers.hpp index 1cbf702891c6..ea016881b0c4 100644 --- a/src/Domain/Structure/OrientationMapHelpers.hpp +++ b/src/Domain/Structure/OrientationMapHelpers.hpp @@ -13,6 +13,8 @@ /// \cond template +class Mesh; +template class OrientationMap; template class Index; @@ -20,6 +22,14 @@ template class Variables; /// \endcond +/// \ingroup ComputationalDomainGroup +/// \brief Orient a sliced Mesh to the logical frame of a neighbor element with +/// the given orientation. +template +Mesh orient_mesh_on_slice( + const Mesh& mesh_on_slice, size_t sliced_dim, + const OrientationMap& orientation_of_neighbor); + /// @{ /// \ingroup ComputationalDomainGroup /// \brief Orient variables to the data-storage order of a neighbor element with diff --git a/src/Evolution/DgSubcell/Actions/ReconstructionCommunication.hpp b/src/Evolution/DgSubcell/Actions/ReconstructionCommunication.hpp index 80c55475096c..df63fbe77334 100644 --- a/src/Evolution/DgSubcell/Actions/ReconstructionCommunication.hpp +++ b/src/Evolution/DgSubcell/Actions/ReconstructionCommunication.hpp @@ -337,6 +337,7 @@ struct ReceiveDataForReconstruction { inbox.erase(current_time_step_id); const Mesh& subcell_mesh = db::get>(box); + const auto& mortar_meshes = get>(box); db::mutate, Tags::DataForRdmpTci, evolution::dg::Tags::MortarData, @@ -347,7 +348,7 @@ struct ReceiveDataForReconstruction { ghost_zone_size = db::get(box) .ghost_zone_size(), - &received_data, &subcell_mesh]( + &received_data, &subcell_mesh, &mortar_meshes]( const gsl::not_null*> ghost_data_ptr, const gsl::not_null rdmp_tci_data_ptr, @@ -381,6 +382,8 @@ struct ReceiveDataForReconstruction { .has_value()) { mortar_data->at(mortar_id).neighbor().face_mesh = received_mortar_data.second.interface_mesh; + mortar_data->at(mortar_id).neighbor().mortar_mesh = + mortar_meshes.at(mortar_id); mortar_data->at(mortar_id).neighbor().mortar_data = std::move( *received_mortar_data.second.boundary_correction_data); } diff --git a/src/Evolution/DiscontinuousGalerkin/Actions/ApplyBoundaryCorrections.hpp b/src/Evolution/DiscontinuousGalerkin/Actions/ApplyBoundaryCorrections.hpp index a73e89f5da0d..e361375b15ae 100644 --- a/src/Evolution/DiscontinuousGalerkin/Actions/ApplyBoundaryCorrections.hpp +++ b/src/Evolution/DiscontinuousGalerkin/Actions/ApplyBoundaryCorrections.hpp @@ -306,10 +306,12 @@ bool receive_boundary_data_global_time_stepping( make_not_null(&db::as_access(*box)), received_temporal_id_and_data); } + const auto& mortar_meshes = + get>(*box); db::mutate, evolution::dg::Tags::MortarNextTemporalId, domain::Tags::NeighborMesh>( - [&received_temporal_id_and_data]( + [&received_temporal_id_and_data, &mortar_meshes]( const gsl::not_null>*> mortar_data, @@ -338,6 +340,8 @@ bool receive_boundary_data_global_time_stepping( .has_value()) { mortar_data->at(mortar_id).neighbor().face_mesh = received_mortar_data.second.interface_mesh; + mortar_data->at(mortar_id).neighbor().mortar_mesh = + mortar_meshes.at(mortar_id); mortar_data->at(mortar_id).neighbor().mortar_data = std::move( received_mortar_data.second.boundary_correction_data.value()); } @@ -411,13 +415,14 @@ bool receive_boundary_data_local_time_stepping( } ASSERT(inbox_ptr != nullptr, "The inbox pointer should not be null."); InboxMap& inbox = *inbox_ptr; + const auto& mortar_meshes = get>(*box); const bool have_all_intermediate_messages = db::mutate< evolution::dg::Tags::MortarDataHistory, evolution::dg::Tags::MortarNextTemporalId, domain::Tags::NeighborMesh>( - [&inbox, &needed_time]( + [&inbox, &needed_time, &mortar_meshes]( const gsl::not_null, @@ -464,6 +469,7 @@ bool receive_boundary_data_local_time_stepping( received_mortar_data->second.volume_mesh_ghost_cell_data); neighbor_mortar_data.face_mesh = received_mortar_data->second.interface_mesh; + neighbor_mortar_data.mortar_mesh = mortar_meshes.at(mortar_id); neighbor_mortar_data.mortar_data = std::move( received_mortar_data->second.boundary_correction_data.value()); boundary_data_history->at(mortar_id).remote().insert( diff --git a/src/Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivative.hpp b/src/Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivative.hpp index 08851d63da9f..7525d4cde495 100644 --- a/src/Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivative.hpp +++ b/src/Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivative.hpp @@ -844,6 +844,7 @@ void ComputeTimeDerivativeat(mortar_id).local(); local_mortar_data.face_normal_magnitude = face_normal_magnitude; if (using_gauss_points) { + local_mortar_data.volume_mesh = volume_mesh; local_mortar_data.volume_det_inv_jacobian = volume_det_inv_jacobian; local_mortar_data.face_det_jacobian = face_det_jacobian; diff --git a/src/Evolution/DiscontinuousGalerkin/Actions/InternalMortarDataImpl.hpp b/src/Evolution/DiscontinuousGalerkin/Actions/InternalMortarDataImpl.hpp index 3d1e4f54ab09..83bbb4a73fa7 100644 --- a/src/Evolution/DiscontinuousGalerkin/Actions/InternalMortarDataImpl.hpp +++ b/src/Evolution/DiscontinuousGalerkin/Actions/InternalMortarDataImpl.hpp @@ -253,6 +253,7 @@ void internal_mortar_data_impl( // Can use the local_mortar_data auto& local_mortar = mortar_data_ptr->at(mortar_id).local(); local_mortar.face_mesh = face_mesh; + local_mortar.mortar_mesh = mortar_mesh; // If this is the first time, initialize the data. If we don't do this, // then the DataVector will be non-owning which we don't want if (UNLIKELY(not local_mortar.mortar_data.has_value())) { @@ -294,6 +295,7 @@ void internal_mortar_data_impl( if (Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)) { auto& local_mortar = mortar_data_ptr->at(mortar_id).local(); local_mortar.face_mesh = face_mesh; + local_mortar.mortar_mesh = mortar_mesh; // If this is the first time, initialize the data. If we don't do this, // then the DataVector will be non-owning which we don't want if (UNLIKELY(not local_mortar.mortar_data.has_value())) { diff --git a/src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.hpp b/src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.hpp index 6ed3b123f741..ccd0f7e2ee2e 100644 --- a/src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.hpp +++ b/src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.hpp @@ -72,6 +72,98 @@ mortars_apply_impl(const std::vector>& initial_extents, Spectral::Quadrature quadrature, const Element& element, const TimeStepId& next_temporal_id, const Mesh& volume_mesh); + +template +void p_project( + const gsl::not_null< + ::dg::MortarMap>*> + /* mortar_data */, + const gsl::not_null<::dg::MortarMap>*> mortar_mesh, + const gsl::not_null< + ::dg::MortarMap>*> + /* mortar_size */, + const gsl::not_null<::dg::MortarMap*> + /* mortar_next_temporal_id */, + const gsl::not_null< + DirectionMap>>>>*> + normal_covector_and_magnitude, + const gsl::not_null mortar_data_history, + const Mesh& new_mesh, const Element& new_element, + const std::unordered_map, amr::Info>& neighbor_info, + const std::pair, Element>& old_mesh_and_element) { + const auto& [old_mesh, old_element] = old_mesh_and_element; + ASSERT(old_element.id() == new_element.id(), + "p-refinement should not have changed the element id"); + + const bool mesh_changed = old_mesh != new_mesh; + + for (const auto& [direction, neighbors] : new_element.neighbors()) { + const auto sliced_away_dimension = direction.dimension(); + const auto old_face_mesh = old_mesh.slice_away(sliced_away_dimension); + const auto new_face_mesh = new_mesh.slice_away(sliced_away_dimension); + const bool face_mesh_changed = old_face_mesh != new_face_mesh; + if (face_mesh_changed) { + (*normal_covector_and_magnitude)[direction] = std::nullopt; + } + for (const auto& neighbor : neighbors) { + const DirectionalId mortar_id{direction, neighbor}; + if (mortar_mesh->contains(mortar_id)) { + const auto new_neighbor_mesh = neighbors.orientation().inverse_map()( + neighbor_info.at(neighbor).new_mesh); + const auto& old_mortar_mesh = mortar_mesh->at(mortar_id); + auto new_mortar_mesh = ::dg::mortar_mesh( + new_mesh.slice_away(direction.dimension()), + new_neighbor_mesh.slice_away(direction.dimension())); + const bool mortar_mesh_changed = old_mortar_mesh != new_mortar_mesh; + + if (mortar_mesh_changed or mesh_changed) { + // mortar_data does not need projecting as it has already been used + // and will be resized automatically + // mortar_size does not change as the mortar has not changed + // next_temporal_id does not change as the mortar has not changed + if (not mortar_data_history->empty()) { + auto& boundary_history = mortar_data_history->at(mortar_id); + auto local_history = boundary_history.local(); + auto remote_history = boundary_history.remote(); + const auto project_local_boundary_data = + [&new_mortar_mesh, &new_face_mesh, &new_mesh]( + const TimeStepId& /* id */, + const gsl::not_null<::evolution::dg::MortarData*> + mortar_data) { + p_project(mortar_data, new_mortar_mesh, new_face_mesh, + new_mesh); + }; + local_history.for_each(project_local_boundary_data); + const auto project_remote_boundary_data = + [&new_mortar_mesh]( + const TimeStepId& /* id */, + const gsl::not_null<::evolution::dg::MortarData*> + mortar_data) { + p_project_only_mortar_data(mortar_data, new_mortar_mesh); + }; + remote_history.for_each(project_remote_boundary_data); + boundary_history.clear_coupling_cache(); + } + mortar_mesh->at(mortar_id) = std::move(new_mortar_mesh); + } + } else { + ERROR("h-refinement not implemented yet"); + } + } + } + + for (const auto& direction : new_element.external_boundaries()) { + const auto sliced_away_dimension = direction.dimension(); + const auto old_face_mesh = old_mesh.slice_away(sliced_away_dimension); + const auto new_face_mesh = new_mesh.slice_away(sliced_away_dimension); + const bool face_mesh_changed = old_face_mesh != new_face_mesh; + if (face_mesh_changed) { + (*normal_covector_and_magnitude)[direction] = std::nullopt; + } + } +} } // namespace detail /*! @@ -192,48 +284,21 @@ struct ProjectMortars : tt::ConformsTo { const gsl::not_null< ::dg::MortarMap>*> mortar_size, - const gsl::not_null< - ::dg::MortarMap*> /*mortar_next_temporal_id*/, + const gsl::not_null<::dg::MortarMap*> + mortar_next_temporal_id, const gsl::not_null< DirectionMap>>>>*> normal_covector_and_magnitude, - const gsl::not_null - /*mortar_data_history*/, + const gsl::not_null mortar_data_history, const Mesh& new_mesh, const Element& new_element, const std::unordered_map, amr::Info>& neighbor_info, - const std::pair, Element>& /*old_mesh_and_element*/) { - if (Metavariables::local_time_stepping) { - ERROR("AMR with local time-stepping is not yet supported"); - } - - mortar_data->clear(); - mortar_mesh->clear(); - mortar_size->clear(); - // mortar_next_temporal_id is not changed, but this will break when - // h-refinement is enabled and the neighbors are no longer the same - for (const auto& [direction, neighbors] : new_element.neighbors()) { - (*normal_covector_and_magnitude)[direction] = std::nullopt; - for (const auto& neighbor : neighbors) { - const DirectionalId mortar_id{direction, neighbor}; - mortar_data->emplace(mortar_id, MortarDataHolder{}); - const auto new_neighbor_mesh = neighbors.orientation().inverse_map()( - neighbor_info.at(neighbor).new_mesh); - mortar_mesh->emplace( - mortar_id, - ::dg::mortar_mesh( - new_mesh.slice_away(direction.dimension()), - new_neighbor_mesh.slice_away(direction.dimension()))); - mortar_size->emplace( - mortar_id, - ::dg::mortar_size(new_element.id(), neighbor, direction.dimension(), - neighbors.orientation())); - } - } - for (const auto& direction : new_element.external_boundaries()) { - (*normal_covector_and_magnitude)[direction] = std::nullopt; - } + const std::pair, Element>& old_mesh_and_element) { + detail::p_project(mortar_data, mortar_mesh, mortar_size, + mortar_next_temporal_id, normal_covector_and_magnitude, + mortar_data_history, new_mesh, new_element, neighbor_info, + old_mesh_and_element); } template diff --git a/src/Evolution/DiscontinuousGalerkin/MortarData.cpp b/src/Evolution/DiscontinuousGalerkin/MortarData.cpp index 6d7640b3e859..5ca9d8c302fe 100644 --- a/src/Evolution/DiscontinuousGalerkin/MortarData.cpp +++ b/src/Evolution/DiscontinuousGalerkin/MortarData.cpp @@ -12,9 +12,11 @@ #include #include +#include "DataStructures/ApplyMatrices.hpp" #include "DataStructures/DataVector.hpp" #include "DataStructures/Tensor/Tensor.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "NumericalAlgorithms/Spectral/Projection.hpp" #include "Utilities/ErrorHandling/Assert.hpp" #include "Utilities/GenerateInstantiations.hpp" #include "Utilities/Serialization/PupStlCpp17.hpp" @@ -26,7 +28,76 @@ void MortarData::pup(PUP::er& p) { p | face_normal_magnitude; p | face_det_jacobian; p | volume_det_inv_jacobian; + p | mortar_mesh; p | face_mesh; + p | volume_mesh; +} + +template +void p_project( + const gsl::not_null<::evolution::dg::MortarData*> mortar_data, + const Mesh& new_mortar_mesh, const Mesh& new_face_mesh, + const Mesh& new_volume_mesh) { + // nothing needs to be done in 1D as mortars/faces are a single point... + if constexpr (Dim > 1) { + if (mortar_data->mortar_data.has_value()) { + const auto& old_mortar_mesh = mortar_data->mortar_mesh.value(); + if (old_mortar_mesh != new_mortar_mesh) { + const auto mortar_projection_matrices = + Spectral::p_projection_matrices(old_mortar_mesh, new_mortar_mesh); + DataVector& vars = mortar_data->mortar_data.value(); + vars = apply_matrices(mortar_projection_matrices, vars, + old_mortar_mesh.extents()); + mortar_data->mortar_mesh = new_mortar_mesh; + } + } + if (mortar_data->face_normal_magnitude.has_value()) { + const auto& old_face_mesh = mortar_data->face_mesh.value(); + if (old_face_mesh != new_face_mesh) { + const auto face_projection_matrices = + Spectral::p_projection_matrices(old_face_mesh, new_face_mesh); + DataVector& n = get(mortar_data->face_normal_magnitude.value()); + n = apply_matrices(face_projection_matrices, n, + old_face_mesh.extents()); + if (mortar_data->face_det_jacobian.has_value()) { + DataVector& det_j = get(mortar_data->face_det_jacobian.value()); + det_j = apply_matrices(face_projection_matrices, det_j, + old_face_mesh.extents()); + } + mortar_data->face_mesh = new_face_mesh; + } + } + } + if (mortar_data->volume_det_inv_jacobian.has_value()) { + const auto& old_volume_mesh = mortar_data->volume_mesh.value(); + if (old_volume_mesh != new_volume_mesh) { + const auto volume_projection_matrices = + Spectral::p_projection_matrices(old_volume_mesh, new_volume_mesh); + DataVector& det_inv_j = get(mortar_data->volume_det_inv_jacobian.value()); + det_inv_j = apply_matrices(volume_projection_matrices, det_inv_j, + old_volume_mesh.extents()); + mortar_data->volume_mesh = new_volume_mesh; + } + } +} + +template +void p_project_only_mortar_data( + const gsl::not_null<::evolution::dg::MortarData*> mortar_data, + const Mesh& new_mortar_mesh) { + // nothing needs to be done in 1D as mortars are a single point... + if constexpr (Dim > 1) { + const auto& old_mortar_mesh = mortar_data->mortar_mesh.value(); + const auto mortar_projection_matrices = + Spectral::p_projection_matrices(old_mortar_mesh, new_mortar_mesh); + DataVector& vars = mortar_data->mortar_data.value(); + vars = apply_matrices(mortar_projection_matrices, vars, + old_mortar_mesh.extents()); + mortar_data->mortar_mesh = new_mortar_mesh; + } else { + (void)mortar_data; + (void)new_mortar_mesh; + } } template @@ -35,7 +106,8 @@ bool operator==(const MortarData& lhs, const MortarData& rhs) { lhs.face_normal_magnitude == rhs.face_normal_magnitude and lhs.face_det_jacobian == rhs.face_det_jacobian and lhs.volume_det_inv_jacobian == rhs.volume_det_inv_jacobian and - lhs.face_mesh == rhs.face_mesh; + lhs.mortar_mesh == rhs.mortar_mesh and + lhs.face_mesh == rhs.face_mesh and lhs.volume_mesh == rhs.volume_mesh; } template @@ -46,22 +118,34 @@ bool operator!=(const MortarData& lhs, const MortarData& rhs) { template std::ostream& operator<<(std::ostream& os, const MortarData& mortar_data) { os << "Mortar data: " << mortar_data.mortar_data << "\n"; + os << "Mortar mesh: " << mortar_data.mortar_mesh << "\n"; os << "Face normal magnitude: " << mortar_data.face_normal_magnitude << "\n"; os << "Face det(J): " << mortar_data.face_det_jacobian << "\n"; os << "Face mesh: " << mortar_data.face_mesh << "\n"; os << "Volume det(invJ): " << mortar_data.volume_det_inv_jacobian << "\n"; + os << "Volume mesh: " << mortar_data.volume_mesh << "\n"; return os; } #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) -#define INSTANTIATION(r, data) \ - template class MortarData; \ - template bool operator==(const MortarData& lhs, \ - const MortarData& rhs); \ - template bool operator!=(const MortarData& lhs, \ - const MortarData& rhs); \ - template std::ostream& operator<<(std::ostream& os, \ +#define INSTANTIATION(r, data) \ + template class MortarData; \ + template void p_project( \ + const gsl::not_null<::evolution::dg::MortarData*> \ + mortar_data, \ + const Mesh& new_mortar_mesh, \ + const Mesh& new_face_mesh, \ + const Mesh& volume_mesh); \ + template void p_project_only_mortar_data( \ + const gsl::not_null<::evolution::dg::MortarData*> \ + mortar_data, \ + const Mesh& new_mortar_mesh); \ + template bool operator==(const MortarData& lhs, \ + const MortarData& rhs); \ + template bool operator!=(const MortarData& lhs, \ + const MortarData& rhs); \ + template std::ostream& operator<<(std::ostream& os, \ const MortarData& mortar_data); GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) diff --git a/src/Evolution/DiscontinuousGalerkin/MortarData.hpp b/src/Evolution/DiscontinuousGalerkin/MortarData.hpp index 8cbeb5159134..228d175818e6 100644 --- a/src/Evolution/DiscontinuousGalerkin/MortarData.hpp +++ b/src/Evolution/DiscontinuousGalerkin/MortarData.hpp @@ -45,8 +45,10 @@ namespace evolution::dg { * In addition, for local time stepping with Gauss points, the determinants of * the volume inverse Jacobian and the face Jacobian are stored. * - * In addition to the (type-erased) fields on the mortar, the face mesh is - * stored + * In addition to the (type-erased) fields on the mortar, the appropriate meshes + * are stored. When setting the mortar data and mortar mesh, the face mesh + * should also be set as it is used when hybridizing DG with finite difference + * or finite volume schemes (DG-subcell). * * If the element and its neighbor have unaligned logical coordinate * systems then the data and meshes are stored in the local logical @@ -62,12 +64,38 @@ struct MortarData { std::optional> face_normal_magnitude{std::nullopt}; std::optional> face_det_jacobian{std::nullopt}; std::optional> volume_det_inv_jacobian{std::nullopt}; + std::optional> mortar_mesh{std::nullopt}; std::optional> face_mesh{std::nullopt}; + std::optional> volume_mesh{std::nullopt}; // NOLINTNEXTLINE(google-runtime-references) void pup(PUP::er& p); }; +/// \brief Projects the mortar data when p-refined +/// +/// \details only updates the stored mesh if the corresponding data exists +/// +/// \note The DG-subcell code stores the face mesh in the MortarData even when +/// the geometric data are not used. Currently projection of MortarData is only +/// needed for local time-stepping which is not yet supported for DG-subcell. +template +void p_project(gsl::not_null<::evolution::dg::MortarData*> mortar_data, + const Mesh& new_mortar_mesh, + const Mesh& new_face_mesh, + const Mesh& new_volume_mesh); + +/// \brief Projects the mortar data (but not the geometric data) when p-refined +/// +/// \details Used to re-project mortar data when the mortar mesh changes +/// reactively after the neighbor face mesh is received. In this case, the +/// geometric data does not need to be updated as it already used the correct +/// face/volume mesh. +template +void p_project_only_mortar_data( + gsl::not_null<::evolution::dg::MortarData*> mortar_data, + const Mesh& new_mortar_mesh); + template bool operator==(const MortarData& lhs, const MortarData& rhs); diff --git a/src/Evolution/DiscontinuousGalerkin/MortarDataHolder.cpp b/src/Evolution/DiscontinuousGalerkin/MortarDataHolder.cpp index dae9f93cd112..1f11fbc30d5f 100644 --- a/src/Evolution/DiscontinuousGalerkin/MortarDataHolder.cpp +++ b/src/Evolution/DiscontinuousGalerkin/MortarDataHolder.cpp @@ -3,6 +3,8 @@ #include "Evolution/DiscontinuousGalerkin/MortarDataHolder.hpp" +#include +#include #include #include "Utilities/GenerateInstantiations.hpp" @@ -14,9 +16,37 @@ void MortarDataHolder::pup(PUP::er& p) { p | neighbor_data_; } +template +bool operator==(const MortarDataHolder& lhs, + const MortarDataHolder& rhs) { + return lhs.local() == rhs.local() and lhs.neighbor() == rhs.neighbor(); +} + +template +bool operator!=(const MortarDataHolder& lhs, + const MortarDataHolder& rhs) { + return not(lhs == rhs); +} + +template +std::ostream& operator<<(std::ostream& os, + const MortarDataHolder& mortar_data_holder) { + os << "Local mortar data:\n" << mortar_data_holder.local() << "\n"; + os << "Neighbor mortar data:\n" << mortar_data_holder.neighbor() << "\n"; + return os; +} + #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) -#define INSTANTIATION(r, data) template class MortarDataHolder; +#define INSTANTIATION(r, data) \ + template class MortarDataHolder; \ + template bool operator==(const MortarDataHolder& lhs, \ + const MortarDataHolder& rhs); \ + template bool operator!=(const MortarDataHolder& lhs, \ + const MortarDataHolder& rhs); \ + template std::ostream& operator<<( \ + std::ostream& os, \ + const MortarDataHolder& mortar_data_holder); GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) diff --git a/src/Evolution/DiscontinuousGalerkin/MortarDataHolder.hpp b/src/Evolution/DiscontinuousGalerkin/MortarDataHolder.hpp index 1801241fcf37..d5b61ef89943 100644 --- a/src/Evolution/DiscontinuousGalerkin/MortarDataHolder.hpp +++ b/src/Evolution/DiscontinuousGalerkin/MortarDataHolder.hpp @@ -3,6 +3,9 @@ #pragma once +#include +#include + #include "Evolution/DiscontinuousGalerkin/MortarData.hpp" /// \cond @@ -33,7 +36,19 @@ class MortarDataHolder { void pup(PUP::er& p); private: - MortarData local_data_; - MortarData neighbor_data_; + MortarData local_data_{}; + MortarData neighbor_data_{}; }; + +template +bool operator==(const MortarDataHolder& lhs, + const MortarDataHolder& rhs); + +template +bool operator!=(const MortarDataHolder& lhs, + const MortarDataHolder& rhs); + +template +std::ostream& operator<<(std::ostream& os, + const MortarDataHolder& mortar_data_holder); } // namespace evolution::dg diff --git a/tests/Unit/Domain/Structure/Test_OrientationMapHelpers.cpp b/tests/Unit/Domain/Structure/Test_OrientationMapHelpers.cpp index 754e69bd9ca7..a8023a19e21f 100644 --- a/tests/Unit/Domain/Structure/Test_OrientationMapHelpers.cpp +++ b/tests/Unit/Domain/Structure/Test_OrientationMapHelpers.cpp @@ -25,6 +25,7 @@ #include "Domain/Structure/OrientationMap.hpp" #include "Domain/Structure/OrientationMapHelpers.hpp" #include "Domain/Structure/Side.hpp" +#include "Framework/TestHelpers.hpp" #include "NumericalAlgorithms/Spectral/Basis.hpp" #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" @@ -401,10 +402,13 @@ void test_0d_orient_variables_on_slice() { Variables>> vars(slice_extents.product()); get(get(vars)) = DataVector{{-0.5}}; get<0>(get>(vars)) = DataVector{{1.0}}; + const auto slice_mesh = Mesh<0>{}; for (const auto& side : {Side::Lower, Side::Upper}) { CAPTURE(side); const OrientationMap<1> orientation_map( std::array, 1>{{Direction<1>(0, side)}}); + CHECK(slice_mesh == + orient_mesh_on_slice(slice_mesh, 0_st, orientation_map)); const auto oriented_vars = orient_variables_on_slice(vars, slice_extents, 0, orientation_map); // 1D boundary is a point, so no change expected @@ -425,6 +429,8 @@ void test_1d_slice_with_orientation(const OrientationMap<2>& orientation_map) { for (const size_t sliced_dim : {0_st, 1_st}) { CAPTURE(sliced_dim); + CHECK(slice_mesh == + orient_mesh_on_slice(slice_mesh, sliced_dim, orientation_map)); // Because `orientation_map` transforms between directions in the volume, we // make a new OrientationMap that transforms between directions on the // slices. In the case of a 1D slice, this is a 1D OrientationMap. @@ -506,10 +512,14 @@ void test_1d_orient_variables_on_slice() { } void test_2d_slice_with_orientation(const OrientationMap<3>& orientation_map) { - const auto slice_extents = Index<2>{3, 4}; - const auto slice_mesh = - Mesh<2>(slice_extents.indices(), Spectral::Basis::Legendre, - Spectral::Quadrature::GaussLobatto); + const auto volume_extents = Index<3>{3, 4, 5}; + const auto volume_mesh = + Mesh<3>(volume_extents.indices(), + {{Spectral::Basis::Chebyshev, Spectral::Basis::Legendre, + Spectral::Basis::FiniteDifference}}, + {{Spectral::Quadrature::GaussLobatto, Spectral::Quadrature::Gauss, + Spectral::Quadrature::CellCentered}}); + const auto oriented_mesh = orientation_map(volume_mesh); const auto affine = Affine2D{Affine(-1.0, 1.0, 2.3, 4.5), Affine(-1.0, 1.0, 0.8, 3.1)}; const auto map = @@ -518,6 +528,14 @@ void test_2d_slice_with_orientation(const OrientationMap<3>& orientation_map) { for (const size_t sliced_dim : {0_st, 1_st, 2_st}) { CAPTURE(sliced_dim); + const auto slice_mesh = volume_mesh.slice_away(sliced_dim); + const auto& slice_extents = slice_mesh.extents(); + CAPTURE(slice_mesh); + CAPTURE(slice_extents); + const auto oriented_sliced_dim = orientation_map(sliced_dim); + const auto expected_mesh = oriented_mesh.slice_away(oriented_sliced_dim); + CHECK(expected_mesh == + orient_mesh_on_slice(slice_mesh, sliced_dim, orientation_map)); // Because `orientation_map` transforms between directions in the volume, we // make a new OrientationMap that transforms between directions on the @@ -567,7 +585,7 @@ void test_2d_slice_with_orientation(const OrientationMap<3>& orientation_map) { get(get(expected_vars)) = oriented_mapped_coords[0]; get<0>(get>(expected_vars)) = oriented_mapped_coords[0]; get<1>(get>(expected_vars)) = oriented_mapped_coords[1]; - CHECK(oriented_vars == expected_vars); + CHECK_VARIABLES_APPROX(oriented_vars, expected_vars); check_vector(vars, oriented_vars, slice_extents, sliced_dim, orientation_map); diff --git a/tests/Unit/Evolution/DgSubcell/Actions/Test_ReconstructionCommunication.cpp b/tests/Unit/Evolution/DgSubcell/Actions/Test_ReconstructionCommunication.cpp index c9ee2e92080e..bb101d076492 100644 --- a/tests/Unit/Evolution/DgSubcell/Actions/Test_ReconstructionCommunication.cpp +++ b/tests/Unit/Evolution/DgSubcell/Actions/Test_ReconstructionCommunication.cpp @@ -97,7 +97,8 @@ struct component { evolution::dg::subcell::Tags::DataForRdmpTci, evolution::dg::subcell::Tags::TciDecision, evolution::dg::subcell::Tags::NeighborTciDecisions, - ::Tags::Variables>, evolution::dg::Tags::MortarData, + ::Tags::Variables>, evolution::dg::Tags::MortarMesh, + evolution::dg::Tags::MortarData, evolution::dg::Tags::MortarNextTemporalId, domain::Tags::NeighborMesh, evolution::dg::subcell::Tags::CellCenteredFlux, Dim>, @@ -258,16 +259,20 @@ void test(const bool use_cell_centered_flux) { } using MortarData = typename evolution::dg::Tags::MortarData::type; + using MortarMesh = typename evolution::dg::Tags::MortarMesh::type; using MortarNextId = typename evolution::dg::Tags::MortarNextTemporalId::type; MortarData mortar_data{}; + MortarMesh mortar_mesh{}; MortarNextId mortar_next_id{}; mortar_data[east_neighbor_id] = {}; + mortar_mesh[east_neighbor_id] = {}; mortar_next_id[east_neighbor_id] = {}; if constexpr (Dim > 1) { const DirectionalId south_neighbor_id{Direction::lower_eta(), south_id}; mortar_data[south_neighbor_id] = {}; + mortar_mesh[south_neighbor_id] = {}; mortar_next_id[south_neighbor_id] = {}; } @@ -291,9 +296,9 @@ void test(const bool use_cell_centered_flux) { evolution::dg::subcell::RdmpTciData{}, neighbor_tci_decision, typename evolution::dg::subcell::Tags::NeighborTciDecisions< Dim>::type{}, - Variables{}, MortarData{}, MortarNextId{}, - typename domain::Tags::NeighborMesh::type{}, cell_centered_flux, - Interps{}, Interps{}}); + Variables{}, MortarMesh{}, MortarData{}, + MortarNextId{}, typename domain::Tags::NeighborMesh::type{}, + cell_centered_flux, Interps{}, Interps{}}); ++neighbor_tci_decision; } } @@ -307,10 +312,10 @@ void test(const bool use_cell_centered_flux) { // Explicitly set RDMP data since this would be set previously by the TCI evolution::dg::subcell::RdmpTciData{{max(get(get(evolved_vars)))}, {min(get(get(evolved_vars)))}}, - self_tci_decision, neighbor_decision, evolved_vars, mortar_data, - mortar_next_id, typename domain::Tags::NeighborMesh::type{}, - cell_centered_flux, fd_to_neighbor_fd_interpolants, - neighbor_dg_to_fd_interpolants}); + self_tci_decision, neighbor_decision, evolved_vars, mortar_mesh, + mortar_data, mortar_next_id, + typename domain::Tags::NeighborMesh::type{}, cell_centered_flux, + fd_to_neighbor_fd_interpolants, neighbor_dg_to_fd_interpolants}); using ghost_data_tag = evolution::dg::subcell::Tags::GhostDataForReconstruction; diff --git a/tests/Unit/Evolution/DgSubcell/Test_CorrectPackagedData.cpp b/tests/Unit/Evolution/DgSubcell/Test_CorrectPackagedData.cpp index 8d748d2539b8..ce3f2068e93d 100644 --- a/tests/Unit/Evolution/DgSubcell/Test_CorrectPackagedData.cpp +++ b/tests/Unit/Evolution/DgSubcell/Test_CorrectPackagedData.cpp @@ -119,8 +119,10 @@ void test() { // Insert neighbor DG data. upper_mortar_data.neighbor().face_mesh = dg_face_mesh; + upper_mortar_data.neighbor().mortar_mesh = dg_face_mesh; upper_mortar_data.neighbor().mortar_data = upper_neighbor_data; lower_mortar_data.neighbor().face_mesh = dg_face_mesh; + lower_mortar_data.neighbor().mortar_mesh = dg_face_mesh; lower_mortar_data.neighbor().mortar_data = lower_neighbor_data; const Mesh subcell_face_mesh = @@ -225,8 +227,10 @@ void test() { dg_face_mesh.number_of_grid_points()}; std::iota(lower_local_data.begin(), lower_local_data.end(), 1.0e7); upper_mortar_data.local().face_mesh = dg_face_mesh; + upper_mortar_data.local().mortar_mesh = dg_face_mesh; upper_mortar_data.local().mortar_data = upper_local_data; lower_mortar_data.local().face_mesh = dg_face_mesh; + lower_mortar_data.local().mortar_mesh = dg_face_mesh; lower_mortar_data.local().mortar_data = lower_local_data; evolution::dg::subcell::correct_package_data( diff --git a/tests/Unit/Evolution/DgSubcell/Test_NeighborReconstructedFaceSolution.cpp b/tests/Unit/Evolution/DgSubcell/Test_NeighborReconstructedFaceSolution.cpp index eb45d48ff3c4..27fe7ac968c2 100644 --- a/tests/Unit/Evolution/DgSubcell/Test_NeighborReconstructedFaceSolution.cpp +++ b/tests/Unit/Evolution/DgSubcell/Test_NeighborReconstructedFaceSolution.cpp @@ -40,10 +40,10 @@ template using NeighborReconstructionMap = DirectionalIdMap; template -using MortarData = evolution::dg::BoundaryData; +using BoundaryData = evolution::dg::BoundaryData; template -using MortarDataMap = DirectionalIdMap>; +using BoundaryDataMap = DirectionalIdMap>; template struct Metavariables { @@ -85,7 +85,7 @@ void test() { evolution::dg::subcell::Tags::DataForRdmpTci, VolumeDouble>>( std::move(neighbor_data_map), std::move(rdmp_tci_data), 2.5); - std::pair> mortar_data_from_neighbors{}; + std::pair> mortar_data_from_neighbors{}; for (size_t d = 0; d < Dim; ++d) { const Mesh dg_volume_mesh{2 + 2 * Dim, Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto}; @@ -109,21 +109,25 @@ void test() { if (d % 2 == 0) { mortar_data_from_neighbors.second[DirectionalId{ Direction{d, Side::Upper}, ElementId{2 * d}}] = - MortarData{dg_volume_mesh, dg_face_mesh, dg_recons_and_rdmp_data, - dg_flux_data, {}, 1}; + BoundaryData{ + dg_volume_mesh, dg_face_mesh, dg_recons_and_rdmp_data, + dg_flux_data, {}, 1}; mortar_data_from_neighbors.second[DirectionalId{ Direction{d, Side::Lower}, ElementId{2 * d + 1}}] = - MortarData{fd_volume_mesh, fd_face_mesh, fd_recons_and_rdmp_data, - std::nullopt, {}, 2}; + BoundaryData{ + fd_volume_mesh, fd_face_mesh, fd_recons_and_rdmp_data, + std::nullopt, {}, 2}; } else { mortar_data_from_neighbors.second[DirectionalId{ Direction{d, Side::Lower}, ElementId{2 * d}}] = - MortarData{dg_volume_mesh, dg_face_mesh, dg_recons_and_rdmp_data, - dg_flux_data, {}, 3}; + BoundaryData{ + dg_volume_mesh, dg_face_mesh, dg_recons_and_rdmp_data, + dg_flux_data, {}, 3}; mortar_data_from_neighbors.second[DirectionalId{ Direction{d, Side::Upper}, ElementId{2 * d + 1}}] = - MortarData{fd_volume_mesh, fd_face_mesh, fd_recons_and_rdmp_data, - std::nullopt, {}, 4}; + BoundaryData{ + fd_volume_mesh, fd_face_mesh, fd_recons_and_rdmp_data, + std::nullopt, {}, 4}; } } evolution::dg::subcell::neighbor_reconstructed_face_solution< diff --git a/tests/Unit/Evolution/DiscontinuousGalerkin/Actions/Test_ApplyBoundaryCorrections.cpp b/tests/Unit/Evolution/DiscontinuousGalerkin/Actions/Test_ApplyBoundaryCorrections.cpp index 3b22c5485a2b..a340d49e37f1 100644 --- a/tests/Unit/Evolution/DiscontinuousGalerkin/Actions/Test_ApplyBoundaryCorrections.cpp +++ b/tests/Unit/Evolution/DiscontinuousGalerkin/Actions/Test_ApplyBoundaryCorrections.cpp @@ -231,13 +231,14 @@ struct SetLocalMortarData { 100 * count + 1000); db::mutate>( - [&face_mesh, &mortar_id, + [&face_mesh, &mortar_id, &mortar_mesh, &type_erased_boundary_data_on_mortar](const auto mortar_data_ptr) { // when using local time stepping, we reset the local mortar data // at the end of the SetLocalMortarData action since the // ComputeTimeDerivative action would've moved the data into the // boundary history. mortar_data_ptr->at(mortar_id).local().face_mesh = face_mesh; + mortar_data_ptr->at(mortar_id).local().mortar_mesh = mortar_mesh; mortar_data_ptr->at(mortar_id).local().mortar_data = std::move(type_erased_boundary_data_on_mortar); }, @@ -268,6 +269,7 @@ struct SetLocalMortarData { evolution::dg::MortarData past_mortar_data{}; past_mortar_data.face_mesh = face_mesh; + past_mortar_data.mortar_mesh = mortar_mesh; past_mortar_data.mortar_data = std::move(type_erased_boundary_data_on_mortar); Scalar local_face_normal_magnitude{ @@ -292,6 +294,7 @@ struct SetLocalMortarData { 100 * count + 300000); past_mortar_data.volume_det_inv_jacobian = local_volume_det_inv_jacobian; + past_mortar_data.volume_mesh = volume_mesh; past_mortar_data.face_det_jacobian = local_face_det_jacobian; } using dt_variables_tag = @@ -301,7 +304,7 @@ struct SetLocalMortarData { evolution::dg::Tags::MortarData, evolution::dg::Tags::MortarDataHistory< Metavariables::volume_dim, typename dt_variables_tag::type>>( - [&det_inv_jacobian, &mortar_id, &past_mortar_data, + [&det_inv_jacobian, &mortar_id, &volume_mesh, &past_mortar_data, &past_time_step_id, &time_step_id]( const auto mortar_data_ptr, const auto mortar_data_history_ptr, @@ -351,6 +354,7 @@ struct SetLocalMortarData { interpolation_matrices, get(det_jacobian), mesh.extents()); local_mortar_data.volume_det_inv_jacobian = det_inv_jacobian; + local_mortar_data.volume_mesh = volume_mesh; local_mortar_data.face_det_jacobian = face_det_jacobian; } mortar_data_history_ptr->at(mortar_id).local().insert( @@ -694,6 +698,7 @@ void test_impl(const Spectral::Quadrature quadrature, if (neighbor_time_step_id < local_next_time_step_id) { evolution::dg::MortarData nhbr_mortar_data{}; nhbr_mortar_data.face_mesh = face_mesh; + nhbr_mortar_data.mortar_mesh = mortar_mesh; nhbr_mortar_data.mortar_data = flux_data; mortar_data_history.at(mortar_id).remote().insert( neighbor_time_step_id, integration_order, @@ -701,6 +706,7 @@ void test_impl(const Spectral::Quadrature quadrature, } } else { all_mortar_data.at(mortar_id).neighbor().face_mesh = face_mesh; + all_mortar_data.at(mortar_id).neighbor().mortar_mesh = mortar_mesh; all_mortar_data.at(mortar_id).neighbor().mortar_data = flux_data; } ++count; diff --git a/tests/Unit/Evolution/DiscontinuousGalerkin/CMakeLists.txt b/tests/Unit/Evolution/DiscontinuousGalerkin/CMakeLists.txt index ac5fd7ced1df..b79f3de98fc5 100644 --- a/tests/Unit/Evolution/DiscontinuousGalerkin/CMakeLists.txt +++ b/tests/Unit/Evolution/DiscontinuousGalerkin/CMakeLists.txt @@ -18,6 +18,7 @@ set(LIBRARY_SOURCES Test_BoundaryCorrectionsHelper.cpp Test_BoundaryData.cpp Test_MortarData.cpp + Test_MortarDataHolder.cpp Test_MortarTags.cpp Test_NormalVectorTags.cpp Test_UsingSubcell.cpp @@ -28,6 +29,7 @@ add_test_library(${LIBRARY} "${LIBRARY_SOURCES}") target_link_libraries( ${LIBRARY} PRIVATE + Amr Boost::boost CoordinateMaps DataStructures diff --git a/tests/Unit/Evolution/DiscontinuousGalerkin/Initialization/Test_Mortars.cpp b/tests/Unit/Evolution/DiscontinuousGalerkin/Initialization/Test_Mortars.cpp index b57487e705e9..2ae6d03ccb19 100644 --- a/tests/Unit/Evolution/DiscontinuousGalerkin/Initialization/Test_Mortars.cpp +++ b/tests/Unit/Evolution/DiscontinuousGalerkin/Initialization/Test_Mortars.cpp @@ -358,6 +358,38 @@ struct Test<3, LocalTimeStepping> { } }; +template +void check_mortar_data(const MortarData& projected, + const MortarData& expected) { + CHECK(projected.mortar_mesh == expected.mortar_mesh); + CHECK(projected.face_mesh == expected.face_mesh); + CHECK(projected.volume_mesh == expected.volume_mesh); + if (projected.mortar_data.has_value()) { + CHECK_ITERABLE_APPROX(projected.mortar_data.value(), + expected.mortar_data.value()); + } else { + CHECK_FALSE(expected.mortar_data.has_value()); + } + if (projected.face_normal_magnitude.has_value()) { + CHECK_ITERABLE_APPROX(projected.face_normal_magnitude.value(), + expected.face_normal_magnitude.value()); + } else { + CHECK_FALSE(expected.face_normal_magnitude.has_value()); + } + if (projected.face_det_jacobian.has_value()) { + CHECK_ITERABLE_APPROX(projected.face_det_jacobian.value(), + expected.face_det_jacobian.value()); + } else { + CHECK_FALSE(expected.face_det_jacobian.has_value()); + } + if (projected.volume_det_inv_jacobian.has_value()) { + CHECK_ITERABLE_APPROX(projected.volume_det_inv_jacobian.value(), + expected.volume_det_inv_jacobian.value()); + } else { + CHECK_FALSE(expected.volume_det_inv_jacobian.has_value()); + } +} + template void test_p_refine( ::dg::MortarMap>& mortar_data, @@ -373,8 +405,8 @@ void test_p_refine( const Mesh& old_mesh, Mesh& new_mesh, const Element& old_element, Element& new_element, std::unordered_map, amr::Info>& neighbor_info, - const ::dg::MortarMap< - Dim, evolution::dg::MortarDataHolder>& /*expected_mortar_data*/, + const ::dg::MortarMap>& + expected_mortar_data, const ::dg::MortarMap>& expected_mortar_mesh, const ::dg::MortarMap>& expected_mortar_size, @@ -400,16 +432,30 @@ void test_p_refine( Metavariables>>(make_not_null(&box), std::make_pair(old_mesh, old_element)); - // Can't check the state as a default constructed MortarData has a default - // constructed TimeStepId which has nans which always compare as unequal - // CHECK(db::get>(box) == expected_mortar_data); + CHECK(db::get>(box) == expected_mortar_data); CHECK(db::get>(box) == expected_mortar_mesh); CHECK(db::get>(box) == expected_mortar_size); CHECK(db::get>(box) == expected_mortar_next_temporal_id); CHECK(db::get>(box) == expected_normal_covector_and_magnitude); - if (not UsingLts) { + if (UsingLts) { + const auto& projected_mortar_data_history = db::get< + Tags::MortarDataHistory::type>>( + box); + for (const auto& [mortar_id, history] : projected_mortar_data_history) { + for (size_t i = 0; i < 3; ++i) { + check_mortar_data( + history.local().data(i), + expected_mortar_data_history.at(mortar_id).local().data(i)); + } + for (size_t i = 0; i < 2; ++i) { + check_mortar_data( + history.remote().data(i), + expected_mortar_data_history.at(mortar_id).remote().data(i)); + } + } + } else { (void)(expected_mortar_data_history); CHECK( db::get< @@ -478,8 +524,6 @@ void test_p_refine_gts() { auto new_element = make_element(); const TimeStepId next_temporal_id{true, 3, Time{Slab{0.2, 3.4}, {6, 100}}}; - // These quantities are re-allocated after projection, so we can - // just set them to empty maps... ::dg::MortarMap> mortar_data{}; ::dg::MortarMap> mortar_mesh{}; ::dg::MortarMap> mortar_size{}; @@ -492,8 +536,18 @@ void test_p_refine_gts() { ::dg::MortarMap mortar_next_temporal_ids{}; std::unordered_map, amr::Info> neighbor_info{}; for (const auto& [direction, neighbors] : old_element.neighbors()) { + normal_covector_and_magnitude[direction] = std::nullopt; for (const auto& neighbor : neighbors) { const DirectionalId mortar_id{direction, neighbor}; + mortar_data.emplace(mortar_id, MortarDataHolder{}); + mortar_mesh.emplace( + mortar_id, + ::dg::mortar_mesh(old_mesh.slice_away(direction.dimension()), + neighbor_mesh.slice_away(direction.dimension()))); + mortar_size.emplace( + mortar_id, + ::dg::mortar_size(old_element.id(), neighbor, direction.dimension(), + neighbors.orientation())); mortar_next_temporal_ids.emplace(mortar_id, next_temporal_id); neighbor_info.emplace( neighbor, @@ -530,6 +584,7 @@ void test_p_refine_gts() { } } for (const auto& direction : new_element.external_boundaries()) { + normal_covector_and_magnitude[direction] = std::nullopt; expected_normal_covector_and_magnitude[direction] = std::nullopt; } @@ -542,6 +597,158 @@ void test_p_refine_gts() { expected_mortar_data_history); } +template +MortarData make_mortar_data(const Mesh& mortar_mesh, + const Mesh& face_mesh, + const Mesh& volume_mesh, + const bool is_local_side, const double value) { + constexpr size_t number_of_components = 1 + Dim; + MortarData mortar_data; + mortar_data.mortar_data = DataVector{ + mortar_mesh.number_of_grid_points() * number_of_components, value}; + mortar_data.mortar_mesh = mortar_mesh; + if (is_local_side) { + mortar_data.face_normal_magnitude = Scalar{ + DataVector{face_mesh.number_of_grid_points(), 2.0 * value}}; + mortar_data.face_det_jacobian = Scalar{ + DataVector{face_mesh.number_of_grid_points(), 3.0 * value}}; + mortar_data.face_mesh = face_mesh; + mortar_data.volume_det_inv_jacobian = Scalar{ + DataVector{volume_mesh.number_of_grid_points(), 4.0 * value}}; + mortar_data.volume_mesh = volume_mesh; + } + return mortar_data; +} + +template +using boundary_history_type = + typename mortar_data_history_type::mapped_type; + +template +void test_p_refine_lts() { + const Mesh old_mesh{4, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + Mesh new_mesh{5, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + const Mesh neighbor_mesh{3, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + + const auto old_element = make_element(); + auto new_element = make_element(); + const TimeStepId next_temporal_id{true, 3, Time{Slab{0.2, 3.4}, {6, 100}}}; + const std::vector local_past_ids{ + TimeStepId{true, 3, Time{Slab{0.2, 3.4}, {0, 100}}}, + TimeStepId{true, 3, Time{Slab{0.2, 3.4}, {2, 100}}}, + TimeStepId{true, 3, Time{Slab{0.2, 3.4}, {4, 100}}}}; + const std::vector remote_past_ids{ + TimeStepId{true, 3, Time{Slab{0.2, 3.4}, {0, 100}}}, + TimeStepId{true, 3, Time{Slab{0.2, 3.4}, {3, 100}}}}; + + ::dg::MortarMap> mortar_data{}; + ::dg::MortarMap> mortar_mesh{}; + ::dg::MortarMap> mortar_size{}; + DirectionMap>>>> + normal_covector_and_magnitude{}; + mortar_data_history_type mortar_data_history{}; + ::dg::MortarMap mortar_next_temporal_ids{}; + std::unordered_map, amr::Info> neighbor_info{}; + for (const auto& [direction, neighbors] : old_element.neighbors()) { + normal_covector_and_magnitude[direction] = std::nullopt; + for (const auto& neighbor : neighbors) { + const DirectionalId mortar_id{direction, neighbor}; + mortar_data.emplace(mortar_id, MortarDataHolder{}); + mortar_mesh.emplace( + mortar_id, + ::dg::mortar_mesh(old_mesh.slice_away(direction.dimension()), + neighbor_mesh.slice_away(direction.dimension()))); + mortar_size.emplace( + mortar_id, + ::dg::mortar_size(old_element.id(), neighbor, direction.dimension(), + neighbors.orientation())); + mortar_next_temporal_ids.emplace(mortar_id, next_temporal_id); + neighbor_info.emplace( + neighbor, + amr::Info{std::array{}, neighbor_mesh}); + mortar_data_history.emplace(mortar_id, boundary_history_type{}); + for (size_t i = 0; i < 3; ++i) { + mortar_data_history.at(mortar_id).local().insert( + local_past_ids[i], 3, + make_mortar_data(mortar_mesh.at(mortar_id), + old_mesh.slice_away(direction.dimension()), + old_mesh, true, 5.0 + static_cast(i))); + } + for (size_t i = 0; i < 2; ++i) { + mortar_data_history.at(mortar_id).remote().insert( + local_past_ids[i], 3, + make_mortar_data(mortar_mesh.at(mortar_id), + neighbor_mesh.slice_away(direction.dimension()), + neighbor_mesh, false, + 5.0 + static_cast(i))); + } + } + } + + ::dg::MortarMap> + expected_mortar_data{}; + ::dg::MortarMap> expected_mortar_mesh{}; + ::dg::MortarMap> + expected_mortar_size{}; + ::dg::MortarMap expected_mortar_next_temporal_ids{}; + DirectionMap>>>> + expected_normal_covector_and_magnitude{}; + mortar_data_history_type expected_mortar_data_history{}; + for (const auto& [direction, neighbors] : new_element.neighbors()) { + expected_normal_covector_and_magnitude[direction] = std::nullopt; + for (const auto& neighbor : neighbors) { + const DirectionalId mortar_id{direction, neighbor}; + expected_mortar_data.emplace(mortar_id, MortarDataHolder{}); + expected_mortar_mesh.emplace( + mortar_id, + ::dg::mortar_mesh(new_mesh.slice_away(direction.dimension()), + neighbor_info.at(neighbor).new_mesh.slice_away( + direction.dimension()))); + expected_mortar_size.emplace( + mortar_id, + ::dg::mortar_size(new_element.id(), neighbor, direction.dimension(), + neighbors.orientation())); + expected_mortar_next_temporal_ids.emplace(mortar_id, next_temporal_id); + expected_mortar_data_history.emplace(mortar_id, + boundary_history_type{}); + for (size_t i = 0; i < 3; ++i) { + expected_mortar_data_history.at(mortar_id).local().insert( + local_past_ids[i], 3, + make_mortar_data(expected_mortar_mesh.at(mortar_id), + new_mesh.slice_away(direction.dimension()), + new_mesh, true, 5.0 + static_cast(i))); + } + for (size_t i = 0; i < 2; ++i) { + expected_mortar_data_history.at(mortar_id).remote().insert( + local_past_ids[i], 3, + make_mortar_data(expected_mortar_mesh.at(mortar_id), + neighbor_mesh.slice_away(direction.dimension()), + neighbor_mesh, false, + 5.0 + static_cast(i))); + } + } + } + for (const auto& direction : new_element.external_boundaries()) { + normal_covector_and_magnitude[direction] = std::nullopt; + expected_normal_covector_and_magnitude[direction] = std::nullopt; + } + + test_p_refine( + mortar_data, mortar_mesh, mortar_size, mortar_next_temporal_ids, + normal_covector_and_magnitude, mortar_data_history, old_mesh, new_mesh, + old_element, new_element, neighbor_info, expected_mortar_data, + expected_mortar_mesh, expected_mortar_size, + expected_mortar_next_temporal_ids, expected_normal_covector_and_magnitude, + expected_mortar_data_history); +} + } // namespace SPECTRE_TEST_CASE("Unit.Evolution.DG.Initialization.Mortars", @@ -563,5 +770,8 @@ SPECTRE_TEST_CASE("Unit.Evolution.DG.Initialization.Mortars", test_p_refine_gts<1>(); test_p_refine_gts<2>(); test_p_refine_gts<3>(); + test_p_refine_lts<1>(); + test_p_refine_lts<2>(); + test_p_refine_lts<3>(); } } // namespace evolution::dg diff --git a/tests/Unit/Evolution/DiscontinuousGalerkin/Test_MortarData.cpp b/tests/Unit/Evolution/DiscontinuousGalerkin/Test_MortarData.cpp index 97f8c12a96e1..5539934d4426 100644 --- a/tests/Unit/Evolution/DiscontinuousGalerkin/Test_MortarData.cpp +++ b/tests/Unit/Evolution/DiscontinuousGalerkin/Test_MortarData.cpp @@ -28,10 +28,12 @@ namespace evolution::dg { namespace { template void assign_with_reference(const gsl::not_null*> mortar_data, + const Mesh& mortar_mesh, const Mesh& face_mesh, const std::optional& data, const std::string& expected_output) { if (data.has_value()) { + mortar_data->mortar_mesh = mortar_mesh; mortar_data->face_mesh = face_mesh; mortar_data->mortar_data = *data; CHECK(mortar_data->mortar_data.has_value()); @@ -81,26 +83,30 @@ void test_global_time_stepping_usage() { std::string local_expected_output = MakeString{} << "Mortar data: " << local_data << "\n" + << "Mortar mesh: " << mortar_mesh << "\n" << "Face normal magnitude: --\n" << "Face det(J): --\n" << "Face mesh: " << local_mesh << "\n" - << "Volume det(invJ): --\n"; + << "Volume det(invJ): --\n" + << "Volume mesh: --\n"; std::string neighbor_expected_output = MakeString{} << "Mortar data: " << neighbor_data << "\n" + << "Mortar mesh: " << mortar_mesh << "\n" << "Face normal magnitude: --\n" << "Face det(J): --\n" << "Face mesh: " << neighbor_mesh << "\n" - << "Volume det(invJ): --\n"; + << "Volume det(invJ): --\n" + << "Volume mesh: --\n"; CHECK_FALSE(local_mortar_data.mortar_data.has_value()); CHECK_FALSE(neighbor_mortar_data.mortar_data.has_value()); - assign_with_reference(make_not_null(&local_mortar_data), + assign_with_reference(make_not_null(&local_mortar_data), mortar_mesh, local_mesh, std::optional{local_data}, local_expected_output); - assign_with_reference(make_not_null(&neighbor_mortar_data), + assign_with_reference(make_not_null(&neighbor_mortar_data), mortar_mesh, neighbor_mesh, std::optional{neighbor_data}, neighbor_expected_output); @@ -129,12 +135,14 @@ void test_local_time_stepping_usage(const bool use_gauss_points) { std::string expected_output = MakeString{} << "Mortar data: " << local_data << "\n" + << "Mortar mesh: " << mortar_mesh << "\n" << "Face normal magnitude: --\n" << "Face det(J): --\n" << "Face mesh: " << local_mesh << "\n" - << "Volume det(invJ): --\n"; + << "Volume det(invJ): --\n" + << "Volume mesh: --\n"; - assign_with_reference(make_not_null(&mortar_data), local_mesh, + assign_with_reference(make_not_null(&mortar_data), mortar_mesh, local_mesh, std::optional{local_data}, expected_output); const auto local_volume_det_inv_jacobian = @@ -183,11 +191,149 @@ void test_local_time_stepping_usage(const bool use_gauss_points) { CHECK_FALSE(mortar_data != deserialized_mortar_data); } +template +std::array make_extents(const size_t first_extent) { + if constexpr (1 == Dim) { + return {first_extent}; + } else if constexpr (2 == Dim) { + return {first_extent, first_extent + 1}; + } else if constexpr (3 == Dim) { + return {first_extent, first_extent + 1, first_extent + 2}; + } else { + return {}; + } +} + +template +void check_mortar_data(const MortarData& projected, + const MortarData& expected) { + CHECK(projected.mortar_mesh == expected.mortar_mesh); + CHECK(projected.face_mesh == expected.face_mesh); + CHECK(projected.volume_mesh == expected.volume_mesh); + if (projected.mortar_data.has_value()) { + CHECK_ITERABLE_APPROX(projected.mortar_data.value(), + expected.mortar_data.value()); + } else { + CHECK_FALSE(expected.mortar_data.has_value()); + } + if (projected.face_normal_magnitude.has_value()) { + CHECK_ITERABLE_APPROX(projected.face_normal_magnitude.value(), + expected.face_normal_magnitude.value()); + } else { + CHECK_FALSE(expected.face_normal_magnitude.has_value()); + } + if (projected.face_det_jacobian.has_value()) { + CHECK_ITERABLE_APPROX(projected.face_det_jacobian.value(), + expected.face_det_jacobian.value()); + } else { + CHECK_FALSE(expected.face_det_jacobian.has_value()); + } + if (projected.volume_det_inv_jacobian.has_value()) { + CHECK_ITERABLE_APPROX(projected.volume_det_inv_jacobian.value(), + expected.volume_det_inv_jacobian.value()); + } else { + CHECK_FALSE(expected.volume_det_inv_jacobian.has_value()); + } +} + +template +void test_p_project() { + const Mesh initial_volume_mesh{make_extents(3), + Spectral::Basis::Legendre, + Spectral::Quadrature::Gauss}; + const Mesh initial_face_mesh = + initial_volume_mesh.slice_away(Dim - 1); + const Mesh initial_mortar_mesh = initial_volume_mesh.slice_away(0); + const Mesh final_volume_mesh{make_extents(4), + Spectral::Basis::Legendre, + Spectral::Quadrature::Gauss}; + const Mesh final_face_mesh = final_volume_mesh.slice_away(Dim - 1); + const Mesh final_mortar_mesh = final_volume_mesh.slice_away(0); + constexpr size_t number_of_components = 1 + Dim; + const DataVector initial_mortar_data{ + initial_mortar_mesh.number_of_grid_points() * number_of_components, 1.0}; + const Scalar initial_face_normal_magnitude{ + DataVector{initial_face_mesh.number_of_grid_points(), 2.0}}; + const Scalar initial_face_det_jacobian{ + DataVector{initial_face_mesh.number_of_grid_points(), 3.0}}; + const Scalar initial_volume_det_inv_jacobian{ + DataVector{initial_volume_mesh.number_of_grid_points(), 4.0}}; + const DataVector final_mortar_data{ + final_mortar_mesh.number_of_grid_points() * number_of_components, 1.0}; + const Scalar final_face_normal_magnitude{ + DataVector{final_face_mesh.number_of_grid_points(), 2.0}}; + const Scalar final_face_det_jacobian{ + DataVector{final_face_mesh.number_of_grid_points(), 3.0}}; + const Scalar final_volume_det_inv_jacobian{ + DataVector{final_volume_mesh.number_of_grid_points(), 4.0}}; + MortarData only_mortar_data{std::optional(initial_mortar_data), + std::nullopt, + std::nullopt, + std::nullopt, + std::optional(initial_mortar_mesh), + std::nullopt, + std::nullopt}; + p_project(make_not_null(&only_mortar_data), final_mortar_mesh, + final_face_mesh, final_volume_mesh); + check_mortar_data(only_mortar_data, + MortarData{std::optional(final_mortar_data), + std::nullopt, std::nullopt, std::nullopt, + std::optional(final_mortar_mesh), + std::nullopt, std::nullopt}); + MortarData only_mortar_data_2{std::optional(initial_mortar_data), + std::nullopt, + std::nullopt, + std::nullopt, + std::optional(initial_mortar_mesh), + std::nullopt, + std::nullopt}; + p_project_only_mortar_data(make_not_null(&only_mortar_data_2), + final_mortar_mesh); + check_mortar_data(only_mortar_data_2, + MortarData{std::optional(final_mortar_data), + std::nullopt, std::nullopt, std::nullopt, + std::optional(final_mortar_mesh), + std::nullopt, std::nullopt}); + MortarData only_gl_data{std::optional(initial_mortar_data), + std::optional(initial_face_normal_magnitude), + std::nullopt, + std::nullopt, + std::optional(initial_mortar_mesh), + std::optional(initial_face_mesh), + std::nullopt}; + p_project(make_not_null(&only_gl_data), final_mortar_mesh, final_face_mesh, + final_volume_mesh); + check_mortar_data( + only_gl_data, + MortarData{std::optional(final_mortar_data), + std::optional(final_face_normal_magnitude), std::nullopt, + std::nullopt, std::optional(final_mortar_mesh), + std::optional(final_face_mesh), std::nullopt}); + MortarData g_data{std::optional(initial_mortar_data), + std::optional(initial_face_normal_magnitude), + std::optional(initial_face_det_jacobian), + std::optional(initial_volume_det_inv_jacobian), + std::optional(initial_mortar_mesh), + std::optional(initial_face_mesh), + std::optional(initial_volume_mesh)}; + p_project(make_not_null(&g_data), final_mortar_mesh, final_face_mesh, + final_volume_mesh); + check_mortar_data( + g_data, MortarData{std::optional(final_mortar_data), + std::optional(final_face_normal_magnitude), + std::optional(final_face_det_jacobian), + std::optional(final_volume_det_inv_jacobian), + std::optional(final_mortar_mesh), + std::optional(final_face_mesh), + std::optional(final_volume_mesh)}); +} + template void test() { test_global_time_stepping_usage(); test_local_time_stepping_usage(true); test_local_time_stepping_usage(false); + test_p_project(); } } // namespace diff --git a/tests/Unit/Evolution/DiscontinuousGalerkin/Test_MortarDataHolder.cpp b/tests/Unit/Evolution/DiscontinuousGalerkin/Test_MortarDataHolder.cpp new file mode 100644 index 000000000000..15c83cea8b40 --- /dev/null +++ b/tests/Unit/Evolution/DiscontinuousGalerkin/Test_MortarDataHolder.cpp @@ -0,0 +1,75 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include +#include + +#include "Evolution/DiscontinuousGalerkin/MortarData.hpp" +#include "Evolution/DiscontinuousGalerkin/MortarDataHolder.hpp" +#include "Framework/TestHelpers.hpp" +#include "NumericalAlgorithms/Spectral/Basis.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "NumericalAlgorithms/Spectral/Quadrature.hpp" +#include "Utilities/GetOutput.hpp" +#include "Utilities/MakeString.hpp" + +namespace evolution::dg { +namespace { +template +void test() { + const Mesh volume_mesh{4, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + const Mesh face_mesh{4, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + const Mesh mortar_mesh{5, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + constexpr size_t number_of_components = 1 + Dim; + DataVector local_mortar_data{ + mortar_mesh.number_of_grid_points() * number_of_components, 2.0}; + DataVector neighbor_mortar_data{ + mortar_mesh.number_of_grid_points() * number_of_components, 3.0}; + Scalar normal_magnitude{ + DataVector{face_mesh.number_of_grid_points(), 4.0}}; + Scalar det_j{DataVector{face_mesh.number_of_grid_points(), 5.0}}; + Scalar det_inv_j{ + DataVector{volume_mesh.number_of_grid_points(), 6.0}}; + MortarDataHolder holder{}; + holder.local() = MortarData{ + local_mortar_data, normal_magnitude, det_j, det_inv_j, + mortar_mesh, face_mesh, volume_mesh}; + holder.neighbor() = MortarData{ + neighbor_mortar_data, std::nullopt, std::nullopt, std::nullopt, + mortar_mesh, std::nullopt, std::nullopt}; + + std::string expected_output = + MakeString{} << "Local mortar data:\n" + << "Mortar data: " << local_mortar_data << "\n" + << "Mortar mesh: " << mortar_mesh << "\n" + << "Face normal magnitude: " << normal_magnitude << "\n" + << "Face det(J): " << det_j << "\n" + << "Face mesh: " << face_mesh << "\n" + << "Volume det(invJ): " << det_inv_j << "\n" + << "Volume mesh: " << volume_mesh << "\n\n" + << "Neighbor mortar data:\n" + << "Mortar data: " << neighbor_mortar_data << "\n" + << "Mortar mesh: " << mortar_mesh << "\n" + << "Face normal magnitude: --\n" + << "Face det(J): --\n" + << "Face mesh: --\n" + << "Volume det(invJ): --\n" + << "Volume mesh: --\n\n"; + + CHECK(get_output(holder) == expected_output); + const auto deserialized_holder = serialize_and_deserialize(holder); + CHECK(holder == deserialized_holder); +} +} // namespace + +SPECTRE_TEST_CASE("Unit.Evolution.DG.MortarDataHolder", "[Unit][Evolution]") { + test<1>(); + test<2>(); + test<3>(); +} +} // namespace evolution::dg