[Mesa-users] Run fails to stop and save using custom accreted_material_j

Mon May 1 16:33:12 UTC 2023

Hi Ebraheem,

That did the trick! I guess I’m not totally following why this causes a problem only for variable angular momentum accretion and only at the point of the stopping condition?

And thanks for the suggestions for the other controls!

Cheers,
- Ben

On Apr 29, 2023, at 1:21 AM, Ebraheem Farag <ekfarag at asu.edu> wrote:

Hello Benjamin,

Comment out "max_star_mass_for_gain = 0.101", this is telling MESA to turn off the mass_gain if the star_mass increases beyond 0.101 Msun during the last timestep, which causes the solver to struggle and retry if it overshoots the stopping condition with mass-loss on. This should resolve your issue.
-------------
Some extra suggested controls to adopt
-------------
To hit the stopping condition more precisely, I would recommend setting:
! for 5th decimal precision on stopping condition, e.g. 0.1010[1...]
when_to_stop_rtol = 1d-5 ! relative tolerance
when_to_stop_atol = 1d-5 ! absolute tolerance

I would also recommend finding a control that helps limit the timesteps, since the solver is max_increasing until it hits the stopping condition
For your case where you are only acretting 1d-3 Msun, perhaps the following will help numerically resolve your model, especially near the end where the mass loss is greatest:
delta_lg_star_mass_limit = 1d-4 ! limits mass change per timestep, I would lower or increase depending on how much mass you are accreting.
delta_lg_star_mass_hard_limit = 1d-3 ! hard limit for mass_change per timestep

Cheers :)

-EbF

On Fri, Apr 28, 2023 at 8:48 PM Roulston, Benjamin via Mesa-users <mesa-users at lists.mesastar.org<mailto:mesa-users at lists.mesastar.org>> wrote:
Hi all,

I have an odd issue with accreting mass and angular momentum that I can’t seem to understand why it is an issue.

Background is I’m trying to accrete a specific abundance to some grid of models. I’m also trying to keep track of rotation and the angular momentum. So I am using

mass_change = 1.0d-8
max_star_mass_for_gain = 0.101
star_mass_max_limit = 0.101

use_accreted_material_j = .true.
accreted_material_j = 1.0d12 ! this IS specific ang momentum

accretion_dump_missing_metals_into_heaviest = .false.
accrete_given_mass_fractions = .true.

And this works fine! The mass is added as expected, and so is angular momentum. The run stops as expected when total mass hits 0.101Msun. Great!

The issue is I want a more dynamic accreted_j, based on some fraction of the Keplarian value as the mass and radius of the accretor changes. Fine, there is a test case (that does run fine with my install) that does this! So I used that code in the run_star_extras.f90 to adapt to my case.

subroutine extras_controls(id, ierr)
integer, intent(in) :: id
integer, intent(out) :: ierr
type (star_info), pointer :: s
ierr = 0
call star_ptr(id, s, ierr)
if (ierr /= 0) return
! this is the place to set any procedure pointers you want to change
! e.g., other_wind, other_mixing, other_energy (see star_data.inc)

s% other_adjust_mdot => accretor_adjust_mdot

subroutine accretor_adjust_mdot(id, ierr)
integer, intent(in) :: id
integer, intent(out) :: ierr
type(star_info), pointer :: s
call star_ptr(id, s, ierr)

s% accreted_material_j = s% x_ctrl(1)*sqrt(standard_cgrav * s% mstar * s% photosphere_r*Rsun)
write(*,*) "accreted_material_j = ", s% accreted_material_j

!write(*,*) "debug", s% mstar_dot/Msun*secyer, 10**(s% x_ctrl(2))
! s% mstar_dot = s% mstar_dot + pow(10d0, s% x_ctrl(2))*Msun/secyer

end subroutine accretor_adjust_mdot

And

mass_change = 1.0d-8
max_star_mass_for_gain = 0.101
star_mass_max_limit = 0.101

accrete_same_as_surface = .false.
premix_omega = .true.

use_accreted_material_j = .true.
use_other_adjust_mdot = .true.

x_ctrl(1) = 0.1d0 ! ratio of accreted am to keplerian

This is where it gets weird. This DOES run! It accretes the mass fine and the angular momentum. However, once the model reaches the mass limit (in this model 0.101Msun) instead of stopping and saving the model, it keeps trying to accrete and the time steps drop down to log_dt ~ -9 and never stops. It gives the following errors in the terminal output:

 accreted_material_j =    29212373139847368.
         83   6.616617   3687.877  -2.810333  -2.810333   0.100900   0.100900   0.758355   0.000009   0.241465   6.528487    672      0
 4.1761E+00   6.616617  -1.037625 -99.000000  -4.131714  -8.000000   0.000000   0.241465   0.000084   0.000180   0.610681      4
 1.0001E+09   2.997841  -2.853451 -18.663860 -99.000000  -3.239033   0.000000   0.000031   0.000018   0.000180 -0.304E-09  max increase

 accreted_material_j =    29228061448833992.
 retry: set_omega_adjust_mass failed in adjust mass      84
 accreted_material_j =    29228061448833992.
         84   6.617098   3688.133  -2.806102  -2.806102   0.100990   0.100990   0.758349   0.000009   0.241471   6.528698    675      1
 3.9543E+00   6.617098  -1.037737 -99.000000  -4.127401  -8.000000   0.000000   0.241471   0.000084   0.000180   0.610356      4
 1.0001E+09   2.998586  -2.853554 -18.186738 -99.000000  -3.238683   0.000000   0.000031   0.000018   0.000180 -0.335E-09         retry

 accreted_material_j =    29237325035824144.
 retry: set_omega_adjust_mass failed in adjust mass      85
 accreted_material_j =    29237325035824144.
 retry: set_omega_adjust_mass failed in adjust mass      85
 accreted_material_j =    29237325035824144.
 retry: set_omega_adjust_mass failed in adjust mass      85
 accreted_material_j =    29237325035824144.
 retry: set_omega_adjust_mass failed in adjust mass      85
 accreted_material_j =    29237325035824144.

However, if I use the same code and functions, but just use a constant momentum

s% accreted_material_j = 1d12

It does work and the model stops and saves. But if I use a larger value like 1d18 it also doesn’t stop the run. So it seems it has issues with larger specific angular momentum. But the model runs fine until the exact point of the mass limit. And it’s not this specific mass limit, if I change the mass limit to 0.2 (or anything else) it runs fun until 0.2 (so it passes the previous limit of 0.101 fine) and then doesn’t stop and the time step limit drops to very low values again.

I can’t seem to understand why the model can run fine until the exact point it’s supposed to be stopping. But only for “large” (really not that large) specific angular momentum.

I hope I’ve explained this clearly, and I’ve included the necessary files to recreate this as well as my starting models. I’m using MESA version "mesa-r22.11.1”.

Cheers,
- Ben
_______________________________________________
mesa-users at lists.mesastar.org<mailto:mesa-users at lists.mesastar.org>
https://urldefense.com/v3/__https://lists.mesastar.org/mailman/listinfo/mesa-users__;!!IKRxdwAv5BmarQ!fw64IBlmWGjNsXCUvzv5ZXbsKzZ3D_8LCFryzQjfC0CJIcBLPdSkF3TGJtPoLCQB2KVZOYEFbMa8gZxYpb81Q1dK$

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20230501/82536bbe/attachment.htm>