[Mesa-users] Loss of angular moment
Roque Caballero Navarro
rcaballeron at hotmail.com
Thu Oct 11 15:15:31 EDT 2018
Dear all,
I've been trying to redistribute (following Pablo's suggestion) the calculated j_dot among the cells which make up the outermost convective region of my model but without success. The simulation always breaks because the extra_jdot values (below the source code and results) are too high in comparison with the corresponding j_rot(k) ones calculated by MESA which, depending on the cell and for my simulation, are are between 10e11 - 10e15.
The j_dot value obtained after applying the different steps (1)-(10) is -7.30e29 which I've tried to compare with the literature but I couldn't find any similar reference for Sun alike model. My impression is that the value is properly calculated but I'm wondering if it can't be use directly without taking into consideration the elapse time between the current and previous step, could this be my error? If not, could you point me to any other suggestion?
Again, many thanks in advance for your help.
Kind regards,
Roque
!Simulation executed for a 1 solar mass model
(1) r_st = s% r(1)
(2) m_st = s% m(1)
(3) omega_surf = s% omega(1)
(4) v_esc = (618 * ((Rsun/r_st)*(m_st/Msun))**0.5) * 100000 !from km/s to cm/s
(5) v_inf = 1.92 * v_esc
(6) B = s% x_ctrl(6) !1G in this example
(7) m_dot = s% mstar_dot
(8) eta_surf = ((r_st * B)**2)/(abs(m_dot) * v_inf)
(9) j_dot = two_thirds * m_dot * omega_surf * r_st**2 * eta_surf
!Distribute j_dot among the outermost convective region
(10) call distribute_j_dot(j_dot, s% extra_jdot)
!Debug output
Rsun= 69598000000.000000 ! from MESA
Msun= 1.9891999999999999E+033 ! from MESA
r_st= 64770270940.741280 ! from MESA
m_st= 1.9891951638625125E+033 ! from MESA
v_esc= 64061701.595639318 ! calculated
v_inf 122998464.31439751 ! calculated
B 1.0000000000000000 !from inlist
m_dot -1974016097772.9548 ! from MESA
eta_surf 17.278301516747817 ! calculated
omega_surf 7.7296419108289482E-006 ! from MESA
j_dot -7.3734592048942183E+029 ! calculated and to be distributed
i_st= 2.7967919984913558E+021 ! from MESA
On 05/15/2018 09:31 PM, Chris Mankovich wrote:
Hi Roque,
I’m trying to implement a magnetic braking routine according to the MESA school exercise proposed by Matteo Cantiello in 2012. The formulas on it are referred to the stellar surface and for that reason I’m calculating the loss of angular moment at the very outermost zone.
If you're referring to the exercise on magnetic braking, i.e., angular momentum loss in the form of Equation 2.1 from the document "Stellar Rotation II", then Jdot represents the loss in total angular momentum, not a specific angular momentum loss (which we'd probably denote jdot). Dimensionally it would make sense to divide this by s% dm(1) and store the result in s% extra_jdot(1), but is that what you want to model physically? As Pablo pointed out, it is likely to cause the outermost layer to counterrotate and it might be difficult to proceed. If you want to apply an angular momentum loss rate like this, I'd suggest that you calculate Jdot from the surface quantities as you have done, but then divide it by the mass of some region of the star (as Pablo also suggested) and set the resulting jdot throughout that region.
Cheers,
Chris
On Sun, May 13, 2018 at 1:01 PM Roque Caballero Navarro via Mesa-users <mesa-users at lists.mesastar.org<mailto:mesa-users at lists.mesastar.org>> wrote:
Hi Pablo,
Many thanks for your feedback.
I'm trying to implement a magnetic braking routine according to the MESA school exercise proposed by Matteo Cantiello in 2012. The formulas on it are referred to the stellar surface and for that reason I'm calculating the loss of angular moment at the very outermost zone.
Coming back to your proposal for calculating the value for extra_jdot(1), I've notice that you're using dm(1) instead of m(1) which is the value I'm using (code line 3) in my implementation. Now I'm wondering if I had to change lines 2 and 3 as follows in order to restrict the calculation just to the zone 1:
1. r_st = s% r(1) - s% r(2)
2. m_st = s% dm(1)
Does it makes sense to you?
Best regards,
Roque
On 05/13/2018 04:10 PM, Pablo Marchant wrote:
Hi Roque, I think your units are mixed up.
On Sun, May 13, 2018, 7:49 AM Roque Caballero Navarro via Mesa-users <mesa-users at lists.mesastar.org<mailto:mesa-users at lists.mesastar.org>> wrote:
Dear all,
I've ajusted the routine. Basically I don't divide now the star ratio by the Sun one as in the previous version.
1. !Star data
2. r_st = s% r(1)
3. m_st = s% m(1)
4. omega_surf = s% omega_avg_surf
5. m_dot = s% mstar_dot
6. !inlist parameter
7. B = s% x_ctrl(6)
8. v_esc = 618 * (((Rsun/r_st)*(m_st/Msun)))**0.5
9. v_inf = 1.92 * v_esc
10. eta_surf = abs((r_st**2 * B**2) / (m_dot * v_inf))
11.
12. !Lost of J
13. j_dot = two_thirds * m_dot * omega_surf * r_st**2 * eta_surf
14. s% extra_jdot(1) = j_dot
j_dot in here is angular momentum per unit time, s% extra_jdot(1) needs to have units of angular momentum per unit mass per unit time.
Also, it's very likely that you want to remove this angular momentum from a physically defined region of the star, rather than the very outermost cell. Just using the outer layer, to make units consistent, you'd have to do
s% extra_jdot(1) = j_dot/s% dm(1)
This will likely be very unstable, and cause your star to counterrotate at the very outer layer unless you use very small timesteps.
Cheers
Now the j_dot value is much bigger (10**39) than the angular moment at surface (10**16) and the simulation breaks. Here the different partial values obtained:
1. Rsun= 69598000000.000000
2. Msun= 1.9891999999999999E+033
3. r_st= 63847123037.660912
4. m_st= 1.9891759644537749E+033
5. v_esc= 645.22853372743293
6. v_inf= 1238.8387570664665
7. B= 100.00000000000000
8. m_dot= -7650742677799.6270
9. eta_surf= 4300948934.7854061
10. omega_surf= 3.1121669509221367E-005
11. j_dot= -2.7830574679286024E+039
12. s% i_rot(1)= 2.7176367467908070E+021
13. s% omega(1)= 3.1311246486658475E-005
Does any one have a clue about what I'm doing wrong?
Kind regards,
Roque
On 04/23/2018 11:37 AM, Rob Farmer wrote:
>Don't you think the same?
Ah i never ran your model long enough, its always easiest (and faster) to send models that are "near" problem points, so people don't waste time constantly rerunning uninteresting bits of evolution.
I think you have a scale issue in your calculation:
omega(1) ~ 10**-5
irot(1) ~10**21
jdot(1) ~10**-10
So the existing angular momentum of the surface is ~ 10**21 * 10**-5 = 10**16 to which you then try to -10**-10, even your -1 value is small compared to the existing angular momentum.
> could it be that I'm not propagating properly back the loss of angular moment?
if your not sure what the codes does, have a look at the code:
grep -rinI other_torque $MESA_DIR/star/p*/*
The look in star/private/solve_omega_mix.f90 then follow where extra_jdot is used.
Rob
On 21 April 2018 at 16:52, Roque Caballero Navarro <rcaballeron at hotmail.com<mailto:rcaballeron at hotmail.com>> wrote:
Hi Rob,
I agree with you that it's needed a rotating model and I get it using the following configuration:
! Rotation off until near ZAMS
change_rotation_flag = .false.
new_rotation_flag = .true.
new_surface_rotation_v = 20 ! solar (km sec^1)
set_near_zams_surface_rotation_v_steps = 10
The star starts to rotate with it approaches to the ZAMS. Using the "debug" option of my routine, I can confirm the other_torque hook is called as soon as the rotation is activated. I also dump the calculated loss of angular moment (and other partial results) in the history fie and non-zero values are reported. I'd say that the code is properly invoked. Don't you think the same? I'll do also some test with the rotation activated since the beginning.
Regarding s% omega_avg_surf value, I've also checked that isn't null when other_torque is invoked. The same values are reported in the history file as well. Nevertheless, I'll follow your advice and take a deeper look into star_utils.
So Rob, for me the question is, if the model is rotating and other_torque routine is invoked (let's put apart if the calculations done on it are right or not), could it be that I'm not propagating properly back the loss of angular moment? I'm assuming that even if I used an "absurd" value, an impact on s% omega_avg_surf value, right?
Many thanks for your help.
Kind regards,
Roque
On 04/20/2018 05:27 PM, Rob Farmer wrote:
Hi
To use other_torque you need a rotating model, you need to set:
change_rotation_flag = .true.
in your star_job inlist to change the rotation flag, new_rotation_flag only changes the rotation flag if change_rotation_flag = .true.
The other issue is your using s% omega_avg_surf, but that gets set after the other_torque call. When using the other_ routine the best bet is to always "goto the source" of the stellar structure data, rather than any use any shorthand variables as they may not be set. In this case you'll want to look at the s% omega array which stores the rotation information for the star. See set_surf_avg_rotation_info in star/private/star_utils.f90 for how s% omega_avg_surf is implemented, you'll want to copy the relevant bits in your run_star_extras file for computing the surface averaged value.
Rob
On 20 April 2018 at 14:18, Roque Caballero Navarro <rcaballeron at hotmail.com<mailto:rcaballeron at hotmail.com>> wrote:
Hi Rob,
Please find attached the requested files.
Kind regards,
Roque
________________________________
De: Rob Farmer <r.j.farmer at uva.nl<mailto:r.j.farmer at uva.nl>>
Enviado: viernes, 20 de abril de 2018 12:12
Para: Roque Caballero Navarro
Cc: mesa-users at lists.mesastar.org<mailto:mesa-users at lists.mesastar.org>
Asunto: Re: [Mesa-users] Loss of angular moment
Hi
This is not enough information to help you. Please provide your inlists and run_star_extras.f
Rob
On 20 April 2018 at 12:06, Roque Caballero Navarro via Mesa-users <mesa-users at lists.mesastar.org<mailto:mesa-users at lists.mesastar.org>> wrote:
Dear all,
I'm trying to implement a magnetic braking routine based on Matteo Cantiello exercise propose during the Mesa school 2012 (here the details<http://cococubed.asu.edu/mesa_market/ewExternalFiles/2012_cantiello.zip>.) Well, for doing it I'm using the hook other_torque and the code is implemented, as usual, in the run_extras_star.f file. Below you can see the core part of the implementation:
1. if (s% use_other_torque) then
2.
3. !Star data
4. r_st = s% r(1)
5. m_st = s% m(1)
6. omega_surf = s% omega_avg_surf
7. m_dot = s% star_mdot
8.
9. v_esc = 618 * (((Rsun/r_st)*(m_st/Msun)))**0.5
10.
11. v_inf = 1.92 * v_esc
12.
13. B = s% x_ctrl(6)
14.
15. eta_surf = abs(((r_st/Rsun)**2/B**2)/(m_dot * v_inf))
16.
17. j_dot = two_thirds * m_dot * omega_surf * (r_st/Rsun)**2 * eta_surf
18.
19. s% extra_jdot(1) = j_dot
20. !s% extra_jdot(1) = -1
21.
22. if (debug_use_other_torque) then
23. write(*,*) "Rsun=", Rsun, "Msun=", Msun, "r_st=", r_st, "m_st=", m_st, &
24. "v_esc=", v_esc, "v_inf", v_inf, "B", B, "m_dot", m_dot, "eta_surf", eta_surf, &
25. "omega_surf", omega_surf, "j_dot", j_dot
26. end if
27.
28. s% x_ctrl(7) = v_esc
29. s% x_ctrl(8) = v_inf
30. s% x_ctrl(9) = eta_surf
31. s% x_ctrl(10) = j_dot
32. s% x_ctrl(11) = m_dot
33.
34. end if
The input parameters for the calculation are the radius, mass, loss of mass and velocity at surface (lines 4 to 7). After some calculations, I end up with a loss of angular moment j_dot which I pass back using the array extra_jdot (line 20). After running a simulation with this routine activated, I'd expect to have lower values of velocity at surface but I couldn't find any difference. I've even assigned an unreal value (line 20) to the extra_jdot array but nothing changed. So, something must be wrong or missing in the implementation. As you maybe has already noticed, I'm just providing a value to one cell, the outermost one. I'm assuming that only this layer is affected by the magnetic braking effect.
As usual, any help is more than welcome.
Kind regards,
Roque
_______________________________________________
mesa-users at lists.mesastar.org<mailto:mesa-users at lists.mesastar.org>
https://lists.mesastar.org/mailman/listinfo/mesa-users
_______________________________________________
mesa-users at lists.mesastar.org<mailto:mesa-users at lists.mesastar.org>
https://lists.mesastar.org/mailman/listinfo/mesa-users
--
Christopher R. Mankovich
PhD Candidate, Department of Astronomy & Astrophysics
University of California, Santa Cruz
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20181011/05f4d05e/attachment.html>
More information about the Mesa-users
mailing list