From 2f4125641c8cb5433c4141e07a75a4b80b2ff45f Mon Sep 17 00:00:00 2001 From: pajkosmi Date: Mon, 30 Jan 2023 17:08:51 -0800 Subject: [PATCH] Add refinement option to sphere domain wedges Co-authored-by: Nils Vu --- src/Domain/Creators/Python/Sphere.cpp | 2 +- src/Domain/Creators/Sphere.cpp | 158 +++- src/Domain/Creators/Sphere.hpp | 105 ++- src/Domain/DomainHelpers.cpp | 16 +- src/Domain/DomainHelpers.hpp | 41 +- .../Punctures/MultiplePunctures.yaml | 4 +- tests/InputFiles/Xcts/TovStar.yaml | 4 +- .../Domain/Creators/Python/Test_Sphere.py | 2 +- tests/Unit/Domain/Creators/Test_Sphere.cpp | 730 ++++++++---------- tests/Unit/Domain/Test_DomainHelpers.cpp | 23 - .../Helpers/Domain/Creators/TestHelpers.hpp | 14 + 11 files changed, 577 insertions(+), 522 deletions(-) diff --git a/src/Domain/Creators/Python/Sphere.cpp b/src/Domain/Creators/Python/Sphere.cpp index 1917f0974cb8..0dd356130baa 100644 --- a/src/Domain/Creators/Python/Sphere.cpp +++ b/src/Domain/Creators/Python/Sphere.cpp @@ -17,7 +17,7 @@ namespace py = pybind11; namespace domain::creators::py_bindings { void bind_sphere(py::module& m) { py::class_>(m, "Sphere") - .def(py::init, + .def(py::init, bool>(), py::arg("inner_radius"), py::arg("outer_radius"), py::arg("inner_cube_sphericity"), py::arg("initial_refinement"), diff --git a/src/Domain/Creators/Sphere.cpp b/src/Domain/Creators/Sphere.cpp index 0054944afe45..567ba8963efd 100644 --- a/src/Domain/Creators/Sphere.cpp +++ b/src/Domain/Creators/Sphere.cpp @@ -7,7 +7,9 @@ #include #include #include +#include #include +#include #include #include "Domain/Block.hpp" @@ -22,6 +24,7 @@ #include "Domain/CoordinateMaps/ProductMaps.tpp" #include "Domain/CoordinateMaps/Wedge.hpp" #include "Domain/Creators/DomainCreator.hpp" +#include "Domain/Creators/ExpandOverBlocks.hpp" #include "Domain/Creators/TimeDependence/None.hpp" #include "Domain/Domain.hpp" #include "Domain/DomainHelpers.hpp" @@ -39,26 +42,24 @@ struct BlockLogical; namespace domain::creators { Sphere::Sphere( - typename InnerRadius::type inner_radius, - typename OuterRadius::type outer_radius, const double inner_cube_sphericity, - typename InitialRefinement::type initial_refinement, - typename InitialGridPoints::type initial_number_of_grid_points, - typename UseEquiangularMap::type use_equiangular_map, + double inner_radius, double outer_radius, double inner_cube_sphericity, + const typename InitialRefinement::type& initial_refinement, + const typename InitialGridPoints::type& initial_number_of_grid_points, + bool use_equiangular_map, std::vector radial_partitioning, + std::vector radial_distribution, std::unique_ptr> time_dependence, std::unique_ptr boundary_condition, const Options::Context& context) - // clang-tidy: trivially copyable - : inner_radius_(std::move(inner_radius)), // NOLINT - outer_radius_(std::move(outer_radius)), // NOLINT + : inner_radius_(inner_radius), + outer_radius_(outer_radius), inner_cube_sphericity_(inner_cube_sphericity), - initial_refinement_( // NOLINT - std::move(initial_refinement)), // NOLINT - initial_number_of_grid_points_( // NOLINT - std::move(initial_number_of_grid_points)), // NOLINT - use_equiangular_map_(std::move(use_equiangular_map)), // NOLINT - time_dependence_(std::move(time_dependence)), // NOLINT + use_equiangular_map_(use_equiangular_map), + radial_partitioning_(std::move(radial_partitioning)), + radial_distribution_(std::move(radial_distribution)), + time_dependence_(std::move(time_dependence)), + boundary_condition_(std::move(boundary_condition)) { if (inner_cube_sphericity_ < 0.0 or inner_cube_sphericity_ >= 1.0) { PARSE_ERROR( @@ -66,6 +67,96 @@ Sphere::Sphere( "Inner cube sphericity must be >= 0.0 and strictly < 1.0, not " + get_output(inner_cube_sphericity_)); } + + if (inner_radius_ > outer_radius_) { + PARSE_ERROR(context, + "Inner radius must be smaller than outer radius, but inner " + "radius is " + + std::to_string(inner_radius_) + " and outer radius is " + + std::to_string(outer_radius_) + "."); + } + if (not std::is_sorted(radial_partitioning_.begin(), + radial_partitioning_.end())) { + PARSE_ERROR(context, + "Specify radial partitioning in ascending order. Specified " + "radial partitioning is: " + + get_output(radial_partitioning_)); + } + if (not radial_partitioning_.empty()) { + if (radial_partitioning_.front() <= inner_radius_) { + PARSE_ERROR(context, + "First radial partition must be larger than inner " + "radius, but is: " + + std::to_string(inner_radius_)); + } + if (radial_partitioning_.back() >= outer_radius_) { + PARSE_ERROR(context, + "Last radial partition must be smaller than outer " + "radius, but is: " + + std::to_string(outer_radius_)); + } + } + + const size_t num_shells = 1 + radial_partitioning_.size(); + if (radial_distribution_.size() != num_shells) { + PARSE_ERROR(context, + "Specify a 'RadialDistribution' for every spherical shell. You " + "specified " + << radial_distribution_.size() + << " items, but the domain has " << num_shells + << " shells."); + } + if (radial_distribution_.front() != + domain::CoordinateMaps::Distribution::Linear) { + PARSE_ERROR(context, + "The 'RadialDistribution' must be 'Linear' for the innermost " + "shell because it changes in sphericity. Add entries to " + "'RadialPartitioning' to add outer shells for which you can " + "select different radial distributions."); + } + + // Create block names and groups + static std::array wedge_directions{ + "UpperZ", "LowerZ", "UpperY", "LowerY", "UpperX", "LowerX"}; + for (size_t shell = 0; shell < num_shells; ++shell) { + std::string shell_prefix = "Shell" + std::to_string(shell); + for (size_t direction = 0; direction < 6; ++direction) { + const std::string wedge_name = + shell_prefix + gsl::at(wedge_directions, direction); + block_names_.emplace_back(wedge_name); + if (num_shells > 1) { + block_groups_[shell_prefix].insert(wedge_name); + } + block_groups_["Wedges"].insert(wedge_name); + } + } + block_names_.emplace_back("InnerCube"); + + // Expand initial refinement and number of grid points over all blocks + const ExpandOverBlocks expand_over_blocks{block_names_, + block_groups_}; + try { + initial_refinement_ = std::visit(expand_over_blocks, initial_refinement); + } catch (const std::exception& error) { + PARSE_ERROR(context, "Invalid 'InitialRefinement': " << error.what()); + } + try { + initial_number_of_grid_points_ = + std::visit(expand_over_blocks, initial_number_of_grid_points); + } catch (const std::exception& error) { + PARSE_ERROR(context, "Invalid 'InitialGridPoints': " << error.what()); + } + + // The central cube has no notion of a "radial" direction, so we set + // refinement and number of grid points of the central cube z direction to + // its y value, which corresponds to the azimuthal direction of the + // wedges. This keeps the boundaries conforming when the radial direction + // is chosen differently to the angular directions. + auto& central_cube_refinement = initial_refinement_.back(); + auto& central_cube_grid_points = initial_number_of_grid_points_.back(); + central_cube_refinement[2] = central_cube_refinement[1]; + central_cube_grid_points[2] = central_cube_grid_points[1]; + using domain::BoundaryConditions::is_none; if (is_none(boundary_condition_)) { PARSE_ERROR( @@ -85,14 +176,14 @@ Sphere::Sphere( } Domain<3> Sphere::create_domain() const { + const size_t num_shells = 1 + radial_partitioning_.size(); std::vector> corners = - corners_for_radially_layered_domains(1, true); + corners_for_radially_layered_domains(num_shells, true); - auto coord_maps = domain::make_vector_coordinate_map_base( - sph_wedge_coordinate_maps(inner_radius_, outer_radius_, - inner_cube_sphericity_, 1.0, - use_equiangular_map_)); + auto coord_maps = domain::make_vector_coordinate_map_base< + Frame::BlockLogical, Frame::Inertial, 3>(sph_wedge_coordinate_maps( + inner_radius_, outer_radius_, inner_cube_sphericity_, 1.0, + use_equiangular_map_, false, radial_partitioning_, radial_distribution_)); if (inner_cube_sphericity_ == 0.0) { if (use_equiangular_map_) { coord_maps.emplace_back( @@ -121,7 +212,8 @@ Domain<3> Sphere::create_domain() const { use_equiangular_map_})); } - Domain<3> domain{std::move(coord_maps), corners}; + Domain<3> domain{std::move(coord_maps), corners, {}, {}, + block_names_, block_groups_}; if (not time_dependence_->is_none()) { const size_t number_of_blocks = domain.blocks().size(); @@ -148,31 +240,19 @@ Sphere::external_boundary_conditions() const { if (boundary_condition_ == nullptr) { return {}; } + + // number of blocks = 1 inner_block + 6 * (number of shells) + size_t number_of_blocks = 1 + 6 * (radial_partitioning_.size() + 1); + std::vector>> - boundary_conditions{7}; - for (size_t i = 0; i < 6; ++i) { + boundary_conditions{number_of_blocks}; + for (size_t i = number_of_blocks - 7; i < number_of_blocks - 1; ++i) { boundary_conditions[i][Direction<3>::upper_zeta()] = boundary_condition_->get_clone(); } return boundary_conditions; } - -std::vector> Sphere::initial_extents() const { - std::vector> extents{ - 6, - {{initial_number_of_grid_points_[1], initial_number_of_grid_points_[1], - initial_number_of_grid_points_[0]}}}; - extents.push_back( - {{initial_number_of_grid_points_[1], initial_number_of_grid_points_[1], - initial_number_of_grid_points_[1]}}); - return extents; -} - -std::vector> Sphere::initial_refinement_levels() const { - return {7, make_array<3>(initial_refinement_)}; -} - std::unordered_map> Sphere::functions_of_time(const std::unordered_map& diff --git a/src/Domain/Creators/Sphere.hpp b/src/Domain/Creators/Sphere.hpp index 8ab77f0399d7..ab8b0be0b4f6 100644 --- a/src/Domain/Creators/Sphere.hpp +++ b/src/Domain/Creators/Sphere.hpp @@ -6,10 +6,14 @@ #include #include #include +#include +#include +#include #include #include "Domain/BoundaryConditions/BoundaryCondition.hpp" #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp" +#include "Domain/CoordinateMaps/Distribution.hpp" #include "Domain/Creators/DomainCreator.hpp" #include "Domain/Creators/TimeDependence/TimeDependence.hpp" #include "Domain/Domain.hpp" @@ -83,21 +87,53 @@ class Sphere : public DomainCreator<3> { }; struct InitialRefinement { - using type = size_t; + using type = + std::variant, + std::vector>, + std::unordered_map>>; static constexpr Options::String help = { - "Initial refinement level in each dimension."}; + "Initial refinement level. Specify one of: a single number, a " + "list representing [phi, theta, r], or such a list for every block " + "in the domain. The central cube always uses the value for 'theta' " + "in both y- and z-direction."}; }; struct InitialGridPoints { - using type = std::array; + using type = + std::variant, + std::vector>, + std::unordered_map>>; static constexpr Options::String help = { - "Initial number of grid points in [r,angular]."}; + "Initial number of grid points. Specify one of: a single number, a " + "list representing [phi, theta, r], or such a list for every block " + "in the domain. The central cube always uses the value for 'theta' " + "in both y- and z-direction."}; }; struct UseEquiangularMap { using type = bool; static constexpr Options::String help = { - "Use equiangular instead of equidistant coordinates."}; + "Use equiangular instead of equidistant coordinates. Equiangular " + "coordinates give better gridpoint spacings in the angular " + "directions, while equidistant coordinates give better gridpoint " + "spacings in the center block."}; + }; + + struct RadialPartitioning { + using type = std::vector; + static constexpr Options::String help = { + "Radial coordinates of the boundaries splitting the spherical shell " + "between InnerRadius and OuterRadius."}; + }; + + struct RadialDistribution { + using type = std::vector; + static constexpr Options::String help = { + "Select the radial distribution of grid points in each spherical " + "shell. The innermost shell must have a 'Linear' distribution because " + "it changes in sphericity. The 'RadialPartitioning' determines the " + "number of shells."}; + static size_t lower_bound_on_size() { return 1; } }; struct TimeDependence { @@ -118,7 +154,7 @@ class Sphere : public DomainCreator<3> { using basic_options = tmpl::list; + RadialPartitioning, RadialDistribution, TimeDependence>; template using options = tmpl::conditional_t< @@ -132,22 +168,16 @@ class Sphere : public DomainCreator<3> { basic_options>; static constexpr Options::String help{ - "Creates a 3D Sphere with seven Blocks.\n" - "Only one refinement level for all dimensions is currently supported.\n" - "The number of gridpoints in the radial direction can be set\n" - "independently of the number of gridpoints in the angular directions.\n" - "The number of gridpoints along the dimensions of the cube is equal\n" - "to the number of gridpoints along the angular dimensions of the " - "wedges.\n" - "Equiangular coordinates give better gridpoint spacings in the angular\n" - "directions, while equidistant coordinates give better gridpoint\n" - "spacings in the center block."}; - - Sphere(typename InnerRadius::type inner_radius, - typename OuterRadius::type outer_radius, double inner_cube_sphericity, - typename InitialRefinement::type initial_refinement, - typename InitialGridPoints::type initial_number_of_grid_points, - typename UseEquiangularMap::type use_equiangular_map, + "A 3D cubed sphere. Six wedges surround a central cube. Additional " + "spherical shells, each composed of six wedges, can be added with the " + "'RadialPartitioning' option."}; + + Sphere(double inner_radius, double outer_radius, double inner_cube_sphericity, + const typename InitialRefinement::type& initial_refinement, + const typename InitialGridPoints::type& initial_number_of_grid_points, + bool use_equiangular_map, std::vector radial_partitioning = {}, + std::vector radial_distribution = + {domain::CoordinateMaps::Distribution::Linear}, std::unique_ptr> time_dependence = nullptr, std::unique_ptr @@ -167,9 +197,21 @@ class Sphere : public DomainCreator<3> { 3, std::unique_ptr>> external_boundary_conditions() const override; - std::vector> initial_extents() const override; + std::vector> initial_extents() const override { + return initial_number_of_grid_points_; + } + + std::vector> initial_refinement_levels() + const override { + return initial_refinement_; + } + + std::vector block_names() const override { return block_names_; } - std::vector> initial_refinement_levels() const override; + std::unordered_map> + block_groups() const override { + return block_groups_; + } auto functions_of_time(const std::unordered_map& initial_expiration_times = {}) const @@ -178,16 +220,21 @@ class Sphere : public DomainCreator<3> { std::unique_ptr> override; private: - typename InnerRadius::type inner_radius_{}; - typename OuterRadius::type outer_radius_{}; + double inner_radius_{}; + double outer_radius_{}; double inner_cube_sphericity_{}; - typename InitialRefinement::type initial_refinement_{}; - typename InitialGridPoints::type initial_number_of_grid_points_{}; - typename UseEquiangularMap::type use_equiangular_map_ = false; + std::vector> initial_refinement_{}; + std::vector> initial_number_of_grid_points_{}; + bool use_equiangular_map_ = false; + std::vector radial_partitioning_{}; + std::vector radial_distribution_{}; std::unique_ptr> time_dependence_; std::unique_ptr boundary_condition_; + std::vector block_names_{}; + std::unordered_map> + block_groups_{}; }; } // namespace creators } // namespace domain diff --git a/src/Domain/DomainHelpers.cpp b/src/Domain/DomainHelpers.cpp index fea0f6c75be1..34312bccabf3 100644 --- a/src/Domain/DomainHelpers.cpp +++ b/src/Domain/DomainHelpers.cpp @@ -608,9 +608,6 @@ std::vector> sph_wedge_coordinate_maps( const ShellWedges which_wedges) { ASSERT(not use_half_wedges or which_wedges == ShellWedges::All, "If we are using half wedges we must also be using ShellWedges::All."); - ASSERT(radial_partitioning.empty() or inner_sphericity == outer_sphericity, - "If we are using more than one layer the inner and outer sphericities " - "must match."); ASSERT(radial_distribution.size() == 1 + radial_partitioning.size(), "Specify a radial distribution for every spherical shell. You " "specified " @@ -626,6 +623,7 @@ std::vector> sph_wedge_coordinate_maps( const size_t number_of_layers = 1 + radial_partitioning.size(); double temp_inner_radius = inner_radius; double temp_outer_radius{}; + double temp_inner_sphericity = inner_sphericity; for (size_t layer_i = 0; layer_i < number_of_layers; layer_i++) { const auto& radial_distribution_this_layer = radial_distribution.at(layer_i); @@ -640,29 +638,29 @@ std::vector> sph_wedge_coordinate_maps( for (size_t face_j = which_wedge_index(which_wedges); face_j < 6; face_j++) { wedges_for_this_layer.emplace_back( - temp_inner_radius, temp_outer_radius, inner_sphericity, + temp_inner_radius, temp_outer_radius, temp_inner_sphericity, outer_sphericity, gsl::at(wedge_orientations, face_j), use_equiangular_map, Halves::Both, radial_distribution_this_layer); } } else { for (size_t i = 0; i < 4; i++) { wedges_for_this_layer.emplace_back( - temp_inner_radius, temp_outer_radius, inner_sphericity, + temp_inner_radius, temp_outer_radius, temp_inner_sphericity, outer_sphericity, gsl::at(wedge_orientations, i), use_equiangular_map, Halves::LowerOnly, radial_distribution_this_layer); wedges_for_this_layer.emplace_back( - temp_inner_radius, temp_outer_radius, inner_sphericity, + temp_inner_radius, temp_outer_radius, temp_inner_sphericity, outer_sphericity, gsl::at(wedge_orientations, i), use_equiangular_map, Halves::UpperOnly, radial_distribution_this_layer); } wedges_for_this_layer.emplace_back( - temp_inner_radius, temp_outer_radius, inner_sphericity, + temp_inner_radius, temp_outer_radius, temp_inner_sphericity, outer_sphericity, gsl::at(wedge_orientations, 4), use_equiangular_map, Halves::Both, radial_distribution_this_layer); wedges_for_this_layer.emplace_back( - temp_inner_radius, temp_outer_radius, inner_sphericity, + temp_inner_radius, temp_outer_radius, temp_inner_sphericity, outer_sphericity, gsl::at(wedge_orientations, 5), use_equiangular_map, Halves::Both, radial_distribution_this_layer); } @@ -672,6 +670,7 @@ std::vector> sph_wedge_coordinate_maps( if (layer_i != radial_partitioning.size()) { temp_inner_radius = radial_partitioning.at(layer_i); + temp_inner_sphericity = outer_sphericity; } } return wedges_for_all_layers; @@ -1422,7 +1421,6 @@ template class domain::CoordinateMaps::ProductOf3Maps< domain::CoordinateMaps::Affine>; template class domain::CoordinateMaps::ProductOf2Maps< domain::CoordinateMaps::Equiangular, domain::CoordinateMaps::Equiangular>; - template class domain::CoordinateMaps::ProductOf3Maps< domain::CoordinateMaps::Equiangular, domain::CoordinateMaps::Equiangular, domain::CoordinateMaps::Equiangular>; diff --git a/src/Domain/DomainHelpers.hpp b/src/Domain/DomainHelpers.hpp index 4de0e961b53a..98f916397beb 100644 --- a/src/Domain/DomainHelpers.hpp +++ b/src/Domain/DomainHelpers.hpp @@ -133,23 +133,30 @@ enum class ShellWedges { OneAlongMinusX }; -/// \ingroup ComputationalDomainGroup -/// These are the CoordinateMaps of the Wedge<3>s used in the Sphere, Shell, and -/// binary compact object DomainCreators. This function can also be used to -/// wrap the Sphere or Shell in a cube made of six Wedge<3>s. -/// The argument `x_coord_of_shell_center` specifies a translation of the Shell -/// in the x-direction in the TargetFrame. For example, the BBH DomainCreator -/// uses this to set the position of each BH. -/// When the argument `use_half_wedges` is set to `true`, the wedges in the -/// +z,-z,+y,-y directions are cut in half along their xi-axes. The resulting -/// ten CoordinateMaps are used for the outermost Blocks of the BBH Domain. -/// The argument `aspect_ratio` sets the equatorial compression factor, -/// used by the EquatorialCompression maps which get composed with the Wedges. -/// This is done if `aspect_ratio` is set to something other than the default -/// value of one. The `radial_partitioning` specifies the radial boundaries of -/// sub-shells between `inner_radius` and `outer_radius`. Set the -/// `radial_distribution` to select the radial distribution of grid points in -/// the spherical shells. +/*! + * \ingroup ComputationalDomainGroup + * These are the CoordinateMaps of the Wedge<3>s used in the Sphere, Shell, and + * binary compact object DomainCreators. This function can also be used to + * wrap the Sphere or Shell in a cube made of six Wedge<3>s. + * + * \param inner_radius Radius of the inner boundary of the shell, or the + * radius circumscribing the inner cube of a sphere. + * \param outer_radius Outer radius of the shell or sphere. + * \param inner_sphericity Specifies if the wedges form a spherical inner + * boundary (1.0) or a cubical inner boundary (0.0). + * \param outer_sphericity Specifies if the wedges form a spherical outer + * boundary (1.0) or a cubical outer boundary (0.0). + * \param use_equiangular_map Toggles the equiangular map of the Wedge map. + * \param use_half_wedges When `true`, the wedges in the +z,-z,+y,-y directions + * are cut in half along their xi-axes. The resulting ten CoordinateMaps are + * used for the outermost Blocks of the BBH Domain. + * \param radial_partitioning Specifies the radial boundaries of sub-shells + * between `inner_radius` and `outer_radius`. If the inner and outer + * sphericities are different, the innermost shell does the transition. + * \param radial_distribution Select the radial distribution of grid points in + * the spherical shells. + * \param which_wedges Select a subset of wedges. + */ std::vector> sph_wedge_coordinate_maps( double inner_radius, double outer_radius, double inner_sphericity, double outer_sphericity, bool use_equiangular_map, diff --git a/tests/InputFiles/Punctures/MultiplePunctures.yaml b/tests/InputFiles/Punctures/MultiplePunctures.yaml index 9045a291666e..330ea91273ee 100644 --- a/tests/InputFiles/Punctures/MultiplePunctures.yaml +++ b/tests/InputFiles/Punctures/MultiplePunctures.yaml @@ -28,7 +28,9 @@ DomainCreator: OuterRadius: 10. InnerCubeSphericity: 0. InitialRefinement: 0 - InitialGridPoints: [4, 4] + InitialGridPoints: 4 + RadialPartitioning: [] + RadialDistribution: [Linear] UseEquiangularMap: True TimeDependence: None BoundaryCondition: Flatness diff --git a/tests/InputFiles/Xcts/TovStar.yaml b/tests/InputFiles/Xcts/TovStar.yaml index 83b9462f3efd..4825c794d08b 100644 --- a/tests/InputFiles/Xcts/TovStar.yaml +++ b/tests/InputFiles/Xcts/TovStar.yaml @@ -35,7 +35,9 @@ DomainCreator: OuterRadius: 1.24984447898 InnerCubeSphericity: 0. InitialRefinement: 0 - InitialGridPoints: [5, 5] + InitialGridPoints: [5, 5, 5] + RadialPartitioning: [] + RadialDistribution: [Linear] UseEquiangularMap: True TimeDependence: None BoundaryCondition: diff --git a/tests/Unit/Domain/Creators/Python/Test_Sphere.py b/tests/Unit/Domain/Creators/Python/Test_Sphere.py index 74995b8141bb..7bdde4252c74 100644 --- a/tests/Unit/Domain/Creators/Python/Test_Sphere.py +++ b/tests/Unit/Domain/Creators/Python/Test_Sphere.py @@ -11,7 +11,7 @@ def test_construction(self): outer_radius=2., inner_cube_sphericity=0., initial_refinement=1, - initial_number_of_grid_points=[3, 3], + initial_number_of_grid_points=[3, 3, 3], use_equiangular_map=False) self.assertIsInstance(sphere, DomainCreator3D) diff --git a/tests/Unit/Domain/Creators/Test_Sphere.cpp b/tests/Unit/Domain/Creators/Test_Sphere.cpp index d797610bb843..71604cce8cc4 100644 --- a/tests/Unit/Domain/Creators/Test_Sphere.cpp +++ b/tests/Unit/Domain/Creators/Test_Sphere.cpp @@ -9,10 +9,13 @@ #include #include #include +#include #include #include #include +#include "DataStructures/Tensor/EagerMath/DotProduct.hpp" +#include "DataStructures/Tensor/EagerMath/Magnitude.hpp" #include "DataStructures/Tensor/Tensor.hpp" #include "Domain/Block.hpp" #include "Domain/BoundaryConditions/BoundaryCondition.hpp" @@ -45,6 +48,7 @@ #include "Helpers/Domain/Creators/TestHelpers.hpp" #include "Helpers/Domain/DomainTestHelpers.hpp" #include "Parallel/RegisterDerivedClassesWithCharm.hpp" +#include "Utilities/CartesianProduct.hpp" #include "Utilities/CloneUniquePtrs.hpp" #include "Utilities/MakeArray.hpp" @@ -59,288 +63,212 @@ create_boundary_condition() { Direction<3>::upper_zeta(), 50); } -auto make_domain_creator(const std::string& opt_string, - const bool use_boundary_conditions) { - if (use_boundary_conditions) { - return TestHelpers::test_option_tag< - domain::OptionTags::DomainCreator<3>, - TestHelpers::domain::BoundaryConditions:: - MetavariablesWithBoundaryConditions<3, domain::creators::Sphere>>( - opt_string + std::string{" BoundaryCondition:\n" - " TestBoundaryCondition:\n" - " Direction: upper-zeta\n" - " BlockId: 50\n"}); +// Calculate block logical coordinates residing on corners of +// inner block or block faces of wedges +// whose normal vector point radially outward from the +// origin. With this domain, this direction corresponds to +// upper zeta. These coordinates will be used to ensure +// they lie on concentric spheres defined by either the inner +// sphere, outer sphere, or radial partition parameters. +tnsr::I logical_coords( + const gsl::not_null generator, const size_t num_blocks, + const size_t block_id, const bool abuts_inner_block) { + std::uniform_real_distribution<> real_dis(-1, 1); + + const double rand_int_xi = (2.0 * (rand() % 2) - 1.0); + const double rand_int_eta = (2.0 * (rand() % 2) - 1.0); + const double rand_int_zeta = (2.0 * (rand() % 2) - 1.0); + const double rand_real_xi = real_dis(*generator); + const double rand_real_eta = real_dis(*generator); + + double xi_logical_coord; + double eta_logical_coord; + // enforce coordinates either fall on the lower or + // upper zeta face of wedges + const double zeta_logical_coord = rand_int_zeta; + + if (block_id == num_blocks - 1) { + // inner block only uses integer corners + xi_logical_coord = rand_int_xi; + eta_logical_coord = rand_int_eta; + + } else if (abuts_inner_block) { + // next to inner block, + // corners only on lower face b/c of square inner block neighbor + // face + xi_logical_coord = rand_int_xi; + eta_logical_coord = rand_int_eta; + + // anywhere on upper zeta face b/c adjacent with spherical wedge + if (rand_int_zeta == 1) { + xi_logical_coord = rand_real_xi; + eta_logical_coord = rand_real_eta; + } } else { - return TestHelpers::test_option_tag< - domain::OptionTags::DomainCreator<3>, - TestHelpers::domain::BoundaryConditions:: - MetavariablesWithoutBoundaryConditions<3, - domain::creators::Sphere>>( - opt_string); + // adjacent to wedges + // everywhere on low or high face should lie on a sphere b/c + // neighbor with spherical wedge + xi_logical_coord = rand_real_xi; + eta_logical_coord = rand_real_eta; } + + return tnsr::I{ + {{xi_logical_coord, eta_logical_coord, zeta_logical_coord}}}; } template void test_sphere_construction( const creators::Sphere& sphere, const double inner_radius, - const double outer_radius, const double inner_cube_sphericity, - const bool use_equiangular_map, - const std::array& expected_sphere_extents, - const std::vector>& expected_refinement_level, - const bool expect_boundary_conditions = false, - const std::tuple...>& - expected_functions_of_time = {}, - const std::vector>>& expected_grid_to_inertial_maps = {}, - const std::vector< - std::unique_ptr>>& - time_dependencies = {}, - const std::unordered_map& initial_expiration_times = - {}) { - CAPTURE(inner_radius); - CAPTURE(outer_radius); - CAPTURE(inner_cube_sphericity); - CAPTURE(use_equiangular_map); - CAPTURE(expect_boundary_conditions); + const double outer_radius, + const std::vector radial_partitioning = {}, + const bool expect_boundary_conditions = true, + const double current_time = 1.0, + const std::array velocity = {}) { + // check consistency of domain const auto domain = TestHelpers::domain::creators::test_domain_creator( sphere, expect_boundary_conditions); - const OrientationMap<3> aligned_orientation{}; - const OrientationMap<3> quarter_turn_ccw_about_zeta( - std::array, 3>{{Direction<3>::lower_eta(), - Direction<3>::upper_xi(), - Direction<3>::upper_zeta()}}); - const OrientationMap<3> half_turn_about_zeta(std::array, 3>{ - {Direction<3>::lower_xi(), Direction<3>::lower_eta(), - Direction<3>::upper_zeta()}}); - const OrientationMap<3> quarter_turn_cw_about_zeta( - std::array, 3>{{Direction<3>::upper_eta(), - Direction<3>::lower_xi(), - Direction<3>::upper_zeta()}}); - const OrientationMap<3> center_relative_to_minus_z( - std::array, 3>{{Direction<3>::upper_xi(), - Direction<3>::lower_eta(), - Direction<3>::lower_zeta()}}); - const OrientationMap<3> center_relative_to_plus_y(std::array, 3>{ - {Direction<3>::upper_xi(), Direction<3>::lower_zeta(), - Direction<3>::upper_eta()}}); - const OrientationMap<3> center_relative_to_minus_y( - std::array, 3>{{Direction<3>::upper_xi(), - Direction<3>::upper_zeta(), - Direction<3>::lower_eta()}}); - const OrientationMap<3> center_relative_to_plus_x(std::array, 3>{ - {Direction<3>::upper_eta(), Direction<3>::upper_zeta(), - Direction<3>::upper_xi()}}); - const OrientationMap<3> center_relative_to_minus_x( - std::array, 3>{{Direction<3>::lower_eta(), - Direction<3>::upper_zeta(), - Direction<3>::lower_xi()}}); - - const std::vector>> expected_block_neighbors{ - {{Direction<3>::upper_xi(), {4, quarter_turn_ccw_about_zeta}}, - {Direction<3>::upper_eta(), {2, aligned_orientation}}, - {Direction<3>::lower_xi(), {5, quarter_turn_cw_about_zeta}}, - {Direction<3>::lower_eta(), {3, aligned_orientation}}, - {Direction<3>::lower_zeta(), {6, aligned_orientation}}}, - {{Direction<3>::upper_xi(), {4, quarter_turn_cw_about_zeta}}, - {Direction<3>::upper_eta(), {3, aligned_orientation}}, - {Direction<3>::lower_xi(), {5, quarter_turn_ccw_about_zeta}}, - {Direction<3>::lower_eta(), {2, aligned_orientation}}, - {Direction<3>::lower_zeta(), {6, center_relative_to_minus_z}}}, - {{Direction<3>::upper_xi(), {4, half_turn_about_zeta}}, - {Direction<3>::upper_eta(), {1, aligned_orientation}}, - {Direction<3>::lower_xi(), {5, half_turn_about_zeta}}, - {Direction<3>::lower_eta(), {0, aligned_orientation}}, - {Direction<3>::lower_zeta(), {6, center_relative_to_plus_y}}}, - {{Direction<3>::upper_xi(), {4, aligned_orientation}}, - {Direction<3>::upper_eta(), {0, aligned_orientation}}, - {Direction<3>::lower_xi(), {5, aligned_orientation}}, - {Direction<3>::lower_eta(), {1, aligned_orientation}}, - {Direction<3>::lower_zeta(), {6, center_relative_to_minus_y}}}, - {{Direction<3>::upper_xi(), {2, half_turn_about_zeta}}, - {Direction<3>::upper_eta(), {0, quarter_turn_cw_about_zeta}}, - {Direction<3>::lower_xi(), {3, aligned_orientation}}, - {Direction<3>::lower_eta(), {1, quarter_turn_ccw_about_zeta}}, - {Direction<3>::lower_zeta(), {6, center_relative_to_plus_x}}}, - {{Direction<3>::upper_xi(), {3, aligned_orientation}}, - {Direction<3>::upper_eta(), {0, quarter_turn_ccw_about_zeta}}, - {Direction<3>::lower_xi(), {2, half_turn_about_zeta}}, - {Direction<3>::lower_eta(), {1, quarter_turn_cw_about_zeta}}, - {Direction<3>::lower_zeta(), {6, center_relative_to_minus_x}}}, - {{Direction<3>::upper_zeta(), {0, aligned_orientation}}, - {Direction<3>::lower_zeta(), - {1, center_relative_to_minus_z.inverse_map()}}, - {Direction<3>::upper_eta(), - {2, center_relative_to_plus_y.inverse_map()}}, - {Direction<3>::lower_eta(), - {3, center_relative_to_minus_y.inverse_map()}}, - {Direction<3>::upper_xi(), {4, center_relative_to_plus_x.inverse_map()}}, - {Direction<3>::lower_xi(), - {5, center_relative_to_minus_x.inverse_map()}}}}; - - const std::vector>> - expected_external_boundaries{{{Direction<3>::upper_zeta()}}, - {{Direction<3>::upper_zeta()}}, - {{Direction<3>::upper_zeta()}}, - {{Direction<3>::upper_zeta()}}, - {{Direction<3>::upper_zeta()}}, - {{Direction<3>::upper_zeta()}}, - {}}; - - std::vector> expected_extents{ - 6, - {{expected_sphere_extents[1], expected_sphere_extents[1], - expected_sphere_extents[0]}}}; - expected_extents.push_back( - {{expected_sphere_extents[1], expected_sphere_extents[1], - expected_sphere_extents[1]}}); - - CHECK(sphere.initial_extents() == expected_extents); - CHECK(sphere.initial_refinement_levels() == expected_refinement_level); - using Wedge3DMap = CoordinateMaps::Wedge<3>; - using Affine = CoordinateMaps::Affine; - using Affine3D = CoordinateMaps::ProductOf3Maps; - using Equiangular = CoordinateMaps::Equiangular; - using Equiangular3D = - CoordinateMaps::ProductOf3Maps; - using BulgedCube = CoordinateMaps::BulgedCube; - - const auto make_coord_maps = [&inner_radius, &outer_radius, - &inner_cube_sphericity, - &use_equiangular_map](const auto frame) { - using TargetFrame = std::decay_t; - auto local_coord_maps = make_vector_coordinate_map_base( - Wedge3DMap{inner_radius, outer_radius, inner_cube_sphericity, 1.0, - OrientationMap<3>{}, use_equiangular_map}, - Wedge3DMap{inner_radius, outer_radius, inner_cube_sphericity, 1.0, - OrientationMap<3>{std::array, 3>{ - {Direction<3>::upper_xi(), Direction<3>::lower_eta(), - Direction<3>::lower_zeta()}}}, - use_equiangular_map}, - Wedge3DMap{inner_radius, outer_radius, inner_cube_sphericity, 1.0, - OrientationMap<3>{std::array, 3>{ - {Direction<3>::upper_xi(), Direction<3>::upper_zeta(), - Direction<3>::lower_eta()}}}, - use_equiangular_map}, - Wedge3DMap{inner_radius, outer_radius, inner_cube_sphericity, 1.0, - OrientationMap<3>{std::array, 3>{ - {Direction<3>::upper_xi(), Direction<3>::lower_zeta(), - Direction<3>::upper_eta()}}}, - use_equiangular_map}, - Wedge3DMap{inner_radius, outer_radius, inner_cube_sphericity, 1.0, - OrientationMap<3>{std::array, 3>{ - {Direction<3>::upper_zeta(), Direction<3>::upper_xi(), - Direction<3>::upper_eta()}}}, - use_equiangular_map}, - Wedge3DMap{inner_radius, outer_radius, inner_cube_sphericity, 1.0, - OrientationMap<3>{std::array, 3>{ - {Direction<3>::lower_zeta(), Direction<3>::lower_xi(), - Direction<3>::upper_eta()}}}, - use_equiangular_map}); - if (inner_cube_sphericity == 0.0) { - if (use_equiangular_map) { - local_coord_maps.emplace_back( - make_coordinate_map_base( - Equiangular3D{ - Equiangular(-1.0, 1.0, -1.0 * inner_radius / sqrt(3.0), - inner_radius / sqrt(3.0)), - Equiangular(-1.0, 1.0, -1.0 * inner_radius / sqrt(3.0), - inner_radius / sqrt(3.0)), - Equiangular(-1.0, 1.0, -1.0 * inner_radius / sqrt(3.0), - inner_radius / sqrt(3.0))})); + const auto& blocks = domain.blocks(); + const auto block_names = sphere.block_names(); + const size_t num_blocks = blocks.size(); + const auto all_boundary_conditions = sphere.external_boundary_conditions(); + + MAKE_GENERATOR(generator); + + // verify if adjacent to inner block + const auto abuts_inner_block = + [&num_blocks](const auto& direction_and_neighbor) { + return direction_and_neighbor.second.id() == num_blocks - 1; + }; + + for (size_t block_id = 0; block_id < num_blocks; ++block_id) { + const auto& block = blocks[block_id]; + const auto& boundary_conditions = all_boundary_conditions[block_id]; + + // This section tests if the logical coordinates of corners from all blocks + // (and points on upper wedge faces) lie on spherical shells specified + // by inner radius, radial partitions, or outer radius + const auto coords_on_spherical_partition = + logical_coords(make_not_null(&generator), num_blocks, block_id, + alg::any_of(block.neighbors(), abuts_inner_block)); + + const double corner_distance_from_origin = + [&block, &coords_on_spherical_partition, ¤t_time, &sphere, + &velocity]() -> double { + // use stationary map if independent of time + if (not block.is_time_dependent()) { + auto location_time_indep = + block.stationary_map()(coords_on_spherical_partition); + return get(magnitude(location_time_indep)); } else { - local_coord_maps.emplace_back( - make_coordinate_map_base( - Affine3D{Affine(-1.0, 1.0, -1.0 * inner_radius / sqrt(3.0), - inner_radius / sqrt(3.0)), - Affine(-1.0, 1.0, -1.0 * inner_radius / sqrt(3.0), - inner_radius / sqrt(3.0)), - Affine(-1.0, 1.0, -1.0 * inner_radius / sqrt(3.0), - inner_radius / sqrt(3.0))})); + // go from logical to grid coords, then grid to inertial coords + auto inertial_location_time_dep = + block.moving_mesh_grid_to_inertial_map()( + block.moving_mesh_logical_to_grid_map()( + coords_on_spherical_partition, current_time, + sphere.functions_of_time()), + current_time, sphere.functions_of_time()); + + // origin in inertial frame (need to shift inertial + // coord by velocity * (final_time - initial_time)) + return sqrt(square(get<0>(inertial_location_time_dep) - + velocity[0] * (current_time - 1.0)) + + square(get<1>(inertial_location_time_dep) - + velocity[1] * (current_time - 1.0)) + + square(get<2>(inertial_location_time_dep) - + velocity[2] * (current_time - 1.0))); + } // end time-dependent if/else + }(); // end lambda + + // construct vector of inner radius, outer radius, and refinements levels + // where inertial block corners have to be located + std::vector expected_corner_radii = radial_partitioning; + expected_corner_radii.insert(expected_corner_radii.begin(), inner_radius); + expected_corner_radii.emplace_back(outer_radius); + + const auto match_demarcation = + [&corner_distance_from_origin](const double radius) { + return corner_distance_from_origin == approx(radius); + }; + + CHECK(alg::any_of(expected_corner_radii, match_demarcation)); + + // if block has 5 neighbors, 1 face should be external, and that direction + // should be upper_zeta, for the sphere + if (block.neighbors().size() == 5) { + CHECK(size(block.external_boundaries()) == 1); + CHECK(*begin(block.external_boundaries()) == Direction<3>::upper_zeta()); + + // also 5 neighbor blocks should only have 1 boundary condition + if (expect_boundary_conditions) { + CHECK(size(boundary_conditions) == 1); } - } else { - local_coord_maps.emplace_back( - make_coordinate_map_base(BulgedCube{ - inner_radius, inner_cube_sphericity, use_equiangular_map})); - } - return local_coord_maps; - }; - - auto coord_maps = - make_coord_maps(tmpl::conditional_t{}); - test_domain_construction(domain, expected_block_neighbors, - expected_external_boundaries, coord_maps, 10.0, - sphere.functions_of_time(), - expected_grid_to_inertial_maps); - const auto coord_maps_copy = clone_unique_ptrs(coord_maps); - - Domain<3> domain_no_corners = Domain<3>{make_coord_maps(Frame::Inertial{})}; - - if (sizeof...(FuncsOfTime) != 0) { - for (const auto& time_dependence : time_dependencies) { - const size_t number_of_blocks = domain_no_corners.blocks().size(); - auto block_maps_grid_to_inertial = - time_dependence->block_maps_grid_to_inertial(number_of_blocks); - auto block_maps_grid_to_distorted = - time_dependence->block_maps_grid_to_distorted(number_of_blocks); - auto block_maps_distorted_to_inertial = - time_dependence->block_maps_distorted_to_inertial(number_of_blocks); - for (size_t block_id = 0; block_id < number_of_blocks; ++block_id) { - domain_no_corners.inject_time_dependent_map_for_block( - block_id, std::move(block_maps_grid_to_inertial[block_id]), - std::move(block_maps_grid_to_distorted[block_id]), - std::move(block_maps_distorted_to_inertial[block_id])); + // Consistency check for like neighbor block directions: if block has 5 + // neighbors, it's external --> four of the neighbors should have upper + // zeta external boundaries + size_t neighbor_count = 0; + for (auto neighbor : block.neighbors()) { + auto neighbor_id = neighbor.second; + + if (size(blocks[neighbor_id.id()].external_boundaries()) == 1) { + if (*begin(blocks[neighbor_id.id()].external_boundaries()) == + Direction<3>::upper_zeta()) { + neighbor_count++; + } + } + } + + CHECK(neighbor_count == 4); + // if > 5 neighbors, none should have external boundaries + } else if (block.neighbors().size() == 6) { + // internal block case + CHECK(size(block.external_boundaries()) == 0); + // internal blocks should not have boundary conditions + if (expect_boundary_conditions) { + CHECK(size(boundary_conditions) == 0); } + } else { + // If here, something is wrong; should only have 5 or 6 neighbors, so + // throw a guaranteed failure. + const bool block_does_not_have_correct_number_of_neighbors = false; + CHECK(block_does_not_have_correct_number_of_neighbors); } - } + } // block loop - test_domain_construction(domain_no_corners, expected_block_neighbors, - expected_external_boundaries, coord_maps_copy, 10.0, - sphere.functions_of_time(), - expected_grid_to_inertial_maps); - test_initial_domain(domain_no_corners, sphere.initial_refinement_levels()); - test_serialization(domain_no_corners); + // verify there are no gaps between blocks + test_physical_separation(domain.blocks(), current_time, + sphere.functions_of_time()); - TestHelpers::domain::creators::test_functions_of_time( - sphere, expected_functions_of_time, initial_expiration_times); -} +} // test_sphere_construction() -void test_sphere_boundaries_equiangular() { - INFO("Sphere boundaries equiangular"); +// ensure CHECK_THROWS_WITH calls are properly captured +void test_parse_errors() { + INFO("Sphere check throws"); const double inner_radius = 1.0; const double outer_radius = 2.0; const size_t refinement = 2; - const std::array grid_points_r_angular{{4, 4}}; - - for (const auto sphericity : {0.0, 0.2, 0.7}) { - CAPTURE(sphericity); - const creators::Sphere sphere{ - inner_radius, outer_radius, sphericity, - refinement, grid_points_r_angular, true}; - test_sphere_construction(sphere, inner_radius, outer_radius, sphericity, - true, grid_points_r_angular, - {7, make_array<3>(refinement)}); - - const creators::Sphere sphere_boundary_condition{ - inner_radius, - outer_radius, - sphericity, - refinement, - grid_points_r_angular, - true, - nullptr, - create_boundary_condition()}; - test_sphere_construction( - sphere_boundary_condition, inner_radius, outer_radius, sphericity, true, - grid_points_r_angular, {7, make_array<3>(refinement)}, true); - } + const std::array initial_extents{{4, 5, 6}}; + const std::vector radial_partitioning = {}; + const std::vector radial_partitioning_unordered = { + {1.5 * inner_radius, 1.1 * inner_radius}}; + const std::vector radial_partitioning_low = { + {0.5 * inner_radius, 1.1 * inner_radius}}; + const std::vector radial_partitioning_high = { + {2.1 * outer_radius, 2.2 * outer_radius}}; + const std::vector radial_distribution{ + domain::CoordinateMaps::Distribution::Linear}; + const std::vector + radial_distribution_too_many{ + domain::CoordinateMaps::Distribution::Linear, + domain::CoordinateMaps::Distribution::Logarithmic}; + const std::vector + radial_distribution_inner_log{ + domain::CoordinateMaps::Distribution::Logarithmic}; CHECK_THROWS_WITH( creators::Sphere( - inner_radius, outer_radius, -1.0, refinement, grid_points_r_angular, - true, nullptr, + inner_radius, outer_radius, -1.0, refinement, initial_extents, true, + radial_partitioning, radial_distribution, nullptr, std::make_unique>(), Options::Context{false, {}, 1, 1}), @@ -348,8 +276,8 @@ void test_sphere_boundaries_equiangular() { "Inner cube sphericity must be >= 0.0 and strictly < 1.0")); CHECK_THROWS_WITH( creators::Sphere( - inner_radius, outer_radius, 1.0, refinement, grid_points_r_angular, - true, nullptr, + inner_radius, outer_radius, 1.0, refinement, initial_extents, true, + radial_partitioning, radial_distribution, nullptr, std::make_unique>(), Options::Context{false, {}, 1, 1}), @@ -357,17 +285,75 @@ void test_sphere_boundaries_equiangular() { "Inner cube sphericity must be >= 0.0 and strictly < 1.0")); CHECK_THROWS_WITH( creators::Sphere( - inner_radius, outer_radius, 2.0, refinement, grid_points_r_angular, - true, nullptr, + inner_radius, outer_radius, 2.0, refinement, initial_extents, true, + radial_partitioning, radial_distribution, nullptr, std::make_unique>(), Options::Context{false, {}, 1, 1}), Catch::Matchers::Contains( "Inner cube sphericity must be >= 0.0 and strictly < 1.0")); + + CHECK_THROWS_WITH( + creators::Sphere( + inner_radius, 0.5 * inner_radius, 0.5, refinement, initial_extents, + true, radial_partitioning, radial_distribution, nullptr, + std::make_unique>(), + Options::Context{false, {}, 1, 1}), + Catch::Matchers::Contains( + "Inner radius must be smaller than outer radius")); + CHECK_THROWS_WITH( creators::Sphere( - inner_radius, outer_radius, 0.0, refinement, grid_points_r_angular, - true, nullptr, + inner_radius, outer_radius, 0.5, refinement, initial_extents, true, + radial_partitioning_unordered, radial_distribution, nullptr, + std::make_unique>(), + Options::Context{false, {}, 1, 1}), + Catch::Matchers::Contains( + "Specify radial partitioning in ascending order.")); + + CHECK_THROWS_WITH( + creators::Sphere( + inner_radius, outer_radius, 0.5, refinement, initial_extents, true, + radial_partitioning_low, radial_distribution, nullptr, + std::make_unique>(), + Options::Context{false, {}, 1, 1}), + Catch::Matchers::Contains( + "First radial partition must be larger than inner")); + CHECK_THROWS_WITH( + creators::Sphere( + inner_radius, outer_radius, 0.5, refinement, initial_extents, true, + radial_partitioning_high, radial_distribution, nullptr, + std::make_unique>(), + Options::Context{false, {}, 1, 1}), + Catch::Matchers::Contains( + "Last radial partition must be smaller than outer")); + CHECK_THROWS_WITH( + creators::Sphere( + inner_radius, outer_radius, 0.5, refinement, initial_extents, true, + radial_partitioning, radial_distribution_too_many, nullptr, + std::make_unique>(), + Options::Context{false, {}, 1, 1}), + Catch::Matchers::Contains( + "Specify a 'RadialDistribution' for every spherical shell. You")); + CHECK_THROWS_WITH( + creators::Sphere( + inner_radius, outer_radius, 0.5, refinement, initial_extents, true, + radial_partitioning, radial_distribution_inner_log, nullptr, + std::make_unique>(), + Options::Context{false, {}, 1, 1}), + Catch::Matchers::Contains( + "The 'RadialDistribution' must be 'Linear' for the")); + + CHECK_THROWS_WITH( + creators::Sphere( + inner_radius, outer_radius, 0.0, refinement, initial_extents, true, + radial_partitioning, radial_distribution, nullptr, std::make_unique>(), Options::Context{false, {}, 1, 1}), @@ -375,8 +361,8 @@ void test_sphere_boundaries_equiangular() { "Cannot have periodic boundary conditions with a Sphere")); CHECK_THROWS_WITH( creators::Sphere( - inner_radius, outer_radius, 0.0, refinement, grid_points_r_angular, - true, nullptr, + inner_radius, outer_radius, 0.0, refinement, initial_extents, true, + radial_partitioning, radial_distribution, nullptr, std::make_unique>(), Options::Context{false, {}, 1, 1}), @@ -385,163 +371,105 @@ void test_sphere_boundaries_equiangular() { "an outflow-type boundary condition, you must use that.")); } -void test_sphere_factory_equiangular() { - INFO("Sphere factory equiangular"); - for (bool use_boundary_conditions : {false, true}) { - const auto sphere = make_domain_creator( - "Sphere:\n" - " InnerRadius: 1\n" - " OuterRadius: 3\n" - " InnerCubeSphericity: 0.0\n" - " InitialRefinement: 2\n" - " InitialGridPoints: [2,3]\n" - " UseEquiangularMap: true\n" - " TimeDependence: None\n", - use_boundary_conditions); - const double inner_radius = 1.0; - const double outer_radius = 3.0; - const size_t refinement_level = 2; - const std::array grid_points_r_angular{{2, 3}}; - test_sphere_construction( - dynamic_cast(*sphere), inner_radius, - outer_radius, 0.0, true, grid_points_r_angular, - {7, make_array<3>(refinement_level)}, use_boundary_conditions); - } -} - -void test_sphere_boundaries_equidistant() { - INFO("Sphere boundaries equidistant"); +// Check wedge neighbors have consistent directions & orientations +void test_sphere_boundaries() { + INFO("Ensure sphere boundaries are equidistant"); const double inner_radius = 1.0; const double outer_radius = 2.0; - const size_t refinement = 2; - const std::array grid_points_r_angular{{4, 4}}; + const size_t initial_refinement = 3; + const std::array initial_extents{{4, 5, 6}}; + const std::array, 2> radial_partitioning{ + {{}, {0.5 * (inner_radius + outer_radius)}}}; + const std::array, 2> + radial_distribution{ + {{domain::CoordinateMaps::Distribution::Linear}, + {domain::CoordinateMaps::Distribution::Linear, + domain::CoordinateMaps::Distribution::Logarithmic}}}; - for (const auto sphericity : {0.0, 0.2, 0.7}) { + for (const auto& [sphericity, equiangular, array_index] : + cartesian_product(make_array(0.0, 0.7), make_array(false, true), + make_array(0.0, 1.0))) { + CAPTURE(inner_radius); + CAPTURE(outer_radius); CAPTURE(sphericity); - const creators::Sphere sphere{ - inner_radius, outer_radius, sphericity, - refinement, grid_points_r_angular, false}; - test_sphere_construction(sphere, inner_radius, outer_radius, sphericity, - false, grid_points_r_angular, - {7, make_array<3>(refinement)}); + CAPTURE(initial_refinement); + CAPTURE(initial_extents); + CAPTURE(equiangular); + CAPTURE(radial_partitioning[array_index]); + CAPTURE(radial_distribution[array_index]); + + const creators::Sphere sphere{inner_radius, + outer_radius, + sphericity, + initial_refinement, + initial_extents, + equiangular, + radial_partitioning[array_index], + radial_distribution[array_index]}; + + test_sphere_construction(sphere, inner_radius, outer_radius, + radial_partitioning[array_index], false); const creators::Sphere sphere_boundary_condition{ inner_radius, outer_radius, sphericity, - refinement, - grid_points_r_angular, - false, + initial_refinement, + initial_extents, + equiangular, + radial_partitioning[array_index], + radial_distribution[array_index], nullptr, create_boundary_condition()}; - test_sphere_construction( - sphere_boundary_condition, inner_radius, outer_radius, sphericity, - false, grid_points_r_angular, {7, make_array<3>(refinement)}, true); - } -} -void test_sphere_factory_equidistant() { - INFO("Sphere factory equidistant"); - for (bool use_boundary_conditions : {false, true}) { - const auto sphere = make_domain_creator( - "Sphere:\n" - " InnerRadius: 1\n" - " OuterRadius: 3\n" - " InnerCubeSphericity: 0.1\n" - " InitialRefinement: 2\n" - " InitialGridPoints: [2,3]\n" - " UseEquiangularMap: false\n" - " TimeDependence: None\n", - use_boundary_conditions); - const double inner_radius = 1.0; - const double outer_radius = 3.0; - const double sphericity = 0.1; - const size_t refinement_level = 2; - const std::array grid_points_r_angular{{2, 3}}; - test_sphere_construction( - dynamic_cast(*sphere), inner_radius, - outer_radius, sphericity, false, grid_points_r_angular, - {7, make_array<3>(refinement_level)}, use_boundary_conditions); + test_sphere_construction(sphere_boundary_condition, inner_radius, + outer_radius, radial_partitioning[array_index], + true); } } +// Check wedge neighbors have consistent directions & orientations, with time +// dependence void test_sphere_factory_time_dependent() { INFO("Sphere factory time dependent"); - // This name must match the hard coded one in UniformTranslation - const std::string f_of_t_name = "Translation"; - const auto helper = [&f_of_t_name]( - const std::unordered_map& - expiration_times) { - const auto domain_creator = - TestHelpers::test_option_tag, - TestHelpers::domain::BoundaryConditions:: - MetavariablesWithoutBoundaryConditions< - 3, domain::creators::Sphere>>( - "Sphere:\n" - " InnerRadius: 1\n" - " OuterRadius: 3\n" - " InnerCubeSphericity: 0.0\n" - " InitialRefinement: 2\n" - " InitialGridPoints: [2,3]\n" - " UseEquiangularMap: false\n" - " TimeDependence:\n" - " UniformTranslation:\n" - " InitialTime: 1.0\n" - " Velocity: [2.3, -0.3, 0.5]\n"); - const auto* sphere = - dynamic_cast(domain_creator.get()); - const double initial_time = 1.0; - const double inner_radius = 1.0; - const double outer_radius = 3.0; - const double sphericity = 0.0; - const size_t refinement_level = 2; - const std::array grid_points_r_angular{{2, 3}}; - const bool use_equiangular_map = false; - const DataVector velocity{{2.3, -0.3, 0.5}}; - std::vector< - std::unique_ptr>> - time_dependencies{1}; - time_dependencies[0] = std::make_unique< - domain::creators::time_dependence::UniformTranslation<3>>( - initial_time, std::array{velocity[0], velocity[1], velocity[2]}); - auto grid_to_inertial_maps = - make_vector_coordinate_map_base( - Translation3D{f_of_t_name}); - for (size_t i = 1; i < sphere->create_domain().blocks().size(); i++) { - grid_to_inertial_maps.emplace_back( - make_coordinate_map_base( - Translation3D{f_of_t_name})); - } + const auto domain_creator = TestHelpers::test_option_tag< + domain::OptionTags::DomainCreator<3>, + TestHelpers::domain::BoundaryConditions:: + MetavariablesWithoutBoundaryConditions<3, domain::creators::Sphere>>( + "Sphere:\n" + " InnerRadius: 1\n" + " OuterRadius: 3\n" + " InnerCubeSphericity: 0.0\n" + " InitialRefinement: 2\n" + " InitialGridPoints: [3, 3, 4]\n" + " UseEquiangularMap: false\n" + " RadialPartitioning: []\n" + " RadialDistribution: [Linear]\n" + " TimeDependence:\n" + " UniformTranslation:\n" + " InitialTime: 1.0\n" + " Velocity: [2.3, -0.3, 0.5]\n"); + const auto* sphere = + dynamic_cast(domain_creator.get()); - test_sphere_construction( - dynamic_cast(*sphere), inner_radius, - outer_radius, sphericity, use_equiangular_map, grid_points_r_angular, - {7, make_array<3>(refinement_level)}, false, - std::make_tuple( - std::pair>{ - f_of_t_name, - {initial_time, - std::array{{{3, 0.0}, velocity, {3, 0.0}}}, - expiration_times.at(f_of_t_name)}}), - std::move(grid_to_inertial_maps), std::move(time_dependencies), - expiration_times); - }; - - std::unordered_map initial_expiration_times{ - {f_of_t_name, std::numeric_limits::infinity()}}; - helper(initial_expiration_times); - initial_expiration_times.at(f_of_t_name) = 10.0; - helper(initial_expiration_times); + const double inner_radius = 1.0; + const double outer_radius = 3.0; + const std::vector radial_partitioning = {}; + const std::array velocity{{2.3, -0.3, 0.5}}; + + for (const auto& expiration_time : make_array(1.0, 10.0)) { + test_sphere_construction(*sphere, inner_radius, outer_radius, + radial_partitioning, false, expiration_time, + velocity); + } } } // namespace +// [[TimeOut, 15]] SPECTRE_TEST_CASE("Unit.Domain.Creators.Sphere", "[Domain][Unit]") { domain::creators::time_dependence::register_derived_with_charm(); - test_sphere_boundaries_equiangular(); - test_sphere_factory_equiangular(); - test_sphere_boundaries_equidistant(); - test_sphere_factory_equidistant(); + test_parse_errors(); + test_sphere_boundaries(); test_sphere_factory_time_dependent(); } } // namespace domain diff --git a/tests/Unit/Domain/Test_DomainHelpers.cpp b/tests/Unit/Domain/Test_DomainHelpers.cpp index 4398f6dd26a8..d036750f7bfc 100644 --- a/tests/Unit/Domain/Test_DomainHelpers.cpp +++ b/tests/Unit/Domain/Test_DomainHelpers.cpp @@ -222,29 +222,6 @@ void test_wedge_errors() { }()), Catch::Contains("If we are using half wedges we must also be using " "ShellWedges::All.")); - - CHECK_THROWS_WITH( - ([]() { - const double inner_radius = 0.5; - const double outer_radius = 2.0; - const double inner_sphericity = 0.9; - const double outer_sphericity = 1.0; - const bool use_equiangular_map = true; - const bool use_half_wedges = true; - std::vector radial_partitioning{1., 1.5}; - const std::vector - radial_distribution{ - domain::CoordinateMaps::Distribution::Logarithmic, - domain::CoordinateMaps::Distribution::Logarithmic, - domain::CoordinateMaps::Distribution::Logarithmic}; - const ShellWedges which_wedges = ShellWedges::All; - static_cast(sph_wedge_coordinate_maps( - inner_radius, outer_radius, inner_sphericity, outer_sphericity, - use_equiangular_map, use_half_wedges, radial_partitioning, - radial_distribution, which_wedges)); - }()), - Catch::Contains("If we are using more than one layer the inner and outer " - "sphericities must match.")); #endif } diff --git a/tests/Unit/Helpers/Domain/Creators/TestHelpers.hpp b/tests/Unit/Helpers/Domain/Creators/TestHelpers.hpp index 61d1015046f7..f6bed7175d71 100644 --- a/tests/Unit/Helpers/Domain/Creators/TestHelpers.hpp +++ b/tests/Unit/Helpers/Domain/Creators/TestHelpers.hpp @@ -65,6 +65,20 @@ Domain test_domain_creator(const DomainCreator& domain_creator, } } } + { + INFO( + "Test block neighbors are never in the same direction as external " + "boundaries") + for (size_t block_id = 0; block_id < block_names.size(); ++block_id) { + for (const auto& neighbor : blocks[block_id].neighbors()) { + // external and neighbor directions should never match + const auto& external_boundaries = + blocks[block_id].external_boundaries(); + CHECK(external_boundaries.find(neighbor.first) == + external_boundaries.end()); + } + } + } } ::domain::creators::register_derived_with_charm();