Tracking the mass that falls onto the BE sphere

I would like to calculate the mass that falls onto a shell of dr=0.1RBe at r=Rbe from the ambient over the course of a simulation. There may be a few ways to do this, but easiest may be what came to mind today: find M(t) of the sphere with r=1.1Rbe, using Visit's query option for total mass in box.

The method

First, I had to determine the right sum to use in Visit. I see that the 'query' option had 2 sums: weighted and not. Additionally, I found that the sum query is for the current plot. If that plot is a slice of the full simulation, the query will return only a sum of variables from that slice's mesh. Using my knowledge of the BE sphere's total mass from my simulation output —- that is, 151 solar mass, I calculated the mass of the ambient quite easily, given the ambient is of uniform density:

amb = (L^3 - {4\over3} * r^3 * \pi) \rho

Putting this into computational units, I discovered that weighted variable sum was the correct query to get the total rho/cell (and hence mass) in the grid.

The plan was to calculate M(t, 1.1Rbe) using the following formula:

(t, 1.1Rbe) = Mtotal(t)-Mambient(t)

I knew I could query-over-time the Mtotal and perhaps make a clip of the box that would isolate only the ambient, and then query-over-time the Mambient. However, I checked first whether I could just treat Mambient as constant. A movie showing rho(t) of just the ambient indicated I could:

Considering the Mtot=Mtot(t) as given by Visit's query:

I was curious where the change in mass was occurring. Since the density of the ambient, and its volume is staying ~ constant, must be inside of the sphere and surrounding vicinity… Here are some plots of the density inside of the sphere at different times:

So, the density is changing, but is the total mass? Visit can integrate these curves, and indeed the area under the curves were different, indicating a change in the total mass of the sphere.

Thus, it seemed justifiable to make the assumption: 1) mass of ambient stays constant, 2) calculate M(t) of the sphere with r=1.1Rbe using the Mtot query of visit. Also, since the mass of the sphere is changing, what does this mean about the boundaries supplying/removing mass from the box?

Using this data, I found the following curve for M(t) of this sphere:

I see that Mbe crit is below the minima of M(t). This is no doubt due to some SMALL round-off error somewhere in computing the total masses in the grid. If for instance, I used simply the formula above to find the mass of the BE sphere, it does not give me exactly what my code gives. The Visit calculation gives me ~0.78 when I subtract the mass of the ambient from the mass of the sphere, whereas my code says it is ~0.72 — in line with the Crit Mbe as expected. This seems to indicate that these curves may be closer to each other than the are here. It would be interesting to see where these curves intersected. At least for now, it is safe to take away from this plot that M(t) does not fluctuate wildly, and never goes more than ~25% over the crit Mbe.

Now, it should also be noted that the crit Mbe is computed by the formula:

be_{crit} = 1.18 * {Cs^4 \over Po^{1/2} * G^{3/2}}

where Cs = isothermal sound speed, Po = external pressure, and G is gravitational constant. This Mbe is a measure of the maximum mass a BE sphere can hold up against its own gravity. For the values used to compute the crit Mbe, I used Cs for the BE sphere (r=Rbe), whereas the Cs would actually increase 10-fold outside of the sphere, since the boundary at Rbe is held at pressure equilibrium.