diff --git a/include/2ndOrder.h b/include/2ndOrder.h
index c1ce83b..0ba5ce7 100644
--- a/include/2ndOrder.h
+++ b/include/2ndOrder.h
@@ -84,4 +84,19 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = {
   2.0
 };
 
+//------------------------------------------------------------------------------
+// Dissipation Operator
+//------------------------------------------------------------------------------
+
+#define NUM_BOUNDS_HOD 1
+#define STENCIL_WIDTH_HOD 2
+
+REAL constant D_HO[2*NUM_BOUNDS_HOD+1][STENCIL_WIDTH_HOD] = {
+  // left boundary coefficients
+  { 0.0, 1.0},
+  // central coefficients
+  {-1.0, 1.0},
+  // right boundary coefficients
+  {-1.0, 0.0}
+};
 
diff --git a/include/2ndOrderExtended.h b/include/2ndOrderExtended.h
index 013c4cd..37e8da4 100644
--- a/include/2ndOrderExtended.h
+++ b/include/2ndOrderExtended.h
@@ -91,3 +91,24 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = {
   12.0/5.0
 };
 
+
+//------------------------------------------------------------------------------
+// Dissipation Operator
+//------------------------------------------------------------------------------
+
+#define NUM_BOUNDS_HOD 3
+#define STENCIL_WIDTH_HOD 2
+
+REAL constant D_HO[2*NUM_BOUNDS_HOD+1][STENCIL_WIDTH_HOD] = {
+  // left boundary coefficients
+  { 0.0, 1.0},
+  {-1.0, 1.0},
+  {-1.0, 1.0},
+  // central coefficients
+  {-1.0, 1.0},
+  // right boundary coefficients
+  {-1.0, 1.0},
+  {-1.0, 1.0},
+  {-1.0, 0.0}
+};
+
diff --git a/include/4thOrder.h b/include/4thOrder.h
index a4c402b..61b5aa4 100644
--- a/include/4thOrder.h
+++ b/include/4thOrder.h
@@ -94,3 +94,27 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = {
   48.0/17.0
 };
 
+
+//------------------------------------------------------------------------------
+// Coeffcients of the high-order dissipation operator
+//------------------------------------------------------------------------------
+
+#define STENCIL_WIDTH_HOD 3
+#define NUM_BOUNDS_HOD 4
+
+
+REAL constant D_HO[2*NUM_BOUNDS_HOD+1][STENCIL_WIDTH_HOD] = {
+  // left boundary coefficients
+  { 0.0, 0.0, -1.0},
+  { 0.0, 2.0, -1.0},
+  {-1.0, 2.0, -1.0},
+  {-1.0, 2.0, -1.0},
+  // central coefficients
+  {-1.0, 2.0, -1.0},
+  // right boundary coefficients
+  {-1.0, 2.0, -1.0},
+  {-1.0, 2.0, -1.0},
+  {-1.0, 2.0,  0.0},
+  {-1.0, 0.0,  0.0}
+};
+
diff --git a/include/4thOrderExtended.h b/include/4thOrderExtended.h
index 45112f8..0d09257 100644
--- a/include/4thOrderExtended.h
+++ b/include/4thOrderExtended.h
@@ -100,4 +100,30 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = {
   3.131115459882583
 };
 
+//------------------------------------------------------------------------------
+// Coeffcients of the high-order dissipation operator
+//------------------------------------------------------------------------------
+
+#define STENCIL_WIDTH_HOD 3
+#define NUM_BOUNDS_HOD 6
+
+
+REAL constant D_HO[2*NUM_BOUNDS_HOD+1][STENCIL_WIDTH_HOD] = {
+  // left boundary coefficients
+  { 0.0, 0.0, -1.0},
+  { 0.0, 2.0, -1.0},
+  {-1.0, 2.0, -1.0},
+  {-1.0, 2.0, -1.0},
+  {-1.0, 2.0, -1.0},
+  {-1.0, 2.0, -1.0},
+  // central coefficients
+  {-1.0, 2.0, -1.0},
+  // right boundary coefficients
+  {-1.0, 2.0, -1.0},
+  {-1.0, 2.0, -1.0},
+  {-1.0, 2.0, -1.0},
+  {-1.0, 2.0, -1.0},
+  {-1.0, 2.0,  0.0},
+  {-1.0, 0.0,  0.0}
+};
 
