Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allowing the parallel for a more general 6*(2n)^2 blocks #176

Merged
merged 10 commits into from
Aug 29, 2024
Merged

Allowing the parallel for a more general 6*(2n)^2 blocks #176

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
e3c1e5b
Tentative solution to multiprocessing
Aug 25, 2024
c8dc0af
dirty fix
Aug 25, 2024
af6ffbb
fixing a bug
Aug 25, 2024
9ad385f
new strategy
Aug 25, 2024
3307a32
find block id
Aug 27, 2024
c0cebf0
Some additional fixes of find neighbor blocks
Aug 27, 2024
645f6ab
Some missing fixes
Aug 27, 2024
c0d05d1
Fixing some missing parts
Aug 27, 2024
bf156b9
Cleaned up the code, kept as a backdoor solution
Aug 28, 2024
92f243d
fixed some build errors
Aug 28, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 25 additions & 2 deletions src/exo3/cs_find_block_id.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,9 @@
#include "cubed_sphere.hpp"

int CubedSphere::FindBlockID(LogicalLocation const& loc) {
int lv2_lx2 = loc.lx2 >> (loc.level - 2);
int lv2_lx3 = loc.lx3 >> (loc.level - 2);

int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);

// Determine the block number
int block_id;
Expand All @@ -27,6 +28,7 @@ int CubedSphere::FindBlockID(LogicalLocation const& loc) {
std::stringstream msg;
msg << "Error: something wrong, check the geometry setup of the "
"cubed sphere. \n";
msg << "lv2_lx2: " << lv2_lx2 << " lv2_lx3: " << lv2_lx3 << "bound_lim: " << bound_lim << "loc.lx2: " << loc.lx2 << "loc.lx3: " << loc.lx3 << std::endl;
msg << "----------------------------------" << std::endl;
ATHENA_ERROR(msg);
}
Expand Down Expand Up @@ -73,3 +75,24 @@ int CubedSphere::FindBlockID(LogicalLocation const& loc) {

return block_id;
}

void CubedSphere::GetLocalIndex(int *lv2_lx2, int *lv2_lx3, int *local_lx2, int *local_lx3, int *bound_lim, LogicalLocation const& loc) {
#ifdef USE_NBLOCKS
// Updated method, need to manually setup in configure.hpp, allow 6*n^2 blocks
*bound_lim = (int)(sqrt(NBLOCKS / 6) - 0.5);
// Find relative location within block

*lv2_lx2 = loc.lx2 / (*bound_lim + 1);
*lv2_lx3 = loc.lx3 / (*bound_lim + 1);
*local_lx2 = loc.lx2 - (*lv2_lx2 * (*bound_lim + 1));
*local_lx3 = loc.lx3 - (*lv2_lx3 * (*bound_lim + 1));
#else
// Old method, suitable for 6*4^n blocks
*bound_lim = (1 << (loc.level - 2)) - 1;
// Find relative location within block
*lv2_lx2 = loc.lx2 >> (loc.level - 2);
*lv2_lx3 = loc.lx3 >> (loc.level - 2);
*local_lx2 = loc.lx2 - (*lv2_lx2 << (loc.level - 2));
*local_lx3 = loc.lx3 - (*lv2_lx3 << (loc.level - 2));
#endif
}
18 changes: 10 additions & 8 deletions src/exo3/cs_transform_ox.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,16 +12,12 @@ void CubedSphere::TransformOX(int *ox2, int *ox3, int *tox2, int *tox3,
LogicalLocation const &loc) {
// Find the block ID
int block_id = FindBlockID(loc);
int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);

// Find relative location within block
int lv2_lx2 = loc.lx2 >> (loc.level - 2);
int lv2_lx3 = loc.lx3 >> (loc.level - 2);
int local_lx2 = loc.lx2 - (lv2_lx2 << (loc.level - 2));
int local_lx3 = loc.lx3 - (lv2_lx3 << (loc.level - 2));
int bound_lim = (1 << (loc.level - 2)) - 1;

// Hard code the cases...
// No need to consider the corner cases, abandon in reading buffers.
// Hard code the cases...
// The corner cases are not used, abandoned after reading buffers.
int target_block = -1; // Block id of target
int target_loc_2, target_loc_3; // local x2 and x3 in target block

