Changes (by Jordan Nicoules):

For reference, thread on the users mailing list: https://lists.einsteintoolkit.org/pipermail/users/2026-January/009852.html

The issue originates in trying to compute volume integrals in an analysis thorn, by using the volume form computed by Coordinates in a Llama grid.

‌

In simulations using Llama, different patch systems don’t compute the volume_form grid function in the same way (when requested by parameter store_volume_form. More precisely, the default behavior is (inverse-jacobian.F90):

volume_form(i,j,k) = detJ

‌

Meanwhile, patch system Thornburg04 has a dedicated routine because it needs to compute weights to account for overlapping patches. But even in the non-overlapping spherical part of the grid, it computes the volume form as (thornburg04.cc):

// set volume form to deterimant of Jacobian
    const CCTK_REAL det =  fabs((  J11[ijk] * J22[ijk--

Ticket URL: https://bitbucket.org/einsteintoolkit/tickets/issues/2912/inconsistent-computation-of-the-volume