Inconsistent units across multiple simulations

[LiamPattinson]

In Issue https://github.com/pyro-kinetics/pyrokinetics/issues/326, we found that calls to _build_cache in functions such as SimulationNormalisation.set_lref were imposing a significant performance penalty. It's necessary to call this as these functions redefine units, and without _build_cache, pint won't recognise the change and will continue using the old units definitions.

I've found another reason beyond performance concerns to rethink this approach:

def test_pyro_consistent_units():
    """Test that after loading a new geometry, old values are still valid."""

    # Create local geometry from G-EQDSK file
    pyro = Pyro(gk_file=gk_templates["GS2"], gk_type="GS2")
    pyro.load_global_eq(eq_templates["GEQDSK"], eq_type="GEQDSK")
    pyro.load_local_geometry(psi_n=0.5)

    # Take reference copies of some quantity that depends on lref
    R = pyro.local_geometry.R
    R_copy_1 = R.copy()
    R_meter_1 = R.to("meter")

    # Load in a different equilibrium file and generate a new geometry
    pyro.load_global_eq(eq_templates["TRANSP"], eq_type="TRANSP")
    pyro.load_local_geometry(psi_n=0.95)

    # Make new reference copies of the _old_ quantity
    # This should not have changed!
    R_copy_2 = R.copy()
    R_meter_2 = R.to("meter")

    # Test that our quantities are unchanged
    np.testing.assert_array_equal(R_copy_1.magnitude, R_copy_2.magnitude)
    assert R_copy_1.units == R_copy_2.units
    np.testing.assert_array_equal(R_meter_1.magnitude, R_meter_2.magnitude)

What we find is that the first two assertions pass, as the numbers and unit name of R haven't been modified. However, the third assertion fails, as the definition of lref_minor_radius0000 (or whatever it gets called) has changed behind the scenes, and therefore R has been invalidated. I'm not sure if this is currently causing any internal issues with Pyrokinetics, but it is a surprising and unintuitive behaviour.

---

To solve this, we could modify functions such as Pyro.load_* so that they always generate a new norms object, as after calling those functions we're always looking at a brand new simulation anyway. We'd also have to modify LocalGeometry.from_global_eq etc to do this.

I've personally been thinking about a solution where we introduce a new class with a name like LocalGKModel that binds together a LocalGeometry, LocalSpecies, Normalisation, and Numerics, and it would be immutable. It could have functions like .with_new_geometry() that would return a brand new LocalGKModel with all new components, including a new Normalisation. Its job would be to ensure that modifications to an existing simulation would always be self-consistent (if lref were to change, then all three of LocalGeometry, LocalSpecies and Normalisation may need to be updated to match), while ensuring that existing simulations aren't affected. By taking on these roles, it should be possible to reduce the burden of managing units from other classes throughout the project.

Pyro could then manage LocalGKModel instances instead of handling their constituent parts directly, and I think it would make parameter scans a bit easier. It could also be the object returned by GKInput classes instead of them effectively returning themselves (which causes some headaches for the file-reading bits of the code), and be the object that GKOutput uses to figure out the input parameters of a run.

This is all just conceptual at this stage though, and I'm sure an actual implementation would be quite messy. An alternative to implementing the above would be to modify Pyro to achieve the same outcome (becoming immutable, function calls returning new Pyro instances), but I'm wary of this breaking existing scripts.

---

As a new release is imminent (Issue https://github.com/pyro-kinetics/pyrokinetics/issues/316), I think it would be best to hold off on this until the next release. My long-winded solution would take quite a lot of work, so if an easier solution presents itself in the interim, it would probably be better to go for that.