[mesa-users] Large network test case

Mathieu mathren90 at gmail.com
Fri Nov 21 23:19:44 EST 2014


Hi Bill, MESA users,

Thanks for the reply, and sorry it took me so long to get back on this. 
In the future I will try the adaptive network or insert isotopes in 
smaller groups.

For the moment I would be happy to take the runs I have until a certain 
point and push them to core collapse with a larger network. Therefore, I 
am trying to restart with mesa_201 a photo/model computed until "Oxygen 
depletion" (center_o16<0.04) with approx21.
I am setting

       adjust_abundances_for_new_isos = .false.

in my star_job namelist, in order to not change my initial abundances.
After a few retries and backups, I am able to run a few iterations, but 
then all my attempts inevitably crash.

I figured that the problem is already present at the very beginning 
after I change the network: I get a large eps_nuc (>10^{10} erg/g/s) 
everywhere in the star.
I get this behavior both with MESA 6794 and 7184, and both with the 
NACRE and the Jina rates. As further evidence, I attach a print-screen 
of the first TRho pgstar plot obtained with the large network. It shows 
that the energy generation rate is >10^7 erg/g/s everywhere in the star. 
I also attach the plots of log_abs_eps_nuc (labeled as 
log_{10}(|eps_nuc-eps_nu|/ [erg g^{-1} s^{-1}), but its just 
log_abs_eps_nuc from the profile*.data file) as a function of mass for 
the first model computed with mesa_201. These plots refer to the run 
with MESA 7184, but the qualitative behavior is the same for MESA 6794 
(to which I would like to stick to finish my current project).

What follows is a description of my attempts to understand what is going 
on and their results.
All the energy generation comes from the nuclear reaction category 
labeled "other", so I think the problem is in the routine get1_derivs in 
$MESA_DIR/net/private/net_derivs.f . To make sure the problem is there, 
I added the following lines after " call get1_derivs " in net_derivs.f :

                      print *, 'case other -> calling get1_derivs ', &
                               'eps_nuc:', eps_nuc_MeV(i_rate)*Qconv
                      print *, 'reaction id:', reaction_id(i), &
                                reaction_Name(reaction_id(i))

(note eps_nuc = eps_nuc_MeV(i_rate)*Qconv, from net_eval.f), and the 
reaction names printed out correspond to electron or positron captures 
reactions.
What follows is an example of the output I get:

          case other -> calling get1_derivs eps_nuc: 
2988594.50685563281392040427208323697
          reaction id:         548 r_na23_wk-minus_mg23
          case other -> calling get1_derivs eps_nuc: 
3364551693.97504178215533047145691508
          reaction id:        1690 r_fe55_wk_mn55

(note that not all the eps_nuc_MeV * Qconv are high as these, but some 
are and they are - likely - the cause of my problem).
The high values of eps_nuc_MeV(i_rate)*Qconv come both from reaction in 
weaklib and reactions not in weaklib.
I have the impression that the problem could be in the "do_one_one_neu" 
routine in $MESA_DIR/net/private/net_derivs_support.f. I added in that 
routine the following lines just after the do loop that sets all the 
eps_nuc_MeV(:) elements. "tmp" is just the value of eps_nuc_MeV(i_rate) 
before the loop. I would expect it to be zero, since eps_nuc_MeV is set 
to zero at the beginning of get_derivs, and then it is not modified 
before the call to get1_derivs, and it is not modified either in 
get1_derivs before the call to do_one_one_neu (am I wrong?)

          if (eps_nuc_MeV(i_rate) .GT. (10.0/Qconv)) then
             print *, '***************',reaction_Name(n% g% 
reaction_id(i)), '***************'
             print *, 'eps_nuc_MeV=', eps_nuc_MeV(i_rate)
             print *, 'eps_nuc=',  eps_nuc_MeV(i_rate)*Qconv
             print *, 'Qconv=', Qconv
             print *, 'initialized value:', tmp
             print *, 'Q=', Q
             print *, 'rvs(i_rate)=', rvs(i_rate), 'r=', r, &
                      'rate_screened=', n% rate_screened(i)
!read(*,*)
          end if

Here is an example of the output:

          ***************r_zn65_wk_cu65 ***************
          eps_nuc_MeV= 2.14636470432912905281222502591286025E-0008
          eps_nuc=   20709272777.4359010681509838041751379
          Qconv=   9.6485339773180915E+017
          initialized value:   2.1463647043291291E-008
          Q=  0.97653930168084468
          rvs(i_rate)=   0.0000000000000000      r= 
0.0000000000000000      rate_screened=   7.9069687275509169E-010

So, it looks that the "initialized" value is not set to zero and causes 
an extremely high eps_nuc. But, this may (or maybe not?) only be part of 
the problem...
I tried to modify the do loop of do_one_one_neu with the naive:

          do j=1,num_rvs
             eps_nuc_MeV(j) = 0.0d0 !initialize
             eps_nuc_MeV(j) = eps_nuc_MeV(j) + rvs(j)*Q
          end do

This does something, but it doesn't fix the issue: the TRho profile now 
doesn't show the red border anywhere (which is probably not good 
either), and the "other" category in my profile files still shows the 
~10^10 eps_nuc values. Since I am not even completely sure that the 
initialization of the eps_nuc_MeV is the fundamental issue here, I am a 
bit lost. I apologize if my explanation is not clear...
I also attach the inlist (loaded after the inlist_massive_defaults - in 
which I set Zbase = 0.019) I am using for this run. The initial model is 
a large file (1.1 Mb), so I put it here for downloads:

     http://www.tapir.caltech.edu/~mathren90/O_depl.mod

