[FLASH-USERS] bug in eos_mgamma.F90 file

[Klaus Weide]

On Thu, 3 Dec 2015, jeevan dahal wrote:

> Hello,
> 
> I noticed a bug in the EoS subroutine of Flash 4.3. In the file
> eos_mgamma.F90 (source/physics/Eos/EosMain/Multigamma/eos_mgamma)  line
> 377  the code is
> 
> 
>          eosData(det+ilo:det+ihi) = eos_gasConstant /
> eosData(abar+ilo:abar+ihi)* &
>                 (eosData(gamc+ilo:gamc+ihi) - 1.0)
> 
> 
> Again, in line 429,
> 
>  eosData(c_v+ilo:c_v+ihi) = eosData(det+ilo:det+ihi)

Jeevan,

This appears to be wrong indeed. Thank you for catching and reporting this 
error!

As a result of this wrong code, the values of the following quantities 
returned by Eos calls would be wrong by a factor of (gamma-1)**2 :

   derived EOS_DET Derivative of internal energy wrt temperature
   derived EOS_CV  Specific heat at constant volume

This error appears to be present only in the 1T Multigamma implementation 
of Eos. I checked and did not see it in the multiTemp variant of 
Multigamma (or Multitype), nor in other 1T Eos implementations (Gamma, 
Helmholtz).

This error has been present in all versions of FLASH 3 and FLASH 4.
It has probably escaped detection over several years because
 (1) The (1T) Multigamma Eos is not widely used -
     Commonly, for simple testing an even simpler Eos ("Gamma") is used,
     and for more realistic simulations, a more physical Eos is used
     (like "Helmholtz");
 (2) EOS_DET and EOS_CV are optional outputs from the Eos that are usually
     not requested (and not independently tested);
 (3) Most code outside the Eos unit that does use the EOS_DET or EOS_CV 
     values returned from Eos calls requires a different Eos 
     implementation and thus would never be used in combination with the 
     1T Multigamma Eos. 


> This does not give the correct value of CV. There seems to be missing
> bracket in calculation of  eosData(det+ilo:det+ihi) which can be corrected
> by writing code as:
> 
>  eosData(det+ilo:det+ihi) = eos_gasConstant / (eosData(abar+ilo:abar+ihi)*
> &
>                 (eosData(gamc+ilo:gamc+ihi) - 1.0))

This should work as a fix.  We will apply an equivalent fix to our code.

Klaus