[mesa-users] controlling timestep using diffusion timescale

Fri Mar 6 09:46:32 EST 2015

Hi Bill,

> really??  can you provide me with inlist etc so I can see this in action.
> when I look at the code, I find subroutine do_solve_diffusion in star/private/diffusion.f90 doing this to save the number of steps taken with isolve:
>             steps_used = iwork(17) ! number of accepted steps
> and the caller of do_solve_diffusion (in star/private/element_diffusion) does this
>          s% num_diffusion_solver_steps = steps_used
> in star/private/timestep, the routine check_diffusion_steps_limit then uses s% num_diffusion_solver_steps to enforce the limits set by the controls diffusion_steps_limit and diffusion_steps_hard_limit.
>
> so if this isn't working, please give me details so we can fix it.

The value of diffusion_dt_div_timescale seems only to be used in the
do_step_loop routine, when using isolve its not taken into account at
all. The attached inlist shows this, modifying
diffusion_dt_div_timescale between 1d99 and 1d-99 makes no difference at
all.

>
>
>> Also, would it be possible to have the value computed by
>> get_timescale(...) in star/private/mod_diffusion.f90 accesible through a
>> public interface, and provide a timestep limit associated to diffusion
>> using this timescale (something like dt_div_diffusion_timescale_limit)?
>
> sure.  but keep in mind that when you use isolve for this, it can take multiple substeps, each of which can be much smaller than the timescale calculated by this routine.   this routine is for use in case you are not using isolve.  and it is called at each substep taken by the diffusion solver, and the value changes a lot as the initial transients are smoothed out allowing larger substeps to be taken.  so it isn't at all clear what you can usefully do with a particular value of this.    we can of course save the final value of get_timescale in a star_data variable so you can inspect it from your run_star_extras.
>
> You can give that a try for yourself by making a patch.  Just after the call on get_timescale, add this line
>
> 	s% xtra1 = timescale
>
> Then you can access s% xtra1 in your run_star_extras.
>
> If this turns out to be a useful way to set timesteps, let me know and we'll add something so you won't have to do a patch for future releases.

The idea here is to use this  diffusion timescale to reduce the
evolutionary step.  So, if I have the timescale stored in xtra1, I can
add to extras_finish step in run_star_extras.f :   s% dt_next = min(s%
dt_next, s%xtra1*0.01), and that would give me the desired limit. I
modified mod_diffusion to store the computed timescale into xtra1,
however, since for isolve the timescale is not computed (could not find
more calls to get_timescale), I added this extra code to the routine

            call update_coeffs((kmax_rad_accel > 0 .and. steps_used >
1), ierr)
            call get_timescale( &
               s, nz, nzlo, nzhi, m, nc, -s%
diffusion_min_X_hard_limit*0.5d0, v_advection_face, &
               upwind_limit, X, X_face, C_div_X, SIG_face, GT_face,
AD_face, cell_dm, r_face, &
               steps_used, total_num_iters, dbg, iter_dbg, j_dbg, k_dbg, &
               j, k, class_chem_id, total_time, timescale)

             s% xtra1 = timescale

I don't know exactly what update_coeffs does, so I might be breaking
something, but without it no timescale is computed (I essentially copied
this from do_step_loop). I attach run_star_extras.f and also
mod_diffusion.f ; the run_star_extras has the code to limit timestep
with xtra1, and mod_diffusion has the changes to store the timescale in
xtra1.

This seems to work well, the computed timescale for the 1.2 solar mass
star changes more or less continuously, and demanding the timestep to be
around 1d-2 or 1d-3 of this value gives me the value for timestep I had
before when using max_years_for_timestep.

Cheers,

Alina

> Cheers,
> Bill
>
>
>
>
>> Cheers!
>>
>> Alina
>>
>>
>> -- 
>> Alina Istrate 
>> PhD student
>>
>>
>> Argelander Institute fur Astronomie
>> Auf dem Hugel 71
>> D-53121 Bonn
>> Germany
>> Room 2.024
>>
>>
>> ------------------------------------------------------------------------------
>> Dive into the World of Parallel Programming The Go Parallel Website, sponsored
>> by Intel and developed in partnership with Slashdot Media, is your hub for all
>> things parallel software development, from weekly thought leadership blogs to
>> news, videos, case studies, tutorials and more. Take a look and join the 
>> conversation now. http://goparallel.sourceforge.net/
>> _______________________________________________
>> mesa-users mailing list
>> mesa-users at lists.sourceforge.net
>> https://lists.sourceforge.net/lists/listinfo/mesa-users

