diff --git a/src/ExaPF.jl b/src/ExaPF.jl
index c47bd2c8..be99319e 100644
--- a/src/ExaPF.jl
+++ b/src/ExaPF.jl
@@ -15,7 +15,7 @@ const KA = KernelAbstractions
 import Base: show, get
 
 export run_pf
-export State, Control, PolarForm
+export State, Control, AllVariables, PolarForm
 
 # Export KernelAbstractions devices
 export CPU
diff --git a/src/utils.jl b/src/utils.jl
index f777514e..517e3e87 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -140,12 +140,3 @@ function blockmul!(y::AbstractArray, A::AbstractMatrix, x::AbstractArray, alpha,
     end
 end
 
-struct BlockSparseMatrix{SMT}
-    k::Int
-    Js::Vector{SMT}
-end
-
-function BlockSparseMatrix(J::AbstractSparseMatrix, k::Int)
-    return BlockSparseMatrix(k, [copy(J) for _ in 1:k])
-end
-
diff --git a/test/Polar/TestPolarForm.jl b/test/Polar/TestPolarForm.jl
index edb56cfc..2bbacbaf 100644
--- a/test/Polar/TestPolarForm.jl
+++ b/test/Polar/TestPolarForm.jl
@@ -40,7 +40,8 @@ function runtests(datafile, device, AT)
     @testset "PolarForm API" begin
         test_polar_api(polar, device, AT)
         test_polar_stack(polar, device, AT)
-        test_polar_blk_stack(polar, device, AT)
+        test_polar_blockstack(polar, device, AT)
+        test_polar_blk_expressions(polar, device, AT)
         test_polar_constraints(polar, device, AT)
         test_polar_powerflow(polar, device, AT)
     end
@@ -49,7 +50,7 @@ function runtests(datafile, device, AT)
         test_constraints_jacobian(polar, device, AT)
         test_constraints_adjoint(polar, device, AT)
         test_full_space_jacobian(polar, device, AT)
-        test_batch_jacobian(polar, device, AT)
+        test_block_jacobian(polar, device, AT)
         test_reduced_gradient(polar, device, AT)
     end
 
diff --git a/test/Polar/api.jl b/test/Polar/api.jl
index 743d4893..e44a6f15 100644
--- a/test/Polar/api.jl
+++ b/test/Polar/api.jl
@@ -6,6 +6,7 @@ function test_polar_stack(polar, device, M)
     mapx = ExaPF.mapping(polar, State()) |> M
     mapu = ExaPF.mapping(polar, Control()) |> M
     u0 = rand(length(mapu)) |> M
+
     stack = ExaPF.NetworkStack(polar)
 
     # Test display
@@ -40,6 +41,56 @@ function test_polar_stack(polar, device, M)
     return nothing
 end
 
+function test_polar_blockstack(polar, device, M)
+    nblocks = 2
+    ngen = get(polar, PS.NumberOfGenerators())
+    nbus = get(polar, PS.NumberOfBuses())
+    nlines = get(polar, PS.NumberOfLines())
+    nu = ExaPF.number(polar, Control())
+    nx = ExaPF.number(polar, State())
+
+    mapx = ExaPF.mapping(polar, State(), nblocks) |> M
+    @test length(mapx) == nx * nblocks
+    mapu = ExaPF.mapping(polar, Control(), nblocks) |> M
+    mapxu = ExaPF.mapping(polar, AllVariables(), nblocks) |> M
+    u0 = rand(nu) |> M
+
+    x0 = rand(nx * nblocks) |> M
+
+    stack = ExaPF.BlockNetworkStack(polar, nblocks)
+
+    # Test display
+    println(devnull, stack)
+    # By defaut, values are equal to 0
+    @test isa(stack.input, M)
+    @test isa(stack.vmag, M)
+
+    # Copy state x0 inside cache
+    copyto!(stack, mapx, x0)
+    @test stack.input[mapx] == x0
+    # Copy control u0 inside cache
+    ExaPF.blockcopy!(stack, mapu, u0)
+
+    # Test that all attributes have valid length
+    @test length(stack.vang) == length(stack.vmag) == nbus * nblocks
+    @test length(stack.pgen) == ngen * nblocks
+    @test length(stack.input) == (ngen + 2 * nbus) * nblocks
+    @test length(stack.ψ) == (2 * nlines + nbus) * nblocks
+
+    # Bounds
+    b♭, b♯ = ExaPF.bounds(polar, stack)
+    @test myisless(b♭, b♯)
+
+    empty!(stack)
+    @test iszero(stack.input)
+
+    # Test constructor with scenarios
+    pload = rand(nbus, nblocks)
+    qload = rand(nbus, nblocks)
+    stack = ExaPF.BlockNetworkStack(polar, pload, qload)
+    return nothing
+end
+
 function test_polar_api(polar, device, M)
     pf = polar.network
     npv, npq, nref = length(pf.pv), length(pf.pq), length(pf.ref)
@@ -159,7 +210,7 @@ function test_polar_powerflow(polar, device, M)
     end
 end
 
-function test_polar_blk_stack(polar, device, M)
+function test_polar_blk_expressions(polar, device, M)
     nblocks = 2
     stack = ExaPF.NetworkStack(polar)
     blk_stack = ExaPF.BlockNetworkStack(polar, nblocks)
diff --git a/test/Polar/first_order.jl b/test/Polar/first_order.jl
index b9be8ffd..341ae225 100644
--- a/test/Polar/first_order.jl
+++ b/test/Polar/first_order.jl
@@ -200,7 +200,7 @@ function test_reduced_gradient(polar, device, MT)
     @test isapprox(grad_fd[:], grad_adjoint, rtol=1e-4)
 end
 
-function test_batch_jacobian(polar, device, MT)
+function test_block_jacobian(polar, device, MT)
     nblocks = 3
     mapx = ExaPF.mapping(polar, State())
 
@@ -226,5 +226,17 @@ function test_batch_jacobian(polar, device, MT)
         J_cpu = jac.J |> SparseMatrixCSC
         @test blk_J_cpu ≈ blockdiag([J_cpu for i in 1:nblocks]...)
     end
+
+    constraints = [
+        ExaPF.PowerFlowBalance(polar),
+        ExaPF.PowerGenerationBounds(polar),
+        ExaPF.VoltageMagnitudeBounds(polar),
+        ExaPF.LineFlows(polar),
+    ]
+    mycons = ExaPF.MultiExpressions(constraints) ∘ ExaPF.PolarBasis(polar)
+    for X in [State(), Control(), AllVariables()]
+        blk_jac = ExaPF.ArrowheadJacobian(polar, mycons, X, nblocks)
+        ExaPF.jacobian!(blk_jac, blk_stack)
+    end
 end
 
diff --git a/test/Polar/second_order.jl b/test/Polar/second_order.jl
index 5fd249b6..ed3fe19b 100644
--- a/test/Polar/second_order.jl
+++ b/test/Polar/second_order.jl
@@ -120,10 +120,15 @@ function test_batch_hessian(polar, device, MT)
     # Block evaluation
     blk_hess = ExaPF.ArrowheadHessian(polar, mycons, ExaPF.State(), nblocks)
     blk_H = ExaPF.hessian!(blk_hess, blk_stack, blk_y)
-
     blk_H_cpu = blk_H |> SparseMatrixCSC
     H_cpu = H |> SparseMatrixCSC
-
     @test blk_H_cpu ≈ blockdiag([H_cpu for i in 1:nblocks]...)
+
+    # Multivariables
+    for X in [State(), Control(), AllVariables()]
+        blk_hess = ExaPF.ArrowheadHessian(polar, mycons, X, nblocks)
+        blk_H = ExaPF.hessian!(blk_hess, blk_stack, blk_y)
+        @test isa(blk_H, AbstractMatrix)
+    end
 end