[mesa-users] bad angular_momentum_j
Pablo Marchant
pamarca at gmail.com
Fri Aug 12 16:37:53 EDT 2016
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20160812/4fc73f38/attachment.html>
More information about the Mesa-users
mailing list