[aspect-devel] Advice for optimizing particle advection for 2D Spherical Shell Models

[Harsha Lokavarapu]

Hi,

As part of a subgroup in Davis, we have been actively studying active
particles when we began to notice that under certain specific situations
the scaling of our models degraded. This occurs only for 2D spherical
annulus domains. The culprit was the particle advection step. Further
investigation showed that the bottleneck within the advection step lay in
the update to the reference location. As we increase the number of
particles per cell for a fixed grid refinement and fixed MPI processes, we
noticed a linear growth in the wall clock time of
transform_from_real_to_unit_cell(). Sometimes, we were looking at 60-70% of
the wall clock time spent in this portion of the code. After several
discussions, we settled on a possible solution. However, attempts to
implement the solution have proven more difficult that expected.

To reproduce the bottleneck, I have attached a modified annulus.prm
benchmark as the test case. This will need to be linked to the compiled
additional library of aspect/benchmarks/annulus. I have added timers to the
local particle advection function to measure the different bits and pieces.
Feel free to compile the version of the code branch at
https://github.com/hlokavarapu/aspect/tree/time_particle_advection

Assuming that the working directory is aspect/benchmarks/annulus ,  for
bash users, run

$ for i in 2 4 8; do sed -i "s/set Number of particles per cell per
direction.*/set Number of particles per cell per direction = $i/"
annulus.prm; mpirun -np 1 aspect annulus.prm >> advect.log ;done

The results will be in the file advect.log.

Anyways, as a starting point, we decided to pick the simplest case, namely
the Euler integration scheme to test our idea. Currently the physical
location of the particle is updated using Euler. Similarly, we would like
to introduce a second Euler step to update the reference location of the
particle. However, this requires the velocity solution at the support
points in the reference cell. If we are using a Q2 velocity element, we can
get away with 9 transformations.

After, burrowing through the Deal.II documentation, I have found
implemented functionality for transforming vectors from unit to real but
not the other way around. Is there a black box solution that already exists
within Deal.II which converts vectors from real to unit?

We also looked through the implementation of transform_from_real_to_unit_cell()
which uses a newton iteration to transform a point in physical space to a
point in the unit space. Will a similar technique be necessary to transform
a vector from real to unit?

Also, are there any locations within Aspect that involves working with the
velocity in the reference cell?

For a box domain, the transform_from_real_to_unit is orders faster than in
the spherical shell domain. Is this because we have a better initial guess
to the newton iteration?

We (I) have been stumped on this for close to a week now. Any leads are
highly appreciated.

Thank You,
Harsha
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20170905/d5374904/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: annulus.prm
Type: application/octet-stream
Size: 3933 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20170905/d5374904/attachment.obj>