[FLASH-USERS] DustCollapse Weird Artifacts in Gravitational Acceleration
Ryan Farber
rjfarber at umich.edu
Wed Mar 1 16:38:32 EST 2017
Hi Josh,
I put the "Grid_fillGuardCells" call right before my calls to
Gravity_acclListOfBlocsk in Driver_initFlash (and after
Gravity_potentialListOfBlocks) but that had no effect (still the jumps on
the boundary, plot looks identical to the one sent in the last email).
Best,
--------
Ryan
On Wed, Mar 1, 2017 at 2:16 PM, Joshua Wall <joshua.e.wall at gmail.com> wrote:
> 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/66b73abc/attachment.htm>
More information about the flash-users
mailing list