incorporating lbpm_random with fixed saturation
This commit is contained in:
James E McClure
parent
c543176244
commit
0acf9c5146
@ -12,93 +12,6 @@
|
||||
#include <Array.h>
|
||||
#include <Domain.h>
|
||||
|
||||
inline void GenerateResidual(char *ID, int Nx, int Ny, int Nz, double Saturation)
|
||||
{
|
||||
//.......................................................................
|
||||
int i,j,k,n,Number,N;
|
||||
int x,y,z,ii,jj,kk;
|
||||
int sizeX,sizeY,sizeZ;
|
||||
int *SizeX, *SizeY, *SizeZ;
|
||||
|
||||
#ifdef NORANDOM
|
||||
srand(10009);
|
||||
#else
|
||||
srand(time(NULL));
|
||||
#endif
|
||||
// float bin;
|
||||
//.......................................................................
|
||||
N = Nx*Ny*Nz;
|
||||
|
||||
int bin, binCount;
|
||||
ifstream Dist("BlobSize.in");
|
||||
Dist >> binCount;
|
||||
// printf("Number of blob sizes: %i \n",binCount);
|
||||
SizeX = new int [binCount];
|
||||
SizeY = new int [binCount];
|
||||
SizeZ = new int [binCount];
|
||||
for (bin=0; bin<binCount; bin++){
|
||||
Dist >> SizeX[bin];
|
||||
Dist >> SizeY[bin];
|
||||
Dist >> SizeZ[bin];
|
||||
// printf("Blob %i dimension: %i x %i x %i \n",bin, SizeX[bin], SizeY[bin], SizeZ[bin]);
|
||||
}
|
||||
Dist.close();
|
||||
//.......................................................................
|
||||
// cout << "Generating blocks... " << endl;
|
||||
// Count for the total number of oil nodes
|
||||
int count = 0;
|
||||
// Count the total number of non-solid nodes
|
||||
int total = 0;
|
||||
for (i=0;i<N;i++){
|
||||
if (ID[i] != 0) total++;
|
||||
}
|
||||
|
||||
float sat = 0.f;
|
||||
Number = 0; // number of features
|
||||
while (sat < Saturation){
|
||||
Number++;
|
||||
// Randomly generate a point in the domain
|
||||
x = Nx*float(rand())/float(RAND_MAX);
|
||||
y = Ny*float(rand())/float(RAND_MAX);
|
||||
z = Nz*float(rand())/float(RAND_MAX);
|
||||
|
||||
bin = int(floor(binCount*float(rand())/float(RAND_MAX)));
|
||||
sizeX = SizeX[bin];
|
||||
sizeY = SizeY[bin];
|
||||
sizeZ = SizeZ[bin];
|
||||
|
||||
// cout << "Sampling from bin no. " << floor(bin) << endl;
|
||||
// cout << "Feature size is: " << sizeX << "x" << sizeY << "x" << sizeZ << endl;
|
||||
|
||||
for (k=z;k<z+sizeZ;k++){
|
||||
for (j=y;j<y+sizeY;j++){
|
||||
for (i=x;i<x+sizeX;i++){
|
||||
// Identify nodes in the domain (periodic BC)
|
||||
ii = i;
|
||||
jj = j;
|
||||
kk = k;
|
||||
if (ii < 1) ii+=(Nx-2);
|
||||
if (jj < 1) jj+=(Ny-2);
|
||||
if (kk < 1) kk+=(Nz-2);
|
||||
if (!(ii < Nx-1)) ii-=(Nx-2);
|
||||
if (!(jj < Ny-1)) jj-=(Ny-2);
|
||||
if (!(kk < Nz-1)) kk-=(Nz-2);
|
||||
|
||||
n = kk*Nx*Ny+jj*Nx+ii;
|
||||
|
||||
if (ID[n] == 2){
|
||||
ID[n] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
sat = float(count)/total;
|
||||
}
|
||||
//.......................................................................
|
||||
}
|
||||
|
||||
|
||||
inline void FlipID(char *ID, int N)
|
||||
{
|
||||
for (int n=0; n<N; n++){
|
||||
@ -182,6 +95,9 @@ int main(int argc, char **argv)
|
||||
|
||||
char LocalRankFilename[40];
|
||||
|
||||
int BoundaryCondition=0;
|
||||
Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
|
||||
|
||||
nx+=2; ny+=2; nz+=2;
|
||||
int N = nx*ny*nz;
|
||||
char *id;
|
||||
@ -194,22 +110,147 @@ int main(int argc, char **argv)
|
||||
fread(SignDist.get(),8,N,DIST);
|
||||
fclose(DIST);
|
||||
|
||||
int count,countGlobal,totalGlobal;
|
||||
count = 0;
|
||||
for (int k=0; k<nz; k++){
|
||||
for (int j=0; j<ny; j++){
|
||||
for (int i=0; i<nx; i++){
|
||||
n = k*nx*ny+j*nx+i;
|
||||
n = k*nx*ny+j*nx+i;
|
||||
if (SignDist(i,j,k) < 0.0) id[n] = 0;
|
||||
else id[n] = 2;
|
||||
else{
|
||||
id[n] = 2;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// total Global is the number of nodes in the pore-space
|
||||
MPI_Allreduce(&count,&totalGlobal,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
|
||||
|
||||
Dm.CommInit(MPI_COMM_WORLD);
|
||||
int iproc = Dm.iproc;
|
||||
int jproc = Dm.jproc;
|
||||
int kproc = Dm.kproc;
|
||||
|
||||
int bin, binCount;
|
||||
ifstream Dist("BlobSize.in");
|
||||
Dist >> binCount;
|
||||
int *SizeX, *SizeY, *SizeZ;
|
||||
// printf("Number of blob sizes: %i \n",binCount);
|
||||
SizeX = new int [binCount];
|
||||
SizeY = new int [binCount];
|
||||
SizeZ = new int [binCount];
|
||||
for (bin=0; bin<binCount; bin++){
|
||||
Dist >> SizeX[bin];
|
||||
Dist >> SizeY[bin];
|
||||
Dist >> SizeZ[bin];
|
||||
// printf("Blob %i dimension: %i x %i x %i \n",bin, SizeX[bin], SizeY[bin], SizeZ[bin]);
|
||||
}
|
||||
Dist.close();
|
||||
|
||||
// Generate the residual NWP
|
||||
if (rank==0) printf("Initializing with NWP saturation = %f \n",Saturation);
|
||||
GenerateResidual(id,nx,ny,nz,Saturation);
|
||||
// GenerateResidual(id,nx,ny,nz,Saturation);
|
||||
|
||||
int x,y,z;
|
||||
int sizeX,sizeY,sizeZ;
|
||||
int ii,jj,kk;
|
||||
int Nx = nx;
|
||||
int Ny = ny;
|
||||
int Nz = nz;
|
||||
float sat = 0.f;
|
||||
int Number = 0; // number of features
|
||||
while (sat < Saturation){
|
||||
if (rank==0){
|
||||
Number++;
|
||||
// Randomly generate a point in the domain
|
||||
x = (Nx-2)*nprocx*float(rand())/float(RAND_MAX);
|
||||
y = (Ny-2)*nprocy*float(rand())/float(RAND_MAX);
|
||||
z = (Nz-2)*nprocz*float(rand())/float(RAND_MAX);
|
||||
|
||||
bin = int(floor(binCount*float(rand())/float(RAND_MAX)));
|
||||
sizeX = SizeX[bin];
|
||||
sizeY = SizeY[bin];
|
||||
sizeZ = SizeZ[bin];
|
||||
}
|
||||
MPI_Bcast(&x,1,MPI_INT,0,MPI_COMM_WORLD);
|
||||
MPI_Bcast(&y,1,MPI_INT,0,MPI_COMM_WORLD);
|
||||
MPI_Bcast(&z,1,MPI_INT,0,MPI_COMM_WORLD);
|
||||
MPI_Bcast(&sizeX,1,MPI_INT,0,MPI_COMM_WORLD);
|
||||
MPI_Bcast(&sizeY,1,MPI_INT,0,MPI_COMM_WORLD);
|
||||
MPI_Bcast(&sizeZ,1,MPI_INT,0,MPI_COMM_WORLD);
|
||||
|
||||
//if (rank==0) printf("Broadcast block at %i,%i,%i \n",x,y,z);
|
||||
|
||||
for (k=z;k<z+sizeZ;k++){
|
||||
for (j=y;j<y+sizeY;j++){
|
||||
for (i=x;i<x+sizeX;i++){
|
||||
|
||||
// Identify nodes in the domain (periodic BC)
|
||||
ii = i;
|
||||
jj = j;
|
||||
kk = k;
|
||||
|
||||
if (ii>nprocx*(Nx-2)) ii-=nprocx*(Nx-2);
|
||||
if (jj>nprocy*(Ny-2)) jj-=nprocy*(Ny-2);
|
||||
if (kk>nprocz*(Nz-2)) kk-=nprocz*(Nz-2);
|
||||
|
||||
// Check if this is in the subdomain
|
||||
if (ii < (iproc+1)*(Nx-2) && jj < (jproc+1)*(Ny-2) && kk < (kproc+1)*(Nz-2) &&
|
||||
ii > iproc*(Nx-2) && jj > jproc*(Ny-2) && kk > kproc*(Nz-2) ){
|
||||
|
||||
// Map from global to local coordinates
|
||||
ii -= iproc*(Nx-2);
|
||||
jj -= jproc*(Ny-2);
|
||||
kk -= kproc*(Nz-2);
|
||||
|
||||
n = kk*Nx*Ny+jj*Nx+ii;
|
||||
|
||||
if (id[n] == 2){
|
||||
id[n] = 1;
|
||||
//count++;
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
count = 0;
|
||||
for (int k=1; k<Nz-1; k++){
|
||||
for (int j=1; j<Ny-1; j++){
|
||||
for (int i=1; i<Nx-1; i++){
|
||||
n=k*Nx*Ny+j*Nx+i;
|
||||
if (id[n] == 1){
|
||||
count++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
|
||||
sat = float(countGlobal)/totalGlobal;
|
||||
//if (rank==0) printf("New count=%i\n",countGlobal);
|
||||
//if (rank==0) printf("New saturation=%f\n",sat);
|
||||
}
|
||||
|
||||
if (InitialWetting == 1) FlipID(id,nx*ny*nz);
|
||||
|
||||
|
||||
count = 0;
|
||||
for (int k=1; k<Nz-1; k++){
|
||||
for (int j=1; j<Ny-1; j++){
|
||||
for (int i=1; i<Nx-1; i++){
|
||||
n=k*Nx*Ny+j*Nx+i;
|
||||
|
||||
if (id[n] == 1){
|
||||
count++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
|
||||
sat = float(countGlobal)/totalGlobal;
|
||||
if (rank==0) printf("Final saturation=%f\n",sat);
|
||||
|
||||
sprintf(LocalRankFilename,"ID.%05i",rank);
|
||||
FILE *ID = fopen(LocalRankFilename,"wb");
|
||||
fwrite(id,1,N,ID);
|
||||
|
||||
Loading…
Reference in New Issue
Block a user