Re: comoving implementation

From: Volker Springel <volker_at_MPA-Garching.MPG.DE>
Date: Sat, 03 Feb 2007 00:45:15 +0100

Yves Revaz wrote:
> Dear Gadget list,
>
> I'm trying to understand some points in the comoving
> implementation of the hydro part of Gadget2.
>
>
> In "timestep.c" when ComovingIntegrationOn is true,
>
> dt_entr = (tend - tstart) * All.Timebase_interval;
> which in fact corresponds to dt_entr = da/a (where a is the scale factor).
>
> Then, the evolution of the entropy is (as for non-comoving simulations)
> SphP[i].Entropy += SphP[i].DtEntropy * dt_entr;
>
>
> SphP[i].DtEntropy is computed in "hydra.c", where we find in
> /* do final operations on results */ :
> SphP[i].DtEntropy *= GAMMA_MINUS1 / (hubble_a2 * pow(SphP[i].Density,
> GAMMA_MINUS1));
> which means, that :
>
> SphP[i].DtEntropy)comoving = DEntropy/da
> = DEntropy/dt * 1/(H(a) a^2)
> = SphP[i].DtEntropy)non-comoving * 1/(H(a) a^2)
>
> Here, I expected to have :
> DtEntropy/da = DtEntropy/dt * 1/(H(a))

Sorry, I'm not fully sure what you mean here... in particular how you
define "DtEntropy)comoving" and a non-comoving DtEntropy.

>
> in order to have then :
> SphP[i].DtEntropy * dt_entr = DEntropy/dt * 1/(H(a)) * da/a
> = DEntropy
>
> because dt*H(a)*a = da
>
> Probably I missed a point here.
>
> It is also difficult to understand the fac_mu correction in "hydra.c":
>
> 1) 0.0001 * soundspeed_i / SphP[i].Hsml / fac_mu
> with fac_mu = a^(3/2(GAMMA - 1))/a
>
> Here I expected to find
> fac_mu = a^(1/2(GAMMA - 1))/a
> in order to write soundspeed_i / SphP[i].Hsml in non-comoving units.
>
> 2) mu_ij = fac_mu * vdotr2 / r
> is also unclear to me
>

Again, I think the confusion probably arises from not being clear about
how the variables used by the code translate to physical quantities. For
example, the "soundspeed" variable (let me denote it with c_c) in
hydra.c is related to the real physical soundspeed c_p by

         c_c = a^(3(gamma-1)/2) c_p

That's because the pressure variable used by the code, P_c, is defined
as
         P_c = A * rho_c^gamma,

where A is the (physical) entropy (denoted by Entropy in the code), and
rho_c= rho_p * a^3 is the comoving density. This means that P_c is
related to the physical pressure as

         P_c = a^(3*gamma) P_p

To work out the correct prefactors of a in various places of hydra.c, I
recommend to first identify the translation between variables used by
the code (which you may call "comoving" if you like) and the underlying
physical quantities, and then to start putting this into the basic
equations in physical units.

>
>
> Another question conserns the
> dt_hydrokick factor which corresponds to the integral
> Int_a1^a2 da/( H(a)*a*a^3(gamma-1) )
> If, following Quinn et al. 97, the origin of dt_gravkick is clear,
> I fail to understand the origin of dt_hydrokick.
>

This follows from the Hamiltonian of fluid dynamics when expressed in
the canonical coordinates (x,p), and from deriving the leap frog as
alternating solutions of kinetic and potential/thermal parts of the
Hamiltonian. If x are comoving coordinates, p=a^2 dx/dt the conjugate
momenta, and s = ln(a) the time variable, then one finds for the part of
the equation of motion due to gas pressure

   dp/ds = 1/(H(a)*a^(3*gamma-1)) [ grad_x P_c/ rho_c ]

where P_c = A*rho_c^gamma, A is physical entropy and rho_c is comoving
density.

For the "kick" in reversible hydro, all x and A are constant, hence the
term in square brackets is constant. The time integration of the
prefactor in front gives the dt_hydrokick factor.

Volker
Received on 2007-02-03 00:46:39