Integrate on local subsets of domain on distributed triangulation

[Lucas Myers]

Hi folks,

I'd like to coarse-grain a finite element field, where the coarse-graining procedure is just convolving the field with a Gaussian function. More concretely, I'd like to do the following:

For a regular Triangulation, we can calculate <X> at every quadrature point by iterating through the entire Triangulation and calculating the integrand at all of the quadrature points. In fact, because there is a Gaussian in distance in the integrand, we only have to look at cells a distance ~3a_0 away.

However, for a distributed triangulation it could be the case that cells 3a_0 away may be on a different MPI process. My first question is: is there a way to extend the number of ghost cells along the boundary such that they are some total distance away from the locally-owned cells on the edge of the domain? And if not, is there a nice utility to allow me to access neighboring cells on other domains (albeit, at the cost of some MPI messages)?

Otherwise, is there a good way to do this? My guess would be to do all of the local integrations first, and make a vector of cells which end up butting against the boundary. Then one could use the `ghost_owners` function to send information about the edge cells to the adjacent subdomains, which would then send information back. However that seems a little in the weeds of using MPI and also would require some error-prone bookkeeping, which I would like to avoid if possible.

Thanks so much for any help!

- Lucas

[Lucas Myers]

Ack, I realize the pasted screenshot of the relevant equation has disappeared! I've included it as an attachment to this reply.

- Lucas

[Wolfgang Bangerth]

> I'd like to coarse-grain a finite element field, where the coarse-graining
> procedure is just convolving the field with a Gaussian function. More
> concretely, I'd like to do the following:
>
>

> For a regular Triangulation, we can calculate <X> at every quadrature point by
> iterating through the entire Triangulation and calculating the integrand at
> all of the quadrature points. In fact, because there is a Gaussian in distance
> in the integrand, we only have to look at cells a distance ~3a_0 away.
>
> However, for a distributed triangulation it could be the case that cells 3a_0
> away may be on a different MPI process. My first question is: is there a way
> to extend the number of ghost cells along the boundary such that they are some
> total distance away from the locally-owned cells on the edge of the domain?

p4est (our parallel mesh engine) actually allows for "thicker" ghost layers,
but we never tried that in deal.II and I doubt that it would work without much
trouble. Either way, p4est allows having a ghost layer of N cells width
whereas you need one that has a specific Euclidean distance from the locally
owned cells -- so these are different concepts anyway.

> And if not, is there a nice utility to allow me to access neighboring cells on
> other domains (albeit, at the cost of some MPI messages)?

Check out the FEPointEvaluation and Utilities::MPI::RemotePointEvaluation
classes. I think they will do what you need.

> Otherwise, is there a good way to do this? My guess would be to do all of the
> local integrations first, and make a vector of cells which end up butting
> against the boundary. Then one could use the `ghost_owners` function to send
> information about the edge cells to the adjacent subdomains, which would then
> send information back. However that seems a little in the weeds of using MPI
> and also would require some error-prone bookkeeping, which I would like to
> avoid if possible.

Yes, and there's a bug in your algorithm as well: Just because the ghost cells
are owned by a specific set of processed doesn't mean that the cell you need
to evaluate the solution on is actually owned by one of the neighboring processes.

Secondly, what you *really* want to do is assess first which points you need
information on from other processes, send those points around, and while you
are waiting for these messages to be delivered and answered, you want to work
on your local cells -- not the other way around: You want to hide
communication latency behind local work.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Lucas Myers]

Hi Wolfgang,

Thanks for the help! Point taken about neighboring processes and communication latency.

Re: FEPointEvaluation and RemotePointEvaluation, I think I actually need more than that because I need the quadrature points on each of the cells that I'm querying, so I can just do the standard procedure with FEValues and the corresponding Quadrature object.

I'll think harder about how to avoid latency issues.

- Lucas

[Wolfgang Bangerth]

On 4/25/23 09:44, Lucas Myers wrote:
>
> Re: FEPointEvaluation and RemotePointEvaluation, I think I actually need
> more than that because I need the quadrature points on each of the cells
> that I'm querying, so I can just do the standard procedure with FEValues
> and the corresponding Quadrature object.

In that case, the way I would implement this is that every process makes
a list of circles: points at which it wants to evaluate the convolution,
and the radius of the convolution.

In a second step, you then determine which of these circles require
communication. This isn't actually entirely trivial, but a good
criterion is probably whether a circle has overlap with an artificial
cell. So:
std::set<Circle> all_circles = ...;
std::set<Circle> circles_requiring_communication;
for (cell=...)
if (cell->is_artificial())
for (c : all_circles)
if (c not already in circles_requiring_communication)
if (c overlaps with cell)
circles_requiring_communication.insert (c);

In a third step, you can either send all of these circles to all other
processes (via Utilities::MPI::all_gather()), or you can be a bit
smarter and first get a bounding box for each of the other processes'
locally owned cells (see the function in GridTools for that); in that
latter case, you'd only send a circle to a process if it overlaps with
that other process's bounding box, using one of the ConsensusAlgorithms
functions.

At this point, each process has a list of circles that overlap with its
own subdomain. Now you work on that by looping over all of its locally
owned cells, on each cell query all circles, and do what you need to do.
This results in a list of contributions that need to be sent back to the
original process requesting that information. This could probably be
done nicely using the request-answer system of the ConsensusAlgorithms
functions.

You'd need to be careful with cases where a process gets sent a circle
that has overlap with its bounding box, but not with any of its cells.
In that case, the reply is simply a zero contribution to the convolution
integral.

I hope this makes sense!

[Lucas Myers]

Hi Wolfgang,

This is extremely helpful, thanks so much!

- Lucas