diff --git a/cpu/D3Q19.cpp b/cpu/D3Q19.cpp
index 32ba966f..45d1a4d4 100644
--- a/cpu/D3Q19.cpp
+++ b/cpu/D3Q19.cpp
@@ -1053,7 +1053,7 @@ extern "C" void ScaLBL_D3Q19_Momentum(double *dist, double *vel, int Np)
 	}
 }
 
-extern "C" void ScaLBL_D3Q19_Pressure(const double *fq, double *Pressure, int N)
+extern "C" void ScaLBL_D3Q19_Pressure(double *dist, double *Pressure, int N)
 {
 	int n;
 	// distributions
diff --git a/gpu/D3Q19.cu b/gpu/D3Q19.cu
index f55db1bd..93ac0232 100644
--- a/gpu/D3Q19.cu
+++ b/gpu/D3Q19.cu
@@ -1552,7 +1552,7 @@ __global__  void dvc_ScaLBL_D3Q19_Momentum(double *dist, double *vel, int N)
 	}
 }
 
-__global__  void dvc_ScaLBL_D3Q19_Pressure(const double *fq, double *Pressure, int N)
+__global__  void dvc_ScaLBL_D3Q19_Pressure(const double *dist, double *Pressure, int N)
 {
 	int n;
 	// distributions
@@ -2344,7 +2344,7 @@ extern "C" void ScaLBL_D3Q19_Momentum(double *dist, double *vel, int Np){
 }
 
 extern "C" void ScaLBL_D3Q19_Pressure(double *fq, double *Pressure, int Np){
-	dvc_ScaLBL_D3Q19_Pressure<<< NBLOCKS,NTHREADS >>>(disteven, distodd, Pressure, Nx, Ny, Nz);
+	dvc_ScaLBL_D3Q19_Pressure<<< NBLOCKS,NTHREADS >>>(fq, Pressure, Np);
 }
 
 extern "C" void ScaLBL_D3Q19_Velocity_BC_z(double *disteven, double *distodd, double uz,int Nx, int Ny, int Nz){
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index d54b33db..c9940957 100755
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -19,7 +19,6 @@ ADD_LBPM_EXECUTABLE( lbpm_captube_pp )
 ADD_LBPM_EXECUTABLE( lbpm_inkbottle_pp )
 ADD_LBPM_EXECUTABLE( lbpm_plates_pp )
 ADD_LBPM_EXECUTABLE( lbpm_squaretube_pp )
-ADD_LBPM_EXECUTABLE( lbpm_BlobAnalysis )
 #ADD_LBPM_EXECUTABLE( TestBubble )
 ADD_LBPM_EXECUTABLE( ComponentLabel )
 ADD_LBPM_EXECUTABLE( ColorToBinary )
diff --git a/tests/lbpm_BlobAnalysis.cpp b/tests/lbpm_BlobAnalysis.cpp
deleted file mode 100644
index 37a993ec..00000000
--- a/tests/lbpm_BlobAnalysis.cpp
+++ /dev/null
@@ -1,375 +0,0 @@
-/*
-This code computes TCAT averages on a blob-by-blob basis in parallel
-It requires that the blobs be labeled using BlobIdentify.cpp 
-James E. McClure 2015
- */
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <sys/stat.h>
-#include <iostream>
-#include <exception>
-#include <stdexcept>
-#include <fstream>
-
-#include "common/Domain.h"
-#include "common/TwoPhase.h"
-#include "common/MPI_Helpers.h"
-#include "common/Utilities.h"
-
-inline void ReadBlobFile(char *FILENAME, int *Data, int N)
-{
-    int n;
-    int value;
-    ifstream File(FILENAME,ios::binary);
-    if (File.good()){
-        for (n=0; n<N; n++){
-            // Write the two density values
-            File.read((char*) &value, sizeof(value));
-            Data[n] = value;
-            
-        }
-    }
-    else {
-        for (n=0; n<N; n++) Data[n] = 1.2e-34;
-    }
-    File.close();
-    
-}
-
-inline void  WriteBlobStates(TwoPhase TCAT, double D, double porosity);
-
-int main(int argc, char **argv)
-{
-	//*****************************************
-	// ***** MPI STUFF ****************
-	//*****************************************
-	// Initialize MPI
-	int rank,nprocs;
-	MPI_Init(&argc,&argv);
-    MPI_Comm comm = MPI_COMM_WORLD;
-	MPI_Comm_rank(comm,&rank);
-	MPI_Comm_size(comm,&nprocs);
-	// parallel domain size (# of sub-domains)
-	int nprocx,nprocy,nprocz;
-	int iproc,jproc,kproc;
-	int Nx,Ny,Nz,N,nspheres;
-	double Lx,Ly,Lz;
-
-	int BC=0;	// type of boundary condition applied: 0-periodic, 1-pressure/velocity
-	int nblobs_global=0; 	// number of blobs in the global system
-    
-    // Get the global number of blobs from arguments
-	/*    if (argc > 1){
-        nblobs_global = atoi(argv[1]);
-        if (rank==0)    printf("Number of global blobs is: %i \n",nblobs_global);
-    }
-    else{
-        ERROR("Number of blobs was not specified");
-    }
-	*/
-	
-	int *CubeList;
-
-	if (rank==0){
-		//.......................................................................
-		// 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;
-		//.......................................................................
-	}
-	//.................................................
-	MPI_Barrier(comm);
-	// Computational domain
-	MPI_Bcast(&Nx,1,MPI_INT,0,comm);
-	MPI_Bcast(&Ny,1,MPI_INT,0,comm);
-	MPI_Bcast(&Nz,1,MPI_INT,0,comm);
-	MPI_Bcast(&nprocx,1,MPI_INT,0,comm);
-	MPI_Bcast(&nprocy,1,MPI_INT,0,comm);
-	MPI_Bcast(&nprocz,1,MPI_INT,0,comm);
-	MPI_Bcast(&nspheres,1,MPI_INT,0,comm);
-	MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm);
-	MPI_Barrier(comm);
-	//.................................................
-	Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
-	TwoPhase Averages(Dm);
-	//	BlobTwoPhase ComponentAverages_NWP(nblobs_global);
-	Nx+=2;Ny+=2;Nz+=2;
-	N=Nx*Ny*Nz; // number of lattice points
-	//.......................................................................
-	// Filenames used
-	char LocalRankString[8];
-	char LocalRankFilename[40];
-	char LocalRestartFile[40];
-	char tmpstr[10];
-	sprintf(LocalRankString,"%05d",rank);
-	sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
-	sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
-	//...........................................................................
-	if (rank == 0) cout << "Reading in domain from signed distance function..." << endl;
-	//.......................................................................
-	sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
-	ReadBinaryFile(LocalRankFilename, Averages.SDs.data(), N);
-	MPI_Barrier(comm);
-
-	//	sprintf(LocalRankFilename,"%s%s","Pressure.",LocalRankString);
-	//ReadBinaryFile(LocalRankFilename, Averages.Press.data(), N);
-	//MPI_Barrier(comm);
-	if (rank == 0) cout << "Domain set." << endl;
-    //.......................................................................
-    sprintf(LocalRankFilename,"%s%s","Label_NWP.",LocalRankString);
-    ReadBlobFile(LocalRankFilename, Averages.Label_NWP.data(), N);
-    MPI_Barrier(comm);
-    if (rank == 0) cout << "Label_NWP set." << endl;
-
-    //.......................................................................
-	//copies of data needed to perform checkpointing from cpu
-	double *Den, *DistEven, *DistOdd;
-	Den = new double[2*N];
-	DistEven = new double[10*N];
-	DistOdd = new double[9*N];
-	//.........................................................................
-	if (rank==0) printf("Reading restart file! \n");
-	// Read in the restart file to CPU buffers
-	ReadCheckpoint(LocalRestartFile, Den, DistEven, DistOdd, N);
-	MPI_Barrier(comm);
-	//.........................................................................
-	// Populate the arrays needed to perform averaging
-	if (rank==0) printf("Populate arrays \n");
-    // Compute porosity
-	double porosity,sum,sum_global;
-    sum=0.0;
-    for (int n=0; n<Nx*Ny*Nz; n++){
-        double phi,da,db,press,vx,vy,vz;
-        double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
-        da = Den[n];
-        db = 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 velocity...................................
-        press = 0.3333333333333333*(f0+f2+f1+f4+f3+f6+f5+f8+f7+f10+
-                                          f9+f12+f11+f14+f13+f16+f15+f18+f17);
-        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;
-        Averages.Phase(n)=(da-db)/(da+db);
-        Averages.Phase_tplus(n)=(da-db)/(da+db);
-        Averages.Phase_tminus(n)=(da-db)/(da+db);
-        Averages.Press(n)=press;
-        Averages.Vel_x(n)=vx;
-        Averages.Vel_y(n)=vy;
-        Averages.Vel_z(n)=vz;
-	if (Averages.SDs(n) > 0.0){
-	  Dm.id[n]=1;
-	  sum += 1.0;
-	}
-	else Dm.id[n]=0;
-    }
-    delete [] DistEven;
-    delete [] DistOdd;
-    
-    MPI_Allreduce(&sum,&sum_global,1,MPI_DOUBLE,MPI_SUM,comm);
-    porosity = sum_global/Dm.Volume;
-    if (rank==0) printf("Porosity = %f \n",porosity);
-    Dm.CommInit(comm);
-    for (int i=0; i<N; i++) Averages.SDs(i) -= 1.0; // map the distance 
-    
-    double beta = 0.95;
-     
-    Averages.UpdateSolid();
-    Averages.Initialize();
-    Averages.ComputeDelPhi();
-    Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn);
-    Averages.UpdateMeshValues();
-    Averages.ComponentAverages();
-    Averages.Reduce();
-    int b=0;
-
-    //  Blobs.Set(Averages.ComponentAverages_NWP.NBLOBS);
-    int dimx = Averages.ComponentAverages_NWP.size(0);
-    int dimy = Averages.ComponentAverages_NWP.size(1);
-    int TotalBlobInfoSize=dimx*dimy;
-   
-    FILE *BLOBLOG = NULL;
-    if (rank==0){
-      	BLOBLOG=fopen("blobs.tcat","w");	
-        //printf("dimx=%i \n",dimx);
-    }
-    //      BlobContainer Blobs;
-    DoubleArray RecvBuffer(dimx);
-    //    MPI_Allreduce(&Averages.ComponentAverages_NWP.data(),&Blobs.data(),1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
-    MPI_Barrier(comm);
-    if (rank==0) printf("All ranks passed gate \n");
-
-    for (int b=0; b<(int)Averages.ComponentAverages_NWP.size(1); b++){
-      
-      MPI_Allreduce(&Averages.ComponentAverages_NWP(0,b),&RecvBuffer(0),dimx,MPI_DOUBLE,MPI_SUM,comm);
-      for (int idx=0; idx<dimx-1; idx++) Averages.ComponentAverages_NWP(idx,b)=RecvBuffer(idx);
-      MPI_Barrier(comm);
-      if (Averages.ComponentAverages_NWP(0,b) > 0.0){
-	double Vn,pn,awn,ans,Jwn,Kwn,lwns,cwns,trawn,trJwn;
-        NULL_USE(Vn); NULL_USE(ans); NULL_USE(Jwn);
-	Vn = Averages.ComponentAverages_NWP(1,b);
-	pn = Averages.ComponentAverages_NWP(2,b)/Averages.ComponentAverages_NWP(0,b);
-	awn = Averages.ComponentAverages_NWP(3,b);
-	ans = Averages.ComponentAverages_NWP(4,b);
-	if (awn != 0.0){
-	  Jwn = Averages.ComponentAverages_NWP(5,b)/Averages.ComponentAverages_NWP(3,b);
-	  Kwn = Averages.ComponentAverages_NWP(6,b)/Averages.ComponentAverages_NWP(3,b);
-	}
-	else Jwn=Kwn=0.0;
-	
-	trawn = Averages.ComponentAverages_NWP(12,b);
-	if (trawn != 0.0){
-	  trJwn = Averages.ComponentAverages_NWP(13,b)/trawn;
-	}
-	else trJwn=0.0;
-
-	lwns = Averages.ComponentAverages_NWP(7,b);
-	if (lwns != 0.0) cwns = Averages.ComponentAverages_NWP(8,b)/Averages.ComponentAverages_NWP(7,b);
-	else  cwns=0.0;
-	Averages.ComponentAverages_NWP(2,b) = pn;
-	Averages.ComponentAverages_NWP(5,b) = trJwn;
-	Averages.ComponentAverages_NWP(6,b) = Kwn;
-	Averages.ComponentAverages_NWP(8,b) = cwns;
-	//	Averages.ComponentAverages_NWP(13,b) = trJwn;
-	
-      }
-    }
-
-    Averages.SortBlobs();
-
-      if (rank==0){
-	//	printf("Reduced blob %i \n",b);
-	fprintf(BLOBLOG,"%.5g %.5g %.5g\n",Averages.vol_w_global,Averages.paw_global,Averages.aws_global);
-    for (int b=0; b<(int)Averages.ComponentAverages_NWP.size(1); b++){
-      if (Averages.ComponentAverages_NWP(0,b) > 0.0){
-	double Vn,pn,awn,ans,Jwn,Kwn,lwns,cwns;
-	Vn = Averages.ComponentAverages_NWP(1,b);
-	pn = Averages.ComponentAverages_NWP(2,b);
-	awn = Averages.ComponentAverages_NWP(3,b);
-	ans = Averages.ComponentAverages_NWP(4,b);
-	Jwn = Averages.ComponentAverages_NWP(5,b);
-	Kwn = Averages.ComponentAverages_NWP(6,b);
-	lwns = Averages.ComponentAverages_NWP(7,b);
-	cwns = Averages.ComponentAverages_NWP(8,b);
-	
-      fprintf(BLOBLOG,"%.5g ", Vn); //Vn
-      fprintf(BLOBLOG,"%.5g ", pn); //pn
-      fprintf(BLOBLOG,"%.5g ", awn); //awn
-      fprintf(BLOBLOG,"%.5g ", ans); //ans
-      fprintf(BLOBLOG,"%.5g ", Jwn); //Jwn
-      fprintf(BLOBLOG,"%.5g ", Kwn); //Kwn
-      fprintf(BLOBLOG,"%.5g ", lwns); //lwns  
-      fprintf(BLOBLOG,"%.5g\n",cwns); //cwns
-      }
-    }      
-      }    
-    if (rank==0)  fclose(BLOBLOG);
-
-    double Length=1.0;
-    if (rank==0) WriteBlobStates(Averages,Length,porosity);
-
-    //MPI_Barrier(comm);
-    //printf("Exit, rank=%i \n",rank);
-	// ****************************************************
-	MPI_Barrier(comm);
-	MPI_Finalize();
-	// ****************************************************
-}
-
-inline void  WriteBlobStates(TwoPhase TCAT, double D, double porosity){
-	FILE *BLOBSTATES= fopen("blobstates.tcat","w");
-	int a;
-	double iVol=1.0/TCAT.Dm.Volume;
-	double PoreVolume;
-	double nwp_volume,vol_n,pan,pn,pw,pawn,pwn,awn,ans,aws,Jwn,Kwn,lwns,cwns,clwns;
-	double sw=0,awnD,awsD,ansD,lwnsDD,JwnD,pc;
-	nwp_volume=vol_n=pan=awn=ans=Jwn=Kwn=lwns=clwns=pawn=0.0;
-	pw = TCAT.paw_global;
-	aws = TCAT.aws;
-	// Compute the averages over the entire non-wetting phsae
-	for (a=0; a<(int)TCAT.ComponentAverages_NWP.size(1); a++){
-		vol_n += TCAT.ComponentAverages_NWP(0,a);
-		pan += TCAT.ComponentAverages_NWP(2,a)*TCAT.ComponentAverages_NWP(0,a);
-		awn += TCAT.ComponentAverages_NWP(3,a);
-		ans += TCAT.ComponentAverages_NWP(4,a);
-		Jwn += TCAT.ComponentAverages_NWP(5,a)*TCAT.ComponentAverages_NWP(3,a);
-		Kwn += TCAT.ComponentAverages_NWP(6,a)*TCAT.ComponentAverages_NWP(3,a);
-		lwns += TCAT.ComponentAverages_NWP(7,a);
-		clwns += TCAT.ComponentAverages_NWP(8,a)*TCAT.ComponentAverages_NWP(7,a);
-		nwp_volume += TCAT.ComponentAverages_NWP(1,a);
-		pawn += TCAT.ComponentAverages_NWP(2,a)*TCAT.ComponentAverages_NWP(3,a);
-	}	
-	
-	// Compute the pore voume (sum of wetting an non-wetting phase volumes)
-	PoreVolume=TCAT.wp_volume_global + nwp_volume;
-	// Subtract off portions of non-wetting phase in order of size
-	for (a=TCAT.ComponentAverages_NWP.size(1)-1; a>0; a--){
-		// Subtract the features one-by-one
-		vol_n -= TCAT.ComponentAverages_NWP(0,a);
-		pan -= TCAT.ComponentAverages_NWP(2,a)*TCAT.ComponentAverages_NWP(0,a);
-		awn -= TCAT.ComponentAverages_NWP(3,a);
-		ans -= TCAT.ComponentAverages_NWP(4,a);
-		Jwn -= TCAT.ComponentAverages_NWP(5,a)*TCAT.ComponentAverages_NWP(3,a);
-		Kwn -= TCAT.ComponentAverages_NWP(6,a)*TCAT.ComponentAverages_NWP(3,a);
-		lwns -= TCAT.ComponentAverages_NWP(7,a);
-		clwns -= TCAT.ComponentAverages_NWP(8,a)*TCAT.ComponentAverages_NWP(7,a);
-		nwp_volume -= TCAT.ComponentAverages_NWP(1,a);
-		pawn -= TCAT.ComponentAverages_NWP(2,a)*TCAT.ComponentAverages_NWP(3,a);
-		
-		// Update wetting phase averages
-		aws += TCAT.ComponentAverages_NWP(4,a);
-		if (vol_n > 64){	// Only consider systems with "large enough" blobs -- 4^3
-			if (fabs(1.0 - nwp_volume/PoreVolume - sw) > 0.005 || a == 1){
-				sw = 1.0 - nwp_volume/PoreVolume;
-
-				JwnD = Jwn*D/awn;
-				//trJwnD = -trJwn*D/trawn;
-				cwns = clwns / lwns;
-				pwn = (pawn/awn-pw)*D/0.058;
-				pn = pan/vol_n;
-				awnD = awn*D*iVol;
-				awsD = aws*D*iVol;
-				ansD = ans*D*iVol;
-				lwnsDD = lwns*D*D*iVol;
-				pc = (pn-pw)*D/0.058; 	// hard-coded surface tension due to being lazy
-
-				fprintf(BLOBSTATES,"%.5g %.5g %.5g ",sw,pn,pw);
-				fprintf(BLOBSTATES,"%.5g %.5g %.5g %.5g ",awnD,awsD,ansD,lwnsDD);
-				fprintf(BLOBSTATES,"%.5g %.5g %.5g %.5g %i\n",pc,pwn,JwnD,cwns,a);
-			}
-		}
-	}
-	fclose(BLOBSTATES);
-}