[FLASH-USERS] using the split PPM hydro solver in 3D spherical polar coordinates

[Rodrigo Fernandez]

Dear FLASH Users,

I'm writing to document the changes to the public FLASH3 distribution
needed to enable the split PPM hydro solver to work in 3D spherical polar
coordinates. I've been told that a number of people are interested in this
functionality. Since I've done some debugging and testing of this myself, I
thought that the results would be useful to those interested.

Best,

Rodrigo

******************************************

The split PPM solver has most of the functionality for operating in 3D
spherical coordinates already built in. Only minor changes are required.
Below I use line numbers from FLASH3.2, although most of the changes are
also applicable to FLASH4, which has almost the same routines. The relevant
files are in

source/physics/Hydro/HydroMain/split/PPM

Steps:

1) In Hydro_detectShock.F90, comment out the call to Driver_abortFlash for
3D spherical coordinates (line 131)

2) In the same subroutine, add

     else if (hy_geometry == SPHERICAL) then

     Before the endif in line 308, after the block labeled "3-d cylindrical
geometry". Then copy the block directly above which computes the velocity
divergence in 3-d cylindrical coordinates, and modify it for 3D spherical
coordinates.

3) In PPMKernel/hydro_1d.F90, add an extra factor of 1/sin(theta) to dtdx
to correct the inverse line element for the azimuthal coordinate. After the
end of the "if" block that starts with "if (igeom >=PHY_CYL)" in line 361,
add

    if (igeom == PHI_SPH) then
      do i = 1, numIntCells8
        dtdx(i) = dtdx(i) / sin(thirdCoord(kCell))

      enddo
    endif

4) If you ever need to use the Hybrid Riemann solver (e.g. for shocks
parallel the grid), a correction needs to be made in hy_ppm_sweep.F90.
After "if (hy_hybridRiemann) then" in line 401, comment out:

    call Hydro_detectShock(...)

    and replace with:

    if (sweepDir == SWEEP_X) then
      call Hydro_detectShock(solnData, shock,blkLimits,blkLimitsGC,&

 primaryCoord,secondCoord,thirdCoord)

    else if (sweepDir == SWEEP_Y) then
      call Hydro_detectShock(solnData, shock,blkLimits,blkLimitsGC,&

 secondCoord,primaryCoord,thirdCoord)

     else
       call Hydro_detectShock(solnData, shock,blkLimits,blkLimitsGC,&

secondCoord,thirdCoord,primaryCoord)

     endif

    FLASH4 is different here, so you'll need to check whether this is still
necessary.

5) If you want to use the full range of polar angles, you also need to put
a safeguard in PPMKernel/avisco.F90.
     In line 416 (xyzswp == SWEEP_Y section of else if block below "(r_sph,
theta, phi)"), replace

    sinth   = sin( xl(i) )

    with

    sinth   = max(sin( xl(i) ),1e-6)

    Or with any other number that is small compared to PI. This factor
divides the phi contribution to the velocity divergence, and is zero if the
face coordinate xl = 0 or PI.


6) In addition, users should be aware that the split PPM solver
reconstructs in the coordinate, not in volume. This is OK for coordinates
with invariant transverse area (e.g. cartesian), but not for others like
spherical radius. Nevertheless, at high enough resolution the correct
behavior is obtained asymptotically, e.g. for spherical radius:

      R_R^3 - R_L^3 = 3*R_C^2 * DR * ( 1 + [DR/R_C]^2/12 )

Where R_R, R_C, and R_L refer to right, centered, and left radial
coordinate, respectively, and DR = R_R - R_L.

Users should make sure that the resolution is good enough for the code to
converge to the right answer within some tolerance, by e.g. carrying out
tests with known solutions.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://flash.rochester.edu/pipermail/flash-users/attachments/20141222/6cacbce5/attachment.htm>