439 lines
14 KiB
C++
439 lines
14 KiB
C++
// Sequential component labeling for two phase systems
|
|
// Reads parallel simulation data and performs connectivity analysis
|
|
// and averaging on a blob-by-blob basis
|
|
// James E. McClure 2015
|
|
|
|
#include <iostream>
|
|
#include <math.h>
|
|
#include "analysis/analysis.h"
|
|
#include "analysis/TwoPhase.h"
|
|
|
|
#define NUM_AVERAGES 30
|
|
|
|
using namespace std;
|
|
|
|
|
|
inline void ReadFromRank(char *FILENAME, DoubleArray &Phase, DoubleArray &Pressure, DoubleArray &Vel_x,
|
|
DoubleArray &Vel_y, DoubleArray &Vel_z, int nx, int ny, int nz, int iproc, int
|
|
jproc, int kproc)
|
|
{
|
|
int i,j,k,q,n,N;
|
|
int iglobal,jglobal,kglobal;
|
|
double value;
|
|
double denA,denB;
|
|
double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9;
|
|
double f10,f11,f12,f13,f14,f15,f16,f17,f18;
|
|
double vx,vy,vz;
|
|
|
|
N = nx*ny*nz;
|
|
|
|
double *Den, *DistEven, *DistOdd;
|
|
|
|
Den = new double[2*N];
|
|
DistEven = new double[10*N];
|
|
DistOdd = new double[9*N];
|
|
|
|
ifstream File(FILENAME,ios::binary);
|
|
for (n=0; n<N; n++){
|
|
// Write the two density values
|
|
File.read((char*) &value, sizeof(value));
|
|
Den[n] = value;
|
|
// if (n== 66276) printf("Density a = %f \n",value);
|
|
File.read((char*) &value, sizeof(value));
|
|
Den[N+n] = value;
|
|
|
|
// if (n== 66276) printf("Density b = %f \n",value);
|
|
// Read the even distributions
|
|
for (q=0; q<10; q++){
|
|
File.read((char*) &value, sizeof(value));
|
|
DistEven[q*N+n] = value;
|
|
}
|
|
// Read the odd distributions
|
|
for (q=0; q<9; q++){
|
|
File.read((char*) &value, sizeof(value));
|
|
DistOdd[q*N+n] = value;
|
|
}
|
|
}
|
|
File.close();
|
|
|
|
// Compute the phase field, pressure and velocity
|
|
for (k=1; k<nz-1; k++){
|
|
for (j=1; j<ny-1; j++){
|
|
for (i=1; i<nz-1; i++){
|
|
//........................................................................
|
|
n = k*nx*ny+j*nx+i;
|
|
//........................................................................
|
|
denA = Den[n];
|
|
denB = Den[N+n];
|
|
//........................................................................
|
|
f0 = DistEven[n];
|
|
f2 = DistEven[N+n];
|
|
f4 = DistEven[2*N+n];
|
|
f6 = DistEven[3*N+n];
|
|
f8 = DistEven[4*N+n];
|
|
f10 = DistEven[5*N+n];
|
|
f12 = DistEven[6*N+n];
|
|
f14 = DistEven[7*N+n];
|
|
f16 = DistEven[8*N+n];
|
|
f18 = DistEven[9*N+n];
|
|
//........................................................................
|
|
f1 = DistOdd[n];
|
|
f3 = DistOdd[1*N+n];
|
|
f5 = DistOdd[2*N+n];
|
|
f7 = DistOdd[3*N+n];
|
|
f9 = DistOdd[4*N+n];
|
|
f11 = DistOdd[5*N+n];
|
|
f13 = DistOdd[6*N+n];
|
|
f15 = DistOdd[7*N+n];
|
|
f17 = DistOdd[8*N+n];
|
|
//........................................................................
|
|
//.................Compute the pressure....................................
|
|
value = 0.3333333333333333*(f0+f2+f1+f4+f3+f6+f5+f8+f7+f10+f9+f12+f11+f14+f13+f16+f15+f18+f17);
|
|
//........................................................................
|
|
//.................Compute the velocity...................................
|
|
vx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14;
|
|
vy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18;
|
|
vz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18;
|
|
//........................................................................
|
|
// save values in global arrays
|
|
//........................................................................
|
|
iglobal = iproc*(nx-2)+i;
|
|
jglobal = jproc*(ny-2)+j;
|
|
kglobal = kproc*(nz-2)+k;
|
|
//........................................................................
|
|
Phase(iglobal,jglobal,kglobal) = (denA-denB)/(denA+denB);
|
|
Pressure(iglobal,jglobal,kglobal) = value;
|
|
Vel_x(iglobal,jglobal,kglobal) = vx;
|
|
Vel_y(iglobal,jglobal,kglobal) = vy;
|
|
Vel_z(iglobal,jglobal,kglobal) = vz;
|
|
//........................................................................
|
|
}
|
|
}
|
|
}
|
|
|
|
delete Den;
|
|
delete DistEven;
|
|
delete DistOdd;
|
|
}
|
|
|
|
int main(int argc, char **argv)
|
|
{
|
|
// Initialize MPI
|
|
Utilities::startup( argc, argv );
|
|
Utilities::MPI comm( MPI_COMM_WORLD );
|
|
int rank = comm.getRank();
|
|
int nprocs = comm.getSize();
|
|
printf("----------------------------------------------------------\n");
|
|
printf("COMPUTING TCAT ANALYSIS FOR NON-WETTING PHASE FEATURES \n");
|
|
printf("----------------------------------------------------------\n");
|
|
|
|
if (nprocs != 1) INSIST(nprocs == 1,"Error: ComponentLabel --serial case!");
|
|
|
|
//.......................................................................
|
|
int nprocx,nprocy,nprocz;
|
|
int Nx, Ny, Nz;
|
|
int nx,ny,nz;
|
|
int nspheres;
|
|
double Lx,Ly,Lz;
|
|
//.......................................................................
|
|
int i,j,k,n;
|
|
int iproc,jproc,kproc;
|
|
//.......................................................................
|
|
// Reading the domain information file
|
|
//.......................................................................
|
|
ifstream domain("Domain.in");
|
|
domain >> nprocx;
|
|
domain >> nprocy;
|
|
domain >> nprocz;
|
|
domain >> nx;
|
|
domain >> ny;
|
|
domain >> nz;
|
|
domain >> nspheres;
|
|
domain >> Lx;
|
|
domain >> Ly;
|
|
domain >> Lz;
|
|
//.......................................................................
|
|
|
|
nx+=2;
|
|
ny+=2;
|
|
nz+=2;
|
|
|
|
nprocs = nprocx*nprocy*nprocz;
|
|
printf("Number of MPI ranks: %i \n", nprocs);
|
|
int BoundaryCondition=0;
|
|
Nx = (nx-2)*nprocx;
|
|
Ny = (ny-2)*nprocy;
|
|
Nz = (nz-2)*nprocz;
|
|
Domain Dm(Nx,Ny,Nz,rank,1,1,1,Lx,Ly,Lz,BoundaryCondition);
|
|
Nx+=2; Ny+=2; Nz+=2;
|
|
printf("Full domain size: %i x %i x %i \n", Nx,Ny,Nz);
|
|
|
|
|
|
TwoPhase Averages(Dm);
|
|
|
|
// Filenames used
|
|
char LocalRankString[8];
|
|
char LocalRankFilename[40];
|
|
char BaseFilename[20];
|
|
sprintf(BaseFilename,"%s","dPdt.");
|
|
|
|
int proc,iglobal,kglobal,jglobal;
|
|
|
|
double * Temp;
|
|
Temp = new double[nx*ny*nz];
|
|
|
|
// read the files and populate main arrays
|
|
for ( kproc=0; kproc<nprocz; kproc++){
|
|
for ( jproc=0; jproc<nprocy; jproc++){
|
|
for ( iproc=0; iproc<nprocx; iproc++){
|
|
|
|
proc = kproc*nprocx*nprocy + jproc*nprocx + iproc;
|
|
|
|
sprintf(LocalRankString,"%05d",proc);
|
|
// sprintf(LocalRankFilename,"%s%s","dPdt.",LocalRankString);
|
|
// printf("Reading file %s \n",LocalRankFilename);
|
|
ReadBinaryFile(LocalRankFilename, Temp, nx*ny*nz);
|
|
for (k=1; k<nz-1; k++){
|
|
for (j=1; j<ny-1; j++){
|
|
for (i=1; i<nz-1; i++){
|
|
//........................................................................
|
|
n = k*nx*ny+j*nx+i;
|
|
//........................................................................
|
|
iglobal = iproc*(nx-2)+i;
|
|
jglobal = jproc*(ny-2)+j;
|
|
kglobal = kproc*(nz-2)+k;
|
|
//........................................................................
|
|
Averages.dPdt(iglobal,jglobal,kglobal) = 0.0;
|
|
Averages.Phase_tplus(iglobal,jglobal,kglobal) = 0.0;
|
|
Averages.Phase_tminus(iglobal,jglobal,kglobal) = 0.0;
|
|
//........................................................................
|
|
}
|
|
}
|
|
}
|
|
|
|
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
|
|
// printf("Reading file %s \n",LocalRankFilename);
|
|
// printf("Sub-domain size: %i x %i x %i \n", nx,ny,nz);
|
|
ReadBinaryFile(LocalRankFilename, Temp, nx*ny*nz);
|
|
for (k=1; k<nz-1; k++){
|
|
for (j=1; j<ny-1; j++){
|
|
for (i=1; i<nz-1; i++){
|
|
|
|
//........................................................................
|
|
n = k*nx*ny+j*nx+i;
|
|
//........................................................................
|
|
iglobal = iproc*(nx-2)+i;
|
|
jglobal = jproc*(ny-2)+j;
|
|
kglobal = kproc*(nz-2)+k;
|
|
//........................................................................
|
|
Averages.SDs(iglobal,jglobal,kglobal) = Temp[n];
|
|
//........................................................................
|
|
}
|
|
}
|
|
}
|
|
|
|
sprintf(LocalRankFilename,"%s%s","Restart.",LocalRankString);
|
|
|
|
ReadFromRank(LocalRankFilename,Averages.Phase,Averages.Press,
|
|
Averages.Vel_x,Averages.Vel_y,Averages.Vel_z,
|
|
nx,ny,nz,iproc,jproc,kproc);
|
|
|
|
/* sprintf(LocalRankFilename,"%s%s","Pressure.",LocalRankString);
|
|
|
|
ReadBinaryFile(LocalRankFilename, Temp, nx*ny*nz);
|
|
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;
|
|
//........................................................................
|
|
iglobal = iproc*(nx-2)+i;
|
|
jglobal = jproc*(ny-2)+j;
|
|
kglobal = kproc*(nz-2)+k;
|
|
//........................................................................
|
|
Averages.Press(iglobal,jglobal,kglobal) = Temp[n];
|
|
//........................................................................
|
|
}
|
|
}
|
|
}
|
|
*/
|
|
//sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
|
|
//ReadBinaryFile(LocalRankFilename, Temp, nx*ny*nz);
|
|
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;
|
|
//........................................................................
|
|
iglobal = iproc*(nx-2)+i;
|
|
jglobal = jproc*(ny-2)+j;
|
|
kglobal = kproc*(nz-2)+k;
|
|
//........................................................................
|
|
//Averages.Phase(iglobal,jglobal,kglobal) = Temp[n];
|
|
Averages.SDn(iglobal,jglobal,kglobal) = Averages.Phase(iglobal,jglobal,kglobal);//Temp[n];
|
|
//........................................................................
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
}
|
|
}
|
|
}
|
|
printf("Read %i ranks of %s \n",nprocs,BaseFilename);
|
|
|
|
delete Temp;
|
|
|
|
// Initializing the blob ID
|
|
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 (Averages.SDs(i,j,k) < 0.0){
|
|
// Solid phase
|
|
Dm.id[n] = 0;
|
|
//PhaseLabel(i,j,k) = 0;
|
|
//WP(i,j,k) = -2;
|
|
//NWP(i,j,k) = -2;
|
|
}
|
|
else if (Averages.Phase(i,j,k) < 0.0){
|
|
// wetting phase
|
|
Dm.id[n] = 2;
|
|
//PhaseLabel(i,j,k) = 2;
|
|
//WP(i,j,k) = -2;
|
|
//NWP(i,j,k) = -1;
|
|
}
|
|
else {
|
|
// non-wetting phase
|
|
Dm.id[n] = 1;
|
|
//PhaseLabel(i,j,k) = 1;
|
|
//WP(i,j,k) = -1;
|
|
//NWP(i,j,k) = -2;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Compute the porosity
|
|
double porosity=0.0;
|
|
for (k=0; k<Nz; k++){
|
|
for (j=0; j<Ny; j++){
|
|
for (i=0; i<Nx; i++){
|
|
if (Dm.id[k*Nx*Ny+j*Nx+i] != 0){
|
|
porosity += 1.0;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
porosity /= (Nx*Ny*Nz*1.0);
|
|
printf("Media porosity is %f \n",porosity);
|
|
Dm.CommInit();
|
|
|
|
/* ****************************************************************
|
|
IDENTIFY ALL COMPONENTS FOR BOTH PHASES
|
|
****************************************************************** */
|
|
// int number_NWP_components = ComputeLocalPhaseComponent(PhaseLabel,1,NWP,true);
|
|
//int number_WP_components = ComputeLocalPhaseComponent(PhaseLabel,2,WP,true);
|
|
|
|
//printf("Number of WP components = %i \n",number_WP_components);
|
|
//printf("Number of NWP components = %i \n",number_NWP_components);
|
|
|
|
// Map the signed distance for the analysis
|
|
|
|
// Compute the porosity
|
|
porosity=0.0;
|
|
for (k=0; k<Nz; k++){
|
|
for (j=0; j<Ny; j++){
|
|
for (i=0; i<Nx; i++){
|
|
if (Averages.SDs(i,j,k) > 0.0){
|
|
porosity += 1.0;
|
|
}
|
|
Averages.SDs(i,j,k) -= (1.0);
|
|
}
|
|
}
|
|
}
|
|
porosity /= (Nx*Ny*Nz*1.0);
|
|
//printf("Media porosity is %f \n",porosity);
|
|
|
|
double beta=0.95;
|
|
int timestep=5;
|
|
Averages.Initialize();
|
|
Averages.ComputeDelPhi();
|
|
Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn);
|
|
Averages.UpdateMeshValues();
|
|
Averages.ComponentAverages();
|
|
Averages.SortBlobs();
|
|
Averages.PrintComponents(timestep);
|
|
|
|
// 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> SignDistVar( new IO::Variable() );
|
|
std::shared_ptr<IO::Variable> LabelWPVar( new IO::Variable() );
|
|
std::shared_ptr<IO::Variable> LabelNWPVar( new IO::Variable() );
|
|
std::shared_ptr<IO::Variable> PhaseIDVar( new IO::Variable() );
|
|
|
|
PhaseVar->name = "phase";
|
|
PhaseVar->type = IO::VariableType::VolumeVariable;
|
|
PhaseVar->dim = 1;
|
|
PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
|
|
meshData[0].vars.push_back(PhaseVar);
|
|
|
|
SignDistVar->name = "SignDist";
|
|
SignDistVar->type = IO::VariableType::VolumeVariable;
|
|
SignDistVar->dim = 1;
|
|
SignDistVar->data.resize(Nx-2,Ny-2,Nz-2);
|
|
meshData[0].vars.push_back(SignDistVar);
|
|
|
|
LabelNWPVar->name = "LabelNWP";
|
|
LabelNWPVar->type = IO::VariableType::VolumeVariable;
|
|
LabelNWPVar->dim = 1;
|
|
LabelNWPVar->data.resize(Nx-2,Ny-2,Nz-2);
|
|
meshData[0].vars.push_back(LabelNWPVar);
|
|
|
|
LabelWPVar->name = "LabelWP";
|
|
LabelWPVar->type = IO::VariableType::VolumeVariable;
|
|
LabelWPVar->dim = 1;
|
|
LabelWPVar->data.resize(Nx-2,Ny-2,Nz-2);
|
|
meshData[0].vars.push_back(LabelWPVar);
|
|
|
|
PhaseIDVar->name = "PhaseID";
|
|
PhaseIDVar->type = IO::VariableType::VolumeVariable;
|
|
PhaseIDVar->dim = 1;
|
|
PhaseIDVar->data.resize(Nx-2,Ny-2,Nz-2);
|
|
meshData[0].vars.push_back(PhaseIDVar);
|
|
|
|
fillData.copy(Averages.SDn,PhaseVar->data);
|
|
fillData.copy(Averages.SDs,SignDistVar->data);
|
|
fillData.copy(Averages.Label_WP,LabelWPVar->data);
|
|
fillData.copy(Averages.Label_NWP,LabelNWPVar->data);
|
|
fillData.copy(Averages.PhaseID,PhaseIDVar->data);
|
|
IO::writeData( 0, meshData, comm );
|
|
/*
|
|
FILE *NWP_FILE;
|
|
NWP_FILE = fopen("NWP.dat","wb");
|
|
fwrite(Averages.Label_NWP.get(),4,Nx*Ny*Nz,NWP_FILE);
|
|
fclose(NWP_FILE);
|
|
|
|
FILE *WP_FILE;
|
|
WP_FILE = fopen("WP.dat","wb");
|
|
fwrite(Averages.Label_WP.get(),4,Nx*Ny*Nz,WP_FILE);
|
|
fclose(WP_FILE);
|
|
|
|
FILE *DISTANCE;
|
|
DISTANCE = fopen("SignDist.dat","wb");
|
|
fwrite(Averages.SDs.get(),8,Nx*Ny*Nz,DISTANCE);
|
|
fclose(DISTANCE);
|
|
*/
|
|
// ****************************************************
|
|
comm.barrier();
|
|
Utilities::shutdown();
|
|
// ****************************************************
|
|
}
|
|
|