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

[Klaus Weide]

Cole,

Sorry for responding late; we have been busy with preparing for the 
FLASH3.0 release.

On Fri, 18 Jan 2008, Cole Miller wrote:

> 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

Note that you may be able to change this from a bisection method to 
Newton-Raphson, by taking advantage of the derivatives (like 
eosData(EOS_DPD) ) that your Eos call should already compute.

>   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?

You are right, there is no clean way to get at the cell coordinates 
within BC handling in FLASH 3.0beta. We are going to provide ways to do 
this in FLASH3.0, to be released soon.

Even if you could get the "blockId" in Grid_applyBCEdge.F90 there would
be possible complications, because Paramesh may sometimes call this
routine on behave of a data block that belongs to a different PE
and thus does not, strictly speaking, have a valid blockId on the
local PE.  So you would also need a modified Grid_getCellCoords that
can deal with this kind of situation. (These difficulties are 
being addressed in FLASH3.0.)

Your trick of always using blockId=1 should of course only work if the 
cells you are trying to fill are cells of block #1. Or, if you just need 
to get the cell spacing (in which case you should just use 
Grid_getDeltas), it should work iff the block you are working on has the 
same resolution as block #1 -- which normally is not the case for any 
other block, since #1 is the root block and a parent or ancestor of all 
others (unless you use an NblockX / NblockY / NblockZ > 1).

>   (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).

But did the block whose guard cells you are fillling have the same
resolution as your root block, block #1?  If not, maybe that is
how the factor comes in: you resolution is off by a factor of 4?

To understand this better, we would have to know more about your
setup. Especially the number of blocks, lrefine_{min,max}, etc.

I hope this helps,
Klaus