Physics Simulation using Coulomb

[danielrmeyer]

I wanted to see what a simple physics simulation would like when physical dimensions were added to the code. To do so I adopted an example from the book Computer Simulation Methods and the corresponding Java library that goes with the book (osp). The code for managing the animation and plotting is basically the same as in the book; however, the code that represents the physical model and the RK4 solver were re-written to be aware of physical dimensions.

The experience was really interesting. To take full advantage of Coulomb it seemed as though I had to make the ODE solver code aware of Coulomb types. Currently the ODE solver is forced to work with a State type that has fields of meters, meters/seconds, and seconds just like in a physics text book. From the perspective of Physics this seems great since we are now working with the same units we work with in the field. From a software engineering perspective I am wondering how we would make a numerical library that is aware of dimensions like this but also generic enough to be useful. Would it be possible to have an ODE solver interface that understands any type of the form [Abstract scalar quantity / Second, Abstract Scalar quantity, Second] so that we can send in states of [Meter, Meter / Second, Second] or [Kilogram, Kilogram / Second, Second] and have that work out?

There was no maven artifact for the osp library, so I built a fat jar in eclipse and included it in the following repository:

https://github.com/danielrmeyer/TypedPhysicsSimulations

Which is built with sbt. For a quick taste here is the code implementing the physical model:

class SHOModel extends Circle with ODE {
  val k = 1.withUnit[Newton / Meter]; // spring constant
  val b = 0.2.withUnit[Newton * Second / Meter] // damping constant

  var state = State(0.withUnit[Meter], 0.withUnit[Meter / Second], 0.withUnit[Second])

  val ode_solver = new RK4(this)

  def getTime(): Quantity[Double, Second] = state.t

  // ODE interface
  def getState() = state

  // ODE interface
  def getRate(state: State): State =
    val x = state.x
    val v = state.v
    //Here we see the compiler do dimensional analysis for us
    val force: Quantity[Double, Newton] = -k * x - b * v
    val rate = force / 1.withUnit[Kilogram]
    State(v * 1.0.withUnit[Second], rate * 1.0.withUnit[Second], 1.withUnit[Second])

The ODE interface (see also RK4.scala) now needs to understand the dimensions of the state:

case class State(x: Quantity[Double, Meter], v: Quantity[Double, Meter / Second], t: Quantity[Double, Second])

trait ODE {

  def getState(): State

