diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index b4ffb285..42ea8d7c 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -521,8 +521,10 @@ void SubPhase::Full(){
 				double vz = 0.5*(Vel_y(i,j,k+1) - Vel_y(i,j,k-1));
 				double wx = 0.5*(Vel_z(i+1,j,k) - Vel_z(i-1,j,k));
 				double wy = 0.5*(Vel_z(i,j+1,k) - Vel_z(i,j-1,k));
-				double wz = 0.5*(Vel_z(i,j,k+1) - Vel_z(i,j,k-1));				
-				Dissipation(i,j,k) = 2*rho*nu*( ux*ux + vy*vy + wz*wz + 0.5*(vx + uy)*(vx + uy)+ 0.5*(vz + wy)*(vz + wy)+ 0.5*(uz + wx)*(uz + wx));
+				double wz = 0.5*(Vel_z(i,j,k+1) - Vel_z(i,j,k-1));
+				if (SDs(i,j,k) > 2.0){
+					Dissipation(i,j,k) = 2*rho*nu*( ux*ux + vy*vy + wz*wz + 0.5*(vx + uy)*(vx + uy)+ 0.5*(vz + wy)*(vz + wy)+ 0.5*(uz + wx)*(uz + wx));					
+				}
 			}
 		}
 	}
diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index 9f7ff428..86abc7a9 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -303,13 +303,13 @@ public:
         }
 
         if ( vis_db->getWithDefault<bool>( "save_distance", false ) ) {
-            ASSERT( visData[0].vars[5]->name == "SignDist" );
+            ASSERT( visData[0].vars[6]->name == "SignDist" );
             Array<double> &SignData = visData[0].vars[6]->data;
             fillData.copy( Averages.SDs, SignData );
         }
 
         if ( vis_db->getWithDefault<bool>( "save_connected_components", false ) ) {
-            ASSERT( visData[0].vars[6]->name == "BlobID" );
+            ASSERT( visData[0].vars[7]->name == "BlobID" );
             Array<double> &BlobData = visData[0].vars[7]->data;
             fillData.copy( Averages.morph_n->label, BlobData );
         }
@@ -658,6 +658,7 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> input_db, const RankInfoStru
     auto VxVar       = std::make_shared<IO::Variable>();
     auto VyVar       = std::make_shared<IO::Variable>();
     auto VzVar       = std::make_shared<IO::Variable>();
+    auto ViscousDissipationVar = std::make_shared<IO::Variable>();
     auto SignDistVar = std::make_shared<IO::Variable>();
     auto BlobIDVar   = std::make_shared<IO::Variable>();
 
@@ -694,6 +695,14 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> input_db, const RankInfoStru
         VzVar->data.resize( d_n[0], d_n[1], d_n[2] );
         d_meshData[0].vars.push_back( VzVar );
     }
+    
+    if ( vis_db->getWithDefault<bool>( "save_dissipation", false ) ) {
+        ViscousDissipationVar->name = "ViscousDissipation";
+    	ViscousDissipationVar->type = IO::VariableType::VolumeVariable;
+    	ViscousDissipationVar->dim  = 1;
+    	ViscousDissipationVar->data.resize( d_n[0], d_n[1], d_n[2] );
+        d_meshData[0].vars.push_back( ViscousDissipationVar );
+    }
 
     if ( vis_db->getWithDefault<bool>( "save_distance", false ) ) {
         SignDistVar->name = "SignDist";
diff --git a/example/DropletCoalescence/Droplets.py b/example/DropletCoalescence/Droplets.py
new file mode 100644
index 00000000..aa8e33d0
--- /dev/null
+++ b/example/DropletCoalescence/Droplets.py
@@ -0,0 +1,35 @@
+import numpy
+import math
+
+nx=400
+ny=200
+nz=200
+N=nx*ny*nz
+Radius=64
+
+mesh=(nx,ny,nz)
+data=numpy.ones(mesh,dtype=numpy.int8)
+
+#print(data)
+print("Create two droplets")
+print("Mesh size: "+repr(mesh))
+print("Droplet radius: "+repr(Radius))
+
+gap = 6
+c1x = nx/2 - gap/2 - Radius
+c2x = nx/2 + gap/2 + Radius
+
+# assign a bubble on each side
+for x in range(0,200):
+  for y in range(0,ny):
+    for z in range(0,nz):
+       if math.sqrt((x-c1x)*(x-c1x)+(y-ny/2)*(y-ny/2)+(z-nz/2)*(z-nz/2) ) < Radius:
+           data[x,y,z]=2
+
+for x in range(200,nx):
+  for y in range(0,ny):
+    for z in range(0,nz):
+       if math.sqrt((x-c2x)*(x-c2x)+(y-ny/2)*(y-ny/2)+(z-nz/2)*(z-nz/2) ) < Radius:
+           data[x,y,z]=2
+
+data.tofile("Droplets.raw")
diff --git a/example/DropletCoalescence/input.db b/example/DropletCoalescence/input.db
new file mode 100644
index 00000000..ba74be4e
--- /dev/null
+++ b/example/DropletCoalescence/input.db
@@ -0,0 +1,51 @@
+MRT {
+    timestepMax = 10000
+    tau = 0.7
+    F = 1e-05, 0, 0
+    Restart = false
+    din = 1.0
+    dout = 1.0
+    flux = 0.0
+}
+
+Color {
+    tauA = 0.7;  
+    tauB = 1.0;  
+    rhoA   = 1.0;  
+    rhoB   = 1.0;  
+    alpha = 5e-3;
+    beta  = 0.95;
+    F = 0, 0, 0
+    Restart = false
+    timestepMax = 40000
+    ComponentLabels = -2
+    ComponentAffinity = -0.5
+}
+
+Domain {
+    Filename = "Droplets.raw"
+    nproc = 1, 1, 2     // Number of processors (Npx,Npy,Npz)
+    n = 200, 200, 200      // Size of local domain (Nx,Ny,Nz)
+    N = 200, 200, 400      // size of the input image
+    voxel_length = 1.0 
+    BC = 0              // Boundary condition type
+    Sw = 0.15
+    ReadType = "8bit"
+    ReadValues = -2, 1, 2
+    WriteValues = -2, 1, 2
+    ComponentLabels = -2
+    HistoryLabels = -2
+}
+
+Analysis {
+    analysis_interval = 1000    // Frequency to perform analysis
+    subphase_analysis_interval = 5000    // Frequency to perform analysis
+    restart_interval = 60000    // Frequency to write restart data
+    visualization_interval = 100000        // 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"
+}
+
+Visualization {
+}