diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index 61bec3b8e..9207e3b66 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -8205,9 +8205,12 @@ ! use_dPrad_dm_form_of_T_gradient_eqn ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! use_flux_limiting_with_dPrad_dm_form + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! These are for alternatives ways to determine the T gradient. ! The standard form of the equation is @@ -8231,6 +8234,11 @@ !| use L_rad = L*gradT/gradr (see, e.g., Cox&Giuli 14.109) !| else !| use L_rad = L + !| if (use_flux_limiting_with_dPrad_dm_form) + !| flxR = 4 * pi * r^2 * abs(dT^4/dm) / kap + !| flxLambda = (6 + 3*flxR) / (6 + 3*flxR + flxR^2) (see Levermore & Pomraning 1981) + !| L_rad = L_rad / flxLambda + !| ! With the resulting ``L_rad``, determine the expected dT/dm by @@ -8239,6 +8247,7 @@ ! :: use_dPrad_dm_form_of_T_gradient_eqn = .false. + use_flux_limiting_with_dPrad_dm_form = .false. use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = .false. diff --git a/star/defaults/profile_columns.list b/star/defaults/profile_columns.list index 644180091..bf0af4b30 100644 --- a/star/defaults/profile_columns.list +++ b/star/defaults/profile_columns.list @@ -148,6 +148,8 @@ !pgas ! gas pressure at center of zone (electrons and ions) !logPgas ! log10(pgas) !pgas_div_ptotal ! pgas/pressure + !flux_limit_lambda ! flux limiter defined as in Levermore & Pomraning 1981 + !flux_limit_R ! flux ratio defined as in Levermore & Pomraning 1981 !eta ! electron degeneracy parameter (eta >> 1 for significant degeneracy) !mu ! mean molecular weight per gas particle (ions + free electrons) @@ -454,6 +456,9 @@ !lum_conv_div_lum_Edd !lum_conv_div_L !lum_rad_div_L + !Frad_div_cUrad ! Frad/(C*Urad), must be < 1 to not violate c. + !flux_limit_lambda + !flux_limit_R !lum_rad_div_L_Edd_sub_fourPrad_div_PchiT ! density increases outward if this is > 0 ! see Joss, Salpeter, and Ostriker, "Critical Luminosity", ApJ 181:429-438, 1973. diff --git a/star/private/alloc.f90 b/star/private/alloc.f90 index d93bb8089..5910b4fe0 100644 --- a/star/private/alloc.f90 +++ b/star/private/alloc.f90 @@ -766,6 +766,11 @@ subroutine star_info_arrays(s, c_in, action_in, ierr) call do1(s% dr_div_csound, c% dr_div_csound) if (failed('dr_div_csound')) exit + call do1(s% flux_limit_R, c% flux_limit_R) + if (failed('flux_limit_R')) exit + call do1(s% flux_limit_lambda, c% flux_limit_lambda) + if (failed('flux_limit_lambda')) exit + call do1(s% ergs_error, c% ergs_error) if (failed('ergs_error')) exit call do1(s% gradr_factor, c% gradr_factor) diff --git a/star/private/ctrls_io.f90 b/star/private/ctrls_io.f90 index 2fdc823c2..604b8fa70 100644 --- a/star/private/ctrls_io.f90 +++ b/star/private/ctrls_io.f90 @@ -344,7 +344,7 @@ module ctrls_io include_composition_in_eps_grav, no_dedt_form_during_relax, & max_abs_rel_change_surf_lnS, & max_num_surf_revisions, Gamma_lnS_eps_grav_full_off, Gamma_lnS_eps_grav_full_on, & - use_dPrad_dm_form_of_T_gradient_eqn, use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn, dedt_eqn_r_scale, & + use_dPrad_dm_form_of_T_gradient_eqn, use_flux_limiting_with_dPrad_dm_form, use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn, dedt_eqn_r_scale, & RTI_A, RTI_B, RTI_C, RTI_D, RTI_max_alpha, RTI_C_X_factor, RTI_C_X0_frac, steps_before_use_velocity_time_centering, & RTI_dm_for_center_eta_nondecreasing, RTI_min_dm_behind_shock_for_full_on, RTI_energy_floor, & RTI_D_mix_floor, RTI_min_m_for_D_mix_floor, RTI_log_max_boost, RTI_m_full_boost, RTI_m_no_boost, & @@ -1843,6 +1843,7 @@ subroutine store_controls(s, ierr) s% Gamma_lnS_eps_grav_full_on = Gamma_lnS_eps_grav_full_on s% use_dPrad_dm_form_of_T_gradient_eqn = use_dPrad_dm_form_of_T_gradient_eqn + s% use_flux_limiting_with_dPrad_dm_form = use_flux_limiting_with_dPrad_dm_form s% use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn s% include_P_in_velocity_time_centering = include_P_in_velocity_time_centering s% include_L_in_velocity_time_centering = include_L_in_velocity_time_centering @@ -3524,6 +3525,7 @@ subroutine set_controls_for_writing(s, ierr) Gamma_lnS_eps_grav_full_on = s% Gamma_lnS_eps_grav_full_on use_dPrad_dm_form_of_T_gradient_eqn = s% use_dPrad_dm_form_of_T_gradient_eqn + use_flux_limiting_with_dPrad_dm_form = s% use_flux_limiting_with_dPrad_dm_form use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = s% use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn steps_before_use_velocity_time_centering = s% steps_before_use_velocity_time_centering include_P_in_velocity_time_centering = s% include_P_in_velocity_time_centering diff --git a/star/private/hydro_temperature.f90 b/star/private/hydro_temperature.f90 index d61eb625c..07530fef9 100644 --- a/star/private/hydro_temperature.f90 +++ b/star/private/hydro_temperature.f90 @@ -58,6 +58,7 @@ subroutine do1_alt_dlnT_dm_eqn(s, k, nvar, ierr) type(auto_diff_real_star_order1) :: L_ad, r_00, area, area2, Lrad_ad, & kap_00, kap_m1, kap_face, d_P_rad_expected_ad, T_m1, T4_m1, T_00, T4_00, & P_rad_m1, P_rad_00, d_P_rad_actual_ad, resid + type(auto_diff_real_star_order1) :: flxR, flxLambda integer :: i_equL logical :: dbg @@ -109,10 +110,29 @@ subroutine do1_alt_dlnT_dm_eqn(s, k, nvar, ierr) !d_P_rad_expected = d_P_rad_expected*s% gradr_factor(k) !TODO(Pablo): check this - P_rad_m1 = (crad/3d0)*T4_m1 - P_rad_00 = (crad/3d0)*T4_00 + P_rad_m1 = (crad/3._dp)*T4_m1 + P_rad_00 = (crad/3._dp)*T4_00 d_P_rad_actual_ad = P_rad_m1 - P_rad_00 + ! enable flux-limited radiation transport derived by Levermore & Pomraning 1981 + s% flux_limit_R(k) = 0._dp + s% flux_limit_lambda(k) =0._dp + if (s% use_flux_limiting_with_dPrad_dm_form) then + ! calculate the flux ratio R + flxR = area * abs(T4_m1 - T4_00) / dm_bar / & + (kap_face * 0.5_dp * (T4_m1 + T4_00)) + + s% flux_limit_R(k) = flxR%val + + ! calculate the flux limiter lambda + flxLambda = (6._dp + 3._dp*flxR) / (6._dp + (3._dp + flxR)*flxR) + + s% flux_limit_lambda(k) = flxLambda%val + + ! calculate d_P_rad given the flux limiter + d_P_rad_expected_ad = d_P_rad_expected_ad / flxLambda + end if + ! residual resid = (d_P_rad_expected_ad - d_P_rad_actual_ad)/scale s% equ(i_equL, k) = resid%val diff --git a/star/private/profile_getval.f90 b/star/private/profile_getval.f90 index b37dac06f..e9adb7ac0 100644 --- a/star/private/profile_getval.f90 +++ b/star/private/profile_getval.f90 @@ -427,6 +427,9 @@ subroutine getval_for_profile(s, c, k, val, int_flag, int_val) case (p_lum_conv_MLT) val = s% L_conv(k)/Lsun + case(p_Frad_div_cUrad) + val = ((s% L(k) - s% L_conv(k)) / (4._dp*pi*pow2(s%r(k)))) & + /(clight * s% Prad(k) *3._dp) case (p_lum_rad_div_L_Edd_sub_fourPrad_div_PchiT) val = get_Lrad_div_Ledd(s,k) - 4*s% Prad(k)/(s% Peos(k)*s% chiT(k)) case (p_lum_rad_div_L_Edd) @@ -774,7 +777,15 @@ subroutine getval_for_profile(s, c, k, val, int_flag, int_val) val = (s% Prad(k)/s% Pgas(k))/max(1d-12,s% L(k)/get_Ledd(s,k)) case (p_pgas_div_p) val = s% Pgas(k)/s% Peos(k) - + case (p_flux_limit_R) + if (s% use_dPrad_dm_form_of_T_gradient_eqn .and. k > 1) & + val = s% flux_limit_R(k) + case (p_flux_limit_lambda) + if (s% use_dPrad_dm_form_of_T_gradient_eqn .and. k > 1) then + val = s% flux_limit_lambda(k) + else + val = 1d0 + end if case (p_cell_collapse_time) if (s% v_flag) then if (k == s% nz) then diff --git a/star/private/read_model.f90 b/star/private/read_model.f90 index ecc70fdfb..929b5e052 100644 --- a/star/private/read_model.f90 +++ b/star/private/read_model.f90 @@ -150,6 +150,9 @@ subroutine finish_load_model(s, restart, ierr) call fill_ad_with_zeros(s% P_face_ad,1,-1) end if + s% flux_limit_R(1:nz) = 0 + s% flux_limit_lambda(1:nz) = 0 + if (s% RSP_flag) then call RSP_setup_part1(s, restart, ierr) if (ierr /= 0) then diff --git a/star/private/star_profile_def.f90 b/star/private/star_profile_def.f90 index 22d53bdcf..5d399df02 100644 --- a/star/private/star_profile_def.f90 +++ b/star/private/star_profile_def.f90 @@ -708,9 +708,13 @@ module star_profile_def integer, parameter :: p_RTI_du_diffusion_kick = p_dPdr_dRhodr_info + 1 integer, parameter :: p_log_du_kick_div_du = p_RTI_du_diffusion_kick + 1 - integer, parameter :: p_lum_rad_div_L_Edd_sub_fourPrad_div_PchiT = p_log_du_kick_div_du + 1 + integer, parameter :: p_Frad_div_cUrad = p_log_du_kick_div_du + 1 + integer, parameter :: p_lum_rad_div_L_Edd_sub_fourPrad_div_PchiT = p_Frad_div_cUrad + 1 - integer, parameter :: p_col_id_max = p_lum_rad_div_L_Edd_sub_fourPrad_div_PchiT + integer, parameter :: p_flux_limit_R = p_lum_rad_div_L_Edd_sub_fourPrad_div_PchiT + 1 + integer, parameter :: p_flux_limit_lambda = p_flux_limit_R + 1 + + integer, parameter :: p_col_id_max = p_flux_limit_lambda character (len=maxlen_profile_column_name) :: profile_column_name(p_col_id_max) type (integer_dict), pointer :: profile_column_names_dict @@ -1395,8 +1399,12 @@ subroutine profile_column_names_init(ierr) profile_column_name(p_tau_epsnuc) = 'tau_epsnuc' profile_column_name(p_tau_cool) = 'tau_cool' + profile_column_name(p_Frad_div_cUrad) = 'Frad_div_cUrad' profile_column_name(p_lum_rad_div_L_Edd_sub_fourPrad_div_PchiT) = 'lum_rad_div_L_Edd_sub_fourPrad_div_PchiT' + profile_column_name(p_flux_limit_R) = 'flux_limit_R' + profile_column_name(p_flux_limit_lambda) = 'flux_limit_lambda' + cnt = 0 do i=1,p_col_id_max if (len_trim(profile_column_name(i)) == 0) then diff --git a/star/test_suite/ppisn/profile_columns.list b/star/test_suite/ppisn/profile_columns.list index dfa0c3268..dfb83f028 100644 --- a/star/test_suite/ppisn/profile_columns.list +++ b/star/test_suite/ppisn/profile_columns.list @@ -456,7 +456,10 @@ !lum_rad_div_L !lum_rad_div_L_Edd_sub_fourPrad_div_PchiT ! density increases outward if this is > 0 ! see Joss, Salpeter, and Ostriker, "Critical Luminosity", ApJ 181:429-438, 1973. - + Frad_div_cUrad + flux_limit_lambda + flux_limit_R + gradT ! mlt value for required temperature gradient dlnT/dlnP gradr ! dlnT/dlnP required for purely radiative transport diff --git a/star_data/private/star_controls.inc b/star_data/private/star_controls.inc index a5c245bd2..453da3472 100644 --- a/star_data/private/star_controls.inc +++ b/star_data/private/star_controls.inc @@ -825,6 +825,7 @@ logical :: use_gravity_rotation_correction, & use_dPrad_dm_form_of_T_gradient_eqn, & + use_flux_limiting_with_dPrad_dm_form, & use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn, & no_dedt_form_during_relax, & include_P_in_velocity_time_centering, & diff --git a/star_data/public/star_data_step_work.inc b/star_data/public/star_data_step_work.inc index 572518933..fbbd05633 100644 --- a/star_data/public/star_data_step_work.inc +++ b/star_data/public/star_data_step_work.inc @@ -102,6 +102,10 @@ real(dp), pointer :: chiT(:) ! dlnPeos_dlnT at constant Rho real(dp), pointer :: QQ(:) ! thermal expansion coefficient + ! flux-limited radiation transport + real(dp), pointer :: flux_limit_lambda(:) ! flux limiter defined as in Levermore & Pomraning 1981 + real(dp), pointer :: flux_limit_R(:) ! flux ratio defined as in Levermore & Pomraning 1981 + ! eos blend information real(dp), pointer :: eos_frac_OPAL_SCVH(:) real(dp), pointer :: eos_frac_HELM(:)