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

Merge branch 'master' of https://github.com/JamesEMcClure/LBPM-WIA

Browse Source
This commit is contained in:
James E McClure 2017-04-18 19:14:02 -04:00
commit ac797709ce
4 changed files with 432 additions and 420 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
common/Domain.cpp
View File

@ -21,6 +21,7 @@ using namespace std;
// Reading the domain information file
void read_domain( int rank, int nprocs, MPI_Comm comm,
int& nprocx, int& nprocy, int& nprocz, int& nx, int& ny, int& nz,

155
common/Domain.h
View File

@ -24,7 +24,6 @@ void read_domain( int rank, int nprocs, MPI_Comm comm,
int& nprocx, int& nprocy, int& nprocz, int& nx, int& ny, int& nz,
int& nspheres, double& Lx, double& Ly, double& Lz );
//! Class to hold domain info
struct Domain{
// Default constructor
@ -597,6 +596,160 @@ inline void ReadBinaryFile(char *FILENAME, double *Data, int N)
File.close();
}
inline double minmod(double &a, double &b){
double value;
value = a;
if ( a*b < 0.0) value=0.0;
else if (fabs(a) > fabs(b)) value = b;
return value;
}
inline double Eikonal(DoubleArray &Distance, char *ID, Domain &Dm, 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] == 1) 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);
if (fabs(GlobalMax) < 1e-5){
if (Dm.rank==0) printf("Exiting with max tolerance of 1e-5 \n");
count=timesteps;
}
}
return GlobalVar;
}
#endif

540
tests/lbpm_segmented_decomp.cpp
View File

