Feature/reintroduce flux limiter

star/defaults/controls.defaults

@@ -8228,9 +8228,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
 
@@ -8254,6 +8257,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
 
@@ -8262,6 +8270,7 @@
       ! ::
 
     use_dPrad_dm_form_of_T_gradient_eqn = .false.
+    use_flux_limiting_with_dPrad_dm_form = .true.
     use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = .false.    
 
 

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)

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)

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, &
@@ -1847,6 +1847,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
@@ -3531,6 +3532,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

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
@@ -113,6 +114,24 @@ subroutine do1_alt_dlnT_dm_eqn(s, k, nvar, ierr)
          P_rad_00 = (crad/3d0)*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) = 0d0
+         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.5d0 * (T4_m1 + T4_00))
+
+            s% flux_limit_R(k) = flxR%val
+
+            ! calculate the flux limiter lambda
+            flxLambda = (6d0 + 3d0*flxR) / (6d0 + (3d0 + 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

star/private/profile_getval.f90

@@ -767,7 +767,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

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

star/private/star_profile_def.f90

@@ -710,7 +710,10 @@ module star_profile_def
 
       integer, parameter :: p_lum_rad_div_L_Edd_sub_fourPrad_div_PchiT = p_log_du_kick_div_du + 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
@@ -1397,6 +1400,9 @@ subroutine profile_column_names_init(ierr)
 
          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

star_data/private/star_controls.inc

@@ -824,6 +824,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, &

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(:)