# [mesa-users] Brunt Vaisala Frequency

Bill Paxton paxton at kitp.ucsb.edu
Sat Feb 26 15:22:31 EST 2011

```On Feb 26, 2011, at 11:52 AM, Jim Fuller wrote:

> Hey everybody,
>
> I've been using mesa to calculate the Brunt-Vaisala frequency in main sequence stars of various masses. Since the brunt freq. is proportional to the force of gravity, it should go to zero at the center of the star, regardless of mass. However, in the models I am generating, it approaches some finite positive value near the center of my stars (for stars of various masses and ages). Has anybody run into this problem or can point me to the location in the code where the brunt frequency is calculated? Thanks,
>

Hi Jim,

Good question -- I wonder if it might be an artifact of smoothing.
Here's how to track down the place in the code where brunt is calculated.

In star/private, do

grep brunt *

among a bunch of output, you'll find that micro.f has what we want.

subroutine do_brunt_N2_for_cell
g = s% grav(k)
r = s% r(k)
csound = interp_val_to_pt(s% csound,k,s% dq)
s% brunt_N2(k) = (g/r)*(-g*r/csound**2 - s% dlnRho_dlnR(k))

the call on interp_val_to_pt is because mesa calculates csound at cell centers,
but we need a value at the cell boundary for use in Brunt.
The other things, g, r, and dlnRho_dlnR are defined at cell boundaries already.

we next need to see how dlnRho_dlnR is calculated, so back to grep where we find
subroutine do_brunt_N2
dlnd(:) = s% lnd(:)
dlnr(:) = s% lnr(:)
call make_delta(dlnd,s% brunt_nsmooth,nz,.false.)
call make_delta(dlnr,s% brunt_nsmooth,nz,.false.)
s% dlnRho_dlnR(:) = dlnd(:)/dlnr(:)

the routine make_delta is defined in the same place.
it smooths the "raw" values of lnd and lnr before taking differences.
the amount of smoothing is controlled by 'brunt_nsmooth'

checking star_defaults.dek shows that the default for brunt_nsmooth is 50,
so there is quite a bit of smoothing happening.

You might want to try smaller values for brunt_nsmooth and see what happens.

Alternatively, you might want a different smoothing scheme --
the one in the code now is very crude.  If you find something that
works well, please let me know and we'll put it into mesa!

Hope that helps.

Cheers,
Bill

p.s., I also see that there is a bug in the calculation of dlnRho_dlnR!
The code in do_brunt_N2 should actually be using cell centered lnr.
I'll fix that for the next release.   If you'd like to try changing it sooner,
here's a modified version of the routine -- just replace the current
definition in micro.f and then do

cd star/test; ./mk ; ./ck ; ./export

subroutine do_brunt_N2(s,nzlo,nzhi,ierr)
type (star_info), pointer :: s
integer, intent(in) :: nzlo, nzhi
integer, intent(out) :: ierr
logical, parameter :: use_omp = .true.
integer :: nz
double precision, pointer :: dlnd(:), dlnr(:)
ierr = 0
if (dbg) write(*,*) 'do_brunt_N2 call foreach_cell', nzlo, nzhi
nz = s% nz
allocate(dlnd(nz), dlnr(nz))
dlnd(:) = s% lnd(:)
dlnr(:) = s% lnr(:)
call make_delta_lnd(dlnd,s% brunt_nsmooth,nz)
call make_delta_lnr(dlnr,s% brunt_nsmooth,nz)
s% dlnRho_dlnR(:) = dlnd(:)/dlnr(:)
deallocate(dlnd,dlnr)
call make_smooth(s% dlnRho_dlnR,s% brunt_nsmooth,nz)
call foreach_cell(s,nzlo,nzhi,use_omp,do_brunt_N2_for_cell,ierr)

contains

subroutine make_smooth(z,nsmooth,nz)
double precision, pointer :: z(:)
integer, intent(in) :: nz, nsmooth
integer :: j, k
do j=1,nsmooth
z(1) = (3*z(1) + z(2))*0.25d0
do k=2,nz-1
z(k) = (z(k-1) + 2*z(k) + z(k+1))*0.25d0
end do
z(nz) = (z(nz-1) + 3*z(nz))*0.25d0
do k=nz-1,2,-1
z(k) = (z(k-1) + 2*z(k) + z(k+1))*0.25d0
end do
end do
end subroutine make_smooth

subroutine make_delta_lnd(z,nsmooth,nz)
double precision, pointer :: z(:)
integer, intent(in) :: nz, nsmooth
integer :: j, k
call make_smooth(z,nsmooth,nz)
do k=2,nz
z(k-1) = z(k-1) - z(k)
end do
z(nz) = z(nz-1)
end subroutine make_delta_lnd

subroutine make_delta_lnr(z,nsmooth,nz)
double precision, pointer :: z(:)
integer, intent(in) :: nz, nsmooth
integer :: j, k
call make_smooth(z,nsmooth,nz)
do k=2,nz-1
z(k-1) = 0.5d0*(z(k-1) - z(k+1))
end do
z(nz) = z(nz-1)
z(1) = z(2)
end subroutine make_delta_lnr

end subroutine do_brunt_N2

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20110226/bbfc20cd/attachment.html>
```