[mesa-users] questions about the "Eulerian" scheme for eps_grav
Masanobu Kunitomo
kunitomo at geo.titech.ac.jp
Sun Nov 16 02:09:56 EST 2014
Hi Dean,
I am sorry for the late reply. Thanks again for answering my questions!
1) time step control from abundance change
Thank you very much for the explanation and files for the time step control from abundance change. However, I have not succeeded in implementing it with ver. 7184 so far.
I downloaded files from http://mesastar.org/teaching-materials/2014-mesa-summer-school-working-dir/townsley/Long%20Lab%201.zip/view <http://mesastar.org/teaching-materials/2014-mesa-summer-school-working-dir/townsley/Long%20Lab%201.zip/view> and I found that the “run_star_extras.f" file is not in the same format as that in ver.7184, so I modified it. There are three minor modifications:
● modified the format of the functions and subroutines (e.g., subroutine extras_controls(s, ierr) ---> subroutine extras_controls(id, ierr) )
● in the function “extras_finish_step”, I change “eps_grav_time_deriv_separation” to “s% eps_grav_time_deriv_separation"
● commented out “call system_clock(time0,clock_rate)” in the function “extras_startup"
Please check the attached file “run_star_extras.f”. I succeeded in compiling, but just after running (./rn), I got the following message and the calculation terminated:
----------------------------------------------------------------------------------------------------------------------
Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
#0 0x10f030742
#1 0x10f030ed0
#2 0x7fff8afc15a9
#3 0x10c0d0d36
#4 0x10c3943b9
#5 0x10c395ced
#6 0x10c0d3d12
#7 0x10c0d3dae
#8 0x10c0d3de7
./rn: line 2: 21670 Segmentation fault: 11 ./star_vlm
----------------------------------------------------------------------------------------------------------------------
I apologise if my modifications were not appropriate. I have not solved this problem yet.
Thank you very much for the explanation about the advection of species. I understand that we should be careful if there are nuclear reaction in the constant-q region. (Fortunately, I think it never happens in my case.)
2) setting q_below_Eulerian_eps_grav
> For your last question, there is a parameter max_q_for_k_below_const_q, which if you set to 0.99 I think will keep 1% of the star as the constant-q type mesh. If it doesn't do that, I should fix things so that it does do that, as that is what I think that control is supposed to do.
Thank you very much for the suggestion and it worked.
I calculated with the following settings in inlist. The setting "q_for_k_const_mass = 0.95” worked. However, the calculation failed to converge and “q_for_k_below_const_q” and “q_below_Eulerian_eps_grav” are still sometimes very close to 1. Please check the attached figure “t-q.pdf".
--------------------------------------------------------------------------
max_q_for_k_below_const_q = 0.95
min_q_for_k_below_const_q = 0.9
max_q_for_k_const_mass = 0.95
min_q_for_k_const_mass = 0.9
--------------------------------------------------------------------------
Next, I found that the “q_for_k_below_const_q” and “q_below_Eulerian_eps_grav” region can be changed with the “extras_finish_step” function as follows. However, the convergence was worse, although I succeeded in avoid the fluctuation in “q_below_Eulerian_eps_grav” as shown in the attached figure, “t-q2.pdf".
--------------------------------------------------------------------------
! q_below_Eulerian_eps_grav, q_const_mass < 0.9M
do j=1,s% nz
if (s% m(j) < 0.9d0*s% star_mass*msol ) exit
end do
s% k_below_Eulerian_eps_grav = max(j,s% k_below_Eulerian_eps_grav)
s% k_const_mass = max(j,s% k_const_mass)
! q_below_Eulerian_eps_grav < 0.99M
do j=1,s% nz
if (s% m(j) < 0.99d0*s% star_mass*msol ) exit
end do
s% k_below_const_q = max(j,s% k_below_const_q)
! k_Lagrangian_eps_grav
s% k_Lagrangian_eps_grav = s% k_below_Eulerian_eps_grav + 10
--------------------------------------------------------------------------
In the calculation in t-q2.pdf, I used the following default settings:
--------------------------------------------------------------------------
max_q_for_k_below_const_q = 1.0
min_q_for_k_below_const_q = 0.999
max_q_for_k_const_mass = 1.0
min_q_for_k_const_mass = 0.995
--------------------------------------------------------------------------
3) smaller time step in retry and backup
Finally, I found that the convergence is improved by the smaller time steps in retry and backup (“timestep_factor_for_backups=1.0d-1” and “timestep_factor_for_retries=1.0d-5”)! I used the default values, 0.5, for both.
Therefore, the fluctuation in the q or m coordinate regions seems to be irrelevant to the convergence in my study, although I struggled to change the q coordinate region.
Thanks again for your many helps!
Best wishes,
Masanobu
> On Nov 13, 2014, at 15:13, Dean Townsley <Dean.M.Townsley at ua.edu> wrote:
>
> Hi Masanobu,
>
> Sorry to be a little slower getting back to you this time.
>
> The timestep control that we use, implemented by my grad student Broxton Miles, is included in the run_star_extras.f in our Long Lab 1 from the 2014 MESA summerschool:
> http://mesastar.org/teaching-materials/2014-mesa-summer-school-working-dir/townsley/Long%20Lab%201.zip/view <http://mesastar.org/teaching-materials/2014-mesa-summer-school-working-dir/townsley/Long%20Lab%201.zip/view>
> It basically determines the depth of the uniform composition layers based on a gradient check and then reproduces the logic for finding the constant-q/Lagrangian boundary and uses that to determine the timestep necessary to keep the boundary within the uniform abundance region. The logic is maybe a little hard to follow, but it seems to work.
>
>
> My apologies for being vague about what the issue is with the constant-q form in non-uniform abundance regions. I'll try to clarify my experience and possible concerns. There is an option to include the composition gradient terms in the computation of eps_grav (=Ds/Dt) in the constant-q region, however when we tried turning these on and letting the constant-q region extend into non-uniform abundance regions (for large timesteps in accretion) it didn't behave as expected. The actual eps_grav was no longer independent of where the boundary between the mesh regions was, which it should be.
>
> Unfortunately I haven't gotten much further than that, so anything else I say is really speculation so don't hold me to it. (i.e. apologies to Bill if I'm wrong or misunderstanding something) But here goes: I believe that the advection of species in mesa/star is handled by directly remapping the abundances to the grid shells that make up the grid at the end of the timestep, and then diffusion is done "on top of" that. Though I think the composition gradient terms in eps_grav for the constant-q region are right, I haven't worked that out carefully. It might be that there is just some minor bug in their computation that is easily fixed. However, I'm also concerned that the way species transformations (nuclear reactions) are done may need a slight adjustment on the constant-q part of the mesh. In MESA applications so far this is not an issue because nuclear reactions always happen on the Lagrangian part of the mesh. That's probably the cleaner way to do things anyway.
>
>
>
> For your last question, there is a parameter max_q_for_k_below_const_q, which if you set to 0.99 I think will keep 1% of the star as the constant-q type mesh. If it doesn't do that, I should fix things so that it does do that, as that is what I think that control is supposed to do.
>
>
> Convergence issues are always tough. And I have found that figuring out the root cause can be tricky. Having the shifts in the regions with the different types of mesh is one of the drawbacks of using a constant q region that depends on timestep. The above control may help with this. If so let me know, since I would like to know if the change of the mesh regions can cause instability all by itself. As long as the esp_grav_deriv_separation is fairly small, the eps_grav computed at any given physical position in the star should be fairly close to the same regardless of how the mesh (and therefore the appropriate form for Ds/Dt) is chosen. Also falling back to a Lagrangian mesh should be safer.
>
> All that being said, it may be that if some things like the composition gradient terms or the reactions aren't quite fully consistent in the constant-q region that something bad might happen and then the integrator has trouble recovering once time step shortens and the mesh reverts back to mostly Lagrangian because the effective physics has suddenly changed.
>
>
> Again hopefully this is helpful. And do let us know what you find. Thanks for trying things out!
>
> Dean
>
>
> On 11/11/2014 02:49 PM, Masanobu Kunitomo wrote:
>> Hi Dean,
>>
>> Thank you very much for the very clear and prompt reply. I really appreciate your help.
>>
>> I agree not to use the Eulerian terminology and thank you for the explanation. I will use “q coordinate region” instead of “Eulerian region”.
>>
>>
>> I tried the case with "eps_grav_time_deriv_separation = -1” with ver. 7184, and I succeeded in making the entire star the “q coordinate” region! Unfortunately, the calculation failed to converge at least in my setting (intense accretion onto a low-mass protostar). Actually, I wanted to try such a case in order to compare with a previous study which seems to use this setting. As you pointed out, the entirely q coordinate setting is not a good way because of the numerical mixing, so I will use the MESA’s standard setting, namely with three regions: the q coordinate, the m coordinate and the blend region.
>>
>>
>> Thank you very much for the very clear explanation for the q coordinate calculation! I understand that I misunderstood the concept. I thought that the newly added mass is uniformly distributed in the q coordinate region, but in fact it is added to the surface and the mesh is adjusted. This procedure does not underestimate the entropy of accreting material and it would be the same as the surface entropy of the previous time step.
>> I am interested in the stellar surface composition, so your caution is very important for my calculation! If you don't mind, could you please tell me the way you limit the time step? Do you check the element gradient in the q coordinate region in each time steps?
>> (I am sorry but I am not sure I could completely understand your caution. I thought that if the element diffusion is included, the advection from the species gradient is calculated even in the Lagrangian region.)
>>
>>
>> Finally, let me ask one more question. I still have a problem in the case of an accreting protostar.
>> Sometimes the calculation failed to converge. It seems to be triggered by a short time step for some reasons at some point, and then the calculation becomes unstable (i.e., short time steps and many retries). I found that a short time step brings the base of q coordinate region close to the surface suddenly. This is because the short time step results in the small dm (the mass added in the current step). I guess this rapid change of the region trigger the convergence problem… I attached a figure which shows the evolution of the q coordinate region with time. At the end, “q_below_just_added” is very close to the surface, and then the Eulerian region becomes tiny.
>> To avoid this, I was hoping that I could retain a q coordinate region in some finite value (e.g., 1% total mass or more). Is this possible? (I guess it is determined by the subroutine “set_Eulerian_Lagrangian_for_eps_grav” in the file “mesa/star/private/evolve.f”. I tried some modifications such as “dxm_kA = s% xmstar” but a modification in it was not reflected in a calculation...)
>>
>> Best wishes,
>> Masanobu
>>
>>
>>
>>
>>
>>
>>
>>
>>> On Nov 11, 2014, at 03:39, Dean Townsley <Dean.M.Townsley at ua.edu <mailto:Dean.M.Townsley at ua.edu>> wrote:
>>>
>>> Hi Masanobu,
>>>
>>> Before going any further, if you don't mind I'd like to not use the "Eulerian" terminology because, for others on the list not familiar with this usage, the behavior we are talking about is not the usual definition of Eulerian. (Eulerian = coordinates and mesh boundaries fixed in space.) I'm sure you know this, but I want to be clear for others. There is an outer mesh region in which the shell boundaries are at the same coordinate in q=m/Mstar instead of just mass coordinates, m. This is the "Eulerian" region, and the usual term for these kinds of coordinates are homology coordinates.
>>>
>>>
>>> The current strategy implemented in MESA star is to try to make as much as possible of the star Lagrangian in order to avoid numerical mixing between zones. What you want is the opposite of this, in which the entire star is treated with homology coordinates. I think that you can accomplish your goal by setting eps_grav_time_deriv_separation = -1. Though I admit that I haven't tested this case.
>>>
>>> The choice of where the crossover from q coordinates to m coordinates in the grid happens is performed while constructing the "new" grid for the end of the current timestep. For accretion, the procedure starts at the surface placing shells with the same dq as the old grid, giving constant q coordinates. It then searches inward until a point is found where the shells are large enough that the division, in mesh spacing, between points in the old an new mesh that have the same mass coordinate, m, is smaller than eps_grav_time_deriv_separation. This is where it is considered "safe" to evaluate the Lagrangian derivative of entropy, Ds/Dt, directly by finite differencing across the time step at the same m coordinate. So I believe if you set eps_grav_time_deriv_separation to -1, it will not find such a place and will just make the new mesh similar to the existing one in q coordinates, and therefore evaluate Ds/Dt as (partial s /partial t)_q plus the advection term everywhere in the star. (Typically called the "Eulerian" form.)
>>>
>>> Let me know if that doesn't work for you and I can figure out if there needs to be some adjustment to make it work. Please also say if it does work, as we could add that to the comment for this control.
>>>
>>>
>>>
>>> Returning to your second question: In homology coordinates, the mass in any given shell increases slightly when mass is added to the star. i.e. for shell of size dq = dm/Mstar, if Mstar gets bigger while dq is fixed, the dm in that shell is larger. However, the mass "added" to this shell is not taken from outside the star. It is effectively taken (advected) from the neighboring shell. More generally, what is happening is that when mass is added to the star all of the q coordinates move slightly in Lagrangian (=mass) coordinates, and so when one writes down the physics, which is typically phrased in terms of Lagrangian derivatives, one must account for the relative motion of the q coordinates. I've attached the slide from my lecture <http://mesastar.org/teaching-materials/2014-mesa-summer-school-working-dir/townsley/townsley_acc_WD_1.pdf/view> from the mesa summerschool this summer that tries to convey what is happening to the mesh in preparation for computing the next timestep. The movement of the q coordinates appears as an advection term in the gravi-thermal energy equation, and is called the "homologous" part in Sugimoto & Nomoto (1975PASJ...27..197S <http://adsabs.harvard.edu/abs/1975PASJ...27..197S>) (equation 10). This still requires a boundary condition to determine what is "neighboring" to the surface. The usual choice is that material arrives there at the same entropy as is already there, which is naturally implemented by just allowing material to advect from the edge of the grid and not doing anything special there.
>>>
>>> I should place a word of caution that as a result of MESA star favoring Lagrangian coordinates, the advection terms needed for species gradients may have some issues. We'd be happy to have help working out the details if this is important for your problem. In my simulations I currently avoid this issue by choosing timesteps in such a way that the constant-q ("Eulerian") region of the mesh is confined to uniform abundance regions.
>>>
>>>
>>> Hopefully that is helpful.
>>>
>>> Dean
>>>
>>>
>>>
>>>
>>>
>>> On 11/09/2014 10:04 PM, Masanobu Kunitomo wrote:
>>>> Hi Dean,
>>>>
>>>> I appreciate that you gave me the comments so quickly!
>>>>
>>>>
>>>> Thank you for the suggestion. I will update to ver 7184 and use the transition region set by eps_grav_time_deriv_separation. In ver. 6596, I set that the transition region is ten cells or more (namely, “min_cells_for_Eulerian_to_Lagrangian_transition = 10”).
>>>>
>>>> I am sorry for the misleading sentence. I can change the Eulerian region with “min_dxm_Eulerian_div_dxm_added” even in ver. 6596, but sometimes the Eulerian region disappears when the timestep (and then dm) is very small. I wanted to maintain a large Eulerian region even in such a case.
>>>>
>>>> I am glad to have your comments on the CpTMdot controls! I wondered why CpTMdot/L is still used to limit the Eulerian region. Since in our study the accretion is highly adiabatic, I set “min_dxm_Eulerian_div_dxm_CpTMdot_lt_L = -1”.
>>>>
>>>> As for my second question, I am still not sure the two assumptions are physically valid: (1) a parcel of accreting material can instantaneously enter the stellar interior and (2) its entropy is the same as that of the stellar interior. I think your “compressional heating” in Townsley and Bildsten (2004) considers the heat leak of the accreting material, but here we do not take into account such a energy release...
>>>>
>>>>
>>>> Best wishes,
>>>> Masanobu
>>>>
>>>>
>>>>
>>>>> On Nov 10, 2014, at 11:27, Dean Townsley <Dean.M.Townsley at ua.edu <mailto:Dean.M.Townsley at ua.edu>> wrote:
>>>>>
>>>>> Hi Masanobu,
>>>>>
>>>>> Nice to have you using this stuff! Hopefully I can help you a bit with your questions.
>>>>>
>>>>> First, if it's not too much trouble, you might want to use the most recent release (7184). This release contains some changes to how the boundaries for the Eulerian vs. Lagrangian forms for eps_grav are determined, and some updates for consistency with how the mesh behaves (constant dq or constant dm) in these regions. With this new version min_dxm_Eulerian_div_dxm_added should work. Please let me know if it doesn't do what you expect and I can investigate.
>>>>>
>>>>> I think the CpTMdot controls aren't being used anymore because either the "Eulerian" or Lagrangian form of eps_grav is valid in the whole star so that it doesn't make sense to limit on this quantity anymore. I still need to clear this control out cleanly, sorry about that. The control used to adjust for accuracy (or validity) is now eps_grav_time_deriv_separation, but you shouldn't need to adjust it if you push the transition region far into the star anyway. I reproduce the comment about it below, since it may be that you don't need to push the Eulerian region over the whole star and you can just let this control adjust to the mass added each time step for you.
>>>>>
>>>>>
>>>>> !### eps_grav_time_deriv_separation
>>>>>
>>>>> ! Separation (in grid cells) over which eps_grav can be time-differenced when Mstar changes
>>>>> ! The mesh has two major regions - an interior region where the cells are
>>>>> ! Lagrangian and an outer region where they are homologous (constant dq =dm/M).
>>>>> ! There is also a small transition region between these two. In the Lagrangian
>>>>> ! region Ds/dt is evaluated with a Lagrangian finite difference in time, while in the
>>>>> ! homologous region Ds/dt is evaluated with a finite difference in time at constant q
>>>>> ! plus an advection-like term accounting for the movement of the q boundaries in mass.
>>>>> ! In the transition region, these two derivatives are combined. This means that at
>>>>> ! the edges of the transition region, finite differences may cross cell boundaries.
>>>>> ! This control determines how the mesh for the end of the current timestep is
>>>>> ! placed to ensure that finite differences cross no more than this many cell
>>>>> ! boundaries.
>>>>>
>>>>> eps_grav_time_deriv_separation = 1.5
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> On your second question: The new material is added at the outside of the star - not throughout the star. The size of the mesh zones change in the "Eulerian" or "constant q" part of the star (which may be the whole star if you push the boundary that deep), but the entropy and its derivatives are evaluated locally, not with respect to the surface. So a parcel of material enters the star at the surface with the same entropy as the surface, but its entropy then evolves appropriately as it moves down through the star.
>>>>>
>>>>>
>>>>> Hopefully that helps. Like I said, let me know if the newer version doesn't do what you expect.
>>>>>
>>>>> Dean
>>>>>
>>>>>
>>>>> On 11/09/2014 07:06 PM, Masanobu Kunitomo wrote:
>>>>>> Dear all,
>>>>>>
>>>>>>
>>>>>> I would like to ask two questions about (1) changing the region of “Eulerian scheme” (s% k_below_Eulerian_eps_grav) and (2) the validity of the Eulerian scheme. I am using ver. 6596 and working on a protostar’s evolution with intense accretion (the mass accretion rate is 1e-5Msun/yr and CpTMdot/L is at most ~100).
>>>>>>
>>>>>>
>>>>>> (1)
>>>>>> I would like to try the case that eps_grav is calculated by the “Eulerian” scheme in the entire star, namely s% k_below_Eulerian_eps_grav = s% nz (although this may not be good because the Eulerian scheme can cause a numerical diffusion).
>>>>>> I tried to set a large “min_dxm_Eulerian_div_dxm_added” in the inlist file, but this did not work. Next I tried to modify the subroutine "set_Eulerian_Lagrangian_for_eps_grav” in “mesa/star/private/evolve.f”. I tried some modifications such as “dxm_kA = s% xmstar”, but it did not work. (In fact, even if I write “stop”, it was not reflected...) Is it inappropriate to modify the evolve.f file?
>>>>>>
>>>>>> (2)
>>>>>> If my understanding is correct, in the Eulerian scheme, the newly accreted mass is distributed in the entire star uniformly and instantaneously, and its entropy is assumed to match the local value of the stellar interior as shown in the Fig. 6 of Sugimoto et al. 1981 (I have attached this figure). If the star is radiative, the accreting material’s entropy can be much lower than that of the stellar surface. Thus, I think this underestimates the accreting material’s entropy. Is this correct?
>>>>>>
>>>>>>
>>>>>>
>>>>>> I apologise if these are already discussed in this mailing list.
>>>>>> I appreciate any comments. Thanks in advance!
>>>>>>
>>>>>> Best wishes,
>>>>>> Masanobu
>>>>>>
>>>>>> --
>>>>>> Masanobu Kunitomo
>>>>>> Ph.D. student / JSPS Research Fellow
>>>>>> Tokyo Institute of Technology
>>>>>> Dept. of Earth and Planetary Sciences
>>>>>> Email: kunitomo at geo.titech.ac.jp <mailto:kunitomo at geo.titech.ac.jp>
>>>>>>
>>>>>>
>>>>>>
>>>>>> ------------------------------------------------------------------------------
>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> mesa-users mailing list
>>>>>> mesa-users at lists.sourceforge.net <mailto:mesa-users at lists.sourceforge.net>
>>>>>> https://lists.sourceforge.net/lists/listinfo/mesa-users <https://lists.sourceforge.net/lists/listinfo/mesa-users>
>>>>>
>>>>
>>>
>>> <mesa_mesh.pdf>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20141116/d92d369f/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: run_star_extras.f
Type: application/octet-stream
Size: 14039 bytes
Desc: not available
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20141116/d92d369f/attachment.obj>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20141116/d92d369f/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: t-q.pdf
Type: application/pdf
Size: 23692 bytes
Desc: not available
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20141116/d92d369f/attachment.pdf>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20141116/d92d369f/attachment-0002.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: t-q2.pdf
Type: application/pdf
Size: 23894 bytes
Desc: not available
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20141116/d92d369f/attachment-0001.pdf>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20141116/d92d369f/attachment-0003.html>
More information about the Mesa-users
mailing list