diff --git a/intern/cycles/blender/image.cpp b/intern/cycles/blender/image.cpp
index 1777341ecd1..a9b22570e1a 100644
--- a/intern/cycles/blender/image.cpp
+++ b/intern/cycles/blender/image.cpp
@@ -14,6 +14,7 @@
 #include "blender/session.h"
 
 #include "util/half.h"
+#include "util/transform.h"
 #include "util/types_float4.h"
 
 CCL_NAMESPACE_BEGIN
@@ -64,8 +65,6 @@ bool BlenderImageLoader::load_metadata(const ImageDeviceFeatures & /*features*/,
     BKE_image_release_ibuf(b_image, ibuf, lock);
   }
 
-  metadata.depth = 1;
-
   if (is_float) {
     if (metadata.channels == 1) {
       metadata.type = IMAGE_DATA_TYPE_FLOAT;
diff --git a/intern/cycles/blender/image.h b/intern/cycles/blender/image.h
index a52e6bc4838..44fd8e8a1ad 100644
--- a/intern/cycles/blender/image.h
+++ b/intern/cycles/blender/image.h
@@ -9,6 +9,7 @@
 #include "RNA_blender_cpp.hh"
 
 #include "scene/image.h"
+#include "scene/image_vdb.h"
 
 struct Image;
 struct ImageUser;
diff --git a/intern/cycles/blender/volume.cpp b/intern/cycles/blender/volume.cpp
index 43e61e6d46d..7ef83643a5c 100644
--- a/intern/cycles/blender/volume.cpp
+++ b/intern/cycles/blender/volume.cpp
@@ -11,43 +11,44 @@
 #include "blender/util.h"
 
 #include "util/log.h"
+#include "util/vector.h"
 
 #include "BKE_volume_grid.hh"
 
 CCL_NAMESPACE_BEGIN
 
 /* TODO: verify this is not loading unnecessary attributes. */
-class BlenderSmokeLoader : public ImageLoader {
+class BlenderSmokeLoader : public VDBImageLoader {
  public:
-  BlenderSmokeLoader(BL::Object &b_ob, AttributeStandard attribute)
-      : b_domain(object_fluid_gas_domain_find(b_ob)), attribute(attribute)
+  BlenderSmokeLoader(BL::Object &b_ob, AttributeStandard attribute, const float clipping)
+      : VDBImageLoader(Attribute::standard_name(attribute), clipping),
+        b_domain(object_fluid_gas_domain_find(b_ob)),
+        attribute(attribute)
   {
     mesh_texture_space(
         *static_cast<const ::Mesh *>(b_ob.data().ptr.data), texspace_loc, texspace_size);
   }
 
-  bool load_metadata(const ImageDeviceFeatures & /*features*/, ImageMetaData &metadata) override
+  void load_grid() override
   {
     if (!b_domain) {
-      return false;
+      return;
     }
 
+    int channels;
     if (attribute == ATTR_STD_VOLUME_DENSITY || attribute == ATTR_STD_VOLUME_FLAME ||
         attribute == ATTR_STD_VOLUME_HEAT || attribute == ATTR_STD_VOLUME_TEMPERATURE)
     {
-      metadata.type = IMAGE_DATA_TYPE_FLOAT;
-      metadata.channels = 1;
+      channels = 1;
     }
     else if (attribute == ATTR_STD_VOLUME_COLOR) {
-      metadata.type = IMAGE_DATA_TYPE_FLOAT4;
-      metadata.channels = 4;
+      channels = 4;
     }
     else if (attribute == ATTR_STD_VOLUME_VELOCITY) {
-      metadata.type = IMAGE_DATA_TYPE_FLOAT4;
-      metadata.channels = 3;
+      channels = 3;
     }
     else {
-      return false;
+      return;
     }
 
     const int3 resolution = get_int3(b_domain.domain_resolution());
@@ -58,48 +59,41 @@ class BlenderSmokeLoader : public ImageLoader {
       amplify = 1;
     }
 
-    metadata.width = resolution.x * amplify;
-    metadata.height = resolution.y * amplify;
-    metadata.depth = resolution.z * amplify;
+    const size_t width = resolution.x * amplify;
+    const size_t height = resolution.y * amplify;
+    const size_t depth = resolution.z * amplify;
 
     /* Create a matrix to transform from object space to mesh texture space.
      * This does not work with deformations but that can probably only be done
      * well with a volume grid mapping of coordinates. */
-    metadata.transform_3d = transform_translate(-texspace_loc) * transform_scale(texspace_size);
-    metadata.use_transform_3d = true;
+    Transform transform_3d = transform_translate(-texspace_loc) * transform_scale(texspace_size);
 
-    return true;
+    vector<float> voxels;
+    if (!get_voxels(width, height, depth, channels, voxels)) {
+      return;
+    }
+
+    grid_from_dense_voxels(width, height, depth, channels, voxels.data(), transform_3d);
   }
 
-  bool load_pixels(const ImageMetaData & /*metadata*/,
-                   void *pixels,
-                   const size_t /*pixels_size*/,
-                   const bool /*associate_alpha*/) override
+  bool get_voxels(const size_t width,
+                  const size_t height,
+                  const size_t depth,
+                  const int channels,
+                  vector<float> &voxels)
   {
     if (!b_domain) {
       return false;
     }
+
 #ifdef WITH_FLUID
-    const int3 resolution = get_int3(b_domain.domain_resolution());
     int length;
-    int amplify = (b_domain.use_noise()) ? b_domain.noise_scale() : 1;
-
-    /* Velocity and heat data is always low-resolution. */
-    if (attribute == ATTR_STD_VOLUME_VELOCITY || attribute == ATTR_STD_VOLUME_HEAT) {
-      amplify = 1;
-    }
-
-    const int width = resolution.x * amplify;
-    const int height = resolution.y * amplify;
-    const int depth = resolution.z * amplify;
-    const size_t num_pixels = ((size_t)width) * height * depth;
-
-    float *fpixels = (float *)pixels;
+    voxels.resize(width * height * depth * channels);
 
     if (attribute == ATTR_STD_VOLUME_DENSITY) {
       FluidDomainSettings_density_grid_get_length(&b_domain.ptr, &length);
-      if (length == num_pixels) {
-        FluidDomainSettings_density_grid_get(&b_domain.ptr, fpixels);
+      if (length == voxels.size()) {
+        FluidDomainSettings_density_grid_get(&b_domain.ptr, voxels.data());
         return true;
       }
     }
@@ -107,50 +101,53 @@ class BlenderSmokeLoader : public ImageLoader {
       /* this is in range 0..1, and interpreted by the OpenGL smoke viewer
        * as 1500..3000 K with the first part faded to zero density */
       FluidDomainSettings_flame_grid_get_length(&b_domain.ptr, &length);
-      if (length == num_pixels) {
-        FluidDomainSettings_flame_grid_get(&b_domain.ptr, fpixels);
+      if (length == voxels.size()) {
+        FluidDomainSettings_flame_grid_get(&b_domain.ptr, voxels.data());
         return true;
       }
     }
     else if (attribute == ATTR_STD_VOLUME_COLOR) {
       /* the RGB is "premultiplied" by density for better interpolation results */
       FluidDomainSettings_color_grid_get_length(&b_domain.ptr, &length);
-      if (length == num_pixels * 4) {
-        FluidDomainSettings_color_grid_get(&b_domain.ptr, fpixels);
+      if (length == voxels.size()) {
+        FluidDomainSettings_color_grid_get(&b_domain.ptr, voxels.data());
         return true;
       }
     }
     else if (attribute == ATTR_STD_VOLUME_VELOCITY) {
       FluidDomainSettings_velocity_grid_get_length(&b_domain.ptr, &length);
-      if (length == num_pixels * 3) {
-        FluidDomainSettings_velocity_grid_get(&b_domain.ptr, fpixels);
+      if (length == voxels.size()) {
+        FluidDomainSettings_velocity_grid_get(&b_domain.ptr, voxels.data());
         return true;
       }
     }
     else if (attribute == ATTR_STD_VOLUME_HEAT) {
       FluidDomainSettings_heat_grid_get_length(&b_domain.ptr, &length);
-      if (length == num_pixels) {
-        FluidDomainSettings_heat_grid_get(&b_domain.ptr, fpixels);
+      if (length == voxels.size()) {
+        FluidDomainSettings_heat_grid_get(&b_domain.ptr, voxels.data());
         return true;
       }
     }
     else if (attribute == ATTR_STD_VOLUME_TEMPERATURE) {
       FluidDomainSettings_temperature_grid_get_length(&b_domain.ptr, &length);
-      if (length == num_pixels) {
-        FluidDomainSettings_temperature_grid_get(&b_domain.ptr, fpixels);
+      if (length == voxels.size()) {
+        FluidDomainSettings_temperature_grid_get(&b_domain.ptr, voxels.data());
         return true;
       }
     }
     else {
       LOG(ERROR) << "Unknown volume attribute " << Attribute::standard_name(attribute)
                  << "skipping ";
-      fpixels[0] = 0.0f;
+      voxels[0] = 0.0f;
       return false;
     }
-#else
-    (void)pixels;
-#endif
     LOG(ERROR) << "Unexpected smoke volume resolution, skipping";
+#else
+    (void)voxels;
+    (void)width;
+    (void)height;
+    (void)depth;
+#endif
     return false;
   }
 
@@ -209,7 +206,7 @@ static void sync_smoke_volume(
       continue;
     }
 
-    volume->set_clipping(b_domain.clipping());
+    const float clipping = b_domain.clipping();
 
     Attribute *attr = volume->attributes.add(std);
 
@@ -218,7 +215,8 @@ static void sync_smoke_volume(
       continue;
     }
 
-    unique_ptr<ImageLoader> loader = make_unique<BlenderSmokeLoader>(b_ob_info.real_object, std);
+    unique_ptr<ImageLoader> loader = make_unique<BlenderSmokeLoader>(
+        b_ob_info.real_object, std, clipping);
     ImageParams params;
     params.frame = frame;
 
@@ -231,8 +229,9 @@ class BlenderVolumeLoader : public VDBImageLoader {
   BlenderVolumeLoader(BL::BlendData &b_data,
                       BL::Volume &b_volume,
                       const string &grid_name,
-                      BL::VolumeRender::precision_enum precision_)
-      : VDBImageLoader(grid_name), b_volume(b_volume)
+                      BL::VolumeRender::precision_enum precision_,
+                      const float clipping)
+      : VDBImageLoader(grid_name, clipping), b_volume(b_volume)
   {
     b_volume.grids.load(b_data.ptr.data);
 
@@ -285,7 +284,6 @@ static void sync_volume_object(BL::BlendData &b_data,
 
   BL::VolumeRender b_render(b_volume.render());
 
-  volume->set_clipping(b_render.clipping());
   volume->set_step_size(b_render.step_size());
   volume->set_object_space((b_render.space() == BL::VolumeRender::space_OBJECT));
 
@@ -352,7 +350,7 @@ static void sync_volume_object(BL::BlendData &b_data,
                             volume->attributes.add(name, TypeFloat, ATTR_ELEMENT_VOXEL);
 
       unique_ptr<ImageLoader> loader = make_unique<BlenderVolumeLoader>(
-          b_data, b_volume, name.string(), b_render.precision());
+          b_data, b_volume, name.string(), b_render.precision(), b_render.clipping());
       ImageParams params;
       params.frame = b_volume.grids.frame();
 
@@ -377,6 +375,8 @@ void BlenderSync::sync_volume(BObjectInfo &b_ob_info, Volume *volume)
     }
   }
 
+  volume->merge_grids(scene);
+
   /* Tag update. */
   volume->tag_update(scene, true);
 }
diff --git a/intern/cycles/device/cuda/device_impl.cpp b/intern/cycles/device/cuda/device_impl.cpp
index 1c498143120..7914fd936a5 100644
--- a/intern/cycles/device/cuda/device_impl.cpp
+++ b/intern/cycles/device/cuda/device_impl.cpp
@@ -897,11 +897,7 @@ void CUDADevice::tex_alloc(device_texture &mem)
   /* Set Mapping and tag that we need to (re-)upload to device */
   TextureInfo tex_info = mem.info;
 
-  if (mem.info.data_type != IMAGE_DATA_TYPE_NANOVDB_FLOAT &&
-      mem.info.data_type != IMAGE_DATA_TYPE_NANOVDB_FLOAT3 &&
-      mem.info.data_type != IMAGE_DATA_TYPE_NANOVDB_FPN &&
-      mem.info.data_type != IMAGE_DATA_TYPE_NANOVDB_FP16)
-  {
+  if (!is_nanovdb_type(mem.info.data_type)) {
     CUDA_RESOURCE_DESC resDesc;
     memset(&resDesc, 0, sizeof(resDesc));
 
diff --git a/intern/cycles/device/hip/device_impl.cpp b/intern/cycles/device/hip/device_impl.cpp
index f51ab81d305..0bb4bce6c26 100644
--- a/intern/cycles/device/hip/device_impl.cpp
+++ b/intern/cycles/device/hip/device_impl.cpp
@@ -865,11 +865,7 @@ void HIPDevice::tex_alloc(device_texture &mem)
   /* Set Mapping and tag that we need to (re-)upload to device */
   TextureInfo tex_info = mem.info;
 
-  if (mem.info.data_type != IMAGE_DATA_TYPE_NANOVDB_FLOAT &&
-      mem.info.data_type != IMAGE_DATA_TYPE_NANOVDB_FLOAT3 &&
-      mem.info.data_type != IMAGE_DATA_TYPE_NANOVDB_FPN &&
-      mem.info.data_type != IMAGE_DATA_TYPE_NANOVDB_FP16)
-  {
+  if (!is_nanovdb_type(mem.info.data_type)) {
     /* Bindless textures. */
     hipResourceDesc resDesc;
     memset(&resDesc, 0, sizeof(resDesc));
diff --git a/intern/cycles/device/memory.cpp b/intern/cycles/device/memory.cpp
index ec2a46a6c34..86f03669429 100644
--- a/intern/cycles/device/memory.cpp
+++ b/intern/cycles/device/memory.cpp
@@ -180,6 +180,7 @@ device_texture::device_texture(Device *device,
     case IMAGE_DATA_TYPE_BYTE:
     case IMAGE_DATA_TYPE_NANOVDB_FLOAT:
     case IMAGE_DATA_TYPE_NANOVDB_FLOAT3:
+    case IMAGE_DATA_TYPE_NANOVDB_FLOAT4:
     case IMAGE_DATA_TYPE_NANOVDB_FPN:
     case IMAGE_DATA_TYPE_NANOVDB_FP16:
       data_type = TYPE_UCHAR;
diff --git a/intern/cycles/device/metal/device_impl.mm b/intern/cycles/device/metal/device_impl.mm
index 8aa75234cf8..a9d0decc8ab 100644
--- a/intern/cycles/device/metal/device_impl.mm
+++ b/intern/cycles/device/metal/device_impl.mm
@@ -937,11 +937,7 @@ void MetalDevice::tex_alloc_as_buffer(device_texture &mem)
   texture_info[slot] = mem.info;
   texture_slot_map[slot] = mmem->mtlBuffer;
 
-  if (mem.info.data_type == IMAGE_DATA_TYPE_NANOVDB_FLOAT ||
-      mem.info.data_type == IMAGE_DATA_TYPE_NANOVDB_FLOAT3 ||
-      mem.info.data_type == IMAGE_DATA_TYPE_NANOVDB_FPN ||
-      mem.info.data_type == IMAGE_DATA_TYPE_NANOVDB_FP16)
-  {
+  if (is_nanovdb_type(mem.info.data_type)) {
     using_nanovdb = true;
   }
 }
diff --git a/intern/cycles/device/oneapi/device_impl.cpp b/intern/cycles/device/oneapi/device_impl.cpp
index 42fcb4c24ff..7b2bcce8f90 100644
--- a/intern/cycles/device/oneapi/device_impl.cpp
+++ b/intern/cycles/device/oneapi/device_impl.cpp
@@ -842,11 +842,7 @@ void OneapiDevice::tex_alloc(device_texture &mem)
     sycl::ext::oneapi::experimental::bindless_image_sampler samp(
         address_mode, sycl::coordinate_normalization_mode::normalized, filter_mode);
 
-    if (mem.info.data_type != IMAGE_DATA_TYPE_NANOVDB_FLOAT &&
-        mem.info.data_type != IMAGE_DATA_TYPE_NANOVDB_FLOAT3 &&
-        mem.info.data_type != IMAGE_DATA_TYPE_NANOVDB_FPN &&
-        mem.info.data_type != IMAGE_DATA_TYPE_NANOVDB_FP16)
-    {
+    if (!is_nanovdb_type(mem.info.data_type)) {
       sycl::ext::oneapi::experimental::sampled_image_handle imgHandle;
 
       if (memHandle.raw_handle) {
diff --git a/intern/cycles/hydra/volume.cpp b/intern/cycles/hydra/volume.cpp
index c9c37578807..27fecdc187d 100644
--- a/intern/cycles/hydra/volume.cpp
+++ b/intern/cycles/hydra/volume.cpp
@@ -84,6 +84,8 @@ void HdCyclesVolume::Populate(HdSceneDelegate *sceneDelegate, HdDirtyBits dirtyB
       }
     }
 
+    _geom->merge_grids(scene);
+
     rebuild = true;
   }
 }
diff --git a/intern/cycles/kernel/device/cpu/image.h b/intern/cycles/kernel/device/cpu/image.h
index 4e4dd142b30..4d3296a1f83 100644
--- a/intern/cycles/kernel/device/cpu/image.h
+++ b/intern/cycles/kernel/device/cpu/image.h
@@ -713,6 +713,11 @@ template<typename TexT, typename OutT> struct NanoVDBInterpolator {
     return make_float4(r.x, r.y, r.z, 1.0f);
   }
 
+  static ccl_always_inline float4 read(const float4 r)
+  {
+    return r;
+  }
+
   template<typename Acc>
   static ccl_always_inline OutT
   interp_3d_closest(const Acc &acc, const float x, float y, const float z)
@@ -917,6 +922,8 @@ ccl_device float4 kernel_tex_image_interp_3d(KernelGlobals kg,
     }
     case IMAGE_DATA_TYPE_NANOVDB_FLOAT3:
       return NanoVDBInterpolator<packed_float3, float4>::interp_3d(info, P.x, P.y, P.z, interp);
+    case IMAGE_DATA_TYPE_NANOVDB_FLOAT4:
+      return NanoVDBInterpolator<float4, float4>::interp_3d(info, P.x, P.y, P.z, interp);
     case IMAGE_DATA_TYPE_NANOVDB_FPN: {
       const float f = NanoVDBInterpolator<nanovdb::FpN, float>::interp_3d(
           info, P.x, P.y, P.z, interp);
@@ -929,7 +936,6 @@ ccl_device float4 kernel_tex_image_interp_3d(KernelGlobals kg,
     }
 #endif
     default:
-      assert(0);
       return make_float4(
           TEX_IMAGE_MISSING_R, TEX_IMAGE_MISSING_G, TEX_IMAGE_MISSING_B, TEX_IMAGE_MISSING_A);
   }
diff --git a/intern/cycles/kernel/device/gpu/image.h b/intern/cycles/kernel/device/gpu/image.h
index d15b3b4c6ed..dfac8c17224 100644
--- a/intern/cycles/kernel/device/gpu/image.h
+++ b/intern/cycles/kernel/device/gpu/image.h
@@ -321,6 +321,9 @@ ccl_device float4 kernel_tex_image_interp_3d(KernelGlobals kg,
         info, x, y, z, interpolation);
     return make_float4(f, 1.0f);
   }
+  if (texture_type == IMAGE_DATA_TYPE_NANOVDB_FLOAT4) {
+    return kernel_tex_image_interp_nanovdb<float4, float4>(info, x, y, z, interpolation);
+  }
   if (texture_type == IMAGE_DATA_TYPE_NANOVDB_FPN) {
     float f = kernel_tex_image_interp_nanovdb<float, nanovdb::FpN>(info, x, y, z, interpolation);
     return make_float4(f, f, f, 1.0f);
diff --git a/intern/cycles/scene/geometry.cpp b/intern/cycles/scene/geometry.cpp
index bcd0b28ad55..1f7596e3859 100644
--- a/intern/cycles/scene/geometry.cpp
+++ b/intern/cycles/scene/geometry.cpp
@@ -664,13 +664,9 @@ void GeometryManager::device_update_volume_images(Device *device, Scene *scene,
       }
 
       const ImageHandle &handle = attr.data_voxel();
-      /* We can build directly from OpenVDB data structures, no need to
-       * load such images early. */
-      if (!handle.vdb_loader()) {
-        const int slot = handle.svm_slot();
-        if (slot != -1) {
-          volume_images.insert(slot);
-        }
+      const int slot = handle.svm_slot();
+      if (slot != -1) {
+        volume_images.insert(slot);
       }
     }
   }
diff --git a/intern/cycles/scene/image.cpp b/intern/cycles/scene/image.cpp
index 15173eee369..fbddafb516e 100644
--- a/intern/cycles/scene/image.cpp
+++ b/intern/cycles/scene/image.cpp
@@ -48,6 +48,8 @@ const char *name_from_type(ImageDataType type)
       return "nanovdb_float";
     case IMAGE_DATA_TYPE_NANOVDB_FLOAT3:
       return "nanovdb_float3";
+    case IMAGE_DATA_TYPE_NANOVDB_FLOAT4:
+      return "nanovdb_float4";
     case IMAGE_DATA_TYPE_NANOVDB_FPN:
       return "nanovdb_fpn";
     case IMAGE_DATA_TYPE_NANOVDB_FP16:
@@ -367,10 +369,7 @@ void ImageManager::load_image_metadata(Image *img)
 
   metadata.detect_colorspace();
 
-  assert(features.has_nanovdb || (metadata.type != IMAGE_DATA_TYPE_NANOVDB_FLOAT ||
-                                  metadata.type != IMAGE_DATA_TYPE_NANOVDB_FLOAT3 ||
-                                  metadata.type != IMAGE_DATA_TYPE_NANOVDB_FPN ||
-                                  metadata.type != IMAGE_DATA_TYPE_NANOVDB_FP16));
+  assert(features.has_nanovdb || !is_nanovdb_type(metadata.type));
 
   img->need_metadata = false;
 }
@@ -811,9 +810,7 @@ void ImageManager::device_load_image(Device *device,
     }
   }
 #ifdef WITH_NANOVDB
-  else if (type == IMAGE_DATA_TYPE_NANOVDB_FLOAT || type == IMAGE_DATA_TYPE_NANOVDB_FLOAT3 ||
-           type == IMAGE_DATA_TYPE_NANOVDB_FPN || type == IMAGE_DATA_TYPE_NANOVDB_FP16)
-  {
+  else if (is_nanovdb_type(type)) {
     const thread_scoped_lock device_lock(device_mutex);
     void *pixels = img->mem->alloc(img->metadata.byte_size, 0);
 
diff --git a/intern/cycles/scene/image.h b/intern/cycles/scene/image.h
index c1a4642123d..dc47cb47bd2 100644
--- a/intern/cycles/scene/image.h
+++ b/intern/cycles/scene/image.h
@@ -80,7 +80,7 @@ class ImageMetaData {
 /* Information about supported features that Image loaders can use. */
 class ImageDeviceFeatures {
  public:
-  bool has_nanovdb;
+  bool has_nanovdb = true;
 };
 
 /* Image loader base class, that can be subclassed to load image data
diff --git a/intern/cycles/scene/image_oiio.cpp b/intern/cycles/scene/image_oiio.cpp
index 2840f2ca130..d2678861829 100644
--- a/intern/cycles/scene/image_oiio.cpp
+++ b/intern/cycles/scene/image_oiio.cpp
@@ -237,6 +237,7 @@ bool OIIOImageLoader::load_pixels(const ImageMetaData &metadata,
       break;
     case IMAGE_DATA_TYPE_NANOVDB_FLOAT:
     case IMAGE_DATA_TYPE_NANOVDB_FLOAT3:
+    case IMAGE_DATA_TYPE_NANOVDB_FLOAT4:
     case IMAGE_DATA_TYPE_NANOVDB_FPN:
     case IMAGE_DATA_TYPE_NANOVDB_FP16:
     case IMAGE_DATA_NUM_TYPES:
diff --git a/intern/cycles/scene/image_vdb.cpp b/intern/cycles/scene/image_vdb.cpp
index bb1a4eb536c..0eb136a6b5e 100644
--- a/intern/cycles/scene/image_vdb.cpp
+++ b/intern/cycles/scene/image_vdb.cpp
@@ -15,33 +15,30 @@
 CCL_NAMESPACE_BEGIN
 
 #ifdef WITH_OPENVDB
-struct ToDenseOp {
-  openvdb::CoordBBox bbox;
-  void *pixels;
-
-  template<typename GridType, typename FloatDataType, const int channels>
-  bool operator()(const typename GridType::ConstPtr &grid)
-  {
-    openvdb::tools::Dense<FloatDataType, openvdb::tools::LayoutXYZ> dense(bbox,
-                                                                          (FloatDataType *)pixels);
-    openvdb::tools::copyToDense(*grid, dense);
-    return true;
-  }
-};
-
-VDBImageLoader::VDBImageLoader(openvdb::GridBase::ConstPtr grid_, const string &grid_name)
-    : grid_name(grid_name), grid(grid_)
+VDBImageLoader::VDBImageLoader(openvdb::GridBase::ConstPtr grid_,
+                               const string &grid_name,
+                               const float clipping)
+    : grid_name(grid_name), clipping(clipping), grid(grid_)
 {
 }
 #endif
 
-VDBImageLoader::VDBImageLoader(const string &grid_name) : grid_name(grid_name) {}
+VDBImageLoader::VDBImageLoader(const string &grid_name, const float clipping)
+    : grid_name(grid_name), clipping(clipping)
+{
+}
 
 VDBImageLoader::~VDBImageLoader() = default;
 
 bool VDBImageLoader::load_metadata(const ImageDeviceFeatures &features, ImageMetaData &metadata)
 {
-#ifdef WITH_OPENVDB
+  if (!features.has_nanovdb) {
+    return false;
+  }
+
+#ifdef WITH_NANOVDB
+  load_grid();
+
   if (!grid) {
     return false;
   }
@@ -55,56 +52,43 @@ bool VDBImageLoader::load_metadata(const ImageDeviceFeatures &features, ImageMet
   /* Get number of channels from type. */
   metadata.channels = openvdb_num_channels(grid);
 
-  /* Set data type. */
-#  ifdef WITH_NANOVDB
-  if (features.has_nanovdb) {
-    /* Convert OpenVDB to NanoVDB grid. */
-    nanogrid = openvdb_to_nanovdb(grid, precision, 0.0f);
-    if (!nanogrid) {
-      return false;
-    }
+  /* Convert OpenVDB to NanoVDB grid. */
+  nanogrid = openvdb_to_nanovdb(grid, precision, clipping);
+  if (!nanogrid) {
+    grid.reset();
+    return false;
   }
-#  endif
 
   /* Set dimensions. */
   bbox = grid->evalActiveVoxelBoundingBox();
   if (bbox.empty()) {
+    grid.reset();
     return false;
   }
 
-  openvdb::Coord dim = bbox.dim();
-  metadata.width = dim.x();
-  metadata.height = dim.y();
-  metadata.depth = dim.z();
+  if (metadata.channels == 1) {
+    if (precision == 0) {
+      metadata.type = IMAGE_DATA_TYPE_NANOVDB_FPN;
+    }
+    else if (precision == 16) {
+      metadata.type = IMAGE_DATA_TYPE_NANOVDB_FP16;
+    }
+    else {
+      metadata.type = IMAGE_DATA_TYPE_NANOVDB_FLOAT;
+    }
+  }
+  else if (metadata.channels == 3) {
+    metadata.type = IMAGE_DATA_TYPE_NANOVDB_FLOAT3;
+  }
+  else if (metadata.channels == 4) {
+    metadata.type = IMAGE_DATA_TYPE_NANOVDB_FLOAT4;
+  }
+  else {
+    grid.reset();
+    return false;
+  }
 
-#  ifdef WITH_NANOVDB
-  if (nanogrid) {
-    metadata.byte_size = nanogrid.size();
-    if (metadata.channels == 1) {
-      if (precision == 0) {
-        metadata.type = IMAGE_DATA_TYPE_NANOVDB_FPN;
-      }
-      else if (precision == 16) {
-        metadata.type = IMAGE_DATA_TYPE_NANOVDB_FP16;
-      }
-      else {
-        metadata.type = IMAGE_DATA_TYPE_NANOVDB_FLOAT;
-      }
-    }
-    else {
-      metadata.type = IMAGE_DATA_TYPE_NANOVDB_FLOAT3;
-    }
-  }
-  else
-#  endif
-  {
-    if (metadata.channels == 1) {
-      metadata.type = IMAGE_DATA_TYPE_FLOAT;
-    }
-    else {
-      metadata.type = IMAGE_DATA_TYPE_FLOAT4;
-    }
-  }
+  metadata.byte_size = nanogrid.size();
 
   /* Set transform from object space to voxel index. */
   openvdb::math::Mat4f grid_matrix = grid->transform().baseMap()->getAffineMap()->getMat4();
@@ -115,29 +99,15 @@ bool VDBImageLoader::load_metadata(const ImageDeviceFeatures &features, ImageMet
     }
   }
 
-  Transform texture_to_index;
-#  ifdef WITH_NANOVDB
-  if (nanogrid) {
-    texture_to_index = transform_identity();
-  }
-  else
-#  endif
-  {
-    openvdb::Coord min = bbox.min();
-    texture_to_index = transform_translate(min.x(), min.y(), min.z()) *
-                       transform_scale(dim.x(), dim.y(), dim.z());
-  }
-
-  metadata.transform_3d = transform_inverse(index_to_object * texture_to_index);
+  metadata.transform_3d = transform_inverse(index_to_object);
   metadata.use_transform_3d = true;
 
-#  ifndef WITH_NANOVDB
-  (void)features;
-#  endif
+  /* Only NanoGrid needed now, free OpenVDB grid. */
+  grid.reset();
+
   return true;
 #else
   (void)metadata;
-  (void)features;
   return false;
 #endif
 }
@@ -147,24 +117,16 @@ bool VDBImageLoader::load_pixels(const ImageMetaData & /*metadata*/,
                                  const size_t /*pixels_size*/,
                                  const bool /*associate_alpha*/)
 {
-#ifdef WITH_OPENVDB
-#  ifdef WITH_NANOVDB
+#ifdef WITH_NANOVDB
   if (nanogrid) {
     memcpy(pixels, nanogrid.data(), nanogrid.size());
+    return true;
   }
-  else
-#  endif
-  {
-    ToDenseOp op;
-    op.pixels = pixels;
-    op.bbox = bbox;
-    openvdb_grid_type_operation(grid, op);
-  }
-  return true;
 #else
   (void)pixels;
-  return false;
 #endif
+
+  return false;
 }
 
 string VDBImageLoader::name() const
@@ -176,7 +138,7 @@ bool VDBImageLoader::equals(const ImageLoader &other) const
 {
 #ifdef WITH_OPENVDB
   const VDBImageLoader &other_loader = (const VDBImageLoader &)other;
-  return grid == other_loader.grid;
+  return grid && grid == other_loader.grid;
 #else
   (void)other;
   return true;
@@ -204,6 +166,89 @@ openvdb::GridBase::ConstPtr VDBImageLoader::get_grid()
 {
   return grid;
 }
+
+template<typename GridType>
+openvdb::GridBase::ConstPtr create_grid(const float *voxels,
+                                        const size_t width,
+                                        const size_t height,
+                                        const size_t depth,
+                                        Transform transform_3d,
+                                        const float clipping)
+{
+  using ValueType = typename GridType::ValueType;
+  openvdb::GridBase::ConstPtr grid;
+
+  const openvdb::CoordBBox dense_bbox(0, 0, 0, width - 1, height - 1, depth - 1);
+
+  typename GridType::Ptr sparse = GridType::create(ValueType(0.0f));
+  if (dense_bbox.empty()) {
+    return sparse;
+  }
+
+  const openvdb::tools::Dense<const ValueType, openvdb::tools::MemoryLayout::LayoutXYZ> dense(
+      dense_bbox, reinterpret_cast<const ValueType *>(voxels));
+
+  openvdb::tools::copyFromDense(dense, *sparse, ValueType(clipping));
+
+  /* Compute index to world matrix. */
+  const float3 voxel_size = make_float3(1.0f / width, 1.0f / height, 1.0f / depth);
+
+  transform_3d = transform_inverse(transform_3d);
+
+  const openvdb::Mat4R index_to_world_mat((double)(voxel_size.x * transform_3d[0][0]),
+                                          0.0,
+                                          0.0,
+                                          0.0,
+                                          0.0,
+                                          (double)(voxel_size.y * transform_3d[1][1]),
+                                          0.0,
+                                          0.0,
+                                          0.0,
+                                          0.0,
+                                          (double)(voxel_size.z * transform_3d[2][2]),
+                                          0.0,
+                                          (double)transform_3d[0][3],
+                                          (double)transform_3d[1][3],
+                                          (double)transform_3d[2][3],
+                                          1.0);
+
+  const openvdb::math::Transform::Ptr index_to_world_tfm =
+      openvdb::math::Transform::createLinearTransform(index_to_world_mat);
+
+  sparse->setTransform(index_to_world_tfm);
+
+  return sparse;
+}
 #endif
 
+void VDBImageLoader::grid_from_dense_voxels(const size_t width,
+                                            const size_t height,
+                                            const size_t depth,
+                                            const int channels,
+                                            const float *voxels,
+                                            Transform transform_3d)
+{
+#ifdef WITH_OPENVDB
+  /* TODO: Create NanoVDB grid directly? */
+  if (channels == 1) {
+    grid = create_grid<openvdb::FloatGrid>(voxels, width, height, depth, transform_3d, clipping);
+  }
+  else if (channels == 3) {
+    grid = create_grid<openvdb::Vec3fGrid>(voxels, width, height, depth, transform_3d, clipping);
+  }
+  else if (channels == 4) {
+    grid = create_grid<openvdb::Vec4fGrid>(voxels, width, height, depth, transform_3d, clipping);
+  }
+
+  /* Clipping already applied, no need to do it again. */
+  clipping = 0.0f;
+#else
+  (void)width;
+  (void)height;
+  (void)depth;
+  (void)channels;
+  (void)voxels;
+#endif
+}
+
 CCL_NAMESPACE_END
diff --git a/intern/cycles/scene/image_vdb.h b/intern/cycles/scene/image_vdb.h
index e516febf9ab..76a0565256f 100644
--- a/intern/cycles/scene/image_vdb.h
+++ b/intern/cycles/scene/image_vdb.h
@@ -19,22 +19,26 @@
 
 #include "scene/image.h"
 
+#include "util/transform.h"
+
 CCL_NAMESPACE_BEGIN
 
 class VDBImageLoader : public ImageLoader {
  public:
 #ifdef WITH_OPENVDB
-  VDBImageLoader(openvdb::GridBase::ConstPtr grid_, const string &grid_name);
+  VDBImageLoader(openvdb::GridBase::ConstPtr grid_,
+                 const string &grid_name,
+                 const float clipping = 0.001f);
 #endif
-  VDBImageLoader(const string &grid_name);
+  VDBImageLoader(const string &grid_name, const float clipping = 0.001f);
   ~VDBImageLoader() override;
 
-  bool load_metadata(const ImageDeviceFeatures &features, ImageMetaData &metadata) override;
+  bool load_metadata(const ImageDeviceFeatures &features, ImageMetaData &metadata) final;
 
   bool load_pixels(const ImageMetaData &metadata,
                    void *pixels,
                    const size_t pixels_size,
-                   const bool associate_alpha) override;
+                   const bool associate_alpha) final;
 
   string name() const override;
 
@@ -49,14 +53,24 @@ class VDBImageLoader : public ImageLoader {
 #endif
 
  protected:
+  virtual void load_grid() {}
+
+  void grid_from_dense_voxels(const size_t width,
+                              const size_t height,
+                              const size_t depth,
+                              const int channels,
+                              const float *voxels,
+                              Transform transform_3d);
+
   string grid_name;
+  float clipping = 0.001f;
 #ifdef WITH_OPENVDB
   openvdb::GridBase::ConstPtr grid;
   openvdb::CoordBBox bbox;
 #endif
 #ifdef WITH_NANOVDB
   nanovdb::GridHandle<> nanogrid;
-  int precision = 0;
+  int precision = 16;
 #endif
 };
 
diff --git a/intern/cycles/scene/object.cpp b/intern/cycles/scene/object.cpp
index 9077625249a..eac0013e114 100644
--- a/intern/cycles/scene/object.cpp
+++ b/intern/cycles/scene/object.cpp
@@ -327,7 +327,7 @@ float Object::compute_volume_step_size() const
       if (attr.element == ATTR_ELEMENT_VOXEL) {
         ImageHandle &handle = attr.data_voxel();
         const ImageMetaData &metadata = handle.metadata();
-        if (metadata.width == 0 || metadata.height == 0 || metadata.depth == 0) {
+        if (metadata.byte_size == 0) {
           continue;
         }
 
@@ -335,25 +335,13 @@ float Object::compute_volume_step_size() const
         float voxel_step_size = volume->get_step_size();
 
         if (voxel_step_size == 0.0f) {
-          /* Auto detect step size. */
-          float3 size = one_float3();
-#ifdef WITH_NANOVDB
-          /* Dimensions were not applied to image transform with NanoVDB (see image_vdb.cpp) */
-          if (metadata.type != IMAGE_DATA_TYPE_NANOVDB_FLOAT &&
-              metadata.type != IMAGE_DATA_TYPE_NANOVDB_FLOAT3 &&
-              metadata.type != IMAGE_DATA_TYPE_NANOVDB_FPN &&
-              metadata.type != IMAGE_DATA_TYPE_NANOVDB_FP16)
-#endif
-          {
-            size /= make_float3(metadata.width, metadata.height, metadata.depth);
-          }
-
-          /* Step size is transformed from voxel to world space. */
+          /* Auto detect step size.
+           * Step size is transformed from voxel to world space. */
           Transform voxel_tfm = tfm;
           if (metadata.use_transform_3d) {
             voxel_tfm = tfm * transform_inverse(metadata.transform_3d);
           }
-          voxel_step_size = reduce_min(fabs(transform_direction(&voxel_tfm, size)));
+          voxel_step_size = reduce_min(fabs(transform_direction(&voxel_tfm, one_float3())));
         }
         else if (volume->get_object_space()) {
           /* User specified step size in object space. */
diff --git a/intern/cycles/scene/volume.cpp b/intern/cycles/scene/volume.cpp
index 7f7076231a6..db605b583b9 100644
--- a/intern/cycles/scene/volume.cpp
+++ b/intern/cycles/scene/volume.cpp
@@ -8,15 +8,13 @@
 #include "scene/scene.h"
 
 #ifdef WITH_OPENVDB
-#  include <openvdb/tools/Dense.h>
 #  include <openvdb/tools/GridTransformer.h>
 #  include <openvdb/tools/Morphology.h>
-#  include <openvdb/tools/Statistics.h>
 #endif
 
 #include "util/hash.h"
 #include "util/log.h"
-#include "util/openvdb.h"
+#include "util/nanovdb.h"
 #include "util/progress.h"
 #include "util/types.h"
 
@@ -26,7 +24,6 @@ NODE_DEFINE(Volume)
 {
   NodeType *type = NodeType::add("volume", create, NodeType::NONE, Mesh::get_node_type());
 
-  SOCKET_FLOAT(clipping, "Clipping", 0.001f);
   SOCKET_FLOAT(step_size, "Step Size", 0.0f);
   SOCKET_BOOLEAN(object_space, "Object Space", false);
   SOCKET_FLOAT(velocity_scale, "Velocity Scale", 1.0f);
@@ -36,7 +33,6 @@ NODE_DEFINE(Volume)
 
 Volume::Volume() : Mesh(get_node_type(), Geometry::VOLUME)
 {
-  clipping = 0.001f;
   step_size = 0.0f;
   object_space = false;
 }
@@ -61,7 +57,7 @@ enum {
   QUAD_Z_MAX = 5,
 };
 
-#ifdef WITH_OPENVDB
+#if defined(WITH_OPENVDB) && defined(WITH_NANOVDB)
 const int quads_indices[6][4] = {
     /* QUAD_X_MIN */
     {4, 0, 3, 7},
@@ -127,7 +123,6 @@ static void create_quad(const int3 corners[8],
 
   quads.push_back(quad);
 }
-#endif
 
 /* Create a mesh from a volume.
  *
@@ -141,18 +136,14 @@ static void create_quad(const int3 corners[8],
  */
 class VolumeMeshBuilder {
  public:
-#ifdef WITH_OPENVDB
   /* use a MaskGrid to store the topology to save memory */
   openvdb::MaskGrid::Ptr topology_grid;
   openvdb::CoordBBox bbox;
-#endif
   bool first_grid;
 
   VolumeMeshBuilder();
 
-#ifdef WITH_OPENVDB
-  void add_grid(openvdb::GridBase::ConstPtr grid, bool do_clipping, const float volume_clipping);
-#endif
+  void add_grid(const nanovdb::GridHandle<> &grid);
 
   void add_padding(const int pad_size);
 
@@ -169,30 +160,6 @@ class VolumeMeshBuilder {
   void convert_quads_to_tris(const vector<QuadData> &quads, vector<int> &tris);
 
   bool empty_grid() const;
-
-#ifdef WITH_OPENVDB
-  template<typename GridType>
-  void merge_grid(openvdb::GridBase::ConstPtr grid, bool do_clipping, const float volume_clipping)
-  {
-    typename GridType::ConstPtr typed_grid = openvdb::gridConstPtrCast<GridType>(grid);
-
-    if (do_clipping) {
-      using ValueType = typename GridType::ValueType;
-      const typename GridType::Ptr copy = typed_grid->deepCopy();
-      typename GridType::ValueOnIter iter = copy->beginValueOn();
-
-      for (; iter; ++iter) {
-        if (openvdb::math::Abs(iter.getValue()) < ValueType(volume_clipping)) {
-          iter.setValueOff();
-        }
-      }
-
-      typed_grid = copy;
-    }
-
-    topology_grid->topologyUnion(*typed_grid);
-  }
-#endif
 };
 
 VolumeMeshBuilder::VolumeMeshBuilder()
@@ -200,21 +167,21 @@ VolumeMeshBuilder::VolumeMeshBuilder()
   first_grid = true;
 }
 
-#ifdef WITH_OPENVDB
-void VolumeMeshBuilder::add_grid(openvdb::GridBase::ConstPtr grid,
-                                 bool do_clipping,
-                                 const float volume_clipping)
+void VolumeMeshBuilder::add_grid(const nanovdb::GridHandle<> &nanogrid)
 {
   /* set the transform of our grid from the first one */
+  openvdb::MaskGrid::Ptr grid = nanovdb_to_openvdb_mask(nanogrid);
+
   if (first_grid) {
-    topology_grid = openvdb::MaskGrid::create();
-    topology_grid->setTransform(grid->transform().copy());
+    topology_grid = grid;
     first_grid = false;
+    return;
   }
-  /* if the transforms do not match, we need to resample one of the grids so that
+
+  /* If the transforms do not match, we need to resample one of the grids so that
    * its index space registers with that of the other, here we resample our mask
    * grid so memory usage is kept low */
-  else if (topology_grid->transform() != grid->transform()) {
+  if (topology_grid->transform() != grid->transform()) {
     const openvdb::MaskGrid::Ptr temp_grid = topology_grid->copyWithNewTree();
     temp_grid->setTransform(grid->transform().copy());
     openvdb::tools::resampleToMatch<openvdb::tools::BoxSampler>(*topology_grid, *temp_grid);
@@ -222,54 +189,19 @@ void VolumeMeshBuilder::add_grid(openvdb::GridBase::ConstPtr grid,
     topology_grid->setTransform(grid->transform().copy());
   }
 
-  if (grid->isType<openvdb::FloatGrid>()) {
-    merge_grid<openvdb::FloatGrid>(grid, do_clipping, volume_clipping);
-  }
-  else if (grid->isType<openvdb::Vec3fGrid>()) {
-    merge_grid<openvdb::Vec3fGrid>(grid, do_clipping, volume_clipping);
-  }
-  else if (grid->isType<openvdb::Vec4fGrid>()) {
-    merge_grid<openvdb::Vec4fGrid>(grid, do_clipping, volume_clipping);
-  }
-  else if (grid->isType<openvdb::BoolGrid>()) {
-    merge_grid<openvdb::BoolGrid>(grid, do_clipping, volume_clipping);
-  }
-  else if (grid->isType<openvdb::DoubleGrid>()) {
-    merge_grid<openvdb::DoubleGrid>(grid, do_clipping, volume_clipping);
-  }
-  else if (grid->isType<openvdb::Int32Grid>()) {
-    merge_grid<openvdb::Int32Grid>(grid, do_clipping, volume_clipping);
-  }
-  else if (grid->isType<openvdb::Int64Grid>()) {
-    merge_grid<openvdb::Int64Grid>(grid, do_clipping, volume_clipping);
-  }
-  else if (grid->isType<openvdb::Vec3IGrid>()) {
-    merge_grid<openvdb::Vec3IGrid>(grid, do_clipping, volume_clipping);
-  }
-  else if (grid->isType<openvdb::Vec3dGrid>()) {
-    merge_grid<openvdb::Vec3dGrid>(grid, do_clipping, volume_clipping);
-  }
-  else if (grid->isType<openvdb::MaskGrid>()) {
-    topology_grid->topologyUnion(*openvdb::gridConstPtrCast<openvdb::MaskGrid>(grid));
-  }
+  topology_grid->topologyUnion(*grid);
 }
-#endif
 
 void VolumeMeshBuilder::add_padding(const int pad_size)
 {
-#ifdef WITH_OPENVDB
   openvdb::tools::dilateActiveValues(
       topology_grid->tree(), pad_size, openvdb::tools::NN_FACE, openvdb::tools::IGNORE_TILES);
-#else
-  (void)pad_size;
-#endif
 }
 
 void VolumeMeshBuilder::create_mesh(vector<float3> &vertices,
                                     vector<int> &indices,
                                     const float face_overlap_avoidance)
 {
-#ifdef WITH_OPENVDB
   /* We create vertices in index space (is), and only convert them to object
    * space when done. */
   vector<int3> vertices_is;
@@ -284,25 +216,17 @@ void VolumeMeshBuilder::create_mesh(vector<float3> &vertices,
   convert_object_space(vertices_is, vertices, face_overlap_avoidance);
 
   convert_quads_to_tris(quads, indices);
-#else
-  (void)vertices;
-  (void)indices;
-  (void)face_overlap_avoidance;
-#endif
 }
 
-#ifdef WITH_OPENVDB
 static bool is_non_empty_leaf(const openvdb::MaskGrid::TreeType &tree, const openvdb::Coord coord)
 {
   const auto *leaf_node = tree.probeLeaf(coord);
   return (leaf_node && !leaf_node->isEmpty());
 }
-#endif
 
 void VolumeMeshBuilder::generate_vertices_and_quads(vector<ccl::int3> &vertices_is,
                                                     vector<QuadData> &quads)
 {
-#ifdef WITH_OPENVDB
   const openvdb::MaskGrid::TreeType &tree = topology_grid->tree();
   tree.evalLeafBoundingBox(bbox);
 
@@ -366,17 +290,12 @@ void VolumeMeshBuilder::generate_vertices_and_quads(vector<ccl::int3> &vertices_
       create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_Z_MAX);
     }
   }
-#else
-  (void)vertices_is;
-  (void)quads;
-#endif
 }
 
 void VolumeMeshBuilder::convert_object_space(const vector<int3> &vertices,
                                              vector<float3> &out_vertices,
                                              const float face_overlap_avoidance)
 {
-#ifdef WITH_OPENVDB
   /* compute the offset for the face overlap avoidance */
   bbox = topology_grid->evalActiveVoxelBoundingBox();
   openvdb::Coord dim = bbox.dim();
@@ -392,11 +311,6 @@ void VolumeMeshBuilder::convert_object_space(const vector<int3> &vertices,
     const float3 vertex = make_float3((float)p.x(), (float)p.y(), (float)p.z());
     out_vertices.push_back(vertex + point_offset);
   }
-#else
-  (void)vertices;
-  (void)out_vertices;
-  (void)face_overlap_avoidance;
-#endif
 }
 
 void VolumeMeshBuilder::convert_quads_to_tris(const vector<QuadData> &quads, vector<int> &tris)
@@ -417,115 +331,46 @@ void VolumeMeshBuilder::convert_quads_to_tris(const vector<QuadData> &quads, vec
 
 bool VolumeMeshBuilder::empty_grid() const
 {
-#ifdef WITH_OPENVDB
   return !topology_grid ||
          (!topology_grid->tree().hasActiveTiles() && topology_grid->tree().leafCount() == 0);
-#else
-  return true;
-#endif
 }
 
-#ifdef WITH_OPENVDB
-template<typename GridType>
-static openvdb::GridBase::ConstPtr openvdb_grid_from_device_texture(device_texture *image_memory,
-                                                                    const float volume_clipping,
-                                                                    Transform transform_3d)
-{
-  using ValueType = typename GridType::ValueType;
-
-  const openvdb::CoordBBox dense_bbox(0,
-                                      0,
-                                      0,
-                                      image_memory->data_width - 1,
-                                      image_memory->data_height - 1,
-                                      image_memory->data_depth - 1);
-
-  typename GridType::Ptr sparse = GridType::create(ValueType(0.0f));
-  if (dense_bbox.empty()) {
-    return sparse;
-  }
-
-  const openvdb::tools::Dense<ValueType, openvdb::tools::MemoryLayout::LayoutXYZ> dense(
-      dense_bbox, static_cast<ValueType *>(image_memory->host_pointer));
-
-  openvdb::tools::copyFromDense(dense, *sparse, ValueType(volume_clipping));
-
-  /* #copyFromDense will remove any leaf node that contains constant data and replace it with a
-   * tile, however, we need to preserve the leaves in order to generate the mesh, so re-voxelize
-   * the leaves that were pruned. This should not affect areas that were skipped due to the
-   * volume_clipping parameter. */
-  sparse->tree().voxelizeActiveTiles();
-
-  /* Compute index to world matrix. */
-  const float3 voxel_size = make_float3(1.0f / image_memory->data_width,
-                                        1.0f / image_memory->data_height,
-                                        1.0f / image_memory->data_depth);
-
-  transform_3d = transform_inverse(transform_3d);
-
-  const openvdb::Mat4R index_to_world_mat((double)(voxel_size.x * transform_3d[0][0]),
-                                          0.0,
-                                          0.0,
-                                          0.0,
-                                          0.0,
-                                          (double)(voxel_size.y * transform_3d[1][1]),
-                                          0.0,
-                                          0.0,
-                                          0.0,
-                                          0.0,
-                                          (double)(voxel_size.z * transform_3d[2][2]),
-                                          0.0,
-                                          (double)transform_3d[0][3],
-                                          (double)transform_3d[1][3],
-                                          (double)transform_3d[2][3],
-                                          1.0);
-
-  const openvdb::math::Transform::Ptr index_to_world_tfm =
-      openvdb::math::Transform::createLinearTransform(index_to_world_mat);
-
-  sparse->setTransform(index_to_world_tfm);
-
-  return sparse;
-}
-
-static int estimate_required_velocity_padding(openvdb::GridBase::ConstPtr grid,
+static int estimate_required_velocity_padding(const nanovdb::GridHandle<> &grid,
                                               const float velocity_scale)
 {
-  /* TODO: we may need to also find outliers and clamp them to avoid adding too much padding. */
-  openvdb::math::Extrema extrema;
-  openvdb::Vec3d voxel_size;
+  const auto *typed_grid = grid.template grid<nanovdb::Vec3f>(0);
 
-  /* External `.vdb` files have a vec3 type for velocity,
-   * but the Blender exporter creates a vec4. */
-  if (grid->isType<openvdb::Vec3fGrid>()) {
-    const openvdb::Vec3fGrid::ConstPtr vel_grid = openvdb::gridConstPtrCast<openvdb::Vec3fGrid>(
-        grid);
-    extrema = openvdb::tools::extrema(vel_grid->cbeginValueOn());
-    voxel_size = vel_grid->voxelSize();
-  }
-  else if (grid->isType<openvdb::Vec4fGrid>()) {
-    const openvdb::Vec4fGrid::ConstPtr vel_grid = openvdb::gridConstPtrCast<openvdb::Vec4fGrid>(
-        grid);
-    extrema = openvdb::tools::extrema(vel_grid->cbeginValueOn());
-    voxel_size = vel_grid->voxelSize();
-  }
-  else {
-    assert(0);
+  if (typed_grid == nullptr) {
     return 0;
   }
 
+  const nanovdb::Vec3d voxel_size = typed_grid->voxelSize();
+
   /* We should only have uniform grids, so x = y = z, but we never know. */
-  const double max_voxel_size = openvdb::math::Max(voxel_size.x(), voxel_size.y(), voxel_size.z());
+  const double max_voxel_size = openvdb::math::Max(voxel_size[0], voxel_size[1], voxel_size[2]);
   if (max_voxel_size == 0.0) {
     return 0;
   }
 
-  const double estimated_padding = extrema.max() * static_cast<double>(velocity_scale) /
+  /* TODO: we may need to also find outliers and clamp them to avoid adding too much padding. */
+  const nanovdb::Vec3f mn = typed_grid->tree().root().minimum();
+  const nanovdb::Vec3f mx = typed_grid->tree().root().maximum();
+  float max_value = 0.0f;
+  max_value = max(max_value, fabsf(mx[0]));
+  max_value = max(max_value, fabsf(mx[1]));
+  max_value = max(max_value, fabsf(mx[2]));
+  max_value = max(max_value, fabsf(mn[0]));
+  max_value = max(max_value, fabsf(mn[1]));
+  max_value = max(max_value, fabsf(mn[2]));
+
+  const double estimated_padding = max_value * static_cast<double>(velocity_scale) /
                                    max_voxel_size;
 
   return static_cast<int>(std::ceil(estimated_padding));
 }
+#endif
 
+#ifdef WITH_OPENVDB
 static openvdb::FloatGrid::ConstPtr get_vdb_for_attribute(Volume *volume, AttributeStandard std)
 {
   Attribute *attr = volume->attributes.find(std);
@@ -619,7 +464,7 @@ static void merge_scalar_grids_for_velocity(const Scene *scene, Volume *volume)
   const ImageParams params;
   attr->data_voxel() = scene->image_manager->add_image(std::move(loader), params);
 }
-#endif
+#endif /* defined(WITH_OPENVDB) && defined(WITH_NANOVDB) */
 
 /* ************************************************************************** */
 
@@ -662,69 +507,44 @@ void GeometryManager::create_volume_mesh(const Scene *scene, Volume *volume, Pro
     return;
   }
 
+#if defined(WITH_OPENVDB) && defined(WITH_NANOVDB)
   /* Create volume mesh builder. */
   VolumeMeshBuilder builder;
 
-#ifdef WITH_OPENVDB
-  merge_scalar_grids_for_velocity(scene, volume);
-
   for (Attribute &attr : volume->attributes.attributes) {
     if (attr.element != ATTR_ELEMENT_VOXEL) {
       continue;
     }
 
-    bool do_clipping = false;
-
     ImageHandle &handle = attr.data_voxel();
 
     if (handle.empty()) {
       continue;
     }
 
-    /* Try building from OpenVDB grid directly. */
-    VDBImageLoader *vdb_loader = handle.vdb_loader();
-    openvdb::GridBase::ConstPtr grid;
-    if (vdb_loader) {
-      grid = vdb_loader->get_grid();
-
-      /* If building from an OpenVDB grid, we need to manually clip the values. */
-      do_clipping = true;
+    /* Create NanoVDB grid handle from texture memory. */
+    device_texture *texture = handle.image_memory();
+    if (texture == nullptr || texture->host_pointer == nullptr ||
+        !is_nanovdb_type(texture->info.data_type))
+    {
+      continue;
     }
 
-    /* Else fall back to creating an OpenVDB grid from the dense volume data. */
-    if (!grid) {
-      device_texture *image_memory = handle.image_memory();
+    nanovdb::GridHandle grid(
+        nanovdb::HostBuffer::createFull(texture->memory_size(), texture->host_pointer));
 
-      if (image_memory->data_elements == 1) {
-        grid = openvdb_grid_from_device_texture<openvdb::FloatGrid>(
-            image_memory, volume->get_clipping(), handle.metadata().transform_3d);
-      }
-      else if (image_memory->data_elements == 3) {
-        grid = openvdb_grid_from_device_texture<openvdb::Vec3fGrid>(
-            image_memory, volume->get_clipping(), handle.metadata().transform_3d);
-      }
-      else if (image_memory->data_elements == 4) {
-        grid = openvdb_grid_from_device_texture<openvdb::Vec4fGrid>(
-            image_memory, volume->get_clipping(), handle.metadata().transform_3d);
-      }
+    /* Add padding based on the maximum velocity vector. */
+    if (attr.std == ATTR_STD_VOLUME_VELOCITY && scene->need_motion() != Scene::MOTION_NONE) {
+      pad_size = max(pad_size,
+                     estimate_required_velocity_padding(grid, volume->get_velocity_scale()));
     }
 
-    if (grid) {
-      /* Add padding based on the maximum velocity vector. */
-      if (attr.std == ATTR_STD_VOLUME_VELOCITY && scene->need_motion() != Scene::MOTION_NONE) {
-        pad_size = max(pad_size,
-                       estimate_required_velocity_padding(grid, volume->get_velocity_scale()));
-      }
-
-      builder.add_grid(grid, do_clipping, volume->get_clipping());
-    }
+    builder.add_grid(grid);
   }
-#else
-  (void)scene;
-#endif
 
   /* If nothing to build, early out. */
   if (builder.empty_grid()) {
+    LOG(WORK) << "Memory usage volume mesh: 0 Mb. (empty grid)";
     return;
   }
 
@@ -759,6 +579,16 @@ void GeometryManager::create_volume_mesh(const Scene *scene, Volume *volume, Pro
             << (vertices.size() * sizeof(float3) + indices.size() * sizeof(int)) /
                    (1024.0 * 1024.0)
             << "Mb.";
+#else
+  (void)scene;
+#endif /* defined(WITH_OPENVDB) && defined(WITH_NANOVDB) */
+}
+
+void Volume::merge_grids(const Scene *scene)
+{
+#ifdef WITH_OPENVDB
+  merge_scalar_grids_for_velocity(scene, this);
+#endif
 }
 
 CCL_NAMESPACE_END
diff --git a/intern/cycles/scene/volume.h b/intern/cycles/scene/volume.h
index 1d15d14a3e0..7e7f542128e 100644
--- a/intern/cycles/scene/volume.h
+++ b/intern/cycles/scene/volume.h
@@ -16,11 +16,13 @@ class Volume : public Mesh {
 
   Volume();
 
-  NODE_SOCKET_API(float, clipping)
   NODE_SOCKET_API(float, step_size)
   NODE_SOCKET_API(bool, object_space)
   NODE_SOCKET_API(float, velocity_scale)
 
+  /* Merge attributes for efficiency, call right after creating them. */
+  void merge_grids(const Scene *scene);
+
   void clear(bool preserve_shaders = false) override;
 };
 
diff --git a/intern/cycles/util/texture.h b/intern/cycles/util/texture.h
index 7f74d2988de..d4a977be4c1 100644
--- a/intern/cycles/util/texture.h
+++ b/intern/cycles/util/texture.h
@@ -4,6 +4,7 @@
 
 #pragma once
 
+#include "util/defines.h"
 #include "util/transform.h"
 
 CCL_NAMESPACE_BEGIN
@@ -39,12 +40,18 @@ enum ImageDataType {
   IMAGE_DATA_TYPE_USHORT = 7,
   IMAGE_DATA_TYPE_NANOVDB_FLOAT = 8,
   IMAGE_DATA_TYPE_NANOVDB_FLOAT3 = 9,
-  IMAGE_DATA_TYPE_NANOVDB_FPN = 10,
-  IMAGE_DATA_TYPE_NANOVDB_FP16 = 11,
+  IMAGE_DATA_TYPE_NANOVDB_FLOAT4 = 10,
+  IMAGE_DATA_TYPE_NANOVDB_FPN = 11,
+  IMAGE_DATA_TYPE_NANOVDB_FP16 = 12,
 
   IMAGE_DATA_NUM_TYPES
 };
 
+ccl_device_inline bool is_nanovdb_type(int type)
+{
+  return (type >= IMAGE_DATA_TYPE_NANOVDB_FLOAT && type <= IMAGE_DATA_TYPE_NANOVDB_FP16);
+}
+
 /* Alpha types
  * How to treat alpha in images. */
 enum ImageAlphaType {
diff --git a/tests/python/cycles_render_tests.py b/tests/python/cycles_render_tests.py
index aa3e9e7fd05..fd4dfad2783 100644
--- a/tests/python/cycles_render_tests.py
+++ b/tests/python/cycles_render_tests.py
@@ -108,7 +108,6 @@ BLOCKLIST_GPU = [
     'image_log.blend',
     'glass_mix_40964.blend',
     'filter_glossy_refraction_45609.blend',
-    'smoke_color.blend',
     'bevel_mblur.blend',
     # Inconsistency between Embree and Hair primitive on GPU.
     'denoise_hair.blend',