Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Support domain deformations in elliptic solver #5575

Merged
merged 6 commits into from
Oct 31, 2023
Merged

Support domain deformations in elliptic solver #5575

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
5b51973
Increase bracket depth for fold expressions in clang 12
nilsvu Oct 18, 2023
ad18822
Add ElementMap constructor from Block
nilsvu Oct 18, 2023
7852114
Forward functions of time in ElementMap
nilsvu Oct 18, 2023
9d7ec58
Only check FoT are ready if in mutable cache
nilsvu Oct 18, 2023
8b912dc
Support domain deformations in elliptic solver
nilsvu Oct 18, 2023
a5e0ae6
Add an instantiation
nilsvu Oct 23, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 10 additions & 0 deletions cmake/SetupCxxFlags.cmake
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,16 @@ set_property(TARGET SpectreFlags
INTERFACE_COMPILE_OPTIONS
$<$<COMPILE_LANGUAGE:CXX>:-ftemplate-backtrace-limit=0>)

# Increase bracket depth for fold expressions
# (see also https://github.com/llvm/llvm-project/issues/48973)
if (CMAKE_CXX_COMPILER_ID MATCHES "Clang"
AND CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 12.0)
set_property(TARGET SpectreFlags
APPEND PROPERTY
INTERFACE_COMPILE_OPTIONS
$<$<COMPILE_LANGUAGE:CXX>:-fbracket-depth=1024>)
endif()

# Disable cmath setting the error flag. This allows the compiler to more
# aggressively vectorize code since it doesn't need to respect some global
# state.
Expand Down
3 changes: 1 addition & 2 deletions src/Domain/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,15 +70,14 @@ target_link_libraries(
DomainBoundaryConditions
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the elliptic code do you actually evaluate the maps outside of compute items? If not, this commit, if not the previous as well could be replaced by adding a grid->inertial map into the elliptic box and using the time-dependent compute tags we use in the evolution code. Those tags don't do any extra work for "static" domains. I guess they currently require a Tags::Time to be in the box, but we could reasonably template on a bool or something to elide the time if it's not needed.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Currently the map is only needed during initialization; in multiple actions though (domain init, subdomain init, faces & mortars). I thought about removing the ElementMap tag from the DataBox, but kept it for now so I don't have to recreate the ElementMap in these places. Could be done though I think. When we update the map parameters we have to re-evaluate the map as well though, which might be easier if the ElementMap is available. I haven't thought that part through yet.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Alright, that sounds good. Thanks for the explanation!

DomainStructure
ErrorHandling
FunctionsOfTime
Options
Spectral
Utilities
PRIVATE
LinearOperators
RootFinding
SphericalHarmonics
INTERFACE
FunctionsOfTime
)

add_subdirectory(Amr)
Expand Down
51 changes: 51 additions & 0 deletions src/Domain/CoordinateMaps/Composition.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -365,6 +365,15 @@ bool Composition<Frames, Dim, std::index_sequence<Is...>>::is_equal_to(
#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)

#define INSTANTIATE(_, data) \
template class Composition< \
tmpl::list<Frame::BlockLogical, Frame::Grid, Frame::Inertial>, \
DIM(data)>; \
template class Composition< \
tmpl::list<Frame::BlockLogical, Frame::Grid, Frame::Distorted>, \
DIM(data)>; \
template class Composition<tmpl::list<Frame::BlockLogical, Frame::Grid, \
Frame::Distorted, Frame::Inertial>, \
DIM(data)>; \
template class Composition< \
tmpl::list<Frame::ElementLogical, Frame::BlockLogical, Frame::Inertial>, \
DIM(data)>; \
Expand All @@ -390,6 +399,48 @@ GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3))

