[FLASH-USERS] Setup of hydrostatic equilibrium in FLASH 3.0beta?

[Cole Miller]

Hi all,

Has anyone set up a good hydrostatic equilibrium routine for FLASH
3.0beta?  I've got something that basically works (in the context of
an X-ray burst problem), but I'm concerned about its robustness and
there is a mysterious factor of 0.25 that I have to put in, so I'd
appreciate greater understanding.  As far as I can tell, the setup
breaks into two tasks: initializing the interior cells, which seems
straightforward, and setting the boundary conditions at the edge,
which is where potential issues are arising.

1. For the interior cells, I simply wrote a little routine in
    Simulation_initBlock.F90 that computes the vertical density structure
    in hydrostatic equilibrium on a fine mesh.  Then, when densities
    are assigned in each block, I obtain the z coordinate at each point
    and do a linear interpolation to find the density.  The routine
    for HSE assumed a user input of density at the minimum z value
    in the simulation (assuming gravity is in the -z direction),
    combined with a temperature (assumed constant throughout), then
    computed the vertical density structure at the next z value up
    using a search by bisection for a density that would yield a
    pressure gradient force that counteracted gravity.  This involved
    calls to the Eos routine, e.g.,

   eosData(EOS_TEMP) = sim_tempAmbient
   dzz = (1.0/(Nhse-1.0))*(sim_zmax-sim_zmin)
   do k = 1, Nhse
     zz = sim_zmin + (k-1)*dzz
     if (k .EQ. 1) then
       density(k) = sim_rhoAmbient
       eosData(EOS_DENS) = density(k)
       call Eos(MODE_DENS_TEMP,1,eosData,massFraction)
       pressure(k)=eosData(EOS_PRES)
     else
       dlogmin=0.0   ! Minimum ln density in cgs
       dlogmax=log(density(k-1))  ! Max.  This is dens of lower cell
       do ibisec=0,30
         dlogmid=0.5*(dlogmin+dlogmax)
         density(k)=exp(dlogmid)
         eosData(EOS_DENS) = density(k)
         call Eos(MODE_DENS_TEMP,1,eosData,massFraction)
         pressure(k)=eosData(EOS_PRES)
!! Now we divide the pressure gradient by the average density.  If the
!! pressure gradient is too small, then the density in the current
!! zone has to be less
         gradient=(pressure(k-1)-pressure(k))/dzz
         gradient=gradient/(0.5*(density(k)+density(k-1)))
         if (gradient<(-1.00*sim_gconst)) then
           dlogmax=dlogmid
         else
           dlogmin=dlogmid
         endif
       enddo

    Here Nhse is a large value (1001 in my case) to make later
    interpolation more accurate.  sim_tempAmbient, sim_rhoAmbient,
    and sim_gconst are input in flash.par.  This seems to work fine.


2. However, when I ran this, although the inner regions seemed in good
    equilibrium, the bottom part of the simulation is not; there is a lot
    of bouncing and thus generation of heat.  Assuming this is because of
    the lower boundary condition, I started looking at
    Grid_applyBCEdge.F90.

    In principle, I'd like to do the same thing here as in
    Simulation_initBlock.F90.  That is, for each of the guard cells,
    I want to be able to figure out the vertical coordinate and then
    do something similar to what I did above.  For example, for the
    guard cells off the low end of the grid, the *upper*
    boundary condition would be a density of sim_rhoAmbient and a
    temperature of sim_tempAmbient.

    I ran into two problems, possibly related to each other:

    (a) I was not able to find a specific routine that would tell me the
        coordinates of the center of a guard cell.  The only routine that
        I saw that gives cell coordinates is Grid_getCellCoords, and one
        can use this to get the coordinates of guard cells.  This,
        however, requires a blockId as input.  I've been using
        blockId=1 and it seems to work for the simple cases I've
        tried, but I'm worried that it will fail after refinement.
        Is Grid_applyBCEdge.F90 called only once, or for every block
        that borders an edge?  Is there a separate routine that gives
        the coordinates of guard cells?

    (b) When I made this choice, I set up a routine similar to the one
        I reproduced above, but in Grid_applyBCEdge.F90.  It gives me
        good equilibrium conditions, but only if I adopt a strange
        ad hoc factor.  What I seem to need to do is that instead of
            if (gradient<(-1.00*sim_gconst)) then
        I have to use
            if (gradient<(-0.25*sim_gconst)) then

        This factor is independent of the vertical resolution for the
        cases I've tried (of course, the equilibrium is better for better
        vertical resolution).

As I say, doing all this does give very nice hydrostatic equilibrium,
but until I understand these issues I'm leery of using this for any
serious runs.  Any idea what's happening?  Does anyone have a routine
for FLASH 3.0beta?  I see that FLASH2.5 had a hydrostatic equilibrium
routine, but I didn't see anything there that addressed guard cells.

Thanks!
Cole