[Mesa-users] Timestep limits by negative abundances

[Robert Willer]

Hello MESA-Community!

My Name is Robert, I am working on models of classical novae with MESA, 
including overshooting. I am encountering the problem that the timestep 
of my models is limited by negative mass fractions being produced 
(particularly c12, after the initial thermonuclear runaway, as the nova 
envelope is expanding). Do you have any suggestions where I should look 
or what I could try to mitigate this? I have attached the inlist I am 
using, its based on the wd_nova_burst test case and I apologize for it 
being a bit messy, its a work in progress.

Many thanks, Robert



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

! inlist_wd_nova_burst


&star_job

      show_log_description_at_start = .false.
      
      load_saved_model = .true.
      load_model_filename = 'ready.mod' ! from running inlist_setup

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

!      save_model_when_terminate = .true.
!      save_model_filename = 'start_of_wd_nova_burst.mod'
!      required_termination_code_string = 'max_model_number'

      set_initial_age = .true.
      initial_age = 0
 
      set_initial_model_number = .true.
      initial_model_number = 0

      !set_tau_factor = .true.
      !set_to_this_tau_factor = 30
    
     ! jina_reaclib_min_T9 = 3d-3


      
      change_v_flag = .true.
      new_v_flag = .true. !.true.

     ! use_suzuki_weak_rates = .true.    ! Suzuki, Toki, and Nomoto (2016)
      !     Electron-capture and $\beta$-decay rates for sd-shell nuclei in stellar environments relevant to high-density O-Ne-Mg cores
      !     http://adsabs.harvard.edu/abs/2016ApJ...817..163S

      ion_coulomb_corrections = 'PCR2009'
      electron_coulomb_corrections = 'Itoh2002'


      set_initial_cumulative_energy_error = .true.
      new_cumulative_energy_error = 0d0

      set_initial_dt = .true.
      years_for_initial_dt = 1d4
      
      !remove_initial_center_by_mass_fraction_q = 0.5

      pgstar_flag = .true.

      pause_before_terminate = .true.

!Set up the color corrections by using the appropriate files

     ! color_num_files = 2
     ! color_file_names(1) = 'lcb98cor.dat'
     ! color_num_colors(1) = 11
      
     ! color_file_names(2) = 'blackbody_johnson.dat'
     ! color_num_colors(2) = 5



/ ! end of star_job namelist


&eos
/ ! end of eos namelist


&kap
      Zbase = 0.01d0

      use_Type2_opacities = .true.

/ ! end of kap namelist


&controls

      ! limit max_model_number as part of test_suite
      !max_model_number = 223
     ! max_model_number = 2250

       max_num_profile_models = 3000
      
! solver

      energy_eqn_option = 'eps_grav' ! 'eps_grav' !'eps_grav'
      
     ! use_gold2_tolerances = .false.

      limit_for_rel_error_in_energy_conservation = 1d-5
      hard_limit_for_rel_error_in_energy_conservation = 1d-3

      warning_limit_for_max_residual = 1d99
      !convergence_ignore_equL_residuals = .true.
      ignore_too_large_correction = .true.

      neg_mass_fraction_hold = 1
      timestep_factor_for_retries = 0.8d0 !0.5d0

!line driven wind

      cool_wind_full_on_T = 0.8d4
      hot_wind_full_on_T = 1.2d4
      cool_wind_RGB_scheme = 'Reimers'
      cool_wind_AGB_scheme = 'Blocker'
      RGB_to_AGB_wind_switch = 1d-4
      Reimers_scaling_factor = 0.5 !0.5d0
      Blocker_scaling_factor = 0.1 !0.5d0 ! 0.8d0

      hot_wind_scheme = 'Dutch'
      Dutch_scaling_factor = 0.8

      mass_change_full_on_dt = 1d-99 ! (seconds)
      mass_change_full_off_dt = 1d-99 ! (seconds)

      ! NOTE: for super eddington wind,
      ! we use Ledd averaged over mass to optical depth tau = 100
      super_eddington_scaling_factor = 1
         ! parameter for mass loss driven by super Eddington luminosity        
     ! super_eddington_wind_Ledd_factor = 1
         ! multiply Ledd by this factor when computing super Eddington wind
         ! e.g., if this is 2, then only get wind when L > 2*Ledd
      
       surf_avg_tau = 200 !130
       surf_avg_tau_min = 1 !20 !31     