#if defined(__clang__) && __clang_major__ >= 15
#define INSTANTIATE2(_, data) \
template domain::CoordinateMaps::Composition< \
brigand::list<Frame::BlockLogical, Frame::Grid, Frame::Distorted>, \
DIM(data), std::integer_sequence<unsigned long, 0ul, 1ul>>:: \
Composition( \
std::unique_ptr<domain::CoordinateMapBase<Frame::BlockLogical, \
Frame::Grid, DIM(data)>, \
std::default_delete<domain::CoordinateMapBase< \
Frame::BlockLogical, Frame::Grid, DIM(data)>>>, \
std::unique_ptr<domain::CoordinateMapBase< \
Frame::Grid, Frame::Distorted, DIM(data)>, \
std::default_delete<domain::CoordinateMapBase< \
Frame::Grid, Frame::Distorted, DIM(data)>>>); \
template domain::CoordinateMaps::Composition< \
brigand::list<Frame::BlockLogical, Frame::Grid, Frame::Inertial>, \
DIM(data), std::integer_sequence<unsigned long, 0ul, 1ul>>:: \
Composition( \
std::unique_ptr<domain::CoordinateMapBase<Frame::BlockLogical, \
Frame::Grid, DIM(data)>, \
std::default_delete<domain::CoordinateMapBase< \
Frame::BlockLogical, Frame::Grid, DIM(data)>>>, \
std::unique_ptr<domain::CoordinateMapBase< \
Frame::Grid, Frame::Inertial, DIM(data)>, \
std::default_delete<domain::CoordinateMapBase< \
Frame::Grid, Frame::Inertial, DIM(data)>>>); \
template domain::CoordinateMaps::Composition< \
brigand::list<Frame::BlockLogical, Frame::Grid, Frame::Distorted, \
Frame::Inertial>, \
DIM(data), std::integer_sequence<unsigned long, 0ul, 1ul, 2ul>>:: \
Composition( \
std::unique_ptr<domain::CoordinateMapBase<Frame::BlockLogical, \
Frame::Grid, DIM(data)>, \
std::default_delete<domain::CoordinateMapBase< \
Frame::BlockLogical, Frame::Grid, DIM(data)>>>, \
std::unique_ptr<domain::CoordinateMapBase< \
Frame::Grid, Frame::Distorted, DIM(data)>, \
std::default_delete<domain::CoordinateMapBase< \
Frame::Grid, Frame::Distorted, DIM(data)>>>, \
std::unique_ptr< \
domain::CoordinateMapBase<Frame::Distorted, Frame::Inertial, \
DIM(data)>, \
std::default_delete<domain::CoordinateMapBase< \
Frame::Distorted, Frame::Inertial, DIM(data)>>>); \
template domain::CoordinateMaps::Composition< \
brigand::list<Frame::ElementLogical, Frame::BlockLogical, \
Frame::Inertial>, \
Expand Down
5 changes: 1 addition & 4 deletions src/Domain/ElementDistribution.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,10 +63,7 @@ double get_num_points_and_grid_spacing_cost(
initial_extents, element_id, quadrature);
Element<Dim> element = ::domain::Initialization::create_initial_element(
element_id, block, initial_refinement_levels);
ElementMap<Dim, Frame::Grid> element_map{
element_id, block.is_time_dependent()
? block.moving_mesh_logical_to_grid_map().get_clone()
: block.stationary_map().get_to_grid_frame()};
ElementMap<Dim, Frame::Grid> element_map{element_id, block};
const tnsr::I<DataVector, Dim, Frame::ElementLogical> logical_coords =
logical_coordinates(mesh);
const tnsr::I<DataVector, Dim, Frame::Grid> grid_coords =
Expand Down
36 changes: 36 additions & 0 deletions src/Domain/ElementMap.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,10 @@

#include "Domain/ElementMap.hpp"

#include "Domain/CoordinateMaps/Composition.hpp"
#include "Domain/CoordinateMaps/CoordinateMap.hpp" // IWYU pragma: keep
#include "Domain/Structure/Side.hpp"
#include "Utilities/ErrorHandling/Assert.hpp"
#include "Utilities/GenerateInstantiations.hpp"
#include "Utilities/OptimizerHacks.hpp"

