diff --git a/src/fluid/boundary/axis.cpp b/src/fluid/boundary/axis.cpp index 47f0c2954..45ea49ea5 100644 --- a/src/fluid/boundary/axis.cpp +++ b/src/fluid/boundary/axis.cpp @@ -393,12 +393,15 @@ void Axis::ExchangeMPI(int side) { #ifdef WITH_MPI // Load the buffers with data int ibeg,iend,jbeg,jend,kbeg,kend,offset; - int nx,ny,nz; - auto bufferSend = this->bufferSend; + int ny; + Buffer bufferSend = this->bufferSend; IdefixArray1D map = this->mapVars; IdefixArray4D Vc = this->Vc; IdefixArray4D Vs = this->Vs; + //reset pointer + bufferSend.ResetPointer(); + // If MPI Persistent, start receiving even before the buffers are filled MPI_Status sendStatus; @@ -409,64 +412,58 @@ void Axis::ExchangeMPI(int side) { idfx::mpiCallsTimer += MPI_Wtime() - tStart; // Coordinates of the ghost region which needs to be transfered - ibeg = 0; - iend = data->np_tot[IDIR]; - nx = data->np_tot[IDIR]; // Number of points in x - jbeg = 0; - jend = data->nghost[JDIR]; offset = data->end[JDIR]; // Distance between beginning of left and right ghosts ny = data->nghost[JDIR]; - kbeg = data->beg[KDIR]; - kend = data->end[KDIR]; - nz = kend - kbeg; + + // Create the base box to be patched by the ops + BoundingBox baseBox; + baseBox[IDIR][0] = 0; + baseBox[IDIR][1] = data->np_tot[IDIR]; + baseBox[JDIR][0] = 0; + baseBox[JDIR][1] = data->nghost[JDIR]; + baseBox[KDIR][0] = data->beg[KDIR]; + baseBox[KDIR][1] = data->end[KDIR]; + if(side==left) { - idefix_for("LoadBufferX2Vc",0,mapNVars,kbeg,kend,jbeg,jend,ibeg,iend, - KOKKOS_LAMBDA (int n, int k, int j, int i) { - bufferSend(i + (j-jbeg)*nx + (k-kbeg)*nx*ny + n*nx*ny*nz ) = - Vc(map(n),k,j+ny,i); - } - ); + //shift by NY + BoundingBox sendBoxVc = baseBox; + sendBoxVc[JDIR][0] += ny; + sendBoxVc[JDIR][1] += ny; + bufferSend.Pack(Vc, map, sendBoxVc); + if (haveMHD) { - int VsIndex = mapNVars*nx*ny*nz; - idefix_for("LoadBufferX2IDIR",kbeg,kend,jbeg,jend,ibeg,iend+1, - KOKKOS_LAMBDA (int k, int j, int i) { - bufferSend(i + (j-jbeg)*(nx+1) + (k-kbeg)*(nx+1)*ny + VsIndex ) = - Vs(IDIR,k,j+ny,i); - } - ); - VsIndex = mapNVars*nx*ny*nz + (nx+1)*ny*nz; - idefix_for("LoadBufferX2KDIR",kbeg,kend+1,jbeg,jend,ibeg,iend, - KOKKOS_LAMBDA (int k, int j, int i) { - bufferSend(i + (j-jbeg)*nx + (k-kbeg)*nx*ny + VsIndex ) = - Vs(KDIR,k,j+ny,i); - } - ); + //extend by one the end on idir + BoundingBox sendBoxVsIdir = sendBoxVc; + sendBoxVsIdir[IDIR][1] += 1; + bufferSend.Pack(Vs, IDIR, sendBoxVsIdir); + + //extend by one the end on kdir + BoundingBox sendBoxVsKdir = sendBoxVc; + sendBoxVsKdir[KDIR][1] += 1; + bufferSend.Pack(Vs, KDIR, sendBoxVsKdir); } // MHD } else if(side==right) { - idefix_for("LoadBufferX2Vc",0,mapNVars,kbeg,kend,jbeg,jend,ibeg,iend, - KOKKOS_LAMBDA (int n, int k, int j, int i) { - bufferSend(i + (j-jbeg)*nx + (k-kbeg)*nx*ny + n*nx*ny*nz ) = - Vc(map(n),k,j+offset-ny,i); - } - ); + //shift by offset - NY (to select the right part) + BoundingBox sendBoxVc = baseBox; + sendBoxVc[JDIR][0] += offset - ny; + sendBoxVc[JDIR][1] += offset - ny; + bufferSend.Pack(Vc, map, sendBoxVc); // Load face-centered field in the buffer if (haveMHD) { - int VsIndex = mapNVars*nx*ny*nz; - idefix_for("LoadBufferX2IDIR",kbeg,kend,jbeg,jend,ibeg,iend+1, - KOKKOS_LAMBDA (int k, int j, int i) { - bufferSend(i + (j-jbeg)*(nx+1) + (k-kbeg)*(nx+1)*ny + VsIndex ) = - Vs(IDIR,k,j+offset-ny,i); - } - ); - VsIndex = mapNVars*nx*ny*nz + (nx+1)*ny*nz; - - idefix_for("LoadBufferX2KDIR",kbeg,kend+1,jbeg,jend,ibeg,iend, - KOKKOS_LAMBDA (int k, int j, int i) { - bufferSend(i + (j-jbeg)*nx + (k-kbeg)*nx*ny + VsIndex ) = - Vs(KDIR,k,j+offset-ny,i); - } - ); + //extend by one the end on idir + take the right part or the mesh on jdir + BoundingBox sendBoxVsIdir = baseBox; + sendBoxVsIdir[IDIR][1] += 1; + sendBoxVsIdir[JDIR][0] += offset - ny; + sendBoxVsIdir[JDIR][1] += offset - ny; + bufferSend.Pack(Vs, IDIR, sendBoxVsIdir); + + //extend by one the end on kdir + take the right part or the mesh on jdir + BoundingBox sendBoxVsKdir = baseBox; + sendBoxVsKdir[KDIR][1] += 1; + sendBoxVsKdir[JDIR][0] += offset - ny; + sendBoxVsKdir[JDIR][1] += offset - ny; + bufferSend.Pack(Vs, KDIR, sendBoxVsKdir); } // MHD } // side==right @@ -478,63 +475,53 @@ void Axis::ExchangeMPI(int side) { idfx::mpiCallsTimer += MPI_Wtime() - tStart; // Unpack - auto bufferRecv=this->bufferRecv; + Buffer bufferRecv=this->bufferRecv; + bufferRecv.ResetPointer(); auto sVc = this->symmetryVc; if(side==left) { - idefix_for("StoreBufferAxis",0,mapNVars,kbeg,kend,jbeg,jend,ibeg,iend, - KOKKOS_LAMBDA (int n, int k, int j, int i) { - Vc(map(n),k,jend-(j-jbeg)-1,i) = - sVc(map(n))*bufferRecv(i + (j-jbeg)*nx + (k-kbeg)*nx*ny + n*nx*ny*nz ); - } - ); + //unpack Vc + BoundingBox recvBoxVc = baseBox; + bufferRecv.UnpackJDirSymetric(Vc, map, sVc, recvBoxVc); // Load face-centered field in the buffer if (haveMHD) { - int VsIndex = mapNVars*nx*ny*nz; auto sVs = this->symmetryVs; - idefix_for("StoreBufferX2IDIR",kbeg,kend,jbeg,jend,ibeg,iend+1, - KOKKOS_LAMBDA (int k, int j, int i) { - Vs(IDIR,k,jend-(j-jbeg)-1,i) = - sVs(IDIR)*bufferRecv(i + (j-jbeg)*(nx+1) + (k-kbeg)*(nx+1)*ny + VsIndex ); - } - ); - VsIndex = mapNVars*nx*ny*nz + (nx+1)*ny*nz; + //unpack Vs face-centered + BoundingBox recvBoxVsIdir = baseBox; + recvBoxVsIdir[IDIR][1] += 1; + bufferRecv.UnpackJDirSymetric(Vs, IDIR, sVs(IDIR), recvBoxVsIdir); - idefix_for("StoreBufferX2KDIR",kbeg,kend+1,jbeg,jend,ibeg,iend, - KOKKOS_LAMBDA ( int k, int j, int i) { - Vs(KDIR,k,jend-(j-jbeg)-1,i) = - sVs(KDIR)*bufferRecv(i + (j-jbeg)*nx + (k-kbeg)*nx*ny + VsIndex ); - } - ); + //unpack Vs face-centered + BoundingBox recvBoxVsKdir = baseBox; + recvBoxVsKdir[KDIR][1] += 1; + bufferRecv.UnpackJDirSymetric(Vs, KDIR, sVs(KDIR), recvBoxVsKdir); } } else if(side==right) { - idefix_for("StoreBufferAxis",0,mapNVars,kbeg,kend,jbeg,jend,ibeg,iend, - KOKKOS_LAMBDA (int n, int k, int j, int i) { - Vc(map(n),k,jend-(j-jbeg)-1+offset,i) = - sVc(map(n))*bufferRecv(i + (j-jbeg)*nx + (k-kbeg)*nx*ny + n*nx*ny*nz ); - } - ); + //unpack Vc on right part + BoundingBox recvBoxVc = baseBox; + recvBoxVc[JDIR][0] += offset; + recvBoxVc[JDIR][1] += offset; + bufferRecv.UnpackJDirSymetric(Vc, map, sVc, recvBoxVc); // Load face-centered field in the buffer if (haveMHD) { - int VsIndex = mapNVars*nx*ny*nz; auto sVs = this->symmetryVs; - idefix_for("StoreBufferX2IDIR",kbeg,kend,jbeg,jend,ibeg,iend+1, - KOKKOS_LAMBDA (int k, int j, int i) { - Vs(IDIR,k,jend-(j-jbeg)-1+offset,i) = - sVs(IDIR)*bufferRecv(i + (j-jbeg)*(nx+1) + (k-kbeg)*(nx+1)*ny + VsIndex ); - } - ); - VsIndex = mapNVars*nx*ny*nz + (nx+1)*ny*nz; - idefix_for("StoreBufferX2KDIR",kbeg,kend+1,jbeg,jend,ibeg,iend, - KOKKOS_LAMBDA ( int k, int j, int i) { - Vs(KDIR,k,jend-(j-jbeg)-1+offset,i) = - sVs(KDIR)*bufferRecv(i + (j-jbeg)*nx + (k-kbeg)*nx*ny + VsIndex ); - } - ); + //unpack Vs face-centered on right part + BoundingBox recvBoxVsIdir = baseBox; + recvBoxVsIdir[IDIR][1] += 1; + recvBoxVsIdir[JDIR][0] += offset; + recvBoxVsIdir[JDIR][1] += offset; + bufferRecv.UnpackJDirSymetric(Vs, IDIR, sVs(IDIR), recvBoxVsIdir); + + //unpack Vs face-centered on right part + BoundingBox recvBoxVsKdir = baseBox; + recvBoxVsKdir[KDIR][1] += 1; + recvBoxVsKdir[JDIR][0] += offset; + recvBoxVsKdir[JDIR][1] += offset; + bufferRecv.UnpackJDirSymetric(Vs, KDIR, sVs(KDIR), recvBoxVsKdir); } // MHD } @@ -594,8 +581,9 @@ void Axis::InitMPI() { #endif // DIMENSIONS } - this->bufferRecv = IdefixArray1D("bufferRecvAxis", bufferSize); - this->bufferSend = IdefixArray1D("bufferSendAxis", bufferSize); + //build buffers + this->bufferRecv = Buffer(bufferSize); + this->bufferSend = Buffer(bufferSize); // init persistent communications // We receive from procRecv, and we send to procSend @@ -605,10 +593,10 @@ void Axis::InitMPI() { MPI_SAFE_CALL(MPI_Cart_shift(data->mygrid->AxisComm,0,data->mygrid->nproc[KDIR]/2, &procRecv,&procSend )); - MPI_SAFE_CALL(MPI_Send_init(bufferSend.data(), bufferSize, realMPI, procSend, + MPI_SAFE_CALL(MPI_Send_init(bufferSend.data(), bufferSend.Size(), realMPI, procSend, 650, data->mygrid->AxisComm, &sendRequest)); - MPI_SAFE_CALL(MPI_Recv_init(bufferRecv.data(), bufferSize, realMPI, procRecv, + MPI_SAFE_CALL(MPI_Recv_init(bufferRecv.data(), bufferRecv.Size(), realMPI, procRecv, 650, data->mygrid->AxisComm, &recvRequest)); #endif diff --git a/src/fluid/boundary/axis.hpp b/src/fluid/boundary/axis.hpp index 229c772f6..390a628ce 100644 --- a/src/fluid/boundary/axis.hpp +++ b/src/fluid/boundary/axis.hpp @@ -11,6 +11,7 @@ #include #include "idefix.hpp" #include "grid.hpp" +#include "buffer.hpp" // Forward class hydro declaration #include "physics.hpp" @@ -59,8 +60,8 @@ class Axis { MPI_Request sendRequest; MPI_Request recvRequest; - IdefixArray1D bufferSend; - IdefixArray1D bufferRecv; + Buffer bufferSend; + Buffer bufferRecv; int bufferSize; diff --git a/src/fluid/constrainedTransport/EMFexchange.hpp b/src/fluid/constrainedTransport/EMFexchange.hpp index f6e90f75a..aad123346 100644 --- a/src/fluid/constrainedTransport/EMFexchange.hpp +++ b/src/fluid/constrainedTransport/EMFexchange.hpp @@ -30,12 +30,10 @@ void ConstrainedTransport::ExchangeX1(IdefixArray3D ey, IdefixArray3 // Load the buffers with data - int ileft,iright,jbeg,jend,kbeg,kend; - int ny; - [[maybe_unused]] int nz; + int ileft,iright,jbeg,jend; - IdefixArray1D BufferLeft=BufferSendX1[faceLeft]; - IdefixArray1D BufferRight=BufferSendX1[faceRight]; + Buffer BufferLeft=BufferSendX1[faceLeft]; + Buffer BufferRight=BufferSendX1[faceRight]; // If MPI Persistent, start receiving even before the buffers are filled @@ -54,24 +52,30 @@ void ConstrainedTransport::ExchangeX1(IdefixArray3D ey, IdefixArray3 iright = data->end[IDIR]; jbeg = data->beg[JDIR]; jend = data->end[JDIR]; - ny = jend - jbeg; - kbeg = data->beg[KDIR]; - kend = data->end[KDIR]; - nz = kend - kbeg; - - idefix_for("LoadBufferX1Emfz",kbeg,kend,jbeg,jend+1, - KOKKOS_LAMBDA (int k, int j) { - BufferRight( (j-jbeg) + (k-kbeg)*(ny+1) ) = ez(k,j,iright); - } - ); + + // Create the base box to be patched by the ops + BoundingBox baseBox; + baseBox[IDIR][0] = data->beg[IDIR]; + baseBox[IDIR][1] = data->end[IDIR]; + baseBox[JDIR][0] = data->beg[JDIR]; + baseBox[JDIR][1] = data->end[JDIR]; + baseBox[KDIR][0] = data->beg[KDIR]; + baseBox[KDIR][1] = data->end[KDIR]; + + //extend by one the end on jdir && take the ghost on i + BoundingBox sendBoxEz = baseBox; + sendBoxEz[JDIR][1] += 1; + sendBoxEz[IDIR][0] = iright; + sendBoxEz[IDIR][1] = iright + 1; + BufferRight.Pack(ez, sendBoxEz); #if DIMENSIONS == 3 - int Vsindex = (ny+1)*nz; - idefix_for("LoadBufferX1Emfy",kbeg,kend+1,jbeg,jend, - KOKKOS_LAMBDA (int k, int j) { - BufferRight( (j-jbeg) + (k-kbeg)*ny + Vsindex ) = ey(k,j,iright); - } - ); + //extend by one the end on kdir && take the ghost on i + BoundingBox sendBoxEy = baseBox; + sendBoxEy[KDIR][1] += 1; + sendBoxEy[IDIR][0] = iright; + sendBoxEy[IDIR][1] = iright + 1; + BufferRight.Pack(ez, sendBoxEy); #endif // Wait for completion before sending out everything @@ -90,22 +94,22 @@ void ConstrainedTransport::ExchangeX1(IdefixArray3D ey, IdefixArray3 BufferRight=BufferRecvX1[faceRight]; // Erase the emf with the one coming from the left process - - idefix_for("StoreBufferX1Emfz",kbeg,kend,jbeg,jend+1, - KOKKOS_LAMBDA (int k, int j) { - if(lbound == internal || lbound == periodic) { - ez(k,j,ileft) = BufferLeft( (j-jbeg) + (k-kbeg)*(ny+1)); - } - }); + //extend by one the end on jdir && take the ghost revc zone on i + BoundingBox recvBoxEz = baseBox; + recvBoxEz[JDIR][1] += 1; + recvBoxEz[IDIR][0] = ileft; + recvBoxEz[IDIR][1] = ileft + 1; + if(lbound == internal || lbound == periodic) + BufferLeft.Unpack(ez, recvBoxEz); #if DIMENSIONS == 3 - Vsindex = (ny+1)*nz; - idefix_for("StoreBufferX1Emfy",kbeg,kend+1,jbeg,jend, - KOKKOS_LAMBDA (int k, int j) { - if(lbound == internal || lbound == periodic) { - ey(k,j,ileft) = BufferLeft( (j-jbeg) + (k-kbeg)*ny +Vsindex); - } - }); + //extend by one the end on kdir && take the ghost revc zone on i + BoundingBox recvBoxEy = baseBox; + recvBoxEy[KDIR][1] += 1; + recvBoxEy[IDIR][0] = ileft; + recvBoxEy[IDIR][1] = ileft + 1; + if(lbound == internal || lbound == periodic) + BufferLeft.Unpack(ey, recvBoxEy); #endif @@ -118,11 +122,9 @@ void ConstrainedTransport::ExchangeX2(IdefixArray3D ex, IdefixArray3 idfx::pushRegion("Emf::ExchangeX2"); // Load the buffers with data - int jleft,jright,ibeg,iend,kbeg,kend; - int nx; - [[maybe_unused]] int nz; - IdefixArray1D BufferLeft=BufferSendX2[faceLeft]; - IdefixArray1D BufferRight=BufferSendX2[faceRight]; + int jleft,jright; + Buffer BufferLeft=BufferSendX2[faceLeft]; + Buffer BufferRight=BufferSendX2[faceRight]; // If MPI Persistent, start receiving even before the buffers are filled double tStart = MPI_Wtime(); @@ -135,28 +137,32 @@ void ConstrainedTransport::ExchangeX2(IdefixArray3D ex, IdefixArray3 BoundaryType rbound = data->rbound[JDIR]; // Coordinates of the ghost region which needs to be transfered - ibeg = data->beg[IDIR]; - iend = data->end[IDIR]; - nx = iend - ibeg; jleft = data->beg[JDIR]; jright = data->end[JDIR]; - kbeg = data->beg[KDIR]; - kend = data->end[KDIR]; - nz = kend - kbeg; - - idefix_for("LoadBufferX2Emfz",kbeg,kend,ibeg,iend+1, - KOKKOS_LAMBDA (int k, int i) { - BufferRight( (i-ibeg) + (k-kbeg)*(nx+1) ) = ez(k,jright,i); - } - ); + + // Create the base box to be patched by the ops + BoundingBox baseBox; + baseBox[IDIR][0] = data->beg[IDIR]; + baseBox[IDIR][1] = data->end[IDIR]; + baseBox[JDIR][0] = data->beg[JDIR]; + baseBox[JDIR][1] = data->end[JDIR]; + baseBox[KDIR][0] = data->beg[KDIR]; + baseBox[KDIR][1] = data->end[KDIR]; + + //extend by one the end on idir && take the ghost on j + BoundingBox sendBoxEz = baseBox; + sendBoxEz[IDIR][1] += 1; + sendBoxEz[JDIR][0] = jright; + sendBoxEz[JDIR][1] = jright + 1; + BufferRight.Pack(ez, sendBoxEz); #if DIMENSIONS == 3 - int Vsindex = (nx+1)*nz; - idefix_for("LoadBufferX1Emfx",kbeg,kend+1,ibeg,iend, - KOKKOS_LAMBDA (int k, int i) { - BufferRight( (i-ibeg) + (k-kbeg)*nx + Vsindex ) = ex(k,jright,i); - } - ); + //extend by one the end on idir && take the ghost on j + BoundingBox sendBoxEx = baseBox; + sendBoxEx[KDIR][1] += 1; + sendBoxEx[JDIR][0] = jright; + sendBoxEx[JDIR][1] = jright + 1; + BufferRight.Pack(ex, sendBoxEx); #endif // Wait for completion before sending out everything @@ -165,7 +171,7 @@ void ConstrainedTransport::ExchangeX2(IdefixArray3D ex, IdefixArray3 tStart = MPI_Wtime(); MPI_SAFE_CALL(MPI_Startall(2, sendRequestX2)); // Wait for buffers to be received - MPI_Waitall(2,recvRequestX2,recvStatus); + MPI_Waitall(2, recvRequestX2, recvStatus); MPI_Waitall(2, sendRequestX2, sendStatus); idfx::mpiCallsTimer += MPI_Wtime() - tStart; @@ -174,20 +180,21 @@ void ConstrainedTransport::ExchangeX2(IdefixArray3D ex, IdefixArray3 BufferRight=BufferRecvX2[faceRight]; // We average the edge emfs zones - idefix_for("StoreBufferX2Emfz",kbeg,kend,ibeg,iend+1, - KOKKOS_LAMBDA (int k, int i) { - if(lbound == internal || lbound == periodic) { - ez(k,jleft,i) = BufferLeft( (i-ibeg) + (k-kbeg)*(nx+1) ); - } - }); + //extend by one the end on idir && take the ghost revc zone on j + BoundingBox recvBoxEz = baseBox; + recvBoxEz[IDIR][1] += 1; + recvBoxEz[JDIR][0] = jleft; + recvBoxEz[JDIR][1] = jleft + 1; + if(lbound == internal || lbound == periodic) + BufferLeft.Unpack(ez, recvBoxEz); #if DIMENSIONS == 3 - Vsindex = (nx+1)*nz; - idefix_for("StoreBufferX1Emfy",kbeg,kend+1,ibeg,iend, - KOKKOS_LAMBDA (int k, int i) { - if(lbound == internal || lbound == periodic) { - ex(k,jleft,i) = BufferLeft( (i-ibeg) + (k-kbeg)*nx +Vsindex); - } - }); + //extend by one the end on kdir && take the ghost revc zone on j + BoundingBox recvBoxEx = baseBox; + recvBoxEx[KDIR][1] += 1; + recvBoxEx[JDIR][0] = jleft; + recvBoxEx[JDIR][1] = jleft + 1; + if(lbound == internal || lbound == periodic) + BufferLeft.Unpack(ex, recvBoxEx); #endif @@ -201,13 +208,9 @@ void ConstrainedTransport::ExchangeX3(IdefixArray3D ex, IdefixArray3 // Load the buffers with data - int kleft,kright,ibeg,iend,jbeg,jend; - int nx,ny; - IdefixArray1D BufferLeft=BufferSendX3[faceLeft]; - IdefixArray1D BufferRight=BufferSendX3[faceRight]; - - int Vsindex = 0; - + int kleft,kright,jbeg,jend; + Buffer BufferLeft=BufferSendX3[faceLeft]; + Buffer BufferRight=BufferSendX3[faceRight]; // If MPI Persistent, start receiving even before the buffers are filled double tStart = MPI_Wtime(); @@ -220,29 +223,34 @@ void ConstrainedTransport::ExchangeX3(IdefixArray3D ex, IdefixArray3 BoundaryType rbound = data->rbound[KDIR]; // Coordinates of the ghost region which needs to be transfered - ibeg = data->beg[IDIR]; - iend = data->end[IDIR]; - nx = iend - ibeg; - jbeg = data->beg[JDIR]; jend = data->end[JDIR]; - ny = jend - jbeg; kleft = data->beg[KDIR]; kright = data->end[KDIR]; - idefix_for("LoadBufferX3Emfx",jbeg,jend+1,ibeg,iend, - KOKKOS_LAMBDA (int j, int i) { - BufferRight( (i-ibeg) + (j-jbeg)*nx ) = ex(kright,j,i); - } - ); - Vsindex = nx*(ny+1); - - idefix_for("LoadBufferX3Emfy",jbeg,jend,ibeg,iend+1, - KOKKOS_LAMBDA (int j, int i) { - BufferRight( (i-ibeg) + (j-jbeg)*(nx+1) + Vsindex ) = ey(kright,j,i); - } - ); + // Create the base box to be patched by the ops + BoundingBox baseBox; + baseBox[IDIR][0] = data->beg[IDIR]; + baseBox[IDIR][1] = data->end[IDIR]; + baseBox[JDIR][0] = data->beg[JDIR]; + baseBox[JDIR][1] = data->end[JDIR]; + baseBox[KDIR][0] = data->beg[KDIR]; + baseBox[KDIR][1] = data->end[KDIR]; + + //extend by one the end on jdir && take the ghost on k + BoundingBox sendBoxEz = baseBox; + sendBoxEz[JDIR][1] += 1; + sendBoxEz[KDIR][0] = kright; + sendBoxEz[KDIR][1] = kright + 1; + BufferRight.Pack(ez, sendBoxEz); + + //extend by one the end on idir && take the ghost on k + BoundingBox sendBoxEy = baseBox; + sendBoxEy[IDIR][1] += 1; + sendBoxEy[KDIR][0] = kright; + sendBoxEy[KDIR][1] = kright + 1; + BufferRight.Pack(ey, sendBoxEy); // Wait for completion before sending out everything Kokkos::fence(); @@ -250,7 +258,7 @@ void ConstrainedTransport::ExchangeX3(IdefixArray3D ex, IdefixArray3 tStart = MPI_Wtime(); MPI_SAFE_CALL(MPI_Startall(2, sendRequestX3)); // Wait for buffers to be received - MPI_Waitall(2,recvRequestX3,recvStatus); + MPI_Waitall(2, recvRequestX3, recvStatus); MPI_Waitall(2, sendRequestX3, sendStatus); idfx::mpiCallsTimer += MPI_Wtime() - tStart; @@ -259,20 +267,21 @@ void ConstrainedTransport::ExchangeX3(IdefixArray3D ex, IdefixArray3 BufferRight=BufferRecvX3[faceRight]; // We average the edge emfs zones - idefix_for("StoreBufferX3Emfx",jbeg,jend+1,ibeg,iend, - KOKKOS_LAMBDA (int j, int i) { - if(lbound == internal || lbound == periodic) { - ex(kleft,j,i) = BufferLeft( (i-ibeg) + (j-jbeg)*nx ); - } - }); - - Vsindex = nx*(ny+1); - idefix_for("StoreBufferX3Emfy",jbeg,jend,ibeg,iend+1, - KOKKOS_LAMBDA (int j, int i) { - if(lbound == internal || lbound == periodic) { - ey(kleft,j,i) = BufferLeft( (i-ibeg) + (j-jbeg)*(nx+1) + Vsindex ); - } - }); + //extend by one the end on jdir && take the ghost revc zone on k + BoundingBox recvBoxEx = baseBox; + recvBoxEx[JDIR][1] += 1; + recvBoxEx[KDIR][0] = kleft; + recvBoxEx[KDIR][1] = kleft + 1; + if(lbound == internal || lbound == periodic) + BufferLeft.Unpack(ex, recvBoxEx); + + //extend by one the end on idir && take the ghost revc zone on k + BoundingBox recvBoxEy = baseBox; + recvBoxEy[IDIR][1] += 1; + recvBoxEy[KDIR][0] = kleft; + recvBoxEy[KDIR][1] = kleft + 1; + if(lbound == internal || lbound == periodic) + BufferLeft.Unpack(ey, recvBoxEy); idfx::popRegion(); } diff --git a/src/fluid/constrainedTransport/constrainedTransport.hpp b/src/fluid/constrainedTransport/constrainedTransport.hpp index 24c317c09..ab83b23bc 100644 --- a/src/fluid/constrainedTransport/constrainedTransport.hpp +++ b/src/fluid/constrainedTransport/constrainedTransport.hpp @@ -11,6 +11,7 @@ #include "idefix.hpp" #include "input.hpp" #include "riemannSolver.hpp" +#include "buffer.hpp" // Forward declarations #include "physics.hpp" @@ -122,12 +123,12 @@ class ConstrainedTransport { enum {faceRight, faceLeft}; // Buffers for MPI calls - IdefixArray1D BufferSendX1[2]; - IdefixArray1D BufferSendX2[2]; - IdefixArray1D BufferSendX3[2]; - IdefixArray1D BufferRecvX1[2]; - IdefixArray1D BufferRecvX2[2]; - IdefixArray1D BufferRecvX3[2]; + Buffer BufferSendX1[2]; + Buffer BufferSendX2[2]; + Buffer BufferSendX3[2]; + Buffer BufferRecvX1[2]; + Buffer BufferRecvX2[2]; + Buffer BufferRecvX3[2]; IdefixArray1D mapVars; @@ -286,10 +287,10 @@ ConstrainedTransport::ConstrainedTransport(Input &input, Fluid *hydr bufferSizeX1 += data->np_int[JDIR] * (data->np_int[KDIR]+KOFFSET); #endif - BufferRecvX1[faceLeft ] = IdefixArray1D("EmfRecvX1Left", bufferSizeX1); - BufferRecvX1[faceRight] = IdefixArray1D("EmfRecvX1Right",bufferSizeX1); - BufferSendX1[faceLeft ] = IdefixArray1D("EmfSendX1Left", bufferSizeX1); - BufferSendX1[faceRight] = IdefixArray1D("EmfSendX1Right",bufferSizeX1); + BufferRecvX1[faceLeft ] = Buffer(bufferSizeX1); + BufferRecvX1[faceRight] = Buffer(bufferSizeX1); + BufferSendX1[faceLeft ] = Buffer(bufferSizeX1); + BufferSendX1[faceRight] = Buffer(bufferSizeX1); // Number of cells in X2 boundary condition (only required when problem >2D): #if DIMENSIONS >= 2 @@ -301,10 +302,10 @@ ConstrainedTransport::ConstrainedTransport(Input &input, Fluid *hydr bufferSizeX2 += data->np_int[IDIR] * (data->np_int[KDIR]+KOFFSET); #endif - BufferRecvX2[faceLeft ] = IdefixArray1D("EmfRecvX2Left", bufferSizeX2); - BufferRecvX2[faceRight] = IdefixArray1D("EmfRecvX2Right",bufferSizeX2); - BufferSendX2[faceLeft ] = IdefixArray1D("EmfSendX2Left", bufferSizeX2); - BufferSendX2[faceRight] = IdefixArray1D("EmfSendX2Right",bufferSizeX2); + BufferRecvX2[faceLeft ] = Buffer(bufferSizeX2); + BufferRecvX2[faceRight] = Buffer(bufferSizeX2); + BufferSendX2[faceLeft ] = Buffer(bufferSizeX2); + BufferSendX2[faceRight] = Buffer(bufferSizeX2); #endif // Number of cells in X3 boundary condition (only required when problem is 3D): @@ -314,10 +315,10 @@ ConstrainedTransport::ConstrainedTransport(Input &input, Fluid *hydr // ey bufferSizeX3 += (data->np_int[IDIR]+IOFFSET) * data->np_int[JDIR]; - BufferRecvX3[faceLeft ] = IdefixArray1D("EmfRecvX3Left", bufferSizeX3); - BufferRecvX3[faceRight] = IdefixArray1D("EmfRecvX3Right",bufferSizeX3); - BufferSendX3[faceLeft ] = IdefixArray1D("EmfSendX3Left", bufferSizeX3); - BufferSendX3[faceRight] = IdefixArray1D("EmfSendX3Right",bufferSizeX3); + BufferRecvX3[faceLeft ] = Buffer(bufferSizeX3); + BufferRecvX3[faceRight] = Buffer(bufferSizeX3); + BufferSendX3[faceLeft ] = Buffer(bufferSizeX3); + BufferSendX3[faceRight] = Buffer(bufferSizeX3); #endif // DIMENSIONS // Init persistent MPI communications diff --git a/src/mpi/buffer.hpp b/src/mpi/buffer.hpp index ff4dff146..35164355b 100644 --- a/src/mpi/buffer.hpp +++ b/src/mpi/buffer.hpp @@ -31,6 +31,10 @@ class Buffer { return(array.data()); } + const IdefixArray1D & getArray(void) const { + return this->array; + } + int Size() { return(array.size()); } @@ -159,6 +163,33 @@ class Buffer { this->pointer += ninjnk; } + void UnpackJDirSymetric(IdefixArray4D& out, + const int var, + const int symMultiplier, + BoundingBox box) { + const int ni = box[IDIR][1]-box[IDIR][0]; + const int ninj = (box[JDIR][1]-box[JDIR][0])*ni; + const int ninjnk = (box[KDIR][1]-box[KDIR][0])*ninj; + const int ibeg = box[IDIR][0]; + const int jbeg = box[JDIR][0]; + const int kbeg = box[KDIR][0]; + const int iend = box[IDIR][1]; + const int jend = box[JDIR][1]; + const int kend = box[KDIR][1]; + const int offset = this->pointer; + + auto arr = this->array; + idefix_for("UnLoadBuffer4D_var_sym",kbeg,kend,jbeg,jend,ibeg,iend, + KOKKOS_LAMBDA (int k, int j, int i) { + const int jinverted = jend-(j-jbeg)-1; + const int arrIndex = i-ibeg + (j-jbeg)*ni + (k-kbeg)*ninj + offset; + out(var,k,jinverted,i) = symMultiplier * arr(arrIndex); + }); + + // Update pointer + this->pointer += ninjnk; + } + void Unpack(IdefixArray4D& out, IdefixArray1D& map, BoundingBox box) { @@ -186,6 +217,36 @@ class Buffer { this->pointer += ninjnk*map.size(); } + void UnpackJDirSymetric(IdefixArray4D& out, + IdefixArray1D& map, + IdefixArray1D& SymMap, + BoundingBox box) { + const int ni = box[IDIR][1]-box[IDIR][0]; + const int ninj = (box[JDIR][1]-box[JDIR][0])*ni; + const int ninjnk = (box[KDIR][1]-box[KDIR][0])*ninj; + const int ibeg = box[IDIR][0]; + const int jbeg = box[JDIR][0]; + const int kbeg = box[KDIR][0]; + const int iend = box[IDIR][1]; + const int jend = box[JDIR][1]; + const int kend = box[KDIR][1]; + const int offset = this->pointer; + + auto arr = this->array; + idefix_for("UnLoadBuffer4D_map_sym",0,map.size(), + kbeg,kend, + jbeg,jend, + ibeg,iend, + KOKKOS_LAMBDA (int n, int k, int j, int i) { + const int jinverted = jend-(j-jbeg)-1; + const int arrIndex = i-ibeg + (j-jbeg)*ni + (k-kbeg)*ninj + n*ninjnk + offset; + out(map(n),k,jinverted,i) = SymMap(map(n)) * arr(arrIndex); + }); + + // Update pointer + this->pointer += ninjnk*map.size(); + } + private: size_t pointer;