From 556af5babe442c5393e1746f33e94fcfaf707c0b Mon Sep 17 00:00:00 2001 From: Ebraheem Farag <63124736+Debraheem@users.noreply.github.com> Date: Wed, 26 Aug 2026 20:55:02 -0400 Subject: [PATCH] 6: Relaxation, mass adjustment, and shell removal --- docs/source/changelog.rst | 15 +++++ star/defaults/star_job.defaults | 29 +++++++++ star/private/adjust_mass.f90 | 85 ++++++++++++++++++++++++- star/private/eps_mdot.f90 | 16 +++-- star/private/evolve.f90 | 6 +- star/private/relax.f90 | 66 +++++++++++++++---- star/private/remove_shells.f90 | 8 ++- star/private/star_job_ctrls_io.f90 | 12 ++++ star_data/private/star_job_controls.inc | 4 ++ 9 files changed, 217 insertions(+), 24 deletions(-) diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 6fc9d9e9ab..74d04d18fb 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -78,6 +78,10 @@ resolved shocks using the default-off ``Riemann_shock_D_mix_reduction_on`` and pressure-jump, shock-strength, and diffusion-factor profile columns are available for diagnostics. +Entropy-profile relaxation can now use an implicit or normalized energy +source and can temporarily select the Newton correction scale and retry hold. +The new controls restore their incoming solver settings when relaxation ends. + .. _Bug Fixes main: Bug Fixes @@ -161,6 +165,17 @@ Added an optional floor on the atmospheric pressure used by the momentum outer boundary. The floor prevents a hydrostatic atmosphere from supplying less than the radiation pressure at its boundary temperature. +Fixed velocity remapping during mass loss. Momentum is now conservatively +remapped on the cell or dual-cell mass coordinates, and unresolved kinetic +energy is returned locally through ``eps_mdot`` instead of being lost. + +Fixed nonpositive surface optical depths after removing surface cells from a +dynamic model. Surface removal now retains the existing optical depth when +the hydrostatic pressure estimate is nonphysical. Split/merge metric zoning +also disables its logarithmic optical-depth component when the optical-depth +coordinate is not positive, preventing invalid logarithms from driving +unbounded mesh refinement. + Important bug fix for ``r26.4.1`` identified by Emily Sandford and Louis Siebenaler: the ``lowT_Freedman11`` opacity option used ``[M/H]`` labels as the metal mass fraction when interpolating in ``Z``, resulting in incorrect opacities. We recommend users who use these low-temperature opacities, such as in planet models, update to the latest MESA version or employ the fixes in :ref:`the known bugs entry ` and `gh-993 `_. The plasmon neutrino cooling rate used a 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 plasmon neutrinos. This affects all previous MESA versions, and was found and fixed by user Garv Chauhan, see :ref:`the known bugs entry ` and `gh-998 `_. Plasmon neutrinos now use the same Weinberg angle as all other processes and changing its value will affect the corresponding cooling rate. Changes to the plasmon neutrino prefactor for MESA's default Weinberg angle result in small numerical differences for stars where plasmon neutrino cooling is significant. diff --git a/star/defaults/star_job.defaults b/star/defaults/star_job.defaults index 354e99893a..068e3abe72 100644 --- a/star/defaults/star_job.defaults +++ b/star/defaults/star_job.defaults @@ -1098,6 +1098,14 @@ ! ~~~~~~~~~~~~~~~~~~~~~~~~ ! num_timescales_for_relax_entropy ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! relax_entropy_use_implicit_source + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! relax_entropy_use_normalized_source + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! relax_entropy_scale_max_correction + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! relax_entropy_retry_hold + ! ~~~~~~~~~~~~~~~~~~~~~~~~ ! relax_entropy_filename ! ~~~~~~~~~~~~~~~~~~~~~~ ! get_entropy_for_relax_from_eos @@ -1116,6 +1124,23 @@ ! using ``min_T_for_acceleration_limited_conv_velocity = 0`` (see controls.defaults), and the timescale ! for relaxation should be very short (less than a second). + ! If ``relax_entropy_use_implicit_source`` is true, the local density and temperature derivatives of the + ! entropy relaxation heating are included in the structure Jacobian. The heating and target profile are + ! unchanged. + + ! If ``relax_entropy_use_normalized_source`` is true, use + + ! :: + + ! s% extra_heat(k) = T(k)*(desired_entropy(k) - entropy(k))/(timescale_for_relax_entropy*secyer) + + ! so that ``timescale_for_relax_entropy`` is the local entropy e-folding time in the absence of other + ! energy sources. This source can be used with either the explicit or implicit energy hook. + + ! A positive ``relax_entropy_scale_max_correction`` temporarily replaces ``scale_max_correction`` during + ! entropy relaxation. A nonnegative ``relax_entropy_retry_hold`` similarly replaces ``retry_hold``. Both + ! controls are restored after relaxation. + ! ``relax_entropy_filename`` holds the desired entropy profile information ! file format for relax entropy @@ -1165,6 +1190,10 @@ timescale_for_relax_entropy = 1d-9 max_dt_for_relax_entropy = 1d-9 num_timescales_for_relax_entropy = 100 + relax_entropy_use_implicit_source = .false. + relax_entropy_use_normalized_source = .false. + relax_entropy_scale_max_correction = 0d0 + relax_entropy_retry_hold = -1 relax_entropy_filename = '' get_entropy_for_relax_from_eos = '' diff --git a/star/private/adjust_mass.f90 b/star/private/adjust_mass.f90 index 1bf13de859..f3c1f6f05b 100644 --- a/star/private/adjust_mass.f90 +++ b/star/private/adjust_mass.f90 @@ -166,7 +166,7 @@ subroutine update_radius(s) end subroutine update_radius - subroutine do_adjust_mass(s, species, ierr) + subroutine do_adjust_mass(s, species, velocity_remap_kinetic_energy, ierr) use adjust_xyz, only: get_xa_for_accretion use star_utils, only: report_xa_bad_nums, & start_time, update_time @@ -175,6 +175,7 @@ subroutine do_adjust_mass(s, species, ierr) type (star_info), pointer :: s integer, intent(in) :: species + real(dp), intent(out) :: velocity_remap_kinetic_energy(:) integer, intent(out) :: ierr real(dp) :: & @@ -186,7 +187,7 @@ subroutine do_adjust_mass(s, species, ierr) real(dp), dimension(:), allocatable :: & rxm_old, rxm_new, old_cell_mass, new_cell_mass, & oldloc, newloc, oldval, newval, xm_old, xm_new, & - xq_center_old, xq_center_new + xq_center_old, xq_center_new, velocity_old real(dp), dimension(:,:), allocatable :: xa_old real(dp), pointer :: work(:) @@ -201,6 +202,7 @@ subroutine do_adjust_mass(s, species, ierr) if (dbg) write(*,*) 'do_adjust_mass' + velocity_remap_kinetic_energy(:) = 0d0 call save_for_eps_mdot(s) ierr = 0 @@ -281,6 +283,11 @@ subroutine do_adjust_mass(s, species, ierr) do k=1,nz old_cell_mass(k) = old_xmstar*s% dq(k) + if (s% u_flag) then + velocity_old(k) = s% xh(s% i_u,k) + else if (s% v_flag) then + velocity_old(k) = s% xh(s% i_v,k) + end if end do xm_old(1) = 0 xq_center_old(1) = 0.5d0*s% dq(1) @@ -418,6 +425,12 @@ subroutine do_adjust_mass(s, species, ierr) end if end do + if ((s% u_flag .or. s% v_flag) .and. delta_m < 0d0) then + call remap_velocity_for_mass_loss( & + s, nz, k_const_mass, velocity_old, old_cell_mass, new_cell_mass, & + velocity_remap_kinetic_energy) + end if + call set_xa(s, nz, k_const_mass, species, xa_old, xaccrete, & rxm_old, rxm_new, mmax, old_cell_mass, new_cell_mass, ierr) if (ierr /= 0) then @@ -524,6 +537,72 @@ subroutine do_adjust_mass(s, species, ierr) contains + subroutine remap_velocity_for_mass_loss( & + s, nz, k_const_mass, velocity_old, old_cell_mass, new_cell_mass, & + remap_kinetic_energy) + use mass_utils, only: integrate_conserved + type (star_info), pointer :: s + integer, intent(in) :: nz, k_const_mass + real(dp), intent(in) :: velocity_old(:), old_cell_mass(:), new_cell_mass(:) + real(dp), intent(out) :: remap_kinetic_energy(:) + + integer :: k + real(dp) :: vals_outside(2), delta_specific_kinetic_energy + real(dp), allocatable :: vals_old(:,:), vals_new(:,:) + real(qp), allocatable :: dm_old(:), dm_new(:) + + allocate(vals_old(nz,2), vals_new(nz,2), dm_old(nz), dm_new(nz)) + + do k=1,nz + if (s% u_flag) then + dm_old(k) = old_cell_mass(k) + dm_new(k) = new_cell_mass(k) + else if (k == 1) then + ! v(1) has only the outer half of cell 1 in its dual cell. + dm_old(k) = 0.5d0*old_cell_mass(k) + dm_new(k) = 0.5d0*new_cell_mass(k) + else + ! Other v points have half of each adjacent cell. + dm_old(k) = 0.5d0*(old_cell_mass(k-1) + old_cell_mass(k)) + dm_new(k) = 0.5d0*(new_cell_mass(k-1) + new_cell_mass(k)) + end if + vals_old(k,1) = velocity_old(k) + vals_old(k,2) = 0.5d0*pow2(velocity_old(k)) + end do + vals_outside(:) = 0d0 + + ! Remap specific momentum and kinetic energy on the same mass overlaps. + call integrate_conserved( & + vals_new, vals_old, vals_outside, dm_new, dm_old, nz, 2) + + remap_kinetic_energy(:) = 0d0 + do k=1,nz + if (k > k_const_mass) then + vals_new(k,1) = velocity_old(k) + delta_specific_kinetic_energy = 0d0 + else + ! Momentum projection removes only unresolved velocity variance. + delta_specific_kinetic_energy = max(0d0, & + vals_new(k,2) - 0.5d0*pow2(vals_new(k,1))) + end if + if (s% u_flag) then + s% u(k) = vals_new(k,1) + s% xh(s% i_u,k) = s% u(k) + remap_kinetic_energy(k) = & + new_cell_mass(k)*delta_specific_kinetic_energy + else + s% v(k) = vals_new(k,1) + s% xh(s% i_v,k) = s% v(k) + remap_kinetic_energy(k) = remap_kinetic_energy(k) + & + 0.5d0*new_cell_mass(k)*delta_specific_kinetic_energy + if (k > 1) remap_kinetic_energy(k-1) = remap_kinetic_energy(k-1) + & + 0.5d0*new_cell_mass(k-1)*delta_specific_kinetic_energy + end if + end do + + end subroutine remap_velocity_for_mass_loss + + real(dp) function angular_momentum_removed(ierr) result(J) ! when call this, s% j_rot is still for old mass integer, intent(out) :: ierr @@ -580,7 +659,7 @@ subroutine do_alloc(ierr) integer, intent(out) :: ierr allocate(rxm_old(nz), rxm_new(nz), old_cell_mass(nz), new_cell_mass(nz), & xa_old(species,nz), oldloc(nz), newloc(nz), oldval(nz), newval(nz), & - xm_old(nz), xm_new(nz), xq_center_old(nz), xq_center_new(nz)) + xm_old(nz), xm_new(nz), xq_center_old(nz), xq_center_new(nz), velocity_old(nz)) call do_work_arrays(.true.,ierr) end subroutine do_alloc diff --git a/star/private/eps_mdot.f90 b/star/private/eps_mdot.f90 index 40166c40d9..08cfbf16ca 100644 --- a/star/private/eps_mdot.f90 +++ b/star/private/eps_mdot.f90 @@ -469,13 +469,14 @@ subroutine leak(nz, i_start, i_end, i_min, i_max, j_min, j_max, pf,& end subroutine leak - subroutine calculate_eps_mdot(s, dt, ierr) + subroutine calculate_eps_mdot(s, dt, velocity_remap_kinetic_energy, ierr) use adjust_mass, only: compute_prev_mesh_dm ! Inputs type (star_info), pointer :: s - real(dp) :: dt - integer :: ierr + real(dp), intent(in) :: dt + real(dp), intent(in) :: velocity_remap_kinetic_energy(:) + integer, intent(out) :: ierr ! Intermediates logical, parameter :: dbg = .false. @@ -491,6 +492,7 @@ subroutine calculate_eps_mdot(s, dt, ierr) integer, dimension(:,:), allocatable :: ranges real(qp), dimension(:), allocatable :: mesh_intersects + ierr = 0 if (s% mstar_dot == 0d0 .or. dt <= 0d0) then s% eps_mdot(1:s%nz) = 0d0 s% mdot_adiabatic_surface = 0d0 @@ -634,6 +636,10 @@ subroutine calculate_eps_mdot(s, dt, ierr) change_sum = change_sum + sum%value() / (dt) + ! Velocity averaging converts unresolved kinetic energy into heat. + ! Keep that local source out of the thermal leakage calculation. + sum = sum - real(velocity_remap_kinetic_energy(j),qp) + ! Multiplicative factors eps_mdot_per_total_mass(j) = sum % value() / (s%dm(j) * dt) / total_mass_through_cell(j) @@ -642,6 +648,7 @@ subroutine calculate_eps_mdot(s, dt, ierr) err = 0d0 do j=1,nz err = err + eps_mdot_per_total_mass(j) * s%dm(j) * dt * total_mass_through_cell(j) + err = err + velocity_remap_kinetic_energy(j) end do err = err - s%mdot_acoustic_surface err = err - te_bar(1) * delta_m @@ -657,7 +664,8 @@ subroutine calculate_eps_mdot(s, dt, ierr) total_mass_through_cell, eps_mdot_per_total_mass,& accumulated, mdot_adiabatic_surface, leak_frac) do j=1,nz - s%eps_mdot(j) = accumulated(j) + s%eps_mdot(j) = accumulated(j) + & + velocity_remap_kinetic_energy(j)/(s%dm(j)*dt) end do s% mdot_adiabatic_surface = -mdot_adiabatic_surface diff --git a/star/private/evolve.f90 b/star/private/evolve.f90 index fe37b19043..ab22e1d423 100644 --- a/star/private/evolve.f90 +++ b/star/private/evolve.f90 @@ -533,6 +533,7 @@ integer function do_step_part2(id, first_try) explicit_mdot, max_wind_mdot, wind_mdot, r_phot, kh_timescale, dmskhf, dmsfac, & too_large_wind_mdot, too_small_wind_mdot, boost, mstar_dot_nxt, & surf_omega_div_omega_crit_limit, dt + real(dp), allocatable :: velocity_remap_kinetic_energy(:) integer :: ph_k, mdot_action real(dp) :: implicit_mdot, ph_L, iwind_tolerance, iwind_lambda @@ -559,6 +560,7 @@ integer function do_step_part2(id, first_try) clock_rate = s% system_clock_rate trace = s% trace_evolve nz = s% nz + allocate(velocity_remap_kinetic_energy(nz)) call setup_for_implicit_mdot_loop @@ -588,13 +590,13 @@ integer function do_step_part2(id, first_try) else - call do_adjust_mass(s, s% species, ierr) + call do_adjust_mass(s, s% species, velocity_remap_kinetic_energy, ierr) if (failed('do_adjust_mass')) return s% star_mdot = s% mstar_dot/(Msun/secyer) ! dm/dt in msolar per year call set_vars_if_needed(s, dt, 'after do_adjust_mass', ierr) if (failed('set_vars_if_needed after do_adjust_mass')) return - call calculate_eps_mdot(s, dt, ierr) + call calculate_eps_mdot(s, dt, velocity_remap_kinetic_energy, ierr) if (failed('calculate_eps_mdot')) return if (s% mstar_dot /= 0d0) then diff --git a/star/private/relax.f90 b/star/private/relax.f90 index fc925dd094..98abdd7dda 100644 --- a/star/private/relax.f90 +++ b/star/private/relax.f90 @@ -56,9 +56,9 @@ module relax real(dp), parameter :: min_dlnz = -12 real(dp), parameter :: min_z = 1d-12 - ! some relax routines depend on things such as other_energy and other_torque - ! to which interpolation parameters cannot be passed directly. So for simplicity - ! use these two global variables instead. + ! some relax routines depend on energy and torque hooks to which interpolation + ! parameters cannot be passed directly. So for simplicity use these two global + ! variables instead. integer :: relax_num_pts real(dp), pointer :: relax_work_array(:) @@ -702,16 +702,19 @@ subroutine do_relax_entropy( & integer, intent(out) :: ierr integer, parameter :: lipar=2 - integer :: lrpar, max_model_number + integer :: lrpar, max_model_number, retry_hold real(dp), pointer :: rpar(:) - real(dp) :: starting_dt_next, mix_factor, dxdt_nuc_factor, max_years_for_timestep, time - logical :: do_element_diffusion, use_other_energy + real(dp) :: starting_dt_next, mix_factor, dxdt_nuc_factor, max_years_for_timestep, & + scale_max_correction, time + logical :: do_element_diffusion, use_other_energy, use_other_energy_implicit type (star_info), pointer :: s real(dp), pointer :: x(:), f1(:), f(:,:) integer, target :: ipar_ary(lipar) integer, pointer :: ipar(:) procedure (other_energy_interface), pointer :: & other_energy => null() + procedure (other_energy_implicit_interface), pointer :: & + other_energy_implicit => null() ipar => ipar_ary @@ -755,9 +758,24 @@ subroutine do_relax_entropy( & s% max_years_for_timestep = s% job% max_dt_for_relax_entropy s% dt_next = min(s% dt_next, s% job% max_dt_for_relax_entropy * secyer) use_other_energy = s% use_other_energy - s% use_other_energy = .true. other_energy => s% other_energy - s% other_energy => entropy_relax_other_energy + use_other_energy_implicit = s% use_other_energy_implicit + other_energy_implicit => s% other_energy_implicit + if (s% job% relax_entropy_use_implicit_source) then + s% use_other_energy = .false. + s% use_other_energy_implicit = .true. + s% other_energy_implicit => entropy_relax_other_energy + else + s% use_other_energy = .true. + s% other_energy => entropy_relax_other_energy + s% use_other_energy_implicit = .false. + end if + scale_max_correction = s% scale_max_correction + if (s% job% relax_entropy_scale_max_correction > 0d0) & + s% scale_max_correction = s% job% relax_entropy_scale_max_correction + retry_hold = s% retry_hold + if (s% job% relax_entropy_retry_hold >= 0) & + s% retry_hold = s% job% relax_entropy_retry_hold time = s% time s% time = 0d0 @@ -774,6 +792,10 @@ subroutine do_relax_entropy( & s% max_years_for_timestep = max_years_for_timestep s% use_other_energy = use_other_energy s% other_energy => other_energy + s% use_other_energy_implicit = use_other_energy_implicit + s% other_energy_implicit => other_energy_implicit + s% scale_max_correction = scale_max_correction + s% retry_hold = retry_hold s% time = time call error_check('relax entropy',ierr) @@ -922,25 +944,43 @@ subroutine entropy_relax_other_energy(id, ierr) integer, intent(out) :: ierr type (star_info), pointer :: s integer :: k, nz, num_pts - real(dp), pointer :: vals(:), xq(:), x(:), f(:) + real(dp), allocatable :: vals(:), xq(:) + real(dp), pointer :: x(:), f(:) ierr = 0 call star_ptr(id, s, ierr) + if (ierr /= 0) return nz = s% nz num_pts = relax_num_pts allocate(vals(nz), xq(nz), stat=ierr) + if (ierr /= 0) return f(1:4*num_pts) => relax_work_array(num_pts+1:5*num_pts) x(1:num_pts) => relax_work_array(1:num_pts) xq(1) = s% dq(1)/2 ! xq for cell center do k = 2, nz xq(k) = xq(k-1) + (s% dq(k) + s% dq(k-1))/2 end do - call interp_values(x, num_pts, f, nz, xq, vals(:), ierr) - if (ierr /= 0) return + call interp_values(x, num_pts, f, nz, xq, vals, ierr) + if (ierr /= 0) then + deallocate(vals, xq) + return + end if do k = 1, s% nz - s% extra_heat(k) = ( 1d0 - exp(s%lnS(k))/vals(k) ) * exp(s%lnE(k)) - s% extra_heat(k) = s% extra_heat(k) / (s% job% timescale_for_relax_entropy * secyer) + if (s% job% relax_entropy_use_normalized_source) then + if (s% job% relax_entropy_use_implicit_source) then + s% extra_heat(k) = wrap_T_00(s,k)*(vals(k) - wrap_s_00(s,k)) + else + s% extra_heat(k) = s% T(k)*(vals(k) - exp(s% lnS(k))) + end if + else if (s% job% relax_entropy_use_implicit_source) then + s% extra_heat(k) = (1d0 - wrap_s_00(s,k)/vals(k))*wrap_e_00(s,k) + else + s% extra_heat(k) = (1d0 - exp(s% lnS(k))/vals(k))*exp(s% lnE(k)) + end if + s% extra_heat(k) = s% extra_heat(k) / & + (s% job% timescale_for_relax_entropy * secyer) end do + deallocate(vals, xq) end subroutine entropy_relax_other_energy subroutine do_relax_angular_momentum( & diff --git a/star/private/remove_shells.f90 b/star/private/remove_shells.f90 index d498a5e36a..e52b34404b 100644 --- a/star/private/remove_shells.f90 +++ b/star/private/remove_shells.f90 @@ -1041,7 +1041,10 @@ subroutine do_remove_surface(id, surface_k, ierr) if (skip < 1 .or. skip >= nz_old) return tau_surf_new = tau_eff(s,1+skip) - tau_factor_new = s% tau_factor*tau_surf_new/s% tau(1) + tau_factor_new = s% tau_factor + ! Keep the existing optical depth when no positive HSE atmosphere exists. + if (tau_surf_new > 0d0 .and. .not. is_bad(tau_surf_new)) & + tau_factor_new = tau_surf_new/s% tau_base if (dbg) write(*,1) 'tau_surf_old', s% tau(1) if (dbg) write(*,1) 'tau_factor_old', s% tau_factor @@ -1122,7 +1125,8 @@ subroutine do_remove_surface(id, surface_k, ierr) if (dbg) write(*,2) 's% dm(nz)/Msun', nz, s% dm(nz)/Msun if (dbg) write(*,2) 's% m(nz)/msun', nz, s% m(nz)/Msun - if (s% use_momentum_outer_bc) then + if (s% use_momentum_outer_bc .and. & + tau_factor_new > 0d0 .and. .not. is_bad(tau_factor_new)) then s% tau_factor = tau_factor_new s% force_tau_factor = s% tau_factor end if diff --git a/star/private/star_job_ctrls_io.f90 b/star/private/star_job_ctrls_io.f90 index a95c209178..6521753515 100644 --- a/star/private/star_job_ctrls_io.f90 +++ b/star/private/star_job_ctrls_io.f90 @@ -382,6 +382,10 @@ module star_job_ctrls_io timescale_for_relax_entropy, & max_dt_for_relax_entropy, & num_timescales_for_relax_entropy, & + relax_entropy_use_implicit_source, & + relax_entropy_use_normalized_source, & + relax_entropy_scale_max_correction, & + relax_entropy_retry_hold, & relax_entropy_filename, & get_entropy_for_relax_from_eos, & report_cell_for_xm, & @@ -916,6 +920,10 @@ subroutine store_star_job_controls(s) s% job% timescale_for_relax_entropy = timescale_for_relax_entropy s% job% max_dt_for_relax_entropy = max_dt_for_relax_entropy s% job% num_timescales_for_relax_entropy = num_timescales_for_relax_entropy + s% job% relax_entropy_use_implicit_source = relax_entropy_use_implicit_source + s% job% relax_entropy_use_normalized_source = relax_entropy_use_normalized_source + s% job% relax_entropy_scale_max_correction = relax_entropy_scale_max_correction + s% job% relax_entropy_retry_hold = relax_entropy_retry_hold s% job% relax_entropy_filename = relax_entropy_filename s% job% get_entropy_for_relax_from_eos = get_entropy_for_relax_from_eos s% job% report_cell_for_xm = report_cell_for_xm @@ -1454,6 +1462,10 @@ subroutine set_star_job_controls_for_writing(s, ierr) timescale_for_relax_entropy = s% job% timescale_for_relax_entropy max_dt_for_relax_entropy = s% job% max_dt_for_relax_entropy num_timescales_for_relax_entropy = s% job% num_timescales_for_relax_entropy + relax_entropy_use_implicit_source = s% job% relax_entropy_use_implicit_source + relax_entropy_use_normalized_source = s% job% relax_entropy_use_normalized_source + relax_entropy_scale_max_correction = s% job% relax_entropy_scale_max_correction + relax_entropy_retry_hold = s% job% relax_entropy_retry_hold relax_entropy_filename = s% job% relax_entropy_filename get_entropy_for_relax_from_eos = s% job% get_entropy_for_relax_from_eos report_cell_for_xm = s% job% report_cell_for_xm diff --git a/star_data/private/star_job_controls.inc b/star_data/private/star_job_controls.inc index 2267f94abf..22caacce4c 100644 --- a/star_data/private/star_job_controls.inc +++ b/star_data/private/star_job_controls.inc @@ -259,6 +259,10 @@ real(dp) :: timescale_for_relax_entropy real(dp) :: max_dt_for_relax_entropy real(dp) :: num_timescales_for_relax_entropy + logical :: relax_entropy_use_implicit_source + logical :: relax_entropy_use_normalized_source + real(dp) :: relax_entropy_scale_max_correction + integer :: relax_entropy_retry_hold character (len=strlen) :: relax_entropy_filename character (len=strlen) :: get_entropy_for_relax_from_eos