[mesa-users] questions about the "Eulerian" scheme for eps_grav

Tue Nov 11 15:49:29 EST 2014

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> 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/20141112/487bb4ba/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: t-q.pdf
Type: application/pdf
Size: 29730 bytes
Desc: not available
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20141112/487bb4ba/attachment.pdf>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20141112/487bb4ba/attachment-0001.html>