Dynamic system of units

[jedbrown]

It's common in computational physics (including fluid dynamics) for a solver to transform all quantities in terms of reference scales in order to improve conditioning of matrices that will arise. This can be accomplished by choosing representative base units and parsing all user input into that system, then further operations will be zero-cost. Fuzzing these scales is also a useful technique when dimensionally-unsafe code is present (such as calls to linear algebra libraries; consider the 5x5 Jacobian matrix that would arise from differentiating the From implementations).

Do you have thoughts about how to handle this? I don't know how to get serde_yaml::from_str to parse into a system of units with dynamic scales, but perhaps we could parse into SI and have conversion methods.

[Tehforsch]

I am not very familiar with this, so I have some questions before I comment on the implementation side of it.

If I understood you correctly, the point is to choose different base representations (one for each dimension), right? So, for example, instead of everything being converted to meters, as is currently the case in the SI definitions, it would be converted to kilometers. This can currently be achieved "statically", i.e. by defining the system as such:

unit_system!(
    Quantity,
    Dimension,
    [
        def Length = { length: 1 },
        unit (meters, "m") = 1e-3 * kilometers,
        unit (kilometers, "km") = Length,
    ]
);

(As a side note, in the upcoming version I want to improve the unit_system macro, so that this will be made more convenient.)

But if I understand you correctly, you'd like to define the scale during runtime i.e. by computing some relevant length scale of the system and then defining the unit system accordingly?

Fuzzing these scales is also a useful technique when dimensionally-unsafe code is present (such as calls to linear algebra libraries; consider the 5x5 Jacobian matrix that would arise from differentiating the From implementations).

What do you mean by fuzzing these scales? Changing the scales and seeing if the output is the same as a method of testing the interface to the dimensionally-unsafe code?

Do you have thoughts about how to handle this? I don't know how to get serde_yaml::from_str to parse into a system of units with dynamic scales, but perhaps we could parse into SI and have conversion methods.

That sounds reasonable to me. I assume that there would be some kind of conversion struct that specifies the factor for each dimension anyways, so that could be used to convert between what I assume will be two different Quantity types.

[jedbrown]

Yeah, we'd use scales of the input (e.g., length, flow velocity) to set internal units and yes, that's what the fuzzing would be for. We use an ad-hoc fuzzing in our current computational mechanics solvers (which are mostly written in C without dimensional analysis libs) to audit dimensional consistency, but it's a blunt tool because it doesn't tell you where the problem is.

I assume we end up with dynamic::Quantity that can only be created by converting from an si::Quantity of the same dimensions while using a specified (non-const) reference scale. What I don't know how to ensure is to distinguish from another dynamic::Quantity created with respect to a different reference scale. Perhaps it can be done with "branding" using phantom lifetimes or some other ZST trick, but I don't know if such constructs can work with adt/const generic features. Pragmatically, I'd be okay even if we can't statically guarantee this -- in the applications I envision, conflating quantities from multiple reference scales would be a very unlikely bug.

[Tehforsch]

I assume we end up with dynamic::Quantity that can only be created by converting from an si::Quantity of the same dimensions while using a specified (non-const) reference scale. What I don't know how to ensure is to distinguish from another dynamic::Quantity created with respect to a different reference scale. Perhaps it can be done with "branding" using phantom lifetimes or some other ZST trick, but I don't know if such constructs can work with adt/const generic features. Pragmatically, I'd be okay even if we can't statically guarantee this -- in the applications I envision, conflating quantities from multiple reference scales would be a very unlikely bug.

I currently can't think of a way to do this statically, but dynamically one ugly implementation strategy that I can think of is introducing some kind of global AtomicUsize which gets incremented every time a new reference scale is created and then gets stored inside any dynamic::Quantity created using that scale. There would be some additional work required to make sure that creating the "same" reference scale twice (which could very easily happen in user code) will still create compatible quantities. Its very messy, but its all I could think of for now.

[Tehforsch]

Your comment in #56 got me thinking about this again:

The gauge idea can perhaps be better expressed as a static case of https://github.com/Tehforsch/diman/issues/35, but I suspect something along those lines will be useful when working with physical sensors.

For the static version of this, I think adding a sub-expression to the unit_system macro could be a very clean solution - something like this modulo syntax

unit_system!(
    quantity_type Quantity1;
    #[base(Length)]
    unit meters;
    ...
    system(
        quantity_type Quantity2;
        #[base(Length)]
        unit kilometers;
    )
);

This would "only" have to add a new quantity type and add the numeric traits and methods on that quantity and probably add a single conversion method such as let q: Quantity2 = quantity1.as_quantity2(). Maybe I am missing some of the gory details here, but I think that should work without too many problems.

If I understand the requirements for dynamic systems correctly, we'd compute some quantities at runtime:

    let (length: Length<f32>, time: Time<f32>) = compute();

For the sake of simplicity, lets say that these are already the base dimensions of the original unit system. Ideally, we'd like to have a quantity type that behaves exactly like the original quantity, except that it's base representation is different. Now, in reality I don't think this will quite work, but I think we could at the very least do something like

struct DynamicQuantity<const D: usize, S>(S);

fn main() {
    let (reference_length, reference_time) = compute_reference_scales();
    let constructor = Constructor { length: reference_length, time: reference_time };
    let length: DynamicQuantity<...> = constructor.new(meter, 100.0);
    let time: DynamicQuantity<...> = constructor.new_from_quantity(10.0 * second);
}

Under the hood, this would be enabled by code which is automatically generated by the unit_system macro that generates the Constructor type as well as the DynamicQuantity type and reads (very, very approximately):

struct Constructor {
    length: Length,
    time: Time,
}

impl Constructor {
    fn new<const D: Dimension, const M: Magnitude, S>(u: Unit<D, M>, val: S) -> DynamicQuantity<S, D>
    where
        S: Mul<Magnitude, Output = S>,
    {
        DynamicQuantity(S * M * self.conversion_factor::<D>())
    }

    fn new_from_quantity<const D: Dimension, S>(q: Quantity<S, D>) -> DynamicQuantity<S, D> {
        DynamicQuantity(q.value_unchecked() * self.conversion_factor(D))
    }

    fn conversion_factor<const D: Dimension>(&self) {
        // or 1 over that, who knows
        self.length.powi(D.length) * self.time.powi(D.time)
    }
}

In this implementation, the DynamicQuantity is basically exactly like a normal Quantity, except that it doesn't have any constructors of its own (or, if we want the API the other way around, all of its constructors require the Constructor/Scale type as argument)

Of course, this does not guarantee any unit safety yet. One can very easily mix two differently generated dynamic quantities with this implementation.

If one needs multiple reference scales at the same time, an alternative could be to extend the Quantity type to contain a third generic argument: struct Quantity<const D: Dimension, S, ZST>(S, core::marker::PhantomData<ZST>); and to implement all traits only on Quantity with matching ZST. (On a side note, I am not entirely sure we want to rely on the the compiler always optimizing PhantomData away for a struct like this that absolutely requires it) This would provide having multiple (that is, n, where n is defined at compile time) reference scales in a unit-safe way, but does not really help with creating thousands of them and preventing mixing those.