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

refactor lbpm_BlobAnalysis

Browse Source
This commit is contained in:
James E McClure
2015-07-13 10:20:44 -04:00
parent b073af56cd
commit 711c9c85b6
1 changed files with 60 additions and 61 deletions
Download Patch File Download Diff File Expand all files Collapse all files

121
tests/lbpm_BlobAnalysis.cpp
View File

@@ -104,7 +104,7 @@ int main(int argc, char **argv)
//.................................................
Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
TwoPhase Averages(Dm);
// BlobTwoPhase BlobAverages(nblobs_global);
// BlobTwoPhase ComponentAverages_NWP(nblobs_global);
Nx+=2;Ny+=2;Nz+=2;
N=Nx*Ny*Nz; // number of lattice points
//.......................................................................
@@ -128,10 +128,10 @@ int main(int argc, char **argv)
//MPI_Barrier(MPI_COMM_WORLD);
if (rank == 0) cout << "Domain set." << endl;
//.......................................................................
sprintf(LocalRankFilename,"%s%s","BlobLabel.",LocalRankString);
ReadBlobFile(LocalRankFilename, Averages.BlobLabel.get(), N);
sprintf(LocalRankFilename,"%s%s","Label_NWP.",LocalRankString);
ReadBlobFile(LocalRankFilename, Averages.Label_NWP.get(), N);
MPI_Barrier(MPI_COMM_WORLD);
if (rank == 0) cout << "BlobLabel set." << endl;
if (rank == 0) cout << "Label_NWP set." << endl;
//.......................................................................
//copies of data needed to perform checkpointing from cpu
@@ -205,7 +205,6 @@ int main(int argc, char **argv)
double beta = 0.95;
Averages.SetupCubes(Dm);
Averages.UpdateSolid();
Averages.Initialize();
Averages.ComputeDelPhi();
@@ -215,9 +214,9 @@ int main(int argc, char **argv)
Averages.Reduce();
int b=0;
// Blobs.Set(Averages.BlobAverages.NBLOBS);
int dimx = Averages.BlobAverages.size(0);
int dimy = Averages.BlobAverages.size(1);
// Blobs.Set(Averages.ComponentAverages_NWP.NBLOBS);
int dimx = Averages.ComponentAverages_NWP.size(0);
int dimy = Averages.ComponentAverages_NWP.size(1);
int TotalBlobInfoSize=dimx*dimy;
FILE *BLOBLOG;
@@ -227,41 +226,41 @@ int main(int argc, char **argv)
}
// BlobContainer Blobs;
DoubleArray RecvBuffer(dimx);
// MPI_Allreduce(&Averages.BlobAverages.get(),&Blobs.get(),1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
// MPI_Allreduce(&Averages.ComponentAverages_NWP.get(),&Blobs.get(),1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Barrier(MPI_COMM_WORLD);
if (rank==0) printf("All ranks passed gate \n");
for (int b=0; b<(int)Averages.BlobAverages.size(1); b++){
for (int b=0; b<(int)Averages.ComponentAverages_NWP.size(1); b++){
MPI_Allreduce(&Averages.BlobAverages(0,b),&RecvBuffer(0),dimx,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
for (int idx=0; idx<dimx-1; idx++) Averages.BlobAverages(idx,b)=RecvBuffer(idx);
MPI_Allreduce(&Averages.ComponentAverages_NWP(0,b),&RecvBuffer(0),dimx,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
for (int idx=0; idx<dimx-1; idx++) Averages.ComponentAverages_NWP(idx,b)=RecvBuffer(idx);
MPI_Barrier(MPI_COMM_WORLD);
if (Averages.BlobAverages(0,b) > 0.0){
if (Averages.ComponentAverages_NWP(0,b) > 0.0){
double Vn,pn,awn,ans,Jwn,Kwn,lwns,cwns,trawn,trJwn;
Vn = Averages.BlobAverages(1,b);
pn = Averages.BlobAverages(2,b)/Averages.BlobAverages(0,b);
awn = Averages.BlobAverages(3,b);
ans = Averages.BlobAverages(4,b);
Vn = Averages.ComponentAverages_NWP(1,b);
pn = Averages.ComponentAverages_NWP(2,b)/Averages.ComponentAverages_NWP(0,b);
awn = Averages.ComponentAverages_NWP(3,b);
ans = Averages.ComponentAverages_NWP(4,b);
if (awn != 0.0){
Jwn = Averages.BlobAverages(5,b)/Averages.BlobAverages(3,b);
Kwn = Averages.BlobAverages(6,b)/Averages.BlobAverages(3,b);
Jwn = Averages.ComponentAverages_NWP(5,b)/Averages.ComponentAverages_NWP(3,b);
Kwn = Averages.ComponentAverages_NWP(6,b)/Averages.ComponentAverages_NWP(3,b);
}
else Jwn=Kwn=0.0;
trawn = Averages.BlobAverages(12,b);
trawn = Averages.ComponentAverages_NWP(12,b);
if (trawn != 0.0){
trJwn = Averages.BlobAverages(13,b)/trawn;
trJwn = Averages.ComponentAverages_NWP(13,b)/trawn;
}
else trJwn=0.0;
lwns = Averages.BlobAverages(7,b);
if (lwns != 0.0) cwns = Averages.BlobAverages(8,b)/Averages.BlobAverages(7,b);
lwns = Averages.ComponentAverages_NWP(7,b);
if (lwns != 0.0) cwns = Averages.ComponentAverages_NWP(8,b)/Averages.ComponentAverages_NWP(7,b);
else cwns=0.0;
Averages.BlobAverages(2,b) = pn;
Averages.BlobAverages(5,b) = trJwn;
Averages.BlobAverages(6,b) = Kwn;
Averages.BlobAverages(8,b) = cwns;
// Averages.BlobAverages(13,b) = trJwn;
Averages.ComponentAverages_NWP(2,b) = pn;
Averages.ComponentAverages_NWP(5,b) = trJwn;
Averages.ComponentAverages_NWP(6,b) = Kwn;
Averages.ComponentAverages_NWP(8,b) = cwns;
// Averages.ComponentAverages_NWP(13,b) = trJwn;
}
}
@@ -271,17 +270,17 @@ int main(int argc, char **argv)
if (rank==0){
// printf("Reduced blob %i \n",b);
fprintf(BLOBLOG,"%.5g %.5g %.5g\n",Averages.vol_w_global,Averages.paw_global,Averages.aws_global);
for (int b=0; b<(int)Averages.BlobAverages.size(1); b++){
if (Averages.BlobAverages(0,b) > 0.0){
for (int b=0; b<(int)Averages.ComponentAverages_NWP.size(1); b++){
if (Averages.ComponentAverages_NWP(0,b) > 0.0){
double Vn,pn,awn,ans,Jwn,Kwn,lwns,cwns;
Vn = Averages.BlobAverages(1,b);
pn = Averages.BlobAverages(2,b);
awn = Averages.BlobAverages(3,b);
ans = Averages.BlobAverages(4,b);
Jwn = Averages.BlobAverages(5,b);
Kwn = Averages.BlobAverages(6,b);
lwns = Averages.BlobAverages(7,b);
cwns = Averages.BlobAverages(8,b);
Vn = Averages.ComponentAverages_NWP(1,b);
pn = Averages.ComponentAverages_NWP(2,b);
awn = Averages.ComponentAverages_NWP(3,b);
ans = Averages.ComponentAverages_NWP(4,b);
Jwn = Averages.ComponentAverages_NWP(5,b);
Kwn = Averages.ComponentAverages_NWP(6,b);
lwns = Averages.ComponentAverages_NWP(7,b);
cwns = Averages.ComponentAverages_NWP(8,b);
fprintf(BLOBLOG,"%.5g ", Vn); //Vn
fprintf(BLOBLOG,"%.5g ", pn); //pn
@@ -318,37 +317,37 @@ inline void WriteBlobStates(TwoPhase TCAT, double D, double porosity){
pw = TCAT.paw_global;
aws = TCAT.aws;
// Compute the averages over the entire non-wetting phsae
for (a=0; a<(int)TCAT.BlobAverages.size(1); a++){
vol_n += TCAT.BlobAverages(0,a);
pan += TCAT.BlobAverages(2,a)*TCAT.BlobAverages(0,a);
awn += TCAT.BlobAverages(3,a);
ans += TCAT.BlobAverages(4,a);
Jwn += TCAT.BlobAverages(5,a)*TCAT.BlobAverages(3,a);
Kwn += TCAT.BlobAverages(6,a)*TCAT.BlobAverages(3,a);
lwns += TCAT.BlobAverages(7,a);
clwns += TCAT.BlobAverages(8,a)*TCAT.BlobAverages(7,a);
nwp_volume += TCAT.BlobAverages(1,a);
pawn += TCAT.BlobAverages(2,a)*TCAT.BlobAverages(3,a);
for (a=0; a<(int)TCAT.ComponentAverages_NWP.size(1); a++){
vol_n += TCAT.ComponentAverages_NWP(0,a);
pan += TCAT.ComponentAverages_NWP(2,a)*TCAT.ComponentAverages_NWP(0,a);
awn += TCAT.ComponentAverages_NWP(3,a);
ans += TCAT.ComponentAverages_NWP(4,a);
Jwn += TCAT.ComponentAverages_NWP(5,a)*TCAT.ComponentAverages_NWP(3,a);
Kwn += TCAT.ComponentAverages_NWP(6,a)*TCAT.ComponentAverages_NWP(3,a);
lwns += TCAT.ComponentAverages_NWP(7,a);
clwns += TCAT.ComponentAverages_NWP(8,a)*TCAT.ComponentAverages_NWP(7,a);
nwp_volume += TCAT.ComponentAverages_NWP(1,a);
pawn += TCAT.ComponentAverages_NWP(2,a)*TCAT.ComponentAverages_NWP(3,a);
}
// Compute the pore voume (sum of wetting an non-wetting phase volumes)
PoreVolume=TCAT.wp_volume_global + nwp_volume;
// Subtract off portions of non-wetting phase in order of size
for (a=TCAT.BlobAverages.size(1)-1; a>0; a--){
for (a=TCAT.ComponentAverages_NWP.size(1)-1; a>0; a--){
// Subtract the features one-by-one
vol_n -= TCAT.BlobAverages(0,a);
pan -= TCAT.BlobAverages(2,a)*TCAT.BlobAverages(0,a);
awn -= TCAT.BlobAverages(3,a);
ans -= TCAT.BlobAverages(4,a);
Jwn -= TCAT.BlobAverages(5,a)*TCAT.BlobAverages(3,a);
Kwn -= TCAT.BlobAverages(6,a)*TCAT.BlobAverages(3,a);
lwns -= TCAT.BlobAverages(7,a);
clwns -= TCAT.BlobAverages(8,a)*TCAT.BlobAverages(7,a);
nwp_volume -= TCAT.BlobAverages(1,a);
pawn -= TCAT.BlobAverages(2,a)*TCAT.BlobAverages(3,a);
vol_n -= TCAT.ComponentAverages_NWP(0,a);
pan -= TCAT.ComponentAverages_NWP(2,a)*TCAT.ComponentAverages_NWP(0,a);
awn -= TCAT.ComponentAverages_NWP(3,a);
ans -= TCAT.ComponentAverages_NWP(4,a);
Jwn -= TCAT.ComponentAverages_NWP(5,a)*TCAT.ComponentAverages_NWP(3,a);
Kwn -= TCAT.ComponentAverages_NWP(6,a)*TCAT.ComponentAverages_NWP(3,a);
lwns -= TCAT.ComponentAverages_NWP(7,a);
clwns -= TCAT.ComponentAverages_NWP(8,a)*TCAT.ComponentAverages_NWP(7,a);
nwp_volume -= TCAT.ComponentAverages_NWP(1,a);
pawn -= TCAT.ComponentAverages_NWP(2,a)*TCAT.ComponentAverages_NWP(3,a);
// Update wetting phase averages
aws += TCAT.BlobAverages(4,a);
aws += TCAT.ComponentAverages_NWP(4,a);
if (vol_n > 64){ // Only consider systems with "large enough" blobs -- 4^3
if (fabs(1.0 - nwp_volume/PoreVolume - sw) > 0.005 || a == 1){
sw = 1.0 - nwp_volume/PoreVolume;