fix the inlet cross-sectional area calculations

This commit is contained in:
Rex Zhe Li
2017-12-20 18:17:00 +11:00
parent a9dda36754
commit 0b2f116fb0

View File

@@ -370,7 +370,6 @@ int main(int argc, char **argv)
// char value;
char *id;
id = new char[N];
int sum = 0;
double sum_local;
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
if (BoundaryCondition > 0) iVol_global = 1.0/(1.0*(Nx-2)*nprocx*(Ny-2)*nprocy*((Nz-2)*nprocz-6));
@@ -387,41 +386,6 @@ int main(int argc, char **argv)
ReadBinaryFile(LocalRankFilename, Averages->SDs.data(), N);
MPI_Barrier(comm);
if (rank == 0) cout << "Domain set." << endl;
//.......................................................................
// Assign the phase ID field based on the signed distance
//.......................................................................
double inlet_area_local=0.f;
double InletArea;
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
id[n] = 0;
}
}
}
sum=0;
pore_vol = 0.0;
for ( k=0;k<Nz;k++){
for ( j=0;j<Ny;j++){
for ( i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
if (Averages->SDs(n) > 0.0){
id[n] = 2;
}
// compute the porosity (actual interface location used)
if (Averages->SDs(n) > 0.0){
sum++;
if (k==1){
inlet_area_local+=1.f;
}
}
}
}
}
MPI_Allreduce(&inlet_area_local,&InletArea,1,MPI_DOUBLE,MPI_SUM,comm);
if (rank==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
sprintf(LocalRankFilename,"ID.%05i",rank);
@@ -430,9 +394,27 @@ int main(int argc, char **argv)
if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx");
readID=fread(id,1,N,IDFILE);
if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID (rank=%i) \n",rank);
fclose(IDFILE);
// Calculate the actual inlet cross-sectional area
double inlet_area_local=0.f;
double InletArea;
if (Dm.kproc==0)
{
for ( k=1;k<2;k++){
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
int n = k*Nx*Ny+j*Nx+i;
if (id[n] > 0){
inlet_area_local+=1.f;
}
}
}
}
}
MPI_Allreduce(&inlet_area_local,&InletArea,1,MPI_DOUBLE,MPI_SUM,comm);
// Set up kstart, kfinish so that the reservoirs are excluded from averaging
int kstart,kfinish;
kstart = 1;