Describe the issue:
The drag law size is coded wrong in drag.hpp
Right now, it is
#if HAVE_ENERGY == 1 cs = std::sqrt(eos.GetGamma(VcGas(PRS,k,j,i),VcGas(RHO,k,j,i) *VcGas(PRS,k,j,i)/VcGas(RHO,k,j,i)));
which amounts to $c_s = \sqrt{ \gamma (P,P) }$.
A parenthesis should be moved, so that it becomes:
#if HAVE_ENERGY == 1 cs = std::sqrt(eos.GetGamma(VcGas(PRS,k,j,i),VcGas(RHO,k,j,i)) *VcGas(PRS,k,j,i)/VcGas(RHO,k,j,i));
which amounts to $c_s = \sqrt{ \gamma (P,\rho)P / \rho }$.
Error message:
runtime information:
To my knowledge, this is present on all versions of the code where this function is implemented.
In practice, for simulations with a constant $\gamma $ the sound speed is taken close or equal to unity.
The size drag law then becomes another gamma drag law, with $\gamma_i = 1/ \beta_i $.
I have made a couple tests with the parenthesis at the right place, such as adding dust to the AmbipolarWind test setup and checking the dust drift against terminal velocity. I found a good agreement between the two..
Describe the issue:
The drag law size is coded wrong in drag.hpp
$c_s = \sqrt{ \gamma (P,P) }$ .
$c_s = \sqrt{ \gamma (P,\rho)P / \rho }$ .
Right now, it is
#if HAVE_ENERGY == 1 cs = std::sqrt(eos.GetGamma(VcGas(PRS,k,j,i),VcGas(RHO,k,j,i) *VcGas(PRS,k,j,i)/VcGas(RHO,k,j,i)));which amounts to
A parenthesis should be moved, so that it becomes:
#if HAVE_ENERGY == 1 cs = std::sqrt(eos.GetGamma(VcGas(PRS,k,j,i),VcGas(RHO,k,j,i)) *VcGas(PRS,k,j,i)/VcGas(RHO,k,j,i));which amounts to
Error message:
runtime information:
To my knowledge, this is present on all versions of the code where this function is implemented.$\gamma $ the sound speed is taken close or equal to unity.$\gamma_i = 1/ \beta_i $ .
In practice, for simulations with a constant
The size drag law then becomes another gamma drag law, with
I have made a couple tests with the parenthesis at the right place, such as adding dust to the AmbipolarWind test setup and checking the dust drift against terminal velocity. I found a good agreement between the two..