Menu
▾
▴
Re: [mesa-users] bad angular_momentum_j
|
From: Pablo M. <pa...@gm...> - 2016-08-12 20:38:01
|
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 |