  def getRate(state:State): State

}

The files of interest are here:
https://github.com/danielrmeyer/TypedPhysicsSimulations/blob/main/src/main/scala/SHOApp.scala
https://github.com/danielrmeyer/TypedPhysicsSimulations/blob/main/src/main/scala/numerics/ODE.scala
https://github.com/danielrmeyer/TypedPhysicsSimulations/blob/main/src/main/scala/numerics/RK4.scala

In any case I would love feed back on this code and how I am using the Coulomb library. I felt a rush of excitement as the compiler worked out the physical dimensions for me.

I also found that there were a bunch of implicit assumptions in the original RK4 code (RK4.java from the osp library) that I was forced to make explicit by adding in physical types. I am planning to create some more examples like this. I am very interested in the question of whether the underlying numerical libraries should be aware of dimensions like I did here in RK4.scala and ODE.scala

When running the simulation you should see something like the following:

[erikerlandson]

@danielrmeyer that is very cool, thank you for sharing it! I will have to read through this more carefully when I have a minute.

One of the things I hoped people could use coulomb for is doing unit aware numeric analysis. I'm happy to see other people such as yourself looking at this idea!

[erikerlandson]

I also found that there were a bunch of implicit assumptions in the original RK4 code (RK4.java from the osp library) that I was forced to make explicit by adding in physical types.

Overall this seems like a very good thing, and in fact "making implicit assumptions explicit" is one of the reasons for strong type systems in general. I tried to expand on this connection more in this blog:
https://erikerlandson.github.io/blog/2020/04/26/your-data-type-is-a-unit/

I'm unsure if an integration with RK4 is a good fit for the coulomb repo itself, but I'd be happy to help advise on a coulomb/RK4 integration if I can.

[erikerlandson]

In any case I would love feed back on this code and how I am using the Coulomb library

@danielrmeyer IMO you are using coulomb exactly the ways I had in mind. I'm curious how your experience was using it. One of my goals was to make it ergonomic and nice to use, since adding unit awareness to code does increase complexity some.

[erikerlandson]

@danielrmeyer when you mention making "underlying numerical libraries" unit aware, which libraries do you have in mind?

[danielrmeyer]

@erikerlandson Thanks for looking at this! For "underlying numerical libraries" I am thinking of how we can interact with things like ODE solvers and not lose the unit consistency. Most of the ODE solvers want arrays of doubles and so we have to convert input values into doubles and then convert them back to meters or seconds--which is unsatisfying. I am not an expert in this area by any means so it is possible I am not thinking about this right; however, after pondering it for a couple weeks I am really leaning towards the idea that differential equations do have units. The very concept of the derivative seems to imply that there is at least a time unit or equivalent and therefore a rigorous mathematically correct ODE solver ought to understand units. Also, in this example I felt I had to re-write the ODE solver for the example to come out compelling.

I really enjoyed using Coulomb and I would love to help develop it or crank out some more examples like above. Currently I do software development for a living; however, I studied Physics in college an I am a thesis short of my graduate degree (which I dropped to get into the software industry). I have a goal of going back and doing Physics at some point in the future and my focus would be to improve how scientific computation is done and unit awareness is forefront on my mind. In Physics they teach us to "check your work" my making sure the units come out right. I remember getting points taken on exams for not checking the units. Unfortunately, when we start trying to use the computer to solve physics problems we have to drop all of that and just work with floats and doubles and the studious might keep units in the comments. In any case I had been thinking scala 3 would be a good fit for scientific computation and then came across the unit libraries in my research--I landed on Coulomb because I saw you were taking advantage of the opaque type, which seems like a really important feature for this kind of work. I would like to do some bench marking at some point to see if we are really getting zero runtime overhead.

I am trying to figure out my next steps here. I am thinking of trying to start converting that entire osp (open source physics) library over to scala 3 with Coulomb integrated. The osp library is a java library used for teaching and I can imagine an argument for releasing a new copy of that book if the scala 3 + Coulomb units prove compelling enough (I am not affiliated with that book in anyway). It sure looks compelling to me.

Another area of interest of mine is symbolic computation and I am planning to see how far I can get implementing a variable type akin to what is in the sympy python library, but I would like to attach a unit to the symbol. Rather than var x = Symbol('x') I'd like to say var x = Symbol('x', Meter) or something like that and have the resulting symbolic expressions type checked on units.

I have admittedly taken a detour to check out the uom library in Rust recently. I would rather stick with scala 3 as I think the syntax gets us closer to actually doing physics on the computer. I've never been that interested in Julia because I think the future is in static typing--I just don't think it is compelling enough to get me off Python for things where a dynamic language is appropriate. I also looked into F#; however, the fact that scala 3 is powerful enough to encode this stuff in a library won me over quickly.

In any case I think Coulomb is a big step towards actually using the computer to do physics the way it really should be done. I am super excited about this library. Thank you for putting in the time to build it!

[erikerlandson]

@danielrmeyer if you are interested in the compilation properties of opaque type, and how expressions are (or aren't) optimized in coulomb, you might find this discussion interesting:
#537 (comment)
Or also this:
typelevel/spire#1310

In principle opaque type seems like the correct abstraction but currently I believe AnyVal has superior compilation optimizations. @benhutchison also is advocating me to stick it out with opaque and I'm willing to let the compiler catch up for now.

[erikerlandson]

I am thinking of trying to start converting that entire osp (open source physics) library over to scala 3 with Coulomb integrated.

This is the kind of community integration that I think would be very valuable, and a great testbed for the ergonomics of coulomb - I spend most of my spare time fooling around with the coulomb core itself, trying to figure out how to make it better, and although I loved freshman physics and chem, it's been many years :)

I agree that unit analysis is a great fit for integral and differential calculus, since the relation between the units of an expression, its integrals, and its derivatives, are all very structured and deterministic.

[erikerlandson]

@danielrmeyer You might also find this interesting, where I work out a unit analysis of the vector algebra for linear regression:
https://erikerlandson.github.io/blog/2020/05/06/unit-analysis-for-linear-regression/

Or these:
https://erikerlandson.github.io/blog/2022/08/12/a-unit-analysis-of-matrix-inversions-and-svd/