diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9311c973eb4bbbba562438169cc5934f1068b427..82376fedcbcb63704ede2a86f723dbb0e4fc44c6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -40,7 +40,7 @@
 * Coupling between transient models for water flow and solute transport !96
 * Initial conditions generated from H5 input data !130
 * Generic Python VTK file reader !143
-* Linear interpolator for initial conditions and scaling fields !145
+* Linear interpolator for initial conditions and scaling fields !145, !156
 * Parameterizations for hydrodynamic dispersion in solute transport !141
 * Generic Python VTK file reader !143, !150
 
diff --git a/dune/dorie/common/interpolator.hh b/dune/dorie/common/interpolator.hh
index 6ff3aa1037091374e4b1327bd0e0a328ebf6a4d7..7441108ffa5fb31ed32c4f649497eedaa6b167d4 100644
--- a/dune/dorie/common/interpolator.hh
+++ b/dune/dorie/common/interpolator.hh
@@ -267,7 +267,7 @@ public:
 
         // compute indices and position vectors in mesh units
         const auto pos_unit = get_pos_unit(pos);
-        const auto [idx_lower, idx_upper] = get_knot_indices(pos);
+        const auto [idx_lower, idx_upper] = get_knot_indices(pos_unit);
         const auto pos_unit_diff = get_pos_unit_diff(pos_unit, idx_lower);
         const auto& data = this->_data;
 
diff --git a/dune/dorie/test/test-interpolators.cc b/dune/dorie/test/test-interpolators.cc
index c34e8c38e3f79905720aa4441142edcce95d5ae0..82724780264665541b50207373274358ecca1913 100644
--- a/dune/dorie/test/test-interpolators.cc
+++ b/dune/dorie/test/test-interpolators.cc
@@ -138,10 +138,12 @@ template<>
 void test_linear<2> ()
 {
     constexpr int dim = 2;
-    std::vector<double> data({0.0, 1.0, 1.0, 2.0});
+    std::vector<double> data({0.0, 0.0, 1.0,
+                              1.0, 1.0, 2.0,
+                              3.0, 3.0, 4.0});
 
     std::vector<size_t> shape(dim);
-    std::fill(begin(shape), end(shape), 2);
+    std::fill(begin(shape), end(shape), 3);
 
     Dune::FieldVector<double, dim> extensions(1.0);
     Dune::FieldVector<double, dim> offset(0.0);
@@ -156,22 +158,25 @@ void test_linear<2> ()
     
     // check without offset
     using Dune::FloatCmp::eq; // floating-point comparison
-    std::vector<Dune::FieldVector<double, dim>> points(5);
+    std::vector<Dune::FieldVector<double, dim>> points(7);
 
     points[0] = {0.0, 0.0};
     points[1] = {1.0, 0.0};
     points[2] = {0.0, 1.0};
     points[3] = {1.0, 1.0};
-    points[4] = {0.75, 0.75};
+    points[4] = {0.5, 0.5};
+    points[5] = {0.25, 0.5};
+    points[6] = {1.0, 0.75};
     assert(eq(interp->evaluate(points[0]), 0.0));
     assert(eq(interp->evaluate(points[1]), 1.0));
-    assert(eq(interp->evaluate(points[2]), 1.0));
-    assert(eq(interp->evaluate(points[3]), 2.0));
-    assert(eq(interp->evaluate(points[4]), 1.5));
-
+    assert(eq(interp->evaluate(points[2]), 3.0));
+    assert(eq(interp->evaluate(points[3]), 4.0));
+    assert(eq(interp->evaluate(points[4]), 1.0));
+    assert(eq(interp->evaluate(points[5]), 1.0));
+    assert(eq(interp->evaluate(points[6]), 3.0));
 
     // check with offset
-    std::fill(offset.begin(), offset.end(), -1.0);
+    std::fill(offset.begin(), offset.end(), -0.5);
     interp = Dune::Dorie::InterpolatorFactory<InterpolatorTraits<dim>>
         ::create("linear",
                  data,
@@ -179,11 +184,8 @@ void test_linear<2> ()
                  extensions,
                  offset);
 
-    points.resize(2);
-    points[0] = {0.0, 0.0};
-    points[1] = {-0.25, -0.25};
-    assert(eq(interp->evaluate(points[0]), 2.0));
-    assert(eq(interp->evaluate(points[1]), 1.5));
+    assert(eq(interp->evaluate(points[0]), 1.0));
+    assert(eq(interp->evaluate(points[5]), 3.5));
 }
 
 /// Test the linear interpolator in 3D
@@ -197,7 +199,7 @@ void test_linear<3> ()
     std::vector<size_t> shape(dim);
     std::fill(begin(shape), end(shape), 2);
 
-    Dune::FieldVector<double, dim> extensions(1.0);
+    Dune::FieldVector<double, dim> extensions(2.0);
     Dune::FieldVector<double, dim> offset(0.0);
 
     // build interpolator
@@ -213,10 +215,10 @@ void test_linear<3> ()
     std::vector<Dune::FieldVector<double, dim>> points(5);
 
     points[0] = {0.0, 0.0, 0.0};
-    points[1] = {1.0, 0.0, 0.0};
-    points[2] = {0.0, 0.0, 1.0};
-    points[3] = {1.0, 1.0, 1.0};
-    points[4] = {0.5, 0.5, 1.0};
+    points[1] = {2.0, 0.0, 0.0};
+    points[2] = {0.0, 0.0, 2.0};
+    points[3] = {2.0, 2.0, 2.0};
+    points[4] = {1.0, 1.0, 2.0};
     assert(eq(interp->evaluate(points[0]), 0.0));
     assert(eq(interp->evaluate(points[1]), 0.0));
     assert(eq(interp->evaluate(points[2]), 1.0));
@@ -224,7 +226,7 @@ void test_linear<3> ()
     assert(eq(interp->evaluate(points[4]), 1.5));
 
     // check with offset
-    std::fill(offset.begin(), offset.end(), -1.0);
+    std::fill(offset.begin(), offset.end(), -2.0);
     interp = Dune::Dorie::InterpolatorFactory<InterpolatorTraits<dim>>
         ::create("linear",
                  data,
@@ -234,7 +236,7 @@ void test_linear<3> ()
 
     points.resize(2);
     points[0] = {0.0, 0.0, 0.0};
-    points[1] = {-0.5, -0.5, 0.0};
+    points[1] = {-1.0, -1.0, 0.0};
     assert(eq(interp->evaluate(points[0]), 2.0));
     assert(eq(interp->evaluate(points[1]), 1.5));
 }