Factored random phase initialization into pre-processor tests/lbpm_random_pp
This commit is contained in:
James E McClure
parent
b551af3764
commit
bb1d2658d6
@ -1127,102 +1127,6 @@ inline void SignedDistance(double *Distance, int nspheres, double *List_cx, doub
|
||||
for (n=0; n<N; n++) Distance[n] = Distance[n]/hx;
|
||||
}
|
||||
|
||||
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++){
|
||||
if (ID[n] == 1) ID[n] = 2;
|
||||
else if (ID[n] == 2) ID[n] = 1;
|
||||
}
|
||||
}
|
||||
|
||||
inline void WriteLocalSolidID(char *FILENAME, char *ID, int N)
|
||||
{
|
||||
char value;
|
||||
|
||||
@ -23,6 +23,104 @@
|
||||
* James E. McClure 2013-2014
|
||||
*/
|
||||
|
||||
|
||||
// The functions GenerateResidual and FlipID may be deprecated
|
||||
// see lbpm_random_pp for latest version
|
||||
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++){
|
||||
if (ID[n] == 1) ID[n] = 2;
|
||||
else if (ID[n] == 2) ID[n] = 1;
|
||||
}
|
||||
}
|
||||
|
||||
#define USE_NEW_WRITER
|
||||
|
||||
using namespace std;
|
||||
@ -487,7 +585,7 @@ int main(int argc, char **argv)
|
||||
id[0] = id[Nx-1] = id[(Ny-1)*Nx] = id[(Ny-1)*Nx + Nx-1] = 0;
|
||||
id[(Nz-1)*Nx*Ny] = id[(Nz-1)*Nx*Ny+Nx-1] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx + Nx-1] = 0;
|
||||
//.........................................................
|
||||
|
||||
/*
|
||||
#ifdef USE_EXP_CONTACT_ANGLE
|
||||
// If negative phi_s is chosen, flip the ID for the wetting and non-wetting phase
|
||||
if (phi_s < 0.0 && !pBC){
|
||||
@ -508,7 +606,7 @@ int main(int argc, char **argv)
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
*/
|
||||
// Initialize communication structures in averaging domain
|
||||
for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = id[i];
|
||||
Dm.CommInit(MPI_COMM_WORLD);
|
||||
|
||||
@ -372,9 +372,9 @@ int main(int argc, char **argv)
|
||||
}
|
||||
|
||||
// Generate the residual NWP
|
||||
if (BoundaryCondition==0 && InitialCondition==0 && rank==0) printf("Initializing with NWP saturation = %f \n",wp_saturation);
|
||||
if (BoundaryCondition==0 && InitialCondition==0) GenerateResidual(id,Nx,Ny,Nz,wp_saturation);
|
||||
if (InitialCondition == 2){
|
||||
// if (BoundaryCondition==0 && InitialCondition==0 && rank==0) printf("Initializing with NWP saturation = %f \n",wp_saturation);
|
||||
// if (BoundaryCondition==0 && InitialCondition==0) GenerateResidual(id,Nx,Ny,Nz,wp_saturation);
|
||||
// if (InitialCondition == 2){
|
||||
if (rank==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
|
||||
sprintf(LocalRankFilename,"ID.%05i",rank);
|
||||
FILE *IDFILE = fopen(LocalRankFilename,"rb");
|
||||
@ -382,7 +382,7 @@ int main(int argc, char **argv)
|
||||
fread(id,1,N,IDFILE);
|
||||
fclose(IDFILE);
|
||||
// CopyToDevice(ID, id, N);
|
||||
}
|
||||
// }
|
||||
|
||||
|
||||
// Set up kstart, kfinish so that the reservoirs are excluded from averaging
|
||||
@ -438,7 +438,7 @@ int main(int argc, char **argv)
|
||||
id[(Nz-1)*Nx*Ny] = id[(Nz-1)*Nx*Ny+Nx-1] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx + Nx-1] = 0;
|
||||
//.........................................................
|
||||
|
||||
// If positive phi_s is chosen, flip the ID for the wetting and non-wetting phase
|
||||
/* // If positive phi_s is chosen, flip the ID for the wetting and non-wetting phase
|
||||
if (phi_s > 0.0 && BoundaryCondition==0){
|
||||
phi_s = -phi_s;
|
||||
das = (phi_s+1.0)*0.5;
|
||||
@ -446,6 +446,7 @@ int main(int argc, char **argv)
|
||||
if (rank == 0) printf("Resetting phi_s = %f, das = %f, dbs = %f \n", phi_s, das, dbs);
|
||||
FlipID(id,Nx*Ny*Nz);
|
||||
}
|
||||
*/
|
||||
// Initialize communication structures in averaging domain
|
||||
for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = id[i];
|
||||
Dm.CommInit(MPI_COMM_WORLD);
|
||||
|
||||
218
tests/lbpm_random_pp.cpp
Normal file
218
tests/lbpm_random_pp.cpp
Normal file
@ -0,0 +1,218 @@
|
||||
/*
|
||||
* Pre-processor to generate signed distance function from segmented data
|
||||
* segmented data should be stored in a raw binary file as 1-byte integer (type char)
|
||||
* will output distance functions for phases
|
||||
*/
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <math.h>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <sstream>
|
||||
#include <Array.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++){
|
||||
if (ID[n] == 1) ID[n] = 2;
|
||||
else if (ID[n] == 2) ID[n] = 1;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
// Initialize MPI
|
||||
int rank, nprocs;
|
||||
MPI_Init(&argc,&argv);
|
||||
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
|
||||
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
|
||||
|
||||
int InitialWetting;
|
||||
double Saturation;
|
||||
if (argc == 3){
|
||||
sscanf(argv[1],"%lf",&Saturation);
|
||||
sscanf(argv[2],"%d",&InitialWetting);
|
||||
if (rank==0){
|
||||
printf("Initializing wetting phase saturation of %f \n",Saturation);
|
||||
if (InitialWetting == 1)
|
||||
printf("Begin from connected wetting phase \n");
|
||||
else
|
||||
printf("Begin from connected non-wetting phase \n");
|
||||
}
|
||||
}
|
||||
|
||||
if (InitialWetting == 1) Saturation=1.0-Saturation;
|
||||
//.......................................................................
|
||||
// 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;
|
||||
|
||||
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;
|
||||
|
||||
}
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
// Computational domain
|
||||
MPI_Bcast(&nx,1,MPI_INT,0,MPI_COMM_WORLD);
|
||||
MPI_Bcast(&ny,1,MPI_INT,0,MPI_COMM_WORLD);
|
||||
MPI_Bcast(&nz,1,MPI_INT,0,MPI_COMM_WORLD);
|
||||
MPI_Bcast(&nprocx,1,MPI_INT,0,MPI_COMM_WORLD);
|
||||
MPI_Bcast(&nprocy,1,MPI_INT,0,MPI_COMM_WORLD);
|
||||
MPI_Bcast(&nprocz,1,MPI_INT,0,MPI_COMM_WORLD);
|
||||
MPI_Bcast(&nspheres,1,MPI_INT,0,MPI_COMM_WORLD);
|
||||
MPI_Bcast(&Lx,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
|
||||
MPI_Bcast(&Ly,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
|
||||
MPI_Bcast(&Lz,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
|
||||
//.................................................
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
|
||||
// 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 LocalRankFilename[40];
|
||||
|
||||
nx+=2; ny+=2; nz+=2;
|
||||
int N = nx*ny*nz;
|
||||
char *id;
|
||||
id = new char[N];
|
||||
|
||||
DoubleArray SignDist(nx,ny,nz);
|
||||
// Read the signed distance from file
|
||||
sprintf(LocalRankFilename,"SignDist.%05i",rank);
|
||||
FILE *DIST = fopen(LocalRankFilename,"rb");
|
||||
fread(SignDist.get(),8,N,DIST);
|
||||
fclose(DIST);
|
||||
|
||||
for (int k=0; k<nz; k++){
|
||||
for (int j=0; j<ny; j++){
|
||||
for (int i=0; i<nx; i++){
|
||||
if (SignDist(i,j,k) < 0.0) id[n] = 0;
|
||||
else id[n] = 2;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Generate the residual NWP
|
||||
if (rank==0) printf("Initializing with NWP saturation = %f \n",wp_saturation);
|
||||
GenerateResidual(id,Nx,Ny,Nz,wp_saturation);
|
||||
|
||||
if (InitialWetting == 1) FlipID(id,Nx*Ny*Nz);
|
||||
|
||||
sprintf(LocalRankFilename,"ID.%05i",rank);
|
||||
FILE *ID = fopen(LocalRankFilename,"wb");
|
||||
fread(id,1,N,ID);
|
||||
fclose(ID);
|
||||
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
MPI_Finalize();
|
||||
return 0;
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user