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

merge vis docs

Browse Source
This commit is contained in:
James McClure 2021-11-03 09:18:25 -04:00
commit 898107e650
23 changed files with 931 additions and 423 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

732
analysis/morphology.cpp
View File

@ -24,6 +24,409 @@ inline void UnpackID(const int *list, int count, signed char *recvbuf, signed ch
}
}
Morphology::Morphology(){
/* MPI tags*/
sendtag = recvtag = 1381;
}
Morphology::~Morphology(){
}
void Morphology::Initialize(std::shared_ptr <Domain> Dm, DoubleArray &Distance){
/* Loop over all faces and determine overlaps */
size_t Nx = Dm->Nx;
size_t Ny = Dm->Ny;
size_t Nz = Dm->Nz;
size_t N = Nx*Ny*Nz;
int *tmpShift_x, *tmpShift_y, *tmpShift_z;
double *tmpDistance;
tmpShift_x = new int [N];
tmpShift_y= new int [N];
tmpShift_z = new int [N];
tmpDistance = new double [N];
double distance, boundary_distance;
/* Loop over the local sub-domain and create overlap lists for each neighboring sub-domain */
int sendLoc = 0; // counter for the local sub-domain send values
int recvLoc = 0; // counter for the local recv
//...................................................
/* x face */
sendCount = recvCount = 0;
for (size_t k=1; k<Nz-1; k++){
for (size_t j=1; j<Ny-1; j++){
for (size_t i=1; i<Nx-1; i++){
distance = Distance(i,j,k);
// Distance to x boundary
boundary_distance = double(i-1);
if (distance > boundary_distance){
tmpShift_x[sendCount] = Nx + i-2;
tmpShift_y[sendCount] = j;
tmpShift_z[sendCount] = k;
tmpDistance[sendCount++] = distance;
int n = k*Nx*Ny + j*Nx + i;
sendID.push_back(n);
}
}
}
}
Dm->Comm.Irecv(&recvCount,1,Dm->rank_X(),recvtag+0);
Dm->Comm.send(&sendCount,1,Dm->rank_x(),sendtag+0);
Dm->Comm.barrier();
sendOffset_x = sendLoc;
recvOffset_X = recvLoc;
sendLoc += sendCount;
recvLoc += recvCount;
sendCount_x = sendCount;
recvCount_X = recvCount;
/* grow the arrays */
xShift.resize(recvLoc);
yShift.resize(recvLoc);
zShift.resize(recvLoc);
morphRadius.resize(recvLoc);
//..............................
/* send the morphological radius */
Dm->Comm.Irecv(&morphRadius[recvOffset_X],recvCount,Dm->rank_X(),recvtag+0);
Dm->Comm.send(&tmpDistance[0],sendCount,Dm->rank_x(),sendtag+0);
/* send the shift values */
Dm->Comm.Irecv(&xShift[recvOffset_X],recvCount,Dm->rank_X(),recvtag+1);
Dm->Comm.send(&tmpShift_x[0],sendCount,Dm->rank_x(),sendtag+1);
Dm->Comm.Irecv(&yShift[recvOffset_X],recvCount,Dm->rank_X(),recvtag+2);
Dm->Comm.send(&tmpShift_y[0],sendCount,Dm->rank_x(),sendtag+2);
Dm->Comm.Irecv(&zShift[recvOffset_X],recvCount,Dm->rank_X(),recvtag+3);
Dm->Comm.send(&tmpShift_z[0],sendCount,Dm->rank_x(),sendtag+3);
Dm->Comm.barrier();
//...................................................
//...................................................
/* X face */
sendCount = recvCount = 0;
for (size_t k=1; k<Nz-1; k++){
for (size_t j=1; j<Ny-1; j++){
for (size_t i=1; i<Nx-1; i++){
distance = Distance(i,j,k);
// Distance to x boundary
boundary_distance = double(Nx-i-1);
if (distance > boundary_distance){
tmpShift_x[sendCount] = (i)-(Nx-2);
tmpShift_y[sendCount] = j;
tmpShift_z[sendCount] = k;
tmpDistance[sendCount++] = distance;
int n = k*Nx*Ny + j*Nx + i;
sendID.push_back(n);
}
}
}
}
Dm->Comm.Irecv(&recvCount,1,Dm->rank_x(),recvtag+0);
Dm->Comm.send(&sendCount,1,Dm->rank_X(),sendtag+0);
Dm->Comm.barrier();
sendOffset_X = sendLoc;
recvOffset_x = recvLoc;
sendLoc += sendCount;
recvLoc += recvCount;
sendCount_X = sendCount;
recvCount_x = recvCount;
/* grow the arrays */
xShift.resize(recvLoc);
yShift.resize(recvLoc);
zShift.resize(recvLoc);
morphRadius.resize(recvLoc);
//..............................
/* send the morphological radius */
Dm->Comm.Irecv(&morphRadius[recvOffset_x],recvCount,Dm->rank_x(),recvtag+0);
Dm->Comm.send(&tmpDistance[0],sendCount,Dm->rank_X(),sendtag+0);
/* send the shift values */
Dm->Comm.Irecv(&xShift[recvOffset_x],recvCount,Dm->rank_x(),recvtag+1);
Dm->Comm.send(&tmpShift_x[0],sendCount,Dm->rank_X(),sendtag+1);
Dm->Comm.Irecv(&yShift[recvOffset_x],recvCount,Dm->rank_x(),recvtag+2);
Dm->Comm.send(&tmpShift_y[0],sendCount,Dm->rank_X(),sendtag+2);
Dm->Comm.Irecv(&zShift[recvOffset_x],recvCount,Dm->rank_x(),recvtag+3);
Dm->Comm.send(&tmpShift_z[0],sendCount,Dm->rank_X(),sendtag+3);
Dm->Comm.barrier();
//...................................................
//...................................................
/* y face */
sendCount = recvCount = 0;
for (size_t k=1; k<Nz-1; k++){
for (size_t j=1; j<Ny-1; j++){
for (size_t i=1; i<Nx-1; i++){
distance = Distance(i,j,k);
// Distance to y boundary
boundary_distance = double(j-1);
if (distance > boundary_distance){
tmpShift_x[sendCount] = i;
tmpShift_y[sendCount] = Ny + j-2;
tmpShift_z[sendCount] = k;
tmpDistance[sendCount++] = distance;
int n = k*Nx*Ny + j*Nx + i;
sendID.push_back(n);
}
}
}
}
Dm->Comm.Irecv(&recvCount,1,Dm->rank_Y(),recvtag+0);
Dm->Comm.send(&sendCount,1,Dm->rank_y(),sendtag+0);
Dm->Comm.barrier();
sendOffset_y = sendLoc;
recvOffset_Y = recvLoc;
sendLoc += sendCount;
recvLoc += recvCount;
sendCount_y = sendCount;
recvCount_Y = recvCount;
/* grow the arrays */
xShift.resize(recvLoc);
yShift.resize(recvLoc);
zShift.resize(recvLoc);
morphRadius.resize(recvLoc);
//..............................
/* send the morphological radius */
Dm->Comm.Irecv(&morphRadius[recvOffset_Y],recvCount,Dm->rank_Y(),recvtag+0);
Dm->Comm.send(&tmpDistance[0],sendCount,Dm->rank_y(),sendtag+0);
/* send the shift values */
Dm->Comm.Irecv(&xShift[recvOffset_Y],recvCount,Dm->rank_Y(),recvtag+1);
Dm->Comm.send(&tmpShift_x[0],sendCount,Dm->rank_y(),sendtag+1);
Dm->Comm.Irecv(&yShift[recvOffset_Y],recvCount,Dm->rank_Y(),recvtag+2);
Dm->Comm.send(&tmpShift_y[0],sendCount,Dm->rank_y(),sendtag+2);
Dm->Comm.Irecv(&zShift[recvOffset_Y],recvCount,Dm->rank_Y(),recvtag+3);
Dm->Comm.send(&tmpShift_z[0],sendCount,Dm->rank_y(),sendtag+3);
Dm->Comm.barrier();
//...................................................
//...................................................
/* X face */
sendCount = recvCount = 0;
for (size_t k=1; k<Nz-1; k++){
for (size_t j=1; j<Ny-1; j++){
for (size_t i=1; i<Nx-1; i++){
distance = Distance(i,j,k);
// Distance to x boundary
boundary_distance = double(Ny-j-1);
if (distance > boundary_distance){
tmpShift_x[sendCount] = i;
tmpShift_y[sendCount] = j-(Ny-2);
tmpShift_z[sendCount] = k;
tmpDistance[sendCount++] = distance;
int n = k*Nx*Ny + j*Nx + i;
sendID.push_back(n);
}
}
}
}
Dm->Comm.Irecv(&recvCount,1,Dm->rank_y(),recvtag+0);
Dm->Comm.send(&sendCount,1,Dm->rank_Y(),sendtag+0);
Dm->Comm.barrier();
sendOffset_Y = sendLoc;
recvOffset_y = recvLoc;
sendLoc += sendCount;
recvLoc += recvCount;
sendCount_Y = sendCount;
recvCount_y = recvCount;
/* grow the arrays */
xShift.resize(recvLoc);
yShift.resize(recvLoc);
zShift.resize(recvLoc);
morphRadius.resize(recvLoc);
//..............................
/* send the morphological radius */
Dm->Comm.Irecv(&morphRadius[recvOffset_y],recvCount,Dm->rank_y(),recvtag+0);
Dm->Comm.send(&tmpDistance[0],sendCount,Dm->rank_Y(),sendtag+0);
/* send the shift values */
Dm->Comm.Irecv(&xShift[recvOffset_y],recvCount,Dm->rank_y(),recvtag+1);
Dm->Comm.send(&tmpShift_x[0],sendCount,Dm->rank_Y(),sendtag+1);
Dm->Comm.Irecv(&yShift[recvOffset_y],recvCount,Dm->rank_y(),recvtag+2);
Dm->Comm.send(&tmpShift_y[0],sendCount,Dm->rank_Y(),sendtag+2);
Dm->Comm.Irecv(&zShift[recvOffset_y],recvCount,Dm->rank_y(),recvtag+3);
Dm->Comm.send(&tmpShift_z[0],sendCount,Dm->rank_Y(),sendtag+3);
Dm->Comm.barrier();
//...................................................
//...................................................
/* z face */
sendCount = recvCount = 0;
for (size_t k=1; k<Nz-1; k++){
for (size_t j=1; j<Ny-1; j++){
for (size_t i=1; i<Nx-1; i++){
distance = Distance(i,j,k);
// Distance to z boundary
boundary_distance = double(k-1);
if (distance > boundary_distance){
tmpShift_x[sendCount] = i;
tmpShift_y[sendCount] = j;
tmpShift_z[sendCount] = (Nz-2) + k;
tmpDistance[sendCount++] = distance;
int n = k*Nx*Ny + j*Nx + i;
sendID.push_back(n);
}
}
}
}
Dm->Comm.Irecv(&recvCount,1,Dm->rank_Z(),recvtag+0);
Dm->Comm.send(&sendCount,1,Dm->rank_z(),sendtag+0);
Dm->Comm.barrier();
sendOffset_z = sendLoc;
recvOffset_Z = recvLoc;
sendLoc += sendCount;
recvLoc += recvCount;
sendCount_z = sendCount;
recvCount_Z = recvCount;
/* grow the arrays */
xShift.resize(recvLoc);
yShift.resize(recvLoc);
zShift.resize(recvLoc);
morphRadius.resize(recvLoc);
//..............................
/* send the morphological radius */
Dm->Comm.Irecv(&morphRadius[recvOffset_Z],recvCount,Dm->rank_Z(),recvtag+0);
Dm->Comm.send(&tmpDistance[0],sendCount,Dm->rank_z(),sendtag+0);
/* send the shift values */
Dm->Comm.Irecv(&xShift[recvOffset_Z],recvCount,Dm->rank_Z(),recvtag+1);
Dm->Comm.send(&tmpShift_x[0],sendCount,Dm->rank_z(),sendtag+1);
Dm->Comm.Irecv(&yShift[recvOffset_Z],recvCount,Dm->rank_Z(),recvtag+2);
Dm->Comm.send(&tmpShift_y[0],sendCount,Dm->rank_z(),sendtag+2);
Dm->Comm.Irecv(&zShift[recvOffset_Z],recvCount,Dm->rank_Z(),recvtag+3);
Dm->Comm.send(&tmpShift_z[0],sendCount,Dm->rank_z(),sendtag+3);
Dm->Comm.barrier();
//...................................................
/* Z face */
sendCount = recvCount = 0;
for (size_t k=1; k<Nz-1; k++){
for (size_t j=1; j<Ny-1; j++){
for (size_t i=1; i<Nx-1; i++){
distance = Distance(i,j,k);
// Distance to x boundary
boundary_distance = double(Nz-k-1);
if (distance > boundary_distance){
tmpShift_x[sendCount] = i;
tmpShift_y[sendCount] = j;
tmpShift_z[sendCount] = k-(Nz-2);
tmpDistance[sendCount++] = distance;
int n = k*Nx*Ny + j*Nx + i;
sendID.push_back(n);
}
}
}
}
Dm->Comm.Irecv(&recvCount,1,Dm->rank_z(),recvtag+0);
Dm->Comm.send(&sendCount,1,Dm->rank_Z(),sendtag+0);
Dm->Comm.barrier();
sendOffset_Z = sendLoc;
recvOffset_z = recvLoc;
sendLoc += sendCount;
recvLoc += recvCount;
sendCount_Z = sendCount;
recvCount_z = recvCount;
/* grow the arrays */
xShift.resize(recvLoc);
yShift.resize(recvLoc);
zShift.resize(recvLoc);
morphRadius.resize(recvLoc);
//..............................
/* send the morphological radius */
Dm->Comm.Irecv(&morphRadius[recvOffset_z],recvCount,Dm->rank_z(),recvtag+0);
Dm->Comm.send(&tmpDistance[0],sendCount,Dm->rank_Z(),sendtag+0);
/* send the shift values */
Dm->Comm.Irecv(&xShift[recvOffset_z],recvCount,Dm->rank_z(),recvtag+1);
Dm->Comm.send(&tmpShift_x[0],sendCount,Dm->rank_Z(),sendtag+1);
Dm->Comm.Irecv(&yShift[recvOffset_z],recvCount,Dm->rank_z(),recvtag+2);
Dm->Comm.send(&tmpShift_y[0],sendCount,Dm->rank_Z(),sendtag+2);
Dm->Comm.Irecv(&zShift[recvOffset_z],recvCount,Dm->rank_z(),recvtag+3);
Dm->Comm.send(&tmpShift_z[0],sendCount,Dm->rank_Z(),sendtag+3);
Dm->Comm.barrier();
//...................................................
/* resize the send / recv lists */
sendCount = sendLoc;
recvCount = recvLoc;
sendList.resize(sendLoc);
recvList.resize(recvLoc);
localID.resize(sendCount);
nonlocalID.resize(recvCount);
/*printf(" offset %i for send (x) %i \n", sendOffset_x, sendCount_x);
printf(" offset %i for send (X) %i \n", sendOffset_X, sendCount_X);
printf(" offset %i for send (y) %i \n", sendOffset_y, sendCount_y);
printf(" offset %i for send (Y) %i \n", sendOffset_Y, sendCount_Y);
printf(" offset %i for send (z) %i \n", sendOffset_z, sendCount_z);
printf(" offset %i for send (Z) %i \n", sendOffset_Z, sendCount_Z);
*/
}
int Morphology::GetOverlaps(std::shared_ptr <Domain> Dm, signed char *id, const signed char ErodeLabel, const signed char NewLabel){
int Nx = Dm->Nx;
int Ny = Dm->Ny;
int Nz = Dm->Nz;
int LocalNumber=0;
int i,j,k,ii,jj,kk;
int imin,jmin,kmin,imax,jmax,kmax;
for (int idx=0; idx<sendCount; idx++){
int n = sendID[idx];
localID[idx] = id[n];
}
//printf("send x -- offset: %i, count: %i \n",sendOffset_x,sendCount_x);
Dm->Comm.Irecv(&nonlocalID[recvOffset_X],recvCount_X,Dm->rank_x(),recvtag+2);
Dm->Comm.send(&localID[sendOffset_x],sendCount_x,Dm->rank_X(),sendtag+2);
//printf("send X \n");
Dm->Comm.Irecv(&nonlocalID[recvOffset_x],recvCount_x,Dm->rank_X(),recvtag+3);
Dm->Comm.send(&localID[sendOffset_X],sendCount_X,Dm->rank_x(),sendtag+3);
//printf("send y \n");
Dm->Comm.Irecv(&nonlocalID[recvOffset_Y],recvCount_Y,Dm->rank_y(),recvtag+4);
Dm->Comm.send(&localID[sendOffset_y],sendCount_y,Dm->rank_Y(),sendtag+4);
//printf("send Y \n");
Dm->Comm.Irecv(&nonlocalID[recvOffset_y],recvCount_y,Dm->rank_Y(),recvtag+5);
Dm->Comm.send(&localID[sendOffset_Y],sendCount_Y,Dm->rank_y(),sendtag+5);
//printf("send z \n");
Dm->Comm.Irecv(&nonlocalID[recvOffset_Z],recvCount_Z,Dm->rank_z(),recvtag+6);
Dm->Comm.send(&localID[sendOffset_z],sendCount_z,Dm->rank_Z(),sendtag+6);
//printf("send Z \n");
Dm->Comm.Irecv(&nonlocalID[recvOffset_z],recvCount_z,Dm->rank_Z(),recvtag+7);
Dm->Comm.send(&localID[sendOffset_Z],sendCount_Z,Dm->rank_z(),sendtag+7);
for (int idx=0; idx<recvCount; idx++){
double radius = morphRadius[idx];
signed char label = nonlocalID[idx];
/* get the neighboring site index */
i = xShift[idx];
j = yShift[idx];
k = zShift[idx];
int Window = int(radius);
// loop over the window and update
if (label == NewLabel){
imin=max(1,i-Window);
jmin=max(1,j-Window);
kmin=max(1,k-Window);
imax=min(Nx-1,i+Window);
jmax=min(Ny-1,j+Window);
kmax=min(Nz-1,k+Window);
for (kk=kmin; kk<kmax; kk++){
for (jj=jmin; jj<jmax; jj++){
for (ii=imin; ii<imax; ii++){
int nn = kk*Nx*Ny+jj*Nx+ii;
double dsq = double((ii-i)*(ii-i)+(jj-j)*(jj-j)+(kk-k)*(kk-k));
if (id[nn] == ErodeLabel && dsq <= radius*radius){
LocalNumber+=1.0;
id[nn]=NewLabel;
}
}
}
}
}
}
Dm->Comm.barrier();
return LocalNumber;
}
//***************************************************************************************
double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain> Dm, double VoidFraction, signed char ErodeLabel, signed char NewLabel){
// SignDist is the distance to the object that you want to constaing the morphological opening
@ -60,6 +463,9 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
}
Dm->Comm.barrier();
Morphology Structure;
Structure.Initialize(Dm,SignDist);
// total Global is the number of nodes in the pore-space
totalGlobal = Dm->Comm.sumReduce( count );
maxdistGlobal = Dm->Comm.sumReduce( maxdist );
@ -69,56 +475,6 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
if (rank==0) printf("Maximum pore size: %f \n",maxdistGlobal);
final_void_fraction = volume_fraction; //initialize
// Communication buffers
signed char *sendID_x, *sendID_y, *sendID_z, *sendID_X, *sendID_Y, *sendID_Z;
signed char *sendID_xy, *sendID_yz, *sendID_xz, *sendID_Xy, *sendID_Yz, *sendID_xZ;
signed char *sendID_xY, *sendID_yZ, *sendID_Xz, *sendID_XY, *sendID_YZ, *sendID_XZ;
signed char *recvID_x, *recvID_y, *recvID_z, *recvID_X, *recvID_Y, *recvID_Z;
signed char *recvID_xy, *recvID_yz, *recvID_xz, *recvID_Xy, *recvID_Yz, *recvID_xZ;
signed char *recvID_xY, *recvID_yZ, *recvID_Xz, *recvID_XY, *recvID_YZ, *recvID_XZ;
// send buffers
sendID_x = new signed char [Dm->sendCount("x")];
sendID_y = new signed char [Dm->sendCount("y")];
sendID_z = new signed char [Dm->sendCount("z")];
sendID_X = new signed char [Dm->sendCount("X")];
sendID_Y = new signed char [Dm->sendCount("Y")];
sendID_Z = new signed char [Dm->sendCount("Z")];
sendID_xy = new signed char [Dm->sendCount("xy")];
sendID_yz = new signed char [Dm->sendCount("yz")];
sendID_xz = new signed char [Dm->sendCount("xz")];
sendID_Xy = new signed char [Dm->sendCount("Xy")];
sendID_Yz = new signed char [Dm->sendCount("Yz")];
sendID_xZ = new signed char [Dm->sendCount("xZ")];
sendID_xY = new signed char [Dm->sendCount("xY")];
sendID_yZ = new signed char [Dm->sendCount("yZ")];
sendID_Xz = new signed char [Dm->sendCount("Xz")];
sendID_XY = new signed char [Dm->sendCount("XY")];
sendID_YZ = new signed char [Dm->sendCount("YZ")];
sendID_XZ = new signed char [Dm->sendCount("XZ")];
//......................................................................................
// recv buffers
recvID_x = new signed char [Dm->recvCount("x")];
recvID_y = new signed char [Dm->recvCount("y")];
recvID_z = new signed char [Dm->recvCount("z")];
recvID_X = new signed char [Dm->recvCount("X")];
recvID_Y = new signed char [Dm->recvCount("Y")];
recvID_Z = new signed char [Dm->recvCount("Z")];
recvID_xy = new signed char [Dm->recvCount("xy")];
recvID_yz = new signed char [Dm->recvCount("yz")];
recvID_xz = new signed char [Dm->recvCount("xz")];
recvID_Xy = new signed char [Dm->recvCount("Xy")];
recvID_xZ = new signed char [Dm->recvCount("xZ")];
recvID_xY = new signed char [Dm->recvCount("xY")];
recvID_yZ = new signed char [Dm->recvCount("yZ")];
recvID_Yz = new signed char [Dm->recvCount("Yz")];
recvID_Xz = new signed char [Dm->recvCount("Xz")];
recvID_XY = new signed char [Dm->recvCount("XY")];
recvID_YZ = new signed char [Dm->recvCount("YZ")];
recvID_XZ = new signed char [Dm->recvCount("XZ")];
//......................................................................................
int sendtag,recvtag;
sendtag = recvtag = 7;
int ii,jj,kk;
int imin,jmin,kmin,imax,jmax,kmax;
int Nx = nx;
@ -138,12 +494,18 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
double Rcrit_old = maxdistGlobal;
double Rcrit_new = maxdistGlobal;
while (void_fraction_new > VoidFraction)
int numTry = 0;
int maxTry = 100;
while (void_fraction_new > VoidFraction && numTry < maxTry)
{
numTry++;
void_fraction_diff_old = void_fraction_diff_new;
void_fraction_old = void_fraction_new;
Rcrit_old = Rcrit_new;
Rcrit_new -= deltaR*Rcrit_old;
if (Rcrit_new < 0.5 ){
numTry = maxTry;
}
int Window=round(Rcrit_new);
if (Window == 0) Window = 1; // If Window = 0 at the begining, after the following process will have sw=1.0
// and sw<Sw will be immediately broken
@ -178,66 +540,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
}
}
}
// Pack and send the updated ID values
PackID(Dm->sendList("x"), Dm->sendCount("x") ,sendID_x, id);
PackID(Dm->sendList("X"), Dm->sendCount("X") ,sendID_X, id);
PackID(Dm->sendList("y"), Dm->sendCount("y") ,sendID_y, id);
PackID(Dm->sendList("Y"), Dm->sendCount("Y") ,sendID_Y, id);
PackID(Dm->sendList("z"), Dm->sendCount("z") ,sendID_z, id);
PackID(Dm->sendList("Z"), Dm->sendCount("Z") ,sendID_Z, id);
PackID(Dm->sendList("xy"), Dm->sendCount("xy") ,sendID_xy, id);
PackID(Dm->sendList("Xy"), Dm->sendCount("Xy") ,sendID_Xy, id);
PackID(Dm->sendList("xY"), Dm->sendCount("xY") ,sendID_xY, id);
PackID(Dm->sendList("XY"), Dm->sendCount("XY") ,sendID_XY, id);
PackID(Dm->sendList("xz"), Dm->sendCount("xz") ,sendID_xz, id);
PackID(Dm->sendList("Xz"), Dm->sendCount("Xz") ,sendID_Xz, id);
PackID(Dm->sendList("xZ"), Dm->sendCount("xZ") ,sendID_xZ, id);
PackID(Dm->sendList("XZ"), Dm->sendCount("XZ") ,sendID_XZ, id);
PackID(Dm->sendList("yz"), Dm->sendCount("yz") ,sendID_yz, id);
PackID(Dm->sendList("Yz"), Dm->sendCount("Yz") ,sendID_Yz, id);
PackID(Dm->sendList("yZ"), Dm->sendCount("yZ") ,sendID_yZ, id);
PackID(Dm->sendList("YZ"), Dm->sendCount("YZ") ,sendID_YZ, id);
//......................................................................................
Dm->Comm.sendrecv(sendID_x,Dm->sendCount("x"),Dm->rank_x(),sendtag,recvID_X,Dm->recvCount("X"),Dm->rank_X(),recvtag);
Dm->Comm.sendrecv(sendID_X,Dm->sendCount("X"),Dm->rank_X(),sendtag,recvID_x,Dm->recvCount("x"),Dm->rank_x(),recvtag);
Dm->Comm.sendrecv(sendID_y,Dm->sendCount("y"),Dm->rank_y(),sendtag,recvID_Y,Dm->recvCount("Y"),Dm->rank_Y(),recvtag);
Dm->Comm.sendrecv(sendID_Y,Dm->sendCount("Y"),Dm->rank_Y(),sendtag,recvID_y,Dm->recvCount("y"),Dm->rank_y(),recvtag);
Dm->Comm.sendrecv(sendID_z,Dm->sendCount("z"),Dm->rank_z(),sendtag,recvID_Z,Dm->recvCount("Z"),Dm->rank_Z(),recvtag);
Dm->Comm.sendrecv(sendID_Z,Dm->sendCount("Z"),Dm->rank_Z(),sendtag,recvID_z,Dm->recvCount("z"),Dm->rank_z(),recvtag);
Dm->Comm.sendrecv(sendID_xy,Dm->sendCount("xy"),Dm->rank_xy(),sendtag,recvID_XY,Dm->recvCount("XY"),Dm->rank_XY(),recvtag);
Dm->Comm.sendrecv(sendID_XY,Dm->sendCount("XY"),Dm->rank_XY(),sendtag,recvID_xy,Dm->recvCount("xy"),Dm->rank_xy(),recvtag);
Dm->Comm.sendrecv(sendID_Xy,Dm->sendCount("Xy"),Dm->rank_Xy(),sendtag,recvID_xY,Dm->recvCount("xY"),Dm->rank_xY(),recvtag);
Dm->Comm.sendrecv(sendID_xY,Dm->sendCount("xY"),Dm->rank_xY(),sendtag,recvID_Xy,Dm->recvCount("Xy"),Dm->rank_Xy(),recvtag);
Dm->Comm.sendrecv(sendID_xz,Dm->sendCount("xz"),Dm->rank_xz(),sendtag,recvID_XZ,Dm->recvCount("XZ"),Dm->rank_XZ(),recvtag);
Dm->Comm.sendrecv(sendID_XZ,Dm->sendCount("XZ"),Dm->rank_XZ(),sendtag,recvID_xz,Dm->recvCount("xz"),Dm->rank_xz(),recvtag);
Dm->Comm.sendrecv(sendID_Xz,Dm->sendCount("Xz"),Dm->rank_Xz(),sendtag,recvID_xZ,Dm->recvCount("xZ"),Dm->rank_xZ(),recvtag);
Dm->Comm.sendrecv(sendID_xZ,Dm->sendCount("xZ"),Dm->rank_xZ(),sendtag,recvID_Xz,Dm->recvCount("Xz"),Dm->rank_Xz(),recvtag);
Dm->Comm.sendrecv(sendID_yz,Dm->sendCount("yz"),Dm->rank_yz(),sendtag,recvID_YZ,Dm->recvCount("YZ"),Dm->rank_YZ(),recvtag);
Dm->Comm.sendrecv(sendID_YZ,Dm->sendCount("YZ"),Dm->rank_YZ(),sendtag,recvID_yz,Dm->recvCount("yz"),Dm->rank_yz(),recvtag);
Dm->Comm.sendrecv(sendID_Yz,Dm->sendCount("Yz"),Dm->rank_Yz(),sendtag,recvID_yZ,Dm->recvCount("yZ"),Dm->rank_yZ(),recvtag);
Dm->Comm.sendrecv(sendID_yZ,Dm->sendCount("yZ"),Dm->rank_yZ(),sendtag,recvID_Yz,Dm->recvCount("Yz"),Dm->rank_Yz(),recvtag);
//......................................................................................
UnpackID(Dm->recvList("x"), Dm->recvCount("x") ,recvID_x, id);
UnpackID(Dm->recvList("X"), Dm->recvCount("X") ,recvID_X, id);
UnpackID(Dm->recvList("y"), Dm->recvCount("y") ,recvID_y, id);
UnpackID(Dm->recvList("Y"), Dm->recvCount("Y") ,recvID_Y, id);
UnpackID(Dm->recvList("z"), Dm->recvCount("z") ,recvID_z, id);
UnpackID(Dm->recvList("Z"), Dm->recvCount("Z") ,recvID_Z, id);
UnpackID(Dm->recvList("xy"), Dm->recvCount("xy") ,recvID_xy, id);
UnpackID(Dm->recvList("Xy"), Dm->recvCount("Xy") ,recvID_Xy, id);
UnpackID(Dm->recvList("xY"), Dm->recvCount("xY") ,recvID_xY, id);
UnpackID(Dm->recvList("XY"), Dm->recvCount("XY") ,recvID_XY, id);
UnpackID(Dm->recvList("xz"), Dm->recvCount("xz") ,recvID_xz, id);
UnpackID(Dm->recvList("Xz"), Dm->recvCount("Xz") ,recvID_Xz, id);
UnpackID(Dm->recvList("xZ"), Dm->recvCount("xZ") ,recvID_xZ, id);
UnpackID(Dm->recvList("XZ"), Dm->recvCount("XZ") ,recvID_XZ, id);
UnpackID(Dm->recvList("yz"), Dm->recvCount("yz") ,recvID_yz, id);
UnpackID(Dm->recvList("Yz"), Dm->recvCount("Yz") ,recvID_Yz, id);
UnpackID(Dm->recvList("yZ"), Dm->recvCount("yZ") ,recvID_yZ, id);
UnpackID(Dm->recvList("YZ"), Dm->recvCount("YZ") ,recvID_YZ, id);
//......................................................................................
//double GlobalNumber = Dm->Comm.sumReduce( LocalNumber );
LocalNumber += Structure.GetOverlaps(Dm,id,ErodeLabel,NewLabel);
count = 0.f;
for (int k=1; k<Nz-1; k++){
@ -253,10 +556,10 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
countGlobal = Dm->Comm.sumReduce( count );
void_fraction_new = countGlobal/totalGlobal;
void_fraction_diff_new = abs(void_fraction_new-VoidFraction);
/* if (rank==0){
if (rank==0){
printf(" %f ",void_fraction_new);
printf(" %f\n",Rcrit_new);
} */
}
}
if (void_fraction_diff_new<void_fraction_diff_old){
@ -275,30 +578,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
}
return final_void_fraction;
}
/*
double morph_open()
{
fillHalo<char> fillChar(Dm->Comm,Dm->rank_info,{Nx-2,Ny-2,Nz-2},{1,1,1},0,1);
GlobalNumber = Dm->Comm.sumReduce( LocalNumber );
count = 0.f;
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] == 2){
count+=1.0;
}
}
}
}
countGlobal = Dm->Comm.sumReduce( count );
return countGlobal;
}
*/
//***************************************************************************************
double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain> Dm, double VoidFraction){
@ -307,6 +587,9 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
// id is a labeled map
// Dm contains information about the domain structure
signed char ErodeLabel = 2;
signed char NewLabel = 1;
int nx = Dm->Nx;
int ny = Dm->Ny;
int nz = Dm->Nz;
@ -319,6 +602,9 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
IntArray phase_label(nx,ny,nz);
Array<char> ID(nx,ny,nz);
fillHalo<char> fillChar(Dm->Comm,Dm->rank_info,{nx-2,ny-2,nz-2},{1,1,1},0,1);
Morphology Structure;
Structure.Initialize(Dm,SignDist);
int n;
double final_void_fraction;
@ -334,7 +620,7 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
if ( SignDist(i,j,k) > maxdist) maxdist=SignDist(i,j,k);
if ( SignDist(i,j,k) > 0.0 ){
count += 1.0;
id[n] = 2;
id[n] = ErodeLabel;
}
ID(i,j,k) = id[n];
}
@ -351,57 +637,6 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
if (rank==0) printf("Volume fraction for morphological opening: %f \n",volume_fraction);
if (rank==0) printf("Maximum pore size: %f \n",maxdistGlobal);
/* // Communication buffers
signed char *sendID_x, *sendID_y, *sendID_z, *sendID_X, *sendID_Y, *sendID_Z;
signed char *sendID_xy, *sendID_yz, *sendID_xz, *sendID_Xy, *sendID_Yz, *sendID_xZ;
signed char *sendID_xY, *sendID_yZ, *sendID_Xz, *sendID_XY, *sendID_YZ, *sendID_XZ;
signed char *recvID_x, *recvID_y, *recvID_z, *recvID_X, *recvID_Y, *recvID_Z;
signed char *recvID_xy, *recvID_yz, *recvID_xz, *recvID_Xy, *recvID_Yz, *recvID_xZ;
signed char *recvID_xY, *recvID_yZ, *recvID_Xz, *recvID_XY, *recvID_YZ, *recvID_XZ;
// send buffers
sendID_x = new signed char [Dm->sendCount("x")];
sendID_y = new signed char [Dm->sendCount("y")];
sendID_z = new signed char [Dm->sendCount("z")];
sendID_X = new signed char [Dm->sendCount("X")];
sendID_Y = new signed char [Dm->sendCount("Y")];
sendID_Z = new signed char [Dm->sendCount("Z")];
sendID_xy = new signed char [Dm->sendCount("xy")];
sendID_yz = new signed char [Dm->sendCount("yz")];
sendID_xz = new signed char [Dm->sendCount("xz")];
sendID_Xy = new signed char [Dm->sendCount("Xy")];
sendID_Yz = new signed char [Dm->sendCount("Yz")];
sendID_xZ = new signed char [Dm->sendCount("xZ")];
sendID_xY = new signed char [Dm->sendCount("xY")];
sendID_yZ = new signed char [Dm->sendCount("yZ")];
sendID_Xz = new signed char [Dm->sendCount("Xz")];
sendID_XY = new signed char [Dm->sendCount("XY")];
sendID_YZ = new signed char [Dm->sendCount("YZ")];
sendID_XZ = new signed char [Dm->sendCount("XZ")];
//......................................................................................
// recv buffers
recvID_x = new signed char [Dm->recvCount("x")];
recvID_y = new signed char [Dm->recvCount("y")];
recvID_z = new signed char [Dm->recvCount("z")];
recvID_X = new signed char [Dm->recvCount("X")];
recvID_Y = new signed char [Dm->recvCount("Y")];
recvID_Z = new signed char [Dm->recvCount("Z")];
recvID_xy = new signed char [Dm->recvCount("xy")];
recvID_yz = new signed char [Dm->recvCount("yz")];
recvID_xz = new signed char [Dm->recvCount("xz")];
recvID_Xy = new signed char [Dm->recvCount("Xy")];
recvID_xZ = new signed char [Dm->recvCount("xZ")];
recvID_xY = new signed char [Dm->recvCount("xY")];
recvID_yZ = new signed char [Dm->recvCount("yZ")];
recvID_Yz = new signed char [Dm->recvCount("Yz")];
recvID_Xz = new signed char [Dm->recvCount("Xz")];
recvID_XY = new signed char [Dm->recvCount("XY")];
recvID_YZ = new signed char [Dm->recvCount("YZ")];
recvID_XZ = new signed char [Dm->recvCount("XZ")];
//......................................................................................
int sendtag,recvtag;
sendtag = recvtag = 7;
*/
int ii,jj,kk;
int imin,jmin,kmin,imax,jmax,kmax;
int Nx = nx;
@ -422,7 +657,6 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
// if (rank==0) printf("Max. distance =%f, Initial critical radius = %f \n",maxdistGlobal,Rcrit_new);
//}
Dm->Comm.barrier();
FILE *DRAIN = fopen("morphdrain.csv","w");
fprintf(DRAIN,"sw radius\n");
@ -453,86 +687,36 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
for (jj=jmin; jj<jmax; jj++){
for (ii=imin; ii<imax; ii++){
double dsq = double((ii-i)*(ii-i)+(jj-j)*(jj-j)+(kk-k)*(kk-k));
if (ID(ii,jj,kk) == 2 && dsq <= (Rcrit_new+1)*(Rcrit_new+1)){
if (ID(ii,jj,kk) == ErodeLabel && dsq <= (Rcrit_new+1)*(Rcrit_new+1)){
LocalNumber+=1.0;
//id[nn]=1;
ID(ii,jj,kk)=1;
ID(ii,jj,kk)=NewLabel;
id[kk*Nx*Ny+jj*Nx+ii] = NewLabel;
}
}
}
}
}
// move on
}
}
}
LocalNumber += Structure.GetOverlaps(Dm,id,ErodeLabel,NewLabel);
for(int k=1; k<Nz-1; k++){
for(int j=1; j<Ny-1; j++){
for(int i=1; i<Nx-1; i++){
ID(i,j,k) = id[k*Nx*Ny+j*Nx+i];
}
}
}
fillChar.fill(ID);
// Pack and send the updated ID values
/* PackID(Dm->sendList("x"), Dm->sendCount("x") ,sendID_x, id);
PackID(Dm->sendList("X"), Dm->sendCount("X") ,sendID_X, id);
PackID(Dm->sendList("y"), Dm->sendCount("y") ,sendID_y, id);
PackID(Dm->sendList("Y"), Dm->sendCount("Y") ,sendID_Y, id);
PackID(Dm->sendList("z"), Dm->sendCount("z") ,sendID_z, id);
PackID(Dm->sendList("Z"), Dm->sendCount("Z") ,sendID_Z, id);
PackID(Dm->sendList("xy"), Dm->sendCount("xy") ,sendID_xy, id);
PackID(Dm->sendList("Xy"), Dm->sendCount("Xy") ,sendID_Xy, id);
PackID(Dm->sendList("xY"), Dm->sendCount("xY") ,sendID_xY, id);
PackID(Dm->sendList("XY"), Dm->sendCount("XY") ,sendID_XY, id);
PackID(Dm->sendList("xz"), Dm->sendCount("xz") ,sendID_xz, id);
PackID(Dm->sendList("Xz"), Dm->sendCount("Xz") ,sendID_Xz, id);
PackID(Dm->sendList("xZ"), Dm->sendCount("xZ") ,sendID_xZ, id);
PackID(Dm->sendList("XZ"), Dm->sendCount("XZ") ,sendID_XZ, id);
PackID(Dm->sendList("yz"), Dm->sendCount("yz") ,sendID_yz, id);
PackID(Dm->sendList("Yz"), Dm->sendCount("Yz") ,sendID_Yz, id);
PackID(Dm->sendList("yZ"), Dm->sendCount("yZ") ,sendID_yZ, id);
PackID(Dm->sendList("YZ"), Dm->sendCount("YZ") ,sendID_YZ, id);
//......................................................................................
Dm->Comm.sendrecv(sendID_x,Dm->sendCount("x"),Dm->rank_x(),sendtag,recvID_X,Dm->recvCount("X"),Dm->rank_X(),recvtag);
Dm->Comm.sendrecv(sendID_X,Dm->sendCount("X"),Dm->rank_X(),sendtag,recvID_x,Dm->recvCount("x"),Dm->rank_x(),recvtag);
Dm->Comm.sendrecv(sendID_y,Dm->sendCount("y"),Dm->rank_y(),sendtag,recvID_Y,Dm->recvCount("Y"),Dm->rank_Y(),recvtag);
Dm->Comm.sendrecv(sendID_Y,Dm->sendCount("Y"),Dm->rank_Y(),sendtag,recvID_y,Dm->recvCount("y"),Dm->rank_y(),recvtag);
Dm->Comm.sendrecv(sendID_z,Dm->sendCount("z"),Dm->rank_z(),sendtag,recvID_Z,Dm->recvCount("Z"),Dm->rank_Z(),recvtag);
Dm->Comm.sendrecv(sendID_Z,Dm->sendCount("Z"),Dm->rank_Z(),sendtag,recvID_z,Dm->recvCount("z"),Dm->rank_z(),recvtag);
Dm->Comm.sendrecv(sendID_xy,Dm->sendCount("xy"),Dm->rank_xy(),sendtag,recvID_XY,Dm->recvCount("XY"),Dm->rank_XY(),recvtag);
Dm->Comm.sendrecv(sendID_XY,Dm->sendCount("XY"),Dm->rank_XY(),sendtag,recvID_xy,Dm->recvCount("xy"),Dm->rank_xy(),recvtag);
Dm->Comm.sendrecv(sendID_Xy,Dm->sendCount("Xy"),Dm->rank_Xy(),sendtag,recvID_xY,Dm->recvCount("xY"),Dm->rank_xY(),recvtag);
Dm->Comm.sendrecv(sendID_xY,Dm->sendCount("xY"),Dm->rank_xY(),sendtag,recvID_Xy,Dm->recvCount("Xy"),Dm->rank_Xy(),recvtag);
Dm->Comm.sendrecv(sendID_xz,Dm->sendCount("xz"),Dm->rank_xz(),sendtag,recvID_XZ,Dm->recvCount("XZ"),Dm->rank_XZ(),recvtag);
Dm->Comm.sendrecv(sendID_XZ,Dm->sendCount("XZ"),Dm->rank_XZ(),sendtag,recvID_xz,Dm->recvCount("xz"),Dm->rank_xz(),recvtag);
Dm->Comm.sendrecv(sendID_Xz,Dm->sendCount("Xz"),Dm->rank_Xz(),sendtag,recvID_xZ,Dm->recvCount("xZ"),Dm->rank_xZ(),recvtag);
Dm->Comm.sendrecv(sendID_xZ,Dm->sendCount("xZ"),Dm->rank_xZ(),sendtag,recvID_Xz,Dm->recvCount("Xz"),Dm->rank_Xz(),recvtag);
Dm->Comm.sendrecv(sendID_yz,Dm->sendCount("yz"),Dm->rank_yz(),sendtag,recvID_YZ,Dm->recvCount("YZ"),Dm->rank_YZ(),recvtag);
Dm->Comm.sendrecv(sendID_YZ,Dm->sendCount("YZ"),Dm->rank_YZ(),sendtag,recvID_yz,Dm->recvCount("yz"),Dm->rank_yz(),recvtag);
Dm->Comm.sendrecv(sendID_Yz,Dm->sendCount("Yz"),Dm->rank_Yz(),sendtag,recvID_yZ,Dm->recvCount("yZ"),Dm->rank_yZ(),recvtag);
Dm->Comm.sendrecv(sendID_yZ,Dm->sendCount("yZ"),Dm->rank_yZ(),sendtag,recvID_Yz,Dm->recvCount("Yz"),Dm->rank_Yz(),recvtag);
//......................................................................................
UnpackID(Dm->recvList("x"), Dm->recvCount("x") ,recvID_x, id);
UnpackID(Dm->recvList("X"), Dm->recvCount("X") ,recvID_X, id);
UnpackID(Dm->recvList("y"), Dm->recvCount("y") ,recvID_y, id);
UnpackID(Dm->recvList("Y"), Dm->recvCount("Y") ,recvID_Y, id);
UnpackID(Dm->recvList("z"), Dm->recvCount("z") ,recvID_z, id);
UnpackID(Dm->recvList("Z"), Dm->recvCount("Z") ,recvID_Z, id);
UnpackID(Dm->recvList("xy"), Dm->recvCount("xy") ,recvID_xy, id);
UnpackID(Dm->recvList("Xy"), Dm->recvCount("Xy") ,recvID_Xy, id);
UnpackID(Dm->recvList("xY"), Dm->recvCount("xY") ,recvID_xY, id);
UnpackID(Dm->recvList("XY"), Dm->recvCount("XY") ,recvID_XY, id);
UnpackID(Dm->recvList("xz"), Dm->recvCount("xz") ,recvID_xz, id);
UnpackID(Dm->recvList("Xz"), Dm->recvCount("Xz") ,recvID_Xz, id);
UnpackID(Dm->recvList("xZ"), Dm->recvCount("xZ") ,recvID_xZ, id);
UnpackID(Dm->recvList("XZ"), Dm->recvCount("XZ") ,recvID_XZ, id);
UnpackID(Dm->recvList("yz"), Dm->recvCount("yz") ,recvID_yz, id);
UnpackID(Dm->recvList("Yz"), Dm->recvCount("Yz") ,recvID_Yz, id);
UnpackID(Dm->recvList("yZ"), Dm->recvCount("yZ") ,recvID_yZ, id);
UnpackID(Dm->recvList("YZ"), Dm->recvCount("YZ") ,recvID_YZ, id);
//......................................................................................
// double GlobalNumber = Dm->Comm.sumReduce( LocalNumber );
*/
Dm->Comm.barrier();
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;
if (ID(i,j,k) == 1){
if (ID(i,j,k) == NewLabel){
phase(i,j,k) = 1.0;
}
else
@ -551,50 +735,12 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
for (int i=0; i<nx; i++){
n=k*nx*ny+j*nx+i;
if (ID(i,j,k) == 1 && phase_label(i,j,k) > 1){
ID(i,j,k) = 2;
ID(i,j,k) = ErodeLabel;
}
id[n] = ID(i,j,k);
}
}
}
/*
* Extract only the connected part of NWP
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] == 2){
phase(i,j,k) = 1.0;
}
else if (id[n] == 1){
// nwp
phase(i,j,k) = -1.0;
}
else{i
// treat solid as WP since films can connect
phase(i,j,k) = 1.0;
}
}
}
}
ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,SignDist,vF,vS,phase_label,Dm->Comm);
Dm->Comm.barrier();
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] == 2 && phase_label(i,j,k) > 1){
id[n] = 20;
}
}
}
}
// done
*/
count = 0.f;
for (int k=1; k<nz-1; k++){

90
analysis/morphology.h
View File

@ -6,3 +6,93 @@
double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain> Dm, double VoidFraction, signed char ErodeLabel, signed char ReplaceLabel);
double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain> Dm, double VoidFraction);
double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id, std::shared_ptr<Domain> Dm, double TargetVol, double WallFactor);
#ifndef MORPHOLOGY_INC
#define MORPHOLOGY_INC
/**
* \class Morphology
* @brief
* The Morphology class supports morphological operations on complex structures
*
*/
class Morphology{
public:
/**
* \brief Create a flow adaptor to operate on the LB model
*/
Morphology();
/**
* \brief Destructor
*/
~Morphology();
/**
* \brief Initialize morphology structure from distance map
* @param Dm Domain structure
* @param Distance Signed distance to boundary of structure
*/
void Initialize(std::shared_ptr <Domain> Dm, DoubleArray &Distance);
/**
* \brief Find all sites such that the reach of the signed distance at the site overlaps with a sub-domain boundary
* @param Dm Domain structure
* @param id image labels
* @param ErodeLabel label to erode based on morphological operation
* @param NewLabel label to assign based on morphological operation
*/
int GetOverlaps(std::shared_ptr <Domain> Dm, signed char *id, const signed char ErodeLabel, const signed char NewLabel);
/*
* data structures to store non-local morphological information
*/
std::vector<int> xShift, yShift, zShift;
std::vector<int> sendID;
std::vector<double> morphRadius;
std::vector<unsigned char> localID;
std::vector<unsigned char> nonlocalID;
private:
int sendtag,recvtag;
//......................................................................................
int sendCount, recvCount;
//......................................................................................
int sendOffset_x, sendOffset_y, sendOffset_z, sendOffset_X, sendOffset_Y, sendOffset_Z;
int sendOffset_xy, sendOffset_yz, sendOffset_xz, sendOffset_Xy, sendOffset_Yz, sendOffset_xZ;
int sendOffset_xY, sendOffset_yZ, sendOffset_Xz, sendOffset_XY, sendOffset_YZ, sendOffset_XZ;
int sendOffset_xyz, sendOffset_XYZ, sendOffset_xYz, sendOffset_XyZ;
int sendOffset_Xyz, sendOffset_xYZ, sendOffset_xyZ, sendOffset_XYz;
//......................................................................................
int recvOffset_x, recvOffset_y, recvOffset_z, recvOffset_X, recvOffset_Y, recvOffset_Z;
int recvOffset_xy, recvOffset_yz, recvOffset_xz, recvOffset_Xy, recvOffset_Yz, recvOffset_xZ;
int recvOffset_xY, recvOffset_yZ, recvOffset_Xz, recvOffset_XY, recvOffset_YZ, recvOffset_XZ;
int recvOffset_xyz, recvOffset_XYZ, recvOffset_xYz, recvOffset_XyZ;
int recvOffset_Xyz, recvOffset_xYZ, recvOffset_xyZ, recvOffset_XYz;
//......................................................................................
int sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z;
int sendCount_xy, sendCount_yz, sendCount_xz, sendCount_Xy, sendCount_Yz, sendCount_xZ;
int sendCount_xY, sendCount_yZ, sendCount_Xz, sendCount_XY, sendCount_YZ, sendCount_XZ;
int sendCount_xyz, sendCount_XYZ, sendCount_xYz, sendCount_XyZ;
int sendCount_Xyz, sendCount_xYZ, sendCount_xyZ, sendCount_XYz;
//......................................................................................
int recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y, recvCount_Z;
int recvCount_xy, recvCount_yz, recvCount_xz, recvCount_Xy, recvCount_Yz, recvCount_xZ;
int recvCount_xY, recvCount_yZ, recvCount_Xz, recvCount_XY, recvCount_YZ, recvCount_XZ;
int recvCount_xyz, recvCount_XYZ, recvCount_xYz, recvCount_XyZ;
int recvCount_Xyz, recvCount_xYZ, recvCount_xyZ, recvCount_XYz;
//......................................................................................
std::vector<char> sendList;
std::vector<char> recvList;
//......................................................................................
// Communication buffers
signed char *sendID_x, *sendID_y, *sendID_z, *sendID_X, *sendID_Y, *sendID_Z;
signed char *sendID_xy, *sendID_yz, *sendID_xz, *sendID_Xy, *sendID_Yz, *sendID_xZ;
signed char *sendID_xY, *sendID_yZ, *sendID_Xz, *sendID_XY, *sendID_YZ, *sendID_XZ;
signed char *recvID_x, *recvID_y, *recvID_z, *recvID_X, *recvID_Y, *recvID_Z;
signed char *recvID_xy, *recvID_yz, *recvID_xz, *recvID_Xy, *recvID_Yz, *recvID_xZ;
signed char *recvID_xY, *recvID_yZ, *recvID_Xz, *recvID_XY, *recvID_YZ, *recvID_XZ;
};
#endif

BIN
docs/source/_static/images/DiscPack-morphdrain.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 9.5 KiB

BIN
docs/source/_static/images/lbpm-visit-workflow-i.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 516 KiB

BIN
docs/source/_static/images/lbpm-visit-workflow-ii.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 371 KiB

BIN
docs/source/_static/images/lbpm-visit-workflow-iii.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 455 KiB

BIN
docs/source/_static/images/lbpm-visit-workflow-iv.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 1.9 MiB

5
docs/source/developerGuide/doxygen/analysis/Morphology.rst Normal file
View File

@ -0,0 +1,5 @@
============================================
Morphology
============================================
.. doxygenfile:: morphology.h
:project: LBPM Doxygen

5
docs/source/examples/color/steadyState.rst Normal file
View File

@ -0,0 +1,5 @@
*******************
Steady-state flow
*******************
In this example we simulate a steady-state flow with a constant driving force.

52
docs/source/examples/domain/domain.rst Normal file
View File

@ -0,0 +1,52 @@
*****************
Input Domain
*****************
LBPM provides a flexible framework to ingest 3D image data.
To illustrate the basic capabilities, this tutorial considers a quasi-2D
flow cell. Source files for the example are included in the LBPM repository
in the directory ``examples/DiscPack``. A simple python code is included
to set up the flow domain.
Based on LBPM convention, external boundary conditions are applied in the
z-direction. This means that the domain should be set up so that the direction
you want to set boundary conditions is aligned with the z-axis. For the quasi-2D
example, a depth of ``3`` voxels is used for the x-direction. *Based on LBPM
internal data structures at least three voxels must be provided in each direction*
The specified domain decomposition must also observe this rule.
Image data is stored internally within LBPM as signed 8-bit binary data. This means that
up to 256 labels can be provided in the input image. LBPM convention takes all
non-positive labels to be immobile (treated as solid). In this example, the solid regions
are assigned a value of ``0``. It is possible to provide up to ``128`` different labels
for the solid. Also, note that python supports only the unsigned 8-bit datatype. For the unsigned data
type, labels assigned values ``128,...255`` in python will correspond to labels
``-127,...-1`` when read in as type ``signed char`` within LBPM.
.. code:: python
import numpy as np
import matplotlib.pylab as plt
import pandas as pd
# Set the size of the domain
Nx=3
Ny=128
Nz=128
D=pd.read_csv("discs.csv",sep=" ")
ID = np.ones(Nx*Ny*Nz,dtype='uint8')
ID.shape = (Nz,Ny,Nx)
# Set the solid labels
for idx in range(len(D)):
cx=D['cx'][idx] / dx
cy=D['cy'][idx] /dx
r=D['r'][idx] /dx
for i in range(0,Nz):
for j in range(0,Ny):
if ( (cx-i)*(cx-i) + (cy-j)*(cy-j) < r*r ):
ID[i,j,0] = 0
ID[i,j,1] = 0
ID[i,j,2] = 0
# write input file to disc
ID.tofile("discs_3x128x128.raw")

110
docs/source/examples/morphology/morphOpen.rst Normal file
View File

@ -0,0 +1,110 @@
*****************************
Morphological pre-processors
*****************************
LBPM includes morphological pre-processing tools as utility functions.
It is often useful to generate initial conditions for a 2-phase flow simulation based on a morphological approach. In particular, morphological tools can be used to provide a physical reasonable initial condition in cases where direct experimental observations are not available. These initial configurations are compatible with any of the 2-phase simulation protocols used by lbpm_color_simulator. These initialization approaches alter the fluid labels within the input files, writing a new file with the new morphologically assigned labels.
There are two main morphological pre-processors in LBPM
* ``lbpm_morphdrain_pp`` -- initialize fluid configuration based on morphological drainage
* ``lbpm_morphopen_pp`` -- initialize fluid configuration based on morphological opening
Here we demonstrate ``lbpm_morphdrain_pp`` because it offers the most information. Although it is not perfect, the morphological drainage operation does a good job of approximating configurations observed along primary drainage processes performed under water-wet conditions. A limitation is that fluid trapped in the corners will not stop the morphological operation from eroding it. This should not discourage you too much -- morphological tools are very practical and can save you a lot of time! It is also a good thing to be skeptical.
Since the morphological operation works on the input domain, associated parameters are added to the ``Domain`` section of the input file. Here we will set a target saturation Sw = 0.20, which will run the morphological drainage operation until the fluid labeled as 2 occupies 20% of the pore space or less. For the case considered in
``example/DiscPack`` we specify the following information in the input file
.. code:: c
Domain {
Filename = "discs_3x128x128.raw"
ReadType = "8bit" // data type
N = 3, 128, 128 // size of original image
nproc = 1, 2, 2 // process grid
n = 3, 64, 64 // sub-domain size
voxel_length = 1.0 // voxel length (in microns)
ReadValues = 0, 1, 2 // labels within the original image
WriteValues = 0, 2, 2 // associated labels to be used by LBPM
BC = 0 // fully periodic BC
Sw = 0.35 // target saturation for morphological tools
}
Once this has been set, we launch lbpm_morphdrain_pp in the same way as other parallel tools
.. code:: bash
mpirun -np 4 $LBPM_BIN/lbpm_morphdrain_pp input.db
Successful output looks like the following
.. code:: bash
Performing morphological opening with target saturation 0.350000
voxel length = 1.000000 micron
voxel length = 1.000000 micron
Input media: discs_3x128x128.raw
Relabeling 3 values
oldvalue=0, newvalue =0
oldvalue=1, newvalue =2
oldvalue=2, newvalue =2
Dimensions of segmented image: 3 x 128 x 128
Reading 8-bit input data
Read segmented data from discs_3x128x128.raw
Label=0, Count=11862
Label=1, Count=37290
Label=2, Count=0
Distributing subdomains across 4 processors
Process grid: 1 x 2 x 2
Subdomain size: 3 x 64 x 64
Size of transition region: 0
Media porosity = 0.758667
Initialized solid phase -- Converting to Signed Distance function
Volume fraction for morphological opening: 0.758667
Maximum pore size: 116.773801
1.000000 110.935111
1.000000 105.388355
1.000000 100.118937
1.000000 95.112990
1.000000 90.357341
1.000000 85.839474
1.000000 81.547500
1.000000 77.470125
1.000000 73.596619
1.000000 69.916788
1.000000 66.420949
1.000000 63.099901
1.000000 59.944906
1.000000 56.947661
1.000000 54.100278
1.000000 51.395264
1.000000 48.825501
1.000000 46.384226
1.000000 44.065014
1.000000 41.861764
1.000000 39.768675
1.000000 37.780242
1.000000 35.891230
1.000000 34.096668
1.000000 32.391835
0.575114 30.772243
0.433119 29.233631
0.291231 27.771949
Final void fraction =0.291231
Final critical radius=27.771949
Writing ID file
Writing file to: discs_3x128x128.raw.morphdrain.raw
The final configuration can be visualized in python by loading the output file
``discs_3x128x128.raw.morphdrain.raw``.
.. figure:: ../../_static/images/DiscPack-morphdrain.png
:width: 600
:alt: morphdrain
Fluid configuration resulting from orphological drainage algorithm applied to a 2D disc pack.

15
docs/source/examples/running.rst
View File

@ -1,10 +1,21 @@
==============
LBPM examples
Examples
==============
There are two main components to running LBPM simulators.
First is understanding how to launch MPI tasks on your system,
which depends on the particular implementation of MPI that you are using,
as well as other details of the local configuration. The second component is
understanding the LBPM input file structure.
understanding the LBPM input file structure. The examples included provide
a basic introduction to working with LBPM.
.. toctree::
:glob:
:maxdepth: 2
domain/*
morphology/*
color/*

2
docs/source/index.rst
View File

@ -12,9 +12,9 @@ LBPM -- Documentation
:caption: Contents:
install
examples/*
userGuide/*
developerGuide/*
examples/*
publications/*
Indices and tables

6
docs/source/userGuide/IO/fileformat.rst
View File

@ -5,9 +5,9 @@ I/O conventions for LBPM
There are three main kinds of output file that are supported by LBPM.
* CSV files --
* CSV files -- space-delimited CSV files are used by the internal analysis framework
* formatted binary files --
* formatted binary files -- SILO and HDF5 formats are supported for visualization data
* unformatted binary files --
* unformatted binary files -- ``.raw`` extension

5
docs/source/userGuide/analysis/analysis.rst
View File

@ -1,5 +0,0 @@
===========================
Internal Analysis Framework
===========================
placeholder for analysis

2
docs/source/userGuide/index.rst
View File

@ -10,8 +10,6 @@ Welcome to the LBPM user guide.
models/*
analysis/*
visualization/*
IO/*

5
docs/source/userGuide/models/color/index.rst
View File

@ -266,3 +266,8 @@ the inlet or outlet, the ``Domain`` section of the database may specify the foll
- ``InletLayerPhase = 2`` -- establish a reservoir of component B at the inlet
- ``OutletLayerPhase = 1`` -- establish a reservoir of component A at the outlet
****************
Example data
****************
Example data can be downloaded from https://www.digitalrocksportal.org/projects/326

2
docs/source/userGuide/models/color/protocols/fractionalFlow.rst
View File

@ -52,7 +52,7 @@ where :math:`\max{|\mathbf{u}_i|}` is the maximum flow speed within fluid :math:
:math:`U_\epsilon` is a threshold speed that is set to minimize the influence of spurious
currents on the mass seeding algorithm. The sum of the weighting function is used to normalize
the local weights so that the added mass will match the value specified by ``mass_fraction_factor``.
If the flow is slower than :math:`\epsilon_m`, the algorithm will tend to add mass evenly to the system.
If the flow is slower than :math:`U_\epsilon`, the algorithm will tend to add mass evenly to the system.
For example, if the water is only present in films that flow very slowly, then mass will
be evenly seeded throughout entire water film. Alternatively, if one or both fluids
flows through distinct channels, the mass will be disproportionately added to these

49
docs/source/userGuide/visualization/visit.rst
View File

@ -2,6 +2,51 @@
Visualizing simulation data with visit
======================================
placeholder for visit
Paraview > v 5.6 also works
(Paraview > v 5.6 also works)
Control over visualization options is set within the ``Visualization`` section of the input file
.. code:: c
Visualization {
write_silo = true // write SILO databases with assigned variables
save_8bit_raw = true // write labeled 8-bit binary files with phase assignments
save_phase_field = true // save phase field within SILO database
save_pressure = false // save pressure field within SILO database
save_velocity = false // save velocity field within SILO database
}
LBPM provides two main options for visualization. The first is the writing of 8-bit raw binary files, which are labeled based on the timestep. For example, if ``visualization_interval = 10000`` (specified within the Analysis section of the input file) the first 8-bit binary file will be written when ``timestep = 1000`` and will be named ``id_t1000.raw``. Additional files will be written subsequently at the specified interval. Similarly, higher fidelity visualization files are written using the SILO format, which are stored within the directories ``vis1000/``. The summary file ``LBM.visit`` enumerates these files so that they can be loaded directly into VisIt or other visualization software. By default, only the phase field will be saved. Visualization for other variables, such as the pressure and velocity fields, can be enabled by setting the associated flags to ``true``.
The VisIt software is able to natively read the SILO format. To import the data fields written by LBPM, open the VisIt GUI and select ``File > Open file`` from the top menu. Then select the LBM.visit file that you would like to read
.. figure:: ../../_static/images/lbpm-visit-workflow-i.png
:width: 600
:alt: VisIt GUI
Opening data in the VisIt GUI.
Once the file has been opened, any database fields within the file will be visible from within the program. Here we construct an 3D countour the phase field to visualize the boundary of the oil. The menu for the ``Contour`` object will show that the data for phase is available. The pressure and velocity fields will only be visible if they have been explicitly enabled within the simulation options (see ``Visualization`` details above)
.. figure:: ../../_static/images/lbpm-visit-workflow-ii.png
:width: 600
:alt: VisIt GUI
Selecting isosurface contour to represent a surface.
Once the contour has been created, double click the object icon to open the countour plot attributes window. In the field for ``Select by`` choose ``values`` and specify ``0`` as the contour to draw. Note that it is possible to change the color and opacity of the contour, which is useful if you want to draw multiple contours on the same plot
.. figure:: ../../_static/images/lbpm-visit-workflow-iii.png
:width: 600
:alt: VisIt GUI
Drawing an isosurface.
Once the attributes have been selected, click the Draw button to render the contour. Depending on the machine where you are rendering and the size of the image, it may take several minutes to render the window
.. figure:: ../../_static/images/lbpm-visit-workflow-iv.png
:width: 600
:alt: VisIt GUI
Rendering an isosurface.

150
example/DiscPack/DiscPack.ipynb
View File

File diff suppressed because one or more lines are too long

6
example/DiscPack/discs.csv Normal file
View File

@ -0,0 +1,6 @@
cx cy r
0.25 0.2 0.12
0.45 0.52 0.15
0.7 0.8 0.12
0.8 0.25 0.13
0.2 0.85 0.1
1 cx cy r
2 0.25 0.2 0.12
3 0.45 0.52 0.15
4 0.7 0.8 0.12
5 0.8 0.25 0.13
6 0.2 0.85 0.1

1
tests/CMakeLists.txt
View File

@ -73,6 +73,7 @@ ADD_LBPM_TEST( TestMomentsD3Q19 )
ADD_LBPM_TEST( TestInterfaceSpeed ../example/Bubble/input.db)
ADD_LBPM_TEST( test_dcel_minkowski )
ADD_LBPM_TEST( test_dcel_tri_normal )
ADD_LBPM_TEST( test_morph )
ADD_LBPM_TEST( TestMassConservationD3Q7 ../example/Bubble/input.db)
#ADD_LBPM_TEST_1_2_4( TestTwoPhase )
ADD_LBPM_TEST_1_2_4( TestBlobIdentify )

117
tests/test_morph.cpp Normal file
View File

@ -0,0 +1,117 @@
#include <iostream>
#include <math.h>
#include "analysis/Minkowski.h"
#include "analysis/morphology.h"
#include "common/Domain.h"
#include "common/SpherePack.h"
#include "ProfilerApp.h"
/*
* Compare the measured and analytical curvature for a sphere
*
*/
std::shared_ptr<Database> loadInputs( )
{
//auto db = std::make_shared<Database>( "Domain.in" );
auto db = std::make_shared<Database>();
db->putScalar<int>( "BC", 0 );
db->putVector<int>( "nproc", { 1, 1, 1 } );
db->putVector<int>( "n", { 64, 64, 64 } );
db->putScalar<int>( "nspheres", 1 );
db->putVector<double>( "L", { 1, 1, 1 } );
return db;
}
int main(int argc, char **argv)
{
Utilities::startup( argc, argv );
Utilities::MPI comm( MPI_COMM_WORLD );
int toReturn = 0;
{
int i,j,k;
// Load inputs
auto db = loadInputs( );
int Nx = db->getVector<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "n" )[2];
auto Dm = std::make_shared<Domain>( db, comm );
Dm->CommInit();
Nx+=2; Ny+=2; Nz+=2;
int N = Nx*Ny*Nz;
signed char *ID;
ID = new signed char[N];
DoubleArray Phase(Nx,Ny,Nz);
DoubleArray Distance(Nx,Ny,Nz);
Minkowski sphere(Dm);
printf("Set distance map for a sphere \n");
for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
int n = k*Nx*Ny + j*Nx + i;
// put sphere with center at (0.8*Nx, 0.5*Ny, 0.5*Nz)
if (i+1 > 0.3*(Nx-2)){
Distance(i,j,k) = 0.4*(Nx-2) - sqrt((1.0*(i-1)-0.8*(Nx-2))*(1.0*i-0.8*(Nx-2))+(1.0*(j-1)-
0.5*(Ny-2))*(1.0*(j-1)-0.5*(Ny-2))+(1.0*(k-1)-0.5*(Nz-2))*(1.0*(k-1)-0.5*(Nz-2)));
}
// reflection sphere center at (-0.2*Nx, 0.5*Ny, 0.5*Nz)
else {
Distance(i,j,k) = 0.4*(Nx-2) - sqrt((1.0*(i-1)+0.2*(Nx-2))*(1.0*(i-1)+0.2*(Nx-2))+(1.0*(j-1)-
0.5*(Ny-2))*(1.0*(j-1)-0.5*(Ny-2))+(1.0*(k-1)-0.5*(Nz-2))*(1.0*(k-1)-0.5*(Nz-2)));
}
/* set the labels */
if (Distance(i,j,k) > 0.0){
ID[n] = 2;
}
else {
ID[n] = 0;
}
}
}
}
printf("Perform morphological opening \n");
signed char ErodeLabel = 2;
signed char NewLabel = 1;
double Saturation = 0.9;
MorphOpen(Distance,ID,Dm,Saturation,ErodeLabel,NewLabel);
for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
int n = k*Nx*Ny + j*Nx + i;
sphere.id(i,j,k) = 1;
if (ID[n] == 1)
sphere.id(i,j,k) = 0;
}
}
}
printf(" Measure the opening \n");
sphere.MeasureObject();
//sphere.ComputeScalar(Distance,0.f);
printf(" Volume = %f (analytical = %f) \n", sphere.Vi,0.256*0.33333333333333*3.14159*double((Nx-2)*(Nx-2)*(Nx-2)));
double error = fabs(sphere.Vi - 0.256*0.33333333333333*3.14159*double((Nx-2)*(Nx-2)*(Nx-2)))/ (0.256*0.33333333333333*3.14159*double((Nx-2)*(Nx-2)*(Nx-2)));
printf(" Relative error %f \n", error);
if (error > 0.03){
toReturn = 10;
printf("ERROR (test_morph): difference from analytical volume is too large\n");
FILE *OUTFILE;
OUTFILE = fopen("test_morph_ID.raw","wb");
fwrite(ID,1,N,OUTFILE);
fclose(OUTFILE);
}
else {
printf("SUCCESS\n");
}
}
PROFILE_SAVE("test_morph");
Utilities::shutdown();
return toReturn;
}