diff --git a/include/NeoN/finiteVolume/cellCentred/interpolation/inlineInterpolationKernels.hpp b/include/NeoN/finiteVolume/cellCentred/interpolation/inlineInterpolationKernels.hpp
new file mode 100644
index 0000000000..ead163f34d
--- /dev/null
+++ b/include/NeoN/finiteVolume/cellCentred/interpolation/inlineInterpolationKernels.hpp
@@ -0,0 +1,69 @@
+// SPDX-FileCopyrightText: 2023 - 2026 NeoN authors
+//
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include <variant>
+
+#include "NeoN/core/parallelAlgorithms.hpp"
+#include "NeoN/core/primitives/scalar.hpp"
+#include "NeoN/core/primitives/label.hpp"
+#include "NeoN/core/view.hpp"
+
+namespace NeoN::finiteVolume::cellCentred
+{
+
+/* @brief Device-callable weight kernel for linear interpolation.
+** Captures the pre-computed geometry weights views (internal + boundary); no virtual dispatch.
+** The boundary view covers both physical and proc boundary faces in a single contiguous array.
+*/
+struct LinearInlineKernel
+{
+    View<const scalar> weights;
+    View<const scalar> bWeights;
+
+    NEON_INLINE_FUNCTION scalar weight(localIdx facei, scalar /*flux*/) const
+    {
+        return weights[facei];
+    }
+
+    NEON_INLINE_FUNCTION scalar boundaryWeight(localIdx bfacei, scalar /*flux*/) const
+    {
+        return bWeights[bfacei];
+    }
+
+    NEON_INLINE_FUNCTION scalar procBoundaryWeight(localIdx bcfacei, scalar /*flux*/) const
+    {
+        return bWeights[bcfacei];
+    }
+};
+
+/* @brief Device-callable weight kernel for upwind interpolation.
+** Internal weight: 1 if flux >= 0 (owner upwind), 0 otherwise.
+** Physical boundary weight: always 1 (value comes from the BC, not interpolated).
+** Proc boundary weight: flux-sign, matching the coupled internal-face convention.
+*/
+struct UpwindInlineKernel
+{
+    NEON_INLINE_FUNCTION scalar weight([[maybe_unused]] localIdx facei, scalar flux) const
+    {
+        return flux >= scalar(0) ? scalar(1) : scalar(0);
+    }
+
+    NEON_INLINE_FUNCTION scalar
+    boundaryWeight([[maybe_unused]] localIdx bfacei, [[maybe_unused]] scalar flux) const
+    {
+        return scalar(1);
+    }
+
+    NEON_INLINE_FUNCTION scalar
+    procBoundaryWeight([[maybe_unused]] localIdx bcfacei, scalar flux) const
+    {
+        return flux >= scalar(0) ? scalar(1) : scalar(0);
+    }
+};
+
+using InlineWeightKernel = std::variant<LinearInlineKernel, UpwindInlineKernel>;
+
+} // namespace NeoN::finiteVolume::cellCentred
diff --git a/include/NeoN/finiteVolume/cellCentred/interpolation/linear.hpp b/include/NeoN/finiteVolume/cellCentred/interpolation/linear.hpp
index 7c4a7648a4..36f6e8707d 100644
--- a/include/NeoN/finiteVolume/cellCentred/interpolation/linear.hpp
+++ b/include/NeoN/finiteVolume/cellCentred/interpolation/linear.hpp
@@ -35,6 +35,7 @@ void computeLinearInterpolation(
     SurfaceField<ValueType>& dst
 );
 
+
 template<typename ValueType>
 class Linear : public SurfaceInterpolationFactory<ValueType>::template Register<Linear<ValueType>>
 {
@@ -84,6 +85,12 @@ class Linear : public SurfaceInterpolationFactory<ValueType>::template Register<
     }
 
 
+    InlineWeightKernel inlineWeightKernel() const override
+    {
+        const SurfaceField<scalar>& w = geometryScheme_->weights();
+        return LinearInlineKernel {w.internalVector().view(), w.boundaryData().value().view()};
+    }
+
     std::unique_ptr<SurfaceInterpolationFactory<ValueType>> clone() const override
     {
         return std::make_unique<Linear>(*this);
diff --git a/include/NeoN/finiteVolume/cellCentred/interpolation/linearUpwind.hpp b/include/NeoN/finiteVolume/cellCentred/interpolation/linearUpwind.hpp
index e68c94898f..408a668298 100644
--- a/include/NeoN/finiteVolume/cellCentred/interpolation/linearUpwind.hpp
+++ b/include/NeoN/finiteVolume/cellCentred/interpolation/linearUpwind.hpp
@@ -159,6 +159,8 @@ class LinearUpwind :
         );
     }
 
