diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index bdbf5b0e..6d914092 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -40,6 +40,8 @@ void ScaLBL_ColorModel::ReadParams(string filename){
 	inletB=0.f;
 	outletA=0.f;
 	outletB=1.f;
+	
+	if (BoundaryCondition==4) flux = din*rhoA; // mass flux must adjust for density (see formulation for details)
 
 	// Read domain parameters
 	auto L = domain_db->getVector<double>( "L" );
@@ -56,19 +58,14 @@ void ScaLBL_ColorModel::ReadParams(string filename){
 	nprocy = nproc[1];
 	nprocz = nproc[2];
 
-	if (BoundaryCondition==4) flux = din*rhoA; // mass flux must adjust for density (see formulation for details)
-
+}
+void ScaLBL_ColorModel::SetDomain(){
 	Dm  = std::shared_ptr<Domain>(new Domain(domain_db,comm));      // full domain for analysis
 	Mask  = std::shared_ptr<Domain>(new Domain(domain_db,comm));    // mask domain removes immobile phases
-
 	Nx+=2; Ny+=2; Nz += 2;
 	N = Nx*Ny*Nz;
 	for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1;               // initialize this way
 	Averages = std::shared_ptr<TwoPhase> ( new TwoPhase(Dm) ); // TwoPhase analysis object
-
-	// local copy of the ids
-	id = new char[N];
-
 	MPI_Barrier(comm);
 	Dm->CommInit();
 	MPI_Barrier(comm);
@@ -80,7 +77,7 @@ void ScaLBL_ColorModel::ReadInput(){
     //.......................................................................
     if (rank == 0)    printf("Read input media... \n");
     //.......................................................................
-    Dm->ReadIDs();
+    Mask->ReadIDs();
     
     sprintf(LocalRankString,"%05d",Dm->rank());
     sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
@@ -194,7 +191,7 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
 		for (int j=0;j<Ny;j++){
 			for (int i=0;i<Nx;i++){
 				int n = k*Nx*Ny+j*Nx+i;
-				VALUE=Dm->id[n];
+				VALUE=Mask->id[n];
 				// Assign the affinity from the paired list
 				for (int idx=0; idx < NLABELS; idx++){
 					//printf("rank=%i, idx=%i, value=%i, %i, \n",rank(),idx, VALUE,LabelList[idx]);
@@ -222,7 +219,7 @@ void ScaLBL_ColorModel::Create(){
 	
 	//.........................................................
 	// Initialize communication structures in averaging domain
-	for (int i=0; i<Nx*Ny*Nz; i++) Mask->id[i] = Dm->id[i];
+	for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i];
 	Mask->CommInit();
 	Np=Mask->PoreCount();
 	//...........................................................................
@@ -285,13 +282,7 @@ void ScaLBL_ColorModel::Create(){
 /********************************************************
  * AssignComponentLabels                                 *
  ********************************************************/
-
-void ScaLBL_ColorModel::Initialize(){
-	/*
-	 * This function initializes model (includes both mobile and immobile components)
-	 */
-       int rank=Dm->rank();
-	// Compute the solid interaction potential and copy result to device
+void ScaLBL_ColorModel::SolidPotential(){
 	if (rank==0) printf("Computing solid interaction potential \n");
 	double *PhaseLabel;
 	PhaseLabel=new double [Nx*Ny*Nz];
@@ -390,8 +381,15 @@ void ScaLBL_ColorModel::Initialize(){
 	PFILE = fopen(LocalRankFilename,"wb");
 	fwrite(Psnorm.data(),8,N,PFILE);
 	fclose(PFILE);
-
+}
+void ScaLBL_ColorModel::Initialize(){
+	/*
+	 * This function initializes model
+	 */
+       int rank=Dm->rank();
 	double count_wet=0.f;
+	double *PhaseLabel;
+	PhaseLabel=new double [Nx*Ny*Nz];
 	for (int k=1; k<Nz-1; k++){
 		for (int j=1; j<Ny-1; j++){
 			for (int i=1; i<Nx-1; i++){
diff --git a/models/ColorModel.h b/models/ColorModel.h
index e045972a..67f7a21d 100644
--- a/models/ColorModel.h
+++ b/models/ColorModel.h
@@ -24,10 +24,11 @@ public:
 	// functions in they should be run
 	void ReadParams(string filename);
 	void ReadParams(std::shared_ptr<Database> db0);
-    void AssignComponentLabels(double *phase);
+	void SetDomain();
 	void ReadInput();
 	void Create();
 	void Initialize();
+	void SolidPotential();
 	void Run();
 	void WriteDebug();
 	
@@ -74,6 +75,7 @@ private:
    
     //int rank,nprocs;
     void LoadParams(std::shared_ptr<Database> db0);
+    void AssignComponentLabels(double *phase);
     	
 };
 
diff --git a/tests/TestColorBubble.cpp b/tests/TestColorBubble.cpp
index 96313ac7..ffbd38d3 100644
--- a/tests/TestColorBubble.cpp
+++ b/tests/TestColorBubble.cpp
@@ -14,8 +14,6 @@ using namespace std;
 inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius){
 	// initialize a bubble
 	int i,j,k,n;
-	int Np=0;
-	double sum=0.f;
 	int rank = ColorModel.Mask->rank();
 	int nprocx = ColorModel.Mask->nprocx();
 	int nprocy = ColorModel.Mask->nprocy();
@@ -24,21 +22,11 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
 	int Ny = ColorModel.Mask->Ny;
 	int Nz = ColorModel.Mask->Nz;
 	if (rank == 0) cout << "Setting up bubble..." << endl;
-	sum=0; Np=0;
 	for (k=0;k<Nz;k++){
 		for (j=0;j<Ny;j++){
 			for (i=0;i<Nx;i++){
 				n = k*Nx*Ny + j*Nz + i;
 				ColorModel.Averages->SDs(i,j,k) = 100.f;
-				// Initialize phase positions field
-				if (ColorModel.Averages->SDs(i,j,k) < 0.0){
-					ColorModel.id[n] = 0;
-					ColorModel.id[n] = 0;
-				}
-				else {
-					sum++;
-					Np++;
-				}
 			}
 		}
 	}
@@ -47,19 +35,19 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
 		for (j=0;j<Ny;j++){
 			for (i=0;i<Nx;i++){
 				n = k*Nx*Ny + j*Nz + i;
-				int iglobal= i+(Nx-2)*ColorModel.Dm->iproc();
-				int jglobal= j+(Ny-2)*ColorModel.Dm->jproc();
-				int kglobal= k+(Nz-2)*ColorModel.Dm->kproc();
+				int iglobal= i+(Nx-2)*ColorModel.Mask->iproc();
+				int jglobal= j+(Ny-2)*ColorModel.Mask->jproc();
+				int kglobal= k+(Nz-2)*ColorModel.Mask->kproc();
 				// Initialize phase position field for parallel bubble test
 				if ((iglobal-0.5*(Nx-2)*nprocx)*(iglobal-0.5*(Nx-2)*nprocx)
 						+(jglobal-0.5*(Ny-2)*nprocy)*(jglobal-0.5*(Ny-2)*nprocy)
 						+(kglobal-0.5*(Nz-2)*nprocz)*(kglobal-0.5*(Nz-2)*nprocz) < BubbleRadius*BubbleRadius){
-					ColorModel.id[n] = 2;
-					ColorModel.id[n] = 2;
+					ColorModel.Mask->id[n] = 2;
+					ColorModel.Mask->id[n] = 2;
 				}
 				else{
-					ColorModel.id[n]=1;
-					ColorModel.id[n]=1;
+					ColorModel.Mask->id[n]=1;
+					ColorModel.Mask->id[n]=1;
 				}
 			}
 		}
@@ -92,12 +80,13 @@ int main(int argc, char **argv)
 		auto filename = argv[1];
 		ScaLBL_ColorModel ColorModel(rank,nprocs,comm);
 		ColorModel.ReadParams(filename);
-		//ColorModel.ReadInput();
+		ColorModel.SetDomain();    
+		//ColorModel.ReadInput(); 
 		double radius=15.5;
 		InitializeBubble(ColorModel,radius);
-		ColorModel.Create();
-		ColorModel.Initialize();
-		ColorModel.Run();
+		ColorModel.Create();       // creating the model will create data structure to match the pore structure and allocate variables
+		ColorModel.Initialize();   // initializing the model will set initial conditions for variables
+		ColorModel.Run();	       
 		ColorModel.WriteDebug();
 	}
 	// ****************************************************