diff --git a/common/Domain.h b/common/Domain.h
index a70181ce..a8284cc4 100755
--- a/common/Domain.h
+++ b/common/Domain.h
@@ -1127,102 +1127,6 @@ inline void SignedDistance(double *Distance, int nspheres, double *List_cx, doub
 	for (n=0; n<N; n++)	Distance[n] = Distance[n]/hx;
 }
 
-inline void GenerateResidual(char *ID, int Nx, int Ny, int Nz, double Saturation)
-{
-	//.......................................................................
-	int i,j,k,n,Number,N;
-	int x,y,z,ii,jj,kk;
-	int sizeX,sizeY,sizeZ;
-	int *SizeX, *SizeY, *SizeZ;
-
-#ifdef NORANDOM
-	srand(10009);
-#else
-	srand(time(NULL));
-#endif
-//	float bin;
-	//.......................................................................
-	N = Nx*Ny*Nz;
-	
-	int bin, binCount;	
-	ifstream Dist("BlobSize.in");
-	Dist >> binCount;
-//	printf("Number of blob sizes: %i \n",binCount);
-	SizeX = new int [binCount];
-	SizeY = new int [binCount];
-	SizeZ = new int [binCount];
-	for (bin=0; bin<binCount; bin++){
-		Dist >> SizeX[bin];
-		Dist >> SizeY[bin];
-		Dist >> SizeZ[bin];
-	//	printf("Blob %i dimension: %i x %i x %i \n",bin, SizeX[bin], SizeY[bin], SizeZ[bin]);
-	}
-	Dist.close();
-	//.......................................................................
-//	cout << "Generating blocks... " << endl;	
-	// Count for the total number of oil nodes 
-	int count = 0;
-	// Count the total number of non-solid nodes
-	int total = 0;
-	for (i=0;i<N;i++){
-		if (ID[i] != 0) total++;
-	}
-	
-	float sat = 0.f;
-	Number = 0;		// number of features
-	while (sat < Saturation){
-		Number++;
-		// Randomly generate a point in the domain
-		x = Nx*float(rand())/float(RAND_MAX);
-		y = Ny*float(rand())/float(RAND_MAX);
-		z = Nz*float(rand())/float(RAND_MAX);
-		
-		bin = int(floor(binCount*float(rand())/float(RAND_MAX)));
-		sizeX = SizeX[bin];
-		sizeY = SizeY[bin];
-		sizeZ = SizeZ[bin];
-		
-//		cout << "Sampling from bin no. " << floor(bin) << endl; 
-//		cout << "Feature size is: " << sizeX << "x" << sizeY << "x" << sizeZ << endl; 
-		
-		for (k=z;k<z+sizeZ;k++){
-			for (j=y;j<y+sizeY;j++){
-				for (i=x;i<x+sizeX;i++){
-					// Identify nodes in the domain (periodic BC)
-					ii = i;
-					jj = j;
-					kk = k;					
-					if (ii < 1)			ii+=(Nx-2);
-					if (jj < 1)			jj+=(Ny-2);
-					if (kk < 1)			kk+=(Nz-2);
-					if (!(ii < Nx-1))		ii-=(Nx-2);
-					if (!(jj < Ny-1))		jj-=(Ny-2);
-					if (!(kk < Nz-1))		kk-=(Nz-2);
-					
-					n = kk*Nx*Ny+jj*Nx+ii;
-					
-					if (ID[n] == 2){
-						ID[n] = 1;
-						count++;
-					}
-				}
-			}
-		}
-		sat = float(count)/total;
-	}
-	//.......................................................................
-}
-
-
-
-inline void FlipID(char *ID, int N)
-{
-	for (int n=0; n<N; n++){
-		if  	 (ID[n] == 1)	ID[n] = 2;
-		else if  (ID[n] == 2)	ID[n] = 1;
-	}
-}
-
 inline void WriteLocalSolidID(char *FILENAME, char *ID, int N)
 {
 	char value;
diff --git a/tests/lb2_Color_wia_mpi.cpp b/tests/lb2_Color_wia_mpi.cpp
index 84f07f56..d7ea3ff4 100644
--- a/tests/lb2_Color_wia_mpi.cpp
+++ b/tests/lb2_Color_wia_mpi.cpp
@@ -23,6 +23,104 @@
  * James E. McClure 2013-2014
  */
 
+
+// The functions GenerateResidual and FlipID may be deprecated
+// see lbpm_random_pp for latest version
+inline void GenerateResidual(char *ID, int Nx, int Ny, int Nz, double Saturation)
+{
+	//.......................................................................
+	int i,j,k,n,Number,N;
+	int x,y,z,ii,jj,kk;
+	int sizeX,sizeY,sizeZ;
+	int *SizeX, *SizeY, *SizeZ;
+
+#ifdef NORANDOM
+	srand(10009);
+#else
+	srand(time(NULL));
+#endif
+//	float bin;
+	//.......................................................................
+	N = Nx*Ny*Nz;
+
+	int bin, binCount;
+	ifstream Dist("BlobSize.in");
+	Dist >> binCount;
+//	printf("Number of blob sizes: %i \n",binCount);
+	SizeX = new int [binCount];
+	SizeY = new int [binCount];
+	SizeZ = new int [binCount];
+	for (bin=0; bin<binCount; bin++){
+		Dist >> SizeX[bin];
+		Dist >> SizeY[bin];
+		Dist >> SizeZ[bin];
+	//	printf("Blob %i dimension: %i x %i x %i \n",bin, SizeX[bin], SizeY[bin], SizeZ[bin]);
+	}
+	Dist.close();
+	//.......................................................................
+//	cout << "Generating blocks... " << endl;
+	// Count for the total number of oil nodes
+	int count = 0;
+	// Count the total number of non-solid nodes
+	int total = 0;
+	for (i=0;i<N;i++){
+		if (ID[i] != 0) total++;
+	}
+
+	float sat = 0.f;
+	Number = 0;		// number of features
+	while (sat < Saturation){
+		Number++;
+		// Randomly generate a point in the domain
+		x = Nx*float(rand())/float(RAND_MAX);
+		y = Ny*float(rand())/float(RAND_MAX);
+		z = Nz*float(rand())/float(RAND_MAX);
+
+		bin = int(floor(binCount*float(rand())/float(RAND_MAX)));
+		sizeX = SizeX[bin];
+		sizeY = SizeY[bin];
+		sizeZ = SizeZ[bin];
+
+//		cout << "Sampling from bin no. " << floor(bin) << endl;
+//		cout << "Feature size is: " << sizeX << "x" << sizeY << "x" << sizeZ << endl;
+
+		for (k=z;k<z+sizeZ;k++){
+			for (j=y;j<y+sizeY;j++){
+				for (i=x;i<x+sizeX;i++){
+					// Identify nodes in the domain (periodic BC)
+					ii = i;
+					jj = j;
+					kk = k;
+					if (ii < 1)			ii+=(Nx-2);
+					if (jj < 1)			jj+=(Ny-2);
+					if (kk < 1)			kk+=(Nz-2);
+					if (!(ii < Nx-1))		ii-=(Nx-2);
+					if (!(jj < Ny-1))		jj-=(Ny-2);
+					if (!(kk < Nz-1))		kk-=(Nz-2);
+
+					n = kk*Nx*Ny+jj*Nx+ii;
+
+					if (ID[n] == 2){
+						ID[n] = 1;
+						count++;
+					}
+				}
+			}
+		}
+		sat = float(count)/total;
+	}
+	//.......................................................................
+}
+
+
+inline void FlipID(char *ID, int N)
+{
+	for (int n=0; n<N; n++){
+		if  	 (ID[n] == 1)	ID[n] = 2;
+		else if  (ID[n] == 2)	ID[n] = 1;
+	}
+}
+
 #define USE_NEW_WRITER
 
 using namespace std;
@@ -487,7 +585,7 @@ int main(int argc, char **argv)
 	id[0] = id[Nx-1] = id[(Ny-1)*Nx] = id[(Ny-1)*Nx + Nx-1] = 0;
 	id[(Nz-1)*Nx*Ny] = id[(Nz-1)*Nx*Ny+Nx-1] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx + Nx-1] = 0;
 	//.........................................................
