AMReX-Codes / amrex

AMReX: Software Framework for Block Structured AMR
https://amrex-codes.github.io/amrex
Other
546 stars 347 forks source link

Request for supporting non-local boundary conditions #1538

Open eschnett opened 3 years ago

eschnett commented 3 years ago

The Einstein Toolkit supports "non-local" boundary conditions. These are similar in spirit to periodic boundaries, but have a different definition. For example, there is a 90-degree rotational symmetry where one simulations one quadrant of the domain, where boundaries are defined by rotating that quadrant by 90 degrees in either direction. Implementing such boundary conditions requires communication, and I thus hope it can be supported by AMReX in some way.

There are three kinds of such boundary conditions we would like to have supported:

I implemented the respective boundary conditions and described them in more detail in serial Fortran 90 in this gist: https://gist.github.com/eschnett/c7852a5c5982ed501cbaa89d9b600b48.

I am looking for functions in AMReX that would handle the communication necessary to support these boundary conditions.

zingale commented 3 years ago

this sounds related to https://github.com/AMReX-Astro/Castro/issues/602

eschnett commented 3 years ago

Yes, this is probably what I called "90 degree rotational" symmetry.

Implementing this in the same way as periodic boundary conditions would be convenient. However (as you say in your issue), one also needs to swap certain tensor components. The implementation I had in mind was a generic API to copy rectangular chunks of data between multifabs, while possibly rotating/mirroring the slabs. Handling signs and exchanging tensor components would then be handled by the caller.

WeiqunZhang commented 3 years ago

In the Fortran example, there are src with the intent of in and dst with the intent of inout. But if you do a copy from src to dst before calling it, it seems that the function only needs to take one array. Is that correct?

eschnett commented 3 years ago

I think one array suffices for the "rotating 180degree" and "polar" boundary conditions. In our code, using one array would not work for the "rotating 90 degree" boundary condition: Let's assume we have three variables vx, vy, and vz for the velocity. Then the boundary of vx is filled by values from vy and vice versa, whereas vz is filled from itself.

Actually, combining src and dst would work, if we later add local loops that exchange data between vx and vy. That would not be too expensive since we need to fix up the signs anyway, unless we add flags to the API describing which signs should be changed while copying values.

WeiqunZhang commented 3 years ago

OK. Then I will work on the one array versions. You can then make local changes.

WeiqunZhang commented 3 years ago

@eschnett I have implemented (1) 90 degree rotation and (2) 180 degree rotation (amrex/Src/Base/AMReX_NonLocalBC.H). Please give them a test. For (3) "across the pole" on a sphere, I want to make sure I understand it. Is the following correct?

// Fill the lo-y and hi-y boundary regions (including the four corners with
// f(x,y) = f(x+Lx/2,y+Ly) and f(x,y) = f(x-Lx/2,y+Ly).  Note that this also
// implies that f(x,y) = f(x+Lx,y) and f(x,y) = f(x,y+2*Ly).
eschnett commented 3 years ago

The "across the pole" boundary conditions come from discretizing the surface of a sphere with coordinates theta and phi (0 <= theta <= pi and 0 <= phi < 2pi). The system is periodic in the phi direction (f(theta, phi+2pi) = f(theta,phi)), as you stated. The non-local boundary conditions are theta=0 and theta=pi boundaries. The conditions there are f(-theta,phi) = f(theta,phi+pi) and f(pi+theta,phi) = f(pi-theta,phi+pi), respectively. In both cases, the expression phi+pi may have to be normalized to be in the range [0, 2pi).

In other words, both the theta=0 and theta=pi boundaries are essentially reflecting boundaries, but with the added complexity that they are also shifted by 1/2 the domain size in the phi direction.

f(1:nghosts, 1:nphi/2, k) = f(2*nghosts:nghosts+1:-1, nphi/2+1:nphi, k)
f(1:nghosts, nphi/2+1:nphi, k) = f(2*nghosts:nghosts+1:-1, 1:nphi/2, k)

f(ntheta-nghosts+1:ntheta, 1:nphi/2, k) = f(ntheta-nghosts:ntheta-2*nghosts:-1,  nphi/2+1:nphi, k)
f(ntheta-nghosts+1:ntheta, nphi/2+1:nphi, k) = f(ntheta-nghosts:ntheta-2*nghosts:-1,  1:nphi/2, k)
WeiqunZhang commented 3 years ago

OK. Got it.

WeiqunZhang commented 3 years ago

@eschnett I have also implemented the polar boundary. Please give them a try. Note that currently they can only run on CPUs. I still need to port them to GPUs. But that should be straightforward once we have verifed the CPU versions work.

eschnett commented 3 years ago

@WeiqunZhang Thanks a lot! That's fine, we'll want to test on CPUs first anyway.

WeiqunZhang commented 3 years ago

GPU support has been added.