+    InlineWeightKernel inlineWeightKernel() const override { return UpwindInlineKernel {}; }
+
     std::unique_ptr<SurfaceInterpolationFactory<ValueType>> clone() const override
     {
         return std::make_unique<LinearUpwind>(*this);
diff --git a/include/NeoN/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp b/include/NeoN/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp
index cfa4f51b73..cae6380b56 100644
--- a/include/NeoN/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp
+++ b/include/NeoN/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp
@@ -16,6 +16,7 @@
 #include "NeoN/finiteVolume/cellCentred/fields/surfaceField.hpp"
 #include "NeoN/finiteVolume/cellCentred/fields/volumeField.hpp"
 #include "NeoN/finiteVolume/cellCentred/boundary.hpp"
+#include "NeoN/finiteVolume/cellCentred/interpolation/inlineInterpolationKernels.hpp"
 
 namespace NeoN::finiteVolume::cellCentred
 {
@@ -91,6 +92,11 @@ class SurfaceInterpolationFactory :
         fill(corr.boundaryData().value(), zero<ValueType>());
     }
 
+    /* @brief Returns a device-callable weight kernel for use inside Kokkos kernels.
+     * The kernel computes the face interpolation weight inline per face without virtual dispatch.
+     */
+    virtual InlineWeightKernel inlineWeightKernel() const = 0;
+
     // Pure virtual function for cloning
     virtual std::unique_ptr<SurfaceInterpolationFactory<ValueType>> clone() const = 0;
 
@@ -159,6 +165,11 @@ class SurfaceInterpolation
 
     bool corrected() const { return interpolationKernel_->corrected(); }
 
+    InlineWeightKernel inlineWeightKernel() const
+    {
+        return interpolationKernel_->inlineWeightKernel();
+    }
+
     void correction(
         const SurfaceField<scalar>& flux,
         const VolumeField<ValueType>& src,
diff --git a/include/NeoN/finiteVolume/cellCentred/interpolation/upwind.hpp b/include/NeoN/finiteVolume/cellCentred/interpolation/upwind.hpp
index 7e940f3f6e..c4832e1dc8 100644
--- a/include/NeoN/finiteVolume/cellCentred/interpolation/upwind.hpp
+++ b/include/NeoN/finiteVolume/cellCentred/interpolation/upwind.hpp
@@ -95,6 +95,8 @@ class Upwind : public SurfaceInterpolationFactory<ValueType>::template Register<
         computeUpwindInterpolationWeights(faceFlux, src, weights);
     }
 
+    InlineWeightKernel inlineWeightKernel() const override { return UpwindInlineKernel {}; }
+
     std::unique_ptr<SurfaceInterpolationFactory<ValueType>> clone() const override
     {
         return std::make_unique<Upwind>(*this);
diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp
index 352ff03019..b0c3e09363 100644
--- a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp
+++ b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp
@@ -87,6 +87,32 @@ void computeDiv(
     );
 }
 
+template<typename WeightKernel, typename ValueType>
+void computeDivExpFusedInternal(
+    const Executor& exec,
+    localIdx nInternalFaces,
+    const WeightKernel& kernel,
+    View<const scalar> fluxV,
+    View<const localIdx> ownV,
+    View<const localIdx> neiV,
+    View<const ValueType> phiV,
+    View<ValueType> resV
+)
+{
+    parallelFor(
+        exec,
+        {0, nInternalFaces},
+        NEON_LAMBDA(const localIdx i) {
+            const auto w = kernel.weight(i, fluxV[i]);
+            const ValueType phiF = w * phiV[ownV[i]] + (scalar(1) - w) * phiV[neiV[i]];
+            const ValueType flux = fluxV[i] * phiF;
+            Kokkos::atomic_add(&resV[ownV[i]], flux);
+            Kokkos::atomic_sub(&resV[neiV[i]], flux);
+        },
+        "sumFluxesInternal"
+    );
+}
+
 template<typename ValueType>
 void computeDivExp(
     const SurfaceField<scalar>& faceFlux,
@@ -98,41 +124,77 @@ void computeDivExp(
 {
     const UnstructuredMesh& mesh = phi.mesh();
     const auto exec = phi.exec();
-    SurfaceField<ValueType> phif(
-        exec, "phif", mesh, createCalculatedBCs<SurfaceBoundary<ValueType>>(mesh)
+    const auto nInternalFaces = mesh.nInternalFaces();
+    const auto nBoundaryFaces = mesh.nBoundaryFaces();
+
+    // Corrected schemes (e.g. linearUpwind) compute interpolate() = weight*phi + correction in one
+    // pass; fall back to allocating phif so correction() is included in the face value.
+    if (surfInterp.corrected())
+    {
+        SurfaceField<ValueType> phif(
+            exec, "phif", mesh, createCalculatedBCs<SurfaceBoundary<ValueType>>(mesh)
+        );
+        surfInterp.interpolate(faceFlux, phi, phif);
+        phif.boundaryData().value() = phi.boundaryData().value();
+        computeDiv<ValueType>(
+            exec,
+            nInternalFaces,
+            nBoundaryFaces,
+            mesh.faceNeighbors().view(),
+            mesh.faceOwners().view(),
+            mesh.boundaryMesh().faceOwners().view(),
+            faceFlux.internalVector().view(),
+            faceFlux.boundaryData().value().view(),
+            phif.internalVector().view(),
+            phif.boundaryData().value().view(),
+            mesh.cellVolumes().view(),
+            divPhi.view(),
+            operatorScaling
+        );
+        return;
+    }
+
+    // Non-corrected schemes: fuse face-value computation into the divergence kernels,
+    // eliminating the temporary phif SurfaceField allocation.
+    const auto [fluxV, ownV, neiV] =
+        views(faceFlux.internalVector(), mesh.faceOwners(), mesh.faceNeighbors());
+    const auto [bFluxV, bFaceOwnersV] =
+        views(faceFlux.boundaryData().value(), mesh.boundaryMesh().faceOwners());
+    const auto phiV = phi.internalVector().view();
+    const auto bPhiV = phi.boundaryData().value().view();
+    const auto volV = mesh.cellVolumes().view();
+    auto resV = divPhi.view();
+
+    std::visit(
+        [&](auto&& kernel) {
+            computeDivExpFusedInternal(exec, nInternalFaces, kernel, fluxV, ownV, neiV, phiV, resV);
+        },
+        surfInterp.inlineWeightKernel()
     );
-    // TODO: remove or implement
-    // fill(phif.internalVector(), NeoN::zero<ValueType>::value);
-    surfInterp.interpolate(faceFlux, phi, phif);
 
-    // TODO: currently we just copy the boundary values over
-    phif.boundaryData().value() = phi.boundaryData().value();
+    parallelFor(
+        exec,
+        {0, nBoundaryFaces},
+        NEON_LAMBDA(const localIdx bfi) {
+            Kokkos::atomic_add(&resV[bFaceOwnersV[bfi]], bFluxV[bfi] * bPhiV[bfi]);
+        },
+        "sumFluxesBoundary"
+    );
 
-    auto nInternalFaces = mesh.nInternalFaces();
-    auto nBoundaryFaces = mesh.nBoundaryFaces();
-    computeDiv<ValueType>(
+    parallelFor(
         exec,
-        nInternalFaces,
-        nBoundaryFaces,
-        mesh.faceNeighbors().view(),
-        mesh.faceOwners().view(),
-        mesh.boundaryMesh().faceOwners().view(),
-        faceFlux.internalVector().view(),
-        faceFlux.boundaryData().value().view(),
-        phif.internalVector().view(),
-        phif.boundaryData().value().view(),
-        mesh.cellVolumes().view(),
-        divPhi.view(),
-        operatorScaling
+        {0, static_cast<localIdx>(volV.size())},
+        NEON_LAMBDA(const localIdx celli) { resV[celli] *= operatorScaling[celli] / volV[celli]; },
+        "normalizeFluxes"
     );
 }
 
-template<typename FieldValueType, typename AssemblyType = FieldValueType>
+template<typename FieldValueType, typename AssemblyType = FieldValueType, typename WeightKernel>
 void computeDivProcBoundImpl(
     la::LinearSystem<AssemblyType, FieldValueType>& ls,
     const SurfaceField<scalar>& faceFlux,
     const VolumeField<FieldValueType>& phi,
-    const SurfaceField<scalar>& weights,
+    WeightKernel weightKernel,
     const dsl::Coeff coeff
 )
 {
@@ -147,7 +209,6 @@ void computeDivProcBoundImpl(
         views(mesh.boundaryMesh().faceOwners(), mesh.boundaryMesh().weights());
 
     const auto bFluxV = faceFlux.boundaryData().value().view();
-    const auto bWeightsV = weights.boundaryData().value().view();
     auto bValues = ls.offDiagonalMatrix().values().view();
     // boundaryMatrix records the diagonal contribution so removeBoundaryContributions can reverse
     // it (proc slots live at [nBoundaryFaces, nBoundaryFaces + nProcBoundaryFaces)).
@@ -175,7 +236,8 @@ void computeDivProcBoundImpl(
             auto isOwnerFace = isOwner[bcfacei] > 0.0;
             auto sign = isOwnerFace ? scalar(-1) : scalar(1);
 
-            auto weight = isOwnerFace ? bWeightsV[bcfacei] : (scalar(1) - bWeightsV[bcfacei]);
+            auto bw = weightKernel.procBoundaryWeight(bcfacei, bFluxV[bcfacei]);
+            auto weight = isOwnerFace ? bw : (scalar(1) - bw);
             auto fluxContrib = sign * weight * bFluxV[bcfacei] * ownCoeff * one<AssemblyType>();
 
             Kokkos::atomic_sub(&values[ma.diagIdx(cell)], fluxContrib);
@@ -189,12 +251,12 @@ void computeDivProcBoundImpl(
 }
 
 
-template<typename FieldValueType, typename AssemblyType = FieldValueType>
+template<typename FieldValueType, typename AssemblyType = FieldValueType, typename WeightKernel>
 void computeDivBoundImpl(
     la::LinearSystem<AssemblyType, FieldValueType>& ls,
     const SurfaceField<scalar>& faceFlux,
     const VolumeField<FieldValueType>& phi,
-    const SurfaceField<scalar>& weights,
+    WeightKernel weightKernel,
     const dsl::Coeff operatorScaling
 )
 {
@@ -208,9 +270,8 @@ void computeDivBoundImpl(
 
     auto values = ls.matrix().values().view();
 
-    auto [bFaceFluxV, bweights, refGradient, valueFraction, refValue] = views(
+    auto [bFaceFluxV, refGradient, valueFraction, refValue] = views(
         faceFlux.boundaryData().value(),
-        weights.boundaryData().value(),
         phi.boundaryData().refGrad(),
         phi.boundaryData().valueFraction(),
         phi.boundaryData().refValue()
@@ -232,8 +293,8 @@ void computeDivBoundImpl(
             auto refValFrac = valueFraction[bfi];
             auto refGradFrac = 1.0 - refValFrac;
 
-            auto flux =
-                bFaceFluxV[bfi] * -bweights[bfi] * ownCoeff * refGradFrac * one<AssemblyType>();
+            auto bw = weightKernel.boundaryWeight(bfi, bFaceFluxV[bfi]);
+            auto flux = bFaceFluxV[bfi] * -bw * ownCoeff * refGradFrac * one<AssemblyType>();
 
             // since upper triangular value is "outside" of system matrix
             // it is stored separately in bMatrix
@@ -245,11 +306,10 @@ void computeDivBoundImpl(
             //   φ_f = refValFrac * refValue               (Dirichlet part)
             //       + refGradFrac * (φ_C + refGradient/δ)  (Neumann part)
             // The implicit valFrac2 * φ_C term is handled via fluxContrib above.
-            // bweights converts the Dirichlet face value to a cell-to-face flux contribution;
+            // bw converts the Dirichlet face value to a cell-to-face flux contribution;
             // the Neumann gradient correction (refGradient/δ) enters directly as a known increment.
-            auto valueRhs =
-                (bweights[bfi] * bFaceFluxV[bfi] * ownCoeff * (refValFrac * refValue[bfi]))
-                + refGradFrac * refGradient[bfi] * (1 / deltaCoeffs[bfi]);
+            auto valueRhs = (bw * bFaceFluxV[bfi] * ownCoeff * (refValFrac * refValue[bfi]))
+                          + refGradFrac * refGradient[bfi] * (1 / deltaCoeffs[bfi]);
             Kokkos::atomic_sub(&rhs[ownRow], valueRhs);
             bRhs[bfi] += valueRhs;
         },
@@ -258,12 +318,12 @@ void computeDivBoundImpl(
 }
 
 
-template<typename FieldValueType, typename AssemblyType = FieldValueType>
+template<typename FieldValueType, typename AssemblyType = FieldValueType, typename WeightKernel>
 void computeDivIntImp(
     la::LinearSystem<AssemblyType, FieldValueType>& ls,
     const SurfaceField<scalar>& faceFlux,
     const VolumeField<FieldValueType>& phi,
-    const SurfaceField<scalar>& weights,
+    WeightKernel weightKernel,
     const dsl::Coeff coeff
 )
 {
@@ -273,9 +333,8 @@ void computeDivIntImp(
 
     const auto ma = ls.faceToMatrixAddress()->view(ls.matrix().sparsity()->rowOffs().view());
 
-    const auto [fluxV, weightsV, ownV, neiV, surfFaceCells] = views(
+    const auto [fluxV, ownV, neiV, surfFaceCells] = views(
         faceFlux.internalVector(),
-        weights.internalVector(),
         mesh.faceOwners(),
         mesh.faceNeighbors(),
         mesh.boundaryMesh().faceOwners()
@@ -301,8 +360,9 @@ void computeDivIntImp(
             // Decompose face flux via linear interpolation:
             //   ownFluxContrib = w * F_f     — part attributed to the owner cell value
             //   neiFluxContrib = (1-w) * F_f — part attributed to the neighbor cell value
-            auto ownFluxContrib = -fluxV[facei] * weightsV[facei] * one<AssemblyType>();
-            auto neiFluxContrib = +fluxV[facei] * (1.0 - weightsV[facei]) * one<AssemblyType>();
+            const auto w = weightKernel.weight(facei, fluxV[facei]);
+            auto ownFluxContrib = -fluxV[facei] * w * one<AssemblyType>();
+            auto neiFluxContrib = +fluxV[facei] * (1.0 - w) * one<AssemblyType>();
 
             // triangular coefficients - neighbor -> lower, owner -> upper
             values[ma.lowerIdx(neiRow, facei)] += ownFluxContrib * neiCoeff;
@@ -316,12 +376,12 @@ void computeDivIntImp(
     );
 };
 
-template<typename FieldValueType, typename AssemblyType = FieldValueType>
+template<typename FieldValueType, typename AssemblyType = FieldValueType, typename WeightKernel>
 void computeDivIntCellBasedImp(
     la::LinearSystem<AssemblyType, FieldValueType>& ls,
     const SurfaceField<scalar>& faceFlux,
     const VolumeField<FieldValueType>& phi,
-    const SurfaceField<scalar>& weights,
+    WeightKernel weightKernel,
     const dsl::Coeff coeff
 )
 {
@@ -330,7 +390,7 @@ void computeDivIntCellBasedImp(
     const auto ma = ls.faceToMatrixAddress()->view(ls.matrix().sparsity()->rowOffs().view());
     auto iterator = std::dynamic_pointer_cast<la::CellBasedIterator>(ls.getMeshIterator()->get());
 
-    const auto [fluxV, weightsV] = views(faceFlux.internalVector(), weights.internalVector());
+    const auto fluxV = faceFlux.internalVector().view();
 
     auto cellBasedData = iterator->getCellBasedData();
     NF_ASSERT(
@@ -358,7 +418,7 @@ void computeDivIntCellBasedImp(
                 const auto faceIdx = cellFacesValues[startIdx + i];
                 const auto sign = faceSignV[startIdx + i];
                 const auto flux = fluxV[faceIdx];
-                const auto w = weightsV[faceIdx];
+                const auto w = weightKernel.weight(faceIdx, flux);
 
                 AssemblyType offDiag;
                 AssemblyType diagContrib;
@@ -510,24 +570,30 @@ void GaussGreenDiv<FieldValueType, AssemblyType>::div(
     const dsl::Coeff operatorScaling
 ) const
 {
-    // TODO boundary weights should be separate so that we can start internal assembly
-    const auto weights = surfaceInterpolation_.weight(faceFlux, phi);
-    if (auto* cellIter = dynamic_cast<la::CellBasedIterator*>(ls.getMeshIterator()->get().get()))
-    {
-        if (!cellIter->getCellBasedData())
+    const auto inlineKernel = surfaceInterpolation_.inlineWeightKernel();
+    std::visit(
+        [&](auto&& kernel)
         {
-            cellIter->setComputeCellBasedData(
-                phi.mesh(), ls.matrix().sparsity(), ls.faceToMatrixAddress()
-            );
-        }
-        computeDivIntCellBasedImp(ls, faceFlux, phi, weights, operatorScaling);
-    }
-    else
-    {
-        computeDivIntImp(ls, faceFlux, phi, weights, operatorScaling);
-    }
-    computeDivBoundImpl(ls, faceFlux, phi, weights, operatorScaling);
-    computeDivProcBoundImpl(ls, faceFlux, phi, weights, operatorScaling);
+            if (auto* cellIter =
+                    dynamic_cast<la::CellBasedIterator*>(ls.getMeshIterator()->get().get()))
+            {
+                if (!cellIter->getCellBasedData())
+                {
+                    cellIter->setComputeCellBasedData(
+                        phi.mesh(), ls.matrix().sparsity(), ls.faceToMatrixAddress()
+                    );
+                }
+                computeDivIntCellBasedImp(ls, faceFlux, phi, kernel, operatorScaling);
+            }
+            else
+            {
+                computeDivIntImp(ls, faceFlux, phi, kernel, operatorScaling);
+            }
+            computeDivBoundImpl(ls, faceFlux, phi, kernel, operatorScaling);
+            computeDivProcBoundImpl(ls, faceFlux, phi, kernel, operatorScaling);
+        },
+        inlineKernel
+    );
 
     // Deferred correction for corrected schemes (e.g. linearUpwind): the implicit matrix uses the
     // upwind weights above, the gradient correction is added explicitly to the rhs.
diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenDivLaplacian.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDivLaplacian.cpp
index 0851d4c0e0..9dff3c104f 100644
--- a/src/finiteVolume/cellCentred/operators/gaussGreenDivLaplacian.cpp
+++ b/src/finiteVolume/cellCentred/operators/gaussGreenDivLaplacian.cpp
@@ -14,14 +14,13 @@
 namespace NeoN::finiteVolume::cellCentred
 {
 
-template<typename FieldValueType, typename AssemblyType>
+template<typename FieldValueType, typename AssemblyType, typename WeightKernel>
 static void computeDivLaplacianIntImpl(
     la::LinearSystem<AssemblyType, FieldValueType>& ls,
-    const VolumeField<FieldValueType>& /*U*/,
     const SurfaceField<scalar>& phi,
     const SurfaceField<scalar>& gamma,
-    const SurfaceInterpolation<FieldValueType>& /*divSurfInterp*/,
     const FaceNormalGradient<FieldValueType>& faceNormalGradient,
+    WeightKernel weightKernel,
     const dsl::Coeff coeffA,
     const dsl::Coeff coeffB
 )
@@ -46,7 +45,7 @@ static void computeDivLaplacianIntImpl(
             auto nei = neiV[facei];
 
             auto fluxDiv = phiV[facei];
-            const auto weight = (phiV[facei] >= 0) ? 1.0 : 0.0;
+            const auto weight = weightKernel.weight(facei, fluxDiv);
             auto fluxLap = deltaV[facei] * gammaV[facei] * magFaceAreaV[facei];
 
             auto coeffNeiA = coeffA[nei];
@@ -70,14 +69,13 @@ static void computeDivLaplacianIntImpl(
     );
 }
 
-template<typename FieldValueType, typename AssemblyType>
+template<typename FieldValueType, typename AssemblyType, typename WeightKernel>
 static void computeDivLaplacianIntCellBasedImpl(
     la::LinearSystem<AssemblyType, FieldValueType>& ls,
-    const VolumeField<FieldValueType>& /*U*/,
     const SurfaceField<scalar>& phi,
     const SurfaceField<scalar>& gamma,
-    const SurfaceInterpolation<FieldValueType>& /*divSurfInterp*/,
     const FaceNormalGradient<FieldValueType>& faceNormalGradient,
+    WeightKernel weightKernel,
     const dsl::Coeff coeffA,
     const dsl::Coeff coeffB
 )
@@ -122,9 +120,8 @@ static void computeDivLaplacianIntCellBasedImpl(
                 const auto faceIdx = cellFacesValues[startIdx + i];
                 const auto sign = faceSignV[startIdx + i];
 
-                // upwind weight: 1 if flux leaves owner (flux > 0), 0 otherwise
                 const auto flux = phiV[faceIdx];
-                const auto w = (flux >= 0) ? scalar(1) : scalar(0);
+                const auto w = weightKernel.weight(faceIdx, flux);
 
                 // Laplacian face coefficient: δ_f · γ_f · |S_f|
                 const auto fluxLap = deltaV[faceIdx] * gammaV[faceIdx] * magFaceAreaV[faceIdx];
@@ -156,14 +153,14 @@ static void computeDivLaplacianIntCellBasedImpl(
     );
 }
 
-template<typename FieldValueType, typename AssemblyType>
+template<typename FieldValueType, typename AssemblyType, typename WeightKernel>
 static void computeDivLaplacianBoundImpl(
     la::LinearSystem<AssemblyType, FieldValueType>& ls,
     const VolumeField<FieldValueType>& u,
     const SurfaceField<scalar>& phi,
     const SurfaceField<scalar>& gamma,
-    const SurfaceInterpolation<FieldValueType>& /*divSurfInterp*/,
     const FaceNormalGradient<FieldValueType>& faceNormalGradient,
+    WeightKernel weightKernel,
     const dsl::Coeff coeffA,
     const dsl::Coeff coeffB
 )
@@ -208,9 +205,9 @@ static void computeDivLaplacianBoundImpl(
             auto valFrac1 = valueFraction[bfi];
             auto valFrac2 = 1.0 - valFrac1;
 
-            const auto bweights = (fluxDiv >= 0) ? scalar(1) : scalar(0);
+            const auto bw = weightKernel.boundaryWeight(bfi, fluxDiv);
 
-            auto valueDiv = -bweights * coeffAOwn * fluxDiv * valFrac2;
+            auto valueDiv = -bw * coeffAOwn * fluxDiv * valFrac2;
             auto valueLap = bDeltaV[bfi] * coeffBOwn * fluxLap * valFrac1;
             auto valueA = (valueDiv + valueLap) * one<AssemblyType>();
 
@@ -233,14 +230,13 @@ static void computeDivLaplacianBoundImpl(
     );
 }
 
-template<typename FieldValueType, typename AssemblyType>
+template<typename FieldValueType, typename AssemblyType, typename WeightKernel>
 static void computeDivLaplacianProcBoundImpl(
     la::LinearSystem<AssemblyType, FieldValueType>& ls,
-    const VolumeField<FieldValueType>& /*U*/,
     const SurfaceField<scalar>& phi,
     const SurfaceField<scalar>& gamma,
-    const SurfaceInterpolation<FieldValueType>& /*divSurfInterp*/,
     const FaceNormalGradient<FieldValueType>& faceNormalGradient,
+    WeightKernel weightKernel,
     const dsl::Coeff coeffA,
     const dsl::Coeff coeffB
 )
@@ -281,12 +277,12 @@ static void computeDivLaplacianProcBoundImpl(
             auto lapFlux = bGammaV[bcfacei] * bMagSf[bcfacei] * bDeltaCoeffs[bcfacei];
             auto lapValue = lapFlux * ownCoeffB * one<AssemblyType>();
 
-            // Div upwind contribution
+            // Div contribution
             auto isOwnerFace = isOwner[bcfacei] > 0.0;
             auto sign = isOwnerFace ? scalar(-1) : scalar(1);
             auto bFlux = bPhiV[bcfacei];
-            auto weight = isOwnerFace ? (bFlux >= 0 ? scalar(1) : scalar(0))
-                                      : (bFlux >= 0 ? scalar(0) : scalar(1));
+            auto bw = weightKernel.procBoundaryWeight(bcfacei, bFlux);
+            auto weight = isOwnerFace ? bw : (scalar(1) - bw);
             auto divDiag = sign * weight * bFlux * ownCoeffA * one<AssemblyType>();
             auto divOffDiag =
                 -sign * (scalar(1) - weight) * bFlux * ownCoeffA * one<AssemblyType>();
@@ -323,59 +319,39 @@ void GaussGreenDivLaplacian<ValueType>::explicitOperation(Vector<ValueType>& /*s
 template<typename ValueType>
 void GaussGreenDivLaplacian<ValueType>::implicitOperation(la::LinearSystem<ValueType>& ls) const
 {
-    if (auto* cellIter = dynamic_cast<la::CellBasedIterator*>(ls.getMeshIterator()->get().get()))
-    {
-        if (!cellIter->getCellBasedData())
+    const auto inlineKernel = divSurfaceInterpolation_->inlineWeightKernel();
+    std::visit(
+        [&](auto&& kernel)
         {
-            cellIter->setComputeCellBasedData(
-                this->getVector().mesh(), ls.matrix().sparsity(), ls.faceToMatrixAddress()
+            if (auto* cellIter =
+                    dynamic_cast<la::CellBasedIterator*>(ls.getMeshIterator()->get().get()))
+            {
+                if (!cellIter->getCellBasedData())
+                {
+                    cellIter->setComputeCellBasedData(
+                        this->getVector().mesh(), ls.matrix().sparsity(), ls.faceToMatrixAddress()
+                    );
+                }
+                computeDivLaplacianIntCellBasedImpl(
+                    ls, flux_, gamma_, *faceNormalGradient_, kernel, coeffA_, coeffB_
+                );
+            }
+            else
+            {
+                computeDivLaplacianIntImpl(
+                    ls, flux_, gamma_, *faceNormalGradient_, kernel, coeffA_, coeffB_
+                );
+            }
+            computeDivLaplacianBoundImpl(
+                ls, this->getVector(), flux_, gamma_, *faceNormalGradient_, kernel, coeffA_, coeffB_
             );
-        }
-        computeDivLaplacianIntCellBasedImpl(
-            ls,
-            this->getVector(),
-            flux_,
-            gamma_,
-            *divSurfaceInterpolation_,
-            *faceNormalGradient_,
-            coeffA_,
-            coeffB_
-        );
-    }
-    else
-    {
-        computeDivLaplacianIntImpl(
-            ls,
-            this->getVector(),
-            flux_,
-            gamma_,
-            *divSurfaceInterpolation_,
-            *faceNormalGradient_,
-            coeffA_,
-            coeffB_
-        );
-    }
-    computeLaplacianNonOrthCorrImpl(ls, gamma_, this->getVector(), coeffB_, *faceNormalGradient_);
-    computeDivLaplacianBoundImpl(
-        ls,
-        this->getVector(),
-        flux_,
-        gamma_,
-        *divSurfaceInterpolation_,
-        *faceNormalGradient_,
-        coeffA_,
-        coeffB_
-    );
-    computeDivLaplacianProcBoundImpl(
-        ls,
-        this->getVector(),
-        flux_,
-        gamma_,
-        *divSurfaceInterpolation_,
-        *faceNormalGradient_,
-        coeffA_,
-        coeffB_
+            computeDivLaplacianProcBoundImpl(
+                ls, flux_, gamma_, *faceNormalGradient_, kernel, coeffA_, coeffB_
+            );
+        },
+        inlineKernel
     );
+    computeLaplacianNonOrthCorrImpl(ls, gamma_, this->getVector(), coeffB_, *faceNormalGradient_);
 
     // Deferred correction for a corrected div scheme (e.g. linearUpwind): the fused matrix uses
     // upwind div weights, the gradient correction is added explicitly to the rhs (internal + proc).
@@ -398,61 +374,41 @@ void GaussGreenDivLaplacian<ValueType>::implicitOperation(la::LinearSystem<scala
 ) const
     requires(!std::is_same_v<ValueType, scalar>)
 {
-    if (auto* cellIter = dynamic_cast<la::CellBasedIterator*>(ls.getMeshIterator()->get().get()))
-    {
-        if (!cellIter->getCellBasedData())
+    const auto inlineKernel = divSurfaceInterpolation_->inlineWeightKernel();
+    std::visit(
+        [&](auto&& kernel)
         {
-            cellIter->setComputeCellBasedData(
-                this->getVector().mesh(), ls.matrix().sparsity(), ls.faceToMatrixAddress()
+            if (auto* cellIter =
+                    dynamic_cast<la::CellBasedIterator*>(ls.getMeshIterator()->get().get()))
+            {
+                if (!cellIter->getCellBasedData())
+                {
+                    cellIter->setComputeCellBasedData(
+                        this->getVector().mesh(), ls.matrix().sparsity(), ls.faceToMatrixAddress()
+                    );
+                }
+                computeDivLaplacianIntCellBasedImpl<ValueType, scalar>(
+                    ls, flux_, gamma_, *faceNormalGradient_, kernel, coeffA_, coeffB_
+                );
+            }
+            else
+            {
+                computeDivLaplacianIntImpl<ValueType, scalar>(
+                    ls, flux_, gamma_, *faceNormalGradient_, kernel, coeffA_, coeffB_
+                );
+            }
+            computeDivLaplacianBoundImpl<ValueType, scalar>(
+                ls, this->getVector(), flux_, gamma_, *faceNormalGradient_, kernel, coeffA_, coeffB_
             );
-        }
-        computeDivLaplacianIntCellBasedImpl<ValueType, scalar>(
-            ls,
-            this->getVector(),
-            flux_,
-            gamma_,
-            *divSurfaceInterpolation_,
-            *faceNormalGradient_,
-            coeffA_,
-            coeffB_
-        );
-    }
-    else
-    {
-        computeDivLaplacianIntImpl<ValueType, scalar>(
-            ls,
-            this->getVector(),
-            flux_,
-            gamma_,
-            *divSurfaceInterpolation_,
-            *faceNormalGradient_,
-            coeffA_,
-            coeffB_
-        );
-    }
+            computeDivLaplacianProcBoundImpl<ValueType, scalar>(
+                ls, flux_, gamma_, *faceNormalGradient_, kernel, coeffA_, coeffB_
+            );
+        },
+        inlineKernel
+    );
     computeLaplacianNonOrthCorrImpl<ValueType, scalar>(
         ls, gamma_, this->getVector(), coeffB_, *faceNormalGradient_
     );
-    computeDivLaplacianBoundImpl<ValueType, scalar>(
-        ls,
-        this->getVector(),
-        flux_,
-        gamma_,
-        *divSurfaceInterpolation_,
-        *faceNormalGradient_,
-        coeffA_,
-        coeffB_
-    );
-    computeDivLaplacianProcBoundImpl<ValueType, scalar>(
-        ls,
-        this->getVector(),
-        flux_,
-        gamma_,
-        *divSurfaceInterpolation_,
-        *faceNormalGradient_,
-        coeffA_,
-        coeffB_
-    );
 
     // Deferred correction for a corrected div scheme (e.g. linearUpwind), scalar-matrix / Vec3-rhs.
     if (divSurfaceInterpolation_->corrected())
diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenGrad.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenGrad.cpp
index 1bb8bbaf74..87cb1a7f41 100644
--- a/src/finiteVolume/cellCentred/operators/gaussGreenGrad.cpp
+++ b/src/finiteVolume/cellCentred/operators/gaussGreenGrad.cpp
@@ -17,19 +17,16 @@ namespace NeoN::finiteVolume::cellCentred
 ** @param[in] in - Vector on which the gradient should be computed
 ** @param[in,out] out - Vector to hold the result
 */
+template<typename WeightKernel>
 void computeGrad(
     const VolumeField<scalar>& in,
-    const SurfaceInterpolation<scalar>& surfInterp,
+    const WeightKernel& wKernel,
     Vector<Vec3>& out,
     const dsl::Coeff operatorScaling
 )
 {
     const UnstructuredMesh& mesh = in.mesh();
     const auto exec = out.exec();
-    SurfaceField<scalar> phif(
-        exec, "phif", mesh, createCalculatedBCs<SurfaceBoundary<scalar>>(mesh)
-    );
-    surfInterp.interpolate(in, phif);
 
     auto surfGradPhi = out.view();
 
@@ -40,54 +37,57 @@ void computeGrad(
         mesh.faceNormals(),
         mesh.cellVolumes()
     );
-    const auto surfPhif = phif.internalVector().view();
-    const auto surfPhifB = phif.boundaryData().value().view();
+    const auto phiV = in.internalVector().view();
+    const auto bPhiV = in.boundaryData().value().view();
 
-    auto nInternalFaces = mesh.nInternalFaces();
-    auto nBoundaryFaces = mesh.nBoundaryFaces();
+    const auto nInternalFaces = mesh.nInternalFaces();
+    const auto nBoundaryFaces = mesh.nBoundaryFaces();
 
     // Green-Gauss gradient theorem: ∇φ_C = (1/V_C) * sum_f S_f * φ_f
     //
-    // S_f points from owner to neighbour by construction (valid for all internal faces).
-    //   owner cell:     S_f is the outward area vector  →  +S_f * φ_f  (add)
-    //   neighbour cell: S_f points inward to neighbour  → −S_f * φ_f  (subtract)
-    // TODO use NeoN::atomic_
+    // S_f points from owner to neighbour by construction.
+    //   owner cell:     +S_f * φ_f  (add)
+    //   neighbour cell: −S_f * φ_f  (subtract)
     parallelFor(
         exec,
         {0, nInternalFaces},
         NEON_LAMBDA(const localIdx i) {
-            Vec3 flux = faceNormals[i] * surfPhif[i];
-            Kokkos::atomic_add(&surfGradPhi[faceOwners[i]], flux);    // +S_f * φ_f
-            Kokkos::atomic_sub(&surfGradPhi[faceNeighbors[i]], flux); // −S_f * φ_f
+            const auto w = wKernel.weight(i, scalar(0));
+            const scalar phiF = w * phiV[faceOwners[i]] + (scalar(1) - w) * phiV[faceNeighbors[i]];
+            const Vec3 flux = faceNormals[i] * phiF;
+            Kokkos::atomic_add(&surfGradPhi[faceOwners[i]], flux);
+            Kokkos::atomic_sub(&surfGradPhi[faceNeighbors[i]], flux);
         },
         "computeGradInternal"
     );
 
-    // Boundary faces: only the owner cell is on this rank.
+    // Physical boundary: linear interpolation gives w_b * phi_bc at the face.
     const auto bFaceNormals = mesh.boundaryMesh().faceNormals().view();
     parallelFor(
         exec,
         {0, nBoundaryFaces},
         NEON_LAMBDA(const localIdx bfi) {
-            auto own = boundaryFaceOwners[bfi];
-            Vec3 valueOwn = bFaceNormals[bfi] * surfPhifB[bfi];
-            Kokkos::atomic_add(&surfGradPhi[own], valueOwn);
+            const auto own = boundaryFaceOwners[bfi];
+            const auto w = wKernel.boundaryWeight(bfi, scalar(0));
+            Kokkos::atomic_add(&surfGradPhi[own], bFaceNormals[bfi] * (w * bPhiV[bfi]));
         },
         "computeGradBoundary"
     );
 
-    auto nProcBoundaryFaces = mesh.nProcBoundaryFaces();
+    const auto nProcBoundaryFaces = mesh.nProcBoundaryFaces();
     if (nProcBoundaryFaces > 0)
     {
+        // Proc boundary: interpolate between local owner and ghost cell value.
         const auto bSf = mesh.boundaryMesh().faceNormals().view();
         parallelFor(
             exec,
             {0, nProcBoundaryFaces},
             NEON_LAMBDA(const localIdx procFacei) {
-                auto bcfacei = nBoundaryFaces + procFacei;
-                auto own = boundaryFaceOwners[bcfacei];
-                Vec3 flux = bSf[bcfacei] * surfPhifB[bcfacei];
-                Kokkos::atomic_add(&surfGradPhi[own], flux);
+                const auto bcfacei = nBoundaryFaces + procFacei;
+                const auto own = boundaryFaceOwners[bcfacei];
+                const auto w = wKernel.procBoundaryWeight(bcfacei, scalar(0));
+                const scalar phiF = w * phiV[own] + (scalar(1) - w) * bPhiV[bcfacei];
+                Kokkos::atomic_add(&surfGradPhi[own], bSf[bcfacei] * phiF);
             },
             "computeProcGradBoundary"
         );
@@ -203,7 +203,10 @@ void GaussGreenGrad::grad(
     const VolumeField<scalar>& phi, const dsl::Coeff operatorScaling, Vector<Vec3>& gradPhi
 ) const
 {
-    computeGrad(phi, surfaceInterpolation_, gradPhi, operatorScaling);
+    std::visit(
+        [&](auto&& kernel) { computeGrad(phi, kernel, gradPhi, operatorScaling); },
+        surfaceInterpolation_.inlineWeightKernel()
+    );
 };
 
 void GaussGreenGrad::grad(
@@ -211,8 +214,10 @@ void GaussGreenGrad::grad(
 ) const
 {
     fill(gradPhi.internalVector(), zero<Vec3>());
-
-    computeGrad(phi, surfaceInterpolation_, gradPhi.internalVector(), operatorScaling);
+    std::visit(
+        [&](auto&& kernel) { computeGrad(phi, kernel, gradPhi.internalVector(), operatorScaling); },
+        surfaceInterpolation_.inlineWeightKernel()
+    );
     computeBoundaryGrad(phi, gradPhi, operatorScaling);
 }
 
@@ -225,7 +230,10 @@ GaussGreenGrad::grad(const VolumeField<scalar>& phi, const dsl::Coeff operatorSc
     auto gradBCs = createCalculatedProcBCs<VolumeBoundary<Vec3>>(phi.mesh());
     VolumeField<Vec3> gradPhi = VolumeField<Vec3>(phi.exec(), "gradPhi", phi.mesh(), gradBCs);
     fill(gradPhi.internalVector(), zero<Vec3>());
-    computeGrad(phi, surfaceInterpolation_, gradPhi.internalVector(), operatorScaling);
+    std::visit(
+        [&](auto&& kernel) { computeGrad(phi, kernel, gradPhi.internalVector(), operatorScaling); },
+        surfaceInterpolation_.inlineWeightKernel()
+    );
     computeBoundaryGrad(phi, gradPhi, operatorScaling);
     return gradPhi;
 }
@@ -244,9 +252,10 @@ void atomicSubTensor(Tensor* target, size_t row, size_t col, scalar value)
     Kokkos::atomic_sub(&(*target)(row, col), value);
 }
 
+template<typename WeightKernel>
 void computeGradTensor(
     const VolumeField<Vec3>& u,
-    const SurfaceInterpolation<Vec3>& surfInterpVec,
+    const WeightKernel& wKernel,
     Vector<Tensor>& gradU,
     const dsl::Coeff operatorScaling
 )
@@ -254,9 +263,6 @@ void computeGradTensor(
     const UnstructuredMesh& mesh = u.mesh();
     const auto exec = gradU.exec();
 
-    SurfaceField<Vec3> uf(exec, "Uf", mesh, createCalculatedBCs<SurfaceBoundary<Vec3>>(mesh));
-    surfInterpVec.interpolate(u, uf);
-
     auto gT = gradU.view();
 
     const auto [owner, nei, SfAll, V, bFaceCells] = views(
@@ -266,8 +272,8 @@ void computeGradTensor(
         mesh.cellVolumes(),
         mesh.boundaryMesh().faceOwners()
     );
-    const auto ufInt = uf.internalVector().view();
-    const auto ufBound = uf.boundaryData().value().view();
+    const auto uInt = u.internalVector().view();
+    const auto uBound = u.boundaryData().value().view();
 
     const localIdx nInt = mesh.nInternalFaces();
     const localIdx nBnd = mesh.nBoundaryFaces();
@@ -277,9 +283,10 @@ void computeGradTensor(
         {0, nInt},
         NEON_LAMBDA(const localIdx f) {
             const Vec3 sf = SfAll[f];
-            const Vec3 faceU = ufInt[f];
             const auto o = owner[f];
             const auto n = nei[f];
+            const auto w = wKernel.weight(f, scalar(0));
+            const Vec3 faceU = w * uInt[o] + (scalar(1) - w) * uInt[n];
             // gradU(row,col) += Sf[col] * U[row]  (Gauss-Green)
             for (size_t row = 0; row < 3; ++row)
             {
@@ -294,6 +301,7 @@ void computeGradTensor(
         "computeGradTensorInternal"
     );
 
+    // Physical boundary: linear interpolation gives w_b * u_bc at the face.
     const auto bFaceNormals = mesh.boundaryMesh().faceNormals().view();
     parallelFor(
         exec,
@@ -301,7 +309,8 @@ void computeGradTensor(
         NEON_LAMBDA(const localIdx bi) {
             const auto o = bFaceCells[bi];
             const Vec3 sf = bFaceNormals[bi];
-            const Vec3 faceU = ufBound[bi];
+            const auto w = wKernel.boundaryWeight(bi, scalar(0));
+            const Vec3 faceU = w * uBound[bi];
             for (size_t row = 0; row < 3; ++row)
             {
                 for (size_t col = 0; col < 3; ++col)
@@ -313,12 +322,9 @@ void computeGradTensor(
         "computeGradTensorBoundary"
     );
 
-    // Processor faces: add the proc-face flux to the owner cell's gradient, exactly as the scalar
-    // computeGradScalar does. This block was MISSING from the tensor path, so the cell gradient at
-    // every processor-adjacent cell omitted its proc-face contribution — corrupting the Vec3
-    // corrected/limited snGrad at proc-adjacent internal faces (the uncorrected scheme, which does
-    // not use the cell gradient, was unaffected). Proc faces are compressed: use the boundary
-    // mesh's own faceNormals indexed by nBoundaryFaces+procFacei, NOT the OF-full SfAll.
+    // Proc boundary: interpolate between local owner and ghost cell value.
+    // Proc faces are compressed: use the boundary mesh's faceNormals indexed by
+    // nBoundaryFaces+procFacei, NOT the OF-full SfAll.
     const auto nProcBoundaryFaces = mesh.nProcBoundaryFaces();
     if (nProcBoundaryFaces > 0)
     {
@@ -330,7 +336,8 @@ void computeGradTensor(
                 const auto bcfacei = nBnd + procFacei;
                 const auto o = bFaceCells[bcfacei];
                 const Vec3 sf = bSf[bcfacei];
-                const Vec3 faceU = ufBound[bcfacei];
+                const auto w = wKernel.procBoundaryWeight(bcfacei, scalar(0));
+                const Vec3 faceU = w * uInt[o] + (scalar(1) - w) * uBound[bcfacei];
                 for (size_t row = 0; row < 3; ++row)
                 {
                     for (size_t col = 0; col < 3; ++col)
@@ -417,7 +424,11 @@ void GaussGreenGrad::gradTensor(
 ) const
 {
     fill(gradU.internalVector(), zero<Tensor>());
-    computeGradTensor(u, surfaceInterpolationVec_, gradU.internalVector(), operatorScaling);
+    std::visit(
+        [&](auto&& kernel)
+        { computeGradTensor(u, kernel, gradU.internalVector(), operatorScaling); },
+        surfaceInterpolationVec_.inlineWeightKernel()
+    );
     computeBoundaryGradTensor(u, gradU);
 }
 
@@ -429,8 +440,11 @@ GaussGreenGrad::gradTensor(const VolumeField<Vec3>& u, const dsl::Coeff operator
     auto calcBC = createCalculatedProcBCs<VolumeBoundary<Tensor>>(u.mesh());
     VolumeField<Tensor> gradU(u.exec(), "gradU", u.mesh(), calcBC);
     fill(gradU.internalVector(), zero<Tensor>());
-
-    computeGradTensor(u, surfaceInterpolationVec_, gradU.internalVector(), operatorScaling);
+    std::visit(
+        [&](auto&& kernel)
+        { computeGradTensor(u, kernel, gradU.internalVector(), operatorScaling); },
+        surfaceInterpolationVec_.inlineWeightKernel()
+    );
     computeBoundaryGradTensor(u, gradU);
     return gradU;
 }
diff --git a/src/linearAlgebra/ginkgo/ginkgo.cpp b/src/linearAlgebra/ginkgo/ginkgo.cpp
index 1ddc3cfb0b..7b45d7db1f 100644
--- a/src/linearAlgebra/ginkgo/ginkgo.cpp
+++ b/src/linearAlgebra/ginkgo/ginkgo.cpp
@@ -294,6 +294,9 @@ createGkoMtxImpl(std::shared_ptr<const gko::Executor> exec, const CSRMatrix<Vec3
     const auto rowsCopy = unpackRowOffs(mtx.rowOffs());
     const auto colsCopy = unpackColIdx(mtx.colIdxs(), rowsCopy, mtx.rowOffs());
     const auto valuesCopy = unpackMtxValues(mtx.values(), mtx.rowOffs(), rowsCopy);
+    // Ensure all Kokkos unpack kernels have completed before Ginkgo reads
+    // the data via gkoCopyArray, which may use a different SYCL queue.
+    Kokkos::fence();
 
     auto nrows = static_cast<gko::size_type>(3 * mtx.nRows());
     return gko::share(gko::matrix::Csr<scalar, IndexType>::create(
@@ -463,6 +466,7 @@ std::shared_ptr<const gko::matrix::Csr<scalar, IndexType>> createGkoMtxImpl(
     const auto rowsCopy = unpackRowOffs(mtx.rowOffs());
     const auto colsCopy = unpackColIdx(mtx.colIdxs(), rowsCopy, mtx.rowOffs());
     const auto valuesCopy = unpackMtxValues(mtx.values(), mtx.rowOffs(), rowsCopy);
+    Kokkos::fence();
 
     auto nrows = static_cast<gko::size_type>(computeNRows(sys));
     return gko::share(gko::matrix::Csr<scalar, IndexType>::create(
@@ -502,6 +506,9 @@ SolverStats GinkgoSolver::solve(
         const auto gkoMtx = createGkoMtx(sys.matrix());
         auto rhsCopy = unpackVecValues(sys.rhs());
         auto xCopy = unpackVecValues(x);
+        // Ensure all Kokkos unpack kernels have completed before Ginkgo reads
+        // rhsCopy/xCopy, which may use a different SYCL queue on Intel GPU.
+        Kokkos::fence();
 
         auto stats =
             solve_impl(gkoExec_, rhsCopy, xCopy, gkoMtx, factory_->generate(gkoMtx), l1Control);