Print diagnostics info

src/Evolution/Particles/MonteCarlo/EvolvePackets.cpp

@@ -5,6 +5,7 @@
 
 #include "Evolution/Particles/MonteCarlo/Packet.hpp"
 #include "Evolution/Particles/MonteCarlo/Scattering.hpp"
+#include "Parallel/Printf/Printf.hpp"
 
 namespace Particles::MonteCarlo {
 
@@ -61,6 +62,13 @@ void evolve_single_packet_on_geodesic(
   // Calculate time derivative of 3-momentum at beginning of time step
   detail::time_derivative_momentum_geodesic(&dpdt, *packet, lapse, d_lapse,
                                             d_shift, d_inv_spatial_metric);
+  Parallel::printf("Initial coord: %.5f %.5f %.5f\n", packet->coordinates[0],
+                   packet->coordinates[1], packet->coordinates[2]);
+  Parallel::printf("Initial momentum: %.5f %.5f %.5f\n", packet->momentum[0],
+                   packet->momentum[1], packet->momentum[2]);
+  Parallel::printf("dpdt: %.5f %.5f %.5f\n", gsl::at(dpdt, 0), gsl::at(dpdt, 1),
+                   gsl::at(dpdt, 2));
+
   // Take half-step (momentum only, as time derivative is independent of
   // position)
   for (size_t i = 0; i < 3; i++) {
@@ -74,6 +82,10 @@ void evolve_single_packet_on_geodesic(
   // Calculate time derivative of 3-momentum and position at half-step
   detail::time_derivative_momentum_geodesic(&dpdt, *packet, lapse, d_lapse,
                                             d_shift, d_inv_spatial_metric);
+  Parallel::printf("Half step momentum: %.5f %.5f %.5f\n", packet->momentum[0],
+                   packet->momentum[1], packet->momentum[2]);
+  Parallel::printf("dpdt: %.5f %.5f %.5f\n", gsl::at(dpdt, 0), gsl::at(dpdt, 1),
+                   gsl::at(dpdt, 2));
   for (size_t i = 0; i < 3; i++) {
     gsl::at(dxdt_inertial, i) = (-1.0) * shift.get(i)[closest_point_index];
     if (mesh_velocity.has_value()) {
@@ -86,19 +98,29 @@ void evolve_single_packet_on_geodesic(
           packet->momentum[j] / packet->momentum_upper_t;
     }
   }
+  Parallel::printf("dxdt_inert: %.5f %.5f %.5f\n", gsl::at(dxdt_inertial, 0),
+                   gsl::at(dxdt_inertial, 1), gsl::at(dxdt_inertial, 2));
+
   for (size_t i = 0; i < 3; i++) {
     for (size_t j = 0; j < 3; j++) {
       gsl::at(dxdt_logical, i) +=
           gsl::at(dxdt_inertial, j) *
           inverse_jacobian_logical_to_inertial.get(i, j)[closest_point_index];
     }
   }
+  Parallel::printf("dxdt_inert: %.5f %.5f %.5f\n", gsl::at(dxdt_logical, 0),
+                   gsl::at(dxdt_logical, 1), gsl::at(dxdt_logical, 2));
+
   // Take full time step
   for (size_t i = 0; i < 3; i++) {
     packet->momentum[i] = gsl::at(p0, i) + gsl::at(dpdt, i) * time_step;
     packet->coordinates[i] += gsl::at(dxdt_logical, i) * time_step;
   }
   packet->time += time_step;
+  Parallel::printf("Final coord: %.5f %.5f %.5f\n", packet->coordinates[0],
+                   packet->coordinates[1], packet->coordinates[2]);
+  Parallel::printf("Final momentum: %.5f %.5f %.5f\n", packet->momentum[0],
+                   packet->momentum[1], packet->momentum[2]);
 }
 
 }  // namespace Particles::MonteCarlo

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

@@ -195,9 +195,9 @@ void test_evolve_kerr() {
   gr::Solutions::KerrSchild solution(mass, spin, center);
 
   // Set domain and coordintes
-  const Mesh<3> mesh(2, Spectral::Basis::FiniteDifference,
+  const Mesh<3> mesh(3, Spectral::Basis::FiniteDifference,
                      Spectral::Quadrature::CellCentered);
-  const DataVector zero_dv(mesh.number_of_grid_points());
+  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);
   // Mesh velocity set to std::null for now
@@ -247,7 +247,7 @@ void test_evolve_kerr() {
       make_with_value<tnsr::i<DataVector, 3, Frame::Inertial>>(lapse, 0.0);
 
   // Initialize packet on the x-axis, moving in the y-direction
-  Particles::MonteCarlo::Packet packet(4, 1.0, 0, 0.0, 6.5, 0.0, 0.0, 1.0, 0.0,
+  Particles::MonteCarlo::Packet packet(0, 1.0, 5, 0.0, 6.5, 0.0, 0.0, 1.0, 0.0,
                                        1.0, 0.0);
   packet.renormalize_momentum(inverse_spatial_metric, lapse);
 
@@ -305,9 +305,9 @@ void test_evolve_kerr() {
       &coupling_rho_ye, final_time, mesh, mesh_coordinates,
       absorption_opacity, scattering_opacity, energy_at_bin_center,
       lorentz_factor, lower_spatial_four_velocity, lapse, shift, deriv_lapse,
-      deriv_shift, deriv_inverse_spatial_metric, inverse_spatial_metric,
-      mesh_velocity, inverse_jacobian, jacobian_inertial_to_fluid,
-      inverse_jacobian_inertial_to_fluid);
+      deriv_shift, deriv_inverse_spatial_metric, spatial_metric,
+      inverse_spatial_metric, mesh_velocity, inverse_jacobian,
+      jacobian_inertial_to_fluid, inverse_jacobian_inertial_to_fluid);
 
   Parallel::printf("%.5f %.5f %.5f %.5f \n", packets[0].time,
                    packets[0].coordinates.get(0), packets[0].coordinates.get(1),