@ -18,294 +18,308 @@ int main(int argc, char **argv)
// Initialize MPI
int rank, nprocs;
MPI_Init(&argc,&argv);
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
{
int SOLID=atoi(argv[1]);
int NWP=atoi(argv[2]);
//char NWP,SOLID;
//SOLID=argv[1][0];
//NWP=argv[2][0];
if (rank==0){
printf("Solid Label: %i \n",SOLID);
printf("NWP Label: %i \n",NWP);
}
//.......................................................................
// Reading the domain information file
//.......................................................................
int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
double Lx, Ly, Lz;
int Nx,Ny,Nz;
int i,j,k,n;
int BC=0;
char Filename[40];
int xStart,yStart,zStart;
// char fluidValue,solidValue;
std::vector<char> solidValues;
std::vector<char> nwpValues;
std::string line;
if (rank==0){
ifstream domain("Domain.in");
domain >> nprocx;
domain >> nprocy;
domain >> nprocz;
domain >> nx;
domain >> ny;
domain >> nz;
domain >> nspheres;
domain >> Lx;
domain >> Ly;
domain >> Lz;
bool MULTINPUT=false;
ifstream image("Segmented.in");
image >> Filename; // Name of data file containing segmented data
image >> Nx; // size of the binary file
image >> Ny;
image >> Nz;
image >> xStart; // offset for the starting voxel
image >> yStart;
image >> zStart;
int NWP,SOLID,rank_offset;
SOLID=atoi(argv[1]);
NWP=atoi(argv[2]);
//char NWP,SOLID;
//SOLID=argv[1][0];
//NWP=argv[2][0];
if (rank==0){
printf("Solid Label: %i \n",SOLID);
printf("NWP Label: %i \n",NWP);
}
if (argc > 3){
rank_offset = atoi(argv[3]);
}
else{
MULTINPUT=true;
rank_offset=0;
}
}
MPI_Barrier(comm);
// Computational domain
//.................................................
MPI_Bcast(&nx,1,MPI_INT,0,comm);
MPI_Bcast(&ny,1,MPI_INT,0,comm);
MPI_Bcast(&nz,1,MPI_INT,0,comm);
MPI_Bcast(&nprocx,1,MPI_INT,0,comm);
MPI_Bcast(&nprocy,1,MPI_INT,0,comm);
MPI_Bcast(&nprocz,1,MPI_INT,0,comm);
MPI_Bcast(&nspheres,1,MPI_INT,0,comm);
MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Lz,1,MPI_DOUBLE,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(&xStart,1,MPI_INT,0,comm);
MPI_Bcast(&yStart,1,MPI_INT,0,comm);
MPI_Bcast(&zStart,1,MPI_INT,0,comm);
//.................................................
MPI_Barrier(comm);
//.......................................................................
// Reading the domain information file
//.......................................................................
int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
double Lx, Ly, Lz;
int Nx,Ny,Nz;
int i,j,k,n;
int BC=0;
char Filename[40];
int xStart,yStart,zStart;
// char fluidValue,solidValue;
// Check that the number of processors >= the number of ranks
if ( rank==0 ) {
printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
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);
}
if ( nprocs < nprocx*nprocy*nprocz ){
ERROR("Insufficient number of processors");
}
char *SegData = NULL;
// Rank=0 reads the entire segmented data and distributes to worker processes
if (rank==0){
printf("Dimensions of segmented image: %i x %i x %i \n",Nx,Ny,Nz);
SegData = new char[Nx*Ny*Nz];
FILE *SEGDAT = fopen(Filename,"rb");
if (SEGDAT==NULL) ERROR("Error reading segmented data");
size_t ReadSeg;
ReadSeg=fread(SegData,1,Nx*Ny*Nz,SEGDAT);
if (ReadSeg != size_t(Nx*Ny*Nz)) printf("lbpm_segmented_decomp: Error reading segmented data (rank=%i)\n",rank);
fclose(SEGDAT);
printf("Read segmented data from %s \n",Filename);
}
MPI_Barrier(comm);
std::vector<char> solidValues;
std::vector<char> nwpValues;
std::string line;
// Get the rank info
int N = (nx+2)*(ny+2)*(nz+2);
Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
n = k*(nx+2)*(ny+2)+j*(nx+2)+i;
Dm.id[n] = 1;
if (rank==0){
ifstream domain("Domain.in");
domain >> nprocx;
domain >> nprocy;
domain >> nprocz;
domain >> nx;
domain >> ny;
domain >> nz;
domain >> nspheres;
domain >> Lx;
domain >> Ly;
domain >> Lz;
ifstream image("Segmented.in");
image >> Filename; // Name of data file containing segmented data
image >> Nx; // size of the binary file
image >> Ny;
image >> Nz;
image >> xStart; // offset for the starting voxel
image >> yStart;
image >> zStart;
}
MPI_Barrier(comm);
// Computational domain
//.................................................
MPI_Bcast(&nx,1,MPI_INT,0,comm);
MPI_Bcast(&ny,1,MPI_INT,0,comm);
MPI_Bcast(&nz,1,MPI_INT,0,comm);
MPI_Bcast(&nprocx,1,MPI_INT,0,comm);
MPI_Bcast(&nprocy,1,MPI_INT,0,comm);
MPI_Bcast(&nprocz,1,MPI_INT,0,comm);
MPI_Bcast(&nspheres,1,MPI_INT,0,comm);
MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Lz,1,MPI_DOUBLE,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(&xStart,1,MPI_INT,0,comm);
MPI_Bcast(&yStart,1,MPI_INT,0,comm);
MPI_Bcast(&zStart,1,MPI_INT,0,comm);
//.................................................
MPI_Barrier(comm);
// Check that the number of processors >= the number of ranks
if ( rank==0 ) {
printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
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);
}
if ( nprocs < nprocx*nprocy*nprocz ){
ERROR("Insufficient number of processors");
}
char *SegData = NULL;
// Rank=0 reads the entire segmented data and distributes to worker processes
if (rank==0){
printf("Dimensions of segmented image: %i x %i x %i \n",Nx,Ny,Nz);
SegData = new char[Nx*Ny*Nz];
FILE *SEGDAT = fopen(Filename,"rb");
if (SEGDAT==NULL) ERROR("Error reading segmented data");
size_t ReadSeg;
ReadSeg=fread(SegData,1,Nx*Ny*Nz,SEGDAT);
if (ReadSeg != size_t(Nx*Ny*Nz)) printf("lbpm_segmented_decomp: Error reading segmented data (rank=%i)\n",rank);
fclose(SEGDAT);
printf("Read segmented data from %s \n",Filename);
}
MPI_Barrier(comm);
// Get the rank info
int N = (nx+2)*(ny+2)*(nz+2);
Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
n = k*(nx+2)*(ny+2)+j*(nx+2)+i;
Dm.id[n] = 1;
}
}
}
}
Dm.CommInit(comm);
Dm.CommInit(comm);
// number of sites to use for periodic boundary condition transition zone
int z_transition_size = (nprocz*nz - (Nz - zStart))/2;
if (z_transition_size < 0) z_transition_size=0;
// number of sites to use for periodic boundary condition transition zone
int z_transition_size = (nprocz*nz - (Nz - zStart))/2;
if (z_transition_size < 0) z_transition_size=0;
// 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("Subdomain size: %i \n",N);
printf("Size of transition region: %i \n", z_transition_size);
char *tmp;
tmp = new char[N];
for (int kp=0; kp<nprocz; kp++){
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;
// Pack and send the subdomain for rnk
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
int x = xStart + ip*nx + i-1;
int y = yStart + jp*ny + j-1;
// int z = zStart + kp*nz + k-1;
int z = zStart + kp*nz + k-1 - z_transition_size;
if (z<zStart) z=zStart;
if (!(z<Nz)) z=Nz-1;
int nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
int nglobal = z*Nx*Ny+y*Nx+x;
tmp[nlocal] = SegData[nglobal];
}
}
}
if (rnk==0){
// 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("Subdomain size: %i \n",N);
printf("Size of transition region: %i \n", z_transition_size);
char *tmp;
tmp = new char[N];
for (int kp=0; kp<nprocz; kp++){
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;
// Pack and send the subdomain for rnk
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
int x = xStart + ip*nx + i-1;
int y = yStart + jp*ny + j-1;
// int z = zStart + kp*nz + k-1;
int z = zStart + kp*nz + k-1 - z_transition_size;
if (z<zStart) z=zStart;
if (!(z<Nz)) z=Nz-1;
int nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
Dm.id[nlocal] = tmp[nlocal];
int nglobal = z*Nx*Ny+y*Nx+x;
tmp[nlocal] = SegData[nglobal];
}
}
}
}
else{
printf("Sending data to process %i \n", rnk);
MPI_Send(tmp,N,MPI_CHAR,rnk,15,comm);
if (rnk==0){
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
int nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
Dm.id[nlocal] = tmp[nlocal];
}
}
}
}
else{
printf("Sending data to process %i \n", rnk);
MPI_Send(tmp,N,MPI_CHAR,rnk,15,comm);
}
}
}
}
}
}
else{
// Recieve the subdomain from rank = 0
printf("Ready to recieve data %i at process %i \n", N,rank);
MPI_Recv(Dm.id,N,MPI_CHAR,0,15,comm,MPI_STATUS_IGNORE);
else{
// Recieve the subdomain from rank = 0
printf("Ready to recieve data %i at process %i \n", N,rank);
MPI_Recv(Dm.id,N,MPI_CHAR,0,15,comm,MPI_STATUS_IGNORE);
}
MPI_Barrier(comm);
nx+=2; ny+=2; nz+=2;
N=nx*ny*nz;
if (rank==0) printf("All sub-domains recieved \n");
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;
if (Dm.id[n]==char(SOLID)) Dm.id[n] = 0;
else if (Dm.id[n]==char(NWP)) Dm.id[n] = 1;
else Dm.id[n] = 2;
}
}
}
if (rank==0) printf("Domain set \n");
int count = 0;
int total = 0;
int countGlobal = 0;
int totalGlobal = 0;
for (k=1;k<nz-1;k++){
for (j=1;j<ny-1;j++){
for (i=1;i<nx-1;i++){
n=k*nx*ny+j*nx+i;
total++;
if (Dm.id[n] == 0){
count++;
}
}
}
}
MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
MPI_Allreduce(&total,&totalGlobal,1,MPI_INT,MPI_SUM,comm);
float porosity = float(totalGlobal-countGlobal)/totalGlobal;
if (rank==0) printf("Porosity=%f\n",porosity);
if (rank==0){
int xstart = xStart; // Is this correct?
int ystart = yStart;
int zstart = zStart;
//totalGlobal=(Nx-xstart)*(Ny-ystart)*(Nz-zstart);
countGlobal = 0;
for (k=zstart; k<zstart+nprocz*(nz-2); k++){
for (j=ystart; j<ystart+nprocy*(ny-2); j++){
for (i=xstart; i<xstart+nprocx*(nx-2); i++){
n=k*Nx*Ny+j*Nx+i;
if (n < Nx*Ny*Nz){
if (SegData[n] == char(SOLID)){
countGlobal++;
}
}
}
}
}
float porosity = float(totalGlobal-countGlobal)/totalGlobal;
printf("Original Porosity=%f\n",porosity);
}
count = 0;
total = 0;
countGlobal = 0;
totalGlobal = 0;
for (k=1;k<nz-1;k++){
for (j=1;j<ny-1;j++){
for (i=1;i<nx-1;i++){
n=k*nx*ny+j*nx+i;
if (Dm.id[n] != 0) total++;
if (Dm.id[n] == 2) count++;
}
}
}
MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
MPI_Allreduce(&total,&totalGlobal,1,MPI_INT,MPI_SUM,comm);
float saturation = float(countGlobal)/totalGlobal;
if (rank==0) printf("wetting phase saturation=%f\n",saturation);
char LocalRankFilename[40];
sprintf(LocalRankFilename,"ID.%05i",rank+rank_offset);
FILE *ID = fopen(LocalRankFilename,"wb");
fwrite(Dm.id,1,N,ID);
// fwrite(Distance.get(),8,Distance.length(),ID);
fclose(ID);
if (!MULTINPUT){
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;
// DoubleArray SymDist(nx,ny,nz);
char *symid;
symid = new char [N];
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;
int nsym=(nz-k-1)*nx*ny+j*nx+i;
symid[nsym] = Dm.id[n];
//SymDist(i,j,nz-k-1)=Distance(i,j,k);
}
}
}
sprintf(LocalRankFilename,"ID.%05i",symrank);
FILE *SYMID = fopen(LocalRankFilename,"wb");
// fwrite(SymDist.get(),8,SymDist.length(),SYMDIST);
fwrite(symid,1,N,SYMID);
fclose(SYMID);
}
}
MPI_Barrier(comm);
nx+=2; ny+=2; nz+=2;
N=nx*ny*nz;
if (rank==0) printf("All sub-domains recieved \n");
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;
if (Dm.id[n]==char(SOLID)) Dm.id[n] = 0;
else if (Dm.id[n]==char(NWP)) Dm.id[n] = 1;
else Dm.id[n] = 2;
}
}
}
if (rank==0) printf("Domain set \n");
int count = 0;
int total = 0;
int countGlobal = 0;
int totalGlobal = 0;
for (k=1;k<nz-1;k++){
for (j=1;j<ny-1;j++){
for (i=1;i<nx-1;i++){
n=k*nx*ny+j*nx+i;
total++;
if (Dm.id[n] == 0){
count++;
}
}
}
}
MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
MPI_Allreduce(&total,&totalGlobal,1,MPI_INT,MPI_SUM,comm);
float porosity = float(totalGlobal-countGlobal)/totalGlobal;
if (rank==0) printf("Porosity=%f\n",porosity);
if (rank==0){
int xstart = xStart; // Is this correct?
int ystart = yStart;
int zstart = zStart;
//totalGlobal=(Nx-xstart)*(Ny-ystart)*(Nz-zstart);
countGlobal = 0;
for (k=zstart; k<zstart+nprocz*(nz-2); k++){
for (j=ystart; j<ystart+nprocy*(ny-2); j++){
for (i=xstart; i<xstart+nprocx*(nx-2); i++){
n=k*Nx*Ny+j*Nx+i;
if (n < Nx*Ny*Nz){
if (SegData[n] == char(SOLID)){
countGlobal++;
}
}
}
}
}
float porosity = float(totalGlobal-countGlobal)/totalGlobal;
printf("Original Porosity=%f\n",porosity);
}
count = 0;
total = 0;
countGlobal = 0;
totalGlobal = 0;
for (k=1;k<nz-1;k++){
for (j=1;j<ny-1;j++){
for (i=1;i<nx-1;i++){
n=k*nx*ny+j*nx+i;
if (Dm.id[n] != 0) total++;
if (Dm.id[n] == 2) count++;
}
}
}
MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
MPI_Allreduce(&total,&totalGlobal,1,MPI_INT,MPI_SUM,comm);
float saturation = float(countGlobal)/totalGlobal;
if (rank==0) printf("wetting phase saturation=%f\n",saturation);
char LocalRankFilename[40];
sprintf(LocalRankFilename,"ID.%05i",rank);
FILE *ID = fopen(LocalRankFilename,"wb");
fwrite(Dm.id,1,N,ID);
// fwrite(Distance.get(),8,Distance.length(),ID);
fclose(ID);
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;
// DoubleArray SymDist(nx,ny,nz);
char *symid;
symid = new char [N];
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;
int nsym=(nz-k-1)*nx*ny+j*nx+i;
symid[nsym] = Dm.id[n];
//SymDist(i,j,nz-k-1)=Distance(i,j,k);
}
}
}
sprintf(LocalRankFilename,"ID.%05i",symrank);
FILE *SYMID = fopen(LocalRankFilename,"wb");
// fwrite(SymDist.get(),8,SymDist.length(),SYMDIST);
fwrite(symid,1,N,SYMID);
fclose(SYMID);
}
MPI_Barrier(comm);
MPI_Finalize();
MPI_Finalize();
}

156
tests/lbpm_segmented_pp.cpp
View File

@ -26,162 +26,6 @@ inline void MeanFilter(DoubleArray &Mesh){
}
}
inline double minmod(double &a, double &b){
double value;
value = a;
if ( a*b < 0.0) value=0.0;
else if (fabs(a) > fabs(b)) value = b;
return value;
}
inline double Eikonal(DoubleArray &Distance, char *ID, Domain &Dm, 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] == 1) 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);
if (fabs(GlobalMax) < 1e-5){
if (Dm.rank==0) printf("Exiting with max tolerance of 1e-5 \n");
count=timesteps;
}
}
return GlobalVar;
}
int main(int argc, char **argv)
{