[Mesa-users] convergence problem

[Pablo Marchant]

Hi Sheyda,

The pressure scale height is set in star/private/mlt_info.f90, assuming
standard hydrostatic equilibrium, so its set as

scale_height = P / (grav*rho)

If you want to modify the pressure scale height then you'll have to adjust
this accordingly. You can reimplement mlt using the other_mlt hook, just
copy all what's neccesary from mlt_info into your run_star_extras and
adjust the line above. Note that you also have to set properly the partial
derivatives for the scale height, which are required by the newton solver,
this is currently done in mlt_info in this way

d_scale_height_dvb = scale_height*(dP_dvb/P - d_grav_dvb/grav -
dRho_dvb/Rho)

Here the variables denoted as *_dvb are arrays containing the partial
derivatives with respect to each variable. You should be able to figure out
how itworks by digging a bit in mlt_info. From your modifications, you'll
need to figure out how to get the partial derivatives for dphi with respext
to the structure variables.

Cheers

On Tue, May 14, 2019, 07:06 sheyda najafi via Mesa-users <
mesa-users at lists.mesastar.org> wrote:

> Dear MESA users
>
> I have changed extra_dpdm to add a new acceleration term to HSE
> equation in MESA code.
> (p(i-1)- p(i))/dm_bar(i)=( -standard_cgrav*m(r))/(4*pi* s% r(i)**4) -
> beta*(clight**2) *rho_face(i) (dphi/dm_bar)
> therefore
> extra_dpdm(i)= - beta*(clight**2) *rho_face(i) dphi/dm_bar
>
> which i have asked about the quantities and the shape of this new term
> in previous emails.
>
>  My question is how does this extra acceleration change the pressure
> scale height in the code. I can not find the relation. I know that
> pressure scale height depends on p_face and rho_face and grav. The
> extra_dpdm has impact on density and pressure but I can not find the
> relation between grav and extra_dpdm in the code.
>
> Thank you very much
> Sheyda Najafi
>
> On 5/7/19, sheyda najafi <sh.najafi8590 at gmail.com> wrote:
> > Dear Mesa users
> >
> > I have another question regarding my previous emails.
> >
> > As i  said I wanted to study the evolution of 1 solar mass star in the
> > presence of
> > new extra acceleration term :(-beta*(clight**2)*(dphi/dr)).
> >  i have changed the
> > extra_dPdm subroutine of MESA code to add
> > the extra term (-beta*rho(clight**2)*(dphi/dm)) ==>( in lagrangian
> > framework) .
> >  where phi depends on local density of each layer(cell density: s% rho).
> >  which is named (xn) in the code. and beta is a constant.
> > i have added the extra term as follows
> >
> >            do k = 1,2
> >              s% extra_dPdm(k) = -( ((   xn(k)-    xn(k+1))/(s%
> > dm_bar(k+1)))*((beta*(clight**(2.0d0)))*(s% rho_face(k))) )
> >
> >          end do
> >
> >          do k = 2, s% nz-1
> >           s% extra_dPdm(k) = -( (((   xn(k)-    xn(k+1))/(s%
> dm_bar(k+1)))+
> > &
> >          ((   xn(k-1) -   xn(k))/(s% dm_bar(k))))
> > *((beta*(clight**(2.0d0)))*(s% rho_face(k))*(1.d0/2.d0)) )
> >
> >          end do
> >
> >          do k = s% nz-1, s% nz
> >           s% extra_dPdm(k) = -( ((   xn(k-1) -    xn(k))/(s%
> > dm_bar(k)))*((beta*(clight**(2.0d0)))*(s% rho_face(k))) )
> >
> >         end do
> >
> > I noticed that  s% dm_bar(s%nz) = s%dm(s% nz) + s%dm(s%nz-1)/2
> > and not (s% dm(s% nz) + s% dm(s%nz-1))/2.
> >
> > i am not sure that dm_bar must be used in extra_dpdm or  (s% dm(s% nz)
> > + s% dm(s%nz-1))/2.
> >
> > Thank you very much
> > Sheyda Najafi
> >
> >
> > On 10/21/18, sheyda najafi <sh.najafi8590 at gmail.com> wrote:
> >> Dear Josiah and MESA users
> >>
> >>
> >> I have another question regarding my previous emails. i have asked it
> >> before but unfortunately no one answered.
> >>
> >> As i  said I wanted to study the evolution of 1 solar mass star in the
> >> presence of
> >> new extra acceleration term :(-beta*(clight**2)*(dphi/dr)).
> >>  i have changed the
> >> extra_dPdm subroutine of MESA code to add
> >> the extra term (-beta*rho(clight**2)*(dphi/dm)) ==>( in lagrangian
> >> framework) .
> >>  where phi depends on local density of each layer(cell density: s% rho).
> >>  which is named (xn) in the code. and beta is a constant.
> >>
> >>            do k = 1,2
> >>              s% extra_dPdm(k) = -( ((   xn(k)-    xn(k+1))/(s%
> >> dm_bar(k+1)))*((beta*(clight**(2.0d0)))*(s% rho_face(k))) )
> >>
> >>          end do
> >>
> >>          do k = 2, s% nz-1
> >>           s% extra_dPdm(k) = -( (((   xn(k)-    xn(k+1))/(s%
> >> dm_bar(k+1)))+
> >> &
> >>          ((   xn(k-1) -   xn(k))/(s% dm_bar(k))))
> >> *((beta*(clight**(2.0d0)))*(s% rho_face(k))*(1.d0/2.d0)) )
> >>
> >>          end do
> >>
> >>          do k = s% nz-1, s% nz
> >>           s% extra_dPdm(k) = -( ((   xn(k-1) -    xn(k))/(s%
> >> dm_bar(k)))*((beta*(clight**(2.0d0)))*(s% rho_face(k))) )
> >>
> >>         end do
> >>
> >> Since i did not know the value of s%rho(0) i did'nt have xn(0) to
> >> write extra_dpdm it in 1 part therefore i have split it in three
> >> parts. how it is possible to write this term in 1 part. i mean if i
> >> want to write the extra term as :
> >>
> >>          do k = 1, s% nz
> >> s% extra_dpdm(k)= -((   xn(k-1)-    xn(k)) / s% dm_bar(k)) * beta*
> >> (clight**2)* ( s%rho_face(k)).
> >>         end do
> >>
> >> I will have problem in finding the value of xn(0). because xn is
> >> related to cell density s%rho(0) via
> >>
> >>    xn(i)= (( n1*(lambda**(4.d0+n1)) )/(((mpl)**(n1))*clight*beta*s%
> >> rho(i)))**(1.d0/(n1+1.d0))
> >>
> >> lambda, mpl, n1 and beta are constant.
> >>
> >>  which density should i use as rho(0) ??
> >>
> >> (for your kindly information my codes are attached in the email)
> >> Thank you in advance
> >>
> >>
> >> Sheyda
> >>
> >
> _______________________________________________
> mesa-users at lists.mesastar.org
> https://lists.mesastar.org/mailman/listinfo/mesa-users
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20190514/299e3e23/attachment.html>