diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index fb98482..bb7136a 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -33,3 +33,4 @@ configure_file(test_blast_large.ini test_blast_large.ini COPYONLY)
 configure_file(test_implode.ini test_implode.ini COPYONLY)
 configure_file(test_implode_big.ini test_implode_big.ini COPYONLY)
 configure_file(test_four_quadrant.ini test_four_quadrant.ini COPYONLY)
+configure_file(test_discontinuity.ini test_discontinuity.ini COPYONLY)
diff --git a/src/HydroParams.cpp b/src/HydroParams.cpp
index 079375f..0cb617f 100644
--- a/src/HydroParams.cpp
+++ b/src/HydroParams.cpp
@@ -94,6 +94,10 @@ HydroParams::setup(ConfigMap & configMap)
   {
     problemType = PROBLEM_FOUR_QUADRANT;
   }
+  else if (!problemStr.compare("discontinuity"))
+  {
+    problemType = PROBLEM_DISCONTINUITY;
+  }
   else
   {
     std::cout << "Problem is invalid\n";
diff --git a/src/HydroParams.h b/src/HydroParams.h
index 13ecf7a..1bb6580 100644
--- a/src/HydroParams.h
+++ b/src/HydroParams.h
@@ -91,7 +91,8 @@ enum ProblemType
 {
   PROBLEM_IMPLODE,
   PROBLEM_BLAST,
-  PROBLEM_FOUR_QUADRANT
+  PROBLEM_FOUR_QUADRANT,
+  PROBLEM_DISCONTINUITY
 };
 
 // variable names in the order as in component index
diff --git a/src/HydroRun.h b/src/HydroRun.h
index f9897a5..1f54cf6 100644
--- a/src/HydroRun.h
+++ b/src/HydroRun.h
@@ -98,6 +98,8 @@ class HydroRun
   init_blast(DataArray_t Udata);
   void
   init_four_quadrant(DataArray_t Udata);
+  void
+  init_discontinuity(DataArray_t Udata);
 
   // host routines (save data to file, device data are copied into host
   // inside this routine)
@@ -181,6 +183,10 @@ HydroRun<device_t>::HydroRun(HydroParams & params, ConfigMap & configMap)
   {
     init_four_quadrant(U);
   }
+  else if (params.problemType == PROBLEM_DISCONTINUITY)
+  {
+    init_discontinuity(U);
+  }
   else
   {
     std::cout << "Problem : " << params.problemType
@@ -391,6 +397,21 @@ HydroRun<device_t>::init_blast(DataArray_t Udata)
 
 } // HydroRun<device_t>::init_blast
 
+// =======================================================
+// =======================================================
+/**
+ * Hydrodynamical discontinuity Test.
+ *
+ */
+template <typename device_t>
+void
+HydroRun<device_t>::init_discontinuity(DataArray_t Udata)
+{
+
+  InitDiscontinuityFunctor<device_t>::apply(params, Udata);
+
+} // HydroRun<device_t>::init_discontinuity
+
 // =======================================================
 // =======================================================
 /**
diff --git a/src/HydroRunFunctors.h b/src/HydroRunFunctors.h
index 23c740a..1155619 100644
--- a/src/HydroRunFunctors.h
+++ b/src/HydroRunFunctors.h
@@ -1329,6 +1329,70 @@ class InitFourQuadrantFunctor : public HydroBaseFunctor
 
 }; // InitFourQuadrantFunctor
 
+/*************************************************/
+/*************************************************/
+/*************************************************/
+template <typename device_t>
+class InitDiscontinuityFunctor : public HydroBaseFunctor
+{
+
+public:
+  using DataArray_t = DataArray<device_t>;
+  using exec_space = typename device_t::execution_space;
+
+  InitDiscontinuityFunctor(HydroParams params, DataArray_t Udata)
+    : HydroBaseFunctor(params)
+    , Udata(Udata){};
+
+  // static method which does it all: create and execute functor
+  static void
+  apply(HydroParams params, DataArray_t Udata)
+  {
+    InitDiscontinuityFunctor functor(params, Udata);
+    Kokkos::parallel_for(
+      "InitDiscontinuity",
+      Kokkos::MDRangePolicy<exec_space, Kokkos::Rank<2>>({ 0, 0 }, { params.isize, params.jsize }),
+      functor);
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  void
+  operator()(const int & i, const int & j) const
+  {
+
+    const int ghostWidth = params.ghostWidth;
+
+    const real_t xmin = params.xmin;
+    const real_t ymin = params.ymin;
+    const real_t dx = params.dx;
+    const real_t dy = params.dy;
+
+    const real_t gamma0 = params.settings.gamma0;
+
+    real_t x = xmin + dx / 2 + (i - ghostWidth) * dx;
+    real_t y = ymin + dy / 2 + (j - ghostWidth) * dy;
+
+    if (x + y < 1)
+    {
+      Udata(i, j, ID) = 1.0 + x * x;
+      Udata(i, j, IP) = 1.0 / (gamma0 - 1.0);
+      Udata(i, j, IU) = 0.0;
+      Udata(i, j, IV) = 0.0;
+    }
+    else
+    {
+      Udata(i, j, ID) = 0.25;
+      Udata(i, j, IP) = 1.0 / (gamma0 - 1.0);
+      Udata(i, j, IU) = 0.0;
+      Udata(i, j, IV) = 0.0;
+    }
+
+  } // end operator ()
+
+  DataArray_t Udata;
+
+}; // InitDiscontinuityFunctor
+
 /*************************************************/
 /*************************************************/
 /*************************************************/
diff --git a/src/test_discontinuity.ini b/src/test_discontinuity.ini
new file mode 100644
index 0000000..7425a9d
--- /dev/null
+++ b/src/test_discontinuity.ini
@@ -0,0 +1,33 @@
+[run]
+tEnd=0.6
+nStepmax=1200
+nOutput=100
+
+[mesh]
+nx=256
+ny=256
+boundary_type_xmin=2
+boundary_type_xmax=2
+boundary_type_ymin=2
+boundary_type_ymax=2
+
+xmin=0.0
+xmax=1.0
+
+ymin=0.0
+ymax=1.0
+
+[hydro]
+gamma0=1.666
+cfl=0.8
+niter_riemann=10
+iorder=2
+slope_type=2
+problem=discontinuity
+riemann=hllc
+
+[output]
+outputPrefix=test_discontinuity
+
+[other]
+implementationVersion=0