Expand Down Expand Up @@ -53,6 +55,40 @@ ElementMap<Dim, TargetFrame>::ElementMap(
jacobian_{map_slope_},
inverse_jacobian_{map_inverse_slope_} {}

// We could refactor the ElementMap class to use a `Composition` internally also
// for the element-to-block-logical map, so the whole map is just one
// composition.
template <size_t Dim, typename TargetFrame>
ElementMap<Dim, TargetFrame>::ElementMap(ElementId<Dim> element_id,
const Block<Dim>& block)
: ElementMap(
std::move(element_id),
[&block]() -> std::unique_ptr<domain::CoordinateMapBase<
Frame::BlockLogical, TargetFrame, Dim>> {
if constexpr (std::is_same_v<TargetFrame, Frame::Inertial>) {
if (block.is_time_dependent()) {
using CompositionType = domain::CoordinateMaps::Composition<
tmpl::list<Frame::BlockLogical, Frame::Grid,
Frame::Inertial>,
Dim>;
return std::make_unique<CompositionType>(
block.moving_mesh_logical_to_grid_map().get_clone(),
block.moving_mesh_grid_to_inertial_map().get_clone());
} else {
return block.stationary_map().get_clone();
}
} else if constexpr (std::is_same_v<TargetFrame, Frame::Grid>) {
if (block.is_time_dependent()) {
return block.moving_mesh_logical_to_grid_map().get_clone();
} else {
return block.stationary_map().get_to_grid_frame();
}
}
}()) {
ASSERT(element_id_.block_id() == block.id(),
"Element " << element_id_ << " is not in block " << block.id() << ".");
}

template <size_t Dim, typename TargetFrame>
void ElementMap<Dim, TargetFrame>::pup(PUP::er& p) {
p | block_map_;
Expand Down
39 changes: 30 additions & 9 deletions src/Domain/ElementMap.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,9 @@
#include <memory>

#include "DataStructures/Tensor/Tensor.hpp"
#include "Domain/Block.hpp"
#include "Domain/CoordinateMaps/CoordinateMap.hpp"
#include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
#include "Domain/Structure/ElementId.hpp"
#include "Utilities/Gsl.hpp"
#include "Utilities/MakeArray.hpp"
Expand Down Expand Up @@ -45,6 +47,14 @@ class ElementMap {
TargetFrame, Dim>>
block_map);

/// Construct from an `element_id` within the `block`. The (affine)
/// ElementLogical to BlockLogical map is determined by the `element_id`. The
/// BlockLogical to TargetFrame map is determined by the `block`:
/// - If the block is time-independent: the `block.stationary_map()` is used.
/// - If the block is time-dependent: The `block.moving_mesh_*_map()` maps
/// are used. Which maps are used depends on the TargetFrame.
ElementMap(ElementId<Dim> element_id, const Block<Dim>& block);

const domain::CoordinateMapBase<Frame::BlockLogical, TargetFrame, Dim>&
block_map() const {
return *block_map_;
Expand All @@ -54,17 +64,23 @@ class ElementMap {

template <typename T>
tnsr::I<T, Dim, TargetFrame> operator()(
const tnsr::I<T, Dim, Frame::ElementLogical>& source_point) const {
const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
const double time = std::numeric_limits<double>::signaling_NaN(),
const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
auto block_source_point =
apply_affine_transformation_to_point(source_point);
return block_map_->operator()(std::move(block_source_point));
return block_map_->operator()(std::move(block_source_point), time,
functions_of_time);
}

template <typename T>
tnsr::I<T, Dim, Frame::ElementLogical> inverse(
tnsr::I<T, Dim, TargetFrame> target_point) const {
tnsr::I<T, Dim, TargetFrame> target_point,
const double time = std::numeric_limits<double>::signaling_NaN(),
const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
auto block_source_point{
block_map_->inverse(std::move(target_point)).value()};
block_map_->inverse(std::move(target_point), time, functions_of_time)
.value()};
// Apply the affine map to the points
tnsr::I<T, Dim, Frame::ElementLogical> source_point;
for (size_t d = 0; d < Dim; ++d) {
Expand All @@ -77,11 +93,13 @@ class ElementMap {

template <typename T>
InverseJacobian<T, Dim, Frame::ElementLogical, TargetFrame> inv_jacobian(
const tnsr::I<T, Dim, Frame::ElementLogical>& source_point) const {
const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
const double time = std::numeric_limits<double>::signaling_NaN(),
const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
auto block_source_point =
apply_affine_transformation_to_point(source_point);
auto block_inv_jac =
block_map_->inv_jacobian(std::move(block_source_point));
auto block_inv_jac = block_map_->inv_jacobian(std::move(block_source_point),
time, functions_of_time);
InverseJacobian<T, Dim, Frame::ElementLogical, TargetFrame> inv_jac;
for (size_t d = 0; d < Dim; ++d) {
for (size_t i = 0; i < Dim; ++i) {
Expand All @@ -94,10 +112,13 @@ class ElementMap {

template <typename T>
Jacobian<T, Dim, Frame::ElementLogical, TargetFrame> jacobian(
const tnsr::I<T, Dim, Frame::ElementLogical>& source_point) const {
const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
const double time = std::numeric_limits<double>::signaling_NaN(),
const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
auto block_source_point =
apply_affine_transformation_to_point(source_point);
auto block_jac = block_map_->jacobian(std::move(block_source_point));
auto block_jac = block_map_->jacobian(std::move(block_source_point), time,
functions_of_time);
Jacobian<T, Dim, Frame::ElementLogical, TargetFrame> jac;
for (size_t d = 0; d < Dim; ++d) {
for (size_t i = 0; i < Dim; ++i) {
Expand Down
60 changes: 34 additions & 26 deletions src/Domain/FaceNormal.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -198,32 +198,40 @@ GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3),

#undef INSTANTIATION

#define INSTANTIATION(_, data) \
template void unnormalized_face_normal( \
const gsl::not_null< \
tnsr::i<DataVector, GET_DIM(data), Frame::Inertial>*> \
result, \
const Mesh<GET_DIM(data) - 1>& interface_mesh, \
const ElementMap<GET_DIM(data), Frame::Grid>& logical_to_grid_map, \
const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, \
GET_DIM(data)>& grid_to_inertial_map, \
const double time, \
const std::unordered_map< \
std::string, \
std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>& \
functions_of_time, \
const Direction<GET_DIM(data)>& direction); \
template tnsr::i<DataVector, GET_DIM(data), Frame::Inertial> \
unnormalized_face_normal( \
const Mesh<GET_DIM(data) - 1>& interface_mesh, \
const ElementMap<GET_DIM(data), Frame::Grid>& logical_to_grid_map, \
const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, \
GET_DIM(data)>& grid_to_inertial_map, \
const double time, \
const std::unordered_map< \
std::string, \
std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>& \
functions_of_time, \
#define INSTANTIATION(_, data) \
template void unnormalized_face_normal( \
const gsl::not_null< \
tnsr::i<DataVector, GET_DIM(data), Frame::Inertial>*> \
result, \
const Mesh<GET_DIM(data) - 1>& interface_mesh, \
const ElementMap<GET_DIM(data), Frame::Grid>& logical_to_grid_map, \
const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, \
GET_DIM(data)>& grid_to_inertial_map, \
const double time, \
const std::unordered_map< \
std::string, \
std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>& \
functions_of_time, \
const Direction<GET_DIM(data)>& direction); \
template void unnormalized_face_normal( \
const gsl::not_null< \
tnsr::i<DataVector, GET_DIM(data), Frame::Inertial>*> \
result, \
const Mesh<GET_DIM(data) - 1>& interface_mesh, \
const InverseJacobian<DataVector, GET_DIM(data), Frame::ElementLogical, \
Frame::Inertial>& inv_jacobian_on_interface, \
const Direction<GET_DIM(data)>& direction); \
template tnsr::i<DataVector, GET_DIM(data), Frame::Inertial> \
unnormalized_face_normal( \
const Mesh<GET_DIM(data) - 1>& interface_mesh, \
const ElementMap<GET_DIM(data), Frame::Grid>& logical_to_grid_map, \
const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, \
GET_DIM(data)>& grid_to_inertial_map, \
const double time, \
const std::unordered_map< \
std::string, \
std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>& \
functions_of_time, \
const Direction<GET_DIM(data)>& direction);

GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3))
Expand Down
4 changes: 4 additions & 0 deletions src/Domain/FunctionsOfTime/FunctionOfTime.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,4 +84,8 @@ class FunctionOfTime : public PUP::able {
WRAPPED_PUPable_abstract(FunctionOfTime); // NOLINT
};
} // namespace FunctionsOfTime

using FunctionsOfTimeMap = std::unordered_map<
std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>;

} // namespace domain
2 changes: 1 addition & 1 deletion src/Domain/FunctionsOfTime/FunctionsOfTimeAreReady.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ bool functions_of_time_are_ready_impl(
const Component* /*meta*/, const double time,
const std::optional<std::unordered_set<std::string>>& functions_to_check,
Args&&... args) {
if constexpr (Parallel::is_in_global_cache<Metavariables, CacheTag>) {
if constexpr (Parallel::is_in_mutable_global_cache<Metavariables, CacheTag>) {
const auto& proxy =
::Parallel::get_parallel_component<Component>(cache)[array_index];
const Parallel::ArrayComponentId array_component_id =
Expand Down
Loading
Loading