A physicist view on gRPC — Simulating the Sun-Earth-Moon three-body system
Leveraging distributed computing to simulate celestial bodies
All codes are available on GitHub.
It is very frequent, in trading system to define proprietary protocols and implement them using proprietary RPC frameworks.
Recently, I was told about gRPC. And upon hearing about gRPC,
(i) I was surprised that the trading systems I worked on until not long ago, not so many were leveraging the technology;
(ii) I just wanted to implement something to try it out.
Why? What is gRPC? If you ask my new best coding partner (chatGPT), he’ll say this
gRPC is a high-performance, open-source framework for RPC (Remote Procedure Call) developed by Google. It enables client and server applications to communicate transparently and efficiently across different environments and languages.
In a nutshell, gRPC allows you to write applications across multiple processes as it provides the communication protocols as well as the server/client codes.
There is often a lot of boilerplate code that you can avoid with the right abstractions and that is what gRPC is: a great abtraction of RPC’s.
For those of you who have read some of my former publications, you may be aware that Rust has taken a lot of my interest lately. There are crates to use gRPC with Rust. I will probably be using these crates in the near future but as they are non official, and to get introduced to gRPC, I thought I would take the opportunity to gain some programming experience in another language that is fully supported by gRPC dev’ teams.
Since Golang is available with an official support in gRPC and I have always been found of the language, I threw myself in without looking back.
As for the project, and as the name of this blog states, it had to be about some form of mathematical modeling. Recently, I had been publishing a lot about Quant Finance, so my heart felt like doing something else.
Thus, a little out of nowhere, I figured I would simulate the whole Solar System: one gRPC client per celestial body. I was thrilled about that idea.
However, while working on the project, I resized it as I realized that I had to carefully design server and clients using gRPC, but also the Ordinary Differential Equation solver. The latter ended up taking more time than initially planned. So I decided to focus on the classical three body system made of the Sun, the Earth and the Moon. I can always come back to the full Solar System later on.
(Hint… Hint…)
So we’ll build a gRPC system with a central server that shares informations amongst clients: each of these representing either the Sun, the Earth or the Moon.
This project was very pleasant as it threw me back in space and the wonderful mechanical physics that goes with it. I should say though, this post is mainly about gRPC rather than celestial mechanics. As you shall see, all I have done here is coherent from a physcial stand point, but the main take away is that you can solve the Sun-Earth-Moon three body system using gRPC. And I find that very cool…
I even think that it could make a lot of sense to systematically use gRPC for complex physical calculations. gRPC offers a very convenient way to define whatever celestial system you like!
First we’ll go through the maths: you don’t really need a proper experience in mechanics to understand the rough idea — if you disagree with me on that point, I am sorry, and please tell me how I failed to make that material accessible.
Second, I’ll get into the details of the implementation. Although I used chatGPT (obviously) to learn about gRPC and debug my Golang doodles, the codes are really mine as I don’t think one can properly learn anything without doing things themselves.
Third, we’ll have a look at the results. Initially, I was having really crazy results and thought something was wrong with my programs… Turns out it was the initial conditions! They are extremely important. It is now obvious why, but initially, I really thought I could just throw some rough numbers, but it turns out that I actually had to be very careful with these initial parameters!
行くぞ！
Background Maths
Let’s go back to undergrad physics and write down the equations of motion for the three body system made of the Sun, the Earth and the Moon.
Center of Mass is the center of the simulation
On the figure below, a non-scaled diagram of the Sun (point S), the Earth (point E) and the Moon (point M).
Each of these points has a mass and classical mechanics defines the center of mass G of these points with the equation
assuming that O is a fixed point in space. Now if we assume that i (either S, E, M) exerts a force on j (either S, E, M, but not i) and we write it
we could also draw a more detailed diagram with the interactions between the bodies as shown on the figure below.
The 2nd law of Newton applied to each of the bodies gives
In English, the 2nd law of Newton states that the rate of change of velocity of each of the bodies (Sun, Earth or Moon) is equal to the forces applied on each of them.
And thus, by taking the second order derivatives in the definition of the center of mass, we find that the second order derivative of OG is
The 3rd law of Newton now states that necessarily for every couples of i and j such that i is not equal to j
In English, the force exerted by the Moon on the Earth is exactly the opposite of the force exerted by the Earth on the Moon. And his statement remains true for any pair in the three body system.
The result of all these statements is that the right hand side of the equation on the second order derivative of OG is equal to 0
In English, the center of Gravity has uniform motion in a right line: it never accelerates.
And that is why we can set the center of mass as the center of the frame used in the calculations — I am voluntarily not mentionning Galilean Frames but for those of you who are familiar with these concepts, that is what the center of mass G defines.
Equations of Motion
Now that we have our center of frame, let us delve deeper into the geometry and state the equations we will use to simulate the motion of each body.
Let’s start with the position vectors of each celestial body. We write them:
which can be shown on the diagram
where i is any of S, E or M. If we introduce the vector
that represents the position of j relative to i, according to Newton’s law of universal gravitation, the gravitational pull of j on i is
which implies that the 2nd of Newton on body i may be stated as
And introducing a cartesian frame that is centered in G, we find that the equations of motions are
where xᵢ and yᵢ are the coordinates of i.
Dimensionless Equations of Motion
A common practice in physics computations, is to render equations dimensionless to avoid over- or underflows. Physics equations can have quantities with very great values (like in this case) or very little values (like in molecular dynamics or quantum physics).
Here G is extremly small as G = 6,67430e−11 and distances are extremely big as the Sun-Earth distance is about 1.5e+08 kilometers… Some of you may notice that there seems to be a compensation between G and the distances. So yes maybe it is not necessary to rely on dimensionless equations, but it has the merit to remove that worry, so let’s do it anyway.
To get to the dimensionless equation, one has to introduce scaling parameters and dimensionless variables. In our Sun-Earth-Moon three body system, I chose to introduce the scaling parameters
where D is the minimal distance between the Sun and the Earth and T is of the order of magnitude of the orbital period of the Earth around the Sun — there should be a 2π for the full formula.
I then naturally chose the dimensionless variables as
If you use the scaled masses
you eventually find that the dimensionless equations are (I dropped the ~’s for clarity) in cartesian coordinates
or more formally (dropping this time the i’s for clarity)
Before moving on to the implementation section, it will be convenient to state the above equation in an even more reduced form using the vector
Indeed, with the right definition for the function f, the equations of motion become
(assuming that all other bodies’ positions are paremeters of f).
Implementing the system
Let’s get real now, shall we? 🤓
What we want — An implementation of a central server that shares across the clients the computations done by all.
In English, the computational system will be made of clients running the mechanical computations. One client per celestial body. Each client will then send the update to the central server that will in turn broadcast the new information to every client.
First I will talk about the use of gRPC to generate the Golang files necessary for the implementations, second I will give an overview of the implementation of the server, and third, an overview of the client.
The gRPC protocol
Defined in a .proto file, the protocol states what the transmitted data is like, as well as the services available in the server.
syntax = "proto3";
package taiyoukei;
option go_package = "./proto";
message CelestialBody {
uint64 sequence = 1; // Seqwuence in stream
string name = 2; // Name of the celestial body
double mass = 3; // Mass of the celestial body
double x = 4; // x-position
double y = 5; // y-position
double vx = 6; // x-speed
double vy = 7; // x-speed
}
message CelestialBodiesPositionRequest {}
service CelestialService {
rpc CelestialUpdate(stream CelestialBody) returns (stream Data) {}
rpc CelestialBodiesPositions(CelestialBodiesPositionRequest) returns (stream Data) {}
}
message Data {
bool success = 1;
map<string, CelestialBody> content = 2;
}There is only one message being sent around: it contains the kinematic information of the body as well as its mass, its name and a sequence number. The sequence number is there to ensure that the flow of messages remains consistent.
As you can see, there are two services defined with the keyword rpc:
CelestialUpdate, to be used between the server and the computing clientsCelestialBodiesPositions, to be used between the server and an outside watcher used for the analysis of the results.
We will concentrate on the first one as it is the most interesting feature of this server. The second one was straight forward and does not bear any interesting aspects.
To generate the gRPC code I used the command
protoc --go_out=. --go-grpc_out=. --proto_path=proto/ proto/celestial.protowhich produced the files
% tree proto
proto
├── celestial.pb.go # GENERATED
├── celestial.proto # Source file - NOT GENERATED
└── celestial_grpc.pb.go # GENERATEDcelestial.pb.go and celestial_grpc.pb.go contain the boilerplate code for the server and clients infrastructures. The specialized logic for the application we are writting is, however, completely absent.
And that is what I present below.
The Server
Here are the requirements for the server
(i) — it waits for all updates from the clients
(ii) — once all update data is received, it broadcasts it to all clients
When the clients receive the broadcast, they have the positions of all other bodies and they can, in turn, update their position using the equations of motion and resend an update to the server, that resends a broadcast once all updates are received again, etc.
The mental image I have of the data flow between the server and each client is similar to this
The general structure of the CelestialUpdate program is given below.
/*
SERVICE -- bidirectional streaming with clients: reception of updates
and broadcast of these updates to all clients.
*/
func (s *server) CelestialUpdate(stream pb.CelestialService_CelestialUpdateServer) error {
/* Connection object definition and storage in map */
stop := make(chan bool) // -> Channel to ensure smooth stop
if n_connections == ONE {
go func() {
/*
Handles the emission of broadcasts
- only started when there is at least one connected client
- killed when the number of clients drops to zero
*/
err := s.sendData(&data)
if err != nil {
stop <- true // Stops the loop forloop in the parent routine
break
}
}()
}
log.Printf("Initiating data reception for id %v", id)
forloop:
for {
select {
case <-stop:
break forloop
default:
/*
Handles the reception of client updates
*/
}
}
return nil
}CelestialUpdate is a method of an interface generated from the .proto file. And in the generic code shown above, I basically implement the interface for a server struct that I defined myself.
The most important feature of this method is the go-routine started under the condition if n_connection == ONE. This go-routine will be in charge with the broadcast of the received data. Without this simple design, I was literally banging my head trying to find an easy way to ensure the synchronicity between client threads and the broadcast…
Having a dedicated go-routine made everything a lot easier, even if I have to admit that, in terms of design, it is not ideal.
Now that I am writing these lines, I wonder if it would not have been better to have both client- and server-streams instead a single bi-directional stream*…
* I’ll let you tell me in the comments what you think would have been a better solution.
The Client
As for the client, there is nothing particularly difficult with regards to the gRPC communications. The interesting part was in the solver that I implemented to simulate the orbits.
Numerical Solver — Runge-Kutta of order 4 (RK4), a go-to solver stable enough to have a coherent system and easy enough to not loose your hair with endless stability algebra.
Using the previously introduced generic equation
the RK4 algorithm states that for a discrete step Δt the sequence of discrete vectors σₙ is given by the following recursion
Now here is the fantastic part: using a strategy design pattern in Go, I was able to very easily implement the RK4. I did not even have to debug the algorithm, from the get-go, it was correct — which is rare! Not that I am particularly amazing at writing all of this, nor that ChatGPT gave me the answer (I made a point of honor to implement things myself), but by staying as close as possible to the mathematical statements, chances of mistake are greatly reduced. So,
(i) — I introduced an implementation of vector σ, hamiltonVector below, with all the algebraic necessary tools (addition and scalar multiplications)
(ii) — I introduced a rateFunction type in the signature of my solver, providing a strategy like design pattern for the rate function
(iii) — I introduced a Solver interface providing another strategy-like design pattern for the solver this time.
/*
STRATEGY PATTERN 1 TO SOLVE THE ODE ON RATE FUNCTION
*/
type rateFunction func(hamiltonVector) hamiltonVector
/*
STRATEGY PATTERN 2 TO SOLVE THE ODE ON SOLVER
*/
type Solver interface {
step(sigma hamiltonVector, f rateFunction, dt float64) hamiltonVector
}
type RungeKutta4Solver struct {
}
/*
Implementation of RK4
*/
func (rks *RungeKutta4Solver) step(sigma hamiltonVector, f rateFunction, dt float64) hamiltonVector {
F1 := f(sigma)
F2 := f(sigma.Add(F1.scalarMultiply(dt / 2.0)))
F3 := f(sigma.Add(F2.scalarMultiply(dt / 2.0)))
F4 := f(sigma.Add(F3.scalarMultiply(dt)))
sum_F := F1.Add(F2.scalarMultiply(2.0))
sum_F = sum_F.Add(F3.scalarMultiply(2.0))
sum_F = sum_F.Add(F4)
return sigma.Add(sum_F.scalarMultiply(dt / 6.0))
}
type celestialRateFunctionHandler struct {
bodies []hamiltonVector
masses []float64
}
/*
Equations of motion -- specific implementation of the rate function
*/
func (fh *celestialRateFunctionHandler) evaluate(sigmai hamiltonVector) hamiltonVector {
rate := hamiltonVector{
x: sigmai.vx,
y: sigmai.vy,
vx: 0.,
vy: 0.,
}
for j, sigmaj := range fh.bodies {
muj := fh.masses[j]
d_ij := math.Sqrt((sigmaj.x-sigmai.x)*(sigmaj.x-sigmai.x) + (sigmaj.y-sigmai.y)*(sigmaj.y-sigmai.y))
Fx := muj * (sigmaj.x - sigmai.x) / (d_ij * d_ij * d_ij)
Fy := muj * (sigmaj.y - sigmai.y) / (d_ij * d_ij * d_ij)
rate.vx += Fx
rate.vy += Fy
log.Printf("mu = %v d = %v, F = %v", muj, d_ij, math.Sqrt(Fx*Fx+Fy*Fy))
}
return rate
}Going further — You ought to know that RK4 eventually diverges. State of the art orbital computations are now performed with the Gauss-Jackson algorithm.
But…
- It is a lot harder and it would have to be an blogpost by itself.
- RK4 is perfect for a first version of the computing system.
- With the strategy pattern, it will be a breeze to modify the algorithm used to propagate the orbits in the clients.
Results
Now, after all this hardworking, the results. At first, here are the initial conditions I gave to the the system.
- I assumed that the Sun, Earth and Moon are all in the same plane. This assumption was there from the beginning. It is actually wrong because there is a 5 degree incline between the Earth Orbital Plane and the Moon orbital plane. But we (secretly) neglected that.
- There has to be a starting point when all three are aligned — e.g. during an eclipse. And I chose the alignment SME. And this allowed me to set all the y values to 0.
- Because of the masses, the Sun is about a million times more massive than the Earth which is in turn about a hundred times more massive than the Moon. So it seemed fair that G would be the same as S, leading to the values
- As for the velocities, the masses were again such that in first approximations
These initial conditions were plain wrong. Way too approximative! I would always get a divergence of the Moon as it would get to close to Earth. Which would happen everytime because of their initial position and speed. Eventually the Moon would even have its own orbit around the Sun. 😇
So I took my pen, stopped being lazy and rederived appropriate initial conditions. Here is what I got (in adimensional form).
- For the initial positions:
- For the initial velocities
As the Earth and the Moon both orbit counter-clockwise from a North Pole view point. The Moon’s velocity above is set to the orbital velocity of the Earth minus the orbital velocity of the Moon around the Earth. This is because the initial state is in SME configuration.
Had it been in SEM configuration (the Sun, then the Earth, and finally the Moon), the velocities would have been added and not subtracted!
Back when I was preparing for the Ecole Polytechnique, getting Numerical Values of Theoretical Formulae was always a nightmare. I would always mess up something. But thank god, ChatGPT is here to the rescue.
Here is a screenshot of my request and his response.
Using these values as input, I enter these commands (in different tabs of my terminal)
# To start the server
go run server/run.go
# To start the capture client -- not detailed in the text but use it
# to generate the orbit graphs
go run capture/main.go
# To start the Sun
go run client/run.go -name S -mass 1 -x -0.0000030393 -y 0 -vx 0 -vy 0
# To start the Earth
go run client/run.go -name E -mass 0.0000030025 -x 0.999997 -y 0 -vx 0 -vy 1
# To start the Moon
go run client/run.go -name M -mass 0.000000036938 -x 0.997427 -y 0 -vx 0 -vy 0.965816The computations will only start once all the expected clients are connected — here 3 clients are expected.
Now let’s look ate the plots generated by capture/main.go to assess whether or not the simulation makes sense.
Scroll down a little more….
ET VOILÀ!
This output looked extremly promising as you could see that the Earth E orbits around the Sun S and so does the Moon M. In my first trials the Moon was just wrongly being thrown into a different orbit around the Sun than the Earth’s one.
Also, note how the Moon orbit (in green) seems to be crossing back and forth the Earth Orbit (in red). This is a possible indication that the Moon is, as expected from observations, orbiting around the Earth as well.
So let’s check this by centering the graphs on the Earth.
Scroll down a little more….
(Yes I am that annoying)
やったあああああああああああああああああ！
