[mesa-users] Hydrostatic equilibrium in MESA

[Bill Paxton]

On Apr 2, 2012, at 8:03 AM, Manos Chatzopoulos wrote:

> Hi Bill,
> I am trying to figure out a way to independently measure the level of hydrostatic equilibrium in particular MESA output models. 


Hi,

First thing to keep in mind is that r(k) and m(k) are radius and mass at the OUTER EDGE of cell k,
while P(k) and rho(k) are pressure and density at the CENTER of cell k.    That's important!
You cannot do something like (P(k-1) - P(k))/(r(k-1) - r(k)) <<< that's wrong. 
You need to use differences that are both at the center or both at the edge, but don't mix the two.

Next, go look at the code.  Open star/private/hydro_eqns.f and find the routine do1_pressure_eqn.
Note that the equation has been changed from d/dr to d/dm.  That's because we have accurate
values for dm at each cell that are constant during the solution (Lagrangian) whereas dr would
have to be calculated by subtracting the adjacent r's which can lead to nasty truncation errors.

So your base equation is dP/dm = -Gm/(4 pi r^4)
Since m and r are defined at cell boundaries, and P is at cell center, 
it is natural to set up the differences with respect to the boundary.
Remember that k=1 is at the surface, so outward dP = P(k-1) - P(k).
For dm we want the mass between cell centers which is avg of the cell dm's.
So the difference equation for hydrostatic equilibrium is
	(P(k-1) - P(k))/(dm(k-1) + dm(k))/2 = -G*m(k)/(4*pi*r(k)**4)

As a measure of how well this is satisfied, you might take
	abs(lhs - rhs) / (atol + rtol*max(abs(lhs),abs(rhs)))
where lhs and rhs are the sides of the equation and atol and rtol are
absolute and relative tolerances.

Finally, you need to add a few lines of code to your run_star_extras to output this to the profile.
Everything you need is in the structure you access using the s% idiom.
Edit how_many_extra_profile_columns to add another column.
Edit data_for_extra_profile_columns to provide the data for the profile.
For example, it might look like this:
         names(1) = 'measure_of_hydrostatic_balance'
         do k = 2, nz
            lhs = (s% P(k-1) - s% P(k))/(s% dm(k-1) + s% dm(k))/2
            rhs = -G*s% m(k)/(4*pi*s% r(k)**4)
            vals(k,1) = abs(lhs - rhs) / (atol + rtol*max(abs(lhs),abs(rhs)))
         end do
	 vals(1,1) = 0 

Add the name of your new column to your profile_columns.list and give it a try.

Cheers,
Bill





On Apr 2, 2012, at 8:03 AM, Manos Chatzopoulos wrote:

> Hi Bill,
> 
> I am trying to figure out a way to independently measure the level of hydrostatic equilibrium in particular MESA output models. The straightforward way to do so
> is to take the equation of hydrostatic equilibrium
> 
> dP/dr = rho*g
> 
> This equation is altered when rotation is considered with an extra term to account for the change of the shape of the equi-potential surfaces. That term is basically
> 1/2 * Omega^2 grad (r*sintheta)^2
> 
> In the case of a spherical MESA model of course sin(theta) = 1 and we know Omega so that term should be easy to evaluate, for a particular output model.
> 
> The question I have for you, is there a way that in the logXX.data output files we have columns for dP/dr? This would speed up my calculation.
> 
> I know that MESA as a stellar evolution code "forces" every model to be hydrostatic, however I do want to see departures from that when rotation is considered because
> I map some of my MESA output models to a hydro code to calculate explosions, pulses etc and those hydro codes do not correct for rotation, therefore once you map
> the model in them, its not going to be in hydrostatic equilibrium.
> 
> Thanks,
> 
> Manos