Report error in the trace of particle density matrices

[dwillcox]

We can report the L2 and max norm of the error in the trace of the particles' density matrices.

Currently, for each particle the flavor vector length L is calculated from N_ab, the matrix storing the number of neutrinos in each state. We also have the particle weight N_p which is just the number of neutrinos in each particle.

Error in Tr(rho) = 1 can significantly affect the number of neutrinos that transform if L ~ N_p, since the size of L controls the maximum amount of flavor transformation. We can define an error E = |L(rho) * N_p - L_0| where L_0 is the initial length stored in the particles and L(rho) * N_p is the flavor vector length calculated from rho, scaled by the particle weight. Then take the L2 and max norms of E by computing Ep^2 for each particle and depositing it onto the grid E, then use sum and max to get L2 and max norms on the grid.

We can then add an input parameter for the error reporting frequency, like report_error_every.

The idea is that this will tell us when our simulations would benefit from symbolically making the assumption that Tr(rho_ab) = 1 => Tr(N_ab) = N_p, eliminating one of the diagonal components algebraically, and evolving a traceless matrix N_ab_p for particles.

[srichers]

I think there will be a problem when the error in Tr(rho) is comparable to the flavor vector length of rho, rather than when Tr(rho) is comparable to the flavor vector length of rho.

[dwillcox]

Ah, got it, that's simpler to calculate :)

[dwillcox]

So that would just be E_p = |Tr(rho) - 1| then, the way rho is currently defined, I think ...

[srichers]

Yes, that's right. There is also the error in the length, so E2_p = |length - old_length|

[srichers]

The realtime_reduced_data branch reports a sum of Tr(rho) over the entire domain. Initially it is equal to the number of particles and should remain so throughout the calculation.