[FLASH-USERS] Bug in Ionize module?

[Klaus Weide]

On Mon, 17 Oct 2016, Suoqing Ji wrote:

> Hi,
> 
> I’m looking at the Ionize module and suspect there might be a bug in IonizeMain/Eqi/nei.F90 (perhaps this file is out of date). I think the line
> 
> xout(jj2) = xout2(i)/selnel
> 
> should be
> 
> xout(jj2) = xout2(i) * selnel(jj2)
> 
> Even from the view of code grammar, selnel is defined as a single real number, but it refers to an array ion_xfrac when nei.F90 is called from Ioinze.F90. 

There are two versions of the file nei.F90 in the Ionize unit:

  IonizeMain/Eqi/nei.F90
  IonizeMain/Nei/nei.F90

Each contains an implementation of

  subroutine neimn

but their interfaces are incompatible with each other. In particular, 
compare declarations for the last dummy argument:

  IonizeMain/Nei/nei.F90:  real, intent(in) :: xfrac(*)     ! an array
  IonizeMain/Eqi/nei.F90:  real,    intent(in) ::  selnel   ! a scalar

The declaration in IonizeMain/Nei matches the call in Ionize.F90, but the 
one in IonizeMain/Eqi doesn't (as you have found). 

The combination of IonizeMain/Ionize.F90 and IonizeMain/Eqi/nei.F90 does 
not result in a valid Fortran program; but compilers will not usually 
detect this because the code does not use an explicit interface for neimn.

We are running regular regression tests using the NeiTest simulation, but 
those only use the Nei variant.

It appears that IonizeMain/Eqi/nei.F90 is severely out of date;
it also appears that two incompatible versions of nei.F90 were already 
part of FLASH version 2.  So it seems the directory IonizeMain/Eqi is 
mostly of historical value, and the Users Guide is wrong as far as it 
gives the impression that Eqi is a viable alternative to Eqi.

I am not going to discuss under which physical conditions the approach in 
the Ionize unit is applicable / valid / useful; I hope somebody else can 
clarify if needed, and answer your question:

> And I have an additional question: if I understood this correctly, the 
> electron population fraction is set to constant in the code (by electron 
> to proton ratio). But even if we can ignore the contribution of metals, 
> helium might not be fully ionized and may make a significant change to 
> electron density. Then why is it valid that electron density remain 
> unchanged throughout ionization?


Klaus