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

Cleaning up lbpm_segmented_pp, remove comments, re-indent,etc.

Browse Source
This commit is contained in:
James E McClure 2016-11-16 16:14:56 -05:00
parent 4fb098c8c0
commit 84d9fed6e7
1 changed files with 238 additions and 344 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

582
tests/lbpm_segmented_pp.cpp
View File

@ -40,140 +40,141 @@ inline double minmod(double &a, double &b){
inline double Eikonal(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ inline double Eikonal(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){
/* /*
* This routine converts the data in the Distance array to a signed distance * 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 * by solving the equation df/dt = sign(1-|grad f|), where Distance provides
* the values of f on the mesh associated with domain Dm * the values of f on the mesh associated with domain Dm
* It has been tested with segmented data initialized to values [-1,1] * 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 * and will converge toward the signed distance to the surface bounding the associated phases
* *
* Reference: * Reference:
* Min C (2010) On reinitializing level set functions, Journal of Computational Physics 229 * Min C (2010) On reinitializing level set functions, Journal of Computational Physics 229
*/ */
int i,j,k; int i,j,k;
double dt=0.1; double dt=0.1;
double Dx,Dy,Dz; double Dx,Dy,Dz;
double Dxp,Dxm,Dyp,Dym,Dzp,Dzm; double Dxp,Dxm,Dyp,Dym,Dzp,Dzm;
double Dxxp,Dxxm,Dyyp,Dyym,Dzzp,Dzzm; double Dxxp,Dxxm,Dyyp,Dyym,Dzzp,Dzzm;
double sign,norm; double sign,norm;
double LocalVar,GlobalVar,LocalMax,GlobalMax; double LocalVar,GlobalVar,LocalMax,GlobalMax;
int xdim,ydim,zdim; int xdim,ydim,zdim;
xdim=Dm.Nx-2; xdim=Dm.Nx-2;
ydim=Dm.Ny-2; ydim=Dm.Ny-2;
zdim=Dm.Nz-2; zdim=Dm.Nz-2;
fillHalo<double> fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1); fillHalo<double> fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1);
// Arrays to store the second derivatives // Arrays to store the second derivatives
DoubleArray Dxx(Dm.Nx,Dm.Ny,Dm.Nz); DoubleArray Dxx(Dm.Nx,Dm.Ny,Dm.Nz);
DoubleArray Dyy(Dm.Nx,Dm.Ny,Dm.Nz); DoubleArray Dyy(Dm.Nx,Dm.Ny,Dm.Nz);
DoubleArray Dzz(Dm.Nx,Dm.Ny,Dm.Nz); DoubleArray Dzz(Dm.Nx,Dm.Ny,Dm.Nz);
int count = 0; int count = 0;
while (count < timesteps){ while (count < timesteps){
// Communicate the halo of values // Communicate the halo of values
fillData.fill(Distance); fillData.fill(Distance);
// Compute second order derivatives // Compute second order derivatives
for (k=1;k<Dm.Nz-1;k++){ for (k=1;k<Dm.Nz-1;k++){
for (j=1;j<Dm.Ny-1;j++){ for (j=1;j<Dm.Ny-1;j++){
for (i=1;i<Dm.Nx-1;i++){ 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); 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); 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); Dzz(i,j,k) = Distance(i,j,k+1) + Distance(i,j,k-1) - 2*Distance(i,j,k);
} }
} }
} }
fillData.fill(Dxx); fillData.fill(Dxx);
fillData.fill(Dyy); fillData.fill(Dyy);
fillData.fill(Dzz); fillData.fill(Dzz);
LocalMax=LocalVar=0.0; LocalMax=LocalVar=0.0;
// Execute the next timestep // Execute the next timestep
for (k=1;k<Dm.Nz-1;k++){ for (k=1;k<Dm.Nz-1;k++){
for (j=1;j<Dm.Ny-1;j++){ for (j=1;j<Dm.Ny-1;j++){
for (i=1;i<Dm.Nx-1;i++){ for (i=1;i<Dm.Nx-1;i++){
int n = k*Dm.Nx*Dm.Ny + j*Dm.Nx + i; int n = k*Dm.Nx*Dm.Ny + j*Dm.Nx + i;
sign = -1; sign = -1;
if (ID[n] == 1) sign = 1; if (ID[n] == 1) sign = 1;
// local second derivative terms // local second derivative terms
Dxxp = minmod(Dxx(i,j,k),Dxx(i+1,j,k)); Dxxp = minmod(Dxx(i,j,k),Dxx(i+1,j,k));
Dyyp = minmod(Dyy(i,j,k),Dyy(i,j+1,k)); Dyyp = minmod(Dyy(i,j,k),Dyy(i,j+1,k));
Dzzp = minmod(Dzz(i,j,k),Dzz(i,j,k+1)); Dzzp = minmod(Dzz(i,j,k),Dzz(i,j,k+1));
Dxxm = minmod(Dxx(i,j,k),Dxx(i-1,j,k)); Dxxm = minmod(Dxx(i,j,k),Dxx(i-1,j,k));
Dyym = minmod(Dyy(i,j,k),Dyy(i,j-1,k)); Dyym = minmod(Dyy(i,j,k),Dyy(i,j-1,k));
Dzzm = minmod(Dzz(i,j,k),Dzz(i,j,k-1)); Dzzm = minmod(Dzz(i,j,k),Dzz(i,j,k-1));
/* //............Compute upwind derivatives ................... /* //............Compute upwind derivatives ...................
Dxp = Distance(i+1,j,k) - Distance(i,j,k) + 0.5*Dxxp; 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; 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; 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; 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; 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; Dzm = Distance(i,j,k) - Distance(i,j,k-1) + 0.5*Dzzm;
*/ */
Dxp = Distance(i+1,j,k); Dxp = Distance(i+1,j,k);
Dyp = Distance(i,j+1,k); Dyp = Distance(i,j+1,k);
Dzp = Distance(i,j,k+1); Dzp = Distance(i,j,k+1);
Dxm = Distance(i-1,j,k); Dxm = Distance(i-1,j,k);
Dym = Distance(i,j-1,k); Dym = Distance(i,j-1,k);
Dzm = Distance(i,j,k-1); Dzm = Distance(i,j,k-1);
// Compute upwind derivatives for Godunov Hamiltonian // Compute upwind derivatives for Godunov Hamiltonian
if (sign < 0.0){ if (sign < 0.0){
if (Dxp > Dxm) Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp; if (Dxp > Dxm) Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp;
else Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm; else Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm;
if (Dyp > Dym) Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp; if (Dyp > Dym) Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp;
else Dy = Distance(i,j,k) - Dym + 0.5*Dyym; else Dy = Distance(i,j,k) - Dym + 0.5*Dyym;
if (Dzp > Dzm) Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp; if (Dzp > Dzm) Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp;
else Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm; else Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm;
} }
else{ else{
if (Dxp < Dxm) Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp; if (Dxp < Dxm) Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp;
else Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm; else Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm;
if (Dyp < Dym) Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp; if (Dyp < Dym) Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp;
else Dy = Distance(i,j,k) - Dym + 0.5*Dyym; else Dy = Distance(i,j,k) - Dym + 0.5*Dyym;
if (Dzp < Dzm) Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp; if (Dzp < Dzm) Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp;
else Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm; else Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm;
} }
norm=sqrt(Dx*Dx+Dy*Dy+Dz*Dz); norm=sqrt(Dx*Dx+Dy*Dy+Dz*Dz);
if (norm > 1.0) norm=1.0; 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) Distance(i,j,k) += dt*sign*(1.0 - norm);
LocalMax = fabs(dt*sign*(1.0 - norm)); LocalVar += dt*sign*(1.0 - norm);
}
}
}
MPI_Allreduce(&LocalVar,&GlobalVar,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); if (fabs(dt*sign*(1.0 - norm)) > LocalMax)
MPI_Allreduce(&LocalMax,&GlobalMax,1,MPI_DOUBLE,MPI_MAX,Dm.Comm); LocalMax = fabs(dt*sign*(1.0 - norm));
GlobalVar /= (Dm.Nx-2)*(Dm.Ny-2)*(Dm.Nz-2)*Dm.nprocx*Dm.nprocy*Dm.nprocz; }
count++; }
}
if (count%50 == 0 && Dm.rank==0 ) MPI_Allreduce(&LocalVar,&GlobalVar,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
printf("Time=%i, Max variation=%f, Global variation=%f \n",count,GlobalMax,GlobalVar); 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 (fabs(GlobalMax) < 1e-5){ if (count%50 == 0 && Dm.rank==0 )
if (Dm.rank==0) printf("Exiting with max tolerance of 1e-5 \n"); printf("Time=%i, Max variation=%f, Global variation=%f \n",count,GlobalMax,GlobalVar);
count=timesteps;
if (fabs(GlobalMax) < 1e-5){
if (Dm.rank==0) printf("Exiting with max tolerance of 1e-5 \n");
count=timesteps;
}
} }
} return GlobalVar;
return GlobalVar;
} }
@ -182,270 +183,163 @@ int main(int argc, char **argv)
// Initialize MPI // Initialize MPI
int rank, nprocs; int rank, nprocs;
MPI_Init(&argc,&argv); MPI_Init(&argc,&argv);
MPI_Comm comm = MPI_COMM_WORLD; MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank); MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs); MPI_Comm_size(comm,&nprocs);
{ {
//....................................................................... //.......................................................................
// Reading the domain information file // Reading the domain information file
//....................................................................... //.......................................................................
int nprocx, nprocy, nprocz, nx, ny, nz, nspheres; int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
double Lx, Ly, Lz; double Lx, Ly, Lz;
int Nx,Ny,Nz; int Nx,Ny,Nz;
int i,j,k,n; int i,j,k,n;
int BC=0; int BC=0;
char Filename[40]; char Filename[40];
int xStart,yStart,zStart; int xStart,yStart,zStart;
// char fluidValue,solidValue; // char fluidValue,solidValue;
std::vector<char> solidValues; std::vector<char> solidValues;
std::vector<char> nwpValues; std::vector<char> nwpValues;
std::string line; std::string line;
if (rank==0){ if (rank==0){
ifstream domain("Domain.in"); ifstream domain("Domain.in");
domain >> nprocx; domain >> nprocx;
domain >> nprocy; domain >> nprocy;
domain >> nprocz; domain >> nprocz;
domain >> nx; domain >> nx;
domain >> ny; domain >> ny;
domain >> nz; domain >> nz;
domain >> nspheres; domain >> nspheres;
domain >> Lx; domain >> Lx;
domain >> Ly; domain >> Ly;
domain >> Lz; domain >> Lz;
ifstream image("Segmented.in"); ifstream image("Segmented.in");
image >> Filename; // Name of data file containing segmented data image >> Filename; // Name of data file containing segmented data
image >> Nx; // size of the binary file image >> Nx; // size of the binary file
image >> Ny; image >> Ny;
image >> Nz; image >> Nz;
image >> xStart; // offset for the starting voxel image >> xStart; // offset for the starting voxel
image >> yStart; image >> yStart;
image >> zStart; image >> zStart;
} }
MPI_Barrier(comm); MPI_Barrier(comm);
// Computational domain // Computational domain
MPI_Bcast(&nx,1,MPI_INT,0,comm); MPI_Bcast(&nx,1,MPI_INT,0,comm);
MPI_Bcast(&ny,1,MPI_INT,0,comm); MPI_Bcast(&ny,1,MPI_INT,0,comm);
MPI_Bcast(&nz,1,MPI_INT,0,comm); MPI_Bcast(&nz,1,MPI_INT,0,comm);
MPI_Bcast(&nprocx,1,MPI_INT,0,comm); MPI_Bcast(&nprocx,1,MPI_INT,0,comm);
MPI_Bcast(&nprocy,1,MPI_INT,0,comm); MPI_Bcast(&nprocy,1,MPI_INT,0,comm);
MPI_Bcast(&nprocz,1,MPI_INT,0,comm); MPI_Bcast(&nprocz,1,MPI_INT,0,comm);
MPI_Bcast(&nspheres,1,MPI_INT,0,comm); MPI_Bcast(&nspheres,1,MPI_INT,0,comm);
MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm);
//................................................. //.................................................
MPI_Barrier(comm); MPI_Barrier(comm);
// Check that the number of processors >= the number of ranks // Check that the number of processors >= the number of ranks
if ( rank==0 ) { if ( rank==0 ) {
printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz); printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
printf("Number of MPI ranks used: %i \n", nprocs); printf("Number of MPI ranks used: %i \n", nprocs);
printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz); printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz);
} }
if ( nprocs < nprocx*nprocy*nprocz ){ if ( nprocs < nprocx*nprocy*nprocz ){
ERROR("Insufficient number of processors"); ERROR("Insufficient number of processors");
} }
char LocalRankFilename[40]; char LocalRankFilename[40];
int N = (nx+2)*(ny+2)*(nz+2); int N = (nx+2)*(ny+2)*(nz+2);
Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC); Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
for (n=0; n<N; n++) Dm.id[n]=1; for (n=0; n<N; n++) Dm.id[n]=1;
Dm.CommInit(comm); Dm.CommInit(comm);
// Read the phase ID // Read the phase ID
size_t readID; size_t readID;
sprintf(LocalRankFilename,"ID.%05i",rank); sprintf(LocalRankFilename,"ID.%05i",rank);
FILE *ID = fopen(LocalRankFilename,"rb"); FILE *ID = fopen(LocalRankFilename,"rb");
readID=fread(Dm.id,1,N,ID); readID=fread(Dm.id,1,N,ID);
if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID \n"); if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID \n");
fclose(ID); fclose(ID);
// make sure communication // make sure communication
// Set up layers in x direction // Set up layers in x direction
for (k=0; k<nz; k++){ for (k=0; k<nz; k++){
for (j=0; j<ny; j++){ for (j=0; j<ny; j++){
Dm.id[k*nx*ny+j*nx]=1; Dm.id[k*nx*ny+j*nx]=1;
Dm.id[k*nx*ny+j*nx+nx-1] = 1; Dm.id[k*nx*ny+j*nx+nx-1] = 1;
}
}
for (k=0; k<nz; k++){
for (i=0; i<nx; i++){
Dm.id[k*nx*ny+i]=1;
Dm.id[k*nx*ny+(ny-1)*nx+i] = 1;
}
}
for (j=0; j<ny; j++){
for (i=0; i<nx; i++){
Dm.id[j*nx+i]=1;
Dm.id[nx*ny*(nz-1)+j*nx+i] = 1;
}
}
// Initialize the domain and communication
nx+=2; ny+=2; nz+=2;
int count = 0;
N=nx*ny*nz;
char *id;
id = new char [N];
TwoPhase Averages(Dm);
// DoubleArray Distance(nx,ny,nz);
// DoubleArray Phase(nx,ny,nz);
// Solve for the position of the solid phase
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;
// Initialize the solid phase
if (Dm.id[n] == 0) id[n] = 0;
else id[n] = 1;
} }
} }
}
// Initialize the signed distance function for (k=0; k<nz; k++){
for (k=0;k<nz;k++){ for (i=0; i<nx; i++){
for (j=0;j<ny;j++){ Dm.id[k*nx*ny+i]=1;
for (i=0;i<nx;i++){ Dm.id[k*nx*ny+(ny-1)*nx+i] = 1;
n=k*nx*ny+j*nx+i;
// Initialize distance to +/- 1
Averages.SDs(i,j,k) = 2.0*id[n]-1.0;
} }
} }
}
MeanFilter(Averages.SDs);
double LocalVar, TotalVar; for (j=0; j<ny; j++){
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n"); for (i=0; i<nx; i++){
int Maxtime=10*max(max(Dm.Nx*Dm.nprocx,Dm.Ny*Dm.nprocy),Dm.Nz*Dm.nprocz); Dm.id[j*nx+i]=1;
LocalVar = Eikonal(Averages.SDs,id,Dm,Maxtime); Dm.id[nx*ny*(nz-1)+j*nx+i] = 1;
MPI_Allreduce(&LocalVar,&TotalVar,1,MPI_DOUBLE,MPI_SUM,comm);
TotalVar /= nprocs;
if (rank==0) printf("Final variation in signed distance function %f \n",TotalVar);
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");
fwrite(Averages.SDs.data(),8,Averages.SDs.length(),DIST);
fclose(DIST);
/* // Solve for the position of the non-wetting phase
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;
// Initialize the non-wetting phase
if (Dm.id[n] == 1) id[n] = 1;
else id[n] = 0;
} }
} }
}
// 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;
// Initialize distance to +/- 1
Averages.Phase(i,j,k) = 2.0*id[n]-1.0;
}
}
}
MeanFilter(Averages.Phase);
if (rank==0) printf("Initialized non-wetting phase -- Converting to Signed Distance function \n"); // Initialize the domain and communication
SSO(Averages.Phase,id,Dm,100);
for (k=0;k<nz;k++){ nx+=2; ny+=2; nz+=2;
for (j=0;j<ny;j++){ int count = 0;
for (i=0;i<nx;i++){ N=nx*ny*nz;
n=k*nx*ny+j*nx+i;
Averages.Phase(i,j,k) -= 1.0; char *id;
// Initialize distance to +/- 1 id = new char [N];
// Dilation of the non-wetting phase TwoPhase Averages(Dm);
Averages.SDn(i,j,k) = -Averages.Phase(i,j,k); // DoubleArray Distance(nx,ny,nz);
Averages.Phase(i,j,k) = Averages.SDn(i,j,k); // DoubleArray Phase(nx,ny,nz);
Averages.Phase_tplus(i,j,k) = Averages.SDn(i,j,k);
Averages.Phase_tminus(i,j,k) = Averages.SDn(i,j,k); // Solve for the position of the solid phase
Averages.DelPhi(i,j,k) = 0.0; for (k=0;k<nz;k++){
Averages.Press(i,j,k) = 0.0; for (j=0;j<ny;j++){
Averages.Vel_x(i,j,k) = 0.0; for (i=0;i<nx;i++){
Averages.Vel_y(i,j,k) = 0.0; n = k*nx*ny+j*nx+i;
Averages.Vel_z(i,j,k) = 0.0; // Initialize the solid phase
if (Averages.SDs(i,j,k) > 0.0){ if (Dm.id[n] == 0) id[n] = 0;
if (Averages.Phase(i,j,k) > 0.0){ else id[n] = 1;
Dm.id[n] = 2;
}
else{
Dm.id[n] = 1;
}
}
else{
Dm.id[n] = 0;
} }
} }
} }
// 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;
// Initialize distance to +/- 1
Averages.SDs(i,j,k) = 2.0*id[n]-1.0;
}
}
}
MeanFilter(Averages.SDs);
double LocalVar, TotalVar;
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
int Maxtime=10*max(max(Dm.Nx*Dm.nprocx,Dm.Ny*Dm.nprocy),Dm.Nz*Dm.nprocz);
LocalVar = Eikonal(Averages.SDs,id,Dm,Maxtime);
MPI_Allreduce(&LocalVar,&TotalVar,1,MPI_DOUBLE,MPI_SUM,comm);
TotalVar /= nprocs;
if (rank==0) printf("Final variation in signed distance function %f \n",TotalVar);
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");
fwrite(Averages.SDs.data(),8,Averages.SDs.length(),DIST);
fclose(DIST);
} }
MPI_Barrier(comm);
// Create the MeshDataStruct
fillHalo<double> fillData(Dm.Comm,Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
std::vector<IO::MeshDataStruct> meshData(1);
meshData[0].meshName = "domain";
meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) );
std::shared_ptr<IO::Variable> PhaseVar( new IO::Variable() );
std::shared_ptr<IO::Variable> SolidVar( new IO::Variable() );
std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );
PhaseVar->name = "Fluid";
PhaseVar->type = IO::VolumeVariable;
PhaseVar->dim = 1;
PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(PhaseVar);
SolidVar->name = "Solid";
SolidVar->type = IO::VolumeVariable;
SolidVar->dim = 1;
SolidVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(SignDistVar);
BlobIDVar->name = "BlobID";
BlobIDVar->type = IO::VolumeVariable;
BlobIDVar->dim = 1;
BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(BlobIDVar);
fillData.copy(Averages.SDn,PhaseVar->data);
fillData.copy(Averages.SDs,SolidVar->data);
fillData.copy(Averages.Label_NWP,BlobIDVar->data);
IO::writeData( 0, meshData, 2, comm );
// sprintf(LocalRankFilename,"Phase.%05i",rank);
// FILE *PHASE = fopen(LocalRankFilename,"wb");
// fwrite(Averages.Phase.get(),8,Averages.Phase.length(),PHASE);
// fclose(PHASE);
double beta = 0.95;
if (rank==0) printf("initializing the system \n");
Averages.UpdateSolid();
Averages.UpdateMeshValues();
Dm.CommunicateMeshHalo(Averages.Phase);
Dm.CommunicateMeshHalo(Averages.SDn);
Dm.CommunicateMeshHalo(Averages.SDs);
int timestep=5;
Averages.Initialize();
if (rank==0) printf("computing phase components \n");
Averages.ComponentAverages();
if (rank==0) printf("sorting phase components \n");
Averages.SortBlobs();
Averages.PrintComponents(timestep);
*/
}
MPI_Barrier(comm);
MPI_Finalize(); MPI_Finalize();
return 0; return 0;
} }