[Mesa-users] Stability Problems with use_other_energy

Fri Mar 16 11:59:46 EDT 2018

Hi Jose, as a toy example, consider a heating term that depends on the
second derivative of density with mass (don't take it seriously, just an
example):

eps_extra(k) = C*d^2 rho/ dm^2 =
C*(rho(k-1)-2*rho(k)+rho(k+1))/(m(k-1)-m(k+1))

where C is a constant and the last part is just the finite difference
expression of the derivative. Note that this uses the MESA definition where
the surface of the star is at k=1. The partials need to provide information
on how eps_extra changes when different properties change in this and
neighboring cells. The meaning of the 00, m1 and p1 terms is as follows
(the lnd in these terms means natural logarithm of density):

- d_extra_heat_dlndm1(k) is the partial derivative of extra_heat(k) with
respect to ln(rho) at cell (k-1)
- d_extra_heat_dlnd00(k) is the partial derivative of extra_heat(k) with
respect to ln(rho) at cell (k)
- d_extra_heat_dlndp1(k) is the partial derivative of extra_heat(k) with
respect to ln(rho) at cell (k+1)

So, using that d/d(ln rho)=rho*d/d(rho), you can get expressions for the
partial derivatives from the expression for eps_extra. For example:

d_extra_heat_dlndm1(k)= d(eps_extra(k))/d(ln rho(k-1)) = C*
rho(k-1)/(m(k-1)-m(k+1))
d_extra_heat_dlnd00(k)= d(eps_extra(k))/d(ln rho(k)) = -2*C*
rho(k)/(m(k-1)-m(k+1))
d_extra_heat_dlndp1(k)= d(eps_extra(k))/d(ln rho(k+1)) = C*
rho(k+1)/(m(k-1)-m(k+1))

lnT and lnR for the other variables is, as you can guess, natural logarithm
of the temperature and radial coordinates.

Hope this clears things.

Cheers

On Fri, Mar 16, 2018 at 10:40 AM, José Lopes <josevlopes at tecnico.ulisboa.pt>
wrote:

> Hello Pablo,
>
> As a follow-up to your e-mail:
>
> After inspecting the other_energy_implicit.f90 and the required partial
> derivatives for the implicit energy computation, I would like to ask if you
> could identify the variables for the derivatives (e.g.: dlnr00, dlnrm1,
> dlndm1, etc.), or point me to a reference/file where that identification is
> explained, since I am having difficulties in understanding what each
> variable means.
>
> Thank you,
>
> José Vargas Lopes
> PhD Student
>
> CENTRA, Centre for Astrophysics and Gravitation
> IST, University of Lisbon, Portugal
>
> Ligado 16 de março de 2018 no 12:15:37, José Lopes (
> josevlopes at tecnico.ulisboa.pt) escreveu:
>
> Hello Pablo,
>
> Thank you very much for you suggestion!
>
> I will change the implementation and report back once I have some results.
>
> Cheers,
>
> José Vargas Lopes
> PhD Student
>
> CENTRA, Centre for Astrophysics and Gravitation
> IST, University of Lisbon, Portugal
>
>
> No dia 16/03/2018, às 02:39, Pablo Marchant <pamarca at gmail.com> escreveu:
>
> Hi Jose,
>
> my guess is that your implementation needs to be done through the implicit
> other_energy option described in $MESA_DIR/star/other/other_energy_implicit.f90.
> other_energy is only called once at the beginning of the step (an explicit
> calculation), and this can induce large step-to-step variations. The
> implicit routines compute the energy through the Newton iterations and can
> prevent these large variations. One detail though is that you need to
> compute as well partial derivatives. You can see instructions for this in
> the source file I mentioned above.
>
> Cheers
>
> On Thu, Mar 15, 2018 at 9:24 PM, José Lopes via Mesa-users <
> mesa-users at lists.mesastar.org> wrote:
>
>> Dear MESA users,
>>
>> I am developing an extension to MESA to include a Dark Matter Capture
>> framework and interplay with the code. This is all done through the
>> run_star_extras.f and the ‘use_other_energy’ hook. An extra energy
>> contribution is computed for each step (this energy term depends on the
>> stellar structure properties, such has the temperature profile and its
>> second derivative w.r.t. the radius) and is injected as an extra heat term
>> (it can be either positive or negative).
>>
>> So far most results are consistent with what has already been studied for
>> this phenomenology. The problems however start when the extra energy term
>> grows to the same order of magnitude of either the nuclear energy or
>> gravity energy terms. In these cases, the code shows an irregular
>> oscillatory behavior, almost as if it is unable to converge to a stable
>> solution.
>>
>> This oscillatory behavior can be seen in the plots attached to this
>> e-mail (plot1 and plot2). There, you can see the nuclear (black), gravity
>> (red) and extra (blue) energy terms in the first zone (r ~ 0) of the Star,
>> as a function of Model Number. Plot1 shows a case where the extra energy
>> term grows to a magnitude comparable to the nuclear energy term, and plot2
>> shows a similar situation for the gravity term. I also attach (plot3) a
>> case where the extra energy term is smaller and thus the model is able to
>> converge throughout the evolution, although a glitch can be seen in the
>> gravity
>> energy term.
>>
>> I have tried tuning some options in the controls inlist (such has time
>> step control, variable control, mesh-refinement, etc.). However, so far, I
>> have not been able to improve the stability of these models.
>>
>> As anyone encountered a similar behavior? Is there any control option I
>> should focus on to try and improve stability?
>>
>> Thank you all,
>>
>> José Vargas Lopes
>> PhD Student
>>
>> CENTRA, Centre for Astrophysics and Gravitation
>> IST, University of Lisbon, Portugal
>>
>> _______________________________________________
>> mesa-users at lists.mesastar.org
>> https://lists.mesastar.org/mailman/listinfo/mesa-users
>>
>>
>>
>
>
> --
> Pablo Marchant Campos
> M.Sc on Astrophysics, Universidad Católica de Chile
> PhD on Astrophysics, Argelander-Institut für Astronomie, Universität Bonn
>
>

-- 
Pablo Marchant Campos
M.Sc on Astrophysics, Universidad Católica de Chile
PhD on Astrophysics, Argelander-Institut für Astronomie, Universität Bonn
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20180316/7e354b7b/attachment.html>