[FLASH-USERS] DustCollapse Weird Artifacts in Gravitational Acceleration

Tue Feb 28 23:57:54 EST 2017

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/Gravity
> Main/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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://flash.rochester.edu/pipermail/flash-users/attachments/20170228/f8822a9b/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: dustcoll_bhtreeAcc1_hdf5_chk_0000_Slice_x_gacy.png
Type: image/png
Size: 35572 bytes
Desc: not available
URL: <http://flash.rochester.edu/pipermail/flash-users/attachments/20170228/f8822a9b/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: dustcoll_hdf5_chk_0000_Slice_x_gacy.png
Type: image/png
Size: 36012 bytes
Desc: not available
URL: <http://flash.rochester.edu/pipermail/flash-users/attachments/20170228/f8822a9b/attachment-0001.png>