Simulation of gravity from density distribution, help?

[The victorious truther]

Dear forum, 

So I was curious about methods through which to simulate gravitational interactions and while I know there are numerous methods through which to numerically deal with multiple individual bodies I decided to find one to deal with general mass distributions. The model in question having ignore their interactions through non-gravitational forces or add in any accumulation of mass (formation of planetesimals) but still account for the change in the density distribution as well as the distribution of velocities. 

Firstly the primary equation in question that would be used is gauss law of gravitation to calculate the current gravitational potential distribution from the density distribution. Then i'd need an equation to encapsulate the conservation of mass to help with designating what the next distribution of matter would be. Finally, i'll need one more to determine the next velocity distribution and then repeat the whole process. I sort of guessed at this and found this set of PDE's that would suit my interest; gathered from continuum mechanics. 

1) \( \nabla^{2} \phi = 4 \pi G \rho \)

2) \( \frac{\partial \vec{u}}{\partial t} + \vec{u} \cdot \nabla \vec{u} = - \nabla \phi \)

3) \( \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{u}) =0 \)

The first is gauss law for gravity then the second is basically newtons second law but with a total derivative with the mass canceled out because of the weak equivalence principle (also because of mass conservation it takes on that form). The third and final equation is basically just a statement of mass conservation. 

I attempted to go on through and write these out full sale to then later discretize them or apply some other finite difference scheme to put into some programmable format using scipy, python, numpy, etc. These forms are given below, 

1) \( \frac{\partial^{2} \phi}{\partial x^{2}} + \frac{\partial^{2} \phi}{\partial y^{2}} + \frac{\partial^{2} \phi}{\partial z^{2}} = 4 \pi G \rho \)

2) 

\[ \begin{bmatrix} \frac{\partial u_{x}}{\partial t} \\ \frac{\partial u_{y}}{\partial t} \\ \frac{\partial u_{z}}{\partial t} \end{bmatrix} + \begin{bmatrix} u_{x} \frac{\partial u_{x}}{\partial x} + u_{y} \frac{\partial u_{y}}{\partial x} + u_{z} \frac{\partial u_{z}}{\partial x} \\  u_{x} \frac{\partial u_{x}}{\partial y} + u_{y} \frac{\partial u_{y}}{\partial y} + u_{z} \frac{\partial u_{z}}{\partial y} \\  u_{x} \frac{\partial u_{x}}{\partial z} + u_{y} \frac{\partial u_{y}}{\partial z}  + u_{z} \frac{\partial u_{z}}{\partial z}  \end{bmatrix} + \begin{bmatrix} \frac{\partial \phi}{\partial x} \\ \frac{\partial \phi}{\partial y} \\ \frac{\partial \phi}{\partial z}  \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} \]

3) \( \frac{\partial \rho}{\partial t} + u_{x} \frac{\partial \rho}{\partial x} + u_{y} \frac{\partial \rho}{\partial y} + u_{z} \frac{\partial \rho}{\partial z} + \rho \frac{\partial u_{x}}{\partial x} + \rho \frac{\partial u_{y}}{\partial y} + \rho \frac{\partial u_{z}}{\partial z}= 0  \)

Any assistance in terms of implementing these or stable numerical methods would be highly appreciated. 

Edited February 14, 2021 by The victorious truther
Latex was messed up.

[The victorious truther]

I'm most importantly wondering if this system of PDE's correctly describe exactly what i'm trying to model. Simulation of matter moving around while not interacting with itself aside from being gravitationally attractive. No mass accumulation or other non-gravitational forces present. No mass sinks or sources. 

[Markus Hanke]

8 hours ago, The victorious truther said:

I'm most importantly wondering if this system of PDE's correctly describe exactly what i'm trying to model. Simulation of matter moving around while not interacting with itself aside from being gravitationally attractive. No mass accumulation or other non-gravitational forces present. No mass sinks or sources. 

Is your density distribution a continuous function, or are you dealing with a collection of point masses?
At first glance your system of equations looks ok, though I’m not immediately sure whether (3) is actually needed at all, as (1)+(2) should already imply (3) - I haven’t explicitly checked though. As for numerical methods, that’s not my area of expertise, so I can’t offer any suggestions; I’d say it would be very difficult to implement that in code without constraining the form of the density and velocity functions in some way. 

How will the user of the software input the functions? Will there be an interface where the user types the functions symbolically (in which case you need to implement a suitable parser), or do they go directly into the source code, or what did you have in mind?

[The victorious truther]

14 hours ago, Markus Hanke said:

Is your density distribution a continuous function, or are you dealing with a collection of point masses?
At first glance your system of equations looks ok, though I’m not immediately sure whether (3) is actually needed at all, as (1)+(2) should already imply (3) - I haven’t explicitly checked though. As for numerical methods, that’s not my area of expertise, so I can’t offer any suggestions; I’d say it would be very difficult to implement that in code without constraining the form of the density and velocity functions in some way. 

How will the user of the software input the functions? Will there be an interface where the user types the functions symbolically (in which case you need to implement a suitable parser), or do they go directly into the source code, or what did you have in mind?

Well, the density distribution could be discontinuous and the idea here is that there are so many innumerable masses i'm thinking you could model it as an evolving non-interacting continuum of sorts (aside from gravitational interactions). Sort of like how we have equations to deal with individual H2O molecules and perhaps also a handful of them but as you increase the number of them then at some point you could fairly well approximate it by throwing out the assumption of a discretized fluid. Then model the fluid continuously with the Naviar-Stokes equations. I'm hoping the same can be done with a density distribution of sorts. 

Further, at the moment i'm just thinking of directly doing it in the source code. 

[Markus Hanke]

You should be able to do this using the equation system you have written down; density and flow are then vector functions of r, and you need to ensure that the proper boundary conditions (at the surface of the mass distribution) are imposed. This will leave you with a system of PDEs along with boundary and initial conditions. You then need to find a suitable numerical algorithm to generate solutions (not my expertise, so can’t help with that), and some suitable way to plot them graphically.

All in all this is quite a formidable task, both mathematically and in terms of coding - so best of luck  

[sethoflagos]

On 2/15/2021 at 12:41 AM, The victorious truther said:

I'm most importantly wondering if this system of PDE's correctly describe exactly what i'm trying to model. Simulation of matter moving around while not interacting with itself aside from being gravitationally attractive. No mass accumulation or other non-gravitational forces present. No mass sinks or sources. 

 In principle, your equations look okay given your stated constraints, but I would strongly recommend converting to polar coordinates to get as much help as possible from its symmetries. 

Given that your gravitational field terms are directly analogous to the pressure gradient terms in Navier-Stokes applications that I'm somewhat familiar with, you may find that at least for some simple starting conditions, the method of characteristics may help convert your PDEs into ODEs which would then be amenable to numerical integration. 

However, you're using Gauss' Law where I would normally be inserting a (simpler) equation of state, and that may complicate matters significantly. 

[sethoflagos]

On further consideration, by omitting the pressure terms from the Navier-Stokes, haven't you lost control of the conservation of energy? For a rotating system, I can visualise material falling into an equatorial disk from both sides, but there appears to be no mechanism in your system for it to cross the boundary. The momenta will just cancel (conserving momentum) but then so will the corresponding kinetic energy! Restoring the pressure term (which in context, is an expression of the system internal energy) will keep the 1st Law books straight and avoid some embarrassing infinite densities.

[The victorious truther]

I got the second equation wrong when I expanded it out.  It should look like that. 🙁