OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'dfh' of github.com:JamesEMcClure/LBPM-WIA into dfh

Browse Source
This commit is contained in:
James E McClure 2018-05-16 14:41:30 -04:00
commit 9c7eee30f7
28 changed files with 646 additions and 782 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

92
analysis/TwoPhase.cpp
View File

@ -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<kmax; k++){
for (j=1; j<Ny-1; j++){
@ -585,15 +585,15 @@ void TwoPhase::ComponentAverages()
ComponentAverages_WP.fill(0.0);
ComponentAverages_NWP.fill(0.0);
if (Dm.rank==0){
if (Dm.rank()==0){
printf("Number of wetting phase components is %i \n",NumberComponents_WP);
printf("Number of non-wetting phase components is %i \n",NumberComponents_NWP);
}
// 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<kmax; k++){
for (j=1; j<Ny-1; j++){
@ -636,9 +636,9 @@ void TwoPhase::ComponentAverages()
ComponentAverages_NWP(VY,LabelNWP) += 0.125*Vel_y(n);
ComponentAverages_NWP(VZ,LabelNWP) += 0.125*Vel_z(n);
// center of mass
ComponentAverages_NWP(CMX,LabelNWP) += 0.125*(i+cube[p][0]+Dm.iproc*Nx);
ComponentAverages_NWP(CMY,LabelNWP) += 0.125*(j+cube[p][1]+Dm.jproc*Ny);
ComponentAverages_NWP(CMZ,LabelNWP) += 0.125*(k+cube[p][2]+Dm.kproc*Nz);
ComponentAverages_NWP(CMX,LabelNWP) += 0.125*(i+cube[p][0]+Dm.iproc()*Nx);
ComponentAverages_NWP(CMY,LabelNWP) += 0.125*(j+cube[p][1]+Dm.jproc()*Ny);
ComponentAverages_NWP(CMZ,LabelNWP) += 0.125*(k+cube[p][2]+Dm.kproc()*Nz);
// twice the kinetic energy
ComponentAverages_NWP(VSQ,LabelNWP) += 0.125*(Vel_x(n)*Vel_x(n)+Vel_y(n)*Vel_y(n)+Vel_z(n)*Vel_z(n));
@ -656,9 +656,9 @@ void TwoPhase::ComponentAverages()
ComponentAverages_WP(VY,LabelWP)+= 0.125*Vel_y(n);
ComponentAverages_WP(VZ,LabelWP) += 0.125*Vel_z(n);
// Center of mass
ComponentAverages_WP(CMX,LabelWP) += 0.125*(i+cube[p][0]+Dm.iproc*Nx);
ComponentAverages_WP(CMY,LabelWP) += 0.125*(j+cube[p][1]+Dm.jproc*Ny);
ComponentAverages_WP(CMZ,LabelWP) += 0.125*(k+cube[p][2]+Dm.kproc*Nz);
ComponentAverages_WP(CMX,LabelWP) += 0.125*(i+cube[p][0]+Dm.iproc()*Nx);
ComponentAverages_WP(CMY,LabelWP) += 0.125*(j+cube[p][1]+Dm.jproc()*Ny);
ComponentAverages_WP(CMZ,LabelWP) += 0.125*(k+cube[p][2]+Dm.kproc()*Nz);
// twice the kinetic energy
ComponentAverages_WP(VSQ,LabelWP) += 0.125*(Vel_x(n)*Vel_x(n)+Vel_y(n)*Vel_y(n)+Vel_z(n)*Vel_z(n));
@ -803,7 +803,7 @@ void TwoPhase::ComponentAverages()
}
MPI_Barrier(Dm.Comm);
if (Dm.rank==0){
if (Dm.rank()==0){
printf("Component averages computed locally -- reducing result... \n");
}
// Globally reduce the non-wetting phase averages
@ -819,7 +819,7 @@ void TwoPhase::ComponentAverages()
MPI_Allreduce(ComponentAverages_NWP.data(),RecvBuffer.data(),BLOB_AVG_COUNT*NumberComponents_NWP, MPI_DOUBLE,MPI_SUM,Dm.Comm);
// MPI_Reduce(ComponentAverages_NWP.data(),RecvBuffer.data(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,0,Dm.Comm);
if (Dm.rank==0){
if (Dm.rank()==0){
printf("rescaling... \n");
}
@ -907,7 +907,7 @@ void TwoPhase::ComponentAverages()
}
}
if (Dm.rank==0){
if (Dm.rank()==0){
printf("reduce WP averages... \n");
}
@ -1052,15 +1052,15 @@ void TwoPhase::WriteSurfaces(int logcount)
C = P;
}
// 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);
wn_mesh->A.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<NumberComponents_NWP; b++){
//if (ComponentAverages_NWP(TRIMVOL,b) > 0.0){

336
analysis/eikonal.cpp
View File

@ -27,199 +27,55 @@ static inline double minmod(double &a, double &b){
}
template<class TYPE>
static inline MPI_Datatype getType( );
template<> inline MPI_Datatype getType<float>() { return MPI_FLOAT; }
template<> inline MPI_Datatype getType<double>() { return MPI_DOUBLE; }
/******************************************************************
* Solve the eikonal equation *
******************************************************************/
double Eikonal(DoubleArray &Distance, const char *ID, const Domain &Dm, int timesteps)
template<class TYPE>
TYPE Eikonal( Array<TYPE> &Distance, const Array<char> &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<double> 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<Dm.Nz-1;k++){
for (j=1;j<Dm.Ny-1;j++){
for (i=1;i<Dm.Nx-1;i++){
Dxx(i,j,k) = Distance(i+1,j,k) + Distance(i-1,j,k) - 2*Distance(i,j,k);
Dyy(i,j,k) = Distance(i,j+1,k) + Distance(i,j-1,k) - 2*Distance(i,j,k);
Dzz(i,j,k) = Distance(i,j,k+1) + Distance(i,j,k-1) - 2*Distance(i,j,k);
}
}
}
fillData.fill(Dxx);
fillData.fill(Dyy);
fillData.fill(Dzz);
LocalMax=LocalVar=0.0;
// Execute the next timestep
for (k=1;k<Dm.Nz-1;k++){
for (j=1;j<Dm.Ny-1;j++){
for (i=1;i<Dm.Nx-1;i++){
int n = k*Dm.Nx*Dm.Ny + j*Dm.Nx + i;
sign = 1;
if (ID[n] == 0) sign = -1;
// local second derivative terms
Dxxp = minmod(Dxx(i,j,k),Dxx(i+1,j,k));
Dyyp = minmod(Dyy(i,j,k),Dyy(i,j+1,k));
Dzzp = minmod(Dzz(i,j,k),Dzz(i,j,k+1));
Dxxm = minmod(Dxx(i,j,k),Dxx(i-1,j,k));
Dyym = minmod(Dyy(i,j,k),Dyy(i,j-1,k));
Dzzm = minmod(Dzz(i,j,k),Dzz(i,j,k-1));
/* //............Compute upwind derivatives ...................
Dxp = Distance(i+1,j,k) - Distance(i,j,k) + 0.5*Dxxp;
Dyp = Distance(i,j+1,k) - Distance(i,j,k) + 0.5*Dyyp;
Dzp = Distance(i,j,k+1) - Distance(i,j,k) + 0.5*Dzzp;
Dxm = Distance(i,j,k) - Distance(i-1,j,k) + 0.5*Dxxm;
Dym = Distance(i,j,k) - Distance(i,j-1,k) + 0.5*Dyym;
Dzm = Distance(i,j,k) - Distance(i,j,k-1) + 0.5*Dzzm;
*/
Dxp = Distance(i+1,j,k)- Distance(i,j,k) - 0.5*Dxxp;
Dyp = Distance(i,j+1,k)- Distance(i,j,k) - 0.5*Dyyp;
Dzp = Distance(i,j,k+1)- Distance(i,j,k) - 0.5*Dzzp;
Dxm = Distance(i,j,k) - Distance(i-1,j,k) + 0.5*Dxxm;
Dym = Distance(i,j,k) - Distance(i,j-1,k) + 0.5*Dyym;
Dzm = Distance(i,j,k) - Distance(i,j,k-1) + 0.5*Dzzm;
// Compute upwind derivatives for Godunov Hamiltonian
if (sign < 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 < 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<float> &Distance, const Array<char> &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<float> 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<TYPE> fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1);
// Arrays to store the second derivatives
Array<float> Dxx(Dm.Nx,Dm.Ny,Dm.Nz);
Array<float> Dyy(Dm.Nx,Dm.Ny,Dm.Nz);
Array<float> Dzz(Dm.Nx,Dm.Ny,Dm.Nz);
Array<TYPE> Dxx(Dm.Nx,Dm.Ny,Dm.Nz);
Array<TYPE> Dyy(Dm.Nx,Dm.Ny,Dm.Nz);
Array<TYPE> Dzz(Dm.Nx,Dm.Ny,Dm.Nz);
Array<int8_t> sign(ID.size());
for (size_t i=0; i<sign.length(); i++)
sign(i) = ID(i) == 0 ? -1:1;
int count = 0;
double dt = 0.1;
double GlobalVar = 0.0;
while (count < timesteps){
// Communicate the halo of values
fillData.fill(Distance);
// Compute second order derivatives
for (k=1;k<Dm.Nz-1;k++){
for (j=1;j<Dm.Ny-1;j++){
for (i=1;i<Dm.Nx-1;i++){
for (int k=1;k<Dm.Nz-1;k++){
for (int j=1;j<Dm.Ny-1;j++){
for (int i=1;i<Dm.Nx-1;i++){
Dxx(i,j,k) = Distance(i+1,j,k) + Distance(i-1,j,k) - 2*Distance(i,j,k);
Dyy(i,j,k) = Distance(i,j+1,k) + Distance(i,j-1,k) - 2*Distance(i,j,k);
Dzz(i,j,k) = Distance(i,j,k+1) + Distance(i,j,k-1) - 2*Distance(i,j,k);
@ -230,100 +86,112 @@ float Eikonal3D( Array<float> &Distance, const Array<char> &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<Dm.Nz-1;k++){
for (j=1;j<Dm.Ny-1;j++){
for (i=1;i<Dm.Nx-1;i++){
int n = k*Dm.Nx*Dm.Ny + j*Dm.Nx + i;
sign = -1;
if (ID(i,j,k) == 1) sign = 1;
for (int k=1;k<Dm.Nz-1;k++){
for (int j=1;j<Dm.Ny-1;j++){
for (int i=1;i<Dm.Nx-1;i++){
double s = sign(i,j,k);
// local second derivative terms
Dxxp = minmod(Dxx(i,j,k),Dxx(i+1,j,k));
Dyyp = minmod(Dyy(i,j,k),Dyy(i,j+1,k));
Dzzp = minmod(Dzz(i,j,k),Dzz(i,j,k+1));
Dxxm = minmod(Dxx(i,j,k),Dxx(i-1,j,k));
Dyym = minmod(Dyy(i,j,k),Dyy(i,j-1,k));
Dzzm = minmod(Dzz(i,j,k),Dzz(i,j,k-1));
double Dxxp = minmod(Dxx(i,j,k),Dxx(i+1,j,k));
double Dyyp = minmod(Dyy(i,j,k),Dyy(i,j+1,k));
double Dzzp = minmod(Dzz(i,j,k),Dzz(i,j,k+1));
double Dxxm = minmod(Dxx(i,j,k),Dxx(i-1,j,k));
double Dyym = minmod(Dyy(i,j,k),Dyy(i,j-1,k));
double Dzzm = minmod(Dzz(i,j,k),Dzz(i,j,k-1));
/* //............Compute upwind derivatives ...................
Dxp = Distance(i+1,j,k) - Distance(i,j,k) + 0.5*Dxxp;
Dyp = Distance(i,j+1,k) - Distance(i,j,k) + 0.5*Dyyp;
Dzp = Distance(i,j,k+1) - Distance(i,j,k) + 0.5*Dzzp;
Dxm = Distance(i,j,k) - Distance(i-1,j,k) + 0.5*Dxxm;
Dym = Distance(i,j,k) - Distance(i,j-1,k) + 0.5*Dyym;
Dzm = Distance(i,j,k) - Distance(i,j,k-1) + 0.5*Dzzm;
*/
Dxp = Distance(i+1,j,k);
Dyp = Distance(i,j+1,k);
Dzp = Distance(i,j,k+1);
Dxm = Distance(i-1,j,k);
Dym = Distance(i,j-1,k);
Dzm = Distance(i,j,k-1);
*/
double Dxp = Distance(i+1,j,k)- Distance(i,j,k) - 0.5*Dxxp;
double Dyp = Distance(i,j+1,k)- Distance(i,j,k) - 0.5*Dyyp;
double Dzp = Distance(i,j,k+1)- Distance(i,j,k) - 0.5*Dzzp;
double Dxm = Distance(i,j,k) - Distance(i-1,j,k) + 0.5*Dxxm;
double Dym = Distance(i,j,k) - Distance(i,j-1,k) + 0.5*Dyym;
double Dzm = Distance(i,j,k) - Distance(i,j,k-1) + 0.5*Dzzm;
// Compute upwind derivatives for Godunov Hamiltonian
if (sign < 0.0){
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;
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<float>( Array<float>&, const Array<char>&, const Domain&, int );
template double Eikonal<double>( Array<double>&, const Array<char>&, const Domain&, int );
/******************************************************************
* A fast distance calculation *
******************************************************************/
bool CalcDist3DIteration( Array<float> &Distance, const Domain &Dm )
bool CalcDist3DIteration( Array<float> &Distance, const Domain & )
{
const float sq2 = sqrt(2.0f);
const float sq3 = sqrt(3.0f);

7
analysis/eikonal.h
View File

@ -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<float> &Distance, const Array<char> &ID, const Domain &Dm, const int timesteps);
template<class TYPE>
TYPE Eikonal( Array<TYPE> &Distance, const Array<char> &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

2
analysis/runAnalysis.cpp
View File

@ -294,7 +294,7 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> 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;

2
analysis/uCT.cpp
View File

@ -149,7 +149,7 @@ void solve( const Array<float>& VOL, Array<float>& Mean, Array<char>& 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

240
common/Domain.cpp
View File

@ -28,7 +28,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),
@ -65,8 +65,8 @@ Domain::Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz,
initialize( db );
}
Domain::Domain( std::shared_ptr<Database> 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),
@ -121,13 +121,13 @@ void Domain::initialize( std::shared_ptr<Database> 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<int>("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()
{
@ -194,7 +194,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);
@ -314,42 +316,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);
@ -374,42 +376,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);
//......................................................................................
@ -508,7 +510,7 @@ void Domain::AssignComponentLabels(double *phase)
vector <char> Label;
vector <double> Affinity;
// Read the labels
if (rank==0){
if (rank()==0){
printf("Component labels:\n");
ifstream iFILE("ComponentLabels.csv");
if (iFILE.good()){
@ -544,7 +546,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;
@ -552,11 +554,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;
}
@ -574,7 +576,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;
@ -612,42 +614,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);

152
common/Domain.h
View File

@ -64,7 +64,7 @@ private:
class Domain{
public:
//! Default constructor
Domain( std::shared_ptr<Database> db );
Domain( std::shared_ptr<Database> 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<const Database> 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;
};

50
common/ScaLBL.cpp
View File

@ -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;
//......................................................................................

86
tests/GenerateSphereTest.cpp
View File

@ -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<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
@ -136,7 +133,6 @@ inline void MorphOpen(DoubleArray SignDist, char *id, Domain &Dm, int nx, int ny
int sendtag,recvtag;
sendtag = recvtag = 7;
int x,y,z;
int ii,jj,kk;
int Nx = nx;
int Ny = ny;
@ -218,42 +214,42 @@ inline void MorphOpen(DoubleArray SignDist, char *id, Domain &Dm, int nx, int ny
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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.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,Dm.Comm,MPI_STATUS_IGNORE);
//......................................................................................
UnpackID(Dm.recvList_x, Dm.recvCount_x ,recvID_x, id);
UnpackID(Dm.recvList_X, Dm.recvCount_X ,recvID_X, id);
@ -331,7 +327,6 @@ int main(int argc, char **argv)
// parallel domain size (# of sub-domains)
int nprocx,nprocy,nprocz;
int iproc,jproc,kproc;
int sendtag,recvtag;
//*****************************************
// MPI ranks for all 18 neighbors
//**********************************
@ -340,8 +335,6 @@ int main(int argc, char **argv)
int rank_xz,rank_XZ,rank_xZ,rank_Xz;
int rank_yz,rank_YZ,rank_yZ,rank_Yz;
//**********************************
MPI_Request req1[18],req2[18];
MPI_Status stat1[18],stat2[18];
if (rank == 0){
printf("********************************************************\n");
@ -351,7 +344,6 @@ int main(int argc, char **argv)
// Variables that specify the computational domain
string FILENAME;
unsigned int nBlocks, nthreads;
int Nx,Ny,Nz; // local sub-domain size
int nspheres; // number of spheres in the packing
double Lx,Ly,Lz; // Domain length

6
tests/TestBlobAnalyze.cpp
View File

@ -219,7 +219,7 @@ int main(int argc, char **argv)
Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
// const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
TwoPhase Averages(Dm);
int N = (nx+2)*(ny+2)*(nz+2);
Nx = nx+2;
Ny = ny+2;
Nz = nz+2;
@ -264,7 +264,7 @@ int main(int argc, char **argv)
MPI_Barrier(comm);
//.......................................................................
SignedDistance(Averages.Phase.data(),nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,
Dm.iproc,Dm.jproc,Dm.kproc,Dm.nprocx,Dm.nprocy,Dm.nprocz);
Dm.iproc(),Dm.jproc(),Dm.kproc(),Dm.nprocx(),Dm.nprocy(),Dm.nprocz());
//.......................................................................
// Assign the phase ID field based on the signed distance
//.......................................................................
@ -294,8 +294,6 @@ int main(int argc, char **argv)
}
}
}
double beta = 0.95;
if (rank==0) printf("initializing the system \n");
Averages.UpdateSolid();

6
tests/TestBubbleDFH.cpp
View File

@ -317,9 +317,9 @@ int main(int argc, char **argv)
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nz + i;
int iglobal= i+(Nx-2)*Dm.iproc;
int jglobal= j+(Ny-2)*Dm.jproc;
int kglobal= k+(Nz-2)*Dm.kproc;
int iglobal= i+(Nx-2)*Dm.iproc();
int jglobal= j+(Ny-2)*Dm.jproc();
int kglobal= k+(Nz-2)*Dm.kproc();
// Initialize phase position field for parallel bubble test
if ((iglobal-0.5*(Nx-2)*nprocx)*(iglobal-0.5*(Nx-2)*nprocx)
+(jglobal-0.5*(Ny-2)*nprocy)*(jglobal-0.5*(Ny-2)*nprocy)

8
tests/TestColorBubble.cpp
View File

@ -257,7 +257,7 @@ int main(int argc, char **argv)
}
}
}
printf("rank=%i, %i,%i,%i \n",rank,Dm.iproc,Dm.jproc,Dm.jproc);
printf("rank=%i, %i,%i,%i \n",rank,Dm.iproc(),Dm.jproc(),Dm.jproc());
// Initialize a bubble
int BubbleRadius=Nx/3;
int center_x = (Nx-2)*nprocx/2;
@ -268,9 +268,9 @@ int main(int argc, char **argv)
for (j=1;j<Ny-1;j++){
for (i=1;i<Nx-1;i++){
n = k*Nx*Ny + j*Nx + i;
int iglobal= i+(Nx-2)*Dm.iproc;
int jglobal= j+(Ny-2)*Dm.jproc;
int kglobal= k+(Nz-2)*Dm.kproc;
int iglobal= i+(Nx-2)*Dm.iproc();
int jglobal= j+(Ny-2)*Dm.jproc();
int kglobal= k+(Nz-2)*Dm.kproc();
// Initialize phase position field for parallel bubble test
if ((iglobal-center_x)*(iglobal-center_x)

4
tests/TestColorSquareTube.cpp
View File

@ -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);

7
tests/TestPoiseuille.cpp
View File

@ -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<Database>( FILENAME );
auto domain_db = db->getDatabase( "Domain" );
// Load inputs
if (rank==0) printf("Loading input database \n");
auto db = std::make_shared<Database>(FILENAME);
auto domain_db= db-> getDatabase("Domain");
int Nx = domain_db->getVector<int>( "n" )[0];
int Ny = domain_db->getVector<int>( "n" )[1];
int Nz = domain_db->getVector<int>( "n" )[2];

106
tests/TestSegDist.cpp
View File

@ -12,6 +12,7 @@
#include "common/Array.h"
#include "common/Domain.h"
#include "analysis/eikonal.h"
#include "IO/Writer.h"
std::shared_ptr<Database> 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<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "n" )[2];
int nprocx = db->getVector<int>( "nproc" )[0];
int nprocy = db->getVector<int>( "nproc" )[1];
int nprocz = db->getVector<int>( "nproc" )[2];
// Get the rank info
Domain Dm(db);
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
Dm.id[n] = 1;
}
}
@ -67,40 +62,34 @@ int main(int argc, char **argv)
int nx = Nx+2;
int ny = Ny+2;
int nz = Nz+2;
int N = nx*ny*nz;
int count = 0;
char *id;
id = new char [N];
double BubbleRadius = 25;
// Initialize the bubble
double x,y,z , Cx,Cy,Cz;
Cx = 1.0*nx;
Cy = 1.0*ny;
Cz = 1.0*nz;
// Initialize the bubble
double Cx = 1.0*nx;
double Cy = 1.0*ny;
double Cz = 1.0*nz;
DoubleArray Distance(nx,ny,nz);
DoubleArray TrueDist(nx,ny,nz);
Array<char> id(nx,ny,nz);
id.fill(0);
for (k=1;k<nz-1;k++){
for (j=1;j<ny-1;j++){
for (i=1;i<nx-1;i++){
for (int k=1;k<nz-1;k++){
for (int j=1;j<ny-1;j++){
for (int i=1;i<nx-1;i++){
// True signed distance
x = (nx-2)*Dm.iproc+i-1;
y = (ny-2)*Dm.jproc+j-1;
z = (nz-2)*Dm.kproc+k-1;
double x = (nx-2)*Dm.iproc()+i-1;
double y = (ny-2)*Dm.jproc()+j-1;
double z = (nz-2)*Dm.kproc()+k-1;
TrueDist(i,j,k) = sqrt((x-Cx)*(x-Cx)+(y-Cy)*(y-Cy)+(z-Cz)*(z-Cz)) - BubbleRadius;
n = k*nx*ny+j*nx+i;
// Initialize phase positions
if (TrueDist(i,j,k) < 0.0){
id[n] = 0;
}
else{
id[n]=1;
id(i,j,k) = 0;
} else{
id(i,j,k)=1;
}
@ -109,12 +98,11 @@ int main(int argc, char **argv)
}
// Initialize the signed distance function
for (k=0;k<nz;k++){
for (j=0;j<ny;j++){
for (i=0;i<nx;i++){
n=k*nx*ny+j*nx+i;
for (int k=0;k<nz;k++){
for (int j=0;j<ny;j++){
for (int i=0;i<nx;i++){
// Initialize distance to +/- 1
Distance(i,j,k) = 1.0*id[n]-0.5;
Distance(i,j,k) = 2.0*id(i,j,k)-1.0;
}
}
}
@ -128,31 +116,55 @@ int main(int argc, char **argv)
double starttime,stoptime,cputime;
starttime = MPI_Wtime();
Eikonal(Distance,id,Dm,10);
double err1 = Eikonal(Distance,id,Dm,1000);
stoptime = MPI_Wtime();
cputime = (stoptime - starttime);
if (rank==0) printf("Total time: %f seconds \n",cputime);
double Error=0.0;
int Count = 0;
for (k=0;k<nz;k++){
for (j=0;j<ny;j++){
for (i=0;i<nx;i++){
double localError=0.0;
int localCount = 0;
for (int k=0;k<nz;k++){
for (int j=0;j<ny;j++){
for (int i=0;i<nx;i++){
if (fabs(TrueDist(i,j,k)) < 3.0){
Error += (Distance(i,j,k)-TrueDist(i,j,k))*(Distance(i,j,k)-TrueDist(i,j,k));
Count += 1;
localError += (Distance(i,j,k)-TrueDist(i,j,k))*(Distance(i,j,k)-TrueDist(i,j,k));
localCount++;
}
}
}
}
Error = sqrt(Error)/(double (Count));
if (rank==0) printf("Mean error %f \n", Error);
double globalError;
int globalCount;
MPI_Allreduce(&localError,&globalError,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&localCount,&globalCount,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
double err2 = sqrt(globalError)/(double (globalCount));
if (rank==0) printf("global variation %f \n", err1);
if (rank==0) printf("Mean error %f \n", err2);
// Write the results to visit
Array<int> ID0(id.size());
ID0.copy( id );
Array<double> ID(Nx,Ny,Nz);
Array<double> dist1(Nx,Ny,Nz);
Array<double> dist2(Nx,Ny,Nz);
fillHalo<double> 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<IO::MeshDataStruct> 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;

14
tests/TestTorus.cpp
View File

@ -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);
*///..............................................................................

4
tests/lbpm_captube_pp.cpp
View File

@ -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<BubbleBottom){
else if (Dm.kproc()*Nz+k<BubbleBottom){
id[n] = 2;
sum++;
}
else if (Dm.kproc*Nz+k<BubbleTop){
else if (Dm.kproc()*Nz+k<BubbleTop){
id[n] = 1;
sum++;
}

8
tests/lbpm_color_macro_simulator.cpp
View File

@ -367,7 +367,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 (k=0; k<3; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
@ -378,7 +378,7 @@ int main(int argc, char **argv)
}
}
}
if (BoundaryCondition > 0 && Dm.kproc == nprocz-1){
if (BoundaryCondition > 0 && Dm.kproc() == nprocz-1){
for (k=Nz-3; k<Nz; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
@ -488,12 +488,12 @@ int main(int argc, char **argv)
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np);
if (BoundaryCondition >0 ){
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);

4
tests/lbpm_dfh_simulator.cpp
View File

@ -286,7 +286,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<Ny;j++){
for (int i=0;i<Nx;i++){
@ -297,7 +297,7 @@ int main(int argc, char **argv)
}
}
}
if (BoundaryCondition > 0 && Dm.kproc == nprocz-1){
if (BoundaryCondition > 0 && Dm.kproc() == nprocz-1){
for (int k=Nz-3; k<Nz; k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){

4
tests/lbpm_inkbottle_pp.cpp
View File

@ -177,11 +177,11 @@ int main(int argc, char **argv)
if (Averages.SDs(i,j,k) < 0.0){
id[n] = 0;
}
else if (Dm.kproc*Nz+k<BubbleBottom){
else if (Dm.kproc()*Nz+k<BubbleBottom){
id[n] = 2;
sum++;
}
else if (Dm.kproc*Nz+k<BubbleTop){
else if (Dm.kproc()*Nz+k<BubbleTop){
id[n] = 1;
sum++;
}

6
tests/lbpm_morph_pp.cpp
View File

@ -168,9 +168,9 @@ int main(int argc, char **argv)
Dm.CommInit(comm);
int iproc = Dm.iproc;
int jproc = Dm.jproc;
int kproc = Dm.kproc;
int iproc = Dm.iproc();
int jproc = Dm.jproc();
int kproc = Dm.kproc();
// Generate the NWP configuration
//if (rank==0) printf("Performing morphological drainage with critical radius %f \n", Rcrit);

78
tests/lbpm_morphdrain_pp.cpp
View File

@ -178,9 +178,9 @@ int main(int argc, char **argv)
Dm.CommInit(comm);
int iproc = Dm.iproc;
int jproc = Dm.jproc;
int kproc = Dm.kproc;
int iproc = Dm.iproc();
int jproc = Dm.jproc();
int kproc = Dm.kproc();
// Generate the NWP configuration
//if (rank==0) printf("Performing morphological drainage with critical radius %f \n", Rcrit);
@ -360,42 +360,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);

78
tests/lbpm_morphopen_pp.cpp
View File

@ -271,9 +271,9 @@ int main(int argc, char **argv)
*/
Dm.CommInit(comm);
int iproc = Dm.iproc;
int jproc = Dm.jproc;
int kproc = Dm.kproc;
int iproc = Dm.iproc();
int jproc = Dm.jproc();
int kproc = Dm.kproc();
// Generate the NWP configuration
//if (rank==0) printf("Initializing morphological distribution with critical radius %f \n", Rcrit);
@ -415,42 +415,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);

78
tests/lbpm_random_pp.cpp
View File

@ -172,9 +172,9 @@ int main(int argc, char **argv)
if (rank==0) printf("Media Porosity: %f \n",porosity);
Dm.CommInit(comm);
int iproc = Dm.iproc;
int jproc = Dm.jproc;
int kproc = Dm.kproc;
int iproc = Dm.iproc();
int jproc = Dm.jproc();
int kproc = Dm.kproc();
int bin, binCount;
ifstream Dist("BlobSize.in");
@ -346,42 +346,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);

17
tests/lbpm_refine_pp.cpp
View File

@ -170,7 +170,7 @@ int main(int argc, char **argv)
// Write output blocks with the same sub-domain size as origina
// refinement increases the size of the process grid
writerank = 8*Dm.kproc*nprocx*nprocy + 4*Dm.jproc*nprocx + 2*Dm.iproc;
writerank = 8*Dm.kproc()*nprocx*nprocy + 4*Dm.jproc()*nprocx + 2*Dm.iproc();
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
@ -198,7 +198,7 @@ int main(int argc, char **argv)
fwrite(id,1,nx*ny*nz,WRITEID);
fclose(WRITEID);
writerank = 8*Dm.kproc*nprocx*nprocy + 4*Dm.jproc*nprocx + 2*Dm.iproc+1;
writerank = 8*Dm.kproc()*nprocx*nprocy + 4*Dm.jproc()*nprocx + 2*Dm.iproc()+1;
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
@ -227,7 +227,7 @@ int main(int argc, char **argv)
fclose(WRITEID);
writerank = (2*Dm.kproc)*4*nprocx*nprocy + (2*Dm.jproc+1)*2*nprocx + 2*Dm.iproc+1;
writerank = (2*Dm.kproc())*4*nprocx*nprocy + (2*Dm.jproc()+1)*2*nprocx + 2*Dm.iproc()+1;
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
@ -255,7 +255,7 @@ int main(int argc, char **argv)
fwrite(id,1,nx*ny*nz,WRITEID);
fclose(WRITEID);
writerank = (2*Dm.kproc)*4*nprocx*nprocy + (2*Dm.jproc+1)*2*nprocx + 2*Dm.iproc;
writerank = (2*Dm.kproc())*4*nprocx*nprocy + (2*Dm.jproc()+1)*2*nprocx + 2*Dm.iproc();
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
@ -283,8 +283,7 @@ int main(int argc, char **argv)
fwrite(id,1,nx*ny*nz,WRITEID);
fclose(WRITEID);
writerank = (2*Dm.kproc+1)*4*nprocx*nprocy + (2*Dm.jproc)*2*nprocx + 2*Dm.iproc;
writerank = (2*Dm.kproc()+1)*4*nprocx*nprocy + (2*Dm.jproc())*2*nprocx + 2*Dm.iproc();
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
@ -312,7 +311,7 @@ int main(int argc, char **argv)
fwrite(id,1,nx*ny*nz,WRITEID);
fclose(WRITEID);
writerank = (2*Dm.kproc+1)*4*nprocx*nprocy + (2*Dm.jproc)*2*nprocx + 2*Dm.iproc+1;
writerank = (2*Dm.kproc()+1)*4*nprocx*nprocy + (2*Dm.jproc())*2*nprocx + 2*Dm.iproc()+1;
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
@ -340,7 +339,7 @@ int main(int argc, char **argv)
fwrite(id,1,nx*ny*nz,WRITEID);
fclose(WRITEID);
writerank = (2*Dm.kproc+1)*4*nprocx*nprocy + (2*Dm.jproc+1)*2*nprocx + 2*Dm.iproc;
writerank = (2*Dm.kproc()+1)*4*nprocx*nprocy + (2*Dm.jproc()+1)*2*nprocx + 2*Dm.iproc();
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
@ -368,7 +367,7 @@ int main(int argc, char **argv)
fwrite(id,1,nx*ny*nz,WRITEID);
fclose(WRITEID);
writerank = (2*Dm.kproc+1)*4*nprocx*nprocy + (2*Dm.jproc+1)*2*nprocx + 2*Dm.iproc+1;
writerank = (2*Dm.kproc()+1)*4*nprocx*nprocy + (2*Dm.jproc()+1)*2*nprocx + 2*Dm.iproc()+1;
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){

8
tests/lbpm_segmented_decomp.cpp
View File

@ -151,7 +151,7 @@ int main(int argc, char **argv)
// Set up the sub-domains
if (rank==0){
printf("Distributing subdomains across %i processors \n",nprocs);
printf("Process grid: %i x %i x %i \n",Dm.nprocx,Dm.nprocy,Dm.nprocz);
printf("Process grid: %i x %i x %i \n",Dm.nprocx(),Dm.nprocy(),Dm.nprocz());
printf("Subdomain size: %i \n",N);
printf("Size of transition region: %i \n", z_transition_size);
char *tmp;
@ -160,7 +160,7 @@ int main(int argc, char **argv)
for (int jp=0; jp<nprocy; jp++){
for (int ip=0; ip<nprocx; ip++){
// rank of the process that gets this subdomain
int rnk = kp*Dm.nprocx*Dm.nprocy + jp*Dm.nprocx + ip;
int rnk = kp*Dm.nprocx()*Dm.nprocy() + jp*Dm.nprocx() + ip;
// Pack and send the subdomain for rnk
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
@ -342,8 +342,8 @@ int main(int argc, char **argv)
if (rank==0) printf("Writing symmetric domain reflection\n");
MPI_Barrier(comm);
int symrank,sympz;
sympz = 2*nprocz - Dm.kproc -1;
symrank = sympz*nprocx*nprocy + Dm.jproc*nprocx + Dm.iproc;
sympz = 2*nprocz - Dm.kproc() -1;
symrank = sympz*nprocx*nprocy + Dm.jproc()*nprocx + Dm.iproc();
// DoubleArray SymDist(nx,ny,nz);
char *symid;

13
tests/lbpm_segmented_pp.cpp
View File

@ -208,7 +208,7 @@ int main(int argc, char **argv)
// Read from the large block and write the local ID file
if (rank==0) printf("Reading ID file from blocks \n");
fflush(stdout);
porosity = ReadFromBlock(Dm.id,Dm.iproc,Dm.jproc,Dm.kproc,nx,ny,nz);
porosity = ReadFromBlock(Dm.id,Dm.iproc(),Dm.jproc(),Dm.kproc(),nx,ny,nz);
MPI_Barrier(MPI_COMM_WORLD);
if (rank==0) printf("Writing local ID files (poros=%f) \n",porosity);
@ -222,8 +222,7 @@ int main(int argc, char **argv)
}
// Initialize the domain and communication
char *id;
id = new char [N];
Array<char> 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<nx;i++){
n = k*nx*ny+j*nx+i;
// Initialize the solid phase
if (Dm.id[n] == 0) id[n] = 0;
else id[n] = 1;
if (Dm.id[n] == 0) id(i,j,k) = 0;
else id(i,j,k) = 1;
}
}
}
@ -246,7 +245,7 @@ int main(int argc, char **argv)
for (i=0;i<nx;i++){
n=k*nx*ny+j*nx+i;
// Initialize distance to +/- 1
Averages.SDs(i,j,k) = 2.0*double(id[n])-1.0;
Averages.SDs(i,j,k) = 2.0*double(id(i,j,k))-1.0;
}
}
}
@ -254,7 +253,7 @@ int main(int argc, char **argv)
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);
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);

12
tests/lbpm_squaretube_pp.cpp
View File

@ -169,11 +169,11 @@ int main(int argc, char **argv)
if (Averages.SDs(i,j,k) < 0.0){
id[n] = 0;
}
else if (Dm.iproc*Nx+k<BubbleBottom){
else if (Dm.iproc()*Nx+k<BubbleBottom){
id[n] = 2;
sum++;
}
else if (Dm.iproc*Nx+k<BubbleTop){
else if (Dm.iproc()*Nx+k<BubbleTop){
id[n] = 1;
sum++;
}
@ -190,11 +190,11 @@ int main(int argc, char **argv)
if (Averages.SDs(i,j,k) < 0.0){
id[n] = 0;
}
else if (Dm.jproc*Ny+k<BubbleBottom){
else if (Dm.jproc()*Ny+k<BubbleBottom){
id[n] = 2;
sum++;
}
else if (Dm.jproc*Ny+k<BubbleTop){
else if (Dm.jproc()*Ny+k<BubbleTop){
id[n] = 1;
sum++;
}
@ -211,11 +211,11 @@ int main(int argc, char **argv)
if (Averages.SDs(i,j,k) < 0.0){
id[n] = 0;
}
else if (Dm.kproc*Nz+k<BubbleBottom){
else if (Dm.kproc()*Nz+k<BubbleBottom){
id[n] = 2;
sum++;
}
else if (Dm.kproc*Nz+k<BubbleTop){
else if (Dm.kproc()*Nz+k<BubbleTop){
id[n] = 1;
sum++;
}