SUBROUTINE UniformGravity(REAL[] q, REAL[] dqdt)

Defined In:

contrib/astro/physics/i_i_ODEs.f90

Inputs:

REAL[] q. The 1D array that holds the field values for the cell.

REAL[] dqdt. A 1D INTENT(INOUT) array that is returned with UniformGravity()'s contribution to q's time-derivatives.

Description:

Applies time-derivatives for gravity source terms. This particular gravity subroutine orients gravitational forces towards the bottom edge of the problem domain. The mass of the central object is given by mCentral, read in from physics.data. This function is only called in 2D problems.

• dqdt = dqdt + f

Only f(3) & f(iE) are updated, with:

• f(3) = -q(1)*mCentral (~ scaled force density in the negative y-direction?)
• f(iE) = -q(3)*mCentral (~ scaled change in energy density wrt time ?)

In physics.data, mCentral is supposedly defined in units of mSolar (= 1.99e33 g).
In i_setprob.f90 (in subroutine CalcGlobalScale), for iGravity=2, mCentral is defined as:

• mCentral = mCentral/(velScale²/lscale)
  

It seems this puts mCentral as a scaled acceleration, but only works correctly if mCentral in physics.data is actually entered as an acceleration in cgs.

I'd like to change it to: gCentral = (mCentral * mSolar * G / rSolar²) / (velScale²/lscale)

With mCentral in units of mSolar, and rSolar defined in globaldeclarations.f90, this would give gCentral as a scaled gravitational acceleration.
Then I'd have here:

• f(3) = -q(1)*gCentral
• f(iE) = -q(3)*gCentral

Note that this subroutine treats dqdt as an accumulator—it will add its terms to whatever is already in dqdt without overwriting it.

Called In:

ODEs :: Derivatives2D()

Modules Used:

None.

Files Included:

None.