diff --git a/CHANGELOG.md b/CHANGELOG.md
index 509183e0..96e95bc6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,10 @@
 Version History
 ---------------
 
+### Open VKL 2.0.3
+
+-   Fix artifacts when sampling AMR volume
+
 ### Open VKL 2.0.2
 
 -   Fix used element size in copyDeviceBufferToHost
diff --git a/openvkl/devices/cpu/volume/amr/AMRData.cpp b/openvkl/devices/cpu/volume/amr/AMRData.cpp
index 5468b1dc..4e2f4047 100644
--- a/openvkl/devices/cpu/volume/amr/AMRData.cpp
+++ b/openvkl/devices/cpu/volume/amr/AMRData.cpp
@@ -25,9 +25,7 @@ namespace openvkl {
         this->cellWidth = info.cellWidth;
         this->value     = *ispc(data);
         this->dims      = this->box.size() + vec3i(1);
-        // make it an ULP smaller such that floor(relBrickPos*f_dims) yields
-        // correct index even for relBrickPos=1.0f
-        this->f_dims    = nextafter(this->dims, -1);
+        this->f_dims    = vec3f(this->dims);
 
         this->worldBounds =
             box3f(vec3f(this->box.lower) * this->cellWidth,
diff --git a/openvkl/devices/cpu/volume/amr/CellRef.ispc b/openvkl/devices/cpu/volume/amr/CellRef.ispc
index 67e78cb0..49ebcff9 100644
--- a/openvkl/devices/cpu/volume/amr/CellRef.ispc
+++ b/openvkl/devices/cpu/volume/amr/CellRef.ispc
@@ -28,8 +28,10 @@ extern CellRef findCell(const AMR *uniform self,
         for (uniform int i=0;any(true);i++) {
           const AMRBrick *uniform brick = leaf->brickList[i];
           if (brick->cellWidth >= minWidth) {
-            const vec3f relBrickPos
+            vec3f relBrickPosUnclamped
               = (worldSpacePos - brick->bounds.lower) * brick->bounds_scale;
+            // clamp to 1ULP-less-1 for f_bc stay in valid index range
+            const vec3f relBrickPos = min(relBrickPosUnclamped, make_vec3f(0x1.fffffep-1f));
             // brick coords: integer cell coordinates inside brick
             // OPT: the same calculations as below, just in
             // floats. this works as long as all values we calculate
@@ -75,8 +77,10 @@ extern CellRef findLeafCell(const AMR *uniform self,
       if (isLeaf(node)) {
         const AMRLeaf *uniform leaf = &self->leaf[getOfs(node)];
         const AMRBrick *uniform brick = leaf->brickList[0];
-        const vec3f relBrickPos
+        const vec3f relBrickPosUnclamped
           = (worldSpacePos - brick->bounds.lower) * brick->bounds_scale;
+        // clamp to 1ULP-less-1 for f_bc stay in valid index range
+        const vec3f relBrickPos = min(relBrickPosUnclamped, make_vec3f(0x1.fffffep-1f));
         // brick coords: integer cell coordinates inside brick
         // OPT: the same calculations as below, just in
         // floats. this works as long as all values we calculate
diff --git a/openvkl/devices/cpu/volume/amr/DualCell.ih b/openvkl/devices/cpu/volume/amr/DualCell.ih
index 2302c687..e555d9d9 100644
--- a/openvkl/devices/cpu/volume/amr/DualCell.ih
+++ b/openvkl/devices/cpu/volume/amr/DualCell.ih
@@ -40,20 +40,6 @@ inline vec3f dualCellLerpWeightsForLevel(const vec3f &P, const float cellWidth)
   return (xfmed - f_idx);
 }
 
-inline void initDualCell(DualCell &D,
-                         const vec3f &P,
-                         const uniform AMRLevel &level)
-{
-  const float cellWidth = level.cellWidth;
-  const float halfCellWidth = 0.5f*cellWidth;
-  const float rcpCellWidth = rcp(cellWidth);
-  const vec3f xfmed = (P-halfCellWidth)*rcpCellWidth;
-  const vec3f f_idx = floor(xfmed);
-  D.cellID.pos   = f_idx * cellWidth + halfCellWidth;
-  D.cellID.width = cellWidth;
-  D.weights      = xfmed - f_idx;
-}
-
 inline void initDualCell(DualCell &D, const vec3f &P, const float cellWidth)
 {
   const float halfCellWidth = cellWidth * 0.5f;
@@ -65,6 +51,13 @@ inline void initDualCell(DualCell &D, const vec3f &P, const float cellWidth)
   D.weights      = xfmed - f_idx;
 }
 
+inline void initDualCell(DualCell &D,
+                         const vec3f &P,
+                         const uniform AMRLevel &level)
+{
+  initDualCell(D, P, level.cellWidth);
+}
+
 inline bool allCornersAreLeaves(const DualCell &D)
 {
   return
@@ -91,9 +84,8 @@ inline bool allCornersArePresent(const DualCell &D)
     (D.actualWidth[7] == D.cellID.width);
 }
 
-inline float lerp(const DualCell &D)
+inline float lerpWithExplicitWeights(const DualCell &D, const vec3f &w)
 {
-  const vec3f &w = D.weights;
   const float f000 = D.value[C000];
   const float f001 = D.value[C001];
   const float f010 = D.value[C010];
@@ -115,27 +107,9 @@ inline float lerp(const DualCell &D)
   return f;
 }
 
-inline float lerpWithExplicitWeights(const DualCell &D, const vec3f &w)
+inline float lerp(const DualCell &D)
 {
-  const float f000 = D.value[C000];
-  const float f001 = D.value[C001];
-  const float f010 = D.value[C010];
-  const float f011 = D.value[C011];
-  const float f100 = D.value[C100];
-  const float f101 = D.value[C101];
-  const float f110 = D.value[C110];
-  const float f111 = D.value[C111];
-
-  const float f00 = (1.f-w.x)*f000 + w.x*f001;
-  const float f01 = (1.f-w.x)*f010 + w.x*f011;
-  const float f10 = (1.f-w.x)*f100 + w.x*f101;
-  const float f11 = (1.f-w.x)*f110 + w.x*f111;
-
-  const float f0 = (1.f-w.y)*f00+w.y*f01;
-  const float f1 = (1.f-w.y)*f10+w.y*f11;
-
-  const float f = (1.f-w.z)*f0+w.z*f1;
-  return f;
+  return lerpWithExplicitWeights(D, D.weights);
 }
 
 inline float lerpAlpha(const DualCell &D)
@@ -174,4 +148,4 @@ extern void findDualCell(const AMR *uniform self,
   corner */
 extern void findMirroredDualCell(const AMR *uniform self,
                                  const vec3i &loID,
-                                 DualCell &dual);
\ No newline at end of file
+                                 DualCell &dual);
diff --git a/openvkl/devices/gpu/compute/amr/CellRef.h b/openvkl/devices/gpu/compute/amr/CellRef.h
index a689dddd..a216928b 100644
--- a/openvkl/devices/gpu/compute/amr/CellRef.h
+++ b/openvkl/devices/gpu/compute/amr/CellRef.h
@@ -55,8 +55,11 @@ namespace ispc {
         if (isLeaf(node)) {
           const AMRLeaf *leaf   = &self->leaf[getOfs(node)];
           const AMRBrick *brick = leaf->brickList[0];
-          const vec3f relBrickPos =
+          const vec3f relBrickPosUnclamped =
               (worldSpacePos - brick->bounds.lower) * brick->bounds_scale;
+          // clamp to 1ULP-less-1 for f_bc stay in valid index range
+          const vec3f relBrickPos =
+              min(relBrickPosUnclamped, vec3f(0x1.fffffep-1f));
           // brick coords: integer cell coordinates inside brick
           // OPT: the same calculations as below, just in
           // floats. this works as long as all values we calculate
diff --git a/openvkl/devices/gpu/compute/amr/DualCell.h b/openvkl/devices/gpu/compute/amr/DualCell.h
index fa3e5749..04ed1c16 100644
--- a/openvkl/devices/gpu/compute/amr/DualCell.h
+++ b/openvkl/devices/gpu/compute/amr/DualCell.h
@@ -7,7 +7,7 @@ namespace ispc {
 
   struct DualCellID
   {
-    // position of LOWER-LEFT COORDINATE (ie, CENTER of lower-left cell
+    // position of LOWER-LEFT COORDINATE (ie, CENTER of lower-left cell)
     vec3f pos;
     // width of dual cell, also doubles as level indicator
     float width;
@@ -38,10 +38,9 @@ namespace ispc {
     return (xfmed - f_idx);
   }
 
-  inline void initDualCell(DualCell &D, const vec3f &P, const AMRLevel &level)
+  inline void initDualCell(DualCell &D, const vec3f &P, const float cellWidth)
   {
-    const float cellWidth     = level.cellWidth;
-    const float halfCellWidth = 0.5f * cellWidth;
+    const float halfCellWidth = cellWidth * 0.5f;
     const float rcpCellWidth  = rcp(cellWidth);
     const vec3f xfmed         = (P - halfCellWidth) * rcpCellWidth;
     const vec3f f_idx         = floor(xfmed);
@@ -50,15 +49,9 @@ namespace ispc {
     D.weights                 = xfmed - f_idx;
   }
 
-  inline void initDualCell(DualCell &D, const vec3f &P, const float cellWidth)
+  inline void initDualCell(DualCell &D, const vec3f &P, const AMRLevel &level)
   {
-    const float halfCellWidth = cellWidth * 0.5f;
-    const float rcpCellWidth  = rcp(cellWidth);
-    const vec3f xfmed         = (P - halfCellWidth) * rcpCellWidth;
-    const vec3f f_idx         = floor(xfmed);
-    D.cellID.pos              = f_idx * cellWidth + halfCellWidth;
-    D.cellID.width            = cellWidth;
-    D.weights                 = xfmed - f_idx;
+    initDualCell(D, P, level.cellWidth);
   }
 
   inline bool allCornersAreLeaves(const DualCell &D)
@@ -79,9 +72,8 @@ namespace ispc {
            (D.actualWidth[7] == D.cellID.width);
   }
 
-  inline float lerp(const DualCell &D)
+  inline float lerpWithExplicitWeights(const DualCell &D, const vec3f &w)
   {
-    const vec3f &w   = D.weights;
     const float f000 = D.value[C000];
     const float f001 = D.value[C001];
     const float f010 = D.value[C010];
@@ -103,27 +95,9 @@ namespace ispc {
     return f;
   }
 
-  inline float lerpWithExplicitWeights(const DualCell &D, const vec3f &w)
+  inline float lerp(const DualCell &D)
   {
-    const float f000 = D.value[C000];
-    const float f001 = D.value[C001];
-    const float f010 = D.value[C010];
-    const float f011 = D.value[C011];
-    const float f100 = D.value[C100];
-    const float f101 = D.value[C101];
-    const float f110 = D.value[C110];
-    const float f111 = D.value[C111];
-
-    const float f00 = (1.f - w.x) * f000 + w.x * f001;
-    const float f01 = (1.f - w.x) * f010 + w.x * f011;
-    const float f10 = (1.f - w.x) * f100 + w.x * f101;
-    const float f11 = (1.f - w.x) * f110 + w.x * f111;
-
-    const float f0 = (1.f - w.y) * f00 + w.y * f01;
-    const float f1 = (1.f - w.y) * f10 + w.y * f11;
-
-    const float f = (1.f - w.z) * f0 + w.z * f1;
-    return f;
+    return lerpWithExplicitWeights(D, D.weights);
   }
 
   inline float lerpAlpha(const DualCell &D)