Bugfix/neu plasmon cooling2

docs/source/changelog.rst

@@ -44,6 +44,11 @@ The `MESA SDK <http://user.astro.wisc.edu/~townsend/static.php?ref=mesasdk>`__ r
 Bug Fixes
 ---------
 
+The plasma neutrino cooling rate was using an hardcoded prefactor calculated with a Weinberg angle of 0.2319, while all other neutrino cooling processes used calculated prefactors taking the Weinberg angle as input, with default value 0.22290.
+Thus, modifying the value of the Weinberg angle resulted in changes to neutrino cooling processes except for the plasma neutrinos.
+This was the case in all previous MESA versions, and was found and fixed by user Garv Chauhan, see `gh-998 https://github.com/MESAHub/mesa/pull/998`_. Plasma neutrinos now use the same Weinberg angle as all other processes and changing its value will affect the corresponding cooling rate.
+The change in default Weinberg angle results in small numerical differences for stars where plasma neutrino cooling is significant.
+
 Fixed a bug where RSP photo restarts did not immediately reconstruct ``s% L``,
 which could leave ``s% L(1)`` with an uninitialized near-zero value and crash
 MESA when the KH timescale was recalculated on restart.

neu/private/mod_neu.f90

@@ -35,6 +35,7 @@ module mod_neu
       real(dp), parameter :: cvp   = 1.0d0 - cv
       real(dp), parameter :: ca    = 0.5d0
       real(dp), parameter :: cap   = 1.0d0 - ca
+      real(dp), parameter :: tfac0 = cv*cv + (num_neu_fam-1.0d0) * (cvp*cvp)
       real(dp), parameter :: tfac1 = cv*cv + ca*ca + (num_neu_fam-1.0d0) * (cvp*cvp+cap*cap)
       real(dp), parameter :: tfac2 = cv*cv - ca*ca + (num_neu_fam-1.0d0) * (cvp*cvp - cap*cap)
       real(dp), parameter :: tfac3 = tfac2/tfac1
@@ -237,13 +238,13 @@ subroutine neutrinos(T, logT, Rho, logRho, abar, zbar, log10_Tlim,  &
          temp = T
       end if
 
-      if (T <= 0) then
+      if (temp <= 0) then
          info = -1
          return
       end if
 
       if (logT == arg_not_provided) then
-         logtemp = log10(T)
+         logtemp = log10(temp)
       else
          logtemp = logT
       end if
@@ -263,13 +264,13 @@ subroutine neutrinos(T, logT, Rho, logRho, abar, zbar, log10_Tlim,  &
          den = Rho
       end if
 
-      if (Rho <= 0) then
+      if (den <= 0) then
          info = -1
          return
       end if
 
       if (logRho == arg_not_provided) then
-         logden = log10(Rho)
+         logden = log10(den)
       else
          logden = logRho
       end if
@@ -1579,8 +1580,8 @@ subroutine plas_neu(splas,splasdt,splasdd,splasda,splasdz, input)
          splasdz = a2*splasdz + a3*gl2dz*a1
 
 
-         a2      = 0.93153d0 * 3.0d21 * input% xl9
-         a3      = 0.93153d0 * 3.0d21 * 9.0d0*input% xl8*input% xldt
+         a2      = tfac0 * 3.0d21 * input% xl9
+         a3      = tfac0 * 3.0d21 * 9.0d0*input% xl8*input% xldt
 
          a1      = splas
          splas   = a2*a1

star/test_suite/20M_z2m2_high_rotation/history_columns.list

@@ -138,7 +138,7 @@
 
       ! you might want to get a more complete list of mixing regions by using the following
 
-      !mixing_regions <integer> ! note: this includes regions where the mixing type is no_mixing.
+      mixing_regions 20 ! note: this includes regions where the mixing type is no_mixing.
 
          ! the <integer> is the number of regions to report
          ! there will be 2*<integer> columns for this in the log file, 2 for each region.
@@ -340,7 +340,7 @@
 
       ! you might want to get a more complete list of burning regions by using the following
 
-      !burning_regions <integer>
+      burning_regions 20
          ! the <integer> is the number of regions to report
          ! there will be 2*<integer> columns for this in the log file, 2 for each region.
          ! the first column for a region gives int(sign(val)*log10(max(1,abs(val))))

star/test_suite/20M_z2m2_high_rotation/inlist_to_end_core_he_burn

@@ -38,8 +38,7 @@
       set_initial_model_number = .true.
       initial_model_number = 0
 
-      pgstar_flag = .false.
-      !pgstar_flag = .true.
+      pgstar_flag = .true.
 
 / ! end of star_job namelist
 
@@ -219,6 +218,7 @@
       warn_rates_for_high_temp = .true.
       max_safe_logT_for_rates = 10.5d0
 
+      use_superad_reduction = .true. 
 
       Pextra_factor = 2 ! for atm
 

star/test_suite/ns_he/inlist_to_flash

@@ -22,7 +22,7 @@
       new_cumulative_energy_error = 0d0
 
       set_initial_dt = .true.
-      seconds_for_initial_dt = 1d-8
+      seconds_for_initial_dt = 1d-4
 
       set_initial_model_number = .true.
       initial_model_number = 0