Skip to content

6: Relaxation, mass adjustment, and shell removal #1053

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
Debraheem wants to merge 1 commit into from
Draft

6: Relaxation, mass adjustment, and shell removal #1053

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 15 additions & 0 deletions docs/source/changelog.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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 <freedman_lowt_z_bug>` and `gh-993 <https://github.com/MESAHub/mesa/pull/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 <plasmon_weinberg_angle_bug>` and `gh-998 <https://github.com/MESAHub/mesa/pull/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.
Expand Down
29 changes: 29 additions & 0 deletions star/defaults/star_job.defaults
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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 = ''

Expand Down
85 changes: 82 additions & 3 deletions star/private/adjust_mass.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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) :: &
Expand All @@ -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(:)

Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down
16 changes: 12 additions & 4 deletions star/private/eps_mdot.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -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)

Expand All @@ -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
Expand All @@ -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

Expand Down
6 changes: 4 additions & 2 deletions star/private/evolve.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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
Expand Down
Loading