[FLASH-BUGS] Bug in iso7 reaction network

[Lynn Reid]

Chris:

Many thanks for your email with the astute discovery of a bug in the 
iso7 network.  I have changed the source code in subroutine diso7 in 
FLASH2/source/source_terms/burn/iso7/net.F90, and made an equivalent
change in FLASH3's version at 
FLASH3/source/physics/sourceTerms/Burn/BurnMain/nuclearBurn/iso7/bn_networkDenseJakob.F90.  
I have noted in the source code your squashing of the bug, and added 
your email below
to the change logs for both FLASH2 and FLASH3. 

Best wishes,

Lynn Reid

Chris Malone wrote:
>
> I was recently trying to modify the iso7 reaction network when I 
> discovered a bug in the calculation of the dense jacobian element 
> pertaining to the triple alpha reaction.  When calculating the change 
> in composition we have:
>
> dydt(ihe4) = -3.0e0 * y(ihe4) * y(ihe4) * y(ihe4) * ratdum(ir3a) + ...
>
> The -3.0e0 is from the fact that we are losing 3 helium nuclei per 
> triple alpha reaction and the reaction rate has already been divided 
> by 3! = 6 to prevent double counting.  Now to calculte the dense 
> jacobian matrix element, dfdy(ihe4, ihe4), in diso7 (in FLASH version 
> 2) we should simply take the above expression and differentiate with 
> respect to helium composition giving:
>
> dfdy(ihe4, ihe4) = -9.0e0 * y(ihe4) * y(ihe4) * ratdum(ir3a) + ...
>
> However, in the FLASH implementation of iso7 the coefficient is -6.0e0 
> not -9.0e0.  I brought this to the attention of Frank Timmes and he 
> said indeed it was a bug and that he was not suprised that the bug 
> propagated through to FLASH because the paper he wrote that describes 
> this reaction network also had the error.
>
> Frank also commented that the bug does not affect the results of the 
> calculation for the change in helium abundance but rather the 
> efficiency of the algorithm.
>
> Chris


-- 
Dr. Lynn B. Reid, Scientific Programmer
ASC Center for Astrophysical Thermonuclear Flashes 
The University of Chicago
5640 S. Ellis, RI 473; Chicago, IL 60637 USA
+1-773-702-0624 phone