[FLASH-USERS] DustCollapse Weird Artifacts in Gravitational Acceleration

[Joshua Wall]

I'm going to habor a guess here, but try adding a call to the following
after you solve for the potential and before you get the acceleration to:

call Grid_fillGuardCells(CENTER, ALLDIR, doEos=.false., &
            selectBlockType=LEAF,unitReadsMeshDataOnly=.true.)

Don't forget you'll have to add a:

use Gravity_interface, only: Grid_fillGuardCells

to the subroutine before the variable declarations. This will make sure the
guard cells all have the proper potential in them before you attempt to
calculate the acceleration, which uses the guard cells to calculate the
finite difference of the potential.

Cordially,

Josh

On Tue, Feb 28, 2017 at 11:58 PM Ryan Farber <rjfarber at umich.edu> wrote:

> Oops, I didn't notice you had sent this to flash-users when I replied to
> you yesterday.
>
> The BHTree solver with bhtreeAcc=1 does indeed work perfectly (y
> gravitational acceleration attached,
> all plots from this email come from debugging runs done with FLASH 4.4).
>
> However, the "production" version of FLASH I'm working with (heavily
> modified) is a bit too old for me to be
> able to use this fix. Moreover, old work used other solvers in which the
> acceleration is computed from the potential
> in Gravity_accelOneRow (which is done the exact same way independent of
> the solver employed to compute
> the potential, as I've verified), and I'm worried how this bug may affect
> that old work. Also, I just gotta know!!!
> I had a dream about debugging this last night...
>
> Any way, if I use BHTree without the "bhtreeAcc=1" then it computes the
> acceleration through the potential.
> But, in this case Gravity_accelOneRow doesn't get called on the first time
> step. So, I made the
> Driver_initFlash edits as described in my first email to get *accelOneRow
> called on the 1st time step.
> And I get the attached y gravitational acceleration, which looks pretty
> nice but still has a weird jump on the
> boundary.
>
> Do you think the boundary jump is supposed to be there? Or is this a bug?
>
> Best,
> --------
> Ryan
>
> On Mon, Feb 27, 2017 at 1:32 PM, Joshua Wall <joshua.e.wall at gmail.com>
> wrote:
>
> Hello Ryan,
>
>       I double checked this morning, and Richard Wunsch's BHTree gravity
> includes variables for storing the gravitation acceleration at each cell if
> you use his solver with the added option bhtreeAcc at setup like this (make
> sure you take out any references to other gravity solvers in the simulation
> Config file):
>
> ./setup DustCollapse +auto +3d
> --with-unit=physics/Gravity/GravityMain/Poisson/BHTree bhtreeAcc=1
>
> If you check in the code in
> source/physics/Gravity/GravityMain/Poisson/BHTree/Gravity_accelOneRow.F90
> you can see how the accelerations are calculated and stored exactly from
> both the grid and any sinks in the code. I've been using his potential
> solver for some time (in Flash 4.2.2, including getting the accelerations)
> and find it to be both robust and fast. Or perhaps you could use his
> implementation as an example of how to structure your own modifications to
> Gravity_accelOneRow.F90.
>
> Note that if you do use his code it already includes the gravity accel
> variables in source/physics/Gravity/GravityMain/Poisson/BHTree/Config,
> which I encourage you to check out as it describes the other options for
> the code. The paper for this method is in prep still I believe.
>
> Hope this helps!
>
> Josh
>
> On Sat, Feb 25, 2017 at 10:03 PM Ryan Farber <rjfarber at umich.edu> wrote:
>
> Dear FLASH users,
>
> I decided to look at the components of the gravitational acceleration for
> the supplied
> DustCollapse problem (FLASH 4.3), and have noticed some weird artifacts on
> the block boundaries.
>
> The issue seems somewhat similar to one Thomas Peters' refinement issue
> <http://flash.uchicago.edu/pipermail/flash-users/2014-May/001440.html> (a
> plotting bug)
> but I chatted with ngoldbaum on the yt IRC channel and the conclusion was
> the artifacts
> are on the FLASH side of things.
>
> The only edits I made to the DustCollapse problem supplied with FLASH 4.3
> were:
>
> * Config: added the following lines
> VARIABLE grac
> VARIABLE poix
> VARIABLE poiy
> VARIABLE poiz
>
> * Driver_initFlash.F90: copied to problem directory and modified the
> original lines:
>   if(.not. dr_restart) then
>      call Grid_getListOfBlocks(LEAF,blockList,blockCount)
>      call Gravity_potentialListOfBlocks(blockCount,blockList)
>      call Particles_initForces()
>   end if
>
>   ! If we want to free any arrays created during simulation
>   ! initialization that are no longer needed, do it here.
>   call Simulation_freeUserArrays()
>
>   call IO_outputInitial(  dr_nbegin, dr_initialSimTime)
>   if(dr_globalMe==MASTER_PE)print*,'Initial plotfile written'
>
>   if(dr_globalMe==MASTER_PE)print*,'Driver init all done'
>
>
> **to**
>   if(.not. dr_restart) then
>      call Grid_getListOfBlocks(LEAF,blockList,blockCount)
>      call Gravity_potentialListOfBlocks(blockCount,blockList)
>      call Particles_initForces()
>
>      call Gravity_accelListOfBlocks(blockCount, blockList, IAXIS, GRAC_VAR)
>      call Gravity_accelListOfBlocks(blockCount, blockList, JAXIS, GRAC_VAR)
>      call Gravity_accelListOfBlocks(blockCount, blockList, KAXIS, GRAC_VAR)
>   end if
>
>   ! If we want to free any arrays created during simulation
>   ! initialization that are no longer needed, do it here.
>   call Simulation_freeUserArrays()
>
>   call IO_outputInitial(  dr_nbegin, dr_initialSimTime)
>   if(dr_globalMe==MASTER_PE)print*,'Initial plotfile written'
>
>   STOP(99)
>
>   if(dr_globalMe==MASTER_PE)print*,'Driver init all done'
>
> * Gravity_accelOneRow: modified the original lines:
>   do ii = iimin+1, iimax-1
>      grav(ii) = grav(ii) + delxinv * (gpot(ii-1) - gpot(ii+1))
>   enddo
>
>   grav(iimin) = grav(iimin+1)     ! this is invalid data - must not be used
>   grav(iimax) = grav(iimax-1)
>
> **to**
>   do ii = iimin+1, iimax-1
>      grav(ii) = grav(ii) + delxinv * (gpot(ii-1) - gpot(ii+1))
>
>      if (sweepDir .eq. SWEEP_X) then
>        solnVec(POIX_VAR, ii, pos(1), pos(2)) = grav(ii)
>      end if
>      if (sweepDir .eq. SWEEP_Y) then
>        solnVec(POIY_VAR, pos(1), ii, pos(2)) = grav(ii)
>      end if
>      if (sweepDir .eq. SWEEP_Z) then
>        solnVec(POIZ_VAR, pos(1), pos(2), ii) = grav(ii)
>      end if
>   enddo
>
>   grav(iimin) = grav(iimin+1)     ! this is invalid data - must not be used
>   grav(iimax) = grav(iimax-1)
>   if (sweepDir .eq. SWEEP_X) then
>     solnVec(POIX_VAR, iimin, pos(1), pos(2)) = grav(iimin+1)
>     solnVec(POIX_VAR, iimax, pos(1), pos(2)) = grav(iimax-1)
>   end if
>   if (sweepDir .eq. SWEEP_Y) then
>     solnVec(POIY_VAR, pos(1), iimin, pos(2)) = grav(iimin+1)
>     solnVec(POIY_VAR, pos(1), iimax, pos(2)) = grav(iimax-1)
>   end if
>   if (sweepDir .eq. SWEEP_Z) then
>     solnVec(POIZ_VAR, pos(1), pos(2), iimin) = grav(iimin+1)
>     solnVec(POIZ_VAR, pos(1), pos(2), iimax) = grav(iimax-1)
>   end if
>
> Slices were taken along the x-axis.
>
> Thanks for your help,
> Ryan
>
> --
> Joshua Wall
> Doctoral Candidate
> Department of Physics
> Drexel University
> 3141 Chestnut Street
> Philadelphia, PA 19104
>
> --
Joshua Wall
Doctoral Candidate
Department of Physics
Drexel University
3141 Chestnut Street
Philadelphia, PA 19104
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://flash.rochester.edu/pipermail/flash-users/attachments/20170301/716d3f0d/attachment.htm>