[mesa-users] questions about the "Eulerian" scheme for eps_grav
Dean Townsley
Dean.M.Townsley at ua.edu
Thu Nov 13 01:13:09 EST 2014
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
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 \93q coordinate region\94 instead of \93Eulerian
> region\94.
>
>
> I tried the case with "eps_grav_time_deriv_separation = -1\94 with ver.
> 7184, and I succeeded in making the entire star the \93q coordinate\94
> 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\92s 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\85 I attached a figure which shows the evolution of the q
> coordinate region with time. At the end, \93q_below_just_added\94 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
> \93set_Eulerian_Lagrangian_for_eps_grav\94 in the file
> \93mesa/star/private/evolve.f\94. I tried some modifications such as
> \93dxm_kA = s% xmstar\94 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,
>>> \93min_cells_for_Eulerian_to_Lagrangian_transition = 10\94).
>>>
>>> I am sorry for the misleading sentence. I can change the Eulerian
>>> region with \93min_dxm_Eulerian_div_dxm_added\94 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
>>> \93min_dxm_Eulerian_div_dxm_CpTMdot_lt_L = -1\94.
>>>
>>> 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
>>> \93compressional heating\94 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
>>>>> \93Eulerian scheme\94 (s% k_below_Eulerian_eps_grav) and (2) the
>>>>> validity of the Eulerian scheme. I am using ver. 6596 and working
>>>>> on a protostar\92s 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
>>>>> \93Eulerian\94 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 \93min_dxm_Eulerian_div_dxm_added\94 in the
>>>>> inlist file, but this did not work. Next I tried to modify the
>>>>> subroutine "set_Eulerian_Lagrangian_for_eps_grav\94 in
>>>>> \93mesa/star/private/evolve.f\94. I tried some modifications such as
>>>>> \93dxm_kA = s% xmstar\94, but it did not work. (In fact, even if I
>>>>> write \93stop\94, 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\92s entropy can be much lower than
>>>>> that of the stellar surface. Thus, I think this underestimates the
>>>>> accreting material\92s 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
>>>>
>>>
>>
>> <mesa_mesh.pdf>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20141113/ffd00f43/attachment.html>
More information about the Mesa-users
mailing list