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.

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.

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}}

This describes blah blah, whereas my sphere is _blah Blah.

I guess what can be taken home safely from this curve is that M(t) does not fluctuate wildly, and never goes more than 10% over the crit Mbe.