Add ghost zones to MC evolution in Kerr test

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

@@ -184,13 +184,13 @@ tnsr::I<DataVector, 3, Frame::ElementLogical> spatial_coords_logical(
       for (cur_extents[2] = 0; cur_extents[2] < extents[2]; cur_extents[2]++) {
         const size_t storage_index = mesh.storage_index(cur_extents);
         x.get(0)[storage_index] =
-            -1.0 + 2.0 * static_cast<double>(cur_extents[0] + 0.5) /
+            -1.0 + 2.0 * (static_cast<double>(cur_extents[0]) + 0.5) /
                        static_cast<double>(extents[0]);
         x.get(1)[storage_index] =
-            -1.0 + 2.0 * static_cast<double>(cur_extents[1] + 0.5) /
+            -1.0 + 2.0 * (static_cast<double>(cur_extents[1]) + 0.5) /
                        static_cast<double>(extents[1]);
         x.get(2)[storage_index] =
-            -1.0 + 2.0 * static_cast<double>(cur_extents[2] + 0.5) /
+            -1.0 + 2.0 * (static_cast<double>(cur_extents[2]) + 0.5) /
                        static_cast<double>(extents[2]);
       }
     }
@@ -225,14 +225,17 @@ void test_evolve_kerr(const size_t mesh_size_1d) {
   // Set domain and coordintes
   const Mesh<3> mesh(mesh_size_1d, Spectral::Basis::FiniteDifference,
                      Spectral::Quadrature::CellCentered);
-  const DataVector zero_dv(mesh.number_of_grid_points(), 0.0);
-  const auto mesh_coordinates = spatial_coords_logical(mesh);
-  const auto inertial_coordinates = spatial_coords_inertial(mesh_coordinates);
   const size_t num_ghost_zones = 1;
   const size_t extended_mesh_1d_size = mesh_size_1d + 2 * num_ghost_zones;
   const Mesh<3> extended_mesh(extended_mesh_1d_size,
                               Spectral::Basis::FiniteDifference,
                               Spectral::Quadrature::CellCentered);
+  const size_t dv_size = mesh.number_of_grid_points();
+  const size_t extended_dv_size = extended_mesh.number_of_grid_points();
+  DataVector zero_dv(dv_size, 0.0);
+  DataVector extended_zero_dv(extended_dv_size, 0.0);
+  const auto mesh_coordinates = spatial_coords_logical(mesh);
+  const auto inertial_coordinates = spatial_coords_inertial(mesh_coordinates);
 
   // Mesh velocity set to std::null for now
   const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>> mesh_velocity =
@@ -330,18 +333,19 @@ void test_evolve_kerr(const size_t mesh_size_1d) {
   // Interaction rates
   const std::array<double, 2> energy_at_bin_center = {2.0, 5.0};
   std::array<std::array<DataVector, 2>, 2> absorption_opacity = {
-      std::array<DataVector, 2>{{zero_dv, zero_dv}},
-      std::array<DataVector, 2>{{zero_dv, zero_dv}}};
+      std::array<DataVector, 2>{{extended_zero_dv, extended_zero_dv}},
+      std::array<DataVector, 2>{{extended_zero_dv, extended_zero_dv}}};
   std::array<std::array<DataVector, 2>, 2> scattering_opacity = {
-      std::array<DataVector, 2>{{zero_dv, zero_dv}},
-      std::array<DataVector, 2>{{zero_dv, zero_dv}}};
+      std::array<DataVector, 2>{{extended_zero_dv, extended_zero_dv}},
+      std::array<DataVector, 2>{{extended_zero_dv, extended_zero_dv}}};
 
   Scalar<DataVector> coupling_tilde_tau =
-      make_with_value<Scalar<DataVector>>(zero_dv, 0.0);
+      make_with_value<Scalar<DataVector>>(extended_zero_dv, 0.0);
   Scalar<DataVector> coupling_rho_ye =
-      make_with_value<Scalar<DataVector>>(zero_dv, 0.0);
+      make_with_value<Scalar<DataVector>>(extended_zero_dv, 0.0);
   tnsr::i<DataVector, 3, Frame::Inertial> coupling_tilde_s =
-      make_with_value<tnsr::i<DataVector, 3, Frame::Inertial>>(zero_dv, 0.0);
+      make_with_value<tnsr::i<DataVector, 3, Frame::Inertial>>(extended_zero_dv,
+                                                               0.0);
 
   std::vector<Particles::MonteCarlo::Packet> packets{packet};
   Particles::MonteCarlo::TemplatedLocalFunctions<2, 2> MonteCarloStruct;
@@ -365,10 +369,10 @@ void test_evolve_kerr(const size_t mesh_size_1d) {
         mesh_velocity, inverse_jacobian, jacobian_inertial_to_fluid,
         inverse_jacobian_inertial_to_fluid);
     current_time += dt;
-    packets[0].renormalize_momentum(inverse_spatial_metric, lapse);
   }
   const double final_r = sqrt(pow(packets[0].coordinates.get(0) + 6.0, 2) +
    pow(packets[0].coordinates.get(1), 2)+pow(packets[0].coordinates.get(2),2));
+
   //Check r against geodesic evolution value obtained from python code
   CHECK(std::abs(final_r - 6.298490) < 1e-2);
 }
@@ -378,5 +382,5 @@ void test_evolve_kerr(const size_t mesh_size_1d) {
 SPECTRE_TEST_CASE("Unit.Evolution.Particles.MonteCarloEvolution",
                   "[Unit][Evolution]") {
   test_evolve_minkowski();
-  test_evolve_kerr(7);
+  test_evolve_kerr(15);
 }