Re: [mesa-users] Brunt Vaisala Frequency

[Bill P. <pa...@ki...>]

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