<div dir="ltr">Dear FLASH Users,<div><br></div><div>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. </div><div><br></div><div>Best,</div><div><br></div><div>Rodrigo</div><div><br></div><div>******************************************</div><div><br></div><div>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 </div><div><br></div><div>source/physics/Hydro/HydroMain/split/PPM</div><div><br></div><div>Steps:<br></div><div><br></div><div>1) In Hydro_detectShock.F90, comment out the call to Driver_abortFlash for 3D spherical coordinates (line 131)</div><div><br></div><div>2) In the same subroutine, add</div><div><br></div><div>     else if (hy_geometry == SPHERICAL) then</div><div><br></div><div>     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.</div><div><br></div><div>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</div><div><br></div><div>    if (igeom == PHI_SPH) then</div><div>      do i = 1, numIntCells8</div><div>        dtdx(i) = dtdx(i) / sin(thirdCoord(kCell))</div><div><br></div><div>      enddo</div><div>    endif</div><div><br></div><div>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:</div><div><br></div><div>    call Hydro_detectShock(...)</div><div><br></div><div>    and replace with:</div><div><br></div><div><div>    if (sweepDir == SWEEP_X) then</div><div>      call Hydro_detectShock(solnData, shock,blkLimits,blkLimitsGC,&</div><div>                                             primaryCoord,secondCoord,thirdCoord)</div><div><br></div><div>    else if (sweepDir == SWEEP_Y) then</div><div>      call Hydro_detectShock(solnData, shock,blkLimits,blkLimitsGC,&</div><div>                                             secondCoord,primaryCoord,thirdCoord)</div><div><br></div><div>     else</div><div>       call Hydro_detectShock(solnData, shock,blkLimits,blkLimitsGC,&</div><div>                                              secondCoord,thirdCoord,primaryCoord)</div><div><br></div><div>     endif</div></div><div><br></div><div>    FLASH4 is different here, so you'll need to check whether this is still necessary.</div><div><br></div><div>5) If you want to use the full range of polar angles, you also need to put a safeguard in PPMKernel/avisco.F90. </div><div>     In line 416 (xyzswp == SWEEP_Y section of else if block below "(r_sph, theta, phi)"), replace</div><div> </div><div>    sinth   = sin( xl(i) )</div><div><br></div><div>    with </div><div><br></div><div>    sinth   = max(sin( xl(i) ),1e-6)</div><div><br></div><div>    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.</div><div><br></div><div><br></div><div>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:</div><div><br></div><div>      R_R^3 - R_L^3 = 3*R_C^2 * DR * ( 1 + [DR/R_C]^2/12 )</div><div><br></div><div>Where R_R, R_C, and R_L refer to right, centered, and left radial coordinate, respectively, and DR = R_R - R_L.</div><div><br></div><div>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. </div><div><br></div><div><br></div><div><br></div></div>