diff --git a/include/6thOrder.h b/include/6thOrder.h
index 46d32fe..7746136 100644
--- a/include/6thOrder.h
+++ b/include/6thOrder.h
@@ -104,4 +104,29 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = {
   3.1650670378782326
 };
 
+//------------------------------------------------------------------------------
+// Coeffcients of the high-order dissipation operator
+//------------------------------------------------------------------------------
+
+#define STENCIL_WIDTH_HOD 4
+#define NUM_BOUNDS_HOD 6
+
+REAL constant D_HO[2*NUM_BOUNDS_HOD+1][STENCIL_WIDTH_HOD] = {
+  // left boundary coefficients
+  {0.0,  0.0, 3.0, -1.0},
+  {0.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  // central coefficients
+  {1.0, -3.0, 3.0, -1.0},
+  // right boundary coefficients
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0,  0.0},
+  {1.0, -3.0, 0.0,  0.0},
+};
 
diff --git a/include/6thOrderExtended.h b/include/6thOrderExtended.h
index 7ea0e33..561dcc8 100644
--- a/include/6thOrderExtended.h
+++ b/include/6thOrderExtended.h
@@ -111,3 +111,34 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = {
 };
 
 
+//------------------------------------------------------------------------------
+// Coeffcients of the high-order dissipation operator
+//------------------------------------------------------------------------------
+
+#define STENCIL_WIDTH_HOD 4
+#define NUM_BOUNDS_HOD 8
+
+
+REAL constant D_HO[2*NUM_BOUNDS_HOD+1][STENCIL_WIDTH_HOD] = {
+  // left boundary coefficients
+  {0.0,  0.0, 3.0, -1.0},
+  {0.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  // central coefficients
+  {1.0, -3.0, 3.0, -1.0},
+  // right boundary coefficients
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0, -1.0},
+  {1.0, -3.0, 3.0,  0.0},
+  {1.0, -3.0, 0.0,  0.0},
+};
+
diff --git a/include_defines/ideal_MHD_defines.h b/include_defines/ideal_MHD_defines.h
index 5365dc1..e2cf417 100644
--- a/include_defines/ideal_MHD_defines.h
+++ b/include_defines/ideal_MHD_defines.h
@@ -64,7 +64,10 @@ TODO: JANHUNEN_CENTRAL, BRACKBILL_BARNES_CENTRAL, cf.
 // Speeds: USE_HLL_SPEED_KUSANO USE_HLL_SPEED_BKW
 #define USE_HLL_SPEED_BKW
 
+// Dissipation in the Volume
 
+#define VOLUME_DISSIPATION_ACTIVE
+#define HO_DISSIPATION_FACTOR {0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01}
 
 enum Fields {
              // Density
diff --git a/kernel/artificial_dissipation.cl b/kernel/artificial_dissipation.cl
new file mode 100644
index 0000000..e82ded0
--- /dev/null
+++ b/kernel/artificial_dissipation.cl
@@ -0,0 +1,93 @@
+// The high order artificial dissipation, ported here from InductionEq
+// TODO: implement per field dissipation factors
+
+#ifndef HO_DISSIPATION_FACTOR
+	define HO_DISSIPATION_FACTOR 0.01
+#endif
+
+// Auxillary function for high order dissipation
+// Used to compute the intermediate result D*a of equation (?)
+void inline diff_HOD_x(uint ix, uint iy, uint iz, int bound_x, int start, global REAL const *d_field, REAL *HOD) {
+
+  REAL um[NUM_TOTAL_VARS] = {0.0};
+  int i, j;
+
+  for (i = 0; i < STENCIL_WIDTH_HOD; i++) {
+    get_field(ix, iy, iz, i + start, 0, 0, d_field, um);
+
+    for (j = 0; i < NUM_CONSERVED_VARS; i++) {
+    	HOD[j] = HOD[j]
+        + M_INV_X[NUM_BOUNDS_HOD + bound_x]
+        * D_HO[NUM_BOUNDS_HOD + bound_x][start + STENCIL_WIDTH_HOD - 1]
+        * D_HO[NUM_BOUNDS_HOD][i] * pown(-1.0, (int)(ORDER / 2 + 1))
+        * um[j];
+    }
+  }
+
+}
+
+void inline diff_HOD_y(uint ix, uint iy, uint iz, int bound_y, int start, global REAL const *d_field, REAL *HOD) {
+
+  REAL um[NUM_TOTAL_VARS] = {0.0};
+  int i, j;
+
+  for (i = 0; i < STENCIL_WIDTH_HOD; i++) {
+    get_field(ix, iy, iz, 0, i + start, 0, d_field, um);
+
+    for (j = 0; j < NUM_CONSERVED_VARS; j++) {
+        HOD[j] = HOD[j]
+        + M_INV_Y[NUM_BOUNDS_HOD + bound_y]
+        * D_HO[NUM_BOUNDS_HOD + bound_y][start + STENCIL_WIDTH_HOD - 1]
+        * D_HO[NUM_BOUNDS_HOD][i] * pown(-1.0, (int)(ORDER / 2 + 1))
+        * um[j];
+    }   
+  }
+  
+}
+
+void inline diff_HOD_z(uint ix, uint iy, uint iz, int bound_z, int start, global REAL const *d_field, REAL *HOD) {
+
+  REAL um[NUM_TOTAL_VARS] = {0.0};
+  int i, j;
+
+  for (i = 0; i < STENCIL_WIDTH_HOD; i++) {
+    get_field(ix, iy, iz, 0,0, i + start, d_field, um);
+
+    for (j = 0; j < NUM_CONSERVED_VARS; j++) {
+        HOD[j] = HOD[j]
+        + M_INV_Y[NUM_BOUNDS_HOD + bound_z]
+        * D_HO[NUM_BOUNDS_HOD + bound_z][start + STENCIL_WIDTH_HOD - 1]
+        * D_HO[NUM_BOUNDS_HOD][i] * pown(-1.0, (int)(ORDER / 2 + 1))
+        * um[j];
+    }
+  }
+
+}
+
+
+
+// Computes the high order dissipation for fields in u at point (ix, iy, iz)
+void inline high_order_dissipation(uint ix, uint iy, uint iz, global REAL const *u, REAL *HOD) {
+  int i, j;
+  uint idx = calc_idx(ix,iy,iz);
+  REAL ho_dissipation_factor[NUM_CONSERVED_VARS] = HO_DISSIPATION_FACTOR; 
+  int bound_x = get_bound_x(ix, NUM_BOUNDS_HOD);
+  int bound_y = get_bound_y(iy, NUM_BOUNDS_HOD);
+  int bound_z = get_bound_z(iz, NUM_BOUNDS_HOD);
+  REAL HOD_X[NUM_TOTAL_VARS] = {0.0};
+  REAL HOD_Y[NUM_TOTAL_VARS] = {0.0};
+  REAL HOD_Z[NUM_TOTAL_VARS] = {0.0};
+
+  for (i = -STENCIL_WIDTH_HOD + 1; i < 1; i++) {
+    diff_HOD_x(ix, iy, iz, bound_x, i, u, HOD_X);
+    diff_HOD_y(ix, iy, iz, bound_y, i, u, HOD_Y);
+    diff_HOD_z(ix, iy, iz, bound_z, i, u, HOD_Z);
+    for(j = 0; j < NUM_CONSERVED_VARS;j++) {
+      HOD[j] = HOD[j] + (REAL)(ho_dissipation_factor[j]) * HOD_X[j]
+                      + (REAL)(ho_dissipation_factor[j]) * HOD_Y[j]
+                      + (REAL)(ho_dissipation_factor[j]) * HOD_Z[j];
+		}
+        }
+
+}
+
diff --git a/kernel/kernel_time_integrator.cl b/kernel/kernel_time_integrator.cl
index 4240095..ca8bb6a 100644
--- a/kernel/kernel_time_integrator.cl
+++ b/kernel/kernel_time_integrator.cl
@@ -85,6 +85,9 @@ inline void compute_du_dt(REAL time, uint ix, uint iy, uint iz, global REAL cons
 
   add_volume_terms(time, ix, iy, iz, u, du_dt);
   add_surface_terms(time, ix, iy, iz, u, du_dt);
+  #ifdef VOLUME_DISSIPATION_ACTIVE
+    high_order_dissipation(ix, iy, iz, u, du_dt);
+  #endif // VOLUME_DISSIPATION_ACTIVE
 }
 
 //--------------------------------------------------------------------------------------------------