From 27898d890ff83b1341c3d1cc4654a89d96720523 Mon Sep 17 00:00:00 2001 From: Mark Berrill Date: Wed, 16 May 2018 14:00:35 -0400 Subject: [PATCH 1/3] Fixing bug in Domain --- analysis/TwoPhase.cpp | 92 ++++---- analysis/eikonal.cpp | 336 ++++++++------------------- analysis/eikonal.h | 7 +- analysis/runAnalysis.cpp | 2 +- analysis/uCT.cpp | 2 +- common/Domain.cpp | 240 +++++++++---------- common/Domain.h | 152 ++++++------ common/ScaLBL.cpp | 50 ++-- tests/TestColorSquareTube.cpp | 4 +- tests/TestPoiseuille.cpp | 7 +- tests/TestSegDist.cpp | 106 +++++---- tests/TestTorus.cpp | 14 +- tests/lbpm_captube_pp.cpp | 4 +- tests/lbpm_color_macro_simulator.cpp | 8 +- tests/lbpm_dfh_simulator.cpp | 4 +- tests/lbpm_inkbottle_pp.cpp | 4 +- tests/lbpm_morph_pp.cpp | 6 +- tests/lbpm_morphdrain_pp.cpp | 78 +++---- tests/lbpm_morphopen_pp.cpp | 78 +++---- tests/lbpm_random_pp.cpp | 76 +++--- tests/lbpm_segmented_decomp.cpp | 8 +- tests/lbpm_segmented_pp.cpp | 13 +- tests/lbpm_sphere_pp.cpp | 2 - tests/lbpm_squaretube_pp.cpp | 12 +- 24 files changed, 589 insertions(+), 716 deletions(-) diff --git a/analysis/TwoPhase.cpp b/analysis/TwoPhase.cpp index 72acce40..492c295e 100644 --- a/analysis/TwoPhase.cpp +++ b/analysis/TwoPhase.cpp @@ -70,7 +70,7 @@ TwoPhase::TwoPhase(Domain &dm): As_global(0), wwndnw_global(0), wwnsdnwn_global(0), Jwnwwndnw_global(0), dEs(0), dAwn(0), dAns(0) { Nx=dm.Nx; Ny=dm.Ny; Nz=dm.Nz; - Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm.nprocx*Dm.nprocy*Dm.nprocz*1.0; + Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm.nprocx()*Dm.nprocy()*Dm.nprocz()*1.0; TempID = new char[Nx*Ny*Nz]; @@ -135,7 +135,7 @@ TwoPhase::TwoPhase(Domain &dm): Gns_global.resize(6); Gws_global.resize(6); //......................................... - if (Dm.rank==0){ + if (Dm.rank()==0){ TIMELOG = fopen("timelog.tcat","a+"); if (fseek(TIMELOG,0,SEEK_SET) == fseek(TIMELOG,0,SEEK_CUR)) { @@ -165,7 +165,7 @@ TwoPhase::TwoPhase(Domain &dm): } else{ char LocalRankString[8]; - sprintf(LocalRankString,"%05d",Dm.rank); + sprintf(LocalRankString,"%05d",Dm.rank()); char LocalRankFilename[40]; sprintf(LocalRankFilename,"%s%s","timelog.tcat.",LocalRankString); TIMELOG = fopen(LocalRankFilename,"a+"); @@ -418,8 +418,8 @@ void TwoPhase::ComputeLocal() // If external boundary conditions are set, do not average over the inlet kmin=1; kmax=Nz-1; - if (Dm.BoundaryCondition > 0 && Dm.kproc == 0) kmin=4; - if (Dm.BoundaryCondition > 0 && Dm.kproc == Dm.nprocz-1) kmax=Nz-4; + if (Dm.BoundaryCondition > 0 && Dm.kproc() == 0) kmin=4; + if (Dm.BoundaryCondition > 0 && Dm.kproc() == Dm.nprocz()-1) kmax=Nz-4; for (k=kmin; k 0 && Dm.kproc == 0) kmin=4; - if (Dm.BoundaryCondition > 0 && Dm.kproc == Dm.nprocz-1) kmax=Nz-4; + if (Dm.BoundaryCondition > 0 && Dm.kproc() == 0) kmin=4; + if (Dm.BoundaryCondition > 0 && Dm.kproc() == Dm.nprocz()-1) kmax=Nz-4; for (k=kmin; kA.push_back(A); wn_mesh->B.push_back(B); wn_mesh->C.push_back(C); @@ -1070,15 +1070,15 @@ void TwoPhase::WriteSurfaces(int logcount) B = ws_pts(ws_tris(1,r)); C = ws_pts(ws_tris(2,r)); // Remap the points - A.x += 1.0*Dm.iproc*(Nx-2); - A.y += 1.0*Dm.jproc*(Nx-2); - A.z += 1.0*Dm.kproc*(Nx-2); - B.x += 1.0*Dm.iproc*(Nx-2); - B.y += 1.0*Dm.jproc*(Nx-2); - B.z += 1.0*Dm.kproc*(Nx-2); - C.x += 1.0*Dm.iproc*(Nx-2); - C.y += 1.0*Dm.jproc*(Nx-2); - C.z += 1.0*Dm.kproc*(Nx-2); + A.x += 1.0*Dm.iproc()*(Nx-2); + A.y += 1.0*Dm.jproc()*(Nx-2); + A.z += 1.0*Dm.kproc()*(Nx-2); + B.x += 1.0*Dm.iproc()*(Nx-2); + B.y += 1.0*Dm.jproc()*(Nx-2); + B.z += 1.0*Dm.kproc()*(Nx-2); + C.x += 1.0*Dm.iproc()*(Nx-2); + C.y += 1.0*Dm.jproc()*(Nx-2); + C.z += 1.0*Dm.kproc()*(Nx-2); ws_mesh->A.push_back(A); ws_mesh->B.push_back(B); ws_mesh->C.push_back(C); @@ -1088,15 +1088,15 @@ void TwoPhase::WriteSurfaces(int logcount) B = ns_pts(ns_tris(1,r)); C = ns_pts(ns_tris(2,r)); // Remap the points - A.x += 1.0*Dm.iproc*(Nx-2); - A.y += 1.0*Dm.jproc*(Nx-2); - A.z += 1.0*Dm.kproc*(Nx-2); - B.x += 1.0*Dm.iproc*(Nx-2); - B.y += 1.0*Dm.jproc*(Nx-2); - B.z += 1.0*Dm.kproc*(Nx-2); - C.x += 1.0*Dm.iproc*(Nx-2); - C.y += 1.0*Dm.jproc*(Nx-2); - C.z += 1.0*Dm.kproc*(Nx-2); + A.x += 1.0*Dm.iproc()*(Nx-2); + A.y += 1.0*Dm.jproc()*(Nx-2); + A.z += 1.0*Dm.kproc()*(Nx-2); + B.x += 1.0*Dm.iproc()*(Nx-2); + B.y += 1.0*Dm.jproc()*(Nx-2); + B.z += 1.0*Dm.kproc()*(Nx-2); + C.x += 1.0*Dm.iproc()*(Nx-2); + C.y += 1.0*Dm.jproc()*(Nx-2); + C.z += 1.0*Dm.kproc()*(Nx-2); ns_mesh->A.push_back(A); ns_mesh->B.push_back(B); ns_mesh->C.push_back(C); @@ -1224,7 +1224,7 @@ void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT) void TwoPhase::PrintAll(int timestep) { - if (Dm.rank==0){ + if (Dm.rank()==0){ fprintf(TIMELOG,"%i %.5g ",timestep,dEs); // change in surface energy fprintf(TIMELOG,"%.5g %.5g %.5g ",sat_w,paw_global,pan_global); // saturation and pressure fprintf(TIMELOG,"%.5g %.5g %.5g ",awn_global,ans_global,aws_global); // interfacial areas @@ -1277,7 +1277,7 @@ void TwoPhase::PrintAll(int timestep) void TwoPhase::PrintComponents(int timestep) { - if (Dm.rank==0){ + if (Dm.rank()==0){ printf("PRINT %i COMPONENT AVEREAGES: time = %i \n",(int)ComponentAverages_NWP.size(1),timestep); for (int b=0; b 0.0){ diff --git a/analysis/eikonal.cpp b/analysis/eikonal.cpp index 82d46246..01e00405 100644 --- a/analysis/eikonal.cpp +++ b/analysis/eikonal.cpp @@ -27,199 +27,55 @@ static inline double minmod(double &a, double &b){ } +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 * ******************************************************************/ -double Eikonal(DoubleArray &Distance, const char *ID, const Domain &Dm, int timesteps) +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 i,j,k; - double dt=0.1; - double Dx,Dy,Dz; - double Dxp,Dxm,Dyp,Dym,Dzp,Dzm; - double Dxxp,Dxxm,Dyyp,Dyym,Dzzp,Dzzm; - double sign,norm; - double LocalVar,GlobalVar,LocalMax,GlobalMax; - - int xdim,ydim,zdim; - 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); - - // Arrays to store the second derivatives - DoubleArray Dxx(Dm.Nx,Dm.Ny,Dm.Nz); - DoubleArray Dyy(Dm.Nx,Dm.Ny,Dm.Nz); - DoubleArray Dzz(Dm.Nx,Dm.Ny,Dm.Nz); - - int count = 0; - while (count < timesteps){ - - // Communicate the halo of values - fillData.fill(Distance); - - // Compute second order derivatives - for (k=1;k 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); - - norm=sqrt(Dx + Dy + Dz); - if (norm > 1.0) norm=1.0; - - Distance(i,j,k) += dt*sign*(1.0 - norm); - LocalVar += dt*sign*(1.0 - norm); - - if (fabs(dt*sign*(1.0 - norm)) > LocalMax) - LocalMax = fabs(dt*sign*(1.0 - norm)); - } - } - } - - 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; -} - -float Eikonal3D( Array &Distance, const Array &ID, const Domain &Dm, const int timesteps) -{ - PROFILE_START("Eikonal3D"); - /* - * 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 - */ + * 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 i,j,k; - float dt=0.1; - float Dx,Dy,Dz; - float Dxp,Dxm,Dyp,Dym,Dzp,Dzm; - float Dxxp,Dxxm,Dyyp,Dyym,Dzzp,Dzzm; - float sign,norm; - float LocalVar,GlobalVar,LocalMax,GlobalMax; - - int xdim,ydim,zdim; - 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); + 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 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 &Distance, const Array &ID, const Domain &Dm fillData.fill(Dyy); fillData.fill(Dzz); - LocalMax=LocalVar=0.0; + double LocalMax = 0.0; + double LocalVar = 0.0; // Execute the next timestep - // f(n+1) = f(n) + dt*sign(1-|grad f|) - for (k=1;k Dxm) Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp; - else Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm; - - if (Dyp > Dym) Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp; - else Dy = Distance(i,j,k) - Dym + 0.5*Dyym; - - if (Dzp > Dzm) Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp; - else Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm; + double Dx, Dy, Dz; + if ( s < 0.0){ + 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; } else{ - if (Dxp < Dxm) Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp; - else Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm; - - if (Dyp < Dym) Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp; - else Dy = Distance(i,j,k) - Dym + 0.5*Dyym; - - if (Dzp < Dzm) Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp; - else Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm; + 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; } - norm=sqrt(Dx*Dx+Dy*Dy+Dz*Dz); - if (norm > 1.0) norm=1.0; + //Dx = max(Dxp*Dxp,Dxm*Dxm); + //Dy = max(Dyp*Dyp,Dym*Dym); + //Dz = max(Dzp*Dzp,Dzm*Dzm); - Distance(i,j,k) += dt*sign*(1.0 - norm); - LocalVar += dt*sign*(1.0 - norm); + double norm = sqrt(Dx + Dy + Dz); + if (norm > 1.0) + norm = 1.0; - if (fabs(dt*sign*(1.0 - norm)) > LocalMax) - LocalMax = fabs(dt*sign*(1.0 - norm)); + 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)); } } } - MPI_Allreduce(&LocalVar,&GlobalVar,1,MPI_FLOAT,MPI_SUM,Dm.Comm); - MPI_Allreduce(&LocalMax,&GlobalMax,1,MPI_FLOAT,MPI_MAX,Dm.Comm); - GlobalVar /= (Dm.Nx-2)*(Dm.Ny-2)*(Dm.Nz-2)*Dm.nprocx*Dm.nprocy*Dm.nprocz; + 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); + + 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"); + if (Dm.rank()==0) + printf("Exiting with max tolerance of 1e-5 \n"); count=timesteps; } } - PROFILE_STOP("Eikonal3D"); 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 &Dm ) +bool CalcDist3DIteration( Array &Distance, const Domain & ) { const float sq2 = sqrt(2.0f); const float sq3 = sqrt(3.0f); diff --git a/analysis/eikonal.h b/analysis/eikonal.h index 5157e1c5..d89a7c0e 100644 --- a/analysis/eikonal.h +++ b/analysis/eikonal.h @@ -15,23 +15,20 @@ * Reference: * Min C (2010) On reinitializing level set functions, Journal of Computational Physics 229 * - * @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 * @param[in] timesteps Maximum number of timesteps to process * @return Returns the global variation */ -double Eikonal(DoubleArray &Distance, const char *ID, const Domain &Dm, int timesteps); - -float Eikonal3D( Array &Distance, const Array &ID, const Domain &Dm, const int timesteps); +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. * - * @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 diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index 8f9c05e8..2f76b38c 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -294,7 +294,7 @@ runAnalysis::runAnalysis( std::shared_ptr db, NULL_USE( pBC ); INSIST( db, "Input database is empty" ); char rankString[20]; - sprintf(rankString,"%05d",Dm.rank); + sprintf(rankString,"%05d",Dm.rank()); d_N[0] = Dm.Nx; d_N[1] = Dm.Ny; d_N[2] = Dm.Nz; diff --git a/analysis/uCT.cpp b/analysis/uCT.cpp index efa43d47..4ee879af 100644 --- a/analysis/uCT.cpp +++ b/analysis/uCT.cpp @@ -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 - Eikonal3D( Dist, ID, Dm, ID.size(0)*nprocx ); + Eikonal( Dist, ID, Dm, ID.size(0)*nprocx ); fillFloat.fill(Dist); smooth( VOL, Dist, 2.0, MultiScaleSmooth, fillFloat ); // Compute non-local mean diff --git a/common/Domain.cpp b/common/Domain.cpp index 13487c51..66b958ac 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -33,7 +33,7 @@ static inline void fgetl( char * str, int num, FILE * stream ) Domain::Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz, double lx, double ly, double lz, int BC): Nx(0), Ny(0), Nz(0), - Lx(0), Ly(0), Lz(0), Volume(0), rank(0), BoundaryCondition(0), + Lx(0), Ly(0), Lz(0), Volume(0), BoundaryCondition(0), Comm(MPI_COMM_NULL), sendCount_x(0), sendCount_y(0), sendCount_z(0), sendCount_X(0), sendCount_Y(0), sendCount_Z(0), sendCount_xy(0), sendCount_yz(0), sendCount_xz(0), sendCount_Xy(0), sendCount_Yz(0), sendCount_xZ(0), @@ -70,8 +70,8 @@ Domain::Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz, initialize( db ); } Domain::Domain( std::shared_ptr db ): - Nx(0), Ny(0), Nz(0), iproc(0), jproc(0), - Lx(0), Ly(0), Lz(0), Volume(0), rank(0), BoundaryCondition(0), + Nx(0), Ny(0), Nz(0), + Lx(0), Ly(0), Lz(0), Volume(0), BoundaryCondition(0), Comm(MPI_COMM_NULL), sendCount_x(0), sendCount_y(0), sendCount_z(0), sendCount_X(0), sendCount_Y(0), sendCount_Z(0), sendCount_xy(0), sendCount_yz(0), sendCount_xz(0), sendCount_Xy(0), sendCount_Yz(0), sendCount_xZ(0), @@ -126,13 +126,13 @@ void Domain::initialize( std::shared_ptr db ) rank_info = RankInfoStruct(myrank,nproc[0],nproc[1],nproc[2]); // Fill remaining variables N = Nx*Ny*Nz; - Volume = nx*ny*nx*nprocx*nprocy*nprocz*1.0; + Volume = nx*ny*nx*nproc[0]*nproc[1]*nproc[2]*1.0; id = new char[N]; memset(id,0,N); BoundaryCondition = d_db->getScalar("BC"); int nprocs; MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); - INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!"); + INSIST(nprocs == nproc[0]*nproc[1]*nproc[2],"Fatal error in processor count!"); } Domain::~Domain() { @@ -199,7 +199,9 @@ void Domain::CommInit(MPI_Comm Communicator) MPI_Comm_dup(Communicator,&Comm); // set up the neighbor ranks - rank_info = RankInfoStruct( rank, nprocx, nprocy, nprocz ); + int myrank; + MPI_Comm_rank( Comm, &myrank ); + rank_info = RankInfoStruct( myrank, rank_info.nx, rank_info.ny, rank_info.nz ); MPI_Barrier(Communicator); @@ -319,42 +321,42 @@ void Domain::CommInit(MPI_Comm Communicator) sendBuf_YZ = new int [sendCount_YZ]; sendBuf_XZ = new int [sendCount_XZ]; //...................................................................................... - MPI_Isend(&sendCount_x, 1,MPI_INT,rank_x,sendtag+0,Communicator,&req1[0]); - MPI_Irecv(&recvCount_X, 1,MPI_INT,rank_X,recvtag+0,Communicator,&req2[0]); - MPI_Isend(&sendCount_X, 1,MPI_INT,rank_X,sendtag+1,Communicator,&req1[1]); - MPI_Irecv(&recvCount_x, 1,MPI_INT,rank_x,recvtag+1,Communicator,&req2[1]); - MPI_Isend(&sendCount_y, 1,MPI_INT,rank_y,sendtag+2,Communicator,&req1[2]); - MPI_Irecv(&recvCount_Y, 1,MPI_INT,rank_Y,recvtag+2,Communicator,&req2[2]); - MPI_Isend(&sendCount_Y, 1,MPI_INT,rank_Y,sendtag+3,Communicator,&req1[3]); - MPI_Irecv(&recvCount_y, 1,MPI_INT,rank_y,recvtag+3,Communicator,&req2[3]); - MPI_Isend(&sendCount_z, 1,MPI_INT,rank_z,sendtag+4,Communicator,&req1[4]); - MPI_Irecv(&recvCount_Z, 1,MPI_INT,rank_Z,recvtag+4,Communicator,&req2[4]); - MPI_Isend(&sendCount_Z, 1,MPI_INT,rank_Z,sendtag+5,Communicator,&req1[5]); - MPI_Irecv(&recvCount_z, 1,MPI_INT,rank_z,recvtag+5,Communicator,&req2[5]); - MPI_Isend(&sendCount_xy, 1,MPI_INT,rank_xy,sendtag+6,Communicator,&req1[6]); - MPI_Irecv(&recvCount_XY, 1,MPI_INT,rank_XY,recvtag+6,Communicator,&req2[6]); - MPI_Isend(&sendCount_XY, 1,MPI_INT,rank_XY,sendtag+7,Communicator,&req1[7]); - MPI_Irecv(&recvCount_xy, 1,MPI_INT,rank_xy,recvtag+7,Communicator,&req2[7]); - MPI_Isend(&sendCount_Xy, 1,MPI_INT,rank_Xy,sendtag+8,Communicator,&req1[8]); - MPI_Irecv(&recvCount_xY, 1,MPI_INT,rank_xY,recvtag+8,Communicator,&req2[8]); - MPI_Isend(&sendCount_xY, 1,MPI_INT,rank_xY,sendtag+9,Communicator,&req1[9]); - MPI_Irecv(&recvCount_Xy, 1,MPI_INT,rank_Xy,recvtag+9,Communicator,&req2[9]); - MPI_Isend(&sendCount_xz, 1,MPI_INT,rank_xz,sendtag+10,Communicator,&req1[10]); - MPI_Irecv(&recvCount_XZ, 1,MPI_INT,rank_XZ,recvtag+10,Communicator,&req2[10]); - MPI_Isend(&sendCount_XZ, 1,MPI_INT,rank_XZ,sendtag+11,Communicator,&req1[11]); - MPI_Irecv(&recvCount_xz, 1,MPI_INT,rank_xz,recvtag+11,Communicator,&req2[11]); - MPI_Isend(&sendCount_Xz, 1,MPI_INT,rank_Xz,sendtag+12,Communicator,&req1[12]); - MPI_Irecv(&recvCount_xZ, 1,MPI_INT,rank_xZ,recvtag+12,Communicator,&req2[12]); - MPI_Isend(&sendCount_xZ, 1,MPI_INT,rank_xZ,sendtag+13,Communicator,&req1[13]); - MPI_Irecv(&recvCount_Xz, 1,MPI_INT,rank_Xz,recvtag+13,Communicator,&req2[13]); - MPI_Isend(&sendCount_yz, 1,MPI_INT,rank_yz,sendtag+14,Communicator,&req1[14]); - MPI_Irecv(&recvCount_YZ, 1,MPI_INT,rank_YZ,recvtag+14,Communicator,&req2[14]); - MPI_Isend(&sendCount_YZ, 1,MPI_INT,rank_YZ,sendtag+15,Communicator,&req1[15]); - MPI_Irecv(&recvCount_yz, 1,MPI_INT,rank_yz,recvtag+15,Communicator,&req2[15]); - MPI_Isend(&sendCount_Yz, 1,MPI_INT,rank_Yz,sendtag+16,Communicator,&req1[16]); - MPI_Irecv(&recvCount_yZ, 1,MPI_INT,rank_yZ,recvtag+16,Communicator,&req2[16]); - MPI_Isend(&sendCount_yZ, 1,MPI_INT,rank_yZ,sendtag+17,Communicator,&req1[17]); - MPI_Irecv(&recvCount_Yz, 1,MPI_INT,rank_Yz,recvtag+17,Communicator,&req2[17]); + MPI_Isend(&sendCount_x, 1,MPI_INT,rank_x(),sendtag+0,Communicator,&req1[0]); + MPI_Irecv(&recvCount_X, 1,MPI_INT,rank_X(),recvtag+0,Communicator,&req2[0]); + MPI_Isend(&sendCount_X, 1,MPI_INT,rank_X(),sendtag+1,Communicator,&req1[1]); + MPI_Irecv(&recvCount_x, 1,MPI_INT,rank_x(),recvtag+1,Communicator,&req2[1]); + MPI_Isend(&sendCount_y, 1,MPI_INT,rank_y(),sendtag+2,Communicator,&req1[2]); + MPI_Irecv(&recvCount_Y, 1,MPI_INT,rank_Y(),recvtag+2,Communicator,&req2[2]); + MPI_Isend(&sendCount_Y, 1,MPI_INT,rank_Y(),sendtag+3,Communicator,&req1[3]); + MPI_Irecv(&recvCount_y, 1,MPI_INT,rank_y(),recvtag+3,Communicator,&req2[3]); + MPI_Isend(&sendCount_z, 1,MPI_INT,rank_z(),sendtag+4,Communicator,&req1[4]); + MPI_Irecv(&recvCount_Z, 1,MPI_INT,rank_Z(),recvtag+4,Communicator,&req2[4]); + MPI_Isend(&sendCount_Z, 1,MPI_INT,rank_Z(),sendtag+5,Communicator,&req1[5]); + MPI_Irecv(&recvCount_z, 1,MPI_INT,rank_z(),recvtag+5,Communicator,&req2[5]); + MPI_Isend(&sendCount_xy, 1,MPI_INT,rank_xy(),sendtag+6,Communicator,&req1[6]); + MPI_Irecv(&recvCount_XY, 1,MPI_INT,rank_XY(),recvtag+6,Communicator,&req2[6]); + MPI_Isend(&sendCount_XY, 1,MPI_INT,rank_XY(),sendtag+7,Communicator,&req1[7]); + MPI_Irecv(&recvCount_xy, 1,MPI_INT,rank_xy(),recvtag+7,Communicator,&req2[7]); + MPI_Isend(&sendCount_Xy, 1,MPI_INT,rank_Xy(),sendtag+8,Communicator,&req1[8]); + MPI_Irecv(&recvCount_xY, 1,MPI_INT,rank_xY(),recvtag+8,Communicator,&req2[8]); + MPI_Isend(&sendCount_xY, 1,MPI_INT,rank_xY(),sendtag+9,Communicator,&req1[9]); + MPI_Irecv(&recvCount_Xy, 1,MPI_INT,rank_Xy(),recvtag+9,Communicator,&req2[9]); + MPI_Isend(&sendCount_xz, 1,MPI_INT,rank_xz(),sendtag+10,Communicator,&req1[10]); + MPI_Irecv(&recvCount_XZ, 1,MPI_INT,rank_XZ(),recvtag+10,Communicator,&req2[10]); + MPI_Isend(&sendCount_XZ, 1,MPI_INT,rank_XZ(),sendtag+11,Communicator,&req1[11]); + MPI_Irecv(&recvCount_xz, 1,MPI_INT,rank_xz(),recvtag+11,Communicator,&req2[11]); + MPI_Isend(&sendCount_Xz, 1,MPI_INT,rank_Xz(),sendtag+12,Communicator,&req1[12]); + MPI_Irecv(&recvCount_xZ, 1,MPI_INT,rank_xZ(),recvtag+12,Communicator,&req2[12]); + MPI_Isend(&sendCount_xZ, 1,MPI_INT,rank_xZ(),sendtag+13,Communicator,&req1[13]); + MPI_Irecv(&recvCount_Xz, 1,MPI_INT,rank_Xz(),recvtag+13,Communicator,&req2[13]); + MPI_Isend(&sendCount_yz, 1,MPI_INT,rank_yz(),sendtag+14,Communicator,&req1[14]); + MPI_Irecv(&recvCount_YZ, 1,MPI_INT,rank_YZ(),recvtag+14,Communicator,&req2[14]); + MPI_Isend(&sendCount_YZ, 1,MPI_INT,rank_YZ(),sendtag+15,Communicator,&req1[15]); + MPI_Irecv(&recvCount_yz, 1,MPI_INT,rank_yz(),recvtag+15,Communicator,&req2[15]); + MPI_Isend(&sendCount_Yz, 1,MPI_INT,rank_Yz(),sendtag+16,Communicator,&req1[16]); + MPI_Irecv(&recvCount_yZ, 1,MPI_INT,rank_yZ(),recvtag+16,Communicator,&req2[16]); + MPI_Isend(&sendCount_yZ, 1,MPI_INT,rank_yZ(),sendtag+17,Communicator,&req1[17]); + MPI_Irecv(&recvCount_Yz, 1,MPI_INT,rank_Yz(),recvtag+17,Communicator,&req2[17]); MPI_Waitall(18,req1,stat1); MPI_Waitall(18,req2,stat2); MPI_Barrier(Communicator); @@ -379,42 +381,42 @@ void Domain::CommInit(MPI_Comm Communicator) recvList_YZ = new int [recvCount_YZ]; recvList_XZ = new int [recvCount_XZ]; //...................................................................................... - MPI_Isend(sendList_x, sendCount_x,MPI_INT,rank_x,sendtag,Communicator,&req1[0]); - MPI_Irecv(recvList_X, recvCount_X,MPI_INT,rank_X,recvtag,Communicator,&req2[0]); - MPI_Isend(sendList_X, sendCount_X,MPI_INT,rank_X,sendtag,Communicator,&req1[1]); - MPI_Irecv(recvList_x, recvCount_x,MPI_INT,rank_x,recvtag,Communicator,&req2[1]); - MPI_Isend(sendList_y, sendCount_y,MPI_INT,rank_y,sendtag,Communicator,&req1[2]); - MPI_Irecv(recvList_Y, recvCount_Y,MPI_INT,rank_Y,recvtag,Communicator,&req2[2]); - MPI_Isend(sendList_Y, sendCount_Y,MPI_INT,rank_Y,sendtag,Communicator,&req1[3]); - MPI_Irecv(recvList_y, recvCount_y,MPI_INT,rank_y,recvtag,Communicator,&req2[3]); - MPI_Isend(sendList_z, sendCount_z,MPI_INT,rank_z,sendtag,Communicator,&req1[4]); - MPI_Irecv(recvList_Z, recvCount_Z,MPI_INT,rank_Z,recvtag,Communicator,&req2[4]); - MPI_Isend(sendList_Z, sendCount_Z,MPI_INT,rank_Z,sendtag,Communicator,&req1[5]); - MPI_Irecv(recvList_z, recvCount_z,MPI_INT,rank_z,recvtag,Communicator,&req2[5]); - MPI_Isend(sendList_xy, sendCount_xy,MPI_INT,rank_xy,sendtag,Communicator,&req1[6]); - MPI_Irecv(recvList_XY, recvCount_XY,MPI_INT,rank_XY,recvtag,Communicator,&req2[6]); - MPI_Isend(sendList_XY, sendCount_XY,MPI_INT,rank_XY,sendtag,Communicator,&req1[7]); - MPI_Irecv(recvList_xy, recvCount_xy,MPI_INT,rank_xy,recvtag,Communicator,&req2[7]); - MPI_Isend(sendList_Xy, sendCount_Xy,MPI_INT,rank_Xy,sendtag,Communicator,&req1[8]); - MPI_Irecv(recvList_xY, recvCount_xY,MPI_INT,rank_xY,recvtag,Communicator,&req2[8]); - MPI_Isend(sendList_xY, sendCount_xY,MPI_INT,rank_xY,sendtag,Communicator,&req1[9]); - MPI_Irecv(recvList_Xy, recvCount_Xy,MPI_INT,rank_Xy,recvtag,Communicator,&req2[9]); - MPI_Isend(sendList_xz, sendCount_xz,MPI_INT,rank_xz,sendtag,Communicator,&req1[10]); - MPI_Irecv(recvList_XZ, recvCount_XZ,MPI_INT,rank_XZ,recvtag,Communicator,&req2[10]); - MPI_Isend(sendList_XZ, sendCount_XZ,MPI_INT,rank_XZ,sendtag,Communicator,&req1[11]); - MPI_Irecv(recvList_xz, recvCount_xz,MPI_INT,rank_xz,recvtag,Communicator,&req2[11]); - MPI_Isend(sendList_Xz, sendCount_Xz,MPI_INT,rank_Xz,sendtag,Communicator,&req1[12]); - MPI_Irecv(recvList_xZ, recvCount_xZ,MPI_INT,rank_xZ,recvtag,Communicator,&req2[12]); - MPI_Isend(sendList_xZ, sendCount_xZ,MPI_INT,rank_xZ,sendtag,Communicator,&req1[13]); - MPI_Irecv(recvList_Xz, recvCount_Xz,MPI_INT,rank_Xz,recvtag,Communicator,&req2[13]); - MPI_Isend(sendList_yz, sendCount_yz,MPI_INT,rank_yz,sendtag,Communicator,&req1[14]); - MPI_Irecv(recvList_YZ, recvCount_YZ,MPI_INT,rank_YZ,recvtag,Communicator,&req2[14]); - MPI_Isend(sendList_YZ, sendCount_YZ,MPI_INT,rank_YZ,sendtag,Communicator,&req1[15]); - MPI_Irecv(recvList_yz, recvCount_yz,MPI_INT,rank_yz,recvtag,Communicator,&req2[15]); - MPI_Isend(sendList_Yz, sendCount_Yz,MPI_INT,rank_Yz,sendtag,Communicator,&req1[16]); - MPI_Irecv(recvList_yZ, recvCount_yZ,MPI_INT,rank_yZ,recvtag,Communicator,&req2[16]); - MPI_Isend(sendList_yZ, sendCount_yZ,MPI_INT,rank_yZ,sendtag,Communicator,&req1[17]); - MPI_Irecv(recvList_Yz, recvCount_Yz,MPI_INT,rank_Yz,recvtag,Communicator,&req2[17]); + MPI_Isend(sendList_x, sendCount_x,MPI_INT,rank_x(),sendtag,Communicator,&req1[0]); + MPI_Irecv(recvList_X, recvCount_X,MPI_INT,rank_X(),recvtag,Communicator,&req2[0]); + MPI_Isend(sendList_X, sendCount_X,MPI_INT,rank_X(),sendtag,Communicator,&req1[1]); + MPI_Irecv(recvList_x, recvCount_x,MPI_INT,rank_x(),recvtag,Communicator,&req2[1]); + MPI_Isend(sendList_y, sendCount_y,MPI_INT,rank_y(),sendtag,Communicator,&req1[2]); + MPI_Irecv(recvList_Y, recvCount_Y,MPI_INT,rank_Y(),recvtag,Communicator,&req2[2]); + MPI_Isend(sendList_Y, sendCount_Y,MPI_INT,rank_Y(),sendtag,Communicator,&req1[3]); + MPI_Irecv(recvList_y, recvCount_y,MPI_INT,rank_y(),recvtag,Communicator,&req2[3]); + MPI_Isend(sendList_z, sendCount_z,MPI_INT,rank_z(),sendtag,Communicator,&req1[4]); + MPI_Irecv(recvList_Z, recvCount_Z,MPI_INT,rank_Z(),recvtag,Communicator,&req2[4]); + MPI_Isend(sendList_Z, sendCount_Z,MPI_INT,rank_Z(),sendtag,Communicator,&req1[5]); + MPI_Irecv(recvList_z, recvCount_z,MPI_INT,rank_z(),recvtag,Communicator,&req2[5]); + MPI_Isend(sendList_xy, sendCount_xy,MPI_INT,rank_xy(),sendtag,Communicator,&req1[6]); + MPI_Irecv(recvList_XY, recvCount_XY,MPI_INT,rank_XY(),recvtag,Communicator,&req2[6]); + MPI_Isend(sendList_XY, sendCount_XY,MPI_INT,rank_XY(),sendtag,Communicator,&req1[7]); + MPI_Irecv(recvList_xy, recvCount_xy,MPI_INT,rank_xy(),recvtag,Communicator,&req2[7]); + MPI_Isend(sendList_Xy, sendCount_Xy,MPI_INT,rank_Xy(),sendtag,Communicator,&req1[8]); + MPI_Irecv(recvList_xY, recvCount_xY,MPI_INT,rank_xY(),recvtag,Communicator,&req2[8]); + MPI_Isend(sendList_xY, sendCount_xY,MPI_INT,rank_xY(),sendtag,Communicator,&req1[9]); + MPI_Irecv(recvList_Xy, recvCount_Xy,MPI_INT,rank_Xy(),recvtag,Communicator,&req2[9]); + MPI_Isend(sendList_xz, sendCount_xz,MPI_INT,rank_xz(),sendtag,Communicator,&req1[10]); + MPI_Irecv(recvList_XZ, recvCount_XZ,MPI_INT,rank_XZ(),recvtag,Communicator,&req2[10]); + MPI_Isend(sendList_XZ, sendCount_XZ,MPI_INT,rank_XZ(),sendtag,Communicator,&req1[11]); + MPI_Irecv(recvList_xz, recvCount_xz,MPI_INT,rank_xz(),recvtag,Communicator,&req2[11]); + MPI_Isend(sendList_Xz, sendCount_Xz,MPI_INT,rank_Xz(),sendtag,Communicator,&req1[12]); + MPI_Irecv(recvList_xZ, recvCount_xZ,MPI_INT,rank_xZ(),recvtag,Communicator,&req2[12]); + MPI_Isend(sendList_xZ, sendCount_xZ,MPI_INT,rank_xZ(),sendtag,Communicator,&req1[13]); + MPI_Irecv(recvList_Xz, recvCount_Xz,MPI_INT,rank_Xz(),recvtag,Communicator,&req2[13]); + MPI_Isend(sendList_yz, sendCount_yz,MPI_INT,rank_yz(),sendtag,Communicator,&req1[14]); + MPI_Irecv(recvList_YZ, recvCount_YZ,MPI_INT,rank_YZ(),recvtag,Communicator,&req2[14]); + MPI_Isend(sendList_YZ, sendCount_YZ,MPI_INT,rank_YZ(),sendtag,Communicator,&req1[15]); + MPI_Irecv(recvList_yz, recvCount_yz,MPI_INT,rank_yz(),recvtag,Communicator,&req2[15]); + MPI_Isend(sendList_Yz, sendCount_Yz,MPI_INT,rank_Yz(),sendtag,Communicator,&req1[16]); + MPI_Irecv(recvList_yZ, recvCount_yZ,MPI_INT,rank_yZ(),recvtag,Communicator,&req2[16]); + MPI_Isend(sendList_yZ, sendCount_yZ,MPI_INT,rank_yZ(),sendtag,Communicator,&req1[17]); + MPI_Irecv(recvList_Yz, recvCount_Yz,MPI_INT,rank_Yz(),recvtag,Communicator,&req2[17]); MPI_Waitall(18,req1,stat1); MPI_Waitall(18,req2,stat2); //...................................................................................... @@ -513,7 +515,7 @@ void Domain::AssignComponentLabels(double *phase) vector Label; vector Affinity; // Read the labels - if (rank==0){ + if (rank()==0){ printf("Component labels:\n"); ifstream iFILE("ComponentLabels.csv"); if (iFILE.good()){ @@ -549,7 +551,7 @@ void Domain::AssignComponentLabels(double *phase) // Broadcast the list MPI_Bcast(&NLABELS,1,MPI_INT,0,Comm); - //printf("rank=%i, NLABELS=%i \n ",rank,NLABELS); + //printf("rank=%i, NLABELS=%i \n ",rank(),NLABELS); // Copy into contiguous buffers char *LabelList; @@ -557,11 +559,11 @@ void Domain::AssignComponentLabels(double *phase) LabelList=new char[NLABELS]; AffinityList=new double[NLABELS]; - if (rank==0){ + if (rank()==0){ for (int idx=0; idx < NLABELS; idx++){ VALUE=Label[idx]; AFFINITY=Affinity[idx]; - printf("rank=%i, idx=%i, value=%d, affinity=%f \n",rank,idx,VALUE,AFFINITY); + printf("rank=%i, idx=%i, value=%d, affinity=%f \n",rank(),idx,VALUE,AFFINITY); LabelList[idx]=VALUE; AffinityList[idx]=AFFINITY; } @@ -579,7 +581,7 @@ void Domain::AssignComponentLabels(double *phase) VALUE=id[n]; // Assign the affinity from the paired list for (int idx=0; idx < NLABELS; idx++){ - //printf("rank=%i, idx=%i, value=%i, %i, \n",rank,idx, VALUE,LabelList[idx]); + //printf("rank=%i, idx=%i, value=%i, %i, \n",rank(),idx, VALUE,LabelList[idx]); if (VALUE == LabelList[idx]){ AFFINITY=AffinityList[idx]; idx = NLABELS; @@ -617,42 +619,42 @@ void Domain::CommunicateMeshHalo(DoubleArray &Mesh) PackMeshData(sendList_yZ, sendCount_yZ ,sendData_yZ, MeshData); PackMeshData(sendList_YZ, sendCount_YZ ,sendData_YZ, MeshData); //...................................................................................... - MPI_Sendrecv(sendData_x,sendCount_x,MPI_DOUBLE,rank_x,sendtag, - recvData_X,recvCount_X,MPI_DOUBLE,rank_X,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_X,sendCount_X,MPI_DOUBLE,rank_X,sendtag, - recvData_x,recvCount_x,MPI_DOUBLE,rank_x,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_y,sendCount_y,MPI_DOUBLE,rank_y,sendtag, - recvData_Y,recvCount_Y,MPI_DOUBLE,rank_Y,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_Y,sendCount_Y,MPI_DOUBLE,rank_Y,sendtag, - recvData_y,recvCount_y,MPI_DOUBLE,rank_y,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_z,sendCount_z,MPI_DOUBLE,rank_z,sendtag, - recvData_Z,recvCount_Z,MPI_DOUBLE,rank_Z,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_Z,sendCount_Z,MPI_DOUBLE,rank_Z,sendtag, - recvData_z,recvCount_z,MPI_DOUBLE,rank_z,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_xy,sendCount_xy,MPI_DOUBLE,rank_xy,sendtag, - recvData_XY,recvCount_XY,MPI_DOUBLE,rank_XY,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_XY,sendCount_XY,MPI_DOUBLE,rank_XY,sendtag, - recvData_xy,recvCount_xy,MPI_DOUBLE,rank_xy,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_Xy,sendCount_Xy,MPI_DOUBLE,rank_Xy,sendtag, - recvData_xY,recvCount_xY,MPI_DOUBLE,rank_xY,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_xY,sendCount_xY,MPI_DOUBLE,rank_xY,sendtag, - recvData_Xy,recvCount_Xy,MPI_DOUBLE,rank_Xy,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_xz,sendCount_xz,MPI_DOUBLE,rank_xz,sendtag, - recvData_XZ,recvCount_XZ,MPI_DOUBLE,rank_XZ,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_XZ,sendCount_XZ,MPI_DOUBLE,rank_XZ,sendtag, - recvData_xz,recvCount_xz,MPI_DOUBLE,rank_xz,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_Xz,sendCount_Xz,MPI_DOUBLE,rank_Xz,sendtag, - recvData_xZ,recvCount_xZ,MPI_DOUBLE,rank_xZ,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_xZ,sendCount_xZ,MPI_DOUBLE,rank_xZ,sendtag, - recvData_Xz,recvCount_Xz,MPI_DOUBLE,rank_Xz,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_yz,sendCount_yz,MPI_DOUBLE,rank_yz,sendtag, - recvData_YZ,recvCount_YZ,MPI_DOUBLE,rank_YZ,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_YZ,sendCount_YZ,MPI_DOUBLE,rank_YZ,sendtag, - recvData_yz,recvCount_yz,MPI_DOUBLE,rank_yz,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_Yz,sendCount_Yz,MPI_DOUBLE,rank_Yz,sendtag, - recvData_yZ,recvCount_yZ,MPI_DOUBLE,rank_yZ,recvtag,Comm,MPI_STATUS_IGNORE); - MPI_Sendrecv(sendData_yZ,sendCount_yZ,MPI_DOUBLE,rank_yZ,sendtag, - recvData_Yz,recvCount_Yz,MPI_DOUBLE,rank_Yz,recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_x,sendCount_x,MPI_DOUBLE,rank_x(),sendtag, + recvData_X,recvCount_X,MPI_DOUBLE,rank_X(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_X,sendCount_X,MPI_DOUBLE,rank_X(),sendtag, + recvData_x,recvCount_x,MPI_DOUBLE,rank_x(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_y,sendCount_y,MPI_DOUBLE,rank_y(),sendtag, + recvData_Y,recvCount_Y,MPI_DOUBLE,rank_Y(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_Y,sendCount_Y,MPI_DOUBLE,rank_Y(),sendtag, + recvData_y,recvCount_y,MPI_DOUBLE,rank_y(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_z,sendCount_z,MPI_DOUBLE,rank_z(),sendtag, + recvData_Z,recvCount_Z,MPI_DOUBLE,rank_Z(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_Z,sendCount_Z,MPI_DOUBLE,rank_Z(),sendtag, + recvData_z,recvCount_z,MPI_DOUBLE,rank_z(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_xy,sendCount_xy,MPI_DOUBLE,rank_xy(),sendtag, + recvData_XY,recvCount_XY,MPI_DOUBLE,rank_XY(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_XY,sendCount_XY,MPI_DOUBLE,rank_XY(),sendtag, + recvData_xy,recvCount_xy,MPI_DOUBLE,rank_xy(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_Xy,sendCount_Xy,MPI_DOUBLE,rank_Xy(),sendtag, + recvData_xY,recvCount_xY,MPI_DOUBLE,rank_xY(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_xY,sendCount_xY,MPI_DOUBLE,rank_xY(),sendtag, + recvData_Xy,recvCount_Xy,MPI_DOUBLE,rank_Xy(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_xz,sendCount_xz,MPI_DOUBLE,rank_xz(),sendtag, + recvData_XZ,recvCount_XZ,MPI_DOUBLE,rank_XZ(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_XZ,sendCount_XZ,MPI_DOUBLE,rank_XZ(),sendtag, + recvData_xz,recvCount_xz,MPI_DOUBLE,rank_xz(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_Xz,sendCount_Xz,MPI_DOUBLE,rank_Xz(),sendtag, + recvData_xZ,recvCount_xZ,MPI_DOUBLE,rank_xZ(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_xZ,sendCount_xZ,MPI_DOUBLE,rank_xZ(),sendtag, + recvData_Xz,recvCount_Xz,MPI_DOUBLE,rank_Xz(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_yz,sendCount_yz,MPI_DOUBLE,rank_yz(),sendtag, + recvData_YZ,recvCount_YZ,MPI_DOUBLE,rank_YZ(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_YZ,sendCount_YZ,MPI_DOUBLE,rank_YZ(),sendtag, + recvData_yz,recvCount_yz,MPI_DOUBLE,rank_yz(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_Yz,sendCount_Yz,MPI_DOUBLE,rank_Yz(),sendtag, + recvData_yZ,recvCount_yZ,MPI_DOUBLE,rank_yZ(),recvtag,Comm,MPI_STATUS_IGNORE); + MPI_Sendrecv(sendData_yZ,sendCount_yZ,MPI_DOUBLE,rank_yZ(),sendtag, + recvData_Yz,recvCount_Yz,MPI_DOUBLE,rank_Yz(),recvtag,Comm,MPI_STATUS_IGNORE); //........................................................................................ UnpackMeshData(recvList_x, recvCount_x ,recvData_x, MeshData); UnpackMeshData(recvList_X, recvCount_X ,recvData_X, MeshData); diff --git a/common/Domain.h b/common/Domain.h index 7563b013..c2578855 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -64,7 +64,7 @@ private: class Domain{ public: //! Default constructor - Domain( std::shared_ptr db ); + Domain( std::shared_ptr db ); //! Obsolete constructor Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz, @@ -80,7 +80,7 @@ public: Domain& operator=( const Domain& ) = delete; //! Destructor - ~Domain(); + ~Domain(); //! Get the database inline std::shared_ptr getDatabase() const { return d_db; } @@ -108,91 +108,91 @@ private: public: // Public variables (need to create accessors instead) double Lx,Ly,Lz,Volume; - int Nx,Ny,Nz,N; - RankInfoStruct rank_info; + int Nx,Ny,Nz,N; + RankInfoStruct rank_info; - MPI_Comm Comm; // MPI Communicator for this domain + MPI_Comm Comm; // MPI Communicator for this domain - int BoundaryCondition; + int BoundaryCondition; - MPI_Group Group; // Group of processors associated with this domain + MPI_Group Group; // Group of processors associated with this domain - //********************************** - // MPI ranks for all 18 neighbors - //********************************** - const int& iproc = rank_info.ix; - const int& jproc = rank_info.jy; - const int& kproc = rank_info.kz; - const int& nprocx = rank_info.nx; - const int& nprocy = rank_info.ny; - const int& nprocz = rank_info.nz; - const int& rank = rank_info.rank[1][1][1]; - const int& rank_X = rank_info.rank[2][1][1]; - const int& rank_x = rank_info.rank[0][1][1]; - const int& rank_Y = rank_info.rank[1][2][1]; - const int& rank_y = rank_info.rank[1][0][1]; - const int& rank_Z = rank_info.rank[1][1][2]; - const int& rank_z = rank_info.rank[1][1][0]; - const int& rank_XY = rank_info.rank[2][2][1]; - const int& rank_xy = rank_info.rank[0][0][1]; - const int& rank_Xy = rank_info.rank[2][0][1]; - const int& rank_xY = rank_info.rank[0][2][1]; - const int& rank_XZ = rank_info.rank[2][1][2]; - const int& rank_xz = rank_info.rank[0][1][0]; - const int& rank_Xz = rank_info.rank[2][1][0]; - const int& rank_xZ = rank_info.rank[0][1][2]; - const int& rank_YZ = rank_info.rank[1][2][2]; - const int& rank_yz = rank_info.rank[1][0][0]; - const int& rank_Yz = rank_info.rank[1][2][0]; - const int& rank_yZ = rank_info.rank[1][0][2]; + //********************************** + // MPI ranks for all 18 neighbors + //********************************** + inline int iproc() const { return rank_info.ix; } + inline int jproc() const { return rank_info.jy; } + inline int kproc() const { return rank_info.kz; } + inline int nprocx() const { return rank_info.nx; } + inline int nprocy() const { return rank_info.ny; } + inline int nprocz() const { return rank_info.nz; } + inline int rank() const { return rank_info.rank[1][1][1]; } + inline int rank_X() const { return rank_info.rank[2][1][1]; } + inline int rank_x() const { return rank_info.rank[0][1][1]; } + inline int rank_Y() const { return rank_info.rank[1][2][1]; } + inline int rank_y() const { return rank_info.rank[1][0][1]; } + inline int rank_Z() const { return rank_info.rank[1][1][2]; } + inline int rank_z() const { return rank_info.rank[1][1][0]; } + inline int rank_XY() const { return rank_info.rank[2][2][1]; } + inline int rank_xy() const { return rank_info.rank[0][0][1]; } + inline int rank_Xy() const { return rank_info.rank[2][0][1]; } + inline int rank_xY() const { return rank_info.rank[0][2][1]; } + inline int rank_XZ() const { return rank_info.rank[2][1][2]; } + inline int rank_xz() const { return rank_info.rank[0][1][0]; } + inline int rank_Xz() const { return rank_info.rank[2][1][0]; } + inline int rank_xZ() const { return rank_info.rank[0][1][2]; } + inline int rank_YZ() const { return rank_info.rank[1][2][2]; } + inline int rank_yz() const { return rank_info.rank[1][0][0]; } + inline int rank_Yz() const { return rank_info.rank[1][2][0]; } + inline int rank_yZ() const { return rank_info.rank[1][0][2]; } - //********************************** - //...................................................................................... - // Get the actual D3Q19 communication counts (based on location of solid phase) - // Discrete velocity set symmetry implies the sendcount = recvcount - //...................................................................................... - int sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z; - int sendCount_xy, sendCount_yz, sendCount_xz, sendCount_Xy, sendCount_Yz, sendCount_xZ; - int sendCount_xY, sendCount_yZ, sendCount_Xz, sendCount_XY, sendCount_YZ, sendCount_XZ; - //...................................................................................... - int *sendList_x, *sendList_y, *sendList_z, *sendList_X, *sendList_Y, *sendList_Z; - int *sendList_xy, *sendList_yz, *sendList_xz, *sendList_Xy, *sendList_Yz, *sendList_xZ; - int *sendList_xY, *sendList_yZ, *sendList_Xz, *sendList_XY, *sendList_YZ, *sendList_XZ; - //...................................................................................... - int recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y, recvCount_Z; - int recvCount_xy, recvCount_yz, recvCount_xz, recvCount_Xy, recvCount_Yz, recvCount_xZ; - int recvCount_xY, recvCount_yZ, recvCount_Xz, recvCount_XY, recvCount_YZ, recvCount_XZ; - //...................................................................................... - int *recvList_x, *recvList_y, *recvList_z, *recvList_X, *recvList_Y, *recvList_Z; - int *recvList_xy, *recvList_yz, *recvList_xz, *recvList_Xy, *recvList_Yz, *recvList_xZ; - int *recvList_xY, *recvList_yZ, *recvList_Xz, *recvList_XY, *recvList_YZ, *recvList_XZ; - //...................................................................................... - // Solid indicator function - char *id; + //********************************** + //...................................................................................... + // Get the actual D3Q19 communication counts (based on location of solid phase) + // Discrete velocity set symmetry implies the sendcount = recvcount + //...................................................................................... + int sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z; + int sendCount_xy, sendCount_yz, sendCount_xz, sendCount_Xy, sendCount_Yz, sendCount_xZ; + int sendCount_xY, sendCount_yZ, sendCount_Xz, sendCount_XY, sendCount_YZ, sendCount_XZ; + //...................................................................................... + int *sendList_x, *sendList_y, *sendList_z, *sendList_X, *sendList_Y, *sendList_Z; + int *sendList_xy, *sendList_yz, *sendList_xz, *sendList_Xy, *sendList_Yz, *sendList_xZ; + int *sendList_xY, *sendList_yZ, *sendList_Xz, *sendList_XY, *sendList_YZ, *sendList_XZ; + //...................................................................................... + int recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y, recvCount_Z; + int recvCount_xy, recvCount_yz, recvCount_xz, recvCount_Xy, recvCount_Yz, recvCount_xZ; + int recvCount_xY, recvCount_yZ, recvCount_Xz, recvCount_XY, recvCount_YZ, recvCount_XZ; + //...................................................................................... + int *recvList_x, *recvList_y, *recvList_z, *recvList_X, *recvList_Y, *recvList_Z; + int *recvList_xy, *recvList_yz, *recvList_xz, *recvList_Xy, *recvList_Yz, *recvList_xZ; + int *recvList_xY, *recvList_yZ, *recvList_Xz, *recvList_XY, *recvList_YZ, *recvList_XZ; + //...................................................................................... + // Solid indicator function + char *id; - void CommunicateMeshHalo(DoubleArray &Mesh); - void AssignComponentLabels(double *phase); - void CommInit(MPI_Comm comm); - void TestCommInit(MPI_Comm comm); + void CommunicateMeshHalo(DoubleArray &Mesh); + void AssignComponentLabels(double *phase); + void CommInit(MPI_Comm comm); + void TestCommInit(MPI_Comm comm); //void MemoryOptimizedLayout(IntArray &Map, int *neighborList, int Np); private: - int *sendBuf_x, *sendBuf_y, *sendBuf_z, *sendBuf_X, *sendBuf_Y, *sendBuf_Z; - int *sendBuf_xy, *sendBuf_yz, *sendBuf_xz, *sendBuf_Xy, *sendBuf_Yz, *sendBuf_xZ; - int *sendBuf_xY, *sendBuf_yZ, *sendBuf_Xz, *sendBuf_XY, *sendBuf_YZ, *sendBuf_XZ; - //...................................................................................... - int *recvBuf_x, *recvBuf_y, *recvBuf_z, *recvBuf_X, *recvBuf_Y, *recvBuf_Z; - int *recvBuf_xy, *recvBuf_yz, *recvBuf_xz, *recvBuf_Xy, *recvBuf_Yz, *recvBuf_xZ; - int *recvBuf_xY, *recvBuf_yZ, *recvBuf_Xz, *recvBuf_XY, *recvBuf_YZ, *recvBuf_XZ; - //...................................................................................... - double *sendData_x, *sendData_y, *sendData_z, *sendData_X, *sendData_Y, *sendData_Z; - double *sendData_xy, *sendData_yz, *sendData_xz, *sendData_Xy, *sendData_Yz, *sendData_xZ; - double *sendData_xY, *sendData_yZ, *sendData_Xz, *sendData_XY, *sendData_YZ, *sendData_XZ; - double *recvData_x, *recvData_y, *recvData_z, *recvData_X, *recvData_Y, *recvData_Z; - double *recvData_xy, *recvData_yz, *recvData_xz, *recvData_Xy, *recvData_Yz, *recvData_xZ; - double *recvData_xY, *recvData_yZ, *recvData_Xz, *recvData_XY, *recvData_YZ, *recvData_XZ; + int *sendBuf_x, *sendBuf_y, *sendBuf_z, *sendBuf_X, *sendBuf_Y, *sendBuf_Z; + int *sendBuf_xy, *sendBuf_yz, *sendBuf_xz, *sendBuf_Xy, *sendBuf_Yz, *sendBuf_xZ; + int *sendBuf_xY, *sendBuf_yZ, *sendBuf_Xz, *sendBuf_XY, *sendBuf_YZ, *sendBuf_XZ; + //...................................................................................... + int *recvBuf_x, *recvBuf_y, *recvBuf_z, *recvBuf_X, *recvBuf_Y, *recvBuf_Z; + int *recvBuf_xy, *recvBuf_yz, *recvBuf_xz, *recvBuf_Xy, *recvBuf_Yz, *recvBuf_xZ; + int *recvBuf_xY, *recvBuf_yZ, *recvBuf_Xz, *recvBuf_XY, *recvBuf_YZ, *recvBuf_XZ; + //...................................................................................... + double *sendData_x, *sendData_y, *sendData_z, *sendData_X, *sendData_Y, *sendData_Z; + double *sendData_xy, *sendData_yz, *sendData_xz, *sendData_Xy, *sendData_Yz, *sendData_xZ; + double *sendData_xY, *sendData_yZ, *sendData_Xz, *sendData_XY, *sendData_YZ, *sendData_XZ; + double *recvData_x, *recvData_y, *recvData_z, *recvData_X, *recvData_Y, *recvData_Z; + double *recvData_xy, *recvData_yz, *recvData_xz, *recvData_Xy, *recvData_Yz, *recvData_xZ; + double *recvData_xY, *recvData_yZ, *recvData_Xz, *recvData_XY, *recvData_YZ, *recvData_XZ; }; diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index 6390ad2d..59785a79 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -15,25 +15,25 @@ ScaLBL_Communicator::ScaLBL_Communicator(Domain &Dm){ Nz = Dm.Nz; N = Nx*Ny*Nz; next=0; - rank=Dm.rank; - rank_x=Dm.rank_x; - rank_y=Dm.rank_y; - rank_z=Dm.rank_z; - rank_X=Dm.rank_X; - rank_Y=Dm.rank_Y; - rank_Z=Dm.rank_Z; - rank_xy=Dm.rank_xy; - rank_XY=Dm.rank_XY; - rank_xY=Dm.rank_xY; - rank_Xy=Dm.rank_Xy; - rank_xz=Dm.rank_xz; - rank_XZ=Dm.rank_XZ; - rank_xZ=Dm.rank_xZ; - rank_Xz=Dm.rank_Xz; - rank_yz=Dm.rank_yz; - rank_YZ=Dm.rank_YZ; - rank_yZ=Dm.rank_yZ; - rank_Yz=Dm.rank_Yz; + rank=Dm.rank(); + rank_x=Dm.rank_x(); + rank_y=Dm.rank_y(); + rank_z=Dm.rank_z(); + rank_X=Dm.rank_X(); + rank_Y=Dm.rank_Y(); + rank_Z=Dm.rank_Z(); + rank_xy=Dm.rank_xy(); + rank_XY=Dm.rank_XY(); + rank_xY=Dm.rank_xY(); + rank_Xy=Dm.rank_Xy(); + rank_xz=Dm.rank_xz(); + rank_XZ=Dm.rank_XZ(); + rank_xZ=Dm.rank_xZ(); + rank_Xz=Dm.rank_Xz(); + rank_yz=Dm.rank_yz(); + rank_YZ=Dm.rank_YZ(); + rank_yZ=Dm.rank_yZ(); + rank_Yz=Dm.rank_Yz(); sendCount_x=Dm.sendCount_x; sendCount_y=Dm.sendCount_y; sendCount_z=Dm.sendCount_z; @@ -71,12 +71,12 @@ ScaLBL_Communicator::ScaLBL_Communicator(Domain &Dm){ recvCount_YZ=Dm.recvCount_YZ; recvCount_XZ=Dm.recvCount_XZ; - iproc = Dm.iproc; - jproc = Dm.jproc; - kproc = Dm.kproc; - nprocx = Dm.nprocx; - nprocy = Dm.nprocy; - nprocz = Dm.nprocz; + iproc = Dm.iproc(); + jproc = Dm.jproc(); + kproc = Dm.kproc(); + nprocx = Dm.nprocx(); + nprocy = Dm.nprocy(); + nprocz = Dm.nprocz(); BoundaryCondition = Dm.BoundaryCondition; //...................................................................................... diff --git a/tests/TestColorSquareTube.cpp b/tests/TestColorSquareTube.cpp index c219e9b4..be72c0e7 100644 --- a/tests/TestColorSquareTube.cpp +++ b/tests/TestColorSquareTube.cpp @@ -277,12 +277,12 @@ int main(int argc, char **argv) if (rank==0) printf ("Initializing distributions \n"); // Initialize the phase field and variables ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm.last_interior, Np); - if (Dm.kproc==0){ + if (Dm.kproc()==0){ ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2); } - if (Dm.kproc == nprocz-1){ + if (Dm.kproc() == nprocz-1){ ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3); diff --git a/tests/TestPoiseuille.cpp b/tests/TestPoiseuille.cpp index 1abdaf4d..f1f0ad03 100644 --- a/tests/TestPoiseuille.cpp +++ b/tests/TestPoiseuille.cpp @@ -48,13 +48,10 @@ int main(int argc, char **argv) Fx = 0; Fy = 0; Fz = 1e-3; //1.f; // 1e-3; - auto FILENAME = argv[1]; - auto db = std::make_shared( FILENAME ); - auto domain_db = db->getDatabase( "Domain" ); // Load inputs if (rank==0) printf("Loading input database \n"); - auto db = std::make_shared(FILENAME); - auto domain_db= db-> getDatabase("Domain"); + auto db = std::make_shared(argv[1]); + auto domain_db = db->getDatabase("Domain"); int Nx = domain_db->getVector( "n" )[0]; int Ny = domain_db->getVector( "n" )[1]; int Nz = domain_db->getVector( "n" )[2]; diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index be7d00a4..18405907 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -12,6 +12,7 @@ #include "common/Array.h" #include "common/Domain.h" #include "analysis/eikonal.h" +#include "IO/Writer.h" std::shared_ptr loadInputs( int nprocs ) @@ -37,9 +38,6 @@ int main(int argc, char **argv) MPI_Comm_rank(comm,&rank); MPI_Comm_size(comm,&nprocs); { - int i,j,k,n,nn; - int iproc,jproc,kproc; - // Load inputs @@ -47,17 +45,14 @@ int main(int argc, char **argv) int Nx = db->getVector( "n" )[0]; int Ny = db->getVector( "n" )[1]; int Nz = db->getVector( "n" )[2]; - int nprocx = db->getVector( "nproc" )[0]; - int nprocy = db->getVector( "nproc" )[1]; - int nprocz = db->getVector( "nproc" )[2]; // Get the rank info Domain Dm(db); - for (k=0;k id(nx,ny,nz); + id.fill(0); - for (k=1;k ID0(id.size()); + ID0.copy( id ); + 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); + fillData.copy( ID0, ID ); + fillData.copy( Distance, dist1 ); + fillData.copy( TrueDist, 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 ) ); + IO::initialize( "", "silo", false ); + IO::writeData( "testSegDist", data, MPI_COMM_WORLD ); - MPI_Barrier(comm); } + MPI_Barrier(comm); MPI_Finalize(); return 0; diff --git a/tests/TestTorus.cpp b/tests/TestTorus.cpp index 38b6d45c..5c2fa519 100644 --- a/tests/TestTorus.cpp +++ b/tests/TestTorus.cpp @@ -103,9 +103,9 @@ int main(int argc, char **argv) n = k*Nx*Ny+j*Nx+i; // global position relative to center - x = Dm.iproc*Nx+i - CX; - y = Dm.jproc*Ny+j - CY; - z = Dm.kproc*Nz+k - CZ; + x = Dm.iproc()*Nx+i - CX; + y = Dm.jproc()*Ny+j - CY; + z = Dm.kproc()*Nz+k - CZ; // Shrink the sphere sizes by two voxels to make sure they don't touch Averages.SDs(i,j,k) = 100.0; @@ -113,12 +113,12 @@ int main(int argc, char **argv) // Single torus Averages.Phase(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; // Double torus - /* y = Dm.jproc*Ny+j - CY1; - //z = Dm.kproc*Nz+k - CZ +R1; + /* y = Dm.jproc()*Ny+j - CY1; + //z = Dm.kproc()*Nz+k - CZ +R1; Averages.Phase(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; - y = Dm.jproc*Ny+j - CY2; - //z = Dm.kproc*Nz+k - CZ-R1; + y = Dm.jproc()*Ny+j - CY2; + //z = Dm.kproc()*Nz+k - CZ-R1; Averages.Phase(i,j,k) = min(Averages.Phase(i,j,k), sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2); *///.............................................................................. diff --git a/tests/lbpm_captube_pp.cpp b/tests/lbpm_captube_pp.cpp index 348fabd5..ebe5480f 100644 --- a/tests/lbpm_captube_pp.cpp +++ b/tests/lbpm_captube_pp.cpp @@ -140,11 +140,11 @@ int main(int argc, char **argv) if (Averages.SDs(i,j,k) < 0.0){ id[n] = 0; } - else if (Dm.kproc*Nz+k 0 && Dm.kproc == 0){ + if (BoundaryCondition > 0 && Dm.kproc() == 0){ for (k=0; k<3; k++){ for (j=0;j 0 && Dm.kproc == nprocz-1){ + if (BoundaryCondition > 0 && Dm.kproc() == nprocz-1){ for (k=Nz-3; k0 ){ - if (Dm.kproc==0){ + if (Dm.kproc()==0){ ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2); } - if (Dm.kproc == nprocz-1){ + if (Dm.kproc() == nprocz-1){ ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3); diff --git a/tests/lbpm_dfh_simulator.cpp b/tests/lbpm_dfh_simulator.cpp index fbf5d554..d17cdad8 100644 --- a/tests/lbpm_dfh_simulator.cpp +++ b/tests/lbpm_dfh_simulator.cpp @@ -288,7 +288,7 @@ int main(int argc, char **argv) if (rank==0) printf("Media porosity = %f \n",porosity); //......................................................... // If external boundary conditions are applied remove solid - if (BoundaryCondition > 0 && Dm.kproc == 0){ + if (BoundaryCondition > 0 && Dm.kproc() == 0){ for (int k=0; k<3; k++){ for (int j=0;j 0 && Dm.kproc == nprocz-1){ + if (BoundaryCondition > 0 && Dm.kproc() == nprocz-1){ for (int k=Nz-3; k id(nx,ny,nz); TwoPhase Averages(Dm); // DoubleArray Distance(nx,ny,nz); // DoubleArray Phase(nx,ny,nz); @@ -235,8 +234,8 @@ int main(int argc, char **argv) for (i=0;i Date: Wed, 16 May 2018 14:07:55 -0400 Subject: [PATCH 2/3] uncertain changes --- tests/TestPoiseuille.cpp | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/tests/TestPoiseuille.cpp b/tests/TestPoiseuille.cpp index 1abdaf4d..6c3f7de9 100644 --- a/tests/TestPoiseuille.cpp +++ b/tests/TestPoiseuille.cpp @@ -47,14 +47,11 @@ int main(int argc, char **argv) double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); Fx = 0; Fy = 0; Fz = 1e-3; //1.f; // 1e-3; - + // Load inputs + if (rank==0) printf("Loading input database \n"); auto FILENAME = argv[1]; auto db = std::make_shared( FILENAME ); auto domain_db = db->getDatabase( "Domain" ); - // Load inputs - if (rank==0) printf("Loading input database \n"); - auto db = std::make_shared(FILENAME); - auto domain_db= db-> getDatabase("Domain"); int Nx = domain_db->getVector( "n" )[0]; int Ny = domain_db->getVector( "n" )[1]; int Nz = domain_db->getVector( "n" )[2]; From 1e3a77f8a939519ea8e49369cf29ce0e33589353 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 16 May 2018 14:37:21 -0400 Subject: [PATCH 3/3] refactor --- tests/GenerateSphereTest.cpp | 86 ++++++++++++++++-------------------- tests/TestBlobAnalyze.cpp | 6 +-- tests/TestBubbleDFH.cpp | 6 +-- tests/TestColorBubble.cpp | 8 ++-- tests/lbpm_random_pp.cpp | 2 +- tests/lbpm_refine_pp.cpp | 17 ++++--- 6 files changed, 57 insertions(+), 68 deletions(-) diff --git a/tests/GenerateSphereTest.cpp b/tests/GenerateSphereTest.cpp index 6c9edba8..3293205c 100644 --- a/tests/GenerateSphereTest.cpp +++ b/tests/GenerateSphereTest.cpp @@ -46,17 +46,14 @@ inline void UnpackID(int *list, int count, char *recvbuf, char *ID){ inline void MorphOpen(DoubleArray SignDist, char *id, Domain &Dm, int nx, int ny, int nz, int rank, double SW){ - int iproc = Dm.iproc; - int jproc = Dm.jproc; - int kproc = Dm.kproc; int i,j,k,n; - int nprocx=Dm.nprocx; - int nprocy=Dm.nprocy; - int nprocz=Dm.nprocz; double count,countGlobal,totalGlobal; count = 0.f; double maxdist=0.f; double maxdistGlobal; + int nprocx=Dm.nprocx(); + int nprocy=Dm.nprocy(); + int nprocz=Dm.nprocz(); for (int k=0; k