[mesa-users] bad angular_momentum_j

Kenny Van kvan at ualberta.ca
Fri Aug 12 19:16:42 EDT 2016


Hi Robert,

I fixed the problem you pointed out in your email. Unfortunately this
doesn't fix the sharp drop in angular momentum. It seems like the
issue is the turnover time oscillates. The turnover time which I defined in
my code as:

dr / conv_vel(k)

where dr is the size of the cell and conv_vel is the convective velocity of
said cell.

>From the profile.data outputs it doesn't seem like either one of these
variables should have large jumps from model to model but the turnover
time has very large jumps.

On 12 August 2016 at 14:59, Robert Farmer <rjfarmer at asu.edu> wrote:

> Hi
>
> >if (k == 1) then
> >      dr = s% r(k) - s% R_center
> >else
>  >     dr = s% r(k) - s% r(k + 1)
> >end if
>
> Your assuming s%r(1) is the center, when its in fact the outer grid cell,
> s%r(s%nz) is the center. So either use this:
>
> if (k < s% nz) then
>               dr= s% r(k) - s% r(k+1)
>  else
>                dr= s% r(k) - s% R_center
> end if
>
> Or as dr is a possible output in the profile data you could always do:
>
> use star_lib
>
> dr=star_get_profile_output(s,'dr',k)
>
> This works for any profile output (there is also a version for history
> data too)
>
> What that does for your results i don't know just what i saw from a quick
> scan.
>
> Rob
>
> On Fri, Aug 12, 2016 at 1:37 PM, Pablo Marchant <pamarca at gmail.com> wrote:
>
>> Kenny, I won't have time to check this soon, so I'd recommend you debug
>> this in more detail. My impression would be that your formula for jdot is
>> misbehaving, I don't know exactly what are the quantities you are
>> outputting, but seeing the values that come out in the last step, many of
>> those jump by 3-4 orders of magnitude. Try an find the source of that jump.
>>
>> On Fri, Aug 12, 2016 at 7:50 PM, Kenny Van <kvan at ualberta.ca> wrote:
>>
>>> Hi,
>>>
>>> Unfortunately some of the models that I'm testing my modified MB on are
>>> running into this odd issue. The angular momentum decreases more quickly
>>> than the default
>>> MB implemented in MESA but with a sufficiently small timestep this
>>> doesn't appear to be a problem. What occurs after some time in evolution is
>>> the angular momentum
>>> suddenly goes negative. Now in the final step prior the the angular
>>> momentum becoming negative the jdot does not appear to be large enough to
>>> cause this sharp decrease.
>>>
>>> The angular momentum near the final step is on the order of 10^52, while
>>> the jdot is on order of 10^40. This jdot value is multiplied by the
>>> timestep size converted to seconds
>>> so the units are the same.
>>>
>>> I added the following two lines into binary_evolove.f90 before the
>>> angular momentum check
>>>
>>>          write (*,*) "angular_momentum_j", b% angular_momentum_j
>>>          write (*,*) "jdot", b% jdot * b% s_donor %time_step * secyer
>>>
>>> So the jdot is taken directly from binary_evolve just to be sure.
>>>
>>> I'm not sure if this is a problem that I have accidentally created with
>>> my run binary_extras.f but I'm not sure how to fix this.
>>>
>>>
>>> On 9 August 2016 at 10:37, Kenny Van <kvan at ualberta.ca> wrote:
>>>
>>>> Hi Pablo,
>>>>
>>>> I figured out what the issue was. I didn't realize that the different
>>>> timestep values were in different units so that caused issues.
>>>>
>>>> Thanks for the help
>>>>
>>>> On 8 August 2016 at 17:06, Kenny Van <kvan at ualberta.ca> wrote:
>>>>
>>>>> Hi Pablo,
>>>>>
>>>>> I'll verify the problem I'm running into and get back to you.
>>>>>
>>>>> As for the snippet I added before, that was just an attempt to see if
>>>>> I could get MESA to decrease the timestep in the very first iteration. I
>>>>> include a similar code in extra_binary_check_model and just use a small
>>>>> initial timestep now.
>>>>>
>>>>> On 8 August 2016 at 17:02, Pablo Marchant <pamarca at gmail.com> wrote:
>>>>>
>>>>>> Kenny, if you want to debug this just do some small edits to
>>>>>> binary/private/binary_evolve.f90, in particular around these lines
>>>>>>
>>>>>> 257          ! get jdot and update orbital J
>>>>>> 258          b% jdot = get_jdot(b, b% mtransfer_rate, b%
>>>>>> xfer_fraction)
>>>>>> 259          b% angular_momentum_j = b% angular_momentum_j + b%
>>>>>> jdot*b% s_donor% time_step*secyer
>>>>>> 260
>>>>>> 261          if (b% angular_momentum_j <= 0) then
>>>>>> 262             stop 'bad angular_momentum_j'
>>>>>> 263          end if
>>>>>>
>>>>>> Since the stop condition literally has a check for negative angular
>>>>>> momentum, I have a hard time believing what you say. You can put write
>>>>>> statements in there to verify you are getting the numbers you expect. If
>>>>>> this still gives you no answers, then if you want you can send a minimal
>>>>>> example that fails (truly minimal please), and I can give it a look.
>>>>>>
>>>>>> PS:
>>>>>>
>>>>>> Just saw in more detail this snippet you posted before
>>>>>>
>>>>>>          write(*,*) "starting modified MB"
>>>>>>          if (b% angular_momentum_j <= 0) then
>>>>>>             write(*,*), "angular_momentum_j <= 0", b%
>>>>>> angular_momentum_j
>>>>>>             b% s1% dt_next = min(b% s1% dt * 0.50, b% s1% dt_next)
>>>>>>             extras_binary_startup = retry
>>>>>>          end if
>>>>>>
>>>>>> You said you add this to extras_binary_startup. That will do nothing,
>>>>>> extras_binary_startup is only ran at the beginning of the simulation, and
>>>>>> orbital angular momentum will be positive. You should be doing this in
>>>>>> extras_binary_check_model.
>>>>>>
>>>>>> On Mon, Aug 8, 2016 at 9:12 PM, Kenny Van <kvan at ualberta.ca> wrote:
>>>>>>
>>>>>>> Hi,
>>>>>>>
>>>>>>> So unfortunately after some testing I've started running into some
>>>>>>> problems that I cant quite solve. I've included a scaling factor in the
>>>>>>> timestep check to try and retry the run if we encounter a problem with
>>>>>>> angular momentum.
>>>>>>> The problem that I'm running into now is that the code triggers the
>>>>>>> "bad angular_momentum_j" portion even when the angular_momentum_j > 0. At
>>>>>>> every step I have the extras_binary_check_model output the
>>>>>>>
>>>>>>> b% angular_momentum_j,
>>>>>>>
>>>>>>> as well as
>>>>>>>
>>>>>>> b% angular_momentum_j + b% jdot * b% s1% dt_next
>>>>>>>
>>>>>>> so I have an idea of what the current angular momentum is as well as
>>>>>>> approx what it will be in the next timestep.
>>>>>>>
>>>>>>> At the point of the code breaking both of these values are positive
>>>>>>> and I encounter the termination code
>>>>>>>
>>>>>>> STOP bad angular_momentum_j.
>>>>>>>
>>>>>>> I'm not quite sure what I can do from here.
>>>>>>>
>>>>>>> On 4 August 2016 at 21:40, Pablo Marchant <pamarca at gmail.com> wrote:
>>>>>>>
>>>>>>>> You should be able to do that in extras_binary_startup. Anyhow, I'd
>>>>>>>> just go for the simple approach of setting a small initial timestep, if its
>>>>>>>> too small it will grow quite fast.
>>>>>>>>
>>>>>>>> On Fri, Aug 5, 2016 at 3:36 AM, Kenny Van <kvan at ualberta.ca> wrote:
>>>>>>>>
>>>>>>>>> Yes, but I was also hoping to create an adaptive way for MESA to
>>>>>>>>> adjust the first step if necessary
>>>>>>>>>
>>>>>>>>> On 4 August 2016 at 17:00, Pablo Marchant <pamarca at gmail.com>
>>>>>>>>> wrote:
>>>>>>>>>
>>>>>>>>>> Are you running into these issues at the first timestep? If so,
>>>>>>>>>> you need to ensure the first step is small as well.
>>>>>>>>>>
>>>>>>>>>> On Thu, Aug 4, 2016 at 11:02 PM, Kenny Van <kvan at ualberta.ca>
>>>>>>>>>> wrote:
>>>>>>>>>>
>>>>>>>>>>> Hi Everyone,
>>>>>>>>>>>
>>>>>>>>>>> So I've checked the units for Jdot and they seem fine, Jorb/Jdot
>>>>>>>>>>> is on the order of 1e12. The issue I'm currently running into is that the
>>>>>>>>>>> code doesn't quite seem to be doing what I think it should be doing.
>>>>>>>>>>> For extra checks I also included:
>>>>>>>>>>>
>>>>>>>>>>>          write(*,*) "starting modified MB"
>>>>>>>>>>>          if (b% angular_momentum_j <= 0) then
>>>>>>>>>>>             write(*,*), "angular_momentum_j <= 0", b%
>>>>>>>>>>> angular_momentum_j
>>>>>>>>>>>             b% s1% dt_next = min(b% s1% dt * 0.50, b% s1%
>>>>>>>>>>> dt_next)
>>>>>>>>>>>             extras_binary_startup = retry
>>>>>>>>>>>          end if
>>>>>>>>>>>
>>>>>>>>>>> to extra_binary_startup. Just to check if my initial timestep is
>>>>>>>>>>> sufficiently small and the MESA code doesn't seem to be triggering this. I
>>>>>>>>>>> have a similar portion of code in extras_binary_check_model,
>>>>>>>>>>> but I haven't reached a point where the simulation breaks to
>>>>>>>>>>> determine if this works or not. Does the file binary_evolve.f90 trigger and
>>>>>>>>>>> instantly kill the MESA simulation before the code checks the
>>>>>>>>>>> run_binary_extras.f file?
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> On 3 August 2016 at 21:25, Kenny Van <kvan at ualberta.ca> wrote:
>>>>>>>>>>>
>>>>>>>>>>>> Hey Pablo,
>>>>>>>>>>>>
>>>>>>>>>>>> Thanks for the advice, I'll try out the time step adjustment to
>>>>>>>>>>>> see how it goes. I'll also double check the units to be 100% sure that I
>>>>>>>>>>>> haven't crated any jdots that are a few orders too big.
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> On 3 August 2016 at 20:53, Pablo Marchant <pamarca at gmail.com>
>>>>>>>>>>>> wrote:
>>>>>>>>>>>>
>>>>>>>>>>>>> Also, be sure you are getting your units right. That's an easy
>>>>>>>>>>>>> way to create absurdly large jdots (been there). For the process you are
>>>>>>>>>>>>> considering, what's the decay timescale for Jorb (ie. Jorb/Jdot)? If that
>>>>>>>>>>>>> is absurdly tiny then I'd worry.
>>>>>>>>>>>>>
>>>>>>>>>>>>> On Thu, Aug 4, 2016 at 4:50 AM, Pablo Marchant <
>>>>>>>>>>>>> pamarca at gmail.com> wrote:
>>>>>>>>>>>>>
>>>>>>>>>>>>>> Kenny, I'd recommend you to manually adjust the timestep,
>>>>>>>>>>>>>> this is something we discussed here in mesa-users a while ago, copy from
>>>>>>>>>>>>>> Rob Farmer's answer to create a hard timestep limit
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> (Untested)
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> In extras_startup:
>>>>>>>>>>>>>> s%xtra1=0.d0
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> extras_finish_step:
>>>>>>>>>>>>>> s%xtra1=s%log_R
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> extras_check_model:
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> if( abs( 10**s%xtra1_old - 10**s%log_R) > EPS) then
>>>>>>>>>>>>>>     extras_check_model = retry
>>>>>>>>>>>>>>     s% dt = s%dt * SOME_SCALE_FACTOR
>>>>>>>>>>>>>> end if
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> And to create a soft timestep limit:
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> (also untested)
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> if( abs( 10**s%xtra1_old - 10**s%log_R) > EPS) then
>>>>>>>>>>>>>>     s% dt_next = min(s% dt_next, s%dt * SOME_SCALE_FACTOR)
>>>>>>>>>>>>>> end if
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> If the first timestep is the issue, then set a smaller
>>>>>>>>>>>>>> initial timestep. Check this from defaults/star_job.defaults
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>  435          !### set_initial_dt
>>>>>>>>>>>>>>  436          !### years_for_initial_dt
>>>>>>>>>>>>>>  437          !### seconds_for_initial_dt
>>>>>>>>>>>>>>  438          ! if true, set initial timestep, dt, in years
>>>>>>>>>>>>>>  439
>>>>>>>>>>>>>>  440       set_initial_dt = .false.
>>>>>>>>>>>>>>  441       years_for_initial_dt = 0
>>>>>>>>>>>>>>  442       seconds_for_initial_dt = 0
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> On Wed, Aug 3, 2016 at 4:56 PM, Kenny Van <kvan at ualberta.ca>
>>>>>>>>>>>>>> wrote:
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> Hi,
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> I'm currently working on adapting additional magnetic
>>>>>>>>>>>>>>> braking terms into MESA using the run_binary_extras.f file. I'm currently
>>>>>>>>>>>>>>> running into an issue where the amount of angular momentum being removed is
>>>>>>>>>>>>>>> too great in a single timestep and causing the simulation to break. Looking
>>>>>>>>>>>>>>> at the code it seems like the MESA binary evolution dies immediately if the
>>>>>>>>>>>>>>> angular momentum loss is too great instead of retrying with a smaller
>>>>>>>>>>>>>>> timestep. Is there a way to get MESA to retry with a smaller timestep when
>>>>>>>>>>>>>>> it encounters this issue?
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> Thanks
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> ------------------------------------------------------------
>>>>>>>>>>>>>>> ------------------
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> _______________________________________________
>>>>>>>>>>>>>>> mesa-users mailing list
>>>>>>>>>>>>>>> mesa-users at lists.sourceforge.net
>>>>>>>>>>>>>>> https://lists.sourceforge.net/lists/listinfo/mesa-users
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> --
>>>>>>>>>>>>>> Pablo Marchant Campos
>>>>>>>>>>>>>> M.Sc on Astrophysics, Universidad Católica de Chile
>>>>>>>>>>>>>> PhD student, Argelander-Institut für Astronomie
>>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>> --
>>>>>>>>>>>>> Pablo Marchant Campos
>>>>>>>>>>>>> M.Sc on Astrophysics, Universidad Católica de Chile
>>>>>>>>>>>>> PhD student, Argelander-Institut für Astronomie
>>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>> Pablo Marchant Campos
>>>>>>>>>> M.Sc on Astrophysics, Universidad Católica de Chile
>>>>>>>>>> PhD student, Argelander-Institut für Astronomie
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> --
>>>>>>>> Pablo Marchant Campos
>>>>>>>> M.Sc on Astrophysics, Universidad Católica de Chile
>>>>>>>> PhD student, Argelander-Institut für Astronomie
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>> --
>>>>>> Pablo Marchant Campos
>>>>>> M.Sc on Astrophysics, Universidad Católica de Chile
>>>>>> PhD student, Argelander-Institut für Astronomie
>>>>>>
>>>>>
>>>>>
>>>>
>>>
>>
>>
>> --
>> Pablo Marchant Campos
>> M.Sc on Astrophysics, Universidad Católica de Chile
>> PhD student, Argelander-Institut für Astronomie
>>
>> ------------------------------------------------------------
>> ------------------
>> What NetFlow Analyzer can do for you? Monitors network bandwidth and
>> traffic
>> patterns at an interface-level. Reveals which users, apps, and protocols
>> are
>> consuming the most bandwidth. Provides multi-vendor support for NetFlow,
>> J-Flow, sFlow and other flows. Make informed decisions using capacity
>> planning reports. http://sdm.link/zohodev2dev
>> _______________________________________________
>> mesa-users mailing list
>> mesa-users at lists.sourceforge.net
>> https://lists.sourceforge.net/lists/listinfo/mesa-users
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20160812/2de47785/attachment.html>


More information about the Mesa-users mailing list