-- 
Alina Istrate 
PhD student

Argelander Institute fur Astronomie
Auf dem Hugel 71
D-53121 Bonn
Germany
Room 2.024

-------------- next part --------------

&star_job

      show_log_description_at_start = .false. 

 !   create_pre_main_sequence_model = .true. 

      save_model_when_terminate = .true.
      save_model_filename = 'final.mod'

      write_profile_when_terminate = .true.
      filename_for_profile_when_terminate = 'final_profile.data'

     pgstar_flag = .true.
/ ! end of star_job namelist

&controls

      history_interval = 1
      extra_terminal_output_file = 'log1'
      photo_directory = 'photos1'
      log_directory = 'LOGS1'

      max_model_number = 500000
      max_num_profile_models = 100000  

      initial_mass = 1.2
      initial_z = 0.02
      max_allowed_nz = 10000
      photostep = 50 
      profile_interval = 100 
      history_interval =1 
      terminal_cnt = 1
      write_header_frequency = 10
      varcontrol_target = 1d-5
      mesh_delta_coeff = 1.0

     Teff_upper_limit =8d3
     !max_years_for_timestep= 1d5
      show_diffusion_info = .true.
      diffusion_dt_div_timescale= 1d-99
     ! diffusion_steps_limit = 5 
     ! diffusion_iters_limit = 5  

      newton_itermin = 4  
      !report_hydro_solver_progress =.true.

    ! limit_mixing_length_by_dist_to_bdy = 1.0 
    ! conv_bdy_mix_softening_f0 = 0.001
    !  conv_bdy_mix_softening_f = 0.003
    !  conv_bdy_mix_softening_min_D_mix = 0.00

     xa_function_species(1) = 'h1'
     xa_function_weight(1) =200
     xa_function_param(1) =1d-2

      step_overshoot_f_above_burn_h = 0.045
      step_overshoot_f_below_burn_h = 0
      overshoot_f0_above_burn_h = 0.005
      overshoot_f0_below_burn_h = 0
      step_overshoot_D0_coeff = 1.0

      num_trace_history_values = 2
      trace_history_value_name(1) = 'envelope_mass'
      trace_history_value_name(2) = 'star_mass_h1' 
      mixing_length_alpha = 1.89
      use_Ledoux_criterion =.true.
      alpha_semiconvection = 1d-3

      report_why_dt_limits = .true. 
      diffusion_dt_limit =  3.15d7
      do_element_diffusion = .true. ! determines whether or not we do diffusion
      diffusion_calculates_ionization = .true.
      diffusion_tol_correction_max = 1d-1
      diffusion_tol_correction_norm = 1d-3

      diffusion_num_classes = 8 ! number of classes of species for diffusion calculations
      diffusion_class_representative(1) = 'h1'
      diffusion_class_representative(2) = 'he3'
      diffusion_class_representative(3) = 'he4'
      diffusion_class_representative(4) = 'c12'
      diffusion_class_representative(5) = 'n14'
      diffusion_class_representative(6) = 'o16'
      diffusion_class_representative(7) = 'ne20'
      diffusion_class_representative(8) = 'mg24'

      ! in ascending order.  species goes into 1st class with A_max >= species A
      diffusion_class_A_max(1) = 2
      diffusion_class_A_max(2) = 3
      diffusion_class_A_max(3) = 4
      diffusion_class_A_max(4) = 12
      diffusion_class_A_max(5) = 14
      diffusion_class_A_max(6) = 16
      diffusion_class_A_max(7) = 20
      diffusion_class_A_max(8) = 10000

      diffusion_use_isolve = .true.
      diffusion_rtol_for_isolve = 1d-4
      diffusion_atol_for_isolve = 1d-5
      diffusion_maxsteps_for_isolve = 1000
      diffusion_isolve_solver = 'ros2_solver'

    !  report_hydro_solver_progress = .true.
  !    tol_correction_norm = 1d-5
  !    tol_max_correction = 1d-3

/ ! end of controls namelist

&pgstar

    HR_win_flag = .true.
    !Kipp_win_flag = .true.
/ ! end of pgstar namelist
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mod_diffusion.f
Type: text/x-fortran
Size: 92916 bytes
Desc: not available
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20150306/96fb36e1/attachment.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: run_star_extras.f
Type: text/x-fortran
Size: 8248 bytes
Desc: not available
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20150306/96fb36e1/attachment-0001.bin>