Add InverseJacobianInertialToFluidCompute

src/DataStructures/Tensor/IndexType.hpp

@@ -41,6 +41,7 @@ struct FrameIsPhysical {};
 
 struct BlockLogical {};
 struct ElementLogical {};
+struct Fluid {};
 struct Grid {};
 struct Inertial : FrameIsPhysical {};
 struct Distorted {};
@@ -81,6 +82,10 @@ inline std::ostream& operator<<(std::ostream& os,
                                 const Frame::ElementLogical& /*meta*/) {
   return os << "ElementLogical";
 }
+inline std::ostream& operator<<(std::ostream& os,
+                                const Frame::Fluid& /*meta*/) {
+  return os << "Fluid";
+}
 inline std::ostream& operator<<(std::ostream& os,
                                 const Frame::Grid& /*meta*/) {
   return os << "Grid";
@@ -121,6 +126,11 @@ inline std::string prefix<Frame::ElementLogical>() {
   return "ElementLogical_";
 }
 
+template <>
+inline std::string prefix<Frame::Fluid>() {
+  return "Fluid_";
+}
+
 template <>
 inline std::string prefix<Frame::Grid>() {
   return "Grid_";

src/Evolution/Particles/MonteCarlo/CMakeLists.txt

@@ -9,6 +9,7 @@ spectre_target_sources(
   ${LIBRARY}
   PRIVATE
   EvolveMonteCarloPackets.cpp
+  InverseJacobianInertialToFluidCompute.cpp
   MonteCarloPacket.cpp
   )
 
@@ -17,6 +18,7 @@ spectre_target_headers(
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
   EvolveMonteCarloPackets.hpp
+  InverseJacobianInertialToFluidCompute.hpp
   MonteCarloPacket.hpp
   )
 

src/Evolution/Particles/MonteCarlo/InverseJacobianInertialToFluidCompute.cpp

@@ -0,0 +1,88 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Evolution/Particles/MonteCarlo/InverseJacobianInertialToFluidCompute.hpp"
+
+#include "PointwiseFunctions/GeneralRelativity/SpacetimeMetric.hpp"
+
+namespace Particles::MonteCarlo {
+
+void InverseJacobianInertialToFluidCompute::function(
+    gsl::not_null<return_type*> inv_jacobian,
+    const tnsr::I<DataVector, 3>& spatial_velocity,
+    const Scalar<DataVector>& lorentz_factor, const Scalar<DataVector>& lapse,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
+    const tnsr::ii<DataVector, 3>& spatial_metric) {
+  // The components of the inverse jacobian are just the vector components
+  // of the orthonormal tetrad comoving with the fluid in the inertial frame.
+
+  // First, inv_jacobian(a,0) is just u^a
+  inv_jacobian->get(0, 0) = get(lorentz_factor) / get(lapse);
+  for (size_t d = 0; d < 3; d++) {
+    inv_jacobian->get(d + 1, 0) =
+        get(lorentz_factor) *
+        (spatial_velocity.get(d) - shift.get(d) / get(lapse));
+  }
+
+  // Then, the other members of the tetrad are constructed using Gram-Shmidt
+
+  // Temporary memory allocation
+  DataVector temp_dot_product = get(lapse) * 0.0;
+  std::array<DataVector, 4> tetrad_vector{temp_dot_product, temp_dot_product,
+                                          temp_dot_product, temp_dot_product};
+  tnsr::aa<DataVector,3,Frame::Inertial> spacetime_metric =
+    gr::spacetime_metric(lapse, shift, spatial_metric);
+
+  // d = 1,2,3 for tetrad components built from x,y,z
+  for (size_t d = 1; d < 4; d++) {
+    // Base vector for Gram-Shmidt
+    for (size_t i = 0; i < 4; i++) {
+      gsl::at(tetrad_vector, i) = (i == d ? 1.0 : 0.0);
+    }
+    // Projection orthogonal to u^\mu
+    temp_dot_product = 0.0;
+    for (size_t i = 0; i < 3; i++) {
+      temp_dot_product +=
+          spatial_metric.get(d - 1, i) * spatial_velocity.get(i);
+    }
+    temp_dot_product *= get(lorentz_factor);
+    // Note: + sign when projecting here, because the 0th tetrad vector has norm
+    // -1.
+    for (size_t i = 0; i < 4; i++) {
+      gsl::at(tetrad_vector, i) += temp_dot_product * inv_jacobian->get(i, 0);
+    }
+
+    // Loop over other existing tetrad vectors to get orthogonal projection
+    for (size_t a = 0; a < d; a++) {
+      temp_dot_product = 0.0;
+      for (size_t b = 0; b < 4; b++) {
+        for (size_t c = 0; c < 4; c++) {
+          temp_dot_product += spacetime_metric.get(b,c) *
+                              gsl::at(tetrad_vector, b) *
+                              inv_jacobian->get(c, a);
+        }
+      }
+      for (size_t i = 0; i < 4; i++) {
+        gsl::at(tetrad_vector, i) -= temp_dot_product * inv_jacobian->get(i, a);
+      }
+    }
+
+    // Normalize tetrad vector
+    temp_dot_product = 0.0;
+    for (size_t a = 0; a < 4; a++) {
+      temp_dot_product += spacetime_metric.get(a,a) *
+                          gsl::at(tetrad_vector, a) * gsl::at(tetrad_vector, a);
+      for (size_t b = a + 1; b < 4; b++) {
+        temp_dot_product += 2.0 * spacetime_metric.get(a,b) *
+                            gsl::at(tetrad_vector, a) *
+                            gsl::at(tetrad_vector, b);
+      }
+    }
+    temp_dot_product = sqrt(temp_dot_product);
+    for (size_t a = 0; a < 4; a++) {
+      inv_jacobian->get(a, d) = gsl::at(tetrad_vector, a) / temp_dot_product;
+    }
+  }
+}
+
+}  // namespace Particles::MonteCarlo

src/Evolution/Particles/MonteCarlo/InverseJacobianInertialToFluidCompute.hpp

@@ -0,0 +1,44 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+/// \file
+/// Defines tags related to domain quantities
+
+#pragma once
+
+#include "Domain/Tags.hpp"
+#include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
+#include "PointwiseFunctions/Hydro/Tags.hpp"
+
+namespace Frame {
+struct Fluid;
+struct Inertial;
+}  // namespace Frame
+
+namespace Particles::MonteCarlo {
+
+// Inverse Jacobian of the map from inertial coordinate to an orthonormal frame
+// comoving with the fluid. That frame uses the 4-velocity as its time axis, and
+// constructs the other members of the tetrads using Gram-Shmidt's algorithm.
+struct InverseJacobianInertialToFluidCompute
+    : domain::Tags::InverseJacobian<4, Frame::Inertial, Frame::Fluid>,
+      db::ComputeTag {
+  using base = domain::Tags::InverseJacobian<4, typename Frame::Inertial,
+                                             typename Frame::Fluid>;
+  using return_type = typename base::type;
+  using argument_tags =
+      tmpl::list<hydro::Tags::SpatialVelocity<DataVector, 3, Frame::Inertial>,
+                 hydro::Tags::LorentzFactor<DataVector>,
+                 gr::Tags::Lapse<DataVector>,
+                 gr::Tags::Shift<DataVector, 3, Frame::Inertial>,
+                 gr::Tags::SpatialMetric<DataVector, 3> >;
+
+  static void function(gsl::not_null<return_type*> inv_jacobian,
+                const tnsr::I<DataVector, 3>& spatial_velocity,
+                const Scalar<DataVector>& lorentz_factor,
+                const Scalar<DataVector>& lapse,
+                const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
+                const tnsr::ii<DataVector, 3>& spatial_metric);
+};
+
+}  // namespace Particles::MonteCarlo

tests/Unit/Evolution/Particles/MonteCarlo/CMakeLists.txt

@@ -4,6 +4,7 @@
 set(LIBRARY "Test_MonteCarlo")
 
 set(LIBRARY_SOURCES
+  Test_InverseJacobianInertialToFluid.cpp
   Test_MonteCarloPacket.cpp
   )
 
@@ -15,5 +16,10 @@ add_test_library(
 target_link_libraries(
   ${LIBRARY}
   PRIVATE
+  DataStructures
+  GeneralRelativity
+  GeneralRelativityHelpers
+  Hydro
+  HydroHelpers
   MonteCarlo
   )

tests/Unit/Evolution/Particles/MonteCarlo/Test_InverseJacobianInertialToFluid.cpp

@@ -0,0 +1,86 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Framework/TestingFramework.hpp"
+
+#include <limits>
+
+#include "DataStructures/DataBox/DataBox.hpp"
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Domain/Tags.hpp"
+#include "Evolution/Particles/MonteCarlo/InverseJacobianInertialToFluidCompute.hpp"
+#include "Framework/TestHelpers.hpp"
+#include "Helpers/DataStructures/DataBox/TestHelpers.hpp"
+#include "Helpers/PointwiseFunctions/GeneralRelativity/TestHelpers.hpp"
+#include "Helpers/PointwiseFunctions/Hydro/TestHelpers.hpp"
+#include "PointwiseFunctions/GeneralRelativity/SpacetimeMetric.hpp"
+#include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
+#include "PointwiseFunctions/Hydro/Tags.hpp"
+#include "Utilities/Gsl.hpp"
+
+SPECTRE_TEST_CASE(
+    "Unit.Evolution.Particles.MonteCarlo.InverseJacobianInertialToFluid",
+    "[Unit][Evolution]") {
+  const DataVector used_for_size(5);
+  MAKE_GENERATOR(generator);
+
+  const auto lapse =
+      TestHelpers::gr::random_lapse(make_not_null(&generator), used_for_size);
+  const auto shift = TestHelpers::gr::random_shift<3>(make_not_null(&generator),
+                                                      used_for_size);
+  const auto spatial_metric = TestHelpers::gr::random_spatial_metric<3>(
+      make_not_null(&generator), used_for_size);
+  const auto lorentz_factor = TestHelpers::hydro::random_lorentz_factor(
+      make_not_null(&generator), used_for_size);
+  const auto spatial_velocity =
+      TestHelpers::hydro::random_velocity<DataVector, 3>(
+          make_not_null(&generator), lorentz_factor, spatial_metric);
+
+  TestHelpers::db::test_compute_tag<
+      Particles::MonteCarlo::InverseJacobianInertialToFluidCompute>(
+      "InverseJacobian(Inertial,Fluid)");
+  const auto box = db::create<
+      db::AddSimpleTags<
+          hydro::Tags::LorentzFactor<DataVector>,
+          hydro::Tags::SpatialVelocity<DataVector, 3, Frame::Inertial>,
+          gr::Tags::Lapse<DataVector>,
+          gr::Tags::Shift<DataVector, 3, Frame::Inertial>,
+          gr::Tags::SpatialMetric<DataVector, 3, Frame::Inertial>>,
+      db::AddComputeTags<
+          Particles::MonteCarlo::InverseJacobianInertialToFluidCompute>>(
+      lorentz_factor, spatial_velocity, lapse, shift, spatial_metric);
+
+  const auto& inverse_jacobian =
+      db::get<Particles::MonteCarlo::InverseJacobianInertialToFluidCompute>(
+          box);
+  const tnsr::aa<DataVector, 3, Frame::Inertial> spacetime_metric =
+      gr::spacetime_metric(lapse, shift, spatial_metric);
+
+  // Check that the time vector is u^mu
+  CHECK_ITERABLE_APPROX(inverse_jacobian.get(0, 0),
+                        get(lorentz_factor) / get(lapse));
+  for (size_t d = 0; d < 3; d++) {
+    CHECK_ITERABLE_APPROX(inverse_jacobian.get(d + 1, 0),
+                          get(lorentz_factor) * (spatial_velocity.get(d) -
+                                                 shift.get(d) / get(lapse)));
+  }
+
+  // Test that we have orthonormal tetrads
+  DataVector dot_product(used_for_size);
+  DataVector expected_dot_product(used_for_size);
+  for (size_t a = 0; a < 4; a++) {
+    for (size_t b = a; b < 4; b++) {
+      dot_product = 0.0;
+      for (size_t d = 0; d < 4; d++) {
+        for (size_t dd = 0; dd < 4; dd++) {
+          dot_product += spacetime_metric.get(d, dd) *
+                         inverse_jacobian.get(d, a) *
+                         inverse_jacobian.get(dd, b);
+        }
+      }
+      expected_dot_product = (a == b ? (a == 0 ? -1.0 : 1.0) : 0.0);
+      CHECK_ITERABLE_APPROX(dot_product, expected_dot_product);
+    }
+  }
+}