From 6b915a04420c2a9f022c19095322185334989f16 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Fri, 13 Dec 2024 22:26:23 -0800 Subject: [PATCH] Create spherical shell in Sphere domain --- src/Domain/CMakeLists.txt | 1 + src/Domain/CoordinateMaps/CMakeLists.txt | 2 + src/Domain/CoordinateMaps/ShellType.cpp | 39 +++++ src/Domain/CoordinateMaps/ShellType.hpp | 44 +++++ src/Domain/Creators/BinaryCompactObject.cpp | 62 ++++--- src/Domain/Creators/Sphere.cpp | 133 ++++++++++++--- src/Domain/Creators/Sphere.hpp | 33 +++- src/Domain/Domain.cpp | 10 +- src/Domain/Domain.hpp | 7 +- src/Domain/DomainHelpers.cpp | 139 ---------------- src/Domain/DomainHelpers.hpp | 11 +- src/Domain/DomainHelpers.tpp | 151 ++++++++++++++++++ .../ExportCoordinates/ExportCoordinates.hpp | 62 +++---- .../InputFiles/ExportCoordinates/Input3D.yaml | 24 +-- 14 files changed, 478 insertions(+), 240 deletions(-) create mode 100644 src/Domain/CoordinateMaps/ShellType.cpp create mode 100644 src/Domain/CoordinateMaps/ShellType.hpp create mode 100644 src/Domain/DomainHelpers.tpp diff --git a/src/Domain/CMakeLists.txt b/src/Domain/CMakeLists.txt index 39833400bf6a..693bd08fc694 100644 --- a/src/Domain/CMakeLists.txt +++ b/src/Domain/CMakeLists.txt @@ -40,6 +40,7 @@ spectre_target_headers( CreateInitialElement.hpp Domain.hpp DomainHelpers.hpp + DomainHelpers.tpp ElementDistribution.hpp ElementLogicalCoordinates.hpp ElementMap.hpp diff --git a/src/Domain/CoordinateMaps/CMakeLists.txt b/src/Domain/CoordinateMaps/CMakeLists.txt index 340cc7d147c2..085dabc7167e 100644 --- a/src/Domain/CoordinateMaps/CMakeLists.txt +++ b/src/Domain/CoordinateMaps/CMakeLists.txt @@ -31,6 +31,7 @@ spectre_target_sources( Interval.cpp KerrHorizonConforming.cpp Rotation.cpp + ShellType.cpp SpecialMobius.cpp SphericalToCartesianPfaffian.cpp UniformCylindricalEndcap.cpp @@ -72,6 +73,7 @@ spectre_target_headers( ProductMaps.hpp ProductMaps.tpp Rotation.hpp + ShellType.hpp SpecialMobius.hpp SphericalToCartesianPfaffian.hpp Tags.hpp diff --git a/src/Domain/CoordinateMaps/ShellType.cpp b/src/Domain/CoordinateMaps/ShellType.cpp new file mode 100644 index 000000000000..7da241165741 --- /dev/null +++ b/src/Domain/CoordinateMaps/ShellType.cpp @@ -0,0 +1,39 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Domain/CoordinateMaps/ShellType.hpp" + +#include +#include + +#include "Options/Options.hpp" +#include "Options/ParseOptions.hpp" +#include "Utilities/ErrorHandling/Error.hpp" + +namespace domain::CoordinateMaps { + +std::ostream& operator<<(std::ostream& os, const ShellType shell_type) { + switch (shell_type) { + case ShellType::Cubed: + return os << "Cubed"; + case ShellType::Spherical: + return os << "Spherical"; + default: + ERROR("Unknown domain::CoordinateMaps::ShellType"); + } +} + +} // namespace domain::CoordinateMaps + +template <> +domain::CoordinateMaps::ShellType +Options::create_from_yaml::create( + const Options::Option& options) { + const auto shell_type = options.parse_as(); + if (shell_type == "Cubed") { + return domain::CoordinateMaps::ShellType::Cubed; + } else if (shell_type == "Spherical") { + return domain::CoordinateMaps::ShellType::Spherical; + } + PARSE_ERROR(options.context(), "ShellType must be 'Cubed' or 'Spherical'."); +} diff --git a/src/Domain/CoordinateMaps/ShellType.hpp b/src/Domain/CoordinateMaps/ShellType.hpp new file mode 100644 index 000000000000..a8f42112b75e --- /dev/null +++ b/src/Domain/CoordinateMaps/ShellType.hpp @@ -0,0 +1,44 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +/// \cond +namespace Options { +class Option; +template +struct create_from_yaml; +} // namespace Options +/// \endcond + +namespace domain::CoordinateMaps { + +/*! + * \brief Type of spherical shell: built from four (2D) or six (3D) wedges + * ("cubed") or from a single spherical shell ("spherical"). + * + * Used to select the shell type in the input file. + */ +enum class ShellType { + Cubed, + Spherical, +}; + +std::ostream& operator<<(std::ostream& os, ShellType shell_type); + +} // namespace domain::CoordinateMaps + +template <> +struct Options::create_from_yaml { + template + static domain::CoordinateMaps::ShellType create( + const Options::Option& options) { + return create(options); + } +}; +template <> +domain::CoordinateMaps::ShellType +Options::create_from_yaml::create( + const Options::Option& options); diff --git a/src/Domain/Creators/BinaryCompactObject.cpp b/src/Domain/Creators/BinaryCompactObject.cpp index cda5a42579d2..edb713cc2bd5 100644 --- a/src/Domain/Creators/BinaryCompactObject.cpp +++ b/src/Domain/Creators/BinaryCompactObject.cpp @@ -36,6 +36,7 @@ #include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp" #include "Domain/Domain.hpp" #include "Domain/DomainHelpers.hpp" +#include "Domain/DomainHelpers.tpp" #include "Domain/ExcisionSphere.hpp" #include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp" #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" @@ -459,20 +460,19 @@ Domain<3> BinaryCompactObject::create_domain() const { length_inner_cube_ * 0.5, std::array{{offset_x_coord_a_, 0.0, 0.0}})); Maps maps_center_A = - domain::make_vector_coordinate_map_base( - sph_wedge_coordinate_maps( - object_a.inner_radius, object_a.outer_radius, - inner_sphericity_A, 1.0, use_equiangular_map_, - offset_a_optional, false, {}, object_A_radial_distribution), - translation_A); + spherical_shells_coordinate_maps( + object_a.inner_radius, object_a.outer_radius, inner_sphericity_A, + 1.0, use_equiangular_map_, offset_a_optional, false, {}, + object_A_radial_distribution, + {domain::CoordinateMaps::ShellType::Cubed}, ShellWedges::All, + M_PI_2, translation_A); Maps maps_cube_A = - domain::make_vector_coordinate_map_base( - sph_wedge_coordinate_maps( - object_a.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_, - 1.0, 0.0, use_equiangular_map_, offset_a_optional), - translation_A); + spherical_shells_coordinate_maps( + object_a.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_, 1.0, + 0.0, use_equiangular_map_, offset_a_optional, false, {}, + {domain::CoordinateMaps::Distribution::Linear}, + {domain::CoordinateMaps::ShellType::Cubed}, ShellWedges::All, + M_PI_2, translation_A); std::move(maps_center_A.begin(), maps_center_A.end(), std::back_inserter(maps)); std::move(maps_cube_A.begin(), maps_cube_A.end(), std::back_inserter(maps)); @@ -502,20 +502,19 @@ Domain<3> BinaryCompactObject::create_domain() const { length_inner_cube_ * 0.5, std::array{{offset_x_coord_b_, 0.0, 0.0}})); Maps maps_center_B = - domain::make_vector_coordinate_map_base( - sph_wedge_coordinate_maps( - object_b.inner_radius, object_b.outer_radius, - inner_sphericity_B, 1.0, use_equiangular_map_, - offset_b_optional, false, {}, object_B_radial_distribution), - translation_B); + spherical_shells_coordinate_maps( + object_b.inner_radius, object_b.outer_radius, inner_sphericity_B, + 1.0, use_equiangular_map_, offset_b_optional, false, {}, + object_B_radial_distribution, + {domain::CoordinateMaps::ShellType::Cubed}, ShellWedges::All, + M_PI_2, translation_B); Maps maps_cube_B = - domain::make_vector_coordinate_map_base( - sph_wedge_coordinate_maps( - object_b.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_, - 1.0, 0.0, use_equiangular_map_, offset_b_optional), - translation_B); + spherical_shells_coordinate_maps( + object_b.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_, 1.0, + 0.0, use_equiangular_map_, offset_b_optional, false, {}, + {domain::CoordinateMaps::Distribution::Linear}, + {domain::CoordinateMaps::ShellType::Cubed}, ShellWedges::All, + M_PI_2, translation_B); std::move(maps_center_B.begin(), maps_center_B.end(), std::back_inserter(maps)); std::move(maps_cube_B.begin(), maps_cube_B.end(), std::back_inserter(maps)); @@ -545,12 +544,11 @@ Domain<3> BinaryCompactObject::create_domain() const { // --- Outer spherical shell (10 blocks) --- Maps maps_outer_shell = - domain::make_vector_coordinate_map_base( - sph_wedge_coordinate_maps(envelope_radius_, outer_radius_, 1.0, 1.0, - use_equiangular_map_, std::nullopt, true, - {}, {radial_distribution_outer_shell_}, - ShellWedges::All, opening_angle_)); + spherical_shells_coordinate_maps( + envelope_radius_, outer_radius_, 1.0, 1.0, use_equiangular_map_, + std::nullopt, true, {}, {radial_distribution_outer_shell_}, + {domain::CoordinateMaps::ShellType::Cubed}, ShellWedges::All, + opening_angle_); std::move(maps_outer_shell.begin(), maps_outer_shell.end(), std::back_inserter(maps)); diff --git a/src/Domain/Creators/Sphere.cpp b/src/Domain/Creators/Sphere.cpp index 668afdc5fb22..4d752e2ca330 100644 --- a/src/Domain/Creators/Sphere.cpp +++ b/src/Domain/Creators/Sphere.cpp @@ -30,6 +30,7 @@ #include "Domain/Creators/TimeDependence/None.hpp" #include "Domain/Domain.hpp" #include "Domain/DomainHelpers.hpp" +#include "Domain/DomainHelpers.tpp" #include "Domain/Structure/BlockNeighbor.hpp" #include "Domain/Structure/Direction.hpp" #include "Domain/Structure/DirectionMap.hpp" @@ -37,6 +38,7 @@ #include "Utilities/ErrorHandling/Assert.hpp" #include "Utilities/GetOutput.hpp" #include "Utilities/MakeArray.hpp" +#include "Utilities/Overloader.hpp" namespace Frame { struct Inertial; @@ -79,7 +81,7 @@ Sphere::Sphere( std::optional equatorial_compression, std::vector radial_partitioning, const typename RadialDistribution::type& radial_distribution, - ShellWedges which_wedges, + const typename ShellType::type& shell_type, ShellWedges which_wedges, std::optional time_dependent_options, std::unique_ptr outer_boundary_condition, @@ -153,22 +155,72 @@ Sphere::Sphere( "Add entries to 'RadialPartitioning' to add outer shells for " "which you can select different radial distributions."); } + shell_types_ = std::visit( + Overloader{ + [&num_shells = num_shells_]( + const domain::CoordinateMaps::ShellType& uniform_shell_type) { + return std::vector( + num_shells, uniform_shell_type); + }, + [](const std::vector& + shell_types) { return shell_types; }}, + shell_type); + if (shell_types_.size() != num_shells_) { + PARSE_ERROR(context, + "Specify a 'ShellType' for every spherical shell. You " + "specified " + << shell_types_.size() << " items, but the domain has " + << num_shells_ << " shells."); + } + if (fill_interior_ and + shell_types_.front() != domain::CoordinateMaps::ShellType::Cubed) { + PARSE_ERROR(context, + "The 'ShellType' must be 'Cubed' for the innermost " + "shell filled with a cube. " + "Add entries to 'RadialPartitioning' to add outer shells for " + "which you can select different shell types."); + } + if (which_wedges_ != ShellWedges::All and + shell_types_ != + std::vector{ + num_shells_, domain::CoordinateMaps::ShellType::Cubed}) { + PARSE_ERROR(context, + "To specify 'ShellWedges' other than 'All', you must specify " + "'ShellType' as 'Cubed' for all shells."); + } // Create grid anchors grid_anchors_["Center"] = tnsr::I{std::array{0.0, 0.0, 0.0}}; // Determine number of blocks - num_blocks_per_shell_ = which_wedges_ == ShellWedges::All ? 6 - : which_wedges_ == ShellWedges::FourOnEquator ? 4 + num_blocks_per_cubed_shell_ = which_wedges_ == ShellWedges::All ? 6 + : which_wedges_ == ShellWedges::FourOnEquator + ? 4 : 1; - num_blocks_ = num_blocks_per_shell_ * num_shells_ + (fill_interior_ ? 1 : 0); + num_blocks_ = alg::accumulate( + shell_types_, 0_st, + [num_blocks_per_cubed_shell = num_blocks_per_cubed_shell_]( + const size_t num_blocks, + const domain::CoordinateMaps::ShellType shell_type) { + if (shell_type == domain::CoordinateMaps::ShellType::Cubed) { + return num_blocks + num_blocks_per_cubed_shell; + } else if (shell_type == domain::CoordinateMaps::ShellType::Spherical) { + return num_blocks + 1; + } else { + ERROR("Unknown ShellType"); + } + }); + if (fill_interior_) { + ++num_blocks_; + } // 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); + if (shell_types_[shell] == domain::CoordinateMaps::ShellType::Cubed) { for (size_t direction = which_wedge_index(which_wedges_); direction < 6; ++direction) { const std::string wedge_name = @@ -176,6 +228,11 @@ Sphere::Sphere( block_names_.emplace_back(wedge_name); block_groups_[shell_prefix].insert(wedge_name); block_groups_["Wedges"].insert(wedge_name); + } + } else if (shell_types_[shell] == + domain::CoordinateMaps::ShellType::Spherical) { + block_names_.emplace_back(shell_prefix); + block_groups_["SphericalShells"].insert(shell_prefix); } } if (fill_interior_) { @@ -264,10 +321,9 @@ Sphere::Sphere( } Domain<3> Sphere::create_domain() const { - std::vector> corners = - corners_for_radially_layered_domains(num_shells_, fill_interior_, - {{1, 2, 3, 4, 5, 6, 7, 8}}, - which_wedges_); + std::vector> blocks{}; + blocks.reserve(num_blocks_); + double aspect_ratio = 1.0; size_t index_polar_axis = 2; if (equatorial_compression_.has_value()) { @@ -277,14 +333,48 @@ Domain<3> Sphere::create_domain() const { const domain::CoordinateMaps::EquatorialCompression compression{ aspect_ratio, index_polar_axis}; - auto coord_maps = domain::make_vector_coordinate_map_base( - sph_wedge_coordinate_maps( + auto coord_maps = + spherical_shells_coordinate_maps( inner_radius_, outer_radius_, fill_interior_ ? std::get(interior_).sphericity : 1.0, 1.0, use_equiangular_map_, std::nullopt, false, radial_partitioning_, - radial_distribution_, which_wedges_), - compression); + radial_distribution_, shell_types_, which_wedges_, M_PI, compression); + + // Assuming only the outer shell is a spherical shell + const auto corners_of_wedges = corners_for_radially_layered_domains( + num_shells_ - 1, fill_interior_, {{1, 2, 3, 4, 5, 6, 7, 8}}, + which_wedges_); + std::vector>> neighbors_of_wedges{}; + set_internal_boundaries<3>(make_not_null(&neighbors_of_wedges), + corners_of_wedges); + std::unordered_set> neighbors_of_outer_shell{}; + for (size_t i = 0; i < coord_maps.size() - 1; i++) { + if (fill_interior_ and i < 7) { + // Adjacent to the inner cube + neighbors_of_wedges[i][Direction<3>::lower_zeta()] = BlockNeighbor<3>{ + num_blocks_ - 1, + neighbors_of_wedges[i][Direction<3>::lower_zeta()].orientation()}; + } + if (i >= coord_maps.size() - 7) { + // Adjacent to the outer shell + // Account for different ordering of logical coordinates in wedges [theta, + // phi, r] and spherical shells [r, theta, phi] + const OrientationMap<3> orientation{ + {{Direction<3>::upper_zeta(), Direction<3>::upper_xi(), + Direction<3>::upper_eta()}}, + {{Direction<3>::upper_xi(), Direction<3>::upper_eta(), + Direction<3>::upper_zeta()}}}; + neighbors_of_wedges[i][Direction<3>::upper_zeta()] = + BlockNeighbor<3>{coord_maps.size() - 1, orientation}; + neighbors_of_outer_shell.emplace(i, orientation.inverse_map()); + } + blocks.emplace_back(std::move(coord_maps[i]), i, + std::move(neighbors_of_wedges[i]), block_names_[i]); + } + blocks.emplace_back(Block<3>( + std::move(coord_maps.back()), coord_maps.size() - 1, + {{Direction<3>::lower_xi(), std::move(neighbors_of_outer_shell)}}, + block_names_[coord_maps.size() - 1])); std::unordered_map> excision_spheres{}; @@ -318,6 +408,9 @@ Domain<3> Sphere::create_domain() const { BulgedCube{inner_radius_, inner_cube_sphericity, use_equiangular_map_})); } + blocks.emplace_back(std::move(coord_maps.back()), num_blocks_ - 1, + std::move(neighbors_of_wedges.back()), + block_names_.back()); } else { // Set up excision sphere only for ShellWedges::All // - The first 6 blocks enclose the excised sphere, see @@ -338,14 +431,14 @@ Domain<3> Sphere::create_domain() const { } } - Domain<3> domain{std::move(coord_maps), corners, {}, - std::move(excision_spheres), block_names_, block_groups_}; + Domain<3> domain{std::move(blocks), std::move(excision_spheres), + block_groups_}; ASSERT(domain.blocks().size() == num_blocks_, "Unexpected number of blocks. Expected " << num_blocks_ << " but created " << domain.blocks().size() << "."); - if (time_dependent_options_.has_value()) { + if (time_dependent_options_.has_value()) { // TODO std::vector>> block_maps_grid_to_inertial{num_blocks_}; @@ -362,18 +455,18 @@ Domain<3> Sphere::create_domain() const { time_dependent_options_.value()); const size_t first_block_outer_shell = - (num_shells_ - 1) * num_blocks_per_shell_; + (num_shells_ - 1) * num_blocks_per_cubed_shell_; for (size_t block_id = 0; block_id < num_blocks_; block_id++) { const bool is_outer_shell = block_id >= first_block_outer_shell and - block_id < (first_block_outer_shell + num_blocks_per_shell_); + block_id < (first_block_outer_shell + num_blocks_per_cubed_shell_); const bool is_inner_cube = fill_interior_ and (block_id == num_blocks_ - 1); // Correct for 'which_wedges' option - const size_t shell = block_id / num_blocks_per_shell_; + const size_t shell = block_id / num_blocks_per_cubed_shell_; const size_t block_number = shell * 6 + which_wedge_index(which_wedges_) + - block_id % num_blocks_per_shell_; + block_id % num_blocks_per_cubed_shell_; block_maps_grid_to_distorted[block_id] = hard_coded_options.grid_to_distorted_map(block_number, is_inner_cube); diff --git a/src/Domain/Creators/Sphere.hpp b/src/Domain/Creators/Sphere.hpp index f70a2e2dea5b..0b3d4518e3a9 100644 --- a/src/Domain/Creators/Sphere.hpp +++ b/src/Domain/Creators/Sphere.hpp @@ -15,6 +15,9 @@ #include "Domain/BoundaryConditions/BoundaryCondition.hpp" #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp" #include "Domain/CoordinateMaps/Distribution.hpp" +#include "Domain/CoordinateMaps/Identity.hpp" +#include "Domain/CoordinateMaps/ShellType.hpp" +#include "Domain/CoordinateMaps/SphericalToCartesianPfaffian.hpp" #include "Domain/Creators/DomainCreator.hpp" #include "Domain/Creators/TimeDependence/TimeDependence.hpp" #include "Domain/Creators/TimeDependentOptions/Sphere.hpp" @@ -178,10 +181,17 @@ class Sphere : public DomainCreator<3> { using Equiangular3D = CoordinateMaps::ProductOf3Maps; using BulgedCube = CoordinateMaps::BulgedCube; + using LogicalToSphericalMap = domain::CoordinateMaps::ProductOf3Maps< + domain::CoordinateMaps::Interval, domain::CoordinateMaps::Identity<1>, + domain::CoordinateMaps::Identity<1>>; public: using maps_list = tmpl::append< tmpl::list< + domain::CoordinateMap< + Frame::BlockLogical, Frame::Inertial, LogicalToSphericalMap, + domain::CoordinateMaps::SphericalToCartesianPfaffian, + domain::CoordinateMaps::EquatorialCompression>, // Inner cube domain::CoordinateMap, @@ -311,6 +321,22 @@ class Sphere : public DomainCreator<3> { "partitions."}; }; + struct ShellType { + using type = std::variant>; + static constexpr Options::String help = { + "Select type of each spherical shell. Specify a list of N+1 shell " + "types for N radial partitions, i.e., one for each shell." + "The shell type can be 'Cubed' to create six wedges that compose the " + "shell, or 'Spherical' to create a single spherical shell with " + "spherical harmonics. The 'Spherical' option " + "is not widely supported yet, so it is unlikely to work. " + "If the interior of the sphere is filled with a cube, the innermost " + "shell must have a 'Cubed' shell type. You can also specify just a " + "single shell type (not in a vector) which will use the same shell " + "type for all partitions."}; + }; + struct WhichWedges { using type = ShellWedges; static constexpr Options::String help = { @@ -340,7 +366,7 @@ class Sphere : public DomainCreator<3> { using basic_options = tmpl::list; template @@ -371,6 +397,8 @@ class Sphere : public DomainCreator<3> { std::vector radial_partitioning = {}, const typename RadialDistribution::type& radial_distribution = domain::CoordinateMaps::Distribution::Linear, + const typename ShellType::type& shell_type = + domain::CoordinateMaps::ShellType::Cubed, ShellWedges which_wedges = ShellWedges::All, std::optional time_dependent_options = std::nullopt, std::unique_ptr @@ -428,6 +456,7 @@ class Sphere : public DomainCreator<3> { std::optional equatorial_compression_{}; std::vector radial_partitioning_{}; std::vector radial_distribution_{}; + std::vector shell_types_{}; ShellWedges which_wedges_ = ShellWedges::All; std::optional time_dependent_options_{}; bool use_hard_coded_maps_{false}; @@ -435,7 +464,7 @@ class Sphere : public DomainCreator<3> { outer_boundary_condition_; size_t num_shells_; size_t num_blocks_; - size_t num_blocks_per_shell_; + size_t num_blocks_per_cubed_shell_; std::vector block_names_{}; std::unordered_map> block_groups_{}; diff --git a/src/Domain/Domain.cpp b/src/Domain/Domain.cpp index 1c956f278835..925812b040ee 100644 --- a/src/Domain/Domain.cpp +++ b/src/Domain/Domain.cpp @@ -21,8 +21,14 @@ struct BlockLogical; } // namespace Frame template -Domain::Domain(std::vector> blocks) - : blocks_(std::move(blocks)) {} +Domain::Domain( + std::vector> blocks, + std::unordered_map> excision_spheres, + std::unordered_map> + block_groups) + : blocks_(std::move(blocks)), + excision_spheres_(std::move(excision_spheres)), + block_groups_(std::move(block_groups)) {} template Domain::Domain( diff --git a/src/Domain/Domain.hpp b/src/Domain/Domain.hpp index b2822d62187a..50e48bfc44e5 100644 --- a/src/Domain/Domain.hpp +++ b/src/Domain/Domain.hpp @@ -98,7 +98,12 @@ namespace domain {} template class Domain { public: - explicit Domain(std::vector> blocks); + explicit Domain( + std::vector> blocks, + std::unordered_map> + excision_spheres = {}, + std::unordered_map> + block_groups = {}); /*! * \brief Create a Domain using CoordinateMaps to encode the Orientations. diff --git a/src/Domain/DomainHelpers.cpp b/src/Domain/DomainHelpers.cpp index 8051cf1b281d..e7990043e2c3 100644 --- a/src/Domain/DomainHelpers.cpp +++ b/src/Domain/DomainHelpers.cpp @@ -592,145 +592,6 @@ size_t which_wedge_index(const ShellWedges& which_wedges) { } } -std::vector> sph_wedge_coordinate_maps( - const double inner_radius, const double outer_radius, - const double inner_sphericity, const double outer_sphericity, - const bool use_equiangular_map, - const std::optional>>& - offset_options, - const bool use_half_wedges, const std::vector& radial_partitioning, - const std::vector& - radial_distribution, - const ShellWedges which_wedges, const double opening_angle) { - 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_distribution.size() == 1 + radial_partitioning.size(), - "Specify a radial distribution for every spherical shell. You " - "specified " - << radial_distribution.size() << " items, but the domain has " - << 1 + radial_partitioning.size() << " shells."); - - const auto wedge_orientations = orientations_for_sphere_wrappings(); - - using Wedge3DMap = domain::CoordinateMaps::Wedge<3>; - using Halves = Wedge3DMap::WedgeHalves; - std::vector wedges_for_all_layers{}; - - const size_t number_of_layers = 1 + radial_partitioning.size(); - double temp_inner_radius = inner_radius; - double temp_inner_sphericity = inner_sphericity; - if (offset_options.has_value()) { - for (size_t layer_i = 0; layer_i < number_of_layers; layer_i++) { - const auto& radial_distribution_this_layer = - radial_distribution.at(layer_i); - std::optional optional_outer_radius{}; - if (outer_sphericity != 0.0) { - optional_outer_radius = outer_radius; - } else { - optional_outer_radius = std::nullopt; - } - // Generate wedges/half-wedges a layer at a time. - std::vector wedges_for_this_layer{}; - if (not use_half_wedges) { - for (size_t face_j = which_wedge_index(which_wedges); face_j < 6; - face_j++) { - wedges_for_this_layer.emplace_back( - temp_inner_radius, optional_outer_radius, - offset_options.value().first, offset_options.value().second, - 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, optional_outer_radius, - offset_options.value().first, offset_options.value().second, - gsl::at(wedge_orientations, i), use_equiangular_map, - Halves::LowerOnly, radial_distribution_this_layer); - wedges_for_this_layer.emplace_back( - temp_inner_radius, optional_outer_radius, - offset_options.value().first, offset_options.value().second, - gsl::at(wedge_orientations, i), use_equiangular_map, - Halves::UpperOnly, radial_distribution_this_layer); - } - wedges_for_this_layer.emplace_back( - temp_inner_radius, optional_outer_radius, - offset_options.value().first, offset_options.value().second, - gsl::at(wedge_orientations, 4), use_equiangular_map, Halves::Both, - radial_distribution_this_layer); - wedges_for_this_layer.emplace_back( - temp_inner_radius, optional_outer_radius, - offset_options.value().first, offset_options.value().second, - gsl::at(wedge_orientations, 5), use_equiangular_map, Halves::Both, - radial_distribution_this_layer); - } - for (const auto& wedge : wedges_for_this_layer) { - wedges_for_all_layers.push_back(wedge); - } - } - } else { - double temp_outer_radius{}; - for (size_t layer_i = 0; layer_i < number_of_layers; layer_i++) { - const auto& radial_distribution_this_layer = - radial_distribution.at(layer_i); - if (layer_i != radial_partitioning.size()) { - temp_outer_radius = radial_partitioning.at(layer_i); - } else { - temp_outer_radius = outer_radius; - } - // Generate wedges/half-wedges a layer at a time. - std::vector wedges_for_this_layer{}; - if (not use_half_wedges) { - 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, temp_inner_sphericity, - outer_sphericity, gsl::at(wedge_orientations, face_j), - use_equiangular_map, Halves::Both, radial_distribution_this_layer, - std::array({{M_PI_2, M_PI_2}})); - } - } else { - for (size_t i = 0; i < 4; i++) { - wedges_for_this_layer.emplace_back( - 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, - std::array({{opening_angle, M_PI_2}})); - wedges_for_this_layer.emplace_back( - 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, - std::array({{opening_angle, M_PI_2}})); - } - const double endcap_opening_angle = M_PI - opening_angle; - const std::array endcap_opening_angles = { - {endcap_opening_angle, endcap_opening_angle}}; - wedges_for_this_layer.emplace_back( - 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, - endcap_opening_angles, false); - wedges_for_this_layer.emplace_back( - 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, - endcap_opening_angles, false); - } - for (const auto& wedge : wedges_for_this_layer) { - wedges_for_all_layers.push_back(wedge); - } - - if (layer_i != radial_partitioning.size()) { - temp_inner_radius = radial_partitioning.at(layer_i); - temp_inner_sphericity = outer_sphericity; - } - } - } - return wedges_for_all_layers; -} - std::vector frustum_coordinate_maps( const double length_inner_cube, const double length_outer_cube, const bool equiangular_map_at_outer, const bool equiangular_map_at_inner, diff --git a/src/Domain/DomainHelpers.hpp b/src/Domain/DomainHelpers.hpp index 6c9c802d28f9..40d80befd80d 100644 --- a/src/Domain/DomainHelpers.hpp +++ b/src/Domain/DomainHelpers.hpp @@ -16,6 +16,7 @@ #include "DataStructures/Index.hpp" #include "DataStructures/Tensor/Tensor.hpp" #include "Domain/CoordinateMaps/Distribution.hpp" +#include "Domain/CoordinateMaps/ShellType.hpp" #include "Domain/Structure/Direction.hpp" #include "Domain/Structure/Side.hpp" #include "Utilities/ConstantExpressions.hpp" @@ -179,7 +180,10 @@ size_t which_wedge_index(const ShellWedges& which_wedges); * of pi minus this opening angle. This parameter only has an effect if * `use_half_wedges` is set to `true`. */ -std::vector> sph_wedge_coordinate_maps( +template +std::vector< + std::unique_ptr>> +spherical_shells_coordinate_maps( double inner_radius, double outer_radius, double inner_sphericity, double outer_sphericity, bool use_equiangular_map, const std::optional>>& @@ -188,7 +192,10 @@ std::vector> sph_wedge_coordinate_maps( const std::vector& radial_partitioning = {}, const std::vector& radial_distribution = {domain::CoordinateMaps::Distribution::Linear}, - ShellWedges which_wedges = ShellWedges::All, double opening_angle = M_PI_2); + const std::vector& shell_types = + {domain::CoordinateMaps::ShellType::Cubed}, + ShellWedges which_wedges = ShellWedges::All, double opening_angle = M_PI_2, + const AppendMaps&... append_maps); /// \ingroup ComputationalDomainGroup /// These are the ten Frustums used in the DomainCreators for binary compact diff --git a/src/Domain/DomainHelpers.tpp b/src/Domain/DomainHelpers.tpp new file mode 100644 index 000000000000..f38ddce3d2ba --- /dev/null +++ b/src/Domain/DomainHelpers.tpp @@ -0,0 +1,151 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include "Domain/DomainHelpers.hpp" + +#include +#include +#include + +#include "Domain/CoordinateMaps/Affine.hpp" +#include "Domain/CoordinateMaps/CoordinateMap.tpp" +#include "Domain/CoordinateMaps/Distribution.hpp" +#include "Domain/CoordinateMaps/Identity.hpp" +#include "Domain/CoordinateMaps/Interval.hpp" +#include "Domain/CoordinateMaps/ProductMaps.tpp" +#include "Domain/CoordinateMaps/ShellType.hpp" +#include "Domain/CoordinateMaps/SphericalToCartesianPfaffian.hpp" +#include "Domain/CoordinateMaps/Wedge.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" +#include "Utilities/ErrorHandling/Error.hpp" + +template +std::vector< + std::unique_ptr>> +spherical_shells_coordinate_maps( + const double inner_radius, const double outer_radius, + const double inner_sphericity, const double outer_sphericity, + const bool use_equiangular_map, + const std::optional>>& + offset_options, + const bool use_half_wedges, const std::vector& radial_partitioning, + const std::vector& + radial_distribution, + const std::vector& shell_types, + const ShellWedges which_wedges, const double opening_angle, + const AppendMaps&... append_maps) { + 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_distribution.size() == 1 + radial_partitioning.size(), + "Specify a radial distribution for every spherical shell. You " + "specified " + << radial_distribution.size() << " items, but the domain has " + << 1 + radial_partitioning.size() << " shells."); + ASSERT(shell_types.size() == 1 + radial_partitioning.size(), + "Specify a type for every spherical shell. You specified " + << shell_types.size() << " items, but the domain has " + << 1 + radial_partitioning.size() << " shells."); + + const auto wedge_orientations = orientations_for_sphere_wrappings(); + + std::vector< + std::unique_ptr>> + maps{}; + + const size_t number_of_layers = 1 + radial_partitioning.size(); + double temp_inner_radius = inner_radius; + double temp_inner_sphericity = inner_sphericity; + if (offset_options.has_value()) { + ERROR("TODO"); + } else { + double temp_outer_radius{}; + for (size_t layer_i = 0; layer_i < number_of_layers; layer_i++) { + // Determine inner and outer radius for this shell + if (layer_i != radial_partitioning.size()) { + temp_outer_radius = radial_partitioning.at(layer_i); + } else { + temp_outer_radius = outer_radius; + } + if (shell_types[layer_i] == domain::CoordinateMaps::ShellType::Cubed) { + // Compose the shell of wedges (deformed cubes) + using Wedge3DMap = domain::CoordinateMaps::Wedge<3>; + using Halves = Wedge3DMap::WedgeHalves; + std::vector wedges{}; + if (not use_half_wedges) { + for (size_t face_j = which_wedge_index(which_wedges); face_j < 6; + face_j++) { + wedges.emplace_back( + temp_inner_radius, temp_outer_radius, temp_inner_sphericity, + outer_sphericity, gsl::at(wedge_orientations, face_j), + use_equiangular_map, Halves::Both, radial_distribution[layer_i], + std::array({{M_PI_2, M_PI_2}})); + } + } else { + for (size_t i = 0; i < 4; i++) { + wedges.emplace_back( + temp_inner_radius, temp_outer_radius, temp_inner_sphericity, + outer_sphericity, gsl::at(wedge_orientations, i), + use_equiangular_map, Halves::LowerOnly, + radial_distribution[layer_i], + std::array({{opening_angle, M_PI_2}})); + wedges.emplace_back( + temp_inner_radius, temp_outer_radius, temp_inner_sphericity, + outer_sphericity, gsl::at(wedge_orientations, i), + use_equiangular_map, Halves::UpperOnly, + radial_distribution[layer_i], + std::array({{opening_angle, M_PI_2}})); + } + const double endcap_opening_angle = M_PI - opening_angle; + const std::array endcap_opening_angles = { + {endcap_opening_angle, endcap_opening_angle}}; + wedges.emplace_back( + temp_inner_radius, temp_outer_radius, temp_inner_sphericity, + outer_sphericity, gsl::at(wedge_orientations, 4), + use_equiangular_map, Halves::Both, radial_distribution[layer_i], + endcap_opening_angles, false); + wedges.emplace_back( + temp_inner_radius, temp_outer_radius, temp_inner_sphericity, + outer_sphericity, gsl::at(wedge_orientations, 5), + use_equiangular_map, Halves::Both, radial_distribution[layer_i], + endcap_opening_angles, false); + } + auto wedge_maps = + domain::make_vector_coordinate_map_base(std::move(wedges), + append_maps...); + maps.insert(maps.end(), std::make_move_iterator(wedge_maps.begin()), + std::make_move_iterator(wedge_maps.end())); + } else if (shell_types[layer_i] == + domain::CoordinateMaps::ShellType::Spherical) { + // Use a single spherical shell (with spherical harmonics in angular + // directions) + ASSERT( + inner_sphericity == 1. and outer_sphericity == 1., + "Spherical shells only support inner and outer sphericity == 1."); + domain::CoordinateMaps::ProductOf3Maps< + domain::CoordinateMaps::Interval, + domain::CoordinateMaps::Identity<1>, + domain::CoordinateMaps::Identity<1>> + logical_to_spherical_map{{-1., 1., temp_inner_radius, outer_radius, + radial_distribution[layer_i], 0.}, + {}, + {}}; + maps.push_back( + domain::make_coordinate_map_base( + std::move(logical_to_spherical_map), + domain::CoordinateMaps::SphericalToCartesianPfaffian{}, + append_maps...)); + } else { + ERROR("Unknown ShellType"); + } + + if (layer_i != radial_partitioning.size()) { + temp_inner_radius = radial_partitioning.at(layer_i); + temp_inner_sphericity = outer_sphericity; + } + } // for layer + } // if offset_options + return maps; +} diff --git a/src/Executables/ExportCoordinates/ExportCoordinates.hpp b/src/Executables/ExportCoordinates/ExportCoordinates.hpp index d7853a1e5aae..d2e8f3e3aca3 100644 --- a/src/Executables/ExportCoordinates/ExportCoordinates.hpp +++ b/src/Executables/ExportCoordinates/ExportCoordinates.hpp @@ -130,13 +130,13 @@ struct ExportCoordinates { const double time = get(box); const auto& mesh = get>(box); - const auto& inv_jacobian = - db::get>(box); + // const auto& inv_jacobian = + // db::get>(box); const auto& inertial_coordinates = db::get>(box); - const auto deriv_inertial_coordinates = - partial_derivative(inertial_coordinates, mesh, inv_jacobian); + // const auto deriv_inertial_coordinates = + // partial_derivative(inertial_coordinates, mesh, inv_jacobian); // Collect volume data // Remove tensor types, only storing individual components std::vector components; @@ -148,13 +148,13 @@ struct ExportCoordinates { inertial_coordinates.get(d)); } - for (size_t i = 0; i < deriv_inertial_coordinates.size(); ++i) { - components.emplace_back( - "DerivInertialCoordinates_" + - deriv_inertial_coordinates.component_name( - deriv_inertial_coordinates.get_tensor_index(i)), - deriv_inertial_coordinates[i]); - } + // for (size_t i = 0; i < deriv_inertial_coordinates.size(); ++i) { + // components.emplace_back( + // "DerivInertialCoordinates_" + + // deriv_inertial_coordinates.component_name( + // deriv_inertial_coordinates.get_tensor_index(i)), + // deriv_inertial_coordinates[i]); + // } // Also output the determinant of the inverse jacobian, which measures // the expansion and compression of the grid @@ -169,27 +169,27 @@ struct ExportCoordinates { // Also output the jacobian diagnostic, which compares the analytic // Jacobian (via the CoordinateMap) to the numerical Jacobian // (computed via logical_partial_derivative) - const auto& jacobian = determinant_and_inverse(inv_jacobian).second; - tnsr::i jac_diag{ - mesh.number_of_grid_points(), 0.0}; - domain::jacobian_diagnostic(make_not_null(&jac_diag), jacobian, - inertial_coordinates, mesh); - for (size_t i = 0; i < Dim; ++i) { - components.emplace_back( - "JacobianDiagnostic_" + - jac_diag.component_name(jac_diag.get_tensor_index(i)), - jac_diag.get(i)); - } + // const auto& jacobian = determinant_and_inverse(inv_jacobian).second; + // tnsr::i jac_diag{ + // mesh.number_of_grid_points(), 0.0}; + // domain::jacobian_diagnostic(make_not_null(&jac_diag), jacobian, + // inertial_coordinates, mesh); + // for (size_t i = 0; i < Dim; ++i) { + // components.emplace_back( + // "JacobianDiagnostic_" + + // jac_diag.component_name(jac_diag.get_tensor_index(i)), + // jac_diag.get(i)); + // } // Also output the computation domain metric - const auto& flat_logical_metric = - db::get>(box); - for (size_t i = 0; i < flat_logical_metric.size(); ++i) { - components.emplace_back( - db::tag_name>() + - flat_logical_metric.component_suffix(i), - flat_logical_metric[i]); - } + // const auto& flat_logical_metric = + // db::get>(box); + // for (size_t i = 0; i < flat_logical_metric.size(); ++i) { + // components.emplace_back( + // db::tag_name>() + + // flat_logical_metric.component_suffix(i), + // flat_logical_metric[i]); + // } // Send data to volume observer auto& local_observer = *Parallel::local_branch( diff --git a/tests/InputFiles/ExportCoordinates/Input3D.yaml b/tests/InputFiles/ExportCoordinates/Input3D.yaml index 919a4a658c49..0d83fa565aad 100644 --- a/tests/InputFiles/ExportCoordinates/Input3D.yaml +++ b/tests/InputFiles/ExportCoordinates/Input3D.yaml @@ -30,17 +30,19 @@ ResourceInfo: Singletons: Auto DomainCreator: - Brick: - LowerBound: [-0.5, -0.5, -0.5] - UpperBound: [0.5, 0.5, 0.5] - Distribution: [Linear, Linear, Linear] - IsPeriodicIn: [false, false, false] - InitialRefinement: [0, 0, 0] - InitialGridPoints: [3, 3, 3] - TimeDependence: - UniformTranslation: - InitialTime: 0.0 - Velocity: [0.5, 0.0, 0.0] + Sphere: + InnerRadius: 2. + OuterRadius: 10. + Interior: Excise + InitialRefinement: 0 + InitialGridPoints: 12 + UseEquiangularMap: True + EquatorialCompression: None + WhichWedges: All + RadialPartitioning: [4.] + RadialDistribution: [Linear, Logarithmic] + ShellType: [Cubed, Spherical] + TimeDependentMaps: None SpatialDiscretization: ActiveGrid: Dg