generalize morphopen

analysis/morphology.cpp

@@ -144,7 +144,6 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
 	//	Rcrit_new = strtod(argv[2],NULL);
 	//	if (rank==0) printf("Max. distance =%f, Initial critical radius = %f \n",maxdistGlobal,Rcrit_new);
 	//}
-	
 	while (void_fraction_new > VoidFraction)
 	{
 		void_fraction_diff_old = void_fraction_diff_new;
@@ -172,7 +171,11 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
 								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));
+<<<<<<< HEAD
 									if (id[nn] == ErodeLabel && dsq <= Rcrit_new*Rcrit_new){
+=======
+									if (id[nn] == 2 && dsq <= (Rcrit_new+1)*(Rcrit_new+1)){
+>>>>>>> fe5a1100278e24d46e9fb65fae8c837afa9a19a6
 										LocalNumber+=1.0;
 										id[nn]=NewLabel;
 									}
@@ -461,9 +464,9 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
 		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
 		double LocalNumber=0.f;
-		for(int k=0; k<Nz; k++){
-			for(int j=0; j<Ny; j++){
-				for(int i=0; i<Nx; i++){
+		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 (SignDist(i,j,k) > Rcrit_new){
 						// loop over the window and update
@@ -478,7 +481,7 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
 								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] == 2 && dsq <= Rcrit_new*Rcrit_new){
+									if (id[nn] == 2 && dsq <= (Rcrit_new+1)*(Rcrit_new+1)){
 										LocalNumber+=1.0;
 										id[nn]=1;
 									}

example/CornerFlow/input.db

@@ -0,0 +1,49 @@
+Color {
+    tauA   = 0.7;  
+    tauB   = 0.7;  
+    rhoA   = 1.0;  
+    rhoB   = 1.0;  
+    alpha = 1e-3;
+    beta  = 0.95;
+    F = 0, 0, 0
+    Restart = false
+    pBC = 0
+    din = 1.0
+    dout = 1.0
+    timestepMax = 3000
+    interval = 1000
+    tol = 1e-5;
+    das = 0.1
+    dbs = 0.9
+    flux = 0.0
+    ComponentLabels = 0, -1        // labels for immobile components
+    ComponentAffinity = -1.0, 0.5  // wetting condition for immobile components (-1 = water wet / +1 oil wet)
+}
+
+Domain {
+    Filename = "CornerFlow.raw"
+    nproc = 1, 1, 2     // Number of processors (Npx,Npy,Npz)
+    n = 64, 64, 64      // Size of local domain (Nx,Ny,Nz)
+    N = 64, 64, 128      // size of the input image
+    n_spheres = 1       // Number of spheres
+    L = 1, 1, 1         // Length of domain (x,y,z)
+    BC = 0              // Boundary condition type
+    ReadType = "8bit"
+    ReadValues = 0, 1, 2  // labels within the original input image
+    WriteValues = 0, 1, 2 // labels to write (if they aren't the same)
+    ComponentLabels = 0   // labels that are immobile
+    HistoryLabels = -1    // new labels to assign to each immobile component based on fluid history
+    Sw = 0.3              // target saturation for morphological routines
+} 
+
+Analysis {
+    blobid_interval = 1000      // Frequency to perform blob identification
+    analysis_interval = 1000    // Frequency to perform analysis
+    restart_interval = 1000    // Frequency to write restart data
+    visualization_interval = 2000        // Frequency to write visualization data
+    restart_file = "Restart"    // Filename to use for restart file (will append rank)
+    N_threads    = 4            // Number of threads to use
+    load_balance = "independent" // Load balance method to use: "none", "default", "independent"
+}
+
+

gpu/Extras.cu

@@ -9,6 +9,7 @@ extern "C" int ScaLBL_SetDevice(int rank){
 	//int device = local_rank % n_devices; 
 	int device = rank % n_devices; 
 	cudaSetDevice(device); 
+	if (rank < n_devices) printf("MPI rank=%i will use GPU ID %i / %i \n",rank,device,n_devices);
 	return device;
 }
 

tests/lbpm_color_simulator.cpp

@@ -41,7 +41,6 @@ int main(int argc, char **argv)
 	}
 	// Initialize compute device
 	int device=ScaLBL_SetDevice(rank);
-	if (rank<4)	printf("Using GPU ID %i for rank %i \n",device,rank);
 	ScaLBL_DeviceBarrier();
 	MPI_Barrier(comm);
 	

tests/lbpm_morphdrain_pp.cpp

@@ -173,7 +173,7 @@ int main(int argc, char **argv)
 							signed char NEWVALUE=HistoryLabels[idx];
 							if (LOCVAL == VALUE){
 								idx = NLABELS;
-								if (SignDist(i,j,k) < 1.0){
+								if (SignDist(i,j,k) < 2.0){
 									id[n] = NEWVALUE;
 								}
 							}

tests/lbpm_morphopen_pp.cpp

@@ -173,7 +173,7 @@ int main(int argc, char **argv)
 							signed char NEWVALUE=HistoryLabels[idx];
 							if (LOCVAL == VALUE){
 								idx = NLABELS;
-								if (SignDist(i,j,k) < 1.0){
+								if (SignDist(i,j,k) < 2.0){
 									id[n] = NEWVALUE;
 								}
 							}

tests/lbpm_permeability_simulator.cpp

@@ -44,7 +44,6 @@ int main(int argc, char **argv)
 		}
 		// Initialize compute device
 		int device=ScaLBL_SetDevice(rank);
-		if (rank<4)	printf("Using GPU ID %i for rank %i \n",device,rank);
 		ScaLBL_DeviceBarrier();
 		MPI_Barrier(comm);