From 195481f615ce49eb7a788f97aa6d60655a1b7fac Mon Sep 17 00:00:00 2001 From: Mark Berrill Date: Thu, 17 May 2018 12:07:20 -0400 Subject: [PATCH] Moving to a new signed distance function --- analysis/TwoPhase.cpp | 2 +- analysis/analysis.cpp | 4 +- analysis/distance.cpp | 213 ++++++++++++++++++ analysis/distance.h | 40 ++++ analysis/eikonal.cpp | 384 --------------------------------- analysis/eikonal.h | 53 ----- analysis/runAnalysis.cpp | 2 +- analysis/uCT.cpp | 22 +- common/Communication.h | 25 ++- common/Communication.hpp | 157 +++++++++----- tests/BlobIdentifyParallel.cpp | 2 +- tests/ComponentLabel.cpp | 2 +- tests/TestBlobIdentify.cpp | 8 +- tests/TestSegDist.cpp | 156 ++++++-------- tests/lbpm_morph_pp.cpp | 72 +++---- tests/lbpm_morphdrain_pp.cpp | 2 +- tests/lbpm_random_pp.cpp | 2 +- tests/lbpm_refine_pp.cpp | 2 +- tests/lbpm_segmented_pp.cpp | 11 +- tests/testCommunication.cpp | 2 +- 20 files changed, 502 insertions(+), 659 deletions(-) create mode 100644 analysis/distance.cpp create mode 100644 analysis/distance.h delete mode 100644 analysis/eikonal.cpp delete mode 100644 analysis/eikonal.h diff --git a/analysis/TwoPhase.cpp b/analysis/TwoPhase.cpp index 492c295e..b2abda0e 100644 --- a/analysis/TwoPhase.cpp +++ b/analysis/TwoPhase.cpp @@ -232,7 +232,7 @@ void TwoPhase::ColorToSignedDistance(double Beta, DoubleArray &ColorData, Double } } - Eikonal(DistData,TempID,Dm,50); + CalcDist(DistData,TempID,Dm); for (int k=0; k fillData(comm,rank_info,nx,ny,nz,1,1,1,0,1,true,true,true); + fillHalo fillData(comm,rank_info,{nx,ny,nz},{1,1,1},0,1,{true,true,true}); fillData.fill(IDs); // Create a list of all neighbor ranks (excluding self) std::vector neighbors; @@ -420,7 +420,7 @@ static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_ } } // Fill the ghosts - fillHalo fillData2(comm,rank_info,nx,ny,nz,1,1,1,0,1,true,true,true); + fillHalo fillData2(comm,rank_info,{nx,ny,nz},{1,1,1},0,1,{true,true,true}); fillData2.fill(IDs); // Reorder based on size (and compress the id space int N_blobs_global = ReorderBlobIDs2(IDs,N_blobs_tot,ngx,ngy,ngz,comm); diff --git a/analysis/distance.cpp b/analysis/distance.cpp new file mode 100644 index 00000000..c5aa4d1e --- /dev/null +++ b/analysis/distance.cpp @@ -0,0 +1,213 @@ +#include "analysis/distance.h" + + +// Check if we need to recompute distance after updating ghose values +static bool checkUpdate( const Array &d, double dx, double dy, double dz ) +{ + auto s = d.size(); + bool test[3] = { false, false, false }; + // Check x-direction + Vec v1, v2; + for (size_t k=1; k +void CalcDist( Array &Distance, const Array &ID, const Domain &Dm, const std::array& periodic ) +{ + ASSERT( Distance.size() == ID.size() ); + std::array n = { Dm.Nx-2, Dm.Ny-2, Dm.Nz-2 }; + fillHalo fillData( Dm.Comm, Dm.rank_info, n, {1,1,1}, 50, 1, {true,false,false}, periodic ); + Array id(ID.size()); + Array vecDist(Distance.size()); + for (size_t i=0; i &d, const Array &ID0, const Domain &Dm, const std::array& periodic ) +{ + const double dx = 1.0; + const double dy = 1.0; + const double dz = 1.0; + std::array N = { Dm.Nx, Dm.Ny, Dm.Nz }; + std::array n = { Dm.Nx-2, Dm.Ny-2, Dm.Nz-2 }; + // Create ID with ghosts + Array ID(N[0],N[1],N[2]); + fillHalo fillDataID( Dm.Comm, Dm.rank_info, n, {1,1,1}, 50, 1, {true,true,true}, periodic ); + fillDataID.copy( ID0, ID ); + // Fill ghosts with nearest neighbor + for (int k=1; k fillData( Dm.Comm, Dm.rank_info, n, {1,1,1}, 50, 1, {true,false,false}, periodic ); + // Initialize the vector distance + d.fill( Vec( 1e50, 1e50, 1e50 ) ); + const double dx0 = 0.5*dx; + const double dy0 = 0.5*dy; + const double dz0 = 0.5*dz; + //const double dxy0 = 0.25*sqrt( dx*dx + dy*dy ); + //const double dxz0 = 0.25*sqrt( dx*dx + dz*dz ); + //const double dyz0 = 0.25*sqrt( dy*dy + dz*dz ); + //const double dxyz0 = sqrt( dx*dx + dy*dy + dz*dz ); + for (int k=1; k=0; i--) { + auto v1 = d(i,j,k); + auto v2 = d(i+1,j,k); + v2.x -= dx; + if ( v2 < v1 ) + d(i,j,k) = v2; + } + } + } + // Propagate +/- y-direction + for (int k=0; k=0; j--) { + auto v1 = d(i,j,k); + auto v2 = d(i,j+1,k); + v2.y -= dy; + if ( v2 < v1 ) + d(i,j,k) = v2; + } + } + } + // Propagate +/- z-direction + for (int j=0; j=0; k--) { + auto v1 = d(i,j,k); + auto v2 = d(i,j,k+1); + v2.z -= dz; + if ( v2 < v1 ) + d(i,j,k) = v2; + } + } + } + fillData.fill( d ); + bool test = checkUpdate( d, dx, dy, dz ); + test = sumReduce( Dm.Comm, test ); + if ( !test ) + break; + } +} +template void CalcDist( Array&, const Array&, const Domain&, const std::array& ); +template void CalcDist( Array&, const Array&, const Domain&, const std::array& ); + + diff --git a/analysis/distance.h b/analysis/distance.h new file mode 100644 index 00000000..b1bb9af5 --- /dev/null +++ b/analysis/distance.h @@ -0,0 +1,40 @@ +#ifndef Distance_H_INC +#define Distance_H_INC + +#include "common/Domain.h" + + +struct Vec { + double x; + double y; + double z; + inline Vec(): x(0), y(0), z(0) {} + inline Vec( double x_, double y_, double z_ ): x(x_), y(y_), z(z_) {} + inline double norm() const { return sqrt(x*x+y*y+z*z); } +}; +inline bool operator<(const Vec& l, const Vec& r){ return l.x*l.x+l.y*l.y+l.z*l.z < r.x*r.x+r.y*r.y+r.z*r.z; } + + +/*! + * @brief Calculate the distance using a simple method + * @details This routine calculates the vector distance to the nearest domain surface. + * @param[out] Distance Distance function + * @param[in] ID Segmentation id + * @param[in] Dm Domain information + * @param[in] periodic Directions that are periodic + */ +template +void CalcDist( Array &Distance, const Array &ID, const Domain &Dm, const std::array& periodic = {true,true,true} ); + +/*! + * @brief Calculate the distance using a simple method + * @details This routine calculates the vector distance to the nearest domain surface. + * @param[out] Distance Distance function + * @param[in] ID Domain id + * @param[in] Dm Domain information + * @param[in] periodic Directions that are periodic + */ +void CalcVecDist( Array &Distance, const Array &ID, const Domain &Dm, const std::array& periodic ); + + +#endif diff --git a/analysis/eikonal.cpp b/analysis/eikonal.cpp deleted file mode 100644 index 01e00405..00000000 --- a/analysis/eikonal.cpp +++ /dev/null @@ -1,384 +0,0 @@ -#include "analysis/eikonal.h" -#include "analysis/imfilter.h" - - - - -static inline float minmod(float &a, float &b) -{ - float value = a; - if ( a*b < 0.0) - value=0.0; - else if (fabs(a) > fabs(b)) - value = b; - return value; -} - - -static inline double minmod(double &a, double &b){ - - double value; - - value = a; - if ( a*b < 0.0) value=0.0; - else if (fabs(a) > fabs(b)) value = b; - - return value; -} - - -template -static inline MPI_Datatype getType( ); -template<> inline MPI_Datatype getType() { return MPI_FLOAT; } -template<> inline MPI_Datatype getType() { return MPI_DOUBLE; } - - -/****************************************************************** -* Solve the eikonal equation * -******************************************************************/ -template -TYPE Eikonal( Array &Distance, const Array &ID, const Domain &Dm, const int timesteps) -{ - - /* - * This routine converts the data in the Distance array to a signed distance - * by solving the equation df/dt = sign(1-|grad f|), where Distance provides - * the values of f on the mesh associated with domain Dm - * It has been tested with segmented data initialized to values [-1,1] - * and will converge toward the signed distance to the surface bounding the associated phases - * - * Reference: - * Min C (2010) On reinitializing level set functions, Journal of Computational Physics 229 - */ - - int xdim = Dm.Nx-2; - int ydim = Dm.Ny-2; - int zdim = Dm.Nz-2; - fillHalo fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1); - - // Arrays to store the second derivatives - Array Dxx(Dm.Nx,Dm.Ny,Dm.Nz); - Array Dyy(Dm.Nx,Dm.Ny,Dm.Nz); - Array Dzz(Dm.Nx,Dm.Ny,Dm.Nz); - Array sign(ID.size()); - for (size_t i=0; i 0.f) - Dx = Dxp*Dxp; - else - Dx = Dxm*Dxm; - if (Dyp + Dym > 0.f) - Dy = Dyp*Dyp; - else - Dy = Dym*Dym; - if (Dzp + Dzm > 0.f) - Dz = Dzp*Dzp; - else - Dz = Dzm*Dzm; - } - else{ - if (Dxp + Dxm < 0.f) - Dx = Dxp*Dxp; - else - Dx = Dxm*Dxm; - if (Dyp + Dym < 0.f) - Dy = Dyp*Dyp; - else - Dy = Dym*Dym; - if (Dzp + Dzm < 0.f) - Dz = Dzp*Dzp; - else - Dz = Dzm*Dzm; - } - - //Dx = max(Dxp*Dxp,Dxm*Dxm); - //Dy = max(Dyp*Dyp,Dym*Dym); - //Dz = max(Dzp*Dzp,Dzm*Dzm); - - double norm = sqrt(Dx + Dy + Dz); - if (norm > 1.0) - norm = 1.0; - - Distance(i,j,k) += dt*s*(1.0 - norm); - LocalVar += dt*s*(1.0 - norm); - - if (fabs(dt*s*(1.0 - norm)) > LocalMax) - LocalMax = fabs(dt*s*(1.0 - norm)); - } - } - } - - double GlobalMax; - MPI_Allreduce(&LocalVar,&GlobalVar,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&LocalMax,&GlobalMax,1,MPI_DOUBLE,MPI_MAX,Dm.Comm); - GlobalVar /= (Dm.Nx-2)*(Dm.Ny-2)*(Dm.Nz-2)*Dm.nprocx()*Dm.nprocy()*Dm.nprocz(); - count++; - - - if (count%50 == 0 && Dm.rank()==0 ){ - printf("Time=%i, Max variation=%f, Global variation=%f \n",count,GlobalMax,GlobalVar); - fflush(stdout); - } - - if (fabs(GlobalMax) < 1e-5){ - if (Dm.rank()==0) - printf("Exiting with max tolerance of 1e-5 \n"); - count=timesteps; - } - } - return GlobalVar; -} -template float Eikonal( Array&, const Array&, const Domain&, int ); -template double Eikonal( Array&, const Array&, const Domain&, int ); - - -/****************************************************************** -* A fast distance calculation * -******************************************************************/ -bool CalcDist3DIteration( Array &Distance, const Domain & ) -{ - const float sq2 = sqrt(2.0f); - const float sq3 = sqrt(3.0f); - float dist0[27]; - dist0[0] = sq3; dist0[1] = sq2; dist0[2] = sq3; - dist0[3] = sq2; dist0[4] = 1; dist0[5] = sq2; - dist0[6] = sq3; dist0[7] = sq2; dist0[8] = sq3; - dist0[9] = sq2; dist0[10] = 1; dist0[11] = sq2; - dist0[12] = 1; dist0[13] = 0; dist0[14] = 1; - dist0[15] = sq2; dist0[16] = 1; dist0[17] = sq2; - dist0[18] = sq3; dist0[19] = sq2; dist0[20] = sq3; - dist0[21] = sq2; dist0[22] = 1; dist0[23] = sq2; - dist0[24] = sq3; dist0[25] = sq2; dist0[26] = sq3; - bool changed = false; - for (size_t k=1; k &Distance, const Array &ID, const Domain &Dm ) -{ - PROFILE_START("Calc Distance"); - // Initialize the distance to be 0 fore the cells adjacent to the interface - Distance.fill(1e100); - for (size_t k=1; k fillData(Dm.Comm, Dm.rank_info,Dm.Nx,Dm.Ny,Dm.Nz,1,1,1,0,1); - while ( true ) { - // Communicate the halo of values - fillData.fill(Distance); - // The distance of the cell is the minimum of the distance of the neighbors plus the distance to that node - bool changed = CalcDist3DIteration( Distance, Dm ); - changed = sumReduce(Dm.Comm,changed); - if ( !changed ) - break; - } - // Update the sign of the distance - for (size_t i=0; i0 ? 1:-1; - PROFILE_STOP("Calc Distance"); -} - - -/****************************************************************** -* A fast distance calculation * -******************************************************************/ -void CalcDistMultiLevelHelper( Array &Distance, const Domain &Dm ) -{ - size_t ratio = 4; - std::function&)> coarsen = [ratio]( const Array& data ) - { - float tmp = 1e100; - int nx = data.size(0); - int ny = data.size(1); - int nz = data.size(2); - for (int k=0; k fillData(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1); - if ( Nx%ratio==0 && Nx>8 && Ny%ratio==0 && Ny>8 && Nz%ratio==0 && Nz>8 ) { - // Use recursive version - int Nr = std::max(std::max(ratio,ratio),ratio); - // Run Nr iterations, communicate, run Nr iterations - for (int i=0; i dist(Nx,Ny,Nz); - fillData.copy(Distance,dist); - auto db = Dm.getDatabase()->cloneDatabase(); - auto n = db->getVector( "n" ); - db->putVector( "n", { n[0]/ratio, n[1]/ratio, n[2]/ratio } ); - Domain Dm2(db); - Dm2.CommInit(Dm.Comm); - fillHalo fillData2(Dm2.Comm,Dm2.rank_info,Nx/ratio,Ny/ratio,Nz/ratio,1,1,1,0,1); - auto dist2 = dist.coarsen( {ratio,ratio,ratio}, coarsen ); - Array Distance2(Nx/ratio+2,Ny/ratio+2,Nz/ratio+2); - fillData2.copy(dist2,Distance2); - // Solve - CalcDistMultiLevelHelper( Distance2, Dm2 ); - // Interpolate the coarse grid to the fine grid - fillData2.copy(Distance2,dist2); - for (int k=0; k &Distance, const Array &ID, const Domain &Dm ) -{ - PROFILE_START("Calc Distance Multilevel"); - int Nx = Dm.Nx-2; - int Ny = Dm.Ny-2; - int Nz = Dm.Nz-2; - ASSERT(int(Distance.size(0))==Nx+2&&int(Distance.size(1))==Ny+2&&int(Distance.size(2))==Nz+2); - fillHalo fillData(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1); - // Initialize the distance to be 0 fore the cells adjacent to the interface - Distance.fill(1e100); - for (size_t k=1; k0 ? 1:-1; - fillData.fill(Distance); - // Run a quick filter to smooth the data - float sigma = 0.6; - Array H = imfilter::create_filter( { 1 }, "gaussian", &sigma ); - std::vector BC(3,imfilter::BC::replicate); - Distance = imfilter::imfilter_separable( Distance, {H,H,H}, BC ); - PROFILE_STOP("Calc Distance Multilevel"); -} diff --git a/analysis/eikonal.h b/analysis/eikonal.h deleted file mode 100644 index d89a7c0e..00000000 --- a/analysis/eikonal.h +++ /dev/null @@ -1,53 +0,0 @@ -#ifndef Eikonal_H_INC -#define Eikonal_H_INC - -#include "common/Domain.h" - - -/*! - * @brief Calculate the distance solving the Eikonal equation - * @details This routine converts the data in the Distance array to a signed distance - * by solving the equation df/dt = sign*(1-|grad f|), where Distance provides - * the values of f on the mesh associated with domain Dm - * It has been tested with segmented data initialized to values [-1,1] - * and will converge toward the signed distance to the surface bounding the associated phases - * - * Reference: - * Min C (2010) On reinitializing level set functions, Journal of Computational Physics 229 - * - * @param[in/out] Distance Distance function - * @param[in] ID Segmentation id - * @param[in] DM Domain information - * @param[in] timesteps Maximum number of timesteps to process - * @return Returns the global variation - */ -template -TYPE Eikonal( Array &Distance, const Array &ID, const Domain &Dm, int timesteps); - - -/*! - * @brief Calculate the distance using a simple method - * @details This routine calculates the distance using a very simple method working off the segmentation id. - * - * @param[in/out] Distance Distance function - * @param[in] ID Segmentation id - * @param[in] DM Domain information - * @return Returns the global variation - */ -void CalcDist3D( Array &Distance, const Array &ID, const Domain &Dm ); - - -/*! - * @brief Calculate the distance using a multi-level method - * @details This routine calculates the distance using a multi-grid method - * - * @return Returns the number of cubes in the blob - * @param[in/out] Distance Distance function - * @param[in] ID Segmentation id - * @param[in] DM Domain information - * @return Returns the global variation - */ -void CalcDistMultiLevel( Array &Distance, const Array &ID, const Domain &Dm ); - - -#endif diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index 2f76b38c..f28697ff 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -289,7 +289,7 @@ runAnalysis::runAnalysis( std::shared_ptr db, d_ScaLBL_Comm( ScaLBL_Comm ), d_rank_info( rank_info ), d_Map( Map ), - d_fillData(Dm.Comm,Dm.rank_info,Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,1,1,1,0,1) + d_fillData(Dm.Comm,Dm.rank_info,{Dm.Nx-2,Dm.Ny-2,Dm.Nz-2},{1,1,1},0,1) { NULL_USE( pBC ); INSIST( db, "Input database is empty" ); diff --git a/analysis/uCT.cpp b/analysis/uCT.cpp index 4ee879af..172de263 100644 --- a/analysis/uCT.cpp +++ b/analysis/uCT.cpp @@ -1,6 +1,6 @@ #include "analysis/uCT.h" #include "analysis/analysis.h" -#include "analysis/eikonal.h" +#include "analysis/distance.h" #include "analysis/filters.h" #include "analysis/imfilter.h" @@ -149,7 +149,7 @@ void solve( const Array& VOL, Array& Mean, Array& ID, fillFloat.fill( Mean ); segment( Mean, ID, 0.01 ); // Compute the distance using the segmented volume - Eikonal( Dist, ID, Dm, ID.size(0)*nprocx ); + CalcDist( Dist, ID, Dm ); fillFloat.fill(Dist); smooth( VOL, Dist, 2.0, MultiScaleSmooth, fillFloat ); // Compute non-local mean @@ -210,9 +210,7 @@ void refine( const Array& Dist_coarse, //removeDisconnected( ID, Dm ); // Compute the distance using the segmented volume if ( level > 0 ) { - //Eikonal3D( Dist, ID, Dm, ID.size(0)*nprocx ); - //CalcDist3D( Dist, ID, Dm ); - CalcDistMultiLevel( Dist, ID, Dm ); + CalcDist( Dist, ID, Dm ); fillFloat.fill(Dist); } } @@ -232,7 +230,7 @@ void filter_final( Array& ID, Array& Dist, int Ny = Dm.Ny-2; int Nz = Dm.Nz-2; // Calculate the distance - CalcDistMultiLevel( Dist, ID, Dm ); + CalcDist( Dist, ID, Dm ); fillFloat.fill(Dist); // Compute the range to shrink the volume based on the L2 norm of the distance Array Dist0(Nx,Ny,Nz); @@ -256,8 +254,8 @@ void filter_final( Array& ID, Array& Dist, } //Array Dist1 = Dist; //Array Dist2 = Dist; - CalcDistMultiLevel( Dist1, ID1, Dm ); - CalcDistMultiLevel( Dist2, ID2, Dm ); + CalcDist( Dist1, ID1, Dm ); + CalcDist( Dist2, ID2, Dm ); fillFloat.fill(Dist1); fillFloat.fill(Dist2); // Keep those regions that are within dx2 of the new volumes @@ -272,8 +270,8 @@ void filter_final( Array& ID, Array& Dist, } } // Find regions of uncertainty that are entirely contained within another region - fillHalo fillDouble(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1); - fillHalo fillInt(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1); + fillHalo fillDouble(Dm.Comm,Dm.rank_info,{Nx,Ny,Nz},{1,1,1},0,1); + fillHalo fillInt(Dm.Comm,Dm.rank_info,{Nx,Ny,Nz},{1,1,1},0,1); BlobIDArray GlobalBlobID; DoubleArray SignDist(ID.size()); for (size_t i=0; i& ID, Array& Dist, // Perform the final segmentation and update the distance fillFloat.fill(Mean); segment( Mean, ID, 0.01 ); - CalcDistMultiLevel( Dist, ID, Dm ); + CalcDist( Dist, ID, Dm ); fillFloat.fill(Dist); } @@ -356,7 +354,7 @@ void filter_src( const Domain& Dm, Array& src ) int Nx = Dm.Nx-2; int Ny = Dm.Ny-2; int Nz = Dm.Nz-2; - fillHalo fillFloat(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1); + fillHalo fillFloat(Dm.Comm,Dm.rank_info,{Nx,Ny,Nz},{1,1,1},0,1); // Perform a hot-spot filter on the data std::vector BC = { imfilter::BC::replicate, imfilter::BC::replicate, imfilter::BC::replicate }; std::function&)> filter_3D = []( const Array& data ) diff --git a/common/Communication.h b/common/Communication.h index 12241be3..17a6f042 100644 --- a/common/Communication.h +++ b/common/Communication.h @@ -5,6 +5,8 @@ #include "common/Utilities.h" #include "common/Array.h" +#include + // ********** COMMUNICTION ************************************** /* //.............................................................. @@ -43,18 +45,17 @@ public: /*! * @brief Default constructor * @param[in] info Rank and neighbor rank info - * @param[in] nx Number of local cells in the x direction - * @param[in] ny Number of local cells in the y direction - * @param[in] nz Number of local cells in the z direction - * @param[in] ngx Number of ghost cells in the x direction - * @param[in] ngy Number of ghost cells in the y direction - * @param[in] ngz Number of ghost cells in the z direction + * @param[in] n Number of local cells + * @param[in] ng Number of ghost cells * @param[in] tag Initial tag to use for the communication (we will require tag:tag+26) * @param[in] depth Maximum depth to support + * @param[in] fill Fill {faces,edges,corners} + * @param[in] periodic Periodic dimensions */ - fillHalo( MPI_Comm comm, const RankInfoStruct& info, int nx, int ny, int nz, - int ngx, int ngy, int ngz, int tag, int depth, - bool fill_face=true, bool fill_edge=true, bool fill_corner=true ); + fillHalo( MPI_Comm comm, const RankInfoStruct& info, + std::array n, std::array ng, int tag, int depth, + std::array fill = {true,true,true}, + std::array periodic = {true,true,true} ); //! Destructor ~fillHalo( ); @@ -75,15 +76,17 @@ public: private: + MPI_Comm comm; RankInfoStruct info; - int nx, ny, nz, ngx, ngy, ngz, depth; + std::array n, ng; + int depth; bool fill_pattern[3][3][3]; int tag[3][3][3]; int N_send_recv[3][3][3]; TYPE *mem; TYPE *send[3][3][3], *recv[3][3][3]; MPI_Request send_req[3][3][3], recv_req[3][3][3]; - MPI_Comm comm; + size_t N_type; MPI_Datatype datatype; fillHalo(); // Private empty constructor fillHalo(const fillHalo&); // Private copy constructor diff --git a/common/Communication.hpp b/common/Communication.hpp index 26139c66..a8eff9ee 100644 --- a/common/Communication.hpp +++ b/common/Communication.hpp @@ -11,16 +11,30 @@ * Structure to store the rank info * ********************************************************/ template -fillHalo::fillHalo( MPI_Comm comm0, const RankInfoStruct& info0, int nx0, int ny0, int nz0, - int ngx0, int ngy0, int ngz0, int tag0, int depth0, - bool fill_face, bool fill_edge, bool fill_corner ): - info(info0), nx(nx0), ny(ny0), nz(nz0), ngx(ngx0), ngy(ngy0), ngz(ngz0), depth(depth0) +fillHalo::fillHalo( MPI_Comm comm_, const RankInfoStruct& info_, + std::array n_, std::array ng_, int tag0, int depth_, + std::array fill, std::array periodic ): + comm(comm_), info(info_), n(n_), ng(ng_), depth(depth_) { - comm = comm0; - datatype = getMPItype(); + if ( std::is_same() ) { + N_type = 1; + datatype = MPI_DOUBLE; + } else if ( std::is_same() ) { + N_type = 1; + datatype = MPI_FLOAT; + } else if ( sizeof(TYPE)%sizeof(double)==0 ) { + N_type = sizeof(TYPE) / sizeof(double); + datatype = MPI_DOUBLE; + } else if ( sizeof(TYPE)%sizeof(float)==0 ) { + N_type = sizeof(TYPE) / sizeof(float); + datatype = MPI_FLOAT; + } else { + N_type = sizeof(TYPE); + datatype = MPI_BYTE; + } // Set the fill pattern memset(fill_pattern,0,sizeof(fill_pattern)); - if ( fill_face ) { + if ( fill[0] ) { fill_pattern[0][1][1] = true; fill_pattern[2][1][1] = true; fill_pattern[1][0][1] = true; @@ -28,7 +42,7 @@ fillHalo::fillHalo( MPI_Comm comm0, const RankInfoStruct& info0, int nx0, fill_pattern[1][1][0] = true; fill_pattern[1][1][2] = true; } - if ( fill_edge ) { + if ( fill[1] ) { fill_pattern[0][0][1] = true; fill_pattern[0][2][1] = true; fill_pattern[2][0][1] = true; @@ -42,7 +56,7 @@ fillHalo::fillHalo( MPI_Comm comm0, const RankInfoStruct& info0, int nx0, fill_pattern[1][2][0] = true; fill_pattern[1][2][2] = true; } - if ( fill_corner ) { + if ( fill[2] ) { fill_pattern[0][0][0] = true; fill_pattern[0][0][2] = true; fill_pattern[0][2][0] = true; @@ -52,13 +66,50 @@ fillHalo::fillHalo( MPI_Comm comm0, const RankInfoStruct& info0, int nx0, fill_pattern[2][2][0] = true; fill_pattern[2][2][2] = true; } + // Remove communication for non-perioidic directions + if ( !periodic[0] && info.ix==0 ) { + for (int j=0; j<3; j++) { + for (int k=0; k<3; k++) + fill_pattern[0][j][k] = false; + } + } + if ( !periodic[0] && info.ix==info.nx-1 ) { + for (int j=0; j<3; j++) { + for (int k=0; k<3; k++) + fill_pattern[2][j][k] = false; + } + } + if ( !periodic[1] && info.jy==0 ) { + for (int i=0; i<3; i++) { + for (int k=0; k<3; k++) + fill_pattern[i][0][k] = false; + } + } + if ( !periodic[1] && info.jy==info.ny-1 ) { + for (int i=0; i<3; i++) { + for (int k=0; k<3; k++) + fill_pattern[i][2][k] = false; + } + } + if ( !periodic[2] && info.kz==0 ) { + for (int i=0; i<3; i++) { + for (int j=0; j<3; j++) + fill_pattern[i][j][0] = false; + } + } + if ( !periodic[2] && info.kz==info.nz-1 ) { + for (int i=0; i<3; i++) { + for (int j=0; j<3; j++) + fill_pattern[i][j][2] = false; + } + } // Determine the number of elements for each send/recv for (int i=0; i<3; i++) { - int ni = (i-1)==0 ? nx:ngx; + int ni = (i-1)==0 ? n[0]:ng[0]; for (int j=0; j<3; j++) { - int nj = (j-1)==0 ? ny:ngy; + int nj = (j-1)==0 ? n[1]:ng[1]; for (int k=0; k<3; k++) { - int nk = (k-1)==0 ? nz:ngz; + int nk = (k-1)==0 ? n[2]:ng[2]; if ( fill_pattern[i][j][k] ) N_send_recv[i][j][k] = ni*nj*nk; else @@ -106,9 +157,9 @@ void fillHalo::fill( Array& data ) { //PROFILE_START("fillHalo::fill",1); int depth2 = data.size(3); - ASSERT((int)data.size(0)==nx+2*ngx); - ASSERT((int)data.size(1)==ny+2*ngy); - ASSERT((int)data.size(2)==nz+2*ngz); + ASSERT((int)data.size(0)==n[0]+2*ng[0]); + ASSERT((int)data.size(1)==n[1]+2*ng[1]); + ASSERT((int)data.size(2)==n[2]+2*ng[2]); ASSERT(depth2<=depth); ASSERT(data.ndim()==3||data.ndim()==4); // Start the recieves @@ -117,7 +168,7 @@ void fillHalo::fill( Array& data ) for (int k=0; k<3; k++) { if ( !fill_pattern[i][j][k] ) continue; - MPI_Irecv( recv[i][j][k], depth2*N_send_recv[i][j][k], datatype, + MPI_Irecv( recv[i][j][k], N_type*depth2*N_send_recv[i][j][k], datatype, info.rank[i][j][k], tag[2-i][2-j][2-k], comm, &recv_req[i][j][k] ); } } @@ -129,7 +180,7 @@ void fillHalo::fill( Array& data ) if ( !fill_pattern[i][j][k] ) continue; pack( data, i-1, j-1, k-1, send[i][j][k] ); - MPI_Isend( send[i][j][k], depth2*N_send_recv[i][j][k], datatype, + MPI_Isend( send[i][j][k], N_type*depth2*N_send_recv[i][j][k], datatype, info.rank[i][j][k], tag[i][j][k], comm, &send_req[i][j][k] ); } } @@ -162,12 +213,12 @@ template void fillHalo::pack( const Array& data, int i0, int j0, int k0, TYPE *buffer ) { int depth2 = data.size(3); - int ni = i0==0 ? nx:ngx; - int nj = j0==0 ? ny:ngy; - int nk = k0==0 ? nz:ngz; - int is = i0==0 ? ngx:((i0==-1)?ngx:nx); - int js = j0==0 ? ngy:((j0==-1)?ngy:ny); - int ks = k0==0 ? ngz:((k0==-1)?ngz:nz); + int ni = i0==0 ? n[0]:ng[0]; + int nj = j0==0 ? n[1]:ng[1]; + int nk = k0==0 ? n[2]:ng[2]; + int is = i0==0 ? ng[0]:((i0==-1)?ng[0]:n[0]); + int js = j0==0 ? ng[1]:((j0==-1)?ng[1]:n[1]); + int ks = k0==0 ? ng[2]:((k0==-1)?ng[2]:n[2]); for (int d=0; d void fillHalo::unpack( Array& data, int i0, int j0, int k0, const TYPE *buffer ) { int depth2 = data.size(3); - int ni = i0==0 ? nx:ngx; - int nj = j0==0 ? ny:ngy; - int nk = k0==0 ? nz:ngz; - int is = i0==0 ? ngx:((i0==-1)?0:nx+ngx); - int js = j0==0 ? ngy:((j0==-1)?0:ny+ngy); - int ks = k0==0 ? ngz:((k0==-1)?0:nz+ngz); + int ni = i0==0 ? n[0]:ng[0]; + int nj = j0==0 ? n[1]:ng[1]; + int nk = k0==0 ? n[2]:ng[2]; + int is = i0==0 ? ng[0]:((i0==-1)?0:n[0]+ng[0]); + int js = j0==0 ? ng[1]:((j0==-1)?0:n[1]+ng[1]); + int ks = k0==0 ? ng[2]:((k0==-1)?0:n[2]+ng[2]); for (int d=0; d void fillHalo::copy( const Array& src, Array& dst ) { //PROFILE_START("fillHalo::copy",1); - ASSERT( (int)src.size(0)==nx || (int)src.size(0)==nx+2*ngx ); - ASSERT( (int)dst.size(0)==nx || (int)dst.size(0)==nx+2*ngx ); - bool src_halo = (int)src.size(0)==nx+2*ngx; - bool dst_halo = (int)dst.size(0)==nx+2*ngx; + ASSERT( (int)src.size(0)==n[0] || (int)src.size(0)==n[0]+2*ng[0] ); + ASSERT( (int)dst.size(0)==n[0] || (int)dst.size(0)==n[0]+2*ng[0] ); + bool src_halo = (int)src.size(0)==n[0]+2*ng[0]; + bool dst_halo = (int)dst.size(0)==n[0]+2*ng[0]; if ( src_halo ) { - ASSERT((int)src.size(0)==nx+2*ngx); - ASSERT((int)src.size(1)==ny+2*ngy); - ASSERT((int)src.size(2)==nz+2*ngz); + ASSERT((int)src.size(0)==n[0]+2*ng[0]); + ASSERT((int)src.size(1)==n[1]+2*ng[1]); + ASSERT((int)src.size(2)==n[2]+2*ng[2]); } else { - ASSERT((int)src.size(0)==nx); - ASSERT((int)src.size(1)==ny); - ASSERT((int)src.size(2)==nz); + ASSERT((int)src.size(0)==n[0]); + ASSERT((int)src.size(1)==n[1]); + ASSERT((int)src.size(2)==n[2]); } if ( dst_halo ) { - ASSERT((int)dst.size(0)==nx+2*ngx); - ASSERT((int)dst.size(1)==ny+2*ngy); - ASSERT((int)dst.size(2)==nz+2*ngz); + ASSERT((int)dst.size(0)==n[0]+2*ng[0]); + ASSERT((int)dst.size(1)==n[1]+2*ng[1]); + ASSERT((int)dst.size(2)==n[2]+2*ng[2]); } else { - ASSERT((int)dst.size(0)==nx); - ASSERT((int)dst.size(1)==ny); - ASSERT((int)dst.size(2)==nz); + ASSERT((int)dst.size(0)==n[0]); + ASSERT((int)dst.size(1)==n[1]); + ASSERT((int)dst.size(2)==n[2]); } if ( src_halo == dst_halo ) { // Src and dst halos match @@ -236,19 +287,19 @@ void fillHalo::copy( const Array& src, Array& dst ) dst(i) = src(i); } else if ( src_halo && !dst_halo ) { // Src has halos - for (int k=0; k fillData(comm,rank_info,nx,ny,nz,1,1,1,0,1); + fillHalo fillData(comm,rank_info,{nx,ny,nz},{1,1,1},0,1); fillData.fill(Phase); fillData.fill(SignDist); diff --git a/tests/ComponentLabel.cpp b/tests/ComponentLabel.cpp index 16587699..f350a5be 100644 --- a/tests/ComponentLabel.cpp +++ b/tests/ComponentLabel.cpp @@ -431,7 +431,7 @@ int main(int argc, char **argv) Averages.PrintComponents(timestep); // Create the MeshDataStruct - fillHalo fillData(Dm.Comm,Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1); + fillHalo fillData(Dm.Comm,Dm.rank_info,{Nx-2,Ny-2,Nz-2},{1,1,1},0,1); std::vector meshData(1); meshData[0].meshName = "domain"; meshData[0].mesh = std::shared_ptr( new IO::DomainMesh(Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) ); diff --git a/tests/TestBlobIdentify.cpp b/tests/TestBlobIdentify.cpp index 71e6f1cf..ccfc6afc 100644 --- a/tests/TestBlobIdentify.cpp +++ b/tests/TestBlobIdentify.cpp @@ -133,8 +133,8 @@ void shift_data( DoubleArray& data, int sx, int sy, int sz, const RankInfoStruct int ngy = ny+2*abs(sy); int ngz = nz+2*abs(sz); Array tmp1(nx,ny,nz), tmp2(ngx,ngy,ngz), tmp3(ngx,ngy,ngz); - fillHalo fillData1(comm,rank_info,nx,ny,nz,1,1,1,0,1); - fillHalo fillData2(comm,rank_info,nx,ny,nz,abs(sx),abs(sy),abs(sz),0,1); + fillHalo fillData1(comm,rank_info,{nx,ny,nz},{1,1,1},0,1); + fillHalo fillData2(comm,rank_info,{nx,ny,nz},{abs(sx),abs(sy),abs(sz)},0,1); fillData1.copy(data,tmp1); fillData2.copy(tmp1,tmp2); fillData2.fill(tmp2); @@ -189,7 +189,7 @@ int main(int argc, char **argv) fillBubbleData( bubbles, Phase, SignDist, Lx, Ly, Lz, rank_info ); // Communication the halos - fillHalo fillData(comm,rank_info,nx,ny,nz,1,1,1,0,1); + fillHalo fillData(comm,rank_info,{nx,ny,nz},{1,1,1},0,1); fillData.fill(Phase); fillData.fill(SignDist); @@ -205,7 +205,7 @@ int main(int argc, char **argv) // Create the MeshDataStruct std::vector meshData(1); meshData[0].meshName = "domain"; - meshData[0].mesh = std::shared_ptr( new IO::DomainMesh(rank_info,nx,ny,nz,Lx,Ly,Lz) ); + meshData[0].mesh = std::make_shared(rank_info,nx,ny,nz,Lx,Ly,Lz); std::shared_ptr PhaseVar( new IO::Variable() ); std::shared_ptr SignDistVar( new IO::Variable() ); std::shared_ptr BlobIDVar( new IO::Variable() ); diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index 18405907..f2b8b631 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -11,8 +11,8 @@ #include #include "common/Array.h" #include "common/Domain.h" -#include "analysis/eikonal.h" #include "IO/Writer.h" +#include "analysis/distance.h" std::shared_ptr loadInputs( int nprocs ) @@ -21,7 +21,7 @@ std::shared_ptr loadInputs( int nprocs ) auto db = std::make_shared( ); db->putScalar( "BC", 0 ); db->putVector( "nproc", { 2, 2, 2 } ); - db->putVector( "n", { 50, 50, 50 } ); + db->putVector( "n", { 100, 100, 100 } ); db->putScalar( "nspheres", 0 ); db->putVector( "L", { 1, 1, 1 } ); return db; @@ -31,13 +31,13 @@ std::shared_ptr loadInputs( int nprocs ) //*************************************************************************************** int main(int argc, char **argv) { - // Initialize MPI - int rank, nprocs; - MPI_Init(&argc,&argv); + // Initialize MPI + int rank, nprocs; + MPI_Init(&argc,&argv); MPI_Comm comm = MPI_COMM_WORLD; - MPI_Comm_rank(comm,&rank); - MPI_Comm_size(comm,&nprocs); - { + MPI_Comm_rank(comm,&rank); + MPI_Comm_size(comm,&nprocs); + { // Load inputs @@ -49,99 +49,81 @@ int main(int argc, char **argv) // Get the rank info Domain Dm(db); - for (int k=0;k id(nx,ny,nz); id.fill(0); - for (int k=1;k ID0(id.size()); @@ -149,17 +131,17 @@ int main(int argc, char **argv) Array ID(Nx,Ny,Nz); Array dist1(Nx,Ny,Nz); Array dist2(Nx,Ny,Nz); - fillHalo fillData(Dm.Comm, Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1); + fillHalo fillData(Dm.Comm, Dm.rank_info,{Nx,Ny,Nz},{1,1,1},0,1); fillData.copy( ID0, ID ); - fillData.copy( Distance, dist1 ); - fillData.copy( TrueDist, dist2 ); + fillData.copy( TrueDist, dist1 ); + fillData.copy( Distance, dist2 ); std::vector data(1); data[0].meshName = "mesh"; data[0].mesh.reset( new IO::DomainMesh( Dm.rank_info, Nx, Ny, Nz, Dm.Lx, Dm.Ly, Dm.Lz ) ); data[0].vars.emplace_back( new IO::Variable( 1, IO::VariableType::VolumeVariable, "ID", ID ) ); - data[0].vars.emplace_back( new IO::Variable( 1, IO::VariableType::VolumeVariable, "Distance", dist1 ) ); - data[0].vars.emplace_back( new IO::Variable( 1, IO::VariableType::VolumeVariable, "TrueDist", dist2 ) ); - data[0].vars.emplace_back( new IO::Variable( 1, IO::VariableType::VolumeVariable, "error", dist1-dist2 ) ); + data[0].vars.emplace_back( new IO::Variable( 1, IO::VariableType::VolumeVariable, "TrueDist", dist1 ) ); + data[0].vars.emplace_back( new IO::Variable( 1, IO::VariableType::VolumeVariable, "Distance", dist2 ) ); + data[0].vars.emplace_back( new IO::Variable( 1, IO::VariableType::VolumeVariable, "error", dist2-dist1 ) ); IO::initialize( "", "silo", false ); IO::writeData( "testSegDist", data, MPI_COMM_WORLD ); diff --git a/tests/lbpm_morph_pp.cpp b/tests/lbpm_morph_pp.cpp index 515e3afd..270e0142 100644 --- a/tests/lbpm_morph_pp.cpp +++ b/tests/lbpm_morph_pp.cpp @@ -343,42 +343,42 @@ int main(int argc, char **argv) PackID(Dm.sendList_yZ, Dm.sendCount_yZ ,sendID_yZ, id); PackID(Dm.sendList_YZ, Dm.sendCount_YZ ,sendID_YZ, id); //...................................................................................... - MPI_Sendrecv(sendID_x,Dm.sendCount_x,MPI_CHAR,Dm.rank_x,sendtag, - recvID_X,Dm.recvCount_X,MPI_CHAR,Dm.rank_X,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_X,Dm.sendCount_X,MPI_CHAR,Dm.rank_X,sendtag, - recvID_x,Dm.recvCount_x,MPI_CHAR,Dm.rank_x,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_y,Dm.sendCount_y,MPI_CHAR,Dm.rank_y,sendtag, - recvID_Y,Dm.recvCount_Y,MPI_CHAR,Dm.rank_Y,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_Y,Dm.sendCount_Y,MPI_CHAR,Dm.rank_Y,sendtag, - recvID_y,Dm.recvCount_y,MPI_CHAR,Dm.rank_y,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_z,Dm.sendCount_z,MPI_CHAR,Dm.rank_z,sendtag, - recvID_Z,Dm.recvCount_Z,MPI_CHAR,Dm.rank_Z,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_Z,Dm.sendCount_Z,MPI_CHAR,Dm.rank_Z,sendtag, - recvID_z,Dm.recvCount_z,MPI_CHAR,Dm.rank_z,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_xy,Dm.sendCount_xy,MPI_CHAR,Dm.rank_xy,sendtag, - recvID_XY,Dm.recvCount_XY,MPI_CHAR,Dm.rank_XY,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_XY,Dm.sendCount_XY,MPI_CHAR,Dm.rank_XY,sendtag, - recvID_xy,Dm.recvCount_xy,MPI_CHAR,Dm.rank_xy,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_Xy,Dm.sendCount_Xy,MPI_CHAR,Dm.rank_Xy,sendtag, - recvID_xY,Dm.recvCount_xY,MPI_CHAR,Dm.rank_xY,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_xY,Dm.sendCount_xY,MPI_CHAR,Dm.rank_xY,sendtag, - recvID_Xy,Dm.recvCount_Xy,MPI_CHAR,Dm.rank_Xy,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_xz,Dm.sendCount_xz,MPI_CHAR,Dm.rank_xz,sendtag, - recvID_XZ,Dm.recvCount_XZ,MPI_CHAR,Dm.rank_XZ,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_XZ,Dm.sendCount_XZ,MPI_CHAR,Dm.rank_XZ,sendtag, - recvID_xz,Dm.recvCount_xz,MPI_CHAR,Dm.rank_xz,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_Xz,Dm.sendCount_Xz,MPI_CHAR,Dm.rank_Xz,sendtag, - recvID_xZ,Dm.recvCount_xZ,MPI_CHAR,Dm.rank_xZ,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_xZ,Dm.sendCount_xZ,MPI_CHAR,Dm.rank_xZ,sendtag, - recvID_Xz,Dm.recvCount_Xz,MPI_CHAR,Dm.rank_Xz,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_yz,Dm.sendCount_yz,MPI_CHAR,Dm.rank_yz,sendtag, - recvID_YZ,Dm.recvCount_YZ,MPI_CHAR,Dm.rank_YZ,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_YZ,Dm.sendCount_YZ,MPI_CHAR,Dm.rank_YZ,sendtag, - recvID_yz,Dm.recvCount_yz,MPI_CHAR,Dm.rank_yz,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_Yz,Dm.sendCount_Yz,MPI_CHAR,Dm.rank_Yz,sendtag, - recvID_yZ,Dm.recvCount_yZ,MPI_CHAR,Dm.rank_yZ,recvtag,comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendID_yZ,Dm.sendCount_yZ,MPI_CHAR,Dm.rank_yZ,sendtag, - recvID_Yz,Dm.recvCount_Yz,MPI_CHAR,Dm.rank_Yz,recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_x,Dm.sendCount_x,MPI_CHAR,Dm.rank_x(),sendtag, + recvID_X,Dm.recvCount_X,MPI_CHAR,Dm.rank_X(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_X,Dm.sendCount_X,MPI_CHAR,Dm.rank_X(),sendtag, + recvID_x,Dm.recvCount_x,MPI_CHAR,Dm.rank_x(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_y,Dm.sendCount_y,MPI_CHAR,Dm.rank_y(),sendtag, + recvID_Y,Dm.recvCount_Y,MPI_CHAR,Dm.rank_Y(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_Y,Dm.sendCount_Y,MPI_CHAR,Dm.rank_Y(),sendtag, + recvID_y,Dm.recvCount_y,MPI_CHAR,Dm.rank_y(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_z,Dm.sendCount_z,MPI_CHAR,Dm.rank_z(),sendtag, + recvID_Z,Dm.recvCount_Z,MPI_CHAR,Dm.rank_Z(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_Z,Dm.sendCount_Z,MPI_CHAR,Dm.rank_Z(),sendtag, + recvID_z,Dm.recvCount_z,MPI_CHAR,Dm.rank_z(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_xy,Dm.sendCount_xy,MPI_CHAR,Dm.rank_xy(),sendtag, + recvID_XY,Dm.recvCount_XY,MPI_CHAR,Dm.rank_XY(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_XY,Dm.sendCount_XY,MPI_CHAR,Dm.rank_XY(),sendtag, + recvID_xy,Dm.recvCount_xy,MPI_CHAR,Dm.rank_xy(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_Xy,Dm.sendCount_Xy,MPI_CHAR,Dm.rank_Xy(),sendtag, + recvID_xY,Dm.recvCount_xY,MPI_CHAR,Dm.rank_xY(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_xY,Dm.sendCount_xY,MPI_CHAR,Dm.rank_xY(),sendtag, + recvID_Xy,Dm.recvCount_Xy,MPI_CHAR,Dm.rank_Xy(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_xz,Dm.sendCount_xz,MPI_CHAR,Dm.rank_xz(),sendtag, + recvID_XZ,Dm.recvCount_XZ,MPI_CHAR,Dm.rank_XZ(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_XZ,Dm.sendCount_XZ,MPI_CHAR,Dm.rank_XZ(),sendtag, + recvID_xz,Dm.recvCount_xz,MPI_CHAR,Dm.rank_xz(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_Xz,Dm.sendCount_Xz,MPI_CHAR,Dm.rank_Xz(),sendtag, + recvID_xZ,Dm.recvCount_xZ,MPI_CHAR,Dm.rank_xZ(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_xZ,Dm.sendCount_xZ,MPI_CHAR,Dm.rank_xZ(),sendtag, + recvID_Xz,Dm.recvCount_Xz,MPI_CHAR,Dm.rank_Xz(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_yz,Dm.sendCount_yz,MPI_CHAR,Dm.rank_yz(),sendtag, + recvID_YZ,Dm.recvCount_YZ,MPI_CHAR,Dm.rank_YZ(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_YZ,Dm.sendCount_YZ,MPI_CHAR,Dm.rank_YZ(),sendtag, + recvID_yz,Dm.recvCount_yz,MPI_CHAR,Dm.rank_yz(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_Yz,Dm.sendCount_Yz,MPI_CHAR,Dm.rank_Yz(),sendtag, + recvID_yZ,Dm.recvCount_yZ,MPI_CHAR,Dm.rank_yZ(),recvtag,comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendID_yZ,Dm.sendCount_yZ,MPI_CHAR,Dm.rank_yZ(),sendtag, + recvID_Yz,Dm.recvCount_Yz,MPI_CHAR,Dm.rank_Yz(),recvtag,comm,MPI_STATUS_IGNORE); UnpackID(Dm.recvList_x, Dm.recvCount_x ,recvID_x, id); UnpackID(Dm.recvList_X, Dm.recvCount_X ,recvID_X, id); diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index 9a9f2f20..42250dcf 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -155,7 +155,7 @@ int main(int argc, char **argv) xdim=Dm.Nx-2; ydim=Dm.Ny-2; zdim=Dm.Nz-2; - fillHalo fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1); + fillHalo fillData(Dm.Comm, Dm.rank_info,{xdim,ydim,zdim},{1,1,1},0,1); DoubleArray SignDist(nx,ny,nz); // Read the signed distance from file diff --git a/tests/lbpm_random_pp.cpp b/tests/lbpm_random_pp.cpp index 8426396b..b6b75e0a 100644 --- a/tests/lbpm_random_pp.cpp +++ b/tests/lbpm_random_pp.cpp @@ -367,7 +367,7 @@ int main(int argc, char **argv) MPI_Sendrecv(sendID_xY,Dm.sendCount_xY,MPI_CHAR,Dm.rank_xY(),sendtag, recvID_Xy,Dm.recvCount_Xy,MPI_CHAR,Dm.rank_Xy(),recvtag,comm,MPI_STATUS_IGNORE); MPI_Sendrecv(sendID_xz,Dm.sendCount_xz,MPI_CHAR,Dm.rank_xz(),sendtag, - recvID_XZ,Dm.recvCount_XZ,MPI_CHAR,Dm.rank_XZ(),recvtag,comm,MPI_STATUS_IGNORE); + recvID_XZ,Dm.recvCount_XZ,MPI_CHAR,Dm.rank_XZ(),recvtag,comm,MPI_STATUS_IGNORE); MPI_Sendrecv(sendID_XZ,Dm.sendCount_XZ,MPI_CHAR,Dm.rank_XZ(),sendtag, recvID_xz,Dm.recvCount_xz,MPI_CHAR,Dm.rank_xz(),recvtag,comm,MPI_STATUS_IGNORE); MPI_Sendrecv(sendID_Xz,Dm.sendCount_Xz,MPI_CHAR,Dm.rank_Xz(),sendtag, diff --git a/tests/lbpm_refine_pp.cpp b/tests/lbpm_refine_pp.cpp index 951d708e..01c12464 100644 --- a/tests/lbpm_refine_pp.cpp +++ b/tests/lbpm_refine_pp.cpp @@ -87,7 +87,7 @@ int main(int argc, char **argv) // Communication the halos const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); - fillHalo fillData(comm,rank_info,rnx,rny,rnz,1,1,1,0,1); + fillHalo fillData(comm,rank_info,{rnx,rny,rnz},{1,1,1},0,1); nx+=2; ny+=2; nz+=2; rnx+=2; rny+=2; rnz+=2; diff --git a/tests/lbpm_segmented_pp.cpp b/tests/lbpm_segmented_pp.cpp index 92f3100d..a14b0c32 100644 --- a/tests/lbpm_segmented_pp.cpp +++ b/tests/lbpm_segmented_pp.cpp @@ -12,7 +12,7 @@ #include "common/Array.h" #include "common/Domain.h" #include "analysis/TwoPhase.h" -#include "analysis/eikonal.h" +#include "analysis/distance.h" inline void MeanFilter(DoubleArray &Mesh){ for (int k=1; k<(int)Mesh.size(2)-1; k++){ @@ -251,15 +251,8 @@ int main(int argc, char **argv) } MeanFilter(Averages.SDs); - double LocalVar, TotalVar; if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n"); - int Maxtime=10*max(max(Dm.Nx*Dm.nprocx(),Dm.Ny*Dm.nprocy()),Dm.Nz*Dm.nprocz()); - Maxtime=min(Maxtime,MAXTIME); - LocalVar = Eikonal(Averages.SDs,id,Dm,Maxtime); - - MPI_Allreduce(&LocalVar,&TotalVar,1,MPI_DOUBLE,MPI_SUM,comm); - TotalVar /= nprocs; - if (rank==0) printf("Final variation in signed distance function %f \n",TotalVar); + CalcDist(Averages.SDs,id,Dm); sprintf(LocalRankFilename,"SignDist.%05i",rank); FILE *DIST = fopen(LocalRankFilename,"wb"); diff --git a/tests/testCommunication.cpp b/tests/testCommunication.cpp index 04b83719..57ce0959 100644 --- a/tests/testCommunication.cpp +++ b/tests/testCommunication.cpp @@ -221,7 +221,7 @@ int testHalo( MPI_Comm comm, int nprocx, int nprocy, int nprocz, int depth ) } // Communicate the halo - fillHalo fillData(comm,rank_info,nx,ny,nz,1,1,1,0,depth); + fillHalo fillData(comm,rank_info,{nx,ny,nz},{1,1,1},0,depth); fillData.fill(array); // Check the results