I have open-sourced a very useful set of diagnostic routines for GRMHD, GRHD, and even vacuum simulations.
The thorns are available here: https://bitbucket.org/zach_etienne/wvuthorns_diagnostics
Here's a brief description of each thorn: Brief description of included thorns
'''Seed_Magnetic_Fields-modified''' Extended Seed_Magnetic_Fields thorn for binary neutron stars. Supercedes Seed_Magnetic_Fields thorn.
'''Meudon_Bin_NS-modified''' Modifications to Meudon BNS initial data thorn to disable the overwriting of initial lapse/shift, which acts to significantly reduce coordinate eccentricity. Supercedes Meudon_Bin_NS thorn.
'''VolumeIntegrals_GRMHD''' Nice GRMHD volume integration thorn, currently depends on IllinoisGRMHD and Carpet. Performs volume integrals on arbitrary "Swiss-cheese"-like topologies, and even interoperates with Carpet to track NS centers of mass.
'''VolumeIntegrals_vacuum''' Nice GRMHD volume integration thorn, currently depends on ML_BSSN. There is a bit of code duplication and duplicated functionality between VI_GRMHD and VI_vacuum, to ensure that VI_vacuum can be used without enabling a GRMHD code. There is probably a better way of doing this, but I haven't had the time to think deeply about this.
'''particle_tracerET''' Solves the ODE \partial_t x^i = v^i for typically thousands of tracer particles, using an RK4 integration atop the current timestepping. E.g., one RK4 substep in the particle integration might occur every 16 RK4 substeps in the GRMHD evolution. These tracer particle positions are quite useful for visualizing magnetic field lines in a consistent way from frame-to-frame in a movie (recall that in the GRMHD approximation, the magnetic field lines stay attached the fluid elements they thread ["Flux Freezing"]). Note that the velocity must be consistent with the velocity appearing in the GRMHD induction equation. This thorn reads in the HydroBase vel vector gridfunction, which assumes the Valencia formalism, and converts it into the induction equation velocity.
'''smallbPoynET''' Computes b^i^, b^2^, and three spatial components of Poynting flux. It also computes (-1-u,,0,,), which is useful for tracking unbound matter.