I really could use some help to understand what is going on, why and 
possibly how to fix it.
Best regards,

Mathieu



On 11/18/2014 01:05 PM, Bill Paxton wrote:
> Hi Mathieu,
>
> The main thing is to add isotopes before they become active in nuclear reactions.
> Ideally, everything would run so fast, you'd just dump them all in from the beginning.
> But that's still a distant dream for now, so you'll want to creep up on it instead.
> Start with a small-ish net and add isos as you move towards advanced burning.
>
> You can just do this by hand (that's what I did for the test cases).
> Or you can try using an "adaptive network" based on Alex Heger's scheme.
> For that option, take a look at the "adjust_net" test case
>
> Cheers,
> Bill
>
>
>
> On Nov 18, 2014, at 12:48 PM, Mathieu wrote:
>
>> Hi everyone,
>>
>> I have a question regarding the initial model for all the test cases
>> using the large mesa_*.net network
>> ($MESA_DIR/star/test_suite/*_big_net). How where the initial model
>> computed? Where they computed from ZAMS using the large network, or
>> where they computed until some point with a small network and then the
>> network was changed? And if it was changed, how did this happen?
>>
>> I am trying to understand what is the best way to change nuclear
>> network, to avoid running in problems because of the initialization of
>> the abundances of the new isotope introduces (especially if there are a
>> lot of new isotopes.)
>>
>> Thanks!
>>
>> Mathieu
>>
>>
>>
>>
>> ------------------------------------------------------------------------------
>> Download BIRT iHub F-Type - The Free Enterprise-Grade BIRT Server
>> from Actuate! Instantly Supercharge Your Business Reports and Dashboards
>> with Interactivity, Sharing, Native Excel Exports, App Integration & more
>> Get technology previously reserved for billion-dollar corporations, FREE
>> http://pubads.g.doubleclick.net/gampad/clk?id=157005751&iu=/4140/ostg.clktrk
>> _______________________________________________
>> mesa-users mailing list
>> mesa-users at lists.sourceforge.net
>> https://lists.sourceforge.net/lists/listinfo/mesa-users



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20141121/55dd12cc/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: eps_first_large_net.png
Type: image/png
Size: 48265 bytes
Desc: not available
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20141121/55dd12cc/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: TRho_pgstar_screenshot.png
Type: image/png
Size: 13205 bytes
Desc: not available
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20141121/55dd12cc/attachment-0001.png>
-------------- next part --------------
!after inlist_massive_defaults

&star_job

      mesa_dir = '' ! empty string defaults to environment variable MESA_DIR
!      history_columns_file = './history_columns.list' ! if null string, use default.
!      profile_columns_file = './profile_columns.list' ! if null string, use default.      

      show_log_description_at_start = .true.
      show_net_species_info = .true.
      show_net_reactions_info = .true.	  
      
      create_pre_main_sequence_model = .false.

      change_net = .true.      
      new_net_name = 'mesa_201.net'
      adjust_abundances_for_new_isos = .false.    

      set_rates_preference = .true. ! for use by net + rates modules
      new_rates_preference = 2 ! only used if set_rates_preference is true
      ! 1 = NACRE rates -- this is the default
      ! 2 = jina reaclib rates

      eos_file_prefix = 'mesa'
      kappa_file_prefix = 'gs98'

      write_profile_when_terminate = .true. ! write to a given name
      filename_for_profile_when_terminate = 'onset_of_cc_big_net.data'
      save_photo_when_terminate = .true.
      pause_before_terminate = .false.
      pgstar_flag =.true.
      
      load_saved_model = .true.
      saved_model_name = 'O_depl.mod'

/ !end of star_job

&controls

    report_all_dt_limits = .true.
     
    ! output to files and terminal
    photostep = 1
    profile_interval = 1
    history_interval = 1
    terminal_interval = 1
    write_header_frequency = 1
    max_num_profile_models = 100
    
   !Mass Loss
!    use_other_wind = .true.
!    RGB_wind_scheme = 'other'
!    AGB_wind_scheme = 'other'
!    Dutch_wind_eta = 1.0
   
    mixing_length_alpha = 2.0
    MLT_option = 'Henyey'
    
    alpha_semiconvection = 0.1
    
    overshoot_f_above_nonburn = 0.025
    overshoot_f_above_burn_h = 0.025
    overshoot_f_above_burn_he = 0.025
    overshoot_f_above_burn_z = 0.025
    overshoot_f_below_nonburn = 0.025
    overshoot_f_below_burn_h = 0.025
    overshoot_f_below_burn_he = 0.025
    overshoot_f_below_burn_z =  0.025
    
    ! max_allowed_nz = 16000
    ! mesh_delta_coeff_for_highT = 0.7
    ! logT_max_for_standard_mesh_delta_coeff = 9.0 
    ! logT_min_for_highT_mesh_delta_coeff = 9.5 
    
    ! varcontrol_target = 1d-3
    ! mesh_delta_coeff= 0.7
    ! mesh_max_allowed_ratio=3.5
    ! min_timestep_limit = 1d-12 !in sec
    ! when to stop
    fe_core_infall_limit = 1d8 
    
/
    



More information about the Mesa-users mailing list