[mesa-users] approx13

[Bill Paxton]

Good job!   I think that's it.  That code wasn't ever intended for anything but approx20/21 --- so it never hit that particular bug.

i'll fix it.

thanks,
bill


On May 16, 2014, at 2:26 PM, Ken Shen wrote:

> Hi Bill, just realized that my "understanding" last night wasn't quite complete (well, obviously).  Here are the lines in get1_derivs in net_derivs.f where I think the high temperature factors are getting used:
> 
>          if (temp > 1d9) then
>             has_neut = (i1 == neut .or. i2 == neut .or. i3 == neut .or. &
>                o1 == neut .or. o2 == neut .or. o3 == neut)       
>             if (has_neut) then
>                call limit_rate_at_very_high_T( &
>                    ir, rates, rtab, n% factor_neut, n% d_factor_neut_dT, .false.) 
>             else if (prot /= 0) then
>                has_prot = (i1 == prot .or. i2 == prot .or. i3 == prot .or. &
>                   o1 == prot .or. o2 == prot .or. o3 == prot)                  
>                if (has_prot) then
>                   call limit_rate_at_very_high_T( &
>                      ir, rates, rtab, n% factor_prot, n% d_factor_prot_dT, .false.) 
>                end if
>             else
>                has_h1 = (i1 == h1 .or. i2 == h1 .or. i3 == h1 .or. &
>                   o1 == h1 .or. o2 == h1 .or. o3 == h1)                  
>                if (has_h1) then
>                   call limit_rate_at_very_high_T( &
>                      ir, rates, rtab, n% factor_prot, n% d_factor_prot_dT, .false.) 
>                end if
>             end if
>          end if
> 
> If neutrons are not in the isotopes list, "neut" is 0.  Then, has_neut should always be false, and factor_neut should never be used.  However, when neutrons are not in the list, factor_neut ends up being used for ALL reactions, the exact opposite of what I think you intended.
> 
> I think what is happening is that i3 and o3 are set to 0 at the beginning of get1_derivs and stay that way (because no reactions are 3 body in these nets).  So when neut == 0, has_neut is always true, regardless of the reaction (which should have has_neut = false because there are no neutrons).
> 
> I assume that's not what's supposed to happen in this piece of code, right?  Maybe there's supposed to be an "if (neut /= 0)", just like the "if (prot /= 0)" a little further down.  And maybe the same thing for h1 (although I'm confused about prot vs. h1 further up in net_derivs.f).
> 
> Or am I totally misinterpreting the code??

it is now changed to this

         if (temp > 1d9) then
            has_neut = (neut /= 0) .and. (i1 == neut .or. i2 == neut .or. i3 == neut .or. &
               o1 == neut .or. o2 == neut .or. o3 == neut)       
            if (has_neut) then
               call limit_rate_at_very_high_T( &
                   ir, rates, rtab, n% factor_neut, n% d_factor_neut_dT, .false.) 
            else if (prot /= 0) then
               has_prot = (i1 == prot .or. i2 == prot .or. i3 == prot .or. &
                  o1 == prot .or. o2 == prot .or. o3 == prot)                  
               if (has_prot) then
                  call limit_rate_at_very_high_T( &
                     ir, rates, rtab, n% factor_prot, n% d_factor_prot_dT, .false.) 
               end if
            else
               has_h1 = (h1 /= 0) .and. (i1 == h1 .or. i2 == h1 .or. i3 == h1 .or. &
                  o1 == h1 .or. o2 == h1 .or. o3 == h1)                  
               if (has_h1) then
                  call limit_rate_at_very_high_T( &
                     ir, rates, rtab, n% factor_prot, n% d_factor_prot_dT, .false.) 
               end if
            end if
         end if