Enable time dependent translation with transition for Sphere domain #5576

guilara · 2023-10-19T13:51:01Z

Proposed changes

Add a translation map acting on the outer most shell of the spherical domain. If there is no radial partition, apply a uniform translation (since the inner boundary of the translation needs to be spherical and the shape map could also be active).

Note: The test Unit.Domain.Creators.Sphere failed because of roundoff error in assets

ASSERT(max(magnitude(result)) <= outer_radius_.value(), "Coordinates translated outside outer radius, This should not " "happen");

in

spectre/src/Domain/CoordinateMaps/TimeDependent/Translation.cpp

Line 376 in 06b2395

ASSERT(max(magnitude(result)) <= outer_radius_.value(),

When I printed the values during a test, the assert failed by a factor of order 10e-15 of the outer_radius_ and seemed to remain of the same order in time. Changed the assert to allow for roundoff error.

Update: issue above fixed in #5601

Upgrade instructions

Code review checklist

The code is documented and the documentation renders correctly. Run
make doc to generate the documentation locally into BUILD_DIR/docs/html.
Then open index.html.
The code follows the stylistic and code quality guidelines listed in the
code review guide.
The PR lists upgrade instructions and is labeled bugfix or
new feature if appropriate.

Further comments

nilsvu · 2023-10-19T16:59:01Z

src/Domain/Creators/SphereTimeDependentMaps.cpp

    : initial_time_(initial_time),
      initial_l_max_(shape_map_options.l_max),
-      initial_shape_values_(shape_map_options.initial_values) {}
+      initial_shape_values_(shape_map_options.initial_values),
+      translation_velocity_(translation_velocity),


initial_translation_velocity

nilsvu · 2023-10-19T17:33:41Z

src/Domain/Creators/SphereTimeDependentMaps.hpp

@@ -131,15 +139,47 @@ struct TimeDependentMapOptions {
    std::optional<std::variant<KerrSchildFromBoyerLindquist>> initial_values{};
  };

-  using options = tmpl::list<InitialTime, ShapeMapOptions>;
+  struct TranslationMapVelocity {
+    static std::string name() { return "Velocity"; }


InitialVelocity

nilsvu · 2023-10-19T17:54:36Z

src/Domain/Creators/SphereTimeDependentMaps.hpp

+    using group = TranslationMapOptions;
+  };
+
+  struct TranslationMapInnerRadius {


Don't take the inner and outer radii as options. Just apply the transition in the outer shell of the sphere, so inner radius is the last number in radial_partitions and outer radius is the outer radius of the sphere. Add this translation map with the transition only to the blocks that make up the outer shell. Apply a rigid translation (no inner/outer radii) to all other blocks. If there is no partition (just one shell) and the interior of the sphere is filled with a cube then either error out or apply a rigid translation (because the transition only works with spherical inner and outer radii).

nilsvu · 2023-10-20T18:15:16Z

src/Domain/Creators/SphereTimeDependentMaps.cpp

@@ -115,17 +133,21 @@ void TimeDependentMapOptions::build_maps(const std::array<double, 3>& center,
  shape_map_ = ShapeMap{center,         initial_l_max_,
                        initial_l_max_, std::move(transition_func),
                        shape_name,     size_name};
+  if (translation_inner_radius == translation_outer_radius) {


Take the translation radii as std::optional<std::pair<double, double>> then you can either pass std::nullopt or both radii.

You also have to build a second rigid translation map with no radii, which you can apply to the inner blocks. Then you don't need to remove asserts in the Translation map. I think it is necessary to have these two different maps because their Jacobian is discontinuous at the translation_inner_radius.

nilsvu · 2023-10-20T18:18:12Z

src/Domain/Creators/SphereTimeDependentMaps.cpp

+      std::make_unique<FunctionsOfTime::PiecewisePolynomial<2>>(
+          initial_time_,
+          std::array<DataVector, 3>{
+              {{3, 0.0}, initial_velocity_temp, {3, 0.0}}},


std::move(initial_velocity_temp)

nilsvu · 2023-10-20T18:23:28Z

src/Domain/Creators/SphereTimeDependentMaps.cpp

-    return std::make_unique<
-        IdentityForComposition<Frame::Distorted, Frame::Inertial>>(
-        IdentityMap{});
+    return std::make_unique<DistortedToInertialComposition>(translation_map_);


Here you add the rigid translation (not the transition)

nilsvu · 2023-10-20T18:29:36Z

src/Domain/Creators/SphereTimeDependentMaps.cpp

  } else {
-    return std::make_unique<
-        IdentityForComposition<Frame::Grid, Frame::Inertial>>(IdentityMap{});
+    return std::make_unique<GridToInertialSimple>(translation_map_);


Take another bool like transition_translation_map to determine if this block should use the transition or the rigid translation.

guilara · 2023-10-21T13:00:03Z

@nilsvu I think addressed all the comments

knelli2 · 2023-11-01T18:39:40Z

src/Domain/Creators/SphereTimeDependentMaps.hpp

+  struct InitialVelocity {
+    static std::string name() { return "InitialVelocity"; }


We probably also want to be able to specify the initial translation itself, not just the initial velocity, so we'll need an option for that. To keep things organized, make the translation options like the shape options. One struct that has

struct TranslationMapOptions { using type = TranslationMapOptions; static std::string name() { return "TranslationMap"; } struct InitialValues { using type = std::array<std::array<double, 3>, 2>; ... }; ... using options = tmpl::list<InitialValues>; std::array<std::array<double, 3>, 2> initial_values{}; };

knelli2 · 2023-11-01T18:41:18Z

src/Domain/Creators/SphereTimeDependentMaps.hpp

+  TranslationMap rigid_translation_map_{};
+  TranslationMap translation_map_{};


Use linear_falloff_translation_map_ instead of translation_map_ so it's clear. Or just make it a std::pair of translation maps.

knelli2 · 2023-11-01T18:46:23Z

src/Domain/Creators/Sphere.cpp

+      // The translation map linear transition occurs only in the outer-most
+      // shell, such that the outer boundary is not translated. If there
+      // is no radial partition, a uniform translation is applied instead, as
+      // the transition requires spherical shape.


I don't think it makes sense to allow only a single spherical shell while also having a translation. I don't think we have boundary conditions that can handle that. We should probably just error if that happens and require we have two spherical shells (one radial partition).

knelli2 · 2023-11-01T18:46:40Z

src/Domain/Creators/Sphere.cpp

+                      radial_partitioning_.empty()
+                          ? std::nullopt
+                          : std::optional<std::pair<double, double>>{
+                                {radial_partitioning_.back(), outer_radius_}});


Then these don't have to be optional, but can always be passed.

knelli2 · 2023-11-01T18:47:27Z

src/Domain/Creators/SphereTimeDependentMaps.cpp

+
+  if (translation_transition_radii.has_value()) {


And then here you just always have a rigid translation map and a falloff translation map.

knelli2 · 2023-11-01T18:49:41Z

tests/Unit/Domain/Creators/Test_Sphere.cpp

          "an outflow-type boundary condition, you must use that."));
  using TDMO = domain::creators::sphere::TimeDependentMapOptions;


The sphere test failed for one of the builds still.

guilara · 2023-11-08T17:01:35Z

For some reason I am still getting the asserts (in the description) triggered with the Xcts Kerr input file. In particular for the inner shell. Maybe there is still a bug somewhere here or in Translation.

As for the rest, I implemented the changes requested @knelli2

knelli2

The fixups look good. You can squash the existing changes plus these new ones. As for the error in the translation map, if the values are really off by 1e-14 then we just need to account for roundoff. Basically do this in two places in the translation map

diff --git a/src/Domain/CoordinateMaps/TimeDependent/Translation.cpp b/src/Domain/CoordinateMaps/TimeDependent/Translation.cpp
index de3f73ed2d..debb09ef7a 100644
--- a/src/Domain/CoordinateMaps/TimeDependent/Translation.cpp
+++ b/src/Domain/CoordinateMaps/TimeDependent/Translation.cpp
@@ -373,7 +373,8 @@ std::array<tt::remove_cvref_wrap_t<T>, Dim> Translation<Dim>::piecewise_helper(
         get_element(gsl::at(result, i), k) +=
             gsl::at(func_or_deriv_of_time, i) * radial_falloff_factor;
       }
-      ASSERT(max(magnitude(result)) <= outer_radius_.value(),
+      const double eps = std::numeric_limits<double>::epsilon() * 100.0;
+      ASSERT(max(magnitude(result)) - eps <= outer_radius_.value(),
              "Coordinates translated outside outer radius, This should not "
              "happen");
     }

If that doesn't work, then we'll need to investigate more. While you're at it, please also print the numbers in that ASSERT message so we can see how off they are.

Also, please fix all clang-tidy errors when squashing.

knelli2 · 2023-11-08T21:18:45Z

src/Domain/Creators/Sphere.cpp

+        const bool use_rigid_translation =
+            block_id + num_blocks_per_shell_ < num_blocks_;


I think this is incorrect if fill_interior_ is true since the central cube is append at the end of the block list

guilara · 2023-11-13T14:04:53Z

The fixups look good. You can squash the existing changes plus these new ones. As for the error in the translation map, if the values are really off by 1e-14 then we just need to account for roundoff. Basically do this in two places in the translation map
diff --git a/src/Domain/CoordinateMaps/TimeDependent/Translation.cpp b/src/Domain/CoordinateMaps/TimeDependent/Translation.cpp
index de3f73ed2d..debb09ef7a 100644
--- a/src/Domain/CoordinateMaps/TimeDependent/Translation.cpp
+++ b/src/Domain/CoordinateMaps/TimeDependent/Translation.cpp
@@ -373,7 +373,8 @@ std::array<tt::remove_cvref_wrap_t<T>, Dim> Translation<Dim>::piecewise_helper(
         get_element(gsl::at(result, i), k) +=
             gsl::at(func_or_deriv_of_time, i) * radial_falloff_factor;
       }
-      ASSERT(max(magnitude(result)) <= outer_radius_.value(),
+      const double eps = std::numeric_limits<double>::epsilon() * 100.0;
+      ASSERT(max(magnitude(result)) - eps <= outer_radius_.value(),
              "Coordinates translated outside outer radius, This should not "
              "happen");
     }
If that doesn't work, then we'll need to investigate more. While you're at it, please also print the numbers in that ASSERT message so we can see how off they are.

Also, please fix all clang-tidy errors when squashing.

Running the Xcts input file with one partition I get:

For the assert checking blocks outside the inner radius I get:
equal within roundoff (max(magnitude(result)), outer_radius_.value())) = 1 1.0e15*(outer_radius_.value() - max(magnitude(result)))/outer_radius_.value() = -0.177636
or -0.177636e-15

whereas
2) for the assert checking blocks that have at least one point as falling within the inner radius I get
< inner_radius_value()... get_element(gsl::at(result, index), k) = -0.752477 , get_element(gsl::at(result, index), k) = -2.053507 , get_element(gsl::at(result, index), k) = 2.053507 , get_element(radius, k) = 3.000000 , equal within roundoff (max(magnitude(result)), outer_radius_.value())) = 1 , 1.0e15*(outer_radius_.value() - max(magnitude(result)))/outer_radius_.value() = -0.177636

which is actually the same number -0.177636e-15

Add options Add rigid translation Fix clang tidy Fix tests Rename translation map member Reorganize translation map options Add initial translation parameter Fix tests Forbid single shell with hard coded time dep maps Fix Fix Allow for roundoff in translation asserts

guilara · 2023-11-13T20:17:49Z

@nilsvu Would it be ok if I set the spin to zero and remove the time dep maps in the Xcts KerrSchild test? In this PR I require at least two shells for hardcoded time dep maps and that test is timing out.

Or how much can I increase the timeout time?

knelli2 · 2023-12-06T01:28:10Z

@nilsvu Can you take another look at this?

knelli2 · 2024-03-19T06:12:20Z

@guilara I believe we can close this because it is superseded by #5758?

knelli2 · 2024-06-11T00:22:32Z

Closing because this feature has already been implemented in #5758.

nilsvu requested changes Oct 19, 2023

View reviewed changes

nilsvu requested changes Oct 20, 2023

View reviewed changes

guilara force-pushed the TimeDepMaps-Translation-Sphere branch from fff292d to f090181 Compare October 21, 2023 10:13

guilara changed the title ~~Add translation to sphere domain time maps~~ Enable time dependent translation with transition for Sphere domain Oct 21, 2023

guilara force-pushed the TimeDepMaps-Translation-Sphere branch from f090181 to 0f3c17f Compare October 21, 2023 12:10

guilara marked this pull request as ready for review October 21, 2023 12:57

knelli2 self-requested a review October 26, 2023 16:30

guilara force-pushed the TimeDepMaps-Translation-Sphere branch from c842318 to e008c93 Compare October 27, 2023 08:45

guilara mentioned this pull request Nov 1, 2023

Add boost to KerrSchild #5610

Merged

3 tasks

knelli2 requested changes Nov 1, 2023

View reviewed changes

guilara force-pushed the TimeDepMaps-Translation-Sphere branch 4 times, most recently from ff9f425 to 4c11784 Compare November 8, 2023 13:10

knelli2 requested changes Nov 8, 2023

View reviewed changes

guilara force-pushed the TimeDepMaps-Translation-Sphere branch from 4c11784 to 046741a Compare November 13, 2023 10:41

guilara force-pushed the TimeDepMaps-Translation-Sphere branch from 046741a to 066f010 Compare November 13, 2023 15:19

guilara force-pushed the TimeDepMaps-Translation-Sphere branch from 066f010 to b781510 Compare November 13, 2023 17:41

Increase timeout for Xcts KerrSchild test

8567967

nilsvu self-requested a review December 9, 2023 01:10

knelli2 closed this Jun 11, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Enable time dependent translation with transition for Sphere domain #5576

Enable time dependent translation with transition for Sphere domain #5576

guilara commented Oct 19, 2023 •

edited

Loading

nilsvu Oct 19, 2023

nilsvu Oct 19, 2023

nilsvu Oct 19, 2023

nilsvu Oct 20, 2023

nilsvu Oct 20, 2023

nilsvu Oct 20, 2023

nilsvu Oct 20, 2023

guilara commented Oct 21, 2023

knelli2 Nov 1, 2023

knelli2 Nov 1, 2023

knelli2 Nov 1, 2023

knelli2 Nov 1, 2023

knelli2 Nov 1, 2023

knelli2 Nov 1, 2023

guilara commented Nov 8, 2023 •

edited

Loading

knelli2 left a comment

knelli2 Nov 8, 2023

guilara commented Nov 13, 2023 •

edited

Loading

guilara commented Nov 13, 2023

knelli2 commented Dec 6, 2023

knelli2 commented Mar 19, 2024

knelli2 commented Jun 11, 2024

		struct InitialVelocity {
		static std::string name() { return "InitialVelocity"; }

		TranslationMap rigid_translation_map_{};
		TranslationMap translation_map_{};

		"an outflow-type boundary condition, you must use that."));
		using TDMO = domain::creators::sphere::TimeDependentMapOptions;

		const bool use_rigid_translation =
		block_id + num_blocks_per_shell_ < num_blocks_;

Enable time dependent translation with transition for Sphere domain #5576

Enable time dependent translation with transition for Sphere domain #5576

Conversation

guilara commented Oct 19, 2023 • edited Loading

Proposed changes

Upgrade instructions

Code review checklist

Further comments

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

guilara commented Oct 21, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

guilara commented Nov 8, 2023 • edited Loading

knelli2 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

guilara commented Nov 13, 2023 • edited Loading

guilara commented Nov 13, 2023

knelli2 commented Dec 6, 2023

knelli2 commented Mar 19, 2024

knelli2 commented Jun 11, 2024

guilara commented Oct 19, 2023 •

edited

Loading

guilara commented Nov 8, 2023 •

edited

Loading

guilara commented Nov 13, 2023 •

edited

Loading