Menu
▾
▴
Re: [mesa-users] bad angular_momentum_j
|
From: Robert F. <rjf...@as...> - 2016-08-12 20:59:47
|
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 <pa...@gm...> 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 <kv...@ua...> 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 <kv...@ua...> 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 <kv...@ua...> 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 <pa...@gm...> 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 <kv...@ua...> 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 <pa...@gm...> 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 <kv...@ua...> 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 <pa...@gm...>
>>>>>>>> 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 <kv...@ua...>
>>>>>>>>> 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 <kv...@ua...> 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 <pa...@gm...>
>>>>>>>>>>> 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 <
>>>>>>>>>>>> pa...@gm...> 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 <kv...@ua...>
>>>>>>>>>>>>> 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
>>>>>>>>>>>>>> mes...@li...
>>>>>>>>>>>>>> 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
> mes...@li...
> https://lists.sourceforge.net/lists/listinfo/mesa-users
>
>
|