!Convection

      overshoot_scheme(1) = 'exponential'
      overshoot_zone_type(1) = 'any'
      overshoot_zone_loc(1) = 'any'
      overshoot_bdy_loc(1) = 'any'
      overshoot_f(1) = 0.010 !0.014    ! How far out from the start of overshooting to go for the exponential decay in units of H_p
      overshoot_f0(1) = 0.006 !0.004   ! How far back into the convective zone to go to start overshooting in units of scale heights
      overshoot_D_min = 1d0 ! 1d2      
      min_overshoot_q = 0.2 ! allows only the outer part (80% by mass if 0.2) to have overshooting to avoid bogus overshooting in the core)

        ! simplify convection
       redo_conv_for_dr_lt_mixing_length = .true.

      ! in units of the scale height
        min_convective_gap = 1d-2 !-1
       min_thermohaline_gap = 1d-2 !-1
       min_thermohaline_dropout = 1d-2 !-1
       min_semiconvection_gap = 1d-2 !-1
       remove_embedded_semiconvection = .true.

      prune_bad_cz_min_Hp_height = 0.01
      
      prune_bad_cz_min_log_eps_nuc = 1
     
    !  remove_small_D_limit = 1d0 !1d-6
      
      clip_D_limit = 1d0

      max_conv_vel_div_csound = 1d0 !1d99 !1d99 !1d0

      mixing_length_alpha = 1.9 !2.0 !1.9
      MLT_option = 'Henyey' !'TDC'
      
   !   alpha_TDC_DAMP = 1d0
   !   alpha_TDC_DAMPR = 0.1d0
   !   alpha_TDC_PtdVdt = 1d0
   !   steps_before_use_TDC = 0

      
       make_gradr_sticky_in_solver_iters = .true. !.true.
      
  !    min_logT_for_make_gradr_sticky_in_solver_iters = 5d6

      !min_logT_for_make_gradr_sticky_in_solver_iters = 1d7
     
      alt_scale_height_flag = .true. !.false. !.true. !.false.     

      use_Ledoux_criterion = .true.
      
      semiconvection_option = 'Langer_85'
      alpha_semiconvection = 0.1

      thermohaline_option =  'Brown_Garaud_Stellmach_13'
      thermohaline_coeff = 1d0

      num_cells_for_smooth_gradL_composition_term = 0 !3 !avoid smoothing composition gradients, which might cause
                                                      !convective cores to unphysically shrink.
                                                       !See section 3.2 in Moore, K., & Garaud, P. 2016, APJ, 817, 54      
      do_element_diffusion = .true. !.true.
      diffusion_use_full_net = .true. !.false.
      
      do_WD_sedimentation_heating = .false.
      do_diffusion_heating = .true.
      
      show_diffusion_info = .true. ! terminal output for diffusion

      !diffusion_T_full_on = 5d5
      !diffusion_T_full_off = 2d4
      
   
     ! D_mix_ignore_diffusion = 5d2

      diffusion_v_max = 1d99 !1d-4

     ! diffusion_min_dq_at_surface = 1d-10
      diffusion_min_T_at_surface = 5d4
       
      diffusion_dt_limit =  3.15d6 !300 !4d4 !86400 !(one day) !3.15d7 !(3.15d7 seconds is 1 year) no element diffusion if dt < this limit (in seconds)
    
      diffusion_use_cgs_solver = .true.
      diffusion_use_isolve = .true.
      diffusion_rtol_for_isolve = 6d-4 ! 1d-4
      diffusion_atol_for_isolve = 6d-5 ! 1d-5
     ! diffusion_min_X_hard_limit = -5d-3
      diffusion_maxsteps_for_isolve = 1000
      diffusion_isolve_solver = 'ros3pl_solver'

      ! Timestep controls to prevent steps that are difficult for diffusion.
     ! diffusion_steps_limit = 250 !500
     ! diffusion_steps_hard_limit = 750 !1750
     ! diffusion_iters_limit = 10
     ! diffusion_iters_hard_limit = 50

    !  diffusion_min_num_substeps = 5
      
    !  diffusion_max_iters_per_substep = 12 !50
      
 
    ! diffusion_dt_div_timescale = 0.2

    diffusion_num_classes = 14

    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) = 'ne22'
    diffusion_class_representative(9) = 'na23'
    diffusion_class_representative(10) = 'mg24'
    diffusion_class_representative(11) = 'al27'
    diffusion_class_representative(12) = 'si28'
    diffusion_class_representative(13) = 's32'
    diffusion_class_representative(14) = 'fe56'


    diffusion_class_A_max(1) = 2
    diffusion_class_A_max(2) = 3
    diffusion_class_A_max(3) = 4
    diffusion_class_A_max(4) = 13
    diffusion_class_A_max(5) = 15
    diffusion_class_A_max(6) = 18
    diffusion_class_A_max(7) = 20
    diffusion_class_A_max(8) = 22
    diffusion_class_A_max(9) = 23
    diffusion_class_A_max(10) = 26
    diffusion_class_A_max(11) = 27
    diffusion_class_A_max(12) = 30
    diffusion_class_A_max(13) = 34
    diffusion_class_A_max(14) = 1000 
                                                
!predictive mixing

   ! do_conv_premix = .true.
   ! conv_premix_time_factor = 0.5

      !any other convective zone
     ! predictive_mix(1) = .true.
     ! predictive_zone_type(1) = 'any'
     ! predictive_zone_loc(1) = 'any'
     ! predictive_bdy_loc(1) = 'any'
     ! predictive_bdy_q_min(1) = 0d0
     ! predictive_bdy_q_max(1) = 1d0
     ! predictive_superad_thresh(1) = 0.001
     ! predictive_avoid_reversal(1) = ''

     ! predictive_mix(2) = .true.
     ! predictive_zone_type(2) = 'any'
     ! predictive_zone_loc(2) = 'surf'
     ! predictive_bdy_loc(2) = 'any'
     ! predictive_bdy_q_min(2) = 0d0
     ! predictive_bdy_q_max(2) = 1d0
     ! predictive_superad_thresh(2) = 0.001             
   
      mass_change = 1d-9
      accrete_same_as_surface = .false. 
      accretion_h1 = 0.7100 ! Estimate Nomoto 2007
      accretion_h2 = 0
      accretion_he3 = 0.292910D-04
      accretion_he4 = 0.279971 ! Estimated Nomoto 2007
      accretion_zfracs = 4 ! Lodders 03
      
      min_timestep_limit = 1d-24

! functions to give extra weight to things when choosing the mesh

      xa_function_species(1) = 'h1'  ! name of nuclide as defined in chem_def
      xa_function_weight(1) = 24  !24
      xa_function_param(1) = 1d-6

      xa_function_species(2) = 'he4'  ! name of nuclide as defined in chem_def
      xa_function_weight(2) = 12 !12
      xa_function_param(2) = 1d-6


    !  xa_function_species(3) = 'n15'  ! name of nuclide as defined in chem_def
    !  xa_function_weight(3) = 24 
    !  xa_function_param(3) = 1d-8

    !  xa_function_species(4) = 'n15'  ! name of nuclide as defined in chem_def
    !  xa_function_weight(4) = 6
    !  xa_function_param(4) = 1d-8

    !  xa_function_species(5) = 'o14'  ! name of nuclide as defined in chem_def
    !  xa_function_weight(5) = 6
    !  xa_function_param(5) = 1d-8
      
    !  xa_function_species(6) = 'c13'  ! name of nuclide as defined in chem_def
    !  xa_function_weight(6) = 72
    !  xa_function_param(6) = 1d-6

     ! xa_scale = 1d-8

      min_dq_for_xa_convective = 1d-14 !1d-6
     
      mesh_min_dlnR = 1d-14 ! 1d-9

     ! merge_if_dlnR_too_small = .true.
 
     ! min_chem_eqn_scale = 1d-1 !1d0 

      log_kap_function_weight = 1.0 !0.5 ! more resolution for changes in opacity

      gam_function_weight = 5.0 !0.5 ! extra resolution at liquid/solid transition

      gradT_function_weight = 1d-4 !5d-5

      log_tau_function_weight = 1d-1

      E_function_weight = 0.1

     ! T_function1_weight = 140 !110

      !R_function_weight = 1.0

     ! R_function3_weight = 1d-1

     ! M_function_weight = 1d-1

     ! T_function2_weight = 1.0

      mesh_dlog_eps_min_for_extra = 1d0

      mesh_dlog_eps_dlogP_full_on = 5
 
      mesh_dlog_eps_dlogP_full_off = 0.5

     ! mesh_dlog_cno_dlogP_extra = 0.2d0
 
     ! delta_lgL_H_limit = 0.1 ! limit for magnitude of change in lgL_H
     ! lgL_H_burn_min = -0.5 ! ignore changes in lgL_H if value is less than this

     ! delta_lgL_He_limit = 0.25 ! limit for magnitude of change in lgL_He
     ! lgL_He_burn_min = 2.5 ! ignore changes in lgL_He if value is less than this
      
     ! delta_lgRho_limit = 1.5 ! limit for magnitude of max change in log10 density at any point
     ! delta_lgRho_hard_limit = -1

     ! delta_lgT_limit = 0.5
     ! delta_lgL_limit = 0.05
     ! delta_lgTeff_limit = 0.10 
      
     ! delta_lgL_nuc_limit = 0.02
     ! delta_lgL_nuc_hard_limit = 0.20
   !   lgL_nuc_burn_min = -3d0 ! 0.5d0
 
  !    dX_nuc_drop_min_X_limit = 1d-6 !1d-4

   !   lgL_nuc_mix_dist_limit = 1d-10  !1d-6

    !  lgL_H_burn_min = 0.5d0 !1.5d0
 
     ! lgL_He_burn_min = 0.5d0 !2.5d0
           
   !   lgL_z_burn_min = 0.5d0  !2.5d0
  
   !    op_split_burn = .true. !.true.
   !    op_split_burn_min_T = 7d7 !3d7
  !    op_split_burn_min_T_for_variable_T_solver = 1d99     
 
      op_split_burn_eps = 5d-6 !1d-5
      op_split_burn_odescal = 5d-6 !1d-5

      burn_steps_limit = 20000 !10000
      burn_steps_hard_limit = 40000 !20000


      atm_option = 'T_tau'
      atm_T_tau_relation = 'Eddington'
      atm_T_tau_opacity = 'iterated'
      atm_T_tau_max_iters = 500

      Pextra_factor = 3.0 ! 4 * tau_factor
      
      use_superad_reduction = .true.
      
      superad_reduction_limit = -2d0
      
      scale_max_correction_for_negative_surf_lum = .true. !.true.

!controls for dynamical situations

	!dynamical gradT (if hydrostatic is not true)
	!use_dPrad_dm_form_of_T_gradient_eqn = .true.
      !  use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = .true.
        
        ! time centering velocity for better energy conservation
      !  steps_before_use_velocity_time_centering = 0
      !  include_P_in_velocity_time_centering = .true.
      !  P_theta_for_velocity_time_centering = 1d0 ! 0.5d0
      !  include_L_in_velocity_time_centering = .true.
      !  L_theta_for_velocity_time_centering = 1d0 !0.5d0
      !  use_P_d_1_div_rho_form_of_work_when_time_centering_velocity = .true.

       use_Pvsc_art_visc = .true. ! .true.
  
      ! Pvsc_cq = 1.235 !4.0d0
      ! Pvsc_zsh = 0.2d0

       use_dPrad_dm_form_of_T_gradient_eqn = .false.

      !okay_to_reduce_gradT_excess = .true.
      
     ! T_mix_limit = 2d4 

      varcontrol_target = 5d-4
      
      !varcontrol_dt_limit_ratio_hard_max = 5d0 !1d99 !ration of actual varcontrol vs target, allows to limit changes in a single step

    !  delta_HR_limit = 5d-2
 
   !   time_delta_coeff = 0.8
     
     ! delta_lgTeff_limit = 0.005d0 !0.01d0       

     ! delta_lgL_limit = 0.1d0
      
      !max_q_for_conv_timescale = 0.999999999d0

      min_dq = 1d-16    ! gets multiplied by mesh delta coeff

      mesh_delta_coeff = 0.6 !0.5 !0.75 !0.5

      max_allowed_nz = 40000
      
     ! xa_clip_limit = 1d-12

      max_dq = 3d-3

      !tolerate larger jumps and residuals in solver tries

      !tolerate larger deviations in the chemical composition
      !min_xa_hard_limit = -1d-4 ! -1d-5

      ignore_species_in_max_correction = .true. ! .false.
      
      min_xa_hard_limit = -3d-3 ! -1d-5
     
      sum_xa_hard_limit = 5d-3 !7.5d-4 ! 5d-4
      !correction_xa_limit = 1d-5 !5d-3

      !tolerate larger jumps
      max_resid_jump_limit = 1d14
      max_corr_jump_limit = 1d12

      !tolerate larger residuals
      
    !  gold_tol_max_residual3 = 1d-3 !1d-4
    !  gold_tol_residual_norm3 = 1d-5 !5d-6 !1d-6     
      
      !correction_xa_limit = 1d-2 !5d-3     

      !force_timestep_min_years = 1d-4 

      report_bad_negative_xa = .true.
 
    !  max_q_for_conv_timescale = 0.999d0

      convergence_ignore_equL_residuals = .true. 

      !bad_max_correction_series_limit = 5 !2
      
     ! tol_bad_max_correction = 1d-3 ! 0d0
      
      !include_L_in_correction_limits = .false. !.true.
      include_v_in_correction_limits = .false !.true.

   
     ! tol_max_residual1 = 5d-9

      !tol_residual_norm3 = 1d95
      !tol_max_residual3 = 1d95
    
     !tolerate larger max corrections

    !  tol_max_correction = 3d-2 ! 3d-3 
    !  tol_correction_norm = 3d-4 !3d-5
  
    !  restore_mesh_on_retry = .true.

      max_years_for_timestep = 1d4

    !  max_allowed_nz = 20000


!solver matrix limits

     hydro_mtx_max_allowed_logT = 14d0  ! 12d0
   !  hydro_mtx_max_allowed_logRho = 14d0 ! 12d0

     hydro_mtx_min_allowed_logT = -1d1 ! 1d0
   ! hydro_mtx_min_allowed_logRho = -1d2 ! -1d2

!solver controls

    !bad_max_correction_series_limit = 3
    !include_v_in_correction_limits = .false.
    !include_L_in_correction_limits = .false.
    
    !solver_itermin = 12
    
     !solver_itermin_until_reduce_min_corr_coeff = 0

     trace_solver_damping = .true.
     do_solver_damping_for_neg_xa = .true. !.true.
  
     fix_eps_grav_transition_to_grid = .true.   
      
    ! max_tries_for_retry = 50 !25

  !   Tcenter_min_for_sig_min_factor_full_on = 1d7 !3.2d9

   !  Tcenter_max_for_sig_min_factor_full_off = 5d6 !2.8d9

   !  max_delta_m_to_bdy_for_sig_min_factor = 5d-4 !0.5d0

   !  delta_m_lower_for_sig_min_factor = 1d-99 !0.1d0    

    !  report_solver_progress = .true.  
      

    !  delta_mdot_atol = 5d-3 !1d-3
    !  delta_mdot_rtol = 1.0d0 !0.5d0


      num_trace_history_values = 2
      trace_history_value_name(1) = 'rel_E_err'
      trace_history_value_name(2) = 'log_rel_run_E_err'

      photo_interval = 100
      profile_interval = 250
      history_interval = 3
      terminal_interval = 1
      write_header_frequency = 10
      max_num_profile_models = 1000

/ ! end of controls namelist



&pgstar

         Grid8_win_flag = .true.
         Grid8_file_flag =.true.
         Grid8_win_width = 7
         Grid8_file_dir = 'png'
         Grid8_file_interval = 20
         Summary_Burn_xaxis_name = 'logxq' 
         Summary_Burn_xaxis_reversed = .true.
         Summary_Burn_xmin = -14 ! -101d0 ! only used if /= -101d0
         Summary_Burn_xmax = -1 ! -101d0 ! only used if /= -101d0
         Abundance_xaxis_name = 'logxq' 
         Abundance_xaxis_reversed = .true.
         Abundance_xmin = -14 ! -101d0 ! only used if /= -101d0
         Abundance_xmax = -1 ! -101d0 ! only used if /= -101d0

    Text_Summary1_name(5,2) = 'time_step_sec'
    Text_Summary1_name(6,2) = 'max_tau_conv'
    !Text_Summary1_name(7,2) = 'cz_bot_radius'
    !Text_Summary1_name(8,2) = 'cz_top_radius'

/ ! end of pgstar namelist