diff --git a/src/Evolution/DgSubcell/CMakeLists.txt b/src/Evolution/DgSubcell/CMakeLists.txt index cdb6fa5bfed9..41b45a09d086 100644 --- a/src/Evolution/DgSubcell/CMakeLists.txt +++ b/src/Evolution/DgSubcell/CMakeLists.txt @@ -25,6 +25,7 @@ spectre_target_headers( Mesh.hpp NeighborRdmpAndVolumeData.hpp NeighborReconstructedFaceSolution.hpp + NeighborReconstructedFaceSolution.tpp NeighborTciDecision.hpp PerssonTci.hpp PrepareNeighborData.hpp @@ -53,6 +54,7 @@ spectre_target_sources( Matrices.cpp Mesh.cpp NeighborRdmpAndVolumeData.cpp + NeighborTciDecision.cpp PerssonTci.cpp Projection.cpp RdmpTci.cpp diff --git a/src/Evolution/DgSubcell/NeighborReconstructedFaceSolution.hpp b/src/Evolution/DgSubcell/NeighborReconstructedFaceSolution.hpp index c6bfe41a0428..3dc2498794d7 100644 --- a/src/Evolution/DgSubcell/NeighborReconstructedFaceSolution.hpp +++ b/src/Evolution/DgSubcell/NeighborReconstructedFaceSolution.hpp @@ -4,26 +4,15 @@ #pragma once #include -#include #include #include #include -#include -#include "DataStructures/DataBox/DataBox.hpp" +#include "DataStructures/DataBox/Access.hpp" #include "DataStructures/DataVector.hpp" -#include "Domain/Structure/Direction.hpp" -#include "Domain/Structure/DirectionalId.hpp" #include "Domain/Structure/DirectionalIdMap.hpp" -#include "Domain/Structure/ElementId.hpp" -#include "Evolution/DgSubcell/GhostData.hpp" -#include "Evolution/DgSubcell/RdmpTciData.hpp" -#include "Evolution/DgSubcell/Tags/DataForRdmpTci.hpp" -#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" -#include "Evolution/DgSubcell/Tags/Mesh.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" #include "Time/TimeStepId.hpp" -#include "Utilities/ErrorHandling/Assert.hpp" #include "Utilities/Gsl.hpp" namespace evolution::dg::subcell { @@ -58,103 +47,15 @@ namespace evolution::dg::subcell { * the last argument to * `Metavariables::SubcellOptions::DgComputeSubcellNeighborPackagedData::apply`. */ -template +template void neighbor_reconstructed_face_solution( - const gsl::not_null*> box, - const gsl::not_null box, + gsl::not_null, Mesh, std::optional, std::optional, ::TimeStepId, int>>>*> - received_temporal_id_and_data) { - db::mutate, - subcell::Tags::DataForRdmpTci>( - [&received_temporal_id_and_data](const auto subcell_ghost_data_ptr, - const auto rdmp_tci_data_ptr) { - const size_t number_of_evolved_vars = - rdmp_tci_data_ptr->max_variables_values.size(); - for (auto& received_mortar_data : - received_temporal_id_and_data->second) { - const auto& mortar_id = received_mortar_data.first; - ASSERT(std::get<2>(received_mortar_data.second).has_value(), - "The subcell mortar data was not sent at TimeStepId " - << received_temporal_id_and_data->first - << " with mortar id " << mortar_id); - const DataVector& neighbor_ghost_and_subcell_data = - *std::get<2>(received_mortar_data.second); - // Compute min and max over neighbors - const size_t offset_for_min = - neighbor_ghost_and_subcell_data.size() - number_of_evolved_vars; - const size_t offset_for_max = offset_for_min - number_of_evolved_vars; - for (size_t var_index = 0; var_index < number_of_evolved_vars; - ++var_index) { - rdmp_tci_data_ptr->max_variables_values[var_index] = std::max( - rdmp_tci_data_ptr->max_variables_values[var_index], - neighbor_ghost_and_subcell_data[offset_for_max + var_index]); - rdmp_tci_data_ptr->min_variables_values[var_index] = std::min( - rdmp_tci_data_ptr->min_variables_values[var_index], - neighbor_ghost_and_subcell_data[offset_for_min + var_index]); - } - - ASSERT(subcell_ghost_data_ptr->find(mortar_id) == - subcell_ghost_data_ptr->end(), - "The subcell neighbor data is already inserted. Direction: " - << mortar_id.direction - << " with ElementId: " << mortar_id.id); - - (*subcell_ghost_data_ptr)[mortar_id] = GhostData{1}; - GhostData& all_ghost_data = subcell_ghost_data_ptr->at(mortar_id); - DataVector& neighbor_data = - all_ghost_data.neighbor_ghost_data_for_reconstruction(); - neighbor_data.destructive_resize( - neighbor_ghost_and_subcell_data.size() - - 2 * number_of_evolved_vars); - - // Copy over the neighbor data for reconstruction. We need this - // since we might be doing a step unwind and the DG algorithm deletes - // the inbox data after lifting the fluxes to the volume. - // The std::prev avoids copying over the data for the RDMP TCI, which - // is both the maximum and minimum of each evolved variable, so - // `2*number_of_evolved_vars` components. - std::copy(neighbor_ghost_and_subcell_data.begin(), - std::prev(neighbor_ghost_and_subcell_data.end(), - 2 * static_cast( - number_of_evolved_vars)), - neighbor_data.begin()); - } - }, - box); - std::vector> mortars_to_reconstruct_to{}; - for (auto& received_mortar_data : received_temporal_id_and_data->second) { - const auto& mortar_id = received_mortar_data.first; - if (not std::get<3>(received_mortar_data.second).has_value()) { - mortars_to_reconstruct_to.push_back(mortar_id); - } - } - DirectionalIdMap neighbor_reconstructed_evolved_vars = - DgComputeSubcellNeighborPackagedData::apply(*box, - mortars_to_reconstruct_to); - ASSERT(neighbor_reconstructed_evolved_vars.size() == - mortars_to_reconstruct_to.size(), - "Should have reconstructed " - << mortars_to_reconstruct_to.size() << " sides but reconstructed " - << neighbor_reconstructed_evolved_vars.size() << " sides."); - // Now copy over the packaged data _into_ the inbox in order to avoid having - // to make other changes to the DG algorithm (code in - // src/Evolution/DiscontinuousGalerkin) - for (auto& received_mortar_data : received_temporal_id_and_data->second) { - const auto& mortar_id = received_mortar_data.first; - if (not std::get<3>(received_mortar_data.second).has_value()) { - ASSERT(neighbor_reconstructed_evolved_vars.find(mortar_id) != - neighbor_reconstructed_evolved_vars.end(), - "Could not find mortar id " << mortar_id - << " in reconstructed data map."); - std::get<3>(received_mortar_data.second) = - std::move(neighbor_reconstructed_evolved_vars.at(mortar_id)); - } - } -} + received_temporal_id_and_data); } // namespace evolution::dg::subcell diff --git a/src/Evolution/DgSubcell/NeighborReconstructedFaceSolution.tpp b/src/Evolution/DgSubcell/NeighborReconstructedFaceSolution.tpp new file mode 100644 index 000000000000..67b5bd9e86f7 --- /dev/null +++ b/src/Evolution/DgSubcell/NeighborReconstructedFaceSolution.tpp @@ -0,0 +1,128 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "DataStructures/DataBox/Access.hpp" +#include "DataStructures/DataVector.hpp" +#include "Domain/Structure/Direction.hpp" +#include "Domain/Structure/DirectionalId.hpp" +#include "Domain/Structure/DirectionalIdMap.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Evolution/DgSubcell/GhostData.hpp" +#include "Evolution/DgSubcell/RdmpTciData.hpp" +#include "Evolution/DgSubcell/Tags/DataForRdmpTci.hpp" +#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" +#include "Evolution/DgSubcell/Tags/Mesh.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "Time/TimeStepId.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" +#include "Utilities/Gsl.hpp" + +namespace evolution::dg::subcell { +template +void neighbor_reconstructed_face_solution( + const gsl::not_null box, + const gsl::not_null, Mesh, + std::optional, std::optional, + ::TimeStepId, int>>>*> + received_temporal_id_and_data) { + db::mutate, + subcell::Tags::DataForRdmpTci>( + [&received_temporal_id_and_data](const auto subcell_ghost_data_ptr, + const auto rdmp_tci_data_ptr) { + const size_t number_of_evolved_vars = + rdmp_tci_data_ptr->max_variables_values.size(); + for (auto& received_mortar_data : + received_temporal_id_and_data->second) { + const auto& mortar_id = received_mortar_data.first; + ASSERT(std::get<2>(received_mortar_data.second).has_value(), + "The subcell mortar data was not sent at TimeStepId " + << received_temporal_id_and_data->first + << " with mortar id " << mortar_id); + const DataVector& neighbor_ghost_and_subcell_data = + *std::get<2>(received_mortar_data.second); + // Compute min and max over neighbors + const size_t offset_for_min = + neighbor_ghost_and_subcell_data.size() - number_of_evolved_vars; + const size_t offset_for_max = offset_for_min - number_of_evolved_vars; + for (size_t var_index = 0; var_index < number_of_evolved_vars; + ++var_index) { + rdmp_tci_data_ptr->max_variables_values[var_index] = std::max( + rdmp_tci_data_ptr->max_variables_values[var_index], + neighbor_ghost_and_subcell_data[offset_for_max + var_index]); + rdmp_tci_data_ptr->min_variables_values[var_index] = std::min( + rdmp_tci_data_ptr->min_variables_values[var_index], + neighbor_ghost_and_subcell_data[offset_for_min + var_index]); + } + + ASSERT(subcell_ghost_data_ptr->find(mortar_id) == + subcell_ghost_data_ptr->end(), + "The subcell neighbor data is already inserted. Direction: " + << mortar_id.direction + << " with ElementId: " << mortar_id.id); + + (*subcell_ghost_data_ptr)[mortar_id] = GhostData{1}; + GhostData& all_ghost_data = subcell_ghost_data_ptr->at(mortar_id); + DataVector& neighbor_data = + all_ghost_data.neighbor_ghost_data_for_reconstruction(); + neighbor_data.destructive_resize( + neighbor_ghost_and_subcell_data.size() - + 2 * number_of_evolved_vars); + + // Copy over the neighbor data for reconstruction. We need this + // since we might be doing a step unwind and the DG algorithm deletes + // the inbox data after lifting the fluxes to the volume. + // The std::prev avoids copying over the data for the RDMP TCI, which + // is both the maximum and minimum of each evolved variable, so + // `2*number_of_evolved_vars` components. + std::copy(neighbor_ghost_and_subcell_data.begin(), + std::prev(neighbor_ghost_and_subcell_data.end(), + 2 * static_cast( + number_of_evolved_vars)), + neighbor_data.begin()); + } + }, + box); + std::vector> mortars_to_reconstruct_to{}; + for (auto& received_mortar_data : received_temporal_id_and_data->second) { + const auto& mortar_id = received_mortar_data.first; + if (not std::get<3>(received_mortar_data.second).has_value()) { + mortars_to_reconstruct_to.push_back(mortar_id); + } + } + DirectionalIdMap neighbor_reconstructed_evolved_vars = + DgComputeSubcellNeighborPackagedData::apply(*box, + mortars_to_reconstruct_to); + ASSERT(neighbor_reconstructed_evolved_vars.size() == + mortars_to_reconstruct_to.size(), + "Should have reconstructed " + << mortars_to_reconstruct_to.size() << " sides but reconstructed " + << neighbor_reconstructed_evolved_vars.size() << " sides."); + // Now copy over the packaged data _into_ the inbox in order to avoid having + // to make other changes to the DG algorithm (code in + // src/Evolution/DiscontinuousGalerkin) + for (auto& received_mortar_data : received_temporal_id_and_data->second) { + const auto& mortar_id = received_mortar_data.first; + if (not std::get<3>(received_mortar_data.second).has_value()) { + ASSERT(neighbor_reconstructed_evolved_vars.find(mortar_id) != + neighbor_reconstructed_evolved_vars.end(), + "Could not find mortar id " << mortar_id + << " in reconstructed data map."); + std::get<3>(received_mortar_data.second) = + std::move(neighbor_reconstructed_evolved_vars.at(mortar_id)); + } + } +} +} // namespace evolution::dg::subcell diff --git a/src/Evolution/DgSubcell/NeighborTciDecision.cpp b/src/Evolution/DgSubcell/NeighborTciDecision.cpp new file mode 100644 index 000000000000..91e4a4636f42 --- /dev/null +++ b/src/Evolution/DgSubcell/NeighborTciDecision.cpp @@ -0,0 +1,65 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/DgSubcell/NeighborTciDecision.hpp" + +#include +#include +#include +#include + +#include "DataStructures/DataBox/Access.hpp" +#include "DataStructures/DataVector.hpp" +#include "Domain/Structure/Direction.hpp" +#include "Domain/Structure/DirectionalId.hpp" +#include "Domain/Structure/DirectionalIdMap.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Evolution/DgSubcell/Tags/TciStatus.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "Time/TimeStepId.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" +#include "Utilities/GenerateInstantiations.hpp" +#include "Utilities/Gsl.hpp" + +namespace evolution::dg::subcell { +template +void neighbor_tci_decision( + const gsl::not_null box, + const std::pair< + const TimeStepId, + DirectionalIdMap< + Dim, std::tuple, Mesh, std::optional, + std::optional, ::TimeStepId, int>>>& + received_temporal_id_and_data) { + db::mutate>( + [&received_temporal_id_and_data](const auto neighbor_tci_decisions_ptr) { + for (const auto& [directional_element_id, neighbor_data] : + received_temporal_id_and_data.second) { + ASSERT(neighbor_tci_decisions_ptr->contains(directional_element_id), + "The NeighborTciDecisions tag does not contain the neighbor " + << directional_element_id); + neighbor_tci_decisions_ptr->at(directional_element_id) = + std::get<5>(neighbor_data); + } + }, + box); +} + +#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) + +#define INSTANTIATION(r, data) \ + template void neighbor_tci_decision( \ + gsl::not_null box, \ + const std::pair< \ + const TimeStepId, \ + DirectionalIdMap< \ + DIM(data), \ + std::tuple, Mesh, \ + std::optional, std::optional, \ + ::TimeStepId, int>>>& received_temporal_id_and_data); + +GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) + +#undef INSTANTIATION +#undef DIM +} // namespace evolution::dg::subcell diff --git a/src/Evolution/DgSubcell/NeighborTciDecision.hpp b/src/Evolution/DgSubcell/NeighborTciDecision.hpp index 0cb35f65aa9a..b3b7418829b3 100644 --- a/src/Evolution/DgSubcell/NeighborTciDecision.hpp +++ b/src/Evolution/DgSubcell/NeighborTciDecision.hpp @@ -8,16 +8,11 @@ #include #include -#include "DataStructures/DataBox/DataBox.hpp" +#include "DataStructures/DataBox/Access.hpp" #include "DataStructures/DataVector.hpp" -#include "Domain/Structure/Direction.hpp" -#include "Domain/Structure/DirectionalId.hpp" #include "Domain/Structure/DirectionalIdMap.hpp" -#include "Domain/Structure/ElementId.hpp" -#include "Evolution/DgSubcell/Tags/TciStatus.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" #include "Time/TimeStepId.hpp" -#include "Utilities/ErrorHandling/Assert.hpp" #include "Utilities/Gsl.hpp" namespace evolution::dg::subcell { @@ -25,26 +20,13 @@ namespace evolution::dg::subcell { * \brief Copies the neighbors' TCI decisions into * `subcell::Tags::NeighborTciDecisions` */ -template +template void neighbor_tci_decision( - const gsl::not_null*> box, + gsl::not_null box, const std::pair< const TimeStepId, DirectionalIdMap< Dim, std::tuple, Mesh, std::optional, std::optional, ::TimeStepId, int>>>& - received_temporal_id_and_data) { - db::mutate>( - [&received_temporal_id_and_data](const auto neighbor_tci_decisions_ptr) { - for (const auto& [directional_element_id, neighbor_data] : - received_temporal_id_and_data.second) { - ASSERT(neighbor_tci_decisions_ptr->contains(directional_element_id), - "The NeighborTciDecisions tag does not contain the neighbor " - << directional_element_id); - neighbor_tci_decisions_ptr->at(directional_element_id) = - std::get<5>(neighbor_data); - } - }, - box); -} + received_temporal_id_and_data); } // namespace evolution::dg::subcell diff --git a/src/Evolution/DiscontinuousGalerkin/Actions/ApplyBoundaryCorrections.hpp b/src/Evolution/DiscontinuousGalerkin/Actions/ApplyBoundaryCorrections.hpp index a60f758451f4..a57e9deb6e64 100644 --- a/src/Evolution/DiscontinuousGalerkin/Actions/ApplyBoundaryCorrections.hpp +++ b/src/Evolution/DiscontinuousGalerkin/Actions/ApplyBoundaryCorrections.hpp @@ -12,7 +12,9 @@ #include #include +#include "DataStructures/DataBox/AsAccess.hpp" #include "DataStructures/DataBox/DataBox.hpp" +#include "DataStructures/DataBox/PrefixHelpers.hpp" #include "DataStructures/DataBox/Prefixes.hpp" #include "DataStructures/Tensor/EagerMath/Magnitude.hpp" #include "Domain/FaceNormal.hpp" @@ -67,10 +69,9 @@ struct TimeStepper; namespace evolution::dg::subcell { // We use a forward declaration instead of including a header file to avoid // coupling to the DG-subcell libraries for executables that don't use subcell. -template +template void neighbor_reconstructed_face_solution( - gsl::not_null*> box, + gsl::not_null box, gsl::not_null, std::optional, ::TimeStepId, int>>>*> received_temporal_id_and_data); -template +template void neighbor_tci_decision( - gsl::not_null*> box, + gsl::not_null box, const std::pair< const TimeStepId, DirectionalIdMap< @@ -92,6 +93,19 @@ void neighbor_tci_decision( /// \endcond namespace evolution::dg { +namespace detail { +template +struct get_dg_boundary_terms { + using type = typename BoundaryCorrectionClass::dg_boundary_terms_volume_tags; +}; + +template >> +struct TemporaryReference { + using tag = Tag; + using type = const Type&; +}; +} // namespace detail + /// Receive boundary data for global time-stepping. Returns true if /// all necessary data has been received. template @@ -132,9 +146,9 @@ bool receive_boundary_data_global_time_stepping( evolution::dg::subcell::neighbor_reconstructed_face_solution< volume_dim, typename Metavariables::SubcellOptions:: DgComputeSubcellNeighborPackagedData>( - box, make_not_null(&*received_temporal_id_and_data)); + &db::as_access(*box), make_not_null(&*received_temporal_id_and_data)); evolution::dg::subcell::neighbor_tci_decision( - box, *received_temporal_id_and_data); + make_not_null(&db::as_access(*box)), *received_temporal_id_and_data); } db::mutate, @@ -340,6 +354,10 @@ struct ApplyBoundaryCorrections { using variables_tag = typename system::variables_tag; using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>; using DtVariables = typename dt_variables_tag::type; + using volume_tags_for_dg_boundary_terms = + tmpl::remove_duplicates>>>; using TimeStepperType = tmpl::conditional_t; @@ -357,19 +375,22 @@ struct ApplyBoundaryCorrections { using return_tags = tmpl::conditional_t, tmpl::list>; - using argument_tags = tmpl::flatten>, - domain::Tags::Mesh, Tags::MortarMesh, - Tags::MortarSize, ::dg::Tags::Formulation, - evolution::dg::Tags::NormalCovectorAndMagnitude, - ::Tags::TimeStepper, - evolution::Tags::BoundaryCorrection, - tmpl::conditional_t, - tmpl::conditional_t, - domain::Tags::DetInvJacobian>>>; + using argument_tags = tmpl::append< + tmpl::flatten>, + domain::Tags::Mesh, Tags::MortarMesh, + Tags::MortarSize, ::dg::Tags::Formulation, + evolution::dg::Tags::NormalCovectorAndMagnitude, + ::Tags::TimeStepper, + evolution::Tags::BoundaryCorrection, + tmpl::conditional_t, + tmpl::conditional_t, + domain::Tags::DetInvJacobian< + Frame::ElementLogical, Frame::Inertial>>>>, + volume_tags_for_dg_boundary_terms>; // full step + template static void apply( const gsl::not_null vars_to_update, const gsl::not_null mortar_data, @@ -385,15 +406,40 @@ struct ApplyBoundaryCorrections { const TimeStepperType& time_stepper, const typename system::boundary_correction_base& boundary_correction, const TimeDelta& time_step, - const Scalar& gts_det_inv_jacobian = {}) { + const Scalar& gts_det_inv_jacobian, + const VolumeArgs&... volume_args) { apply_impl(vars_to_update, mortar_data, volume_mesh, mortar_meshes, mortar_sizes, dg_formulation, face_normal_covector_and_magnitude, time_stepper, boundary_correction, time_step, std::numeric_limits::signaling_NaN(), - gts_det_inv_jacobian); + gts_det_inv_jacobian, volume_args...); + } + + template + static void apply( + const gsl::not_null vars_to_update, + const gsl::not_null mortar_data, + const Mesh& volume_mesh, + const typename Tags::MortarMesh::type& mortar_meshes, + const typename Tags::MortarSize::type& mortar_sizes, + const ::dg::Formulation dg_formulation, + const DirectionMap< + volume_dim, std::optional>>>>& + face_normal_covector_and_magnitude, + const TimeStepperType& time_stepper, + const typename system::boundary_correction_base& boundary_correction, + const TimeDelta& time_step, const VolumeArgs&... volume_args) { + apply_impl(vars_to_update, mortar_data, volume_mesh, mortar_meshes, + mortar_sizes, dg_formulation, face_normal_covector_and_magnitude, + time_stepper, boundary_correction, time_step, + std::numeric_limits::signaling_NaN(), {}, + volume_args...); } // dense output (LTS only) + template static void apply( const gsl::not_null vars_to_update, const MortarDataType& mortar_data, const Mesh& volume_mesh, @@ -407,11 +453,11 @@ struct ApplyBoundaryCorrections { face_normal_covector_and_magnitude, const LtsTimeStepper& time_stepper, const typename system::boundary_correction_base& boundary_correction, - const double dense_output_time) { + const double dense_output_time, const VolumeArgs&... volume_args) { apply_impl(vars_to_update, &mortar_data, volume_mesh, mortar_meshes, mortar_sizes, dg_formulation, face_normal_covector_and_magnitude, time_stepper, boundary_correction, TimeDelta{}, - dense_output_time, {}); + dense_output_time, {}, volume_args...); } template static void apply_impl( const gsl::not_null vars_to_update, const gsl::not_null mortar_data, @@ -448,7 +495,12 @@ struct ApplyBoundaryCorrections { const TimeStepperType& time_stepper, const typename system::boundary_correction_base& boundary_correction, const TimeDelta& time_step, const double dense_output_time, - const Scalar& gts_det_inv_jacobian) { + const Scalar& gts_det_inv_jacobian, + const VolumeArgs&... volume_args) { + tuples::tagged_tuple_from_typelist> + volume_args_tuple{volume_args...}; + // Set up helper lambda that will compute and lift the boundary corrections ASSERT( volume_mesh.quadrature() == @@ -481,13 +533,14 @@ struct ApplyBoundaryCorrections { [&dense_output_time, &dg_formulation, &face_normal_covector_and_magnitude, &mortar_data, &mortar_meshes, &mortar_sizes, &time_step, &time_stepper, using_gauss_lobatto_points, - &vars_to_update, &volume_det_jacobian, &volume_det_inv_jacobian, + &vars_to_update, &volume_args_tuple, &volume_det_jacobian, + &volume_det_inv_jacobian, &volume_mesh](auto* typed_boundary_correction) { + using BcType = std::decay_t; // Compute internal boundary quantities on the mortar for sides of // the element that have neighbors, i.e. they are not an external // side. - using mortar_tags_list = typename std::decay_t::dg_package_field_tags; + using mortar_tags_list = typename BcType::dg_package_field_tags; // Variables for reusing allocations. The actual values are // not reused. @@ -526,8 +579,9 @@ struct ApplyBoundaryCorrections { &face_mesh, &face_normal_covector_and_magnitude, &local_data_on_mortar, &mortar_id, &mortar_meshes, &mortar_sizes, &neighbor_data_on_mortar, - using_gauss_lobatto_points, &volume_det_jacobian, - &volume_det_inv_jacobian, &volume_dt_correction, &volume_mesh]( + using_gauss_lobatto_points, &volume_args_tuple, + &volume_det_jacobian, &volume_det_inv_jacobian, + &volume_dt_correction, &volume_mesh]( const MortarData& local_mortar_data, const MortarData& neighbor_mortar_data) -> DtVariables { @@ -571,7 +625,8 @@ struct ApplyBoundaryCorrections { call_boundary_correction( make_not_null(&dt_boundary_correction_on_mortar), local_data_on_mortar, neighbor_data_on_mortar, - *typed_boundary_correction, dg_formulation); + *typed_boundary_correction, dg_formulation, volume_args_tuple, + typename BcType::dg_boundary_terms_volume_tags{}); const std::array& mortar_size = mortar_sizes.at(mortar_id); @@ -739,19 +794,25 @@ struct ApplyBoundaryCorrections { } template + typename BoundaryCorrection, typename... AllVolumeArgs, + typename... VolumeTagsForCorrection> static void call_boundary_correction( const gsl::not_null>*> boundary_corrections_on_mortar, const Variables>& local_boundary_data, const Variables>& neighbor_boundary_data, const BoundaryCorrection& boundary_correction, - const ::dg::Formulation dg_formulation) { + const ::dg::Formulation dg_formulation, + const tuples::TaggedTuple...>& + volume_args_tuple, + tmpl::list /*meta*/) { boundary_correction.dg_boundary_terms( make_not_null( &get(*boundary_corrections_on_mortar))..., get(local_boundary_data)..., get(neighbor_boundary_data)..., - dg_formulation); + dg_formulation, + tuples::get>( + volume_args_tuple)...); } }; diff --git a/src/Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp b/src/Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp index fef47e6ece10..32c960420183 100644 --- a/src/Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp +++ b/src/Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp @@ -66,7 +66,7 @@ double dg_package_data( template + typename... VolumeTags> double dg_package_data( const gsl::not_null>*> packaged_data, @@ -75,7 +75,7 @@ double dg_package_data( const tnsr::i& unit_normal_covector, const std::optional>& mesh_velocity, - const db::DataBox& box, tmpl::list /*meta*/, + const db::Access& box, tmpl::list /*meta*/, tmpl::list /*meta*/) { return dg_package_data( packaged_data, boundary_correction, projected_fields, diff --git a/src/Evolution/Systems/Burgers/BoundaryCorrections/Hll.hpp b/src/Evolution/Systems/Burgers/BoundaryCorrections/Hll.hpp index 89b07bd78212..f4ae10d9480e 100644 --- a/src/Evolution/Systems/Burgers/BoundaryCorrections/Hll.hpp +++ b/src/Evolution/Systems/Burgers/BoundaryCorrections/Hll.hpp @@ -94,6 +94,7 @@ class Hll final : public BoundaryCorrection { tmpl::list, CharSpeed>; using dg_package_data_temporary_tags = tmpl::list<>; using dg_package_data_volume_tags = tmpl::list<>; + using dg_boundary_terms_volume_tags = tmpl::list<>; static double dg_package_data( gsl::not_null*> packaged_u, diff --git a/src/Evolution/Systems/Burgers/BoundaryCorrections/Rusanov.hpp b/src/Evolution/Systems/Burgers/BoundaryCorrections/Rusanov.hpp index 85ef8bb64d9a..286a6dfe6b55 100644 --- a/src/Evolution/Systems/Burgers/BoundaryCorrections/Rusanov.hpp +++ b/src/Evolution/Systems/Burgers/BoundaryCorrections/Rusanov.hpp @@ -80,6 +80,7 @@ class Rusanov final : public BoundaryCorrection { tmpl::list, AbsCharSpeed>; using dg_package_data_temporary_tags = tmpl::list<>; using dg_package_data_volume_tags = tmpl::list<>; + using dg_boundary_terms_volume_tags = tmpl::list<>; static double dg_package_data( gsl::not_null*> packaged_u, diff --git a/src/Evolution/Systems/Burgers/Subcell/CMakeLists.txt b/src/Evolution/Systems/Burgers/Subcell/CMakeLists.txt index 3cc556797e65..ff42c83de18f 100644 --- a/src/Evolution/Systems/Burgers/Subcell/CMakeLists.txt +++ b/src/Evolution/Systems/Burgers/Subcell/CMakeLists.txt @@ -5,6 +5,7 @@ spectre_target_sources( ${LIBRARY} PRIVATE GhostData.cpp + NeighborPackagedData.cpp SetInitialRdmpData.cpp TciOnDgGrid.cpp TciOnFdGrid.cpp diff --git a/src/Evolution/Systems/Burgers/Subcell/NeighborPackagedData.cpp b/src/Evolution/Systems/Burgers/Subcell/NeighborPackagedData.cpp new file mode 100644 index 000000000000..da60e36a96c4 --- /dev/null +++ b/src/Evolution/Systems/Burgers/Subcell/NeighborPackagedData.cpp @@ -0,0 +1,156 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/Systems/Burgers/Subcell/NeighborPackagedData.hpp" + +#include +#include +#include + +#include "DataStructures/DataBox/Access.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "Domain/Structure/Direction.hpp" +#include "Domain/Structure/DirectionalId.hpp" +#include "Domain/Structure/DirectionalIdMap.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Domain/TagsTimeDependent.hpp" +#include "Evolution/BoundaryCorrectionTags.hpp" +#include "Evolution/DgSubcell/Mesh.hpp" +#include "Evolution/DgSubcell/NeighborReconstructedFaceSolution.tpp" +#include "Evolution/DgSubcell/Projection.hpp" +#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" +#include "Evolution/DgSubcell/Tags/Mesh.hpp" +#include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp" +#include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp" +#include "Evolution/Systems/Burgers/BoundaryCorrections/BoundaryCorrection.hpp" +#include "Evolution/Systems/Burgers/BoundaryCorrections/Factory.hpp" +#include "Evolution/Systems/Burgers/FiniteDifference/Factory.hpp" +#include "Evolution/Systems/Burgers/FiniteDifference/Reconstructor.hpp" +#include "Evolution/Systems/Burgers/FiniteDifference/Tags.hpp" +#include "Evolution/Systems/Burgers/Fluxes.hpp" +#include "Evolution/Systems/Burgers/Subcell/ComputeFluxes.hpp" +#include "Evolution/Systems/Burgers/System.hpp" +#include "Evolution/Systems/Burgers/Tags.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "Utilities/CallWithDynamicType.hpp" +#include "Utilities/GenerateInstantiations.hpp" +#include "Utilities/TMPL.hpp" + +namespace Burgers::subcell { +DirectionalIdMap<1, DataVector> NeighborPackagedData::apply( + const db::Access& box, + const std::vector>& mortars_to_reconstruct_to) { + // The object to return + DirectionalIdMap<1, DataVector> neighbor_package_data{}; + if (mortars_to_reconstruct_to.empty()) { + return neighbor_package_data; + } + + using evolved_vars_tags = typename System::variables_tag::tags_list; + using fluxes_tags = typename Fluxes::return_tags; + + // subcell currently does not support moving mesh + ASSERT(not db::get>(box).has_value(), + "Haven't yet added support for moving mesh to DG-subcell. This " + "should be easy to generalize, but we will want to consider " + "storing the mesh velocity on the faces instead of " + "re-slicing/projecting."); + + // Project volume variables from DG to subcell mesh + const Mesh<1>& dg_mesh = db::get>(box); + const Mesh<1>& subcell_mesh = + db::get>(box); + const auto volume_vars_subcell = evolution::dg::subcell::fd::project( + db::get(box), dg_mesh, + subcell_mesh.extents()); + + const auto& ghost_subcell_data = + db::get>(box); + + const Burgers::fd::Reconstructor& recons = + db::get(box); + + const auto& boundary_correction = + db::get>(box); + using derived_boundary_corrections = + typename std::decay_t::creatable_classes; + + // perform reconstruction + tmpl::for_each([&](auto derived_correction_v) { + using derived_correction = tmpl::type_from; + if (typeid(boundary_correction) == typeid(derived_correction)) { + using dg_package_field_tags = + typename derived_correction::dg_package_field_tags; + using dg_package_data_argument_tags = + tmpl::append; + + const auto& element = db::get>(box); + + // Variables to store packaged data + Variables packaged_data{}; + // Variables to be reconstructed on the shared interfaces + Variables vars_on_face{}; + + for (const auto& mortar_id : mortars_to_reconstruct_to) { + // Note : 1D mortar has only one face point + const size_t num_face_pts = 1; + vars_on_face.initialize(num_face_pts); + + // Reconstruct field variables on faces + call_with_dynamic_type< + void, typename Burgers::fd::Reconstructor::creatable_classes>( + &recons, + [&element, &mortar_id, &ghost_subcell_data, &subcell_mesh, + &vars_on_face, &volume_vars_subcell](const auto& reconstructor) { + reconstructor->reconstruct_fd_neighbor( + make_not_null(&vars_on_face), volume_vars_subcell, element, + ghost_subcell_data, subcell_mesh, mortar_id.direction); + }); + + // Compute fluxes + Burgers::subcell::compute_fluxes(make_not_null(&vars_on_face)); + + tnsr::i normal_covector = + get>( + *db::get>( + box) + .at(mortar_id.direction)); + for (auto& t : normal_covector) { + t *= -1.0; + } + + // Compute the packaged data + packaged_data.initialize(num_face_pts); + evolution::dg::Actions::detail::dg_package_data( + make_not_null(&packaged_data), + dynamic_cast(boundary_correction), + vars_on_face, normal_covector, {std::nullopt}, box, + typename derived_correction::dg_package_data_volume_tags{}, + dg_package_data_argument_tags{}); + + // Make a view so we can use iterators with std::copy + DataVector packaged_data_view{packaged_data.data(), + packaged_data.size()}; + neighbor_package_data[mortar_id] = DataVector{packaged_data.size()}; + std::copy(packaged_data_view.begin(), packaged_data_view.end(), + neighbor_package_data[mortar_id].begin()); + } + } + }); + + return neighbor_package_data; +} +} // namespace Burgers::subcell + +template void evolution::dg::subcell::neighbor_reconstructed_face_solution< + 1, Burgers::subcell::NeighborPackagedData>( + gsl::not_null box, + gsl::not_null, Mesh<0>, std::optional, + std::optional, ::TimeStepId, int>>>*> + received_temporal_id_and_data); diff --git a/src/Evolution/Systems/Burgers/Subcell/NeighborPackagedData.hpp b/src/Evolution/Systems/Burgers/Subcell/NeighborPackagedData.hpp index 3987858da068..94e70b24e62f 100644 --- a/src/Evolution/Systems/Burgers/Subcell/NeighborPackagedData.hpp +++ b/src/Evolution/Systems/Burgers/Subcell/NeighborPackagedData.hpp @@ -4,40 +4,18 @@ #pragma once #include -#include -#include -#include - -#include "DataStructures/DataBox/DataBox.hpp" -#include "DataStructures/DataVector.hpp" -#include "DataStructures/Index.hpp" -#include "DataStructures/Tensor/Slice.hpp" -#include "DataStructures/Tensor/Tensor.hpp" -#include "DataStructures/Variables.hpp" -#include "Domain/Structure/Direction.hpp" -#include "Domain/Structure/DirectionalId.hpp" -#include "Domain/Structure/DirectionalIdMap.hpp" -#include "Domain/Structure/ElementId.hpp" -#include "Domain/TagsTimeDependent.hpp" -#include "Evolution/BoundaryCorrectionTags.hpp" -#include "Evolution/DgSubcell/Mesh.hpp" -#include "Evolution/DgSubcell/Projection.hpp" -#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" -#include "Evolution/DgSubcell/Tags/Mesh.hpp" -#include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp" -#include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp" -#include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp" -#include "Evolution/Systems/Burgers/BoundaryCorrections/BoundaryCorrection.hpp" -#include "Evolution/Systems/Burgers/BoundaryCorrections/Factory.hpp" -#include "Evolution/Systems/Burgers/FiniteDifference/Reconstructor.hpp" -#include "Evolution/Systems/Burgers/FiniteDifference/Tags.hpp" -#include "Evolution/Systems/Burgers/Fluxes.hpp" -#include "Evolution/Systems/Burgers/Subcell/ComputeFluxes.hpp" -#include "Evolution/Systems/Burgers/System.hpp" -#include "Evolution/Systems/Burgers/Tags.hpp" -#include "NumericalAlgorithms/Spectral/Mesh.hpp" -#include "Utilities/CallWithDynamicType.hpp" -#include "Utilities/TMPL.hpp" +#include + +/// \cond +namespace db { +class Access; +} // namespace db +class DataVector; +template +struct DirectionalId; +template +class DirectionalIdMap; +/// \endcond namespace Burgers::subcell { /*! @@ -55,112 +33,8 @@ namespace Burgers::subcell { * and if one element is doing DG, then we aren't at a shock. */ struct NeighborPackagedData { - template static DirectionalIdMap<1, DataVector> apply( - const db::DataBox& box, - const std::vector>& mortars_to_reconstruct_to) { - // The object to return - DirectionalIdMap<1, DataVector> neighbor_package_data{}; - if (mortars_to_reconstruct_to.empty()) { - return neighbor_package_data; - } - - using evolved_vars_tags = typename System::variables_tag::tags_list; - using fluxes_tags = typename Fluxes::return_tags; - - // subcell currently does not support moving mesh - ASSERT(not db::get>(box).has_value(), - "Haven't yet added support for moving mesh to DG-subcell. This " - "should be easy to generalize, but we will want to consider " - "storing the mesh velocity on the faces instead of " - "re-slicing/projecting."); - - // Project volume variables from DG to subcell mesh - const Mesh<1>& dg_mesh = db::get>(box); - const Mesh<1>& subcell_mesh = - db::get>(box); - const auto volume_vars_subcell = evolution::dg::subcell::fd::project( - db::get(box), dg_mesh, - subcell_mesh.extents()); - - const auto& ghost_subcell_data = - db::get>( - box); - - const Burgers::fd::Reconstructor& recons = - db::get(box); - - const auto& boundary_correction = - db::get>(box); - using derived_boundary_corrections = - typename std::decay_t::creatable_classes; - - // perform reconstruction - tmpl::for_each([&](auto - derived_correction_v) { - using derived_correction = - tmpl::type_from; - if (typeid(boundary_correction) == typeid(derived_correction)) { - using dg_package_field_tags = - typename derived_correction::dg_package_field_tags; - using dg_package_data_argument_tags = - tmpl::append; - - const auto& element = db::get>(box); - - // Variables to store packaged data - Variables packaged_data{}; - // Variables to be reconstructed on the shared interfaces - Variables vars_on_face{}; - - for (const auto& mortar_id : mortars_to_reconstruct_to) { - // Note : 1D mortar has only one face point - const size_t num_face_pts = 1; - vars_on_face.initialize(num_face_pts); - - // Reconstruct field variables on faces - call_with_dynamic_type< - void, typename Burgers::fd::Reconstructor::creatable_classes>( - &recons, - [&element, &mortar_id, &ghost_subcell_data, &subcell_mesh, - &vars_on_face, &volume_vars_subcell](const auto& reconstructor) { - reconstructor->reconstruct_fd_neighbor( - make_not_null(&vars_on_face), volume_vars_subcell, element, - ghost_subcell_data, subcell_mesh, mortar_id.direction); - }); - - // Compute fluxes - Burgers::subcell::compute_fluxes(make_not_null(&vars_on_face)); - - tnsr::i normal_covector = - get>( - *db::get>( - box) - .at(mortar_id.direction)); - for (auto& t : normal_covector) { - t *= -1.0; - } - - // Compute the packaged data - packaged_data.initialize(num_face_pts); - evolution::dg::Actions::detail::dg_package_data( - make_not_null(&packaged_data), - dynamic_cast(boundary_correction), - vars_on_face, normal_covector, {std::nullopt}, box, - typename derived_correction::dg_package_data_volume_tags{}, - dg_package_data_argument_tags{}); - - // Make a view so we can use iterators with std::copy - DataVector packaged_data_view{packaged_data.data(), - packaged_data.size()}; - neighbor_package_data[mortar_id] = DataVector{packaged_data.size()}; - std::copy(packaged_data_view.begin(), packaged_data_view.end(), - neighbor_package_data[mortar_id].begin()); - } - } - }); - - return neighbor_package_data; - } + const db::Access& box, + const std::vector>& mortars_to_reconstruct_to); }; } // namespace Burgers::subcell diff --git a/src/Evolution/Systems/CurvedScalarWave/BoundaryCorrections/UpwindPenalty.hpp b/src/Evolution/Systems/CurvedScalarWave/BoundaryCorrections/UpwindPenalty.hpp index 1b10e79810e1..387e1152d330 100644 --- a/src/Evolution/Systems/CurvedScalarWave/BoundaryCorrections/UpwindPenalty.hpp +++ b/src/Evolution/Systems/CurvedScalarWave/BoundaryCorrections/UpwindPenalty.hpp @@ -120,6 +120,7 @@ class UpwindPenalty final : public BoundaryCorrection { tmpl::list, gr::Tags::Shift, Tags::ConstraintGamma1, Tags::ConstraintGamma2>; using dg_package_data_volume_tags = tmpl::list<>; + using dg_boundary_terms_volume_tags = tmpl::list<>; double dg_package_data( gsl::not_null*> packaged_v_psi, diff --git a/src/Evolution/Systems/ForceFree/BoundaryCorrections/Rusanov.hpp b/src/Evolution/Systems/ForceFree/BoundaryCorrections/Rusanov.hpp index a191dc8b3edf..369200b78806 100644 --- a/src/Evolution/Systems/ForceFree/BoundaryCorrections/Rusanov.hpp +++ b/src/Evolution/Systems/ForceFree/BoundaryCorrections/Rusanov.hpp @@ -101,6 +101,7 @@ class Rusanov final : public BoundaryCorrection { tmpl::list, gr::Tags::Shift>; using dg_package_data_primitive_tags = tmpl::list<>; using dg_package_data_volume_tags = tmpl::list<>; + using dg_boundary_terms_volume_tags = tmpl::list<>; static double dg_package_data( gsl::not_null*> packaged_tilde_e, diff --git a/src/Evolution/Systems/ForceFree/Subcell/CMakeLists.txt b/src/Evolution/Systems/ForceFree/Subcell/CMakeLists.txt index 43c3419b82ab..a0d521798108 100644 --- a/src/Evolution/Systems/ForceFree/Subcell/CMakeLists.txt +++ b/src/Evolution/Systems/ForceFree/Subcell/CMakeLists.txt @@ -7,6 +7,7 @@ spectre_target_sources( PRIVATE ComputeFluxes.hpp GhostData.cpp + NeighborPackagedData.cpp SetInitialRdmpData.cpp TciOnDgGrid.cpp TciOnFdGrid.cpp diff --git a/src/Evolution/Systems/ForceFree/Subcell/NeighborPackagedData.cpp b/src/Evolution/Systems/ForceFree/Subcell/NeighborPackagedData.cpp new file mode 100644 index 000000000000..1cd81f5be3d3 --- /dev/null +++ b/src/Evolution/Systems/ForceFree/Subcell/NeighborPackagedData.cpp @@ -0,0 +1,272 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/Systems/ForceFree/Subcell/NeighborPackagedData.hpp" + +#include +#include +#include +#include +#include +#include + +#include "DataStructures/DataBox/Access.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Index.hpp" +#include "DataStructures/Tensor/Slice.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "Domain/Structure/Direction.hpp" +#include "Domain/Structure/DirectionalId.hpp" +#include "Domain/Structure/DirectionalIdMap.hpp" +#include "Domain/Structure/Element.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Domain/Structure/Side.hpp" +#include "Domain/Tags.hpp" +#include "Domain/TagsTimeDependent.hpp" +#include "Evolution/BoundaryCorrectionTags.hpp" +#include "Evolution/DgSubcell/NeighborReconstructedFaceSolution.tpp" +#include "Evolution/DgSubcell/Projection.hpp" +#include "Evolution/DgSubcell/Reconstruction.hpp" +#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" +#include "Evolution/DgSubcell/Tags/Mesh.hpp" +#include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp" +#include "Evolution/DgSubcell/Tags/SubcellOptions.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp" +#include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp" +#include "Evolution/Systems/ForceFree/BoundaryCorrections/Factory.hpp" +#include "Evolution/Systems/ForceFree/FiniteDifference/Factory.hpp" +#include "Evolution/Systems/ForceFree/FiniteDifference/ReconstructWork.hpp" +#include "Evolution/Systems/ForceFree/FiniteDifference/Reconstructor.hpp" +#include "Evolution/Systems/ForceFree/FiniteDifference/Tags.hpp" +#include "Evolution/Systems/ForceFree/Fluxes.hpp" +#include "Evolution/Systems/ForceFree/Subcell/ComputeFluxes.hpp" +#include "Evolution/Systems/ForceFree/System.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "Utilities/CallWithDynamicType.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/TMPL.hpp" + +namespace ForceFree::subcell { +DirectionalIdMap<3, DataVector> NeighborPackagedData::apply( + const db::Access& box, + const std::vector>& mortars_to_reconstruct_to) { + using evolved_vars_tags = typename System::variables_tag::tags_list; + using fluxes_tags = typename Fluxes::return_tags; + + // Subcell with moving mesh is not supported in the ForceFree system yet. + ASSERT(not db::get>(box).has_value(), + "Haven't yet added support for moving mesh to DG-subcell. This " + "should be easy to generalize, but we will want to consider " + "storing the mesh velocity on the faces instead of " + "re-slicing/projecting."); + + DirectionalIdMap<3, DataVector> neighbor_package_data{}; + if (mortars_to_reconstruct_to.empty()) { + return neighbor_package_data; + } + + // Project volume variables from DG to subcell mesh + const Mesh<3>& dg_mesh = db::get>(box); + const Mesh<3>& subcell_mesh = + db::get>(box); + const auto volume_vars_subcell = evolution::dg::subcell::fd::project( + db::get(box), dg_mesh, + subcell_mesh.extents()); + + // In addition to evolved vars, compute TildeJ on the subcell grid so that + // it could be sent to neighbor elements doing FD + DataVector buffer{ + Variables>::number_of_independent_components * + dg_mesh.number_of_grid_points()}; + const auto volume_tildej_subcell = [&box, &buffer, &dg_mesh, + &subcell_mesh]() { + Variables> var{buffer.data(), buffer.size()}; + for (size_t i = 0; i < 3; ++i) { + get(var).get(i) = db::get(box).get(i); + } + return get(evolution::dg::subcell::fd::project( + var, dg_mesh, subcell_mesh.extents())); + }(); + + const auto& ghost_subcell_data = + db::get>(box); + + const ForceFree::fd::Reconstructor& recons = + db::get(box); + + const auto& boundary_correction = + db::get>(box); + using derived_boundary_corrections = + typename std::decay_t::creatable_classes; + + const auto& subcell_options = + db::get>(box); + + tmpl::for_each([&box, &boundary_correction, + &dg_mesh, + &mortars_to_reconstruct_to, + &neighbor_package_data, + &ghost_subcell_data, &recons, + &subcell_mesh, &subcell_options, + &volume_vars_subcell, + &volume_tildej_subcell]( + auto derived_correction_v) { + using DerivedCorrection = tmpl::type_from; + + using dg_package_field_tags = + typename DerivedCorrection::dg_package_field_tags; + using dg_package_data_temporary_tags = + typename DerivedCorrection::dg_package_data_temporary_tags; + + if (typeid(boundary_correction) == typeid(DerivedCorrection)) { + Variables packaged_data{0}; + + // Comprehensive tags list for computing fluxes and dg_package_data + using dg_package_data_argument_tags = tmpl::append< + evolved_vars_tags, fluxes_tags, + tmpl::remove_duplicates, + gr::Tags::SqrtDetSpatialMetric, + gr::Tags::InverseSpatialMetric, + evolution::dg::Actions::detail::NormalVector<3>>>>; + + // We allocate a Variables object with extra buffer for storing TildeJ. + // TildeJ is not an argument of dg_package_data() function, but it needs + // to be reconstructed for computing fluxes on faces. + Variables, + dg_package_data_argument_tags>> + vars_on_face_with_tildej_buffer{0}; + + for (const auto& mortar_id : mortars_to_reconstruct_to) { + const Direction<3>& direction = mortar_id.direction; + + // Switch to face-centered points + Index<3> extents = subcell_mesh.extents(); + ++extents[direction.dimension()]; + const size_t num_face_pts = + subcell_mesh.extents().slice_away(direction.dimension()).product(); + vars_on_face_with_tildej_buffer.initialize(num_face_pts); + + // create a `view` of the face vars excluding the TildeJ tag, since + // boundary correction takes the flux F^i(TildeQ) as argument but not + // TildeJ itself. + auto vars_on_face = + vars_on_face_with_tildej_buffer + .template reference_subset(); + auto reconstructed_vars_on_face = + vars_on_face_with_tildej_buffer.template reference_subset< + ForceFree::fd::tags_list_for_reconstruction>(); + + // Copy spacetime vars over from volume. + using spacetime_vars_to_copy = tmpl::list< + gr::Tags::Lapse, + gr::Tags::Shift, + gr::Tags::SpatialMetric, + gr::Tags::SqrtDetSpatialMetric, + gr::Tags::InverseSpatialMetric>; + tmpl::for_each( + [&direction, &extents, &vars_on_face, + &spacetime_vars_on_faces = + db::get>(box)]( + auto tag_v) { + using tag = tmpl::type_from; + data_on_slice(make_not_null(&get(vars_on_face)), + get(gsl::at(spacetime_vars_on_faces, + direction.dimension())), + extents, direction.dimension(), + direction.side() == Side::Lower + ? 0 + : extents[direction.dimension()] - 1); + }); + + // Perform FD reconstruction of variables on cell interfaces. Note + // that we are using the vars object with TildeJ buffer so that + // TildeJ can be reconstructed in this step. + call_with_dynamic_type< + void, typename ForceFree::fd::Reconstructor::creatable_classes>( + &recons, [&element = db::get>(box), + &mortar_id, &ghost_subcell_data, &subcell_mesh, + &reconstructed_vars_on_face, &volume_vars_subcell, + &volume_tildej_subcell](const auto& reconstructor) { + reconstructor->reconstruct_fd_neighbor( + make_not_null(&reconstructed_vars_on_face), + volume_vars_subcell, volume_tildej_subcell, element, + ghost_subcell_data, subcell_mesh, mortar_id.direction); + }); + ForceFree::subcell::compute_fluxes( + make_not_null(&vars_on_face_with_tildej_buffer)); + + // Note: since the spacetime isn't dynamical we don't need to + // worry about different normal vectors on the different sides + // of the element. If the spacetime is dynamical, then the normal + // covector needs to be sent, or at least the inverse spatial metric + // from the neighbor so that the normal covector can be computed + // correctly. + tnsr::i normal_covector = + get>( + *db::get>( + box) + .at(mortar_id.direction)); + for (auto& t : normal_covector) { + t *= -1.0; + } + + // Note: Only need to do the projection in 2d and 3d, but GRFFE is + // always 3d currently. + const auto dg_normal_covector = normal_covector; + for (size_t i = 0; i < 3; ++i) { + normal_covector.get(i) = evolution::dg::subcell::fd::project( + dg_normal_covector.get(i), + dg_mesh.slice_away(mortar_id.direction.dimension()), + subcell_mesh.extents().slice_away( + mortar_id.direction.dimension())); + } + + // Compute the packaged data + packaged_data.initialize(num_face_pts); + using dg_package_data_projected_tags = + tmpl::append; + evolution::dg::Actions::detail::dg_package_data( + make_not_null(&packaged_data), + dynamic_cast(boundary_correction), + vars_on_face, normal_covector, {std::nullopt}, box, + typename DerivedCorrection::dg_package_data_volume_tags{}, + dg_package_data_projected_tags{}); + + // Reconstruct the DG solution. Really we should be solving the + // boundary correction and then reconstructing, but away from a shock + // this doesn't matter. + auto dg_packaged_data = evolution::dg::subcell::fd::reconstruct( + packaged_data, dg_mesh.slice_away(mortar_id.direction.dimension()), + subcell_mesh.extents().slice_away(mortar_id.direction.dimension()), + subcell_options.reconstruction_method()); + + // Make a view so we can use iterators with std::copy + DataVector dg_packaged_data_view{dg_packaged_data.data(), + dg_packaged_data.size()}; + neighbor_package_data[mortar_id] = DataVector{dg_packaged_data.size()}; + std::copy(dg_packaged_data_view.begin(), dg_packaged_data_view.end(), + neighbor_package_data[mortar_id].begin()); + } + } + }); + + return neighbor_package_data; +} +} // namespace ForceFree::subcell + +template void evolution::dg::subcell::neighbor_reconstructed_face_solution< + 3, ForceFree::subcell::NeighborPackagedData>( + gsl::not_null box, + gsl::not_null, Mesh<2>, std::optional, + std::optional, ::TimeStepId, int>>>*> + received_temporal_id_and_data); diff --git a/src/Evolution/Systems/ForceFree/Subcell/NeighborPackagedData.hpp b/src/Evolution/Systems/ForceFree/Subcell/NeighborPackagedData.hpp index 37ee2e191857..259eb55c78f4 100644 --- a/src/Evolution/Systems/ForceFree/Subcell/NeighborPackagedData.hpp +++ b/src/Evolution/Systems/ForceFree/Subcell/NeighborPackagedData.hpp @@ -3,49 +3,19 @@ #pragma once -#include #include -#include -#include -#include #include -#include "DataStructures/DataBox/DataBox.hpp" -#include "DataStructures/DataVector.hpp" -#include "DataStructures/Index.hpp" -#include "DataStructures/Tensor/Slice.hpp" -#include "DataStructures/Tensor/Tensor.hpp" -#include "DataStructures/Variables.hpp" -#include "Domain/Structure/Direction.hpp" -#include "Domain/Structure/DirectionalId.hpp" -#include "Domain/Structure/DirectionalIdMap.hpp" -#include "Domain/Structure/Element.hpp" -#include "Domain/Structure/ElementId.hpp" -#include "Domain/Structure/Side.hpp" -#include "Domain/Tags.hpp" -#include "Domain/TagsTimeDependent.hpp" -#include "Evolution/BoundaryCorrectionTags.hpp" -#include "Evolution/DgSubcell/Projection.hpp" -#include "Evolution/DgSubcell/Reconstruction.hpp" -#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" -#include "Evolution/DgSubcell/Tags/Mesh.hpp" -#include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp" -#include "Evolution/DgSubcell/Tags/SubcellOptions.hpp" -#include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp" -#include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp" -#include "Evolution/Systems/ForceFree/FiniteDifference/Factory.hpp" -#include "Evolution/Systems/ForceFree/FiniteDifference/ReconstructWork.hpp" -#include "Evolution/Systems/ForceFree/FiniteDifference/Reconstructor.hpp" -#include "Evolution/Systems/ForceFree/FiniteDifference/Tags.hpp" -#include "Evolution/Systems/ForceFree/Fluxes.hpp" -#include "Evolution/Systems/ForceFree/Subcell/ComputeFluxes.hpp" -#include "Evolution/Systems/ForceFree/System.hpp" -#include "NumericalAlgorithms/Spectral/Mesh.hpp" -#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" -#include "Utilities/CallWithDynamicType.hpp" -#include "Utilities/ErrorHandling/Assert.hpp" -#include "Utilities/Gsl.hpp" -#include "Utilities/TMPL.hpp" +/// \cond +namespace db { +class Access; +} // namespace db +class DataVector; +template +struct DirectionalId; +template +class DirectionalIdMap; +/// \endcond namespace ForceFree::subcell { @@ -64,221 +34,9 @@ namespace ForceFree::subcell { * if one element is doing DG, then we aren't at a shock. */ struct NeighborPackagedData { - template static DirectionalIdMap<3, DataVector> apply( - const db::DataBox& box, - const std::vector>& mortars_to_reconstruct_to) { - using evolved_vars_tags = typename System::variables_tag::tags_list; - using fluxes_tags = typename Fluxes::return_tags; - - // Subcell with moving mesh is not supported in the ForceFree system yet. - ASSERT(not db::get>(box).has_value(), - "Haven't yet added support for moving mesh to DG-subcell. This " - "should be easy to generalize, but we will want to consider " - "storing the mesh velocity on the faces instead of " - "re-slicing/projecting."); - - DirectionalIdMap<3, DataVector> neighbor_package_data{}; - if (mortars_to_reconstruct_to.empty()) { - return neighbor_package_data; - } - - // Project volume variables from DG to subcell mesh - const Mesh<3>& dg_mesh = db::get>(box); - const Mesh<3>& subcell_mesh = - db::get>(box); - const auto volume_vars_subcell = evolution::dg::subcell::fd::project( - db::get(box), dg_mesh, - subcell_mesh.extents()); - - // In addition to evolved vars, compute TildeJ on the subcell grid so that - // it could be sent to neighbor elements doing FD - DataVector buffer{ - Variables>::number_of_independent_components * - dg_mesh.number_of_grid_points()}; - const auto volume_tildej_subcell = [&box, &buffer, &dg_mesh, - &subcell_mesh]() { - Variables> var{buffer.data(), buffer.size()}; - for (size_t i = 0; i < 3; ++i) { - get(var).get(i) = db::get(box).get(i); - } - return get(evolution::dg::subcell::fd::project( - var, dg_mesh, subcell_mesh.extents())); - }(); - - const auto& ghost_subcell_data = - db::get>( - box); - - const ForceFree::fd::Reconstructor& recons = - db::get(box); - - const auto& boundary_correction = - db::get>(box); - using derived_boundary_corrections = - typename std::decay_t::creatable_classes; - - const auto& subcell_options = - db::get>(box); - - tmpl::for_each< - derived_boundary_corrections>([&box, &boundary_correction, &dg_mesh, - &mortars_to_reconstruct_to, - &neighbor_package_data, - &ghost_subcell_data, &recons, - &subcell_mesh, &subcell_options, - &volume_vars_subcell, - &volume_tildej_subcell]( - auto derived_correction_v) { - using DerivedCorrection = tmpl::type_from; - - using dg_package_field_tags = - typename DerivedCorrection::dg_package_field_tags; - using dg_package_data_temporary_tags = - typename DerivedCorrection::dg_package_data_temporary_tags; - - if (typeid(boundary_correction) == typeid(DerivedCorrection)) { - Variables packaged_data{0}; - - // Comprehensive tags list for computing fluxes and dg_package_data - using dg_package_data_argument_tags = tmpl::append< - evolved_vars_tags, fluxes_tags, - tmpl::remove_duplicates, - gr::Tags::SqrtDetSpatialMetric, - gr::Tags::InverseSpatialMetric, - evolution::dg::Actions::detail::NormalVector<3>>>>; - - // We allocate a Variables object with extra buffer for storing TildeJ. - // TildeJ is not an argument of dg_package_data() function, but it needs - // to be reconstructed for computing fluxes on faces. - Variables, - dg_package_data_argument_tags>> - vars_on_face_with_tildej_buffer{0}; - - for (const auto& mortar_id : mortars_to_reconstruct_to) { - const Direction<3>& direction = mortar_id.direction; - - // Switch to face-centered points - Index<3> extents = subcell_mesh.extents(); - ++extents[direction.dimension()]; - const size_t num_face_pts = subcell_mesh.extents() - .slice_away(direction.dimension()) - .product(); - vars_on_face_with_tildej_buffer.initialize(num_face_pts); - - // create a `view` of the face vars excluding the TildeJ tag, since - // boundary correction takes the flux F^i(TildeQ) as argument but not - // TildeJ itself. - auto vars_on_face = - vars_on_face_with_tildej_buffer - .template reference_subset(); - auto reconstructed_vars_on_face = - vars_on_face_with_tildej_buffer.template reference_subset< - ForceFree::fd::tags_list_for_reconstruction>(); - - // Copy spacetime vars over from volume. - using spacetime_vars_to_copy = tmpl::list< - gr::Tags::Lapse, - gr::Tags::Shift, - gr::Tags::SpatialMetric, - gr::Tags::SqrtDetSpatialMetric, - gr::Tags::InverseSpatialMetric>; - tmpl::for_each( - [&direction, &extents, &vars_on_face, - &spacetime_vars_on_faces = - db::get>(box)]( - auto tag_v) { - using tag = tmpl::type_from; - data_on_slice(make_not_null(&get(vars_on_face)), - get(gsl::at(spacetime_vars_on_faces, - direction.dimension())), - extents, direction.dimension(), - direction.side() == Side::Lower - ? 0 - : extents[direction.dimension()] - 1); - }); - - // Perform FD reconstruction of variables on cell interfaces. Note - // that we are using the vars object with TildeJ buffer so that - // TildeJ can be reconstructed in this step. - call_with_dynamic_type< - void, typename ForceFree::fd::Reconstructor::creatable_classes>( - &recons, [&element = db::get>(box), - &mortar_id, &ghost_subcell_data, &subcell_mesh, - &reconstructed_vars_on_face, &volume_vars_subcell, - &volume_tildej_subcell](const auto& reconstructor) { - reconstructor->reconstruct_fd_neighbor( - make_not_null(&reconstructed_vars_on_face), - volume_vars_subcell, volume_tildej_subcell, element, - ghost_subcell_data, subcell_mesh, mortar_id.direction); - }); - ForceFree::subcell::compute_fluxes( - make_not_null(&vars_on_face_with_tildej_buffer)); - - // Note: since the spacetime isn't dynamical we don't need to - // worry about different normal vectors on the different sides - // of the element. If the spacetime is dynamical, then the normal - // covector needs to be sent, or at least the inverse spatial metric - // from the neighbor so that the normal covector can be computed - // correctly. - tnsr::i normal_covector = - get>( - *db::get>( - box) - .at(mortar_id.direction)); - for (auto& t : normal_covector) { - t *= -1.0; - } - - // Note: Only need to do the projection in 2d and 3d, but GRFFE is - // always 3d currently. - const auto dg_normal_covector = normal_covector; - for (size_t i = 0; i < 3; ++i) { - normal_covector.get(i) = evolution::dg::subcell::fd::project( - dg_normal_covector.get(i), - dg_mesh.slice_away(mortar_id.direction.dimension()), - subcell_mesh.extents().slice_away( - mortar_id.direction.dimension())); - } - - // Compute the packaged data - packaged_data.initialize(num_face_pts); - using dg_package_data_projected_tags = - tmpl::append; - evolution::dg::Actions::detail::dg_package_data( - make_not_null(&packaged_data), - dynamic_cast(boundary_correction), - vars_on_face, normal_covector, {std::nullopt}, box, - typename DerivedCorrection::dg_package_data_volume_tags{}, - dg_package_data_projected_tags{}); - - // Reconstruct the DG solution. Really we should be solving the - // boundary correction and then reconstructing, but away from a shock - // this doesn't matter. - auto dg_packaged_data = evolution::dg::subcell::fd::reconstruct( - packaged_data, - dg_mesh.slice_away(mortar_id.direction.dimension()), - subcell_mesh.extents().slice_away( - mortar_id.direction.dimension()), - subcell_options.reconstruction_method()); - - // Make a view so we can use iterators with std::copy - DataVector dg_packaged_data_view{dg_packaged_data.data(), - dg_packaged_data.size()}; - neighbor_package_data[mortar_id] = - DataVector{dg_packaged_data.size()}; - std::copy(dg_packaged_data_view.begin(), dg_packaged_data_view.end(), - neighbor_package_data[mortar_id].begin()); - } - } - }); - - return neighbor_package_data; - } + const db::Access& box, + const std::vector>& mortars_to_reconstruct_to); }; } // namespace ForceFree::subcell diff --git a/src/Evolution/Systems/GeneralizedHarmonic/BoundaryCorrections/UpwindPenalty.hpp b/src/Evolution/Systems/GeneralizedHarmonic/BoundaryCorrections/UpwindPenalty.hpp index 17fd560ded12..2c82bc939adc 100644 --- a/src/Evolution/Systems/GeneralizedHarmonic/BoundaryCorrections/UpwindPenalty.hpp +++ b/src/Evolution/Systems/GeneralizedHarmonic/BoundaryCorrections/UpwindPenalty.hpp @@ -223,6 +223,7 @@ class UpwindPenalty final : public BoundaryCorrection { gr::Tags::Lapse, gr::Tags::Shift>; using dg_package_data_primitive_tags = tmpl::list<>; using dg_package_data_volume_tags = tmpl::list<>; + using dg_boundary_terms_volume_tags = tmpl::list<>; double dg_package_data( gsl::not_null*> diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryCorrections/ProductOfCorrections.hpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryCorrections/ProductOfCorrections.hpp index f6b702de8d62..85fe74b9a067 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryCorrections/ProductOfCorrections.hpp +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryCorrections/ProductOfCorrections.hpp @@ -152,6 +152,10 @@ class ProductOfCorrections final : public BoundaryCorrection { typename DerivedGhCorrection::dg_package_data_volume_tags, typename DerivedValenciaCorrection::dg_package_data_volume_tags>; + using dg_boundary_terms_volume_tags = tmpl::append< + typename DerivedGhCorrection::dg_boundary_terms_volume_tags, + typename DerivedValenciaCorrection::dg_boundary_terms_volume_tags>; + using derived_product_correction_impl = detail::ProductOfCorrectionsImpl< DerivedGhCorrection, DerivedValenciaCorrection, typename DerivedGhCorrection::dg_package_field_tags, diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/CMakeLists.txt b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/CMakeLists.txt index 174142aba556..6ca5807d7261 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/CMakeLists.txt +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/CMakeLists.txt @@ -5,6 +5,7 @@ spectre_target_sources( ${LIBRARY} PRIVATE FixConservativesAndComputePrims.cpp + NeighborPackagedData.cpp PrimitiveGhostData.cpp PrimsAfterRollback.cpp ResizeAndComputePrimitives.cpp diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/NeighborPackagedData.cpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/NeighborPackagedData.cpp new file mode 100644 index 000000000000..6f383ba3537e --- /dev/null +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/NeighborPackagedData.cpp @@ -0,0 +1,334 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/NeighborPackagedData.hpp" + +#include +#include +#include +#include +#include + +#include "DataStructures/DataBox/Access.hpp" +#include "DataStructures/DataBox/PrefixHelpers.hpp" +#include "DataStructures/DataBox/Prefixes.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Index.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "DataStructures/VariablesTag.hpp" +#include "Domain/Structure/Direction.hpp" +#include "Domain/Structure/DirectionalIdMap.hpp" +#include "Domain/Structure/Element.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Domain/Tags.hpp" +#include "Domain/TagsTimeDependent.hpp" +#include "Evolution/BoundaryCorrectionTags.hpp" +#include "Evolution/DgSubcell/NeighborReconstructedFaceSolution.tpp" +#include "Evolution/DgSubcell/Projection.hpp" +#include "Evolution/DgSubcell/Reconstruction.hpp" +#include "Evolution/DgSubcell/ReconstructionMethod.hpp" +#include "Evolution/DgSubcell/SubcellOptions.hpp" +#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" +#include "Evolution/DgSubcell/Tags/Mesh.hpp" +#include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp" +#include "Evolution/DgSubcell/Tags/SubcellOptions.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryCorrections/BoundaryCorrection.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryCorrections/Factory.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Factory.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Tag.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/System.hpp" +#include "Evolution/Systems/GrMhd/GhValenciaDivClean/Tags.hpp" +#include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/ComputeFluxes.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp" +#include "PointwiseFunctions/Hydro/Tags.hpp" +#include "Utilities/CallWithDynamicType.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/TMPL.hpp" + +namespace grmhd::GhValenciaDivClean::subcell { +DirectionalIdMap<3, DataVector> NeighborPackagedData::apply( + const db::Access& box, + const std::vector>& mortars_to_reconstruct_to) { + using evolved_vars_tag = typename System::variables_tag; + using evolved_vars_tags = typename evolved_vars_tag::tags_list; + using prim_tags = typename System::primitive_variables_tag::tags_list; + using recons_prim_tags = tmpl::push_front>>; + using fluxes_tags = + db::wrap_tags_in<::Tags::Flux, typename System::flux_variables, + tmpl::size_t<3>, Frame::Inertial>; + using grmhd_evolved_vars_tag = + typename grmhd::ValenciaDivClean::System::variables_tag; + using grmhd_evolved_vars_tags = typename grmhd_evolved_vars_tag::tags_list; + + DirectionalIdMap<3, DataVector> neighbor_package_data{}; + if (mortars_to_reconstruct_to.empty()) { + return neighbor_package_data; + } + + const auto& ghost_subcell_data = + db::get>(box); + const Mesh<3>& subcell_mesh = + db::get>(box); + const Mesh<3>& dg_mesh = db::get>(box); + const auto& subcell_options = + db::get>(box); + const auto& evolved_vars = db::get(box); + + const Variables volume_spacetime_vars = + evolution::dg::subcell::fd::project( + Variables{ + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast) + const_cast( + get>( + evolved_vars)[0] + .data()), + Variables:: + number_of_independent_components * + evolved_vars.number_of_grid_points()}, + dg_mesh, subcell_mesh.extents()); + + // Note: we need to compare if projecting the entire mesh or only ghost + // zones needed is faster. This probably depends on the number of neighbors + // we have doing FD. + const auto volume_prims = evolution::dg::subcell::fd::project( + db::get(box), dg_mesh, + subcell_mesh.extents()); + + const auto& recons = + db::get(box); + const auto& base_boundary_correction = + db::get>(box); + using derived_boundary_corrections = typename std::decay_t< + decltype(base_boundary_correction)>::creatable_classes; + call_with_dynamic_type( + &base_boundary_correction, + [&box, &dg_mesh, &mortars_to_reconstruct_to, &neighbor_package_data, + &ghost_subcell_data, &recons, &subcell_mesh, &subcell_options, + &volume_prims, &volume_spacetime_vars](const auto* gh_grmhd_correction) { + using DerivedCorrection = std::decay_t; + const auto& boundary_correction = + dynamic_cast(*gh_grmhd_correction); + + using dg_package_data_temporary_tags = + typename DerivedCorrection::dg_package_data_temporary_tags; + using dg_package_data_argument_tags = + tmpl::append, + gr::Tags::SqrtDetSpatialMetric, + gr::Tags::InverseSpatialMetric, + evolution::dg::Actions::detail::NormalVector<3>>>>; + + const auto& element = db::get>(box); + const auto& eos = get(box); + + using dg_package_field_tags = + typename DerivedCorrection::dg_package_field_tags; + Variables vars_on_face{0}; + Variables packaged_data{0}; + for (const auto& mortar_id : mortars_to_reconstruct_to) { + const Direction<3>& direction = mortar_id.direction; + + const Mesh<2> dg_face_mesh = + dg_mesh.slice_away(direction.dimension()); + Index<3> extents = subcell_mesh.extents(); + // Switch to face-centered instead of cell-centered points on the + // FD. There are num_cell_centered+1 face-centered points. + ++extents[direction.dimension()]; + const Index<2> subcell_face_extents = + extents.slice_away(direction.dimension()); + + // Computed prims and cons on face via reconstruction + const size_t num_face_pts = subcell_mesh.extents() + .slice_away(direction.dimension()) + .product(); + vars_on_face.initialize(num_face_pts); + + call_with_dynamic_type( + &recons, [&element, &eos, &mortar_id, &ghost_subcell_data, + &subcell_mesh, &vars_on_face, &volume_prims, + &volume_spacetime_vars](const auto& reconstructor) { + reconstructor->reconstruct_fd_neighbor( + make_not_null(&vars_on_face), volume_prims, + volume_spacetime_vars, eos, element, ghost_subcell_data, + subcell_mesh, mortar_id.direction); + }); + + // Get the mesh velocity if needed + const std::optional> + mesh_velocity_dg = db::get>(box); + std::optional> + mesh_velocity_on_subcell_face = {}; + if (mesh_velocity_dg.has_value()) { + // Slice data on current face + tnsr::I mesh_velocity_on_dg_face = + data_on_slice( + mesh_velocity_dg.value(), dg_mesh.extents(), + direction.dimension(), + direction.side() == Side::Lower + ? 0 + : (dg_mesh.extents(direction.dimension()) - 1)); + + mesh_velocity_on_subcell_face = + tnsr::I{num_face_pts}; + + for (size_t i = 0; i < 3; i++) { + evolution::dg::subcell::fd::project( + make_not_null(&mesh_velocity_on_subcell_face.value().get(i)), + mesh_velocity_on_dg_face.get(i), dg_face_mesh, + subcell_face_extents); + } + } + + grmhd::ValenciaDivClean::subcell::compute_fluxes( + make_not_null(&vars_on_face)); + + if (mesh_velocity_on_subcell_face.has_value()) { + tmpl::for_each( + [&vars_on_face, &mesh_velocity_on_subcell_face](auto tag_v) { + using tag = tmpl::type_from; + using flux_tag = + ::Tags::Flux, Frame::Inertial>; + using FluxTensor = typename flux_tag::type; + const auto& var = get(vars_on_face); + auto& flux = get(vars_on_face); + for (size_t storage_index = 0; storage_index < var.size(); + ++storage_index) { + const auto tensor_index = + var.get_tensor_index(storage_index); + for (size_t j = 0; j < 3; j++) { + const auto flux_storage_index = + FluxTensor::get_storage_index( + prepend(tensor_index, j)); + flux[flux_storage_index] -= + mesh_velocity_on_subcell_face.value().get(j) * + var[storage_index]; + } + } + }); + } + + // Copy over Gamma1 and Gamma2. Future considerations: + // - For LTS do we instead just recompute it since the times might + // not be aligned? + // - Instead of slicing, can we copy the data from the local + // packaged data? + { + Scalar gamma_on_dg_face = data_on_slice( + get(box), + dg_mesh.extents(), direction.dimension(), + direction.side() == Side::Lower + ? 0 + : (dg_mesh.extents(direction.dimension()) - 1)); + evolution::dg::subcell::fd::project( + make_not_null( + &get(get( + vars_on_face))), + get(gamma_on_dg_face), dg_face_mesh, subcell_face_extents); + data_on_slice( + make_not_null(&gamma_on_dg_face), + get(box), + dg_mesh.extents(), direction.dimension(), + direction.side() == Side::Lower + ? 0 + : (dg_mesh.extents(direction.dimension()) - 1)); + evolution::dg::subcell::fd::project( + make_not_null( + &get(get( + vars_on_face))), + get(gamma_on_dg_face), dg_face_mesh, subcell_face_extents); + } + + // Compute the neighbor's normal covector and normalize it. + tnsr::i normal_covector = + get>( + *db::get>( + box) + .at(mortar_id.direction)); + for (auto& t : normal_covector) { + t *= -1.0; + } + // Note: Only need to do the projection in 2d and 3d, but GRMHD is + // always 3d currently. + for (size_t i = 0; i < 3; ++i) { + normal_covector.get(i) = evolution::dg::subcell::fd::project( + normal_covector.get(i), + dg_mesh.slice_away(mortar_id.direction.dimension()), + subcell_mesh.extents().slice_away( + mortar_id.direction.dimension())); + } + // Need to renormalize the normal vector with the neighbor's + // inverse spatial metric. + const auto& inverse_spatial_metric = + get>(vars_on_face); + const auto normal_magnitude = + magnitude(normal_covector, inverse_spatial_metric); + for (size_t i = 0; i < 3; ++i) { + normal_covector.get(i) /= get(normal_magnitude); + } + + auto& normal_vector = + get>( + vars_on_face); + for (size_t i = 0; i < 3; ++i) { + normal_vector.get(i) = + normal_covector.get(0) * inverse_spatial_metric.get(i, 0) + + normal_covector.get(1) * inverse_spatial_metric.get(i, 1) + + normal_covector.get(2) * inverse_spatial_metric.get(i, 2); + } + + // Compute the packaged data + packaged_data.initialize(num_face_pts); + using dg_package_data_projected_tags = tmpl::append< + evolved_vars_tags, fluxes_tags, dg_package_data_temporary_tags, + typename DerivedCorrection::dg_package_data_primitive_tags>; + evolution::dg::Actions::detail::dg_package_data( + make_not_null(&packaged_data), + dynamic_cast(boundary_correction), + vars_on_face, normal_covector, mesh_velocity_on_subcell_face, box, + typename DerivedCorrection::dg_package_data_volume_tags{}, + dg_package_data_projected_tags{}); + + // Reconstruct the DG solution. + // Really we should be solving the boundary correction and + // then reconstructing, but away from a shock this doesn't + // matter. + DataVector dg_data{ + Variables< + dg_package_field_tags>::number_of_independent_components * + dg_face_mesh.number_of_grid_points()}; + Variables dg_packaged_data{dg_data.data(), + dg_data.size()}; + evolution::dg::subcell::fd::reconstruct( + make_not_null(&dg_packaged_data), packaged_data, dg_face_mesh, + subcell_mesh.extents().slice_away( + mortar_id.direction.dimension()), + subcell_options.reconstruction_method()); + + neighbor_package_data[mortar_id] = std::move(dg_data); + } + }); + + return neighbor_package_data; +} +} // namespace grmhd::GhValenciaDivClean::subcell + +template void evolution::dg::subcell::neighbor_reconstructed_face_solution< + 3, grmhd::GhValenciaDivClean::subcell::NeighborPackagedData>( + gsl::not_null box, + gsl::not_null, Mesh<2>, std::optional, + std::optional, ::TimeStepId, int>>>*> + received_temporal_id_and_data); diff --git a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/NeighborPackagedData.hpp b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/NeighborPackagedData.hpp index bb8c018b619b..a08082ff5cff 100644 --- a/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/NeighborPackagedData.hpp +++ b/src/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/NeighborPackagedData.hpp @@ -3,52 +3,19 @@ #pragma once -#include #include -#include -#include -#include #include -#include "DataStructures/DataBox/DataBox.hpp" -#include "DataStructures/DataBox/PrefixHelpers.hpp" -#include "DataStructures/DataBox/Prefixes.hpp" -#include "DataStructures/DataVector.hpp" -#include "DataStructures/Index.hpp" -#include "DataStructures/Tensor/Tensor.hpp" -#include "DataStructures/Variables.hpp" -#include "DataStructures/VariablesTag.hpp" -#include "Domain/Structure/Direction.hpp" -#include "Domain/Structure/DirectionalIdMap.hpp" -#include "Domain/Structure/Element.hpp" -#include "Domain/Structure/ElementId.hpp" -#include "Domain/Tags.hpp" -#include "Domain/TagsTimeDependent.hpp" -#include "Evolution/BoundaryCorrectionTags.hpp" -#include "Evolution/DgSubcell/Projection.hpp" -#include "Evolution/DgSubcell/Reconstruction.hpp" -#include "Evolution/DgSubcell/ReconstructionMethod.hpp" -#include "Evolution/DgSubcell/SubcellOptions.hpp" -#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" -#include "Evolution/DgSubcell/Tags/Mesh.hpp" -#include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp" -#include "Evolution/DgSubcell/Tags/SubcellOptions.hpp" -#include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp" -#include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp" -#include "Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryCorrections/BoundaryCorrection.hpp" -#include "Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryCorrections/Factory.hpp" -#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp" -#include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Tag.hpp" -#include "Evolution/Systems/GrMhd/GhValenciaDivClean/System.hpp" -#include "Evolution/Systems/GrMhd/GhValenciaDivClean/Tags.hpp" -#include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/ComputeFluxes.hpp" -#include "NumericalAlgorithms/Spectral/Mesh.hpp" -#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" -#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp" -#include "PointwiseFunctions/Hydro/Tags.hpp" -#include "Utilities/CallWithDynamicType.hpp" -#include "Utilities/Gsl.hpp" -#include "Utilities/TMPL.hpp" +/// \cond +namespace db { +class Access; +} // namespace db +class DataVector; +template +struct DirectionalId; +template +class DirectionalIdMap; +/// \endcond namespace grmhd::GhValenciaDivClean::subcell { /*! @@ -69,281 +36,8 @@ namespace grmhd::GhValenciaDivClean::subcell { * velocity on the faces instead of re-slicing/projecting. */ struct NeighborPackagedData { - template static DirectionalIdMap<3, DataVector> apply( - const db::DataBox& box, - const std::vector>& mortars_to_reconstruct_to) { - using evolved_vars_tag = typename System::variables_tag; - using evolved_vars_tags = typename evolved_vars_tag::tags_list; - using prim_tags = typename System::primitive_variables_tag::tags_list; - using recons_prim_tags = tmpl::push_front>>; - using fluxes_tags = - db::wrap_tags_in<::Tags::Flux, typename System::flux_variables, - tmpl::size_t<3>, Frame::Inertial>; - using grmhd_evolved_vars_tag = - typename grmhd::ValenciaDivClean::System::variables_tag; - using grmhd_evolved_vars_tags = typename grmhd_evolved_vars_tag::tags_list; - - DirectionalIdMap<3, DataVector> neighbor_package_data{}; - if (mortars_to_reconstruct_to.empty()) { - return neighbor_package_data; - } - - const auto& ghost_subcell_data = - db::get>( - box); - const Mesh<3>& subcell_mesh = - db::get>(box); - const Mesh<3>& dg_mesh = db::get>(box); - const auto& subcell_options = - db::get>(box); - const auto& evolved_vars = db::get(box); - - const Variables volume_spacetime_vars = - evolution::dg::subcell::fd::project( - Variables{ - // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast) - const_cast( - get>( - evolved_vars)[0] - .data()), - Variables:: - number_of_independent_components * - evolved_vars.number_of_grid_points()}, - dg_mesh, subcell_mesh.extents()); - - // Note: we need to compare if projecting the entire mesh or only ghost - // zones needed is faster. This probably depends on the number of neighbors - // we have doing FD. - const auto volume_prims = evolution::dg::subcell::fd::project( - db::get(box), dg_mesh, - subcell_mesh.extents()); - - const auto& recons = - db::get(box); - const auto& base_boundary_correction = - db::get>(box); - using derived_boundary_corrections = typename std::decay_t::creatable_classes; - call_with_dynamic_type( - &base_boundary_correction, - [&box, &dg_mesh, &mortars_to_reconstruct_to, &neighbor_package_data, - &ghost_subcell_data, &recons, &subcell_mesh, &subcell_options, - &volume_prims, - &volume_spacetime_vars](const auto* gh_grmhd_correction) { - using DerivedCorrection = - std::decay_t; - const auto& boundary_correction = - dynamic_cast(*gh_grmhd_correction); - - using dg_package_data_temporary_tags = - typename DerivedCorrection::dg_package_data_temporary_tags; - using dg_package_data_argument_tags = tmpl::append< - evolved_vars_tags, recons_prim_tags, fluxes_tags, - tmpl::remove_duplicates, - gr::Tags::SqrtDetSpatialMetric, - gr::Tags::InverseSpatialMetric, - evolution::dg::Actions::detail::NormalVector<3>>>>; - - const auto& element = db::get>(box); - const auto& eos = get(box); - - using dg_package_field_tags = - typename DerivedCorrection::dg_package_field_tags; - Variables vars_on_face{0}; - Variables packaged_data{0}; - for (const auto& mortar_id : mortars_to_reconstruct_to) { - const Direction<3>& direction = mortar_id.direction; - - const Mesh<2> dg_face_mesh = - dg_mesh.slice_away(direction.dimension()); - Index<3> extents = subcell_mesh.extents(); - // Switch to face-centered instead of cell-centered points on the - // FD. There are num_cell_centered+1 face-centered points. - ++extents[direction.dimension()]; - const Index<2> subcell_face_extents = - extents.slice_away(direction.dimension()); - - // Computed prims and cons on face via reconstruction - const size_t num_face_pts = subcell_mesh.extents() - .slice_away(direction.dimension()) - .product(); - vars_on_face.initialize(num_face_pts); - - call_with_dynamic_type( - &recons, [&element, &eos, &mortar_id, &ghost_subcell_data, - &subcell_mesh, &vars_on_face, &volume_prims, - &volume_spacetime_vars](const auto& reconstructor) { - reconstructor->reconstruct_fd_neighbor( - make_not_null(&vars_on_face), volume_prims, - volume_spacetime_vars, eos, element, ghost_subcell_data, - subcell_mesh, mortar_id.direction); - }); - - // Get the mesh velocity if needed - const std::optional> - mesh_velocity_dg = db::get>(box); - std::optional> - mesh_velocity_on_subcell_face = {}; - if (mesh_velocity_dg.has_value()) { - // Slice data on current face - tnsr::I mesh_velocity_on_dg_face = - data_on_slice( - mesh_velocity_dg.value(), dg_mesh.extents(), - direction.dimension(), - direction.side() == Side::Lower - ? 0 - : (dg_mesh.extents(direction.dimension()) - 1)); - - mesh_velocity_on_subcell_face = - tnsr::I{num_face_pts}; - - for (size_t i = 0; i < 3; i++) { - evolution::dg::subcell::fd::project( - make_not_null( - &mesh_velocity_on_subcell_face.value().get(i)), - mesh_velocity_on_dg_face.get(i), dg_face_mesh, - subcell_face_extents); - } - } - - grmhd::ValenciaDivClean::subcell::compute_fluxes( - make_not_null(&vars_on_face)); - - if (mesh_velocity_on_subcell_face.has_value()) { - tmpl::for_each( - [&vars_on_face, &mesh_velocity_on_subcell_face](auto tag_v) { - using tag = tmpl::type_from; - using flux_tag = - ::Tags::Flux, Frame::Inertial>; - using FluxTensor = typename flux_tag::type; - const auto& var = get(vars_on_face); - auto& flux = get(vars_on_face); - for (size_t storage_index = 0; storage_index < var.size(); - ++storage_index) { - const auto tensor_index = - var.get_tensor_index(storage_index); - for (size_t j = 0; j < 3; j++) { - const auto flux_storage_index = - FluxTensor::get_storage_index( - prepend(tensor_index, j)); - flux[flux_storage_index] -= - mesh_velocity_on_subcell_face.value().get(j) * - var[storage_index]; - } - } - }); - } - - // Copy over Gamma1 and Gamma2. Future considerations: - // - For LTS do we instead just recompute it since the times might - // not be aligned? - // - Instead of slicing, can we copy the data from the local - // packaged data? - { - Scalar gamma_on_dg_face = data_on_slice( - get(box), - dg_mesh.extents(), direction.dimension(), - direction.side() == Side::Lower - ? 0 - : (dg_mesh.extents(direction.dimension()) - 1)); - evolution::dg::subcell::fd::project( - make_not_null( - &get(get( - vars_on_face))), - get(gamma_on_dg_face), dg_face_mesh, subcell_face_extents); - data_on_slice( - make_not_null(&gamma_on_dg_face), - get(box), - dg_mesh.extents(), direction.dimension(), - direction.side() == Side::Lower - ? 0 - : (dg_mesh.extents(direction.dimension()) - 1)); - evolution::dg::subcell::fd::project( - make_not_null( - &get(get( - vars_on_face))), - get(gamma_on_dg_face), dg_face_mesh, subcell_face_extents); - } - - // Compute the neighbor's normal covector and normalize it. - tnsr::i normal_covector = get< - evolution::dg::Tags::NormalCovector<3>>( - *db::get>( - box) - .at(mortar_id.direction)); - for (auto& t : normal_covector) { - t *= -1.0; - } - // Note: Only need to do the projection in 2d and 3d, but GRMHD is - // always 3d currently. - for (size_t i = 0; i < 3; ++i) { - normal_covector.get(i) = evolution::dg::subcell::fd::project( - normal_covector.get(i), - dg_mesh.slice_away(mortar_id.direction.dimension()), - subcell_mesh.extents().slice_away( - mortar_id.direction.dimension())); - } - // Need to renormalize the normal vector with the neighbor's - // inverse spatial metric. - const auto& inverse_spatial_metric = - get>( - vars_on_face); - const auto normal_magnitude = - magnitude(normal_covector, inverse_spatial_metric); - for (size_t i = 0; i < 3; ++i) { - normal_covector.get(i) /= get(normal_magnitude); - } - - auto& normal_vector = - get>( - vars_on_face); - for (size_t i = 0; i < 3; ++i) { - normal_vector.get(i) = - normal_covector.get(0) * inverse_spatial_metric.get(i, 0) + - normal_covector.get(1) * inverse_spatial_metric.get(i, 1) + - normal_covector.get(2) * inverse_spatial_metric.get(i, 2); - } - - // Compute the packaged data - packaged_data.initialize(num_face_pts); - using dg_package_data_projected_tags = tmpl::append< - evolved_vars_tags, fluxes_tags, dg_package_data_temporary_tags, - typename DerivedCorrection::dg_package_data_primitive_tags>; - evolution::dg::Actions::detail::dg_package_data( - make_not_null(&packaged_data), - dynamic_cast(boundary_correction), - vars_on_face, normal_covector, mesh_velocity_on_subcell_face, - box, typename DerivedCorrection::dg_package_data_volume_tags{}, - dg_package_data_projected_tags{}); - - // Reconstruct the DG solution. - // Really we should be solving the boundary correction and - // then reconstructing, but away from a shock this doesn't - // matter. - DataVector dg_data{ - Variables< - dg_package_field_tags>::number_of_independent_components * - dg_face_mesh.number_of_grid_points()}; - Variables dg_packaged_data{dg_data.data(), - dg_data.size()}; - evolution::dg::subcell::fd::reconstruct( - make_not_null(&dg_packaged_data), packaged_data, dg_face_mesh, - subcell_mesh.extents().slice_away( - mortar_id.direction.dimension()), - subcell_options.reconstruction_method()); - - neighbor_package_data[mortar_id] = std::move(dg_data); - } - }); - - return neighbor_package_data; - } + const db::Access& box, + const std::vector>& mortars_to_reconstruct_to); }; } // namespace grmhd::GhValenciaDivClean::subcell diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/Hll.hpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/Hll.hpp index a1b1d42485f3..273782034a82 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/Hll.hpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/Hll.hpp @@ -127,6 +127,7 @@ class Hll final : public BoundaryCorrection { tmpl::list, gr::Tags::Shift>; using dg_package_data_primitive_tags = tmpl::list<>; using dg_package_data_volume_tags = tmpl::list<>; + using dg_boundary_terms_volume_tags = tmpl::list<>; static double dg_package_data( gsl::not_null*> packaged_tilde_d, diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/Rusanov.hpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/Rusanov.hpp index ff01d3c782b8..b7eff38eb6ba 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/Rusanov.hpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/Rusanov.hpp @@ -107,6 +107,7 @@ class Rusanov final : public BoundaryCorrection { tmpl::list, gr::Tags::Shift>; using dg_package_data_primitive_tags = tmpl::list<>; using dg_package_data_volume_tags = tmpl::list<>; + using dg_boundary_terms_volume_tags = tmpl::list<>; static double dg_package_data( gsl::not_null*> packaged_tilde_d, diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/CMakeLists.txt b/src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/CMakeLists.txt index 099b3b200e98..41d20b6e26c3 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/CMakeLists.txt +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/CMakeLists.txt @@ -5,6 +5,7 @@ spectre_target_sources( ${LIBRARY} PRIVATE FixConservativesAndComputePrims.cpp + NeighborPackagedData.cpp PrimitiveGhostData.cpp PrimsAfterRollback.cpp ResizeAndComputePrimitives.cpp diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/NeighborPackagedData.cpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/NeighborPackagedData.cpp new file mode 100644 index 000000000000..0163f56030d0 --- /dev/null +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/NeighborPackagedData.cpp @@ -0,0 +1,282 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/NeighborPackagedData.hpp" + +#include +#include +#include +#include +#include + +#include "DataStructures/DataBox/Access.hpp" +#include "DataStructures/DataBox/PrefixHelpers.hpp" +#include "DataStructures/DataBox/Prefixes.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Index.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "DataStructures/VariablesTag.hpp" +#include "Domain/Structure/Direction.hpp" +#include "Domain/Structure/DirectionalIdMap.hpp" +#include "Domain/Structure/Element.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Domain/Tags.hpp" +#include "Domain/TagsTimeDependent.hpp" +#include "Evolution/BoundaryCorrectionTags.hpp" +#include "Evolution/DgSubcell/NeighborReconstructedFaceSolution.tpp" +#include "Evolution/DgSubcell/Projection.hpp" +#include "Evolution/DgSubcell/Reconstruction.hpp" +#include "Evolution/DgSubcell/ReconstructionMethod.hpp" +#include "Evolution/DgSubcell/SubcellOptions.hpp" +#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" +#include "Evolution/DgSubcell/Tags/Mesh.hpp" +#include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp" +#include "Evolution/DgSubcell/Tags/SubcellOptions.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp" +#include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/BoundaryCorrection.hpp" +#include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/Factory.hpp" +#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Factory.hpp" +#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.hpp" +#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp" +#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Tag.hpp" +#include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/ComputeFluxes.hpp" +#include "Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp" +#include "PointwiseFunctions/Hydro/Tags.hpp" +#include "Utilities/CallWithDynamicType.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/TMPL.hpp" + +namespace grmhd::ValenciaDivClean::subcell { +DirectionalIdMap<3, DataVector> NeighborPackagedData::apply( + const db::Access& box, + const std::vector>& mortars_to_reconstruct_to) { + using evolved_vars_tag = typename System::variables_tag; + using evolved_vars_tags = typename evolved_vars_tag::tags_list; + using fluxes_tags = db::wrap_tags_in<::Tags::Flux, evolved_vars_tags, + tmpl::size_t<3>, Frame::Inertial>; + + ASSERT(not db::get>(box).has_value(), + "Haven't yet added support for moving mesh to DG-subcell. This " + "should be easy to generalize, but we will want to consider " + "storing the mesh velocity on the faces instead of " + "re-slicing/projecting."); + + DirectionalIdMap<3, DataVector> neighbor_package_data{}; + if (mortars_to_reconstruct_to.empty()) { + return neighbor_package_data; + } + + const auto& ghost_subcell_data = + db::get>(box); + const Mesh<3>& subcell_mesh = + db::get>(box); + const Mesh<3>& dg_mesh = db::get>(box); + const auto& subcell_options = + db::get>(box); + + // Note: we need to compare if projecting the entire mesh or only ghost + // zones needed is faster. This probably depends on the number of neighbors + // we have doing FD. + const auto volume_prims = evolution::dg::subcell::fd::project( + db::get(box), dg_mesh, + subcell_mesh.extents()); + + const auto& recons = + db::get(box); + const auto& boundary_correction = + db::get>(box); + using derived_boundary_corrections = + typename std::decay_t::creatable_classes; + tmpl::for_each([&box, &boundary_correction, + &dg_mesh, + &mortars_to_reconstruct_to, + &neighbor_package_data, + &ghost_subcell_data, &recons, + &subcell_mesh, &subcell_options, + &volume_prims]( + auto derived_correction_v) { + using DerivedCorrection = tmpl::type_from; + if (typeid(boundary_correction) == typeid(DerivedCorrection)) { + using dg_package_data_temporary_tags = + typename DerivedCorrection::dg_package_data_temporary_tags; + using dg_package_data_argument_tags = fd::tags_list_for_reconstruct; + + const auto& element = db::get>(box); + const auto& eos = get(box); + + using dg_package_field_tags = + typename DerivedCorrection::dg_package_field_tags; + Variables vars_on_face{0}; + Variables packaged_data{0}; + for (const auto& mortar_id : mortars_to_reconstruct_to) { + const Direction<3>& direction = mortar_id.direction; + + const Mesh<2> dg_face_mesh = dg_mesh.slice_away(direction.dimension()); + Index<3> extents = subcell_mesh.extents(); + // Switch to face-centered instead of cell-centered points on the FD. + // There are num_cell_centered+1 face-centered points. + ++extents[direction.dimension()]; + const Index<2> subcell_face_extents = + extents.slice_away(direction.dimension()); + + // Computed prims and cons on face via reconstruction + const size_t num_face_pts = + subcell_mesh.extents().slice_away(direction.dimension()).product(); + vars_on_face.initialize(num_face_pts); + // Copy spacetime vars over from volume. + using spacetime_vars_to_copy = + tmpl::list, + gr::Tags::Shift, + gr::Tags::SpatialMetric, + gr::Tags::SqrtDetSpatialMetric, + gr::Tags::InverseSpatialMetric>; + tmpl::for_each( + [&direction, &extents, &vars_on_face, + &spacetime_vars_on_faces = + db::get>(box)]( + auto tag_v) { + using tag = tmpl::type_from; + data_on_slice(make_not_null(&get(vars_on_face)), + get(gsl::at(spacetime_vars_on_faces, + direction.dimension())), + extents, direction.dimension(), + direction.side() == Side::Lower + ? 0 + : extents[direction.dimension()] - 1); + }); + + call_with_dynamic_type( + &recons, + [&element, &eos, &mortar_id, &ghost_subcell_data, &subcell_mesh, + &vars_on_face, &volume_prims](const auto& reconstructor) { + reconstructor->reconstruct_fd_neighbor( + make_not_null(&vars_on_face), volume_prims, eos, element, + ghost_subcell_data, subcell_mesh, mortar_id.direction); + }); + + // Get the mesh velocity if needed + const std::optional> + mesh_velocity_dg = db::get>(box); + std::optional> + mesh_velocity_on_subcell_face = {}; + if (mesh_velocity_dg.has_value()) { + // Slice data on current face + tnsr::I mesh_velocity_on_dg_face = + data_on_slice(mesh_velocity_dg.value(), dg_mesh.extents(), + direction.dimension(), + direction.side() == Side::Lower + ? 0 + : (dg_mesh.extents(direction.dimension()) - 1)); + + mesh_velocity_on_subcell_face = + tnsr::I{num_face_pts}; + + for (size_t i = 0; i < 3; i++) { + evolution::dg::subcell::fd::project( + make_not_null(&mesh_velocity_on_subcell_face.value().get(i)), + mesh_velocity_on_dg_face.get(i), dg_face_mesh, + subcell_face_extents); + } + } + + grmhd::ValenciaDivClean::subcell::compute_fluxes( + make_not_null(&vars_on_face)); + + if (mesh_velocity_on_subcell_face.has_value()) { + tmpl::for_each( + [&vars_on_face, &mesh_velocity_on_subcell_face](auto tag_v) { + using tag = tmpl::type_from; + using flux_tag = + ::Tags::Flux, Frame::Inertial>; + using FluxTensor = typename flux_tag::type; + const auto& var = get(vars_on_face); + auto& flux = get(vars_on_face); + for (size_t storage_index = 0; storage_index < var.size(); + ++storage_index) { + const auto tensor_index = var.get_tensor_index(storage_index); + for (size_t j = 0; j < 3; j++) { + const auto flux_storage_index = + FluxTensor::get_storage_index(prepend(tensor_index, j)); + flux[flux_storage_index] -= + mesh_velocity_on_subcell_face.value().get(j) * + var[storage_index]; + } + } + }); + } + + // Note: since the spacetime isn't dynamical we don't need to + // worry about different normal vectors on the different sides + // of the element. If the spacetime is dynamical, then the normal + // covector needs to be sent, or at least the inverse spatial metric + // from the neighbor so that the normal covector can be computed + // correctly. + tnsr::i normal_covector = + get>( + *db::get>( + box) + .at(mortar_id.direction)); + for (auto& t : normal_covector) { + t *= -1.0; + } + // Note: Only need to do the projection in 2d and 3d, but GRMHD is + // always 3d currently. + const auto dg_normal_covector = normal_covector; + for (size_t i = 0; i < 3; ++i) { + normal_covector.get(i) = evolution::dg::subcell::fd::project( + dg_normal_covector.get(i), + dg_mesh.slice_away(mortar_id.direction.dimension()), + subcell_mesh.extents().slice_away( + mortar_id.direction.dimension())); + } + + // Compute the packaged data + packaged_data.initialize(num_face_pts); + using dg_package_data_projected_tags = tmpl::append< + evolved_vars_tags, fluxes_tags, dg_package_data_temporary_tags, + typename DerivedCorrection::dg_package_data_primitive_tags>; + evolution::dg::Actions::detail::dg_package_data( + make_not_null(&packaged_data), + dynamic_cast(boundary_correction), + vars_on_face, normal_covector, mesh_velocity_on_subcell_face, box, + typename DerivedCorrection::dg_package_data_volume_tags{}, + dg_package_data_projected_tags{}); + + // Reconstruct the DG solution. + // Really we should be solving the boundary correction and + // then reconstructing, but away from a shock this doesn't + // matter. + auto dg_packaged_data = evolution::dg::subcell::fd::reconstruct( + packaged_data, dg_mesh.slice_away(mortar_id.direction.dimension()), + subcell_mesh.extents().slice_away(mortar_id.direction.dimension()), + subcell_options.reconstruction_method()); + // Make a view so we can use iterators with std::copy + DataVector dg_packaged_data_view{dg_packaged_data.data(), + dg_packaged_data.size()}; + neighbor_package_data[mortar_id] = DataVector{dg_packaged_data.size()}; + std::copy(dg_packaged_data_view.begin(), dg_packaged_data_view.end(), + neighbor_package_data[mortar_id].begin()); + } + } + }); + + return neighbor_package_data; +} +} // namespace grmhd::ValenciaDivClean::subcell + +template void evolution::dg::subcell::neighbor_reconstructed_face_solution< + 3, grmhd::ValenciaDivClean::subcell::NeighborPackagedData>( + gsl::not_null box, + gsl::not_null, Mesh<2>, std::optional, + std::optional, ::TimeStepId, int>>>*> + received_temporal_id_and_data); diff --git a/src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/NeighborPackagedData.hpp b/src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/NeighborPackagedData.hpp index c673c790f076..0ba12da448d8 100644 --- a/src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/NeighborPackagedData.hpp +++ b/src/Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/NeighborPackagedData.hpp @@ -4,50 +4,18 @@ #pragma once #include -#include -#include -#include #include -#include "DataStructures/DataBox/DataBox.hpp" -#include "DataStructures/DataBox/PrefixHelpers.hpp" -#include "DataStructures/DataBox/Prefixes.hpp" -#include "DataStructures/DataVector.hpp" -#include "DataStructures/Index.hpp" -#include "DataStructures/Tensor/Tensor.hpp" -#include "DataStructures/Variables.hpp" -#include "DataStructures/VariablesTag.hpp" -#include "Domain/Structure/Direction.hpp" -#include "Domain/Structure/DirectionalIdMap.hpp" -#include "Domain/Structure/Element.hpp" -#include "Domain/Structure/ElementId.hpp" -#include "Domain/Tags.hpp" -#include "Domain/TagsTimeDependent.hpp" -#include "Evolution/BoundaryCorrectionTags.hpp" -#include "Evolution/DgSubcell/Projection.hpp" -#include "Evolution/DgSubcell/Reconstruction.hpp" -#include "Evolution/DgSubcell/ReconstructionMethod.hpp" -#include "Evolution/DgSubcell/SubcellOptions.hpp" -#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" -#include "Evolution/DgSubcell/Tags/Mesh.hpp" -#include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp" -#include "Evolution/DgSubcell/Tags/SubcellOptions.hpp" -#include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp" -#include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp" -#include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/BoundaryCorrection.hpp" -#include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/Factory.hpp" -#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.hpp" -#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp" -#include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Tag.hpp" -#include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/ComputeFluxes.hpp" -#include "Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp" -#include "NumericalAlgorithms/Spectral/Mesh.hpp" -#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" -#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp" -#include "PointwiseFunctions/Hydro/Tags.hpp" -#include "Utilities/CallWithDynamicType.hpp" -#include "Utilities/Gsl.hpp" -#include "Utilities/TMPL.hpp" +/// \cond +namespace db { +class Access; +} // namespace db +class DataVector; +template +struct DirectionalId; +template +class DirectionalIdMap; +/// \endcond namespace grmhd::ValenciaDivClean::subcell { /*! @@ -65,231 +33,8 @@ namespace grmhd::ValenciaDivClean::subcell { * if one element is doing DG, then we aren't at a shock. */ struct NeighborPackagedData { - template static DirectionalIdMap<3, DataVector> apply( - const db::DataBox& box, - const std::vector>& mortars_to_reconstruct_to) { - using evolved_vars_tag = typename System::variables_tag; - using evolved_vars_tags = typename evolved_vars_tag::tags_list; - using fluxes_tags = db::wrap_tags_in<::Tags::Flux, evolved_vars_tags, - tmpl::size_t<3>, Frame::Inertial>; - - ASSERT(not db::get>(box).has_value(), - "Haven't yet added support for moving mesh to DG-subcell. This " - "should be easy to generalize, but we will want to consider " - "storing the mesh velocity on the faces instead of " - "re-slicing/projecting."); - - DirectionalIdMap<3, DataVector> neighbor_package_data{}; - if (mortars_to_reconstruct_to.empty()) { - return neighbor_package_data; - } - - const auto& ghost_subcell_data = - db::get>( - box); - const Mesh<3>& subcell_mesh = - db::get>(box); - const Mesh<3>& dg_mesh = db::get>(box); - const auto& subcell_options = - db::get>(box); - - // Note: we need to compare if projecting the entire mesh or only ghost - // zones needed is faster. This probably depends on the number of neighbors - // we have doing FD. - const auto volume_prims = evolution::dg::subcell::fd::project( - db::get(box), dg_mesh, - subcell_mesh.extents()); - - const auto& recons = - db::get(box); - const auto& boundary_correction = - db::get>(box); - using derived_boundary_corrections = - typename std::decay_t::creatable_classes; - tmpl::for_each< - derived_boundary_corrections>([&box, &boundary_correction, &dg_mesh, - &mortars_to_reconstruct_to, - &neighbor_package_data, - &ghost_subcell_data, &recons, - &subcell_mesh, &subcell_options, - &volume_prims]( - auto derived_correction_v) { - using DerivedCorrection = tmpl::type_from; - if (typeid(boundary_correction) == typeid(DerivedCorrection)) { - using dg_package_data_temporary_tags = - typename DerivedCorrection::dg_package_data_temporary_tags; - using dg_package_data_argument_tags = fd::tags_list_for_reconstruct; - - const auto& element = db::get>(box); - const auto& eos = get(box); - - using dg_package_field_tags = - typename DerivedCorrection::dg_package_field_tags; - Variables vars_on_face{0}; - Variables packaged_data{0}; - for (const auto& mortar_id : mortars_to_reconstruct_to) { - const Direction<3>& direction = mortar_id.direction; - - const Mesh<2> dg_face_mesh = - dg_mesh.slice_away(direction.dimension()); - Index<3> extents = subcell_mesh.extents(); - // Switch to face-centered instead of cell-centered points on the FD. - // There are num_cell_centered+1 face-centered points. - ++extents[direction.dimension()]; - const Index<2> subcell_face_extents = - extents.slice_away(direction.dimension()); - - // Computed prims and cons on face via reconstruction - const size_t num_face_pts = subcell_mesh.extents() - .slice_away(direction.dimension()) - .product(); - vars_on_face.initialize(num_face_pts); - // Copy spacetime vars over from volume. - using spacetime_vars_to_copy = - tmpl::list, - gr::Tags::Shift, - gr::Tags::SpatialMetric, - gr::Tags::SqrtDetSpatialMetric, - gr::Tags::InverseSpatialMetric>; - tmpl::for_each( - [&direction, &extents, &vars_on_face, - &spacetime_vars_on_faces = - db::get>(box)]( - auto tag_v) { - using tag = tmpl::type_from; - data_on_slice(make_not_null(&get(vars_on_face)), - get(gsl::at(spacetime_vars_on_faces, - direction.dimension())), - extents, direction.dimension(), - direction.side() == Side::Lower - ? 0 - : extents[direction.dimension()] - 1); - }); - - call_with_dynamic_type( - &recons, - [&element, &eos, &mortar_id, &ghost_subcell_data, &subcell_mesh, - &vars_on_face, &volume_prims](const auto& reconstructor) { - reconstructor->reconstruct_fd_neighbor( - make_not_null(&vars_on_face), volume_prims, eos, element, - ghost_subcell_data, subcell_mesh, mortar_id.direction); - }); - - // Get the mesh velocity if needed - const std::optional> - mesh_velocity_dg = db::get>(box); - std::optional> - mesh_velocity_on_subcell_face = {}; - if (mesh_velocity_dg.has_value()) { - // Slice data on current face - tnsr::I mesh_velocity_on_dg_face = - data_on_slice( - mesh_velocity_dg.value(), dg_mesh.extents(), - direction.dimension(), - direction.side() == Side::Lower - ? 0 - : (dg_mesh.extents(direction.dimension()) - 1)); - - mesh_velocity_on_subcell_face = - tnsr::I{num_face_pts}; - - for (size_t i = 0; i < 3; i++) { - evolution::dg::subcell::fd::project( - make_not_null(&mesh_velocity_on_subcell_face.value().get(i)), - mesh_velocity_on_dg_face.get(i), dg_face_mesh, - subcell_face_extents); - } - } - - grmhd::ValenciaDivClean::subcell::compute_fluxes( - make_not_null(&vars_on_face)); - - if (mesh_velocity_on_subcell_face.has_value()) { - tmpl::for_each([&vars_on_face, - &mesh_velocity_on_subcell_face]( - auto tag_v) { - using tag = tmpl::type_from; - using flux_tag = - ::Tags::Flux, Frame::Inertial>; - using FluxTensor = typename flux_tag::type; - const auto& var = get(vars_on_face); - auto& flux = get(vars_on_face); - for (size_t storage_index = 0; storage_index < var.size(); - ++storage_index) { - const auto tensor_index = var.get_tensor_index(storage_index); - for (size_t j = 0; j < 3; j++) { - const auto flux_storage_index = - FluxTensor::get_storage_index(prepend(tensor_index, j)); - flux[flux_storage_index] -= - mesh_velocity_on_subcell_face.value().get(j) * - var[storage_index]; - } - } - }); - } - - // Note: since the spacetime isn't dynamical we don't need to - // worry about different normal vectors on the different sides - // of the element. If the spacetime is dynamical, then the normal - // covector needs to be sent, or at least the inverse spatial metric - // from the neighbor so that the normal covector can be computed - // correctly. - tnsr::i normal_covector = - get>( - *db::get>( - box) - .at(mortar_id.direction)); - for (auto& t : normal_covector) { - t *= -1.0; - } - // Note: Only need to do the projection in 2d and 3d, but GRMHD is - // always 3d currently. - const auto dg_normal_covector = normal_covector; - for (size_t i = 0; i < 3; ++i) { - normal_covector.get(i) = evolution::dg::subcell::fd::project( - dg_normal_covector.get(i), - dg_mesh.slice_away(mortar_id.direction.dimension()), - subcell_mesh.extents().slice_away( - mortar_id.direction.dimension())); - } - - // Compute the packaged data - packaged_data.initialize(num_face_pts); - using dg_package_data_projected_tags = tmpl::append< - evolved_vars_tags, fluxes_tags, dg_package_data_temporary_tags, - typename DerivedCorrection::dg_package_data_primitive_tags>; - evolution::dg::Actions::detail::dg_package_data( - make_not_null(&packaged_data), - dynamic_cast(boundary_correction), - vars_on_face, normal_covector, mesh_velocity_on_subcell_face, box, - typename DerivedCorrection::dg_package_data_volume_tags{}, - dg_package_data_projected_tags{}); - - // Reconstruct the DG solution. - // Really we should be solving the boundary correction and - // then reconstructing, but away from a shock this doesn't - // matter. - auto dg_packaged_data = evolution::dg::subcell::fd::reconstruct( - packaged_data, - dg_mesh.slice_away(mortar_id.direction.dimension()), - subcell_mesh.extents().slice_away( - mortar_id.direction.dimension()), - subcell_options.reconstruction_method()); - // Make a view so we can use iterators with std::copy - DataVector dg_packaged_data_view{dg_packaged_data.data(), - dg_packaged_data.size()}; - neighbor_package_data[mortar_id] = - DataVector{dg_packaged_data.size()}; - std::copy(dg_packaged_data_view.begin(), dg_packaged_data_view.end(), - neighbor_package_data[mortar_id].begin()); - } - } - }); - - return neighbor_package_data; - } + const db::Access& box, + const std::vector>& mortars_to_reconstruct_to); }; } // namespace grmhd::ValenciaDivClean::subcell diff --git a/src/Evolution/Systems/NewtonianEuler/BoundaryCorrections/Hll.hpp b/src/Evolution/Systems/NewtonianEuler/BoundaryCorrections/Hll.hpp index 2f4c97dc5dfb..b494d5a4cdf8 100644 --- a/src/Evolution/Systems/NewtonianEuler/BoundaryCorrections/Hll.hpp +++ b/src/Evolution/Systems/NewtonianEuler/BoundaryCorrections/Hll.hpp @@ -123,6 +123,7 @@ class Hll final : public BoundaryCorrection { hydro::Tags::SpecificInternalEnergy>; using dg_package_data_volume_tags = tmpl::list>; + using dg_boundary_terms_volume_tags = tmpl::list<>; double dg_package_data( gsl::not_null*> packaged_mass_density, diff --git a/src/Evolution/Systems/NewtonianEuler/BoundaryCorrections/Hllc.hpp b/src/Evolution/Systems/NewtonianEuler/BoundaryCorrections/Hllc.hpp index 9ff1d26c285b..076f50be429e 100644 --- a/src/Evolution/Systems/NewtonianEuler/BoundaryCorrections/Hllc.hpp +++ b/src/Evolution/Systems/NewtonianEuler/BoundaryCorrections/Hllc.hpp @@ -219,6 +219,7 @@ class Hllc final : public BoundaryCorrection { hydro::Tags::SpecificInternalEnergy>; using dg_package_data_volume_tags = tmpl::list>; + using dg_boundary_terms_volume_tags = tmpl::list<>; double dg_package_data( gsl::not_null*> packaged_mass_density, diff --git a/src/Evolution/Systems/NewtonianEuler/BoundaryCorrections/Rusanov.hpp b/src/Evolution/Systems/NewtonianEuler/BoundaryCorrections/Rusanov.hpp index 6b92d3b62cff..5876aa87b0bf 100644 --- a/src/Evolution/Systems/NewtonianEuler/BoundaryCorrections/Rusanov.hpp +++ b/src/Evolution/Systems/NewtonianEuler/BoundaryCorrections/Rusanov.hpp @@ -102,6 +102,7 @@ class Rusanov final : public BoundaryCorrection { hydro::Tags::SpecificInternalEnergy>; using dg_package_data_volume_tags = tmpl::list>; + using dg_boundary_terms_volume_tags = tmpl::list<>; double dg_package_data( gsl::not_null*> packaged_mass_density, diff --git a/src/Evolution/Systems/NewtonianEuler/Subcell/CMakeLists.txt b/src/Evolution/Systems/NewtonianEuler/Subcell/CMakeLists.txt index c69bd5cc4641..3684305a31c6 100644 --- a/src/Evolution/Systems/NewtonianEuler/Subcell/CMakeLists.txt +++ b/src/Evolution/Systems/NewtonianEuler/Subcell/CMakeLists.txt @@ -4,6 +4,7 @@ spectre_target_sources( ${LIBRARY} PRIVATE + NeighborPackagedData.cpp PrimitiveGhostData.cpp PrimsAfterRollback.cpp ResizeAndComputePrimitives.cpp diff --git a/src/Evolution/Systems/NewtonianEuler/Subcell/NeighborPackagedData.cpp b/src/Evolution/Systems/NewtonianEuler/Subcell/NeighborPackagedData.cpp new file mode 100644 index 000000000000..c804ff46323b --- /dev/null +++ b/src/Evolution/Systems/NewtonianEuler/Subcell/NeighborPackagedData.cpp @@ -0,0 +1,243 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/Systems/NewtonianEuler/Subcell/NeighborPackagedData.hpp" + +#include +#include +#include +#include +#include + +#include "DataStructures/DataBox/Access.hpp" +#include "DataStructures/DataBox/PrefixHelpers.hpp" +#include "DataStructures/DataBox/Prefixes.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Index.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "DataStructures/VariablesTag.hpp" +#include "Domain/Structure/Direction.hpp" +#include "Domain/Structure/DirectionalIdMap.hpp" +#include "Domain/Structure/Element.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Domain/Tags.hpp" +#include "Domain/TagsTimeDependent.hpp" +#include "Evolution/BoundaryCorrectionTags.hpp" +#include "Evolution/DgSubcell/NeighborReconstructedFaceSolution.tpp" +#include "Evolution/DgSubcell/Projection.hpp" +#include "Evolution/DgSubcell/Reconstruction.hpp" +#include "Evolution/DgSubcell/ReconstructionMethod.hpp" +#include "Evolution/DgSubcell/SubcellOptions.hpp" +#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" +#include "Evolution/DgSubcell/Tags/Mesh.hpp" +#include "Evolution/DgSubcell/Tags/SubcellOptions.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp" +#include "Evolution/Systems/NewtonianEuler/BoundaryCorrections/BoundaryCorrection.hpp" +#include "Evolution/Systems/NewtonianEuler/BoundaryCorrections/Factory.hpp" +#include "Evolution/Systems/NewtonianEuler/FiniteDifference/Factory.hpp" +#include "Evolution/Systems/NewtonianEuler/FiniteDifference/Reconstructor.hpp" +#include "Evolution/Systems/NewtonianEuler/FiniteDifference/Tag.hpp" +#include "Evolution/Systems/NewtonianEuler/Subcell/ComputeFluxes.hpp" +#include "Evolution/Systems/NewtonianEuler/System.hpp" +#include "Evolution/Systems/NewtonianEuler/Tags.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "Parallel/Tags/Metavariables.hpp" +#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp" +#include "PointwiseFunctions/Hydro/Tags.hpp" +#include "Utilities/CallWithDynamicType.hpp" +#include "Utilities/GenerateInstantiations.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/TMPL.hpp" + +namespace NewtonianEuler::subcell { +template +DirectionalIdMap NeighborPackagedData::apply( + const db::Access& box, + const std::vector>& mortars_to_reconstruct_to) { + using system = NewtonianEuler::System; + using evolved_vars_tag = typename system::variables_tag; + using evolved_vars_tags = typename evolved_vars_tag::tags_list; + using prim_tags = typename system::primitive_variables_tag::tags_list; + using fluxes_tags = db::wrap_tags_in<::Tags::Flux, evolved_vars_tags, + tmpl::size_t, Frame::Inertial>; + + ASSERT(not db::get>(box).has_value(), + "Haven't yet added support for moving mesh to DG-subcell. This " + "should be easy to generalize, but we will want to consider " + "storing the mesh velocity on the faces instead of " + "re-slicing/projecting."); + + DirectionalIdMap neighbor_package_data{}; + if (mortars_to_reconstruct_to.empty()) { + return neighbor_package_data; + } + + const auto& ghost_subcell_data = + db::get>( + box); + const Mesh& subcell_mesh = + db::get>(box); + const Mesh& dg_mesh = db::get>(box); + const auto& subcell_options = + db::get>(box); + + // Note: we need to compare if projecting the entire mesh or only ghost + // zones needed is faster. This probably depends on the number of neighbors + // we have doing FD. + const auto volume_prims = evolution::dg::subcell::fd::project( + db::get(box), dg_mesh, + subcell_mesh.extents()); + + const auto& recons = + db::get>(box); + const auto& boundary_correction = + db::get>(box); + using derived_boundary_corrections = + typename std::decay_t::creatable_classes; + tmpl::for_each([&box, &boundary_correction, + &dg_mesh, + &mortars_to_reconstruct_to, + &neighbor_package_data, + &ghost_subcell_data, &recons, + &subcell_mesh, &subcell_options, + &volume_prims]( + auto derived_correction_v) { + using DerivedCorrection = tmpl::type_from; + if (typeid(boundary_correction) == typeid(DerivedCorrection)) { + using dg_package_data_temporary_tags = + typename DerivedCorrection::dg_package_data_temporary_tags; + using dg_package_data_argument_tags = + tmpl::append; + + const auto& element = db::get>(box); + const auto& eos = get>(box); + + using dg_package_field_tags = + typename DerivedCorrection::dg_package_field_tags; + Variables vars_on_face; + Variables packaged_data; + for (const auto& mortar_id : mortars_to_reconstruct_to) { + const Direction& direction = mortar_id.direction; + + Index extents = subcell_mesh.extents(); + // Switch to face-centered instead of cell-centered points on the FD. + // There are num_cell_centered+1 face-centered points. + ++extents[direction.dimension()]; + + // Computed prims and cons on face via reconstruction + const size_t num_face_pts = + subcell_mesh.extents().slice_away(direction.dimension()).product(); + vars_on_face.initialize(num_face_pts); + + call_with_dynamic_type::creatable_classes>( + &recons, + [&element, &eos, &mortar_id, &ghost_subcell_data, &subcell_mesh, + &vars_on_face, &volume_prims](const auto& reconstructor) { + reconstructor->reconstruct_fd_neighbor( + make_not_null(&vars_on_face), volume_prims, eos, element, + ghost_subcell_data, subcell_mesh, mortar_id.direction); + }); + + NewtonianEuler::subcell::compute_fluxes( + make_not_null(&vars_on_face)); + + tnsr::i normal_covector = + get>( + *db::get>( + box) + .at(mortar_id.direction)); + for (auto& t : normal_covector) { + t *= -1.0; + } + if constexpr (Dim > 1) { + const auto dg_normal_covector = normal_covector; + for (size_t i = 0; i < Dim; ++i) { + normal_covector.get(i) = evolution::dg::subcell::fd::project( + dg_normal_covector.get(i), + dg_mesh.slice_away(mortar_id.direction.dimension()), + subcell_mesh.extents().slice_away( + mortar_id.direction.dimension())); + } + } + + // Compute the packaged data + packaged_data.initialize(num_face_pts); + using dg_package_data_projected_tags = tmpl::append< + evolved_vars_tags, fluxes_tags, dg_package_data_temporary_tags, + typename DerivedCorrection::dg_package_data_primitive_tags>; + evolution::dg::Actions::detail::dg_package_data( + make_not_null(&packaged_data), + dynamic_cast(boundary_correction), + vars_on_face, normal_covector, {std::nullopt}, box, + typename DerivedCorrection::dg_package_data_volume_tags{}, + dg_package_data_projected_tags{}); + + if constexpr (Dim == 1) { + (void)dg_mesh; + (void)subcell_options; + // Make a view so we can use iterators with std::copy + DataVector packaged_data_view{packaged_data.data(), + packaged_data.size()}; + neighbor_package_data[mortar_id] = DataVector{packaged_data.size()}; + std::copy(packaged_data_view.begin(), packaged_data_view.end(), + neighbor_package_data[mortar_id].begin()); + } else { + // Reconstruct the DG solution. + // Really we should be solving the boundary correction and + // then reconstructing, but away from a shock this doesn't + // matter. + auto dg_packaged_data = evolution::dg::subcell::fd::reconstruct( + packaged_data, + dg_mesh.slice_away(mortar_id.direction.dimension()), + subcell_mesh.extents().slice_away( + mortar_id.direction.dimension()), + subcell_options.reconstruction_method()); + // Make a view so we can use iterators with std::copy + DataVector dg_packaged_data_view{dg_packaged_data.data(), + dg_packaged_data.size()}; + neighbor_package_data[mortar_id] = + DataVector{dg_packaged_data.size()}; + std::copy(dg_packaged_data_view.begin(), dg_packaged_data_view.end(), + neighbor_package_data[mortar_id].begin()); + } + } + } + }); + + return neighbor_package_data; +} + +#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) + +#define INSTANTIATION(r, data) \ + template DirectionalIdMap \ + NeighborPackagedData::apply( \ + const db::Access& box, \ + const std::vector>& mortars_to_reconstruct_to); + +GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) + +#undef INSTANTIATION +} // namespace NewtonianEuler::subcell + +#define INSTANTIATION(r, data) \ + template void evolution::dg::subcell::neighbor_reconstructed_face_solution< \ + DIM(data), NewtonianEuler::subcell::NeighborPackagedData>( \ + gsl::not_null box, \ + gsl::not_null, Mesh, \ + std::optional, std::optional, \ + ::TimeStepId, int>>>*> \ + received_temporal_id_and_data); + +GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) + +#undef INSTANTIATION +#undef DIM diff --git a/src/Evolution/Systems/NewtonianEuler/Subcell/NeighborPackagedData.hpp b/src/Evolution/Systems/NewtonianEuler/Subcell/NeighborPackagedData.hpp index 73e5885729c8..e863be7df60f 100644 --- a/src/Evolution/Systems/NewtonianEuler/Subcell/NeighborPackagedData.hpp +++ b/src/Evolution/Systems/NewtonianEuler/Subcell/NeighborPackagedData.hpp @@ -4,49 +4,18 @@ #pragma once #include -#include -#include -#include #include -#include "DataStructures/DataBox/DataBox.hpp" -#include "DataStructures/DataBox/PrefixHelpers.hpp" -#include "DataStructures/DataBox/Prefixes.hpp" -#include "DataStructures/DataVector.hpp" -#include "DataStructures/Index.hpp" -#include "DataStructures/Tensor/Tensor.hpp" -#include "DataStructures/Variables.hpp" -#include "DataStructures/VariablesTag.hpp" -#include "Domain/Structure/Direction.hpp" -#include "Domain/Structure/DirectionalIdMap.hpp" -#include "Domain/Structure/Element.hpp" -#include "Domain/Structure/ElementId.hpp" -#include "Domain/Tags.hpp" -#include "Domain/TagsTimeDependent.hpp" -#include "Evolution/BoundaryCorrectionTags.hpp" -#include "Evolution/DgSubcell/Projection.hpp" -#include "Evolution/DgSubcell/Reconstruction.hpp" -#include "Evolution/DgSubcell/ReconstructionMethod.hpp" -#include "Evolution/DgSubcell/SubcellOptions.hpp" -#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" -#include "Evolution/DgSubcell/Tags/Mesh.hpp" -#include "Evolution/DgSubcell/Tags/SubcellOptions.hpp" -#include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp" -#include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp" -#include "Evolution/Systems/NewtonianEuler/BoundaryCorrections/BoundaryCorrection.hpp" -#include "Evolution/Systems/NewtonianEuler/BoundaryCorrections/Factory.hpp" -#include "Evolution/Systems/NewtonianEuler/FiniteDifference/Reconstructor.hpp" -#include "Evolution/Systems/NewtonianEuler/FiniteDifference/Tag.hpp" -#include "Evolution/Systems/NewtonianEuler/Subcell/ComputeFluxes.hpp" -#include "Evolution/Systems/NewtonianEuler/System.hpp" -#include "Evolution/Systems/NewtonianEuler/Tags.hpp" -#include "NumericalAlgorithms/Spectral/Mesh.hpp" -#include "Parallel/Tags/Metavariables.hpp" -#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp" -#include "PointwiseFunctions/Hydro/Tags.hpp" -#include "Utilities/CallWithDynamicType.hpp" -#include "Utilities/Gsl.hpp" -#include "Utilities/TMPL.hpp" +/// \cond +namespace db { +class Access; +} // namespace db +class DataVector; +template +struct DirectionalId; +template +class DirectionalIdMap; +/// \endcond namespace NewtonianEuler::subcell { /*! @@ -64,167 +33,9 @@ namespace NewtonianEuler::subcell { * if one element is doing DG, then we aren't at a shock. */ struct NeighborPackagedData { - template + template static DirectionalIdMap apply( - const db::DataBox& box, - const std::vector>& mortars_to_reconstruct_to) { - using system = typename std::decay_t(box))>::system; - using evolved_vars_tag = typename system::variables_tag; - using evolved_vars_tags = typename evolved_vars_tag::tags_list; - using prim_tags = typename system::primitive_variables_tag::tags_list; - using fluxes_tags = db::wrap_tags_in<::Tags::Flux, evolved_vars_tags, - tmpl::size_t, Frame::Inertial>; - - ASSERT(not db::get>(box).has_value(), - "Haven't yet added support for moving mesh to DG-subcell. This " - "should be easy to generalize, but we will want to consider " - "storing the mesh velocity on the faces instead of " - "re-slicing/projecting."); - - DirectionalIdMap neighbor_package_data{}; - if (mortars_to_reconstruct_to.empty()) { - return neighbor_package_data; - } - - const auto& ghost_subcell_data = - db::get>( - box); - const Mesh& subcell_mesh = - db::get>(box); - const Mesh& dg_mesh = db::get>(box); - const auto& subcell_options = - db::get>(box); - - // Note: we need to compare if projecting the entire mesh or only ghost - // zones needed is faster. This probably depends on the number of neighbors - // we have doing FD. - const auto volume_prims = evolution::dg::subcell::fd::project( - db::get(box), dg_mesh, - subcell_mesh.extents()); - - const auto& recons = - db::get>(box); - const auto& boundary_correction = - db::get>(box); - using derived_boundary_corrections = - typename std::decay_t::creatable_classes; - tmpl::for_each< - derived_boundary_corrections>([&box, &boundary_correction, &dg_mesh, - &mortars_to_reconstruct_to, - &neighbor_package_data, - &ghost_subcell_data, &recons, - &subcell_mesh, &subcell_options, - &volume_prims]( - auto derived_correction_v) { - using DerivedCorrection = tmpl::type_from; - if (typeid(boundary_correction) == typeid(DerivedCorrection)) { - using dg_package_data_temporary_tags = - typename DerivedCorrection::dg_package_data_temporary_tags; - using dg_package_data_argument_tags = - tmpl::append; - - const auto& element = db::get>(box); - const auto& eos = get>(box); - - using dg_package_field_tags = - typename DerivedCorrection::dg_package_field_tags; - Variables vars_on_face; - Variables packaged_data; - for (const auto& mortar_id : mortars_to_reconstruct_to) { - const Direction& direction = mortar_id.direction; - - Index extents = subcell_mesh.extents(); - // Switch to face-centered instead of cell-centered points on the FD. - // There are num_cell_centered+1 face-centered points. - ++extents[direction.dimension()]; - - // Computed prims and cons on face via reconstruction - const size_t num_face_pts = subcell_mesh.extents() - .slice_away(direction.dimension()) - .product(); - vars_on_face.initialize(num_face_pts); - - call_with_dynamic_type::creatable_classes>( - &recons, - [&element, &eos, &mortar_id, &ghost_subcell_data, &subcell_mesh, - &vars_on_face, &volume_prims](const auto& reconstructor) { - reconstructor->reconstruct_fd_neighbor( - make_not_null(&vars_on_face), volume_prims, eos, element, - ghost_subcell_data, subcell_mesh, mortar_id.direction); - }); - - NewtonianEuler::subcell::compute_fluxes( - make_not_null(&vars_on_face)); - - tnsr::i normal_covector = get< - evolution::dg::Tags::NormalCovector>( - *db::get>( - box) - .at(mortar_id.direction)); - for (auto& t : normal_covector) { - t *= -1.0; - } - if constexpr (Dim > 1) { - const auto dg_normal_covector = normal_covector; - for (size_t i = 0; i < Dim; ++i) { - normal_covector.get(i) = evolution::dg::subcell::fd::project( - dg_normal_covector.get(i), - dg_mesh.slice_away(mortar_id.direction.dimension()), - subcell_mesh.extents().slice_away( - mortar_id.direction.dimension())); - } - } - - // Compute the packaged data - packaged_data.initialize(num_face_pts); - using dg_package_data_projected_tags = tmpl::append< - evolved_vars_tags, fluxes_tags, dg_package_data_temporary_tags, - typename DerivedCorrection::dg_package_data_primitive_tags>; - evolution::dg::Actions::detail::dg_package_data( - make_not_null(&packaged_data), - dynamic_cast(boundary_correction), - vars_on_face, normal_covector, {std::nullopt}, box, - typename DerivedCorrection::dg_package_data_volume_tags{}, - dg_package_data_projected_tags{}); - - if constexpr (Dim == 1) { - (void)dg_mesh; - (void)subcell_options; - // Make a view so we can use iterators with std::copy - DataVector packaged_data_view{packaged_data.data(), - packaged_data.size()}; - neighbor_package_data[mortar_id] = DataVector{packaged_data.size()}; - std::copy(packaged_data_view.begin(), packaged_data_view.end(), - neighbor_package_data[mortar_id].begin()); - } else { - // Reconstruct the DG solution. - // Really we should be solving the boundary correction and - // then reconstructing, but away from a shock this doesn't - // matter. - auto dg_packaged_data = evolution::dg::subcell::fd::reconstruct( - packaged_data, - dg_mesh.slice_away(mortar_id.direction.dimension()), - subcell_mesh.extents().slice_away( - mortar_id.direction.dimension()), - subcell_options.reconstruction_method()); - // Make a view so we can use iterators with std::copy - DataVector dg_packaged_data_view{dg_packaged_data.data(), - dg_packaged_data.size()}; - neighbor_package_data[mortar_id] = - DataVector{dg_packaged_data.size()}; - std::copy(dg_packaged_data_view.begin(), - dg_packaged_data_view.end(), - neighbor_package_data[mortar_id].begin()); - } - } - } - }); - - return neighbor_package_data; - } + const db::Access& box, + const std::vector>& mortars_to_reconstruct_to); }; } // namespace NewtonianEuler::subcell diff --git a/src/Evolution/Systems/RadiationTransport/M1Grey/BoundaryCorrections/Rusanov.hpp b/src/Evolution/Systems/RadiationTransport/M1Grey/BoundaryCorrections/Rusanov.hpp index ce895766107f..e45a25eba087 100644 --- a/src/Evolution/Systems/RadiationTransport/M1Grey/BoundaryCorrections/Rusanov.hpp +++ b/src/Evolution/Systems/RadiationTransport/M1Grey/BoundaryCorrections/Rusanov.hpp @@ -133,6 +133,7 @@ class Rusanov> final using dg_package_data_temporary_tags = tmpl::list<>; using dg_package_data_primitive_tags = tmpl::list<>; using dg_package_data_volume_tags = tmpl::list<>; + using dg_boundary_terms_volume_tags = tmpl::list<>; double dg_package_data( const gsl::not_null { tmpl::list, AbsCharSpeed>; using dg_package_data_temporary_tags = tmpl::list>; using dg_package_data_volume_tags = tmpl::list<>; + using dg_boundary_terms_volume_tags = tmpl::list<>; static double dg_package_data( gsl::not_null*> packaged_u, diff --git a/src/Evolution/Systems/ScalarAdvection/FiniteDifference/AoWeno.cpp b/src/Evolution/Systems/ScalarAdvection/FiniteDifference/AoWeno.cpp index 54a1976a1be3..6a909b534837 100644 --- a/src/Evolution/Systems/ScalarAdvection/FiniteDifference/AoWeno.cpp +++ b/src/Evolution/Systems/ScalarAdvection/FiniteDifference/AoWeno.cpp @@ -78,7 +78,7 @@ void AoWeno53::reconstruct( vars_on_lower_face, const gsl::not_null, Dim>*> vars_on_upper_face, - const Variables>& volume_vars, + const Variables>& volume_vars, const Element& element, const DirectionalIdMap& ghost_data, const Mesh& subcell_mesh) const { @@ -98,7 +98,7 @@ template template void AoWeno53::reconstruct_fd_neighbor( const gsl::not_null*> vars_on_face, - const Variables>& volume_vars, + const Variables>& volume_vars, const Element& element, const DirectionalIdMap& ghost_data, const Mesh& subcell_mesh, @@ -144,10 +144,11 @@ bool operator!=(const AoWeno53& lhs, const AoWeno53& rhs) { } #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) -#define TAGS_LIST(data) \ - tmpl::list, Frame::Inertial>, \ - Tags::VelocityField> +#define TAGS_LIST(data) \ + tmpl::list, \ + Frame::Inertial>, \ + ScalarAdvection::Tags::VelocityField> #define INSTANTIATION(r, data) \ template class AoWeno53; \ @@ -165,14 +166,14 @@ GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2)) vars_on_lower_face, \ gsl::not_null, DIM(data)>*> \ vars_on_upper_face, \ - const Variables>& volume_vars, \ + const Variables>& volume_vars, \ const Element& element, \ const DirectionalIdMap& \ ghost_data, \ const Mesh& subcell_mesh) const; \ template void AoWeno53::reconstruct_fd_neighbor( \ gsl::not_null*> vars_on_face, \ - const Variables>& volume_vars, \ + const Variables>& volume_vars, \ const Element& element, \ const DirectionalIdMap& \ ghost_data, \ diff --git a/src/Evolution/Systems/ScalarAdvection/FiniteDifference/AoWeno.hpp b/src/Evolution/Systems/ScalarAdvection/FiniteDifference/AoWeno.hpp index 1bbc5baf4068..409c47d0fc71 100644 --- a/src/Evolution/Systems/ScalarAdvection/FiniteDifference/AoWeno.hpp +++ b/src/Evolution/Systems/ScalarAdvection/FiniteDifference/AoWeno.hpp @@ -57,9 +57,10 @@ template class AoWeno53 : public Reconstructor { private: using face_vars_tags = - tmpl::list, Frame::Inertial>>; - using volume_vars_tags = tmpl::list; + tmpl::list, + Frame::Inertial>>; + using volume_vars_tags = tmpl::list; public: struct GammaHi { @@ -119,7 +120,7 @@ class AoWeno53 : public Reconstructor { void reconstruct( gsl::not_null, Dim>*> vars_on_lower_face, gsl::not_null, Dim>*> vars_on_upper_face, - const Variables>& volume_vars, + const Variables>& volume_vars, const Element& element, const DirectionalIdMap& ghost_data, @@ -128,7 +129,7 @@ class AoWeno53 : public Reconstructor { template void reconstruct_fd_neighbor( gsl::not_null*> vars_on_face, - const Variables>& volume_vars, + const Variables>& volume_vars, const Element& element, const DirectionalIdMap& ghost_data, diff --git a/src/Evolution/Systems/ScalarAdvection/FiniteDifference/MonotonisedCentral.cpp b/src/Evolution/Systems/ScalarAdvection/FiniteDifference/MonotonisedCentral.cpp index c851e2083bd0..4dd9654d1d63 100644 --- a/src/Evolution/Systems/ScalarAdvection/FiniteDifference/MonotonisedCentral.cpp +++ b/src/Evolution/Systems/ScalarAdvection/FiniteDifference/MonotonisedCentral.cpp @@ -53,7 +53,7 @@ void MonotonisedCentral::reconstruct( vars_on_lower_face, const gsl::not_null, Dim>*> vars_on_upper_face, - const Variables>& volume_vars, + const Variables>& volume_vars, const Element& element, const DirectionalIdMap& ghost_data, const Mesh& subcell_mesh) const { @@ -73,7 +73,7 @@ template template void MonotonisedCentral::reconstruct_fd_neighbor( const gsl::not_null*> vars_on_face, - const Variables>& volume_vars, + const Variables>& volume_vars, const Element& element, const DirectionalIdMap& ghost_data, const Mesh& subcell_mesh, @@ -121,9 +121,10 @@ bool operator!=(const MonotonisedCentral& lhs, } #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) -#define TAGS_LIST(data) \ - tmpl::list, Frame::Inertial>, \ +#define TAGS_LIST(data) \ + tmpl::list, \ + Frame::Inertial>, \ Tags::VelocityField> #define INSTANTIATION(r, data) \ @@ -144,14 +145,14 @@ GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2)) vars_on_lower_face, \ gsl::not_null, DIM(data)>*> \ vars_on_upper_face, \ - const Variables>& volume_vars, \ + const Variables>& volume_vars, \ const Element& element, \ const DirectionalIdMap& \ ghost_data, \ const Mesh& subcell_mesh) const; \ template void MonotonisedCentral::reconstruct_fd_neighbor( \ gsl::not_null*> vars_on_face, \ - const Variables>& volume_vars, \ + const Variables>& volume_vars, \ const Element& element, \ const DirectionalIdMap& \ ghost_data, \ diff --git a/src/Evolution/Systems/ScalarAdvection/FiniteDifference/MonotonisedCentral.hpp b/src/Evolution/Systems/ScalarAdvection/FiniteDifference/MonotonisedCentral.hpp index 47e9231c8430..040ff7c63b47 100644 --- a/src/Evolution/Systems/ScalarAdvection/FiniteDifference/MonotonisedCentral.hpp +++ b/src/Evolution/Systems/ScalarAdvection/FiniteDifference/MonotonisedCentral.hpp @@ -53,7 +53,7 @@ namespace ScalarAdvection::fd { template class MonotonisedCentral : public Reconstructor { private: - using volume_vars_tags = tmpl::list; + using volume_vars_tags = tmpl::list; public: using options = tmpl::list<>; @@ -88,7 +88,7 @@ class MonotonisedCentral : public Reconstructor { void reconstruct( gsl::not_null, Dim>*> vars_on_lower_face, gsl::not_null, Dim>*> vars_on_upper_face, - const Variables>& volume_vars, + const Variables>& volume_vars, const Element& element, const DirectionalIdMap& ghost_data, @@ -97,7 +97,7 @@ class MonotonisedCentral : public Reconstructor { template void reconstruct_fd_neighbor( gsl::not_null*> vars_on_face, - const Variables>& volume_vars, + const Variables>& volume_vars, const Element& element, const DirectionalIdMap& ghost_data, diff --git a/src/Evolution/Systems/ScalarAdvection/Subcell/CMakeLists.txt b/src/Evolution/Systems/ScalarAdvection/Subcell/CMakeLists.txt index 4ec97e1bedfa..938030321edb 100644 --- a/src/Evolution/Systems/ScalarAdvection/Subcell/CMakeLists.txt +++ b/src/Evolution/Systems/ScalarAdvection/Subcell/CMakeLists.txt @@ -5,6 +5,7 @@ spectre_target_sources( ${LIBRARY} PRIVATE GhostData.cpp + NeighborPackagedData.cpp SetInitialRdmpData.cpp TciOnDgGrid.cpp TciOnFdGrid.cpp diff --git a/src/Evolution/Systems/ScalarAdvection/Subcell/NeighborPackagedData.cpp b/src/Evolution/Systems/ScalarAdvection/Subcell/NeighborPackagedData.cpp new file mode 100644 index 000000000000..f3158e62a50b --- /dev/null +++ b/src/Evolution/Systems/ScalarAdvection/Subcell/NeighborPackagedData.cpp @@ -0,0 +1,240 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/Systems/ScalarAdvection/Subcell/NeighborPackagedData.hpp" + +#include +#include +#include +#include +#include + +#include "DataStructures/DataBox/Access.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Index.hpp" +#include "DataStructures/Tensor/Slice.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "Domain/Structure/Direction.hpp" +#include "Domain/Structure/DirectionalIdMap.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Domain/TagsTimeDependent.hpp" +#include "Evolution/BoundaryCorrectionTags.hpp" +#include "Evolution/DgSubcell/Mesh.hpp" +#include "Evolution/DgSubcell/NeighborReconstructedFaceSolution.tpp" +#include "Evolution/DgSubcell/Projection.hpp" +#include "Evolution/DgSubcell/Reconstruction.hpp" +#include "Evolution/DgSubcell/ReconstructionMethod.hpp" +#include "Evolution/DgSubcell/SubcellOptions.hpp" +#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" +#include "Evolution/DgSubcell/Tags/Mesh.hpp" +#include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp" +#include "Evolution/DgSubcell/Tags/SubcellOptions.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp" +#include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp" +#include "Evolution/Systems/ScalarAdvection/BoundaryCorrections/BoundaryCorrection.hpp" +#include "Evolution/Systems/ScalarAdvection/BoundaryCorrections/Factory.hpp" +#include "Evolution/Systems/ScalarAdvection/FiniteDifference/Factory.hpp" +#include "Evolution/Systems/ScalarAdvection/FiniteDifference/Reconstructor.hpp" +#include "Evolution/Systems/ScalarAdvection/FiniteDifference/Tags.hpp" +#include "Evolution/Systems/ScalarAdvection/Fluxes.hpp" +#include "Evolution/Systems/ScalarAdvection/Subcell/ComputeFluxes.hpp" +#include "Evolution/Systems/ScalarAdvection/Subcell/VelocityAtFace.hpp" +#include "Evolution/Systems/ScalarAdvection/System.hpp" +#include "Evolution/Systems/ScalarAdvection/Tags.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "Utilities/CallWithDynamicType.hpp" +#include "Utilities/GenerateInstantiations.hpp" +#include "Utilities/TMPL.hpp" + +namespace ScalarAdvection::subcell { +template +DirectionalIdMap NeighborPackagedData::apply( + const db::Access& box, + const std::vector>& mortars_to_reconstruct_to) { + // The object to return + DirectionalIdMap neighbor_package_data{}; + if (mortars_to_reconstruct_to.empty()) { + return neighbor_package_data; + } + + using evolved_vars_tags = typename System::variables_tag::tags_list; + using fluxes_tags = typename Fluxes::return_tags; + + // subcell currently does not support moving mesh + ASSERT(not db::get>(box).has_value(), + "Haven't yet added support for moving mesh to DG-subcell. This " + "should be easy to generalize, but we will want to consider " + "storing the mesh velocity on the faces instead of " + "re-slicing/projecting."); + + // Project volume variables from DG to subcell mesh + const Mesh& dg_mesh = db::get>(box); + const Mesh& subcell_mesh = + db::get>(box); + const auto& subcell_options = + db::get>(box); + const auto volume_vars_subcell = evolution::dg::subcell::fd::project( + db::get::variables_tag>(box), dg_mesh, + subcell_mesh.extents()); + + const auto& ghost_subcell_data = + db::get>( + box); + + const ScalarAdvection::fd::Reconstructor& recons = + db::get>(box); + + const auto& boundary_correction = + db::get>>(box); + using derived_boundary_corrections = + typename std::decay_t::creatable_classes; + + // perform reconstruction + tmpl::for_each([&](auto derived_correction_v) { + using derived_correction = tmpl::type_from; + if (typeid(boundary_correction) == typeid(derived_correction)) { + using dg_package_field_tags = + typename derived_correction::dg_package_field_tags; + using dg_package_data_temporary_tags = + typename derived_correction::dg_package_data_temporary_tags; + using dg_package_data_argument_tags = + tmpl::append; + const auto& element = db::get>(box); + + // Variables to store packaged data + Variables packaged_data{}; + // Variables to be reconstructed on the shared interfaces + Variables vars_on_face{}; + + for (const auto& mortar_id : mortars_to_reconstruct_to) { + const Direction& direction = mortar_id.direction; + Index extents = subcell_mesh.extents(); + + // Switch to face-centered instead of cell-centered points on the + // FD. There are num_cell_centered+1 face-centered points. + ++extents[direction.dimension()]; + const size_t num_face_pts = + subcell_mesh.extents().slice_away(direction.dimension()).product(); + vars_on_face.initialize(num_face_pts); + + using velocity_field = ScalarAdvection::Tags::VelocityField; + const auto& velocity_on_face = db::get< + evolution::dg::subcell::Tags::OnSubcellFaces>( + box); + data_on_slice(make_not_null(&get(vars_on_face)), + gsl::at(velocity_on_face, direction.dimension()), extents, + direction.dimension(), + direction.side() == Side::Lower + ? 0 + : extents[direction.dimension()] - 1); + + // Reconstruct field variables on faces + call_with_dynamic_type::creatable_classes>( + &recons, + [&element, &mortar_id, &ghost_subcell_data, &subcell_mesh, + &vars_on_face, &volume_vars_subcell](const auto& reconstructor) { + reconstructor->reconstruct_fd_neighbor( + make_not_null(&vars_on_face), volume_vars_subcell, element, + ghost_subcell_data, subcell_mesh, mortar_id.direction); + }); + + // Compute fluxes + ScalarAdvection::subcell::compute_fluxes( + make_not_null(&vars_on_face)); + + tnsr::i normal_covector = + get>( + *db::get>( + box) + .at(mortar_id.direction)); + for (auto& t : normal_covector) { + t *= -1.0; + } + + // Note: need to do the projection in 2D + if constexpr (Dim > 1) { + const auto dg_normal_covector = normal_covector; + for (size_t i = 0; i < Dim; ++i) { + normal_covector.get(i) = evolution::dg::subcell::fd::project( + dg_normal_covector.get(i), + dg_mesh.slice_away(mortar_id.direction.dimension()), + subcell_mesh.extents().slice_away( + mortar_id.direction.dimension())); + } + } + + // Compute the packaged data + packaged_data.initialize(num_face_pts); + evolution::dg::Actions::detail::dg_package_data>( + make_not_null(&packaged_data), + dynamic_cast(boundary_correction), + vars_on_face, normal_covector, {std::nullopt}, box, + typename derived_correction::dg_package_data_volume_tags{}, + dg_package_data_argument_tags{}); + + // Make a view so we can use iterators with std::copy + DataVector packaged_data_view{packaged_data.data(), + packaged_data.size()}; + neighbor_package_data[mortar_id] = DataVector{packaged_data.size()}; + std::copy(packaged_data_view.begin(), packaged_data_view.end(), + neighbor_package_data[mortar_id].begin()); + + // Note : need to reconstruct from FD to DG grid in 2D + if constexpr (Dim > 1) { + const auto dg_packaged_data = evolution::dg::subcell::fd::reconstruct( + packaged_data, + dg_mesh.slice_away(mortar_id.direction.dimension()), + subcell_mesh.extents().slice_away( + mortar_id.direction.dimension()), + subcell_options.reconstruction_method()); + // Make a view so we can use iterators with std::copy + DataVector dg_packaged_data_view{ + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast) + const_cast(dg_packaged_data.data()), + dg_packaged_data.size()}; + neighbor_package_data[mortar_id] = + DataVector{dg_packaged_data.size()}; + std::copy(dg_packaged_data_view.begin(), dg_packaged_data_view.end(), + neighbor_package_data[mortar_id].begin()); + } + } + } + }); + + return neighbor_package_data; +} + +#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) + +#define INSTANTIATION(r, data) \ + template DirectionalIdMap \ + NeighborPackagedData::apply( \ + const db::Access& box, \ + const std::vector>& mortars_to_reconstruct_to); + +GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2)) + +#undef INSTANTIATION + +} // namespace ScalarAdvection::subcell + +#define INSTANTIATION(r, data) \ + template void evolution::dg::subcell::neighbor_reconstructed_face_solution< \ + DIM(data), ScalarAdvection::subcell::NeighborPackagedData>( \ + gsl::not_null box, \ + gsl::not_null, Mesh, \ + std::optional, std::optional, \ + ::TimeStepId, int>>>*> \ + received_temporal_id_and_data); + +GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2)) + +#undef INSTANTIATION +#undef DIM diff --git a/src/Evolution/Systems/ScalarAdvection/Subcell/NeighborPackagedData.hpp b/src/Evolution/Systems/ScalarAdvection/Subcell/NeighborPackagedData.hpp index 107d65e9e41c..bcdb1c278cdf 100644 --- a/src/Evolution/Systems/ScalarAdvection/Subcell/NeighborPackagedData.hpp +++ b/src/Evolution/Systems/ScalarAdvection/Subcell/NeighborPackagedData.hpp @@ -4,45 +4,18 @@ #pragma once #include -#include -#include -#include #include -#include "DataStructures/DataBox/DataBox.hpp" -#include "DataStructures/DataVector.hpp" -#include "DataStructures/Index.hpp" -#include "DataStructures/Tensor/Slice.hpp" -#include "DataStructures/Tensor/Tensor.hpp" -#include "DataStructures/Variables.hpp" -#include "Domain/Structure/Direction.hpp" -#include "Domain/Structure/DirectionalIdMap.hpp" -#include "Domain/Structure/ElementId.hpp" -#include "Domain/TagsTimeDependent.hpp" -#include "Evolution/BoundaryCorrectionTags.hpp" -#include "Evolution/DgSubcell/Mesh.hpp" -#include "Evolution/DgSubcell/Projection.hpp" -#include "Evolution/DgSubcell/Reconstruction.hpp" -#include "Evolution/DgSubcell/ReconstructionMethod.hpp" -#include "Evolution/DgSubcell/SubcellOptions.hpp" -#include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" -#include "Evolution/DgSubcell/Tags/Mesh.hpp" -#include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp" -#include "Evolution/DgSubcell/Tags/SubcellOptions.hpp" -#include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp" -#include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp" -#include "Evolution/Systems/ScalarAdvection/BoundaryCorrections/BoundaryCorrection.hpp" -#include "Evolution/Systems/ScalarAdvection/BoundaryCorrections/Factory.hpp" -#include "Evolution/Systems/ScalarAdvection/FiniteDifference/Reconstructor.hpp" -#include "Evolution/Systems/ScalarAdvection/FiniteDifference/Tags.hpp" -#include "Evolution/Systems/ScalarAdvection/Fluxes.hpp" -#include "Evolution/Systems/ScalarAdvection/Subcell/ComputeFluxes.hpp" -#include "Evolution/Systems/ScalarAdvection/Subcell/VelocityAtFace.hpp" -#include "Evolution/Systems/ScalarAdvection/System.hpp" -#include "Evolution/Systems/ScalarAdvection/Tags.hpp" -#include "NumericalAlgorithms/Spectral/Mesh.hpp" -#include "Utilities/CallWithDynamicType.hpp" -#include "Utilities/TMPL.hpp" +/// \cond +namespace db { +class Access; +} // namespace db +class DataVector; +template +struct DirectionalId; +template +class DirectionalIdMap; +/// \endcond namespace ScalarAdvection::subcell { /*! @@ -61,169 +34,9 @@ namespace ScalarAdvection::subcell { */ struct NeighborPackagedData { - template + template static DirectionalIdMap apply( - const db::DataBox& box, - const std::vector>& mortars_to_reconstruct_to) { - // The object to return - DirectionalIdMap neighbor_package_data{}; - if (mortars_to_reconstruct_to.empty()) { - return neighbor_package_data; - } - - using evolved_vars_tags = typename System::variables_tag::tags_list; - using fluxes_tags = typename Fluxes::return_tags; - - // subcell currently does not support moving mesh - ASSERT(not db::get>(box).has_value(), - "Haven't yet added support for moving mesh to DG-subcell. This " - "should be easy to generalize, but we will want to consider " - "storing the mesh velocity on the faces instead of " - "re-slicing/projecting."); - - // Project volume variables from DG to subcell mesh - const Mesh& dg_mesh = db::get>(box); - const Mesh& subcell_mesh = - db::get>(box); - const auto& subcell_options = - db::get>(box); - const auto volume_vars_subcell = evolution::dg::subcell::fd::project( - db::get::variables_tag>(box), dg_mesh, - subcell_mesh.extents()); - - const auto& ghost_subcell_data = - db::get>( - box); - - const ScalarAdvection::fd::Reconstructor& recons = - db::get>(box); - - const auto& boundary_correction = - db::get>>(box); - using derived_boundary_corrections = - typename std::decay_t::creatable_classes; - - // perform reconstruction - tmpl::for_each([&](auto - derived_correction_v) { - using derived_correction = - tmpl::type_from; - if (typeid(boundary_correction) == typeid(derived_correction)) { - using dg_package_field_tags = - typename derived_correction::dg_package_field_tags; - using dg_package_data_temporary_tags = - typename derived_correction::dg_package_data_temporary_tags; - using dg_package_data_argument_tags = - tmpl::append; - - const auto& element = db::get>(box); - - // Variables to store packaged data - Variables packaged_data{}; - // Variables to be reconstructed on the shared interfaces - Variables vars_on_face{}; - - for (const auto& mortar_id : mortars_to_reconstruct_to) { - const Direction& direction = mortar_id.direction; - Index extents = subcell_mesh.extents(); - - // Switch to face-centered instead of cell-centered points on the - // FD. There are num_cell_centered+1 face-centered points. - ++extents[direction.dimension()]; - const size_t num_face_pts = subcell_mesh.extents() - .slice_away(direction.dimension()) - .product(); - vars_on_face.initialize(num_face_pts); - - using velocity_field = ScalarAdvection::Tags::VelocityField; - const auto& velocity_on_face = - db::get>(box); - data_on_slice(make_not_null(&get(vars_on_face)), - gsl::at(velocity_on_face, direction.dimension()), - extents, direction.dimension(), - direction.side() == Side::Lower - ? 0 - : extents[direction.dimension()] - 1); - - // Reconstruct field variables on faces - call_with_dynamic_type::creatable_classes>( - &recons, - [&element, &mortar_id, &ghost_subcell_data, &subcell_mesh, - &vars_on_face, &volume_vars_subcell](const auto& reconstructor) { - reconstructor->reconstruct_fd_neighbor( - make_not_null(&vars_on_face), volume_vars_subcell, element, - ghost_subcell_data, subcell_mesh, mortar_id.direction); - }); - - // Compute fluxes - ScalarAdvection::subcell::compute_fluxes( - make_not_null(&vars_on_face)); - - tnsr::i normal_covector = get< - evolution::dg::Tags::NormalCovector>( - *db::get>( - box) - .at(mortar_id.direction)); - for (auto& t : normal_covector) { - t *= -1.0; - } - - // Note: need to do the projection in 2D - if constexpr (Dim > 1) { - const auto dg_normal_covector = normal_covector; - for (size_t i = 0; i < Dim; ++i) { - normal_covector.get(i) = evolution::dg::subcell::fd::project( - dg_normal_covector.get(i), - dg_mesh.slice_away(mortar_id.direction.dimension()), - subcell_mesh.extents().slice_away( - mortar_id.direction.dimension())); - } - } - - // Compute the packaged data - packaged_data.initialize(num_face_pts); - evolution::dg::Actions::detail::dg_package_data>( - make_not_null(&packaged_data), - dynamic_cast(boundary_correction), - vars_on_face, normal_covector, {std::nullopt}, box, - typename derived_correction::dg_package_data_volume_tags{}, - dg_package_data_argument_tags{}); - - // Make a view so we can use iterators with std::copy - DataVector packaged_data_view{packaged_data.data(), - packaged_data.size()}; - neighbor_package_data[mortar_id] = DataVector{packaged_data.size()}; - std::copy(packaged_data_view.begin(), packaged_data_view.end(), - neighbor_package_data[mortar_id].begin()); - - // Note : need to reconstruct from FD to DG grid in 2D - if constexpr (Dim > 1) { - const auto dg_packaged_data = - evolution::dg::subcell::fd::reconstruct( - packaged_data, - dg_mesh.slice_away(mortar_id.direction.dimension()), - subcell_mesh.extents().slice_away( - mortar_id.direction.dimension()), - subcell_options.reconstruction_method()); - // Make a view so we can use iterators with std::copy - DataVector dg_packaged_data_view{ - const_cast(dg_packaged_data.data()), - dg_packaged_data.size()}; - neighbor_package_data[mortar_id] = - DataVector{dg_packaged_data.size()}; - std::copy(dg_packaged_data_view.begin(), - dg_packaged_data_view.end(), - neighbor_package_data[mortar_id].begin()); - } - } - } - }); - - return neighbor_package_data; - } + const db::Access& box, + const std::vector>& mortars_to_reconstruct_to); }; } // namespace ScalarAdvection::subcell diff --git a/src/Evolution/Systems/ScalarTensor/BoundaryCorrections/ProductOfCorrections.hpp b/src/Evolution/Systems/ScalarTensor/BoundaryCorrections/ProductOfCorrections.hpp index 14afdeeae5b1..c05ed75279d1 100644 --- a/src/Evolution/Systems/ScalarTensor/BoundaryCorrections/ProductOfCorrections.hpp +++ b/src/Evolution/Systems/ScalarTensor/BoundaryCorrections/ProductOfCorrections.hpp @@ -52,6 +52,10 @@ class ProductOfCorrections final : public BoundaryCorrection { typename DerivedGhCorrection::dg_package_data_volume_tags, typename DerivedScalarCorrection::dg_package_data_volume_tags>; + using dg_boundary_terms_volume_tags = tmpl::append< + typename DerivedGhCorrection::dg_boundary_terms_volume_tags, + typename DerivedScalarCorrection::dg_boundary_terms_volume_tags>; + static std::string name() { return "Product" + pretty_type::name() + "GH" + "And" + pretty_type::name() + "Scalar"; diff --git a/src/Evolution/Systems/ScalarWave/BoundaryCorrections/UpwindPenalty.hpp b/src/Evolution/Systems/ScalarWave/BoundaryCorrections/UpwindPenalty.hpp index 0917b9a2d6fb..12b01f29f01e 100644 --- a/src/Evolution/Systems/ScalarWave/BoundaryCorrections/UpwindPenalty.hpp +++ b/src/Evolution/Systems/ScalarWave/BoundaryCorrections/UpwindPenalty.hpp @@ -214,6 +214,7 @@ class UpwindPenalty final : public BoundaryCorrection { CharSpeedsTensor>; using dg_package_data_temporary_tags = tmpl::list; using dg_package_data_volume_tags = tmpl::list<>; + using dg_boundary_terms_volume_tags = tmpl::list<>; double dg_package_data( gsl::not_null*> packaged_char_speed_v_psi, diff --git a/tests/Unit/Evolution/DgSubcell/Test_NeighborReconstructedFaceSolution.cpp b/tests/Unit/Evolution/DgSubcell/Test_NeighborReconstructedFaceSolution.cpp index b996bc63d3d3..378e5e5fcf1b 100644 --- a/tests/Unit/Evolution/DgSubcell/Test_NeighborReconstructedFaceSolution.cpp +++ b/tests/Unit/Evolution/DgSubcell/Test_NeighborReconstructedFaceSolution.cpp @@ -9,6 +9,7 @@ #include #include +#include "DataStructures/DataBox/Access.hpp" #include "DataStructures/DataBox/DataBox.hpp" #include "DataStructures/DataBox/Tag.hpp" #include "DataStructures/DataVector.hpp" @@ -17,6 +18,7 @@ #include "Domain/Structure/ElementId.hpp" #include "Evolution/DgSubcell/GhostData.hpp" #include "Evolution/DgSubcell/NeighborReconstructedFaceSolution.hpp" +#include "Evolution/DgSubcell/NeighborReconstructedFaceSolution.tpp" #include "Evolution/DgSubcell/RdmpTciData.hpp" #include "Evolution/DgSubcell/Tags/DataForRdmpTci.hpp" #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" @@ -49,11 +51,9 @@ struct Metavariables { static constexpr size_t volume_dim = Dim; struct SubcellOptions { struct DgComputeSubcellNeighborPackagedData { - template static NeighborReconstructionMap apply( - const db::DataBox& box, - const std::vector>& - mortars_to_reconstruct_to) { + const db::Access& box, const std::vector>& + mortars_to_reconstruct_to) { const GhostDataMap& ghost_data = db::get< evolution::dg::subcell::Tags::GhostDataForReconstruction>(box); diff --git a/tests/Unit/Evolution/DiscontinuousGalerkin/Actions/Test_ApplyBoundaryCorrections.cpp b/tests/Unit/Evolution/DiscontinuousGalerkin/Actions/Test_ApplyBoundaryCorrections.cpp index 3c97221494ba..c5c3b6ed84d8 100644 --- a/tests/Unit/Evolution/DiscontinuousGalerkin/Actions/Test_ApplyBoundaryCorrections.cpp +++ b/tests/Unit/Evolution/DiscontinuousGalerkin/Actions/Test_ApplyBoundaryCorrections.cpp @@ -62,6 +62,10 @@ struct Var2 : db::SimpleTag { using type = tnsr::I; }; +struct VolumeTag : db::SimpleTag { + using type = int; +}; + template struct BoundaryTerms; @@ -110,6 +114,7 @@ struct BoundaryTerms final : public BoundaryCorrection { tmpl::append, variables_tags>, MaxAbsCharSpeed>; + using dg_boundary_terms_volume_tags = tmpl::list; void dg_boundary_terms( const gsl::not_null*> boundary_correction_var1, @@ -127,7 +132,7 @@ struct BoundaryTerms final : public BoundaryCorrection { const Scalar& exterior_var1, const tnsr::I& exterior_var2, const Scalar& exterior_max_abs_char_speed, - const dg::Formulation dg_formulation) const { + const dg::Formulation dg_formulation, const int& volume_tag) const { // extra minus sign on exterior normal dot flux because normal faces // opposite direction get(*boundary_correction_var1) = @@ -150,6 +155,7 @@ struct BoundaryTerms final : public BoundaryCorrection { get(exterior_max_abs_char_speed)) * (exterior_var2.get(i) - interior_var2.get(i)); } + CHECK(volume_tag == 10); } }; @@ -397,8 +403,8 @@ struct component { domain::Tags::BoundaryDirectionsInterior; using simple_tags = tmpl::list< - ::Tags::TimeStepId, ::Tags::Next<::Tags::TimeStepId>, ::Tags::TimeStep, - Tags::ConcreteTimeStepper, + VolumeTag, ::Tags::TimeStepId, ::Tags::Next<::Tags::TimeStepId>, + ::Tags::TimeStep, Tags::ConcreteTimeStepper, db::add_tag_prefix<::Tags::dt, typename Metavariables::system::variables_tag>, typename Metavariables::system::variables_tag, @@ -583,7 +589,7 @@ void test_impl(const Spectral::Quadrature quadrature, ActionTesting::emplace_component_and_initialize( &runner, self_id, - {time_step_id, local_next_time_step_id, time_step, + {10, time_step_id, local_next_time_step_id, time_step, std::make_unique(time_stepper), dt_evolved_vars, evolved_vars, mesh, element, inertial_coords, inv_jac, quadrature, typename domain::Tags::NeighborMesh::type{}}); @@ -824,7 +830,7 @@ void test_impl(const Spectral::Quadrature quadrature, get>(neighbor_data_on_mortar), get::MaxAbsCharSpeed>( neighbor_data_on_mortar), - dg_formulation); + dg_formulation, 10); // Project the boundary terms from the mortar to the face const std::array& mortar_size = diff --git a/tests/Unit/Evolution/Systems/Burgers/Subcell/Test_NeighborPackagedData.cpp b/tests/Unit/Evolution/Systems/Burgers/Subcell/Test_NeighborPackagedData.cpp index 3c04ca667dbf..8e5f093162e5 100644 --- a/tests/Unit/Evolution/Systems/Burgers/Subcell/Test_NeighborPackagedData.cpp +++ b/tests/Unit/Evolution/Systems/Burgers/Subcell/Test_NeighborPackagedData.cpp @@ -35,6 +35,7 @@ #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp" #include "Evolution/DgSubcell/Tags/Mesh.hpp" #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp" #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp" #include "Evolution/Systems/Burgers/BoundaryCorrections/BoundaryCorrection.hpp" #include "Evolution/Systems/Burgers/BoundaryCorrections/Factory.hpp" @@ -42,6 +43,7 @@ #include "Evolution/Systems/Burgers/FiniteDifference/Reconstructor.hpp" #include "Evolution/Systems/Burgers/FiniteDifference/Tags.hpp" #include "Evolution/Systems/Burgers/Fluxes.hpp" +#include "Evolution/Systems/Burgers/Subcell/ComputeFluxes.hpp" #include "Evolution/Systems/Burgers/Subcell/NeighborPackagedData.hpp" #include "Evolution/Systems/Burgers/System.hpp" #include "Framework/TestHelpers.hpp" diff --git a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp index 72798759226c..38aef189c285 100644 --- a/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp +++ b/tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_NeighborPackagedData.cpp @@ -41,6 +41,7 @@ #include "Evolution/DgSubcell/Tags/Mesh.hpp" #include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp" #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp" #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp" #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp" #include "Evolution/Systems/GeneralizedHarmonic/BoundaryCorrections/UpwindPenalty.hpp" diff --git a/tests/Unit/Evolution/Systems/NewtonianEuler/Subcell/Test_NeighborPackagedData.cpp b/tests/Unit/Evolution/Systems/NewtonianEuler/Subcell/Test_NeighborPackagedData.cpp index 1b98ca1d319c..8e4f0d724823 100644 --- a/tests/Unit/Evolution/Systems/NewtonianEuler/Subcell/Test_NeighborPackagedData.cpp +++ b/tests/Unit/Evolution/Systems/NewtonianEuler/Subcell/Test_NeighborPackagedData.cpp @@ -52,9 +52,11 @@ #include "Evolution/Systems/NewtonianEuler/FiniteDifference/MonotonisedCentral.hpp" #include "Evolution/Systems/NewtonianEuler/FiniteDifference/Tag.hpp" #include "Evolution/Systems/NewtonianEuler/Subcell/NeighborPackagedData.hpp" +#include "Evolution/Systems/NewtonianEuler/System.hpp" #include "Evolution/Systems/NewtonianEuler/Tags.hpp" #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "Parallel/Tags/Metavariables.hpp" #include "PointwiseFunctions/AnalyticSolutions/NewtonianEuler/SmoothFlow.hpp" #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp" #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp" diff --git a/tests/Unit/Evolution/Systems/ScalarAdvection/Subcell/Test_NeighborPackagedData.cpp b/tests/Unit/Evolution/Systems/ScalarAdvection/Subcell/Test_NeighborPackagedData.cpp index 95d3fbbf0581..496d993d5867 100644 --- a/tests/Unit/Evolution/Systems/ScalarAdvection/Subcell/Test_NeighborPackagedData.cpp +++ b/tests/Unit/Evolution/Systems/ScalarAdvection/Subcell/Test_NeighborPackagedData.cpp @@ -41,6 +41,7 @@ #include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp" #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp" #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp" +#include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp" #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp" #include "Evolution/Initialization/Tags.hpp" #include "Evolution/Systems/ScalarAdvection/BoundaryCorrections/BoundaryCorrection.hpp" @@ -49,6 +50,7 @@ #include "Evolution/Systems/ScalarAdvection/FiniteDifference/Reconstructor.hpp" #include "Evolution/Systems/ScalarAdvection/FiniteDifference/Tags.hpp" #include "Evolution/Systems/ScalarAdvection/Fluxes.hpp" +#include "Evolution/Systems/ScalarAdvection/Subcell/ComputeFluxes.hpp" #include "Evolution/Systems/ScalarAdvection/Subcell/NeighborPackagedData.hpp" #include "Evolution/Systems/ScalarAdvection/Subcell/VelocityAtFace.hpp" #include "Evolution/Systems/ScalarAdvection/System.hpp"