[FLASH-USERS] Gravity, GPOL_VAR, and GPOT_VAR

[Klaus Weide]

Cole,

Thanks for your message. I am responding here only to the second part:

On Wed, 5 Dec 2007, Cole Miller wrote:

> 2. On a much more complicated note, for a long time I was unsuccessful
>   in adding constant gravity to this problem.  I did all the obvious
>   things, e.g., adding REQUIRES physics/Gravity/ to the Config file,
>   along with gconst and gdirec, making sure that there was a line
>   "useGravity = .true." in flash.par, running it in three dimensions
>   instead  of two, and so on.  No luck.  Eventually, with the help of
>   Sean O'Neill here at Maryland, I tracked down the problem.
>
>   In source/physics/Hydro/HydroMain/split/PPM/PPMKernel/
>   there is a program called hydro_1d.F90.  This program reads the
>   gravity and applies the force, if any, to the fluid.  In it, we have the
>   following lines:
>   #ifdef GPOL_VAR
>   	  call Gravity_accelOneRow (pos, xyzswp, blockID, numIntCells8,ograv,GPOL_VAR)
>   #ifndef GPOT_VAR
>   	  call Driver_abortFlash("Shouldn't have gpol defined without gpot")
>   #endif
>   #endif
>   #ifdef GPOT_VAR
>   	  call Gravity_accelOneRow (pos, xyzswp, blockID, numIntCells8,grav,GPOT_VAR)
>   #endif
>   	  dtfac = dt/dt_old
>
>   	  do i = 1,numIntCells8
>   	     dg       = dtfac*(grav(i) - ograv(i))
>   	     hgrav(i) = grav(i) + 0.5e0*dg
>   	     ngrav(i) = grav(i) +	dg
>   	  enddo
>
>   and then hgrav is passed to the Riemann solver as a force.

I am sorry you had to spend time on this.  We have also found after
the beta release the code did not work properly with constant gravity.

In the current FLASH3 development code, the relevant section in
hydro_1d.F90 loos like this:

    ................................................................................
     if (useGravity) then
        pos(1)=jCell; pos(2)=kCell
#ifdef GPOL_VAR
        call Gravity_accelOneRow (pos, xyzswp, blockID, numIntCells8,ograv,GPOL_VAR)
#ifndef GPOT_VAR
        call Driver_abortFlash("Shouldn't have gpol defined without gpot")
#endif
#endif

#ifdef GPOT_VAR
        call Gravity_accelOneRow (pos, xyzswp, blockID, numIntCells8,grav,GPOT_VAR)
        dtfac = dt/dt_old

        do i = 1,numIntCells8
           dg       = dtfac*(grav(i) - ograv(i))
           hgrav(i) = grav(i) + 0.5e0*dg
           ngrav(i) = grav(i) +       dg
        enddo

#else
        ! Neither GPOT_VAR nor GPOL_VAR defined -> assume time-independent gravity field - KW
        call Gravity_accelOneRow (pos, xyzswp, blockID, numIntCells8,grav)
        hgrav = grav
        ngrav = grav
#endif

     else

        ograv(:) = 0.e0
        hgrav(:) = 0.e0
        ngrav(:) = 0.e0

     end if
    ................................................................................

This should esssentially do the same as your fix (while saving soem
unnecessary operations).

>   The problem is that nowhere that I saw is GPOL_VAR or GPOT_VAR defined.
>   Not in Flash.h, not in any of the files after gmake is run.  I tried
>   several variants of the setup script, with the same result.  Therefore,
>   since prior to the lines quoted above the arrays grav() and ograv() are
>   both set to zero, hgrav is as well, so there is no gravity.

GPOT_VAR should be defined (automatically, in Flash.h) if and only if
there is a "VARIABLE gpot" declaration in a Config file that is
included in the setup. That is currently only the case if one uses
a "Poisson" implementation of gravity, where a variable in UNK to hold
the gravitational potential is actually needed.


Klaus Weide, FLASH Center Code Group