Expand Down Expand Up @@ -168,6 +164,7 @@ void CubedSphere::TransformOX(int *ox2, int *ox3, int *tox2, int *tox3,
std::stringstream msg;
msg << "Error: something wrong, check the geometry setup of the cubed "
"sphere. \n";
msg << "Block ID: " << block_id << std::endl;
msg << "----------------------------------" << std::endl;
ATHENA_ERROR(msg);
}
Expand Down Expand Up @@ -210,8 +207,13 @@ void CubedSphere::TransformOX(int *ox2, int *ox3, int *tox2, int *tox3,
msg << "----------------------------------" << std::endl;
ATHENA_ERROR(msg);
}
#ifdef USE_NBLOCKS
lx3_0 = (lx3_0 * (bound_lim + 1));
lx2_0 = (lx2_0 * (bound_lim + 1));
#else
lx3_0 = (lx3_0 << (loc.level - 2));
lx2_0 = (lx2_0 << (loc.level - 2));
#endif
// Add up first block and local positions
int lx3_t = lx3_0 + target_loc_3;
int lx2_t = lx2_0 + target_loc_2;
Expand Down
85 changes: 41 additions & 44 deletions src/exo3/cubed_sphere.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,9 +39,10 @@ CubedSphere::CubedSphere(MeshBlock *pmb) : pmy_block_(pmb) {

Real CubedSphere::GenerateMeshX2(Real x, LogicalLocation const &loc) {
Real x_l, x_u;
int lx2_lv2 = loc.lx2 >> (loc.level - 2);
int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);

switch (lx2_lv2) {
switch (lv2_lx2) {
case 0:
x_l = -0.5;
x_u = 0.0;
Expand All @@ -56,9 +57,10 @@ Real CubedSphere::GenerateMeshX2(Real x, LogicalLocation const &loc) {

Real CubedSphere::GenerateMeshX3(Real x, LogicalLocation const &loc) {
Real x_l, x_u;
int lx3_lv2 = loc.lx3 >> (loc.level - 2);
int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);

switch (lx3_lv2) {
switch (lv2_lx3) {
case 0:
x_l = -0.5;
x_u = -1.0 / 6.0;
Expand All @@ -81,10 +83,8 @@ Real CubedSphere::GenerateMeshX3(Real x, LogicalLocation const &loc) {
void CubedSphere::GetLatLon(Real *lat, Real *lon, int k, int j, int i) const {
auto pcoord = pmy_block_->pcoord;
auto &loc = pmy_block_->loc;
int blockID = FindBlockID(loc);

int lv2_lx2 = loc.lx2 >> (loc.level - 2);
int lv2_lx3 = loc.lx3 >> (loc.level - 2);
int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
// Calculate the needed parameters
Real dX = tan(pcoord->x2v(j));
Real dY = tan(pcoord->x3v(k));
Expand All @@ -99,10 +99,7 @@ void CubedSphere::GetLatLonFace2(Real *lat, Real *lon, int k, int j,
int i) const {
auto pcoord = pmy_block_->pcoord;
auto &loc = pmy_block_->loc;

int lv2_lx2 = loc.lx2 >> (loc.level - 2);
int lv2_lx3 = loc.lx3 >> (loc.level - 2);
int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
int blockID = FindBlockID(loc);

// Calculate the needed parameters
Real dX = tan(pcoord->x2f(j));
Expand All @@ -118,10 +115,7 @@ void CubedSphere::GetLatLonFace3(Real *lat, Real *lon, int k, int j,
int i) const {
auto pcoord = pmy_block_->pcoord;
auto &loc = pcoord->pmy_block->loc;

int lv2_lx2 = loc.lx2 >> (loc.level - 2);
int lv2_lx3 = loc.lx3 >> (loc.level - 2);
int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
int blockID = FindBlockID(loc);

// Calculate the needed parameters
Real dX = tan(pcoord->x2v(j));
Expand All @@ -136,10 +130,7 @@ void CubedSphere::GetUV(Real *U, Real *V, Real V2, Real V3, int k, int j,
int i) const {
auto pcoord = pmy_block_->pcoord;
auto &loc = pmy_block_->loc;

int lv2_lx2 = loc.lx2 >> (loc.level - 2);
int lv2_lx3 = loc.lx3 >> (loc.level - 2);
int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
int blockID = FindBlockID(loc);

Real X = tan(pcoord->x2v(j));
Real Y = tan(pcoord->x3v(k));
Expand All @@ -153,10 +144,7 @@ void CubedSphere::GetVyVz(Real *V2, Real *V3, Real U, Real V, int k, int j,
int i) const {
auto pcoord = pmy_block_->pcoord;
auto &loc = pmy_block_->loc;

int lv2_lx2 = loc.lx2 >> (loc.level - 2);
int lv2_lx3 = loc.lx3 >> (loc.level - 2);
int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
int blockID = FindBlockID(loc);

// Calculate the needed parameters
Real X = tan(pcoord->x2v(j));
Expand All @@ -169,10 +157,7 @@ void CubedSphere::CalculateCoriolisForce2(int i2, int i3, Real v2, Real v3,
Real *cF3) const {
auto pcoord = pmy_block_->pcoord;
auto &loc = pmy_block_->loc;

int lv2_lx2 = loc.lx2 >> (loc.level - 2);
int lv2_lx3 = loc.lx3 >> (loc.level - 2);
int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
int blockID = FindBlockID(loc);

Real x = tan(pcoord->x2v(i2));
Real y = tan(pcoord->x3v(i3));
Expand Down Expand Up @@ -203,10 +188,7 @@ void CubedSphere::CalculateCoriolisForce3(int i2, int i3, Real v1, Real v2,
Real *cF3) const {
auto pcoord = pmy_block_->pcoord;
auto &loc = pmy_block_->loc;

int lv2_lx2 = loc.lx2 >> (loc.level - 2);
int lv2_lx3 = loc.lx3 >> (loc.level - 2);
int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
int blockID = FindBlockID(loc);

Real x = tan(pcoord->x2v(i2));
Real y = tan(pcoord->x3v(i3));
Expand Down Expand Up @@ -333,14 +315,10 @@ void CubedSphere::sendNeighborBlocks(int ox2, int ox3, int tg_rank,
int tg_gid) {
MeshBlock *pmb = pmy_block_;
auto &loc = pmb->loc;

int lv2_lx2 = loc.lx2 >> (loc.level - 2);
int lv2_lx3 = loc.lx3 >> (loc.level - 2);
int blockID = lv2_lx2 + lv2_lx3 * 2 + 1;
int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);
int blockID = FindBlockID(loc);
// Calculate local ID
int local_lx2 = loc.lx2 - (lv2_lx2 << (loc.level - 2));
int local_lx3 = loc.lx3 - (lv2_lx3 << (loc.level - 2));
int bound_lim = (1 << (loc.level - 2)) - 1;
int tox2, tox3;
int DirTag, DirNum, ownTag; // Tag for target, and the numbering of axis
int ox2_bkp = ox2;
Expand Down Expand Up @@ -560,6 +538,16 @@ void CubedSphere::sendNeighborBlocks(int ox2, int ox3, int tg_rank,
}

// Calculate the tag of destination
#ifdef USE_NBLOCKS
if (tox2 == -1)
DirTag = 0 + 4 * pmb->gid + 24 * (bound_lim + 1) * tg_gid;
if (tox2 == 1)
DirTag = 1 + 4 * pmb->gid + 24 * (bound_lim + 1) * tg_gid;
if (tox3 == -1)
DirTag = 2 + 4 * pmb->gid + 24 * (bound_lim + 1) * tg_gid;
if (tox3 == 1)
DirTag = 3 + 4 * pmb->gid + 24 * (bound_lim + 1) * tg_gid;
#else
if (tox2 == -1)
DirTag = 0 + 4 * pmb->gid + 24 * (1 << (loc.level - 2)) * tg_gid;
if (tox2 == 1)
Expand All @@ -568,6 +556,7 @@ void CubedSphere::sendNeighborBlocks(int ox2, int ox3, int tg_rank,
DirTag = 2 + 4 * pmb->gid + 24 * (1 << (loc.level - 2)) * tg_gid;
if (tox3 == 1)
DirTag = 3 + 4 * pmb->gid + 24 * (1 << (loc.level - 2)) * tg_gid;
#endif
// Send by MPI: we don't care whether it is in the same process for now
if (ox2 == -1) ownTag = 0;
if (ox2 == 1) ownTag = 1;
Expand All @@ -586,13 +575,10 @@ void CubedSphere::recvNeighborBlocks(int ox2, int ox3, int tg_rank,
MeshBlock *pmb = pmy_block_;
auto &loc = pmb->loc;

int lv2_lx2 = loc.lx2 >> (loc.level - 2);
int lv2_lx3 = loc.lx3 >> (loc.level - 2);
int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);
int blockID = FindBlockID(loc);
// Calculate local ID
int local_lx2 = loc.lx2 - (lv2_lx2 << (loc.level - 2));
int local_lx3 = loc.lx3 - (lv2_lx3 << (loc.level - 2));
int bound_lim = (1 << (loc.level - 2)) - 1;

int tox2, tox3;
int DirTag, DirNum, ownTag; // Tag for receiving, and the numbering of axis
// to place the values
Expand Down Expand Up @@ -658,6 +644,16 @@ void CubedSphere::recvNeighborBlocks(int ox2, int ox3, int tg_rank,
int dsize = ((kb2 - kb1 + 1) * (jb2 - jb1 + 1) * (ib2 - ib1 + 1) * NWAVE);
Real *data = new Real[dsize];
// Calculate the tag for receiving
#ifdef USE_NBLOCKS
if (ox2 == -1)
DirTag = 0 + 4 * tg_gid + 24 * (bound_lim + 1) * pmb->gid;
if (ox2 == 1)
DirTag = 1 + 4 * tg_gid + 24 * (bound_lim + 1) * pmb->gid;
if (ox3 == -1)
DirTag = 2 + 4 * tg_gid + 24 * (bound_lim + 1) * pmb->gid;
if (ox3 == 1)
DirTag = 3 + 4 * tg_gid + 24 * (bound_lim + 1) * pmb->gid;
#else
if (ox2 == -1)
DirTag = 0 + 4 * tg_gid + 24 * (1 << (loc.level - 2)) * pmb->gid;
if (ox2 == 1)
Expand All @@ -666,6 +662,7 @@ void CubedSphere::recvNeighborBlocks(int ox2, int ox3, int tg_rank,
DirTag = 2 + 4 * tg_gid + 24 * (1 << (loc.level - 2)) * pmb->gid;
if (ox3 == 1)
DirTag = 3 + 4 * tg_gid + 24 * (1 << (loc.level - 2)) * pmb->gid;
#endif

#ifdef MPI_PARALLEL
MPI_Recv(data, dsize, MPI_DOUBLE, tg_rank, DirTag, MPI_COMM_WORLD,
Expand Down
4 changes: 3 additions & 1 deletion src/exo3/cubed_sphere.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,9 @@ class CubedSphere {
~CubedSphere() {}

static int FindBlockID(LogicalLocation const &loc);
static void GetLocalIndex(int *lv2_lx2, int *lv2_lx3, int *local_lx2,
int *local_lx3, int *bound_lim,
LogicalLocation const &loc);
static void TransformOX(int *ox2, int *ox3, int *tox2, int *tox3,
LogicalLocation const &loc);

Expand Down Expand Up @@ -66,7 +69,6 @@ class CubedSphere {
#endif

std::vector<Real> LRDataBuffer[4];

MeshBlock *pmy_block_;
};

Expand Down
21 changes: 8 additions & 13 deletions src/exo3/cubed_sphere_utility.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,10 +102,9 @@ void InteprolateX2(const AthenaArray<Real> &src, AthenaArray<Real> &tgt,
int sk, int ek, int TypeFlag) {
// Interpolation along X2 (j) axis, used before sending data to X3 (k) axis
// Get the local indices
int lv2_lx2 = loc.lx2 >> (loc.level - 2);
int lv2_lx3 = loc.lx3 >> (loc.level - 2);
int local_lx2 = loc.lx2 - (lv2_lx2 << (loc.level - 2));
int bound_lim = (1 << (loc.level - 2)) - 1;
int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
CubedSphere::GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);

int meshblock_size = ej - sj + 1;
int N_blk = meshblock_size *
(bound_lim + 1); // N in X2 direction for each panel. This value
Expand Down Expand Up @@ -233,10 +232,9 @@ void InteprolateX3(const AthenaArray<Real> &src, AthenaArray<Real> &tgt,
int sk, int ek, int TypeFlag) {
// Interpolation along X3 (k) axis, used before sending data to ghost zone in
// X2 (j) direction Get the local indices
int lv2_lx2 = loc.lx2 >> (loc.level - 2);
int lv2_lx3 = loc.lx3 >> (loc.level - 2);
int local_lx3 = loc.lx3 - (lv2_lx3 << (loc.level - 2));
int bound_lim = (1 << (loc.level - 2)) - 1;
int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
CubedSphere::GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);

int meshblock_size = ek - sk + 1;
int N_blk = meshblock_size *
(bound_lim + 1); // N in X2 direction for each panel. This value
Expand Down Expand Up @@ -403,11 +401,8 @@ void PackData(const AthenaArray<Real> &src, Real *buf, int sn, int en, int si,
}

// Get the local indices
int lv2_lx2 = loc.lx2 >> (loc.level - 2);
int lv2_lx3 = loc.lx3 >> (loc.level - 2);
int local_lx2 = loc.lx2 - (lv2_lx2 << (loc.level - 2));
int local_lx3 = loc.lx3 - (lv2_lx3 << (loc.level - 2));
int bound_lim = (1 << (loc.level - 2)) - 1;
int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim;
CubedSphere::GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc);

// Work on interpolation
AthenaArray<Real> interpolatedSrc;
Expand Down
Loading