diff --git a/source/blender/nodes/geometry/nodes/node_geo_volume_grid.cc b/source/blender/nodes/geometry/nodes/node_geo_volume_grid.cc
index ee1a2c524e0..6d34ad1ca34 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_volume_grid.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_volume_grid.cc
@@ -2,6 +2,7 @@
 
 #ifdef WITH_OPENVDB
 #  include <openvdb/openvdb.h>
+#  include <openvdb/tools/Dense.h>
 #  include <openvdb/tools/LevelSetUtil.h>
 #  include <openvdb/tools/ParticlesToLevelSet.h>
 #endif
@@ -18,6 +19,8 @@
 #include "BKE_volume.h"
 #include "BKE_volume_to_mesh.hh"
 
+#include "BLI_task.hh"
+
 namespace blender::nodes::node_geo_volume_grid_cc {
 
 static float map(float x, float in_min, float in_max, float out_min, float out_max)
@@ -103,8 +106,10 @@ static void node_geo_exec(GeoNodeExecParams params)
     return;
   }
 
-  if (bounds_min.x == bounds_max.x || bounds_min.y == bounds_max.y || bounds_min.z == bounds_max.z) {
-    params.error_message_add(NodeWarningType::Error, TIP_("Bounding box volume must be greater than 0"));
+  if (bounds_min.x == bounds_max.x || bounds_min.y == bounds_max.y ||
+      bounds_min.z == bounds_max.z) {
+    params.error_message_add(NodeWarningType::Error,
+                             TIP_("Bounding box volume must be greater than 0"));
     params.set_default_remaining_outputs();
     return;
   }
@@ -117,7 +122,7 @@ static void node_geo_exec(GeoNodeExecParams params)
   const int64_t domain_size = static_cast<int64_t>(resX) * static_cast<int64_t>(resY) *
                               static_cast<int64_t>(resZ);
   FieldEvaluator evaluator(context, domain_size);
-  int eval_idx = evaluator.add(input_field);
+  evaluator.add(input_field);
   {
     evaluator.evaluate();
   }
@@ -125,19 +130,11 @@ static void node_geo_exec(GeoNodeExecParams params)
   /* Store resulting values in openvdb grid. */
   openvdb::FloatGrid::Ptr grid = openvdb::FloatGrid::create(FLT_MAX);
   {
-    VArray<float> results = evaluator.get_evaluated(eval_idx).typed<float>();
-    openvdb::FloatGrid::Accessor accessor = grid->getAccessor();
-    openvdb::Coord ijk;
-    int &i = ijk[0], &j = ijk[1], &k = ijk[2];
-    int idx = 0;
-    for (i = 0; i < resX; i++) {
-      for (j = 0; j < resY; j++) {
-        for (k = 0; k < resZ; k++) {
-          accessor.setValue(ijk, results[idx]);
-          idx++;
-        }
-      }
-    }
+    const VArray_Span<float> results = evaluator.get_evaluated<float>(0);
+    openvdb::tools::Dense<float, openvdb::tools::LayoutZYX> dense_grid{
+        openvdb::math::CoordBBox({0, 0, 0}, {resX - 1, resY - 1, resZ - 1}),
+        const_cast<float *>(results.data())};
+    openvdb::tools::copyFromDense(dense_grid, *grid, 0.0f);
   }
 
   float3 scale_fac = (bounds_max - bounds_min) / float3(resX, resY, resZ);