-	
+/*
 #ifdef USE_EXP_CONTACT_ANGLE
 	// If negative phi_s is chosen, flip the ID for the wetting and non-wetting phase
 	if (phi_s < 0.0 && !pBC){
@@ -508,7 +606,7 @@ int main(int argc, char **argv)
 	}
 
 #endif
-
+*/
 	// Initialize communication structures in averaging domain
 	for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = id[i];
 	Dm.CommInit(MPI_COMM_WORLD);
diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index 6b25ab8f..ba5e6f7d 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -372,9 +372,9 @@ int main(int argc, char **argv)
 	}
 
 	// Generate the residual NWP 
-	if (BoundaryCondition==0 && InitialCondition==0 && rank==0) printf("Initializing with NWP saturation = %f \n",wp_saturation);
-	if (BoundaryCondition==0 && InitialCondition==0) GenerateResidual(id,Nx,Ny,Nz,wp_saturation);
-	if (InitialCondition == 2){
+//	if (BoundaryCondition==0 && InitialCondition==0 && rank==0) printf("Initializing with NWP saturation = %f \n",wp_saturation);
+//	if (BoundaryCondition==0 && InitialCondition==0) GenerateResidual(id,Nx,Ny,Nz,wp_saturation);
+//	if (InitialCondition == 2){
 	    if (rank==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
 	    sprintf(LocalRankFilename,"ID.%05i",rank);
 	    FILE *IDFILE = fopen(LocalRankFilename,"rb");
@@ -382,7 +382,7 @@ int main(int argc, char **argv)
 	    fread(id,1,N,IDFILE);
 	    fclose(IDFILE);
 	    //	    CopyToDevice(ID, id, N);
-	}
+//	}
 
 	
 	// Set up kstart, kfinish so that the reservoirs are excluded from averaging
@@ -438,7 +438,7 @@ int main(int argc, char **argv)
 	id[(Nz-1)*Nx*Ny] = id[(Nz-1)*Nx*Ny+Nx-1] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx + Nx-1] = 0;
 	//.........................................................
 	
-	// If positive phi_s is chosen, flip the ID for the wetting and non-wetting phase
+/*	// If positive phi_s is chosen, flip the ID for the wetting and non-wetting phase
 	if (phi_s > 0.0 && BoundaryCondition==0){
 		phi_s = -phi_s;
 	 	das = (phi_s+1.0)*0.5;
@@ -446,6 +446,7 @@ int main(int argc, char **argv)
 		if (rank == 0)	printf("Resetting phi_s = %f, das = %f, dbs = %f \n", phi_s, das, dbs);
 		FlipID(id,Nx*Ny*Nz);
 	}
+*/
 	// Initialize communication structures in averaging domain
 	for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = id[i];
 	Dm.CommInit(MPI_COMM_WORLD);
diff --git a/tests/lbpm_random_pp.cpp b/tests/lbpm_random_pp.cpp
new file mode 100644
index 00000000..e4a01927
--- /dev/null
+++ b/tests/lbpm_random_pp.cpp
@@ -0,0 +1,218 @@
+/*
+ * Pre-processor to generate signed distance function from segmented data
+ * segmented data should be stored in a raw binary file as 1-byte integer (type char)
+ * will output distance functions for phases
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <Array.h>
+
+inline void GenerateResidual(char *ID, int Nx, int Ny, int Nz, double Saturation)
+{
+	//.......................................................................
+	int i,j,k,n,Number,N;
+	int x,y,z,ii,jj,kk;
+	int sizeX,sizeY,sizeZ;
+	int *SizeX, *SizeY, *SizeZ;
+
+#ifdef NORANDOM
+	srand(10009);
+#else
+	srand(time(NULL));
+#endif
+//	float bin;
+	//.......................................................................
+	N = Nx*Ny*Nz;
+
+	int bin, binCount;
+	ifstream Dist("BlobSize.in");
+	Dist >> binCount;
+//	printf("Number of blob sizes: %i \n",binCount);
+	SizeX = new int [binCount];
+	SizeY = new int [binCount];
+	SizeZ = new int [binCount];
+	for (bin=0; bin<binCount; bin++){
+		Dist >> SizeX[bin];
+		Dist >> SizeY[bin];
+		Dist >> SizeZ[bin];
+	//	printf("Blob %i dimension: %i x %i x %i \n",bin, SizeX[bin], SizeY[bin], SizeZ[bin]);
+	}
+	Dist.close();
+	//.......................................................................
+//	cout << "Generating blocks... " << endl;
+	// Count for the total number of oil nodes
+	int count = 0;
+	// Count the total number of non-solid nodes
+	int total = 0;
+	for (i=0;i<N;i++){
+		if (ID[i] != 0) total++;
+	}
+
+	float sat = 0.f;
+	Number = 0;		// number of features
+	while (sat < Saturation){
+		Number++;
+		// Randomly generate a point in the domain
+		x = Nx*float(rand())/float(RAND_MAX);
+		y = Ny*float(rand())/float(RAND_MAX);
+		z = Nz*float(rand())/float(RAND_MAX);
+
+		bin = int(floor(binCount*float(rand())/float(RAND_MAX)));
+		sizeX = SizeX[bin];
+		sizeY = SizeY[bin];
+		sizeZ = SizeZ[bin];
+
+//		cout << "Sampling from bin no. " << floor(bin) << endl;
+//		cout << "Feature size is: " << sizeX << "x" << sizeY << "x" << sizeZ << endl;
+
+		for (k=z;k<z+sizeZ;k++){
+			for (j=y;j<y+sizeY;j++){
+				for (i=x;i<x+sizeX;i++){
+					// Identify nodes in the domain (periodic BC)
+					ii = i;
+					jj = j;
+					kk = k;
+					if (ii < 1)			ii+=(Nx-2);
+					if (jj < 1)			jj+=(Ny-2);
+					if (kk < 1)			kk+=(Nz-2);
+					if (!(ii < Nx-1))		ii-=(Nx-2);
+					if (!(jj < Ny-1))		jj-=(Ny-2);
+					if (!(kk < Nz-1))		kk-=(Nz-2);
+
+					n = kk*Nx*Ny+jj*Nx+ii;
+
+					if (ID[n] == 2){
+						ID[n] = 1;
+						count++;
+					}
+				}
+			}
+		}
+		sat = float(count)/total;
+	}
+	//.......................................................................
+}
+
+
+inline void FlipID(char *ID, int N)
+{
+	for (int n=0; n<N; n++){
+		if  	 (ID[n] == 1)	ID[n] = 2;
+		else if  (ID[n] == 2)	ID[n] = 1;
+	}
+}
+
+
+int main(int argc, char **argv)
+{
+	// Initialize MPI
+	int rank, nprocs;
+	MPI_Init(&argc,&argv);
+	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+	MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
+
+	int InitialWetting;
+	double Saturation;
+	if (argc == 3){
+		sscanf(argv[1],"%lf",&Saturation);
+		sscanf(argv[2],"%d",&InitialWetting);
+		if (rank==0){
+			printf("Initializing wetting phase saturation of %f \n",Saturation);
+			if (InitialWetting == 1)
+				printf("Begin from connected wetting phase \n");
+			else
+				printf("Begin from connected non-wetting phase \n");
+		}
+	}
+
+	if (InitialWetting == 1)	Saturation=1.0-Saturation;
+    //.......................................................................
+    // Reading the domain information file
+    //.......................................................................
+    int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
+    double Lx, Ly, Lz;
+    int Nx,Ny,Nz;
+    int i,j,k,n;
+	int BC=0;
+
+    if (rank==0){
+    	ifstream domain("Domain.in");
+    	domain >> nprocx;
+    	domain >> nprocy;
+    	domain >> nprocz;
+    	domain >> nx;
+    	domain >> ny;
+    	domain >> nz;
+    	domain >> nspheres;
+    	domain >> Lx;
+    	domain >> Ly;
+    	domain >> Lz;
+
+    }
+	MPI_Barrier(MPI_COMM_WORLD);
+	// Computational domain
+	MPI_Bcast(&nx,1,MPI_INT,0,MPI_COMM_WORLD);
+	MPI_Bcast(&ny,1,MPI_INT,0,MPI_COMM_WORLD);
+	MPI_Bcast(&nz,1,MPI_INT,0,MPI_COMM_WORLD);
+	MPI_Bcast(&nprocx,1,MPI_INT,0,MPI_COMM_WORLD);
+	MPI_Bcast(&nprocy,1,MPI_INT,0,MPI_COMM_WORLD);
+	MPI_Bcast(&nprocz,1,MPI_INT,0,MPI_COMM_WORLD);
+	MPI_Bcast(&nspheres,1,MPI_INT,0,MPI_COMM_WORLD);
+	MPI_Bcast(&Lx,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
+	MPI_Bcast(&Ly,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
+	MPI_Bcast(&Lz,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
+	//.................................................
+	MPI_Barrier(MPI_COMM_WORLD);
+
+    // Check that the number of processors >= the number of ranks
+    if ( rank==0 ) {
+        printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
+        printf("Number of MPI ranks used: %i \n", nprocs);
+        printf("Full domain size: %i x %i x %i  \n",nx*nprocx,ny*nprocy,nz*nprocz);
+    }
+    if ( nprocs < nprocx*nprocy*nprocz ){
+        ERROR("Insufficient number of processors");
+    }
+
+	char LocalRankFilename[40];
+
+	nx+=2; ny+=2; nz+=2;
+	int N = nx*ny*nz;
+    char *id;
+    id = new char[N];
+
+    DoubleArray SignDist(nx,ny,nz);
+    // Read the signed distance from file
+    sprintf(LocalRankFilename,"SignDist.%05i",rank);
+    FILE *DIST = fopen(LocalRankFilename,"rb");
+    fread(SignDist.get(),8,N,DIST);
+    fclose(DIST);
+
+    for (int k=0; k<nz; k++){
+        for (int j=0; j<ny; j++){
+            for (int i=0; i<nx; i++){
+            	if (SignDist(i,j,k) < 0.0)  id[n] = 0;
+            	else 						id[n] = 2;
+            }
+        }
+    }
+
+	// Generate the residual NWP
+	if (rank==0) printf("Initializing with NWP saturation = %f \n",wp_saturation);
+	GenerateResidual(id,Nx,Ny,Nz,wp_saturation);
+
+	if (InitialWetting == 1)	FlipID(id,Nx*Ny*Nz);
+
+    sprintf(LocalRankFilename,"ID.%05i",rank);
+    FILE *ID = fopen(LocalRankFilename,"wb");
+    fread(id,1,N,ID);
+    fclose(ID);
+
+    MPI_Barrier(MPI_COMM_WORLD);
+	MPI_Finalize();
+    return 0;
+}