diff --git a/Project.toml b/Project.toml
index 4e67d7f..18375ea 100644
--- a/Project.toml
+++ b/Project.toml
@@ -51,7 +51,7 @@ BlockArrays = "1"
 CUDA = "5, 6"
 DocStringExtensions = "0.9"
 DynamicExpressions = "2.5.2"
-Exodus = "0.14"
+Exodus = "0.15"
 ForwardDiff = "1"
 GPUArrays = "11"
 Gmsh = "0.3"
diff --git a/src/Fields.jl b/src/Fields.jl
index a1328bc..102ba16 100644
--- a/src/Fields.jl
+++ b/src/Fields.jl
@@ -22,12 +22,6 @@ KA.get_backend(field::AbstractField) = KA.get_backend(field.data)
 
 abstract type AbstractContinuousField{T, D <: AbstractArray{T, 1}, NF} <: AbstractField{T, 2, D} end
 
-function Adapt.adapt_structure(to, field::AbstractContinuousField{T, D, NF}) where {T, D, NF}
-    data = adapt(to, field.data)
-    type = typeof(field).name.name
-    return eval(type){T, typeof(data), NF}(data)
-end
-
 # minimal abstractarray interface methods below
 
 function Base.axes(field::AbstractContinuousField{T, D, NF}) where {T, D, NF}
@@ -86,6 +80,46 @@ end
 
 abstract type AbstractDiscontinuousField{T, D <: AbstractArray{T, 1}} <: AbstractField{T, 1, D} end
 
+# need to implement num_fields method
+# function num_fields end
+
+# just to have something to fall back to
+Base.getindex(field::AbstractDiscontinuousField, i::Int) = field.data[i]
+Base.IndexStyle(::Type{<:AbstractDiscontinuousField}) = IndexLinear()
+Base.size(field::AbstractDiscontinuousField) = size(field.data)
+
+function Base.show(io::IO, field::AbstractDiscontinuousField)
+    println(io, "$(typeof(field)):")
+    for b in 1:num_blocks(field)
+        nf, nepe, ne = block_size(field, b)
+        println(io, "  Block $b:")
+        println(io, "    Number of fields               = $nf")
+        println(io, "    Number of entities per element = $nepe")
+        println(io, "    Number of elements             = $ne")
+    end
+end
+
+function Base.show(io::IO, ::MIME"text/plain", field::AbstractDiscontinuousField)
+    show(io, field)
+end
+
+function block_size(field::AbstractDiscontinuousField, b::Int)
+    return (num_fields(field, b), field.nepes[b], field.nelems[b])
+end
+
+function block_view(field::AbstractDiscontinuousField, b::Int)
+    nfield = num_fields(field, b)
+    nepe = field.nepes[b]
+    nelem = field.nelems[b]
+    boffset = field.offsets[b]
+    bend = boffset + nfield * nepe * nelem - 1
+    return reshape(view(field.data, boffset:bend), nfield, nepe, nelem)
+end
+
+function num_blocks(field::AbstractDiscontinuousField)
+    return length(field.nelems)
+end
+
 ######################################################################################################
 # Connectivity
 ######################################################################################################
@@ -219,6 +253,11 @@ struct H1Field{T, D, NF} <: AbstractContinuousField{T, D, NF}
     end
 end
 
+function Adapt.adapt_structure(to, field::H1Field{T, D, NF}) where {T, D, NF}
+    data = adapt(to, field.data)
+    return H1Field{T, typeof(data), NF}(data)
+end
+
 ######################################################################################################
 # HcurlField
 ######################################################################################################
@@ -229,15 +268,26 @@ Implementation of fields that live in Hdiv spaces.
 """
 struct HcurlField{T, D, NF} <: AbstractContinuousField{T, D, NF}
     data::D
+
+    function HcurlField{T, D, NF}(data::D) where {T, D, NF}
+        new{T, D, NF}(data)
+    end
+
+    function HcurlField{T, D, NF}(data::AbstractMatrix{T}) where {T, D, NF}
+        data = vec(data)
+        return HcurlField{T, D, NF}(data)
+    end
+
+    function HcurlField(data::M) where M <: AbstractMatrix
+        NF = size(data, 1)
+        data = vec(data)
+        return HcurlField{eltype(data), typeof(data), NF}(data)
+    end
 end
 
-"""
-$(TYPEDSIGNATURES)
-"""
-function HcurlField(data::M) where M <: AbstractMatrix
-    NF = size(data, 1)
-    data = vec(data)
-    return HcurlField{eltype(data), typeof(data), NF}(data)
+function Adapt.adapt_structure(to, field::HcurlField{T, D, NF}) where {T, D, NF}
+    data = adapt(to, field.data)
+    return HcurlField{T, typeof(data), NF}(data)
 end
 
 ######################################################################################################
@@ -250,15 +300,26 @@ Implementation of fields that live in Hdiv spaces.
 """
 struct HdivField{T, D, NF} <: AbstractContinuousField{T, D, NF}
     data::D
+
+    function HdivField{T, D, NF}(data::D) where {T, D, NF}
+        new{T, D, NF}(data)
+    end
+
+    function HdivField{T, D, NF}(data::AbstractMatrix{T}) where {T, D, NF}
+        data = vec(data)
+        return HdivField{T, D, NF}(data)
+    end
+
+    function HdivField(data::M) where M <: AbstractMatrix
+        NF = size(data, 1)
+        data = vec(data)
+        return HdivField{eltype(data), typeof(data), NF}(data)
+    end
 end
 
-"""
-$(TYPEDSIGNATURES)
-"""
-function HdivField(data::M) where M <: AbstractMatrix
-    NF = size(data, 1)
-    data = vec(data)
-    return HdivField{eltype(data), typeof(data), NF}(data)
+function Adapt.adapt_structure(to, field::HdivField{T, D, NF}) where {T, D, NF}
+    data = adapt(to, field.data)
+    return HdivField{T, typeof(data), NF}(data)
 end
 
 ######################################################################################################
@@ -266,89 +327,133 @@ end
 ######################################################################################################
 struct L2Field{
     T, # Let it be anything to allow for structs
-    D <: AbstractVector{T}
+    D  <: AbstractVector{T},
+    NF
 } <: AbstractDiscontinuousField{T, D}
-    data::D                    # flat storage (CPU or GPU)
-    nfields::Vector{Int}
-    nepes::Vector{Int} # num nodes, q points, etc.
+    data::D              # flat storage (CPU or GPU)
+    nepes::Vector{Int}   # num nodes, q points, etc.
     nelems::Vector{Int}
     offsets::Vector{Int}
-end
 
-function L2Field(arrs::Vector{<:AbstractArray{T, 3}}) where T
-    nfields = map(x -> size(x, 1), arrs)
-    nepes = map(x -> size(x, 2), arrs)
-    nelems = map(x -> size(x, 3), arrs)
-    offsets = Vector{eltype(nepes)}(undef, 0)
-    offset = 1
-    for b in 1:length(nepes)
-        push!(offsets, offset)
-        offset += nfields[b] * nepes[b] * nelems[b]
+    function L2Field{T, D, NF}(data, nepes, nelems, offsets) where {T, D, NF}
+        new{T, D, NF}(data, nepes, nelems, offsets)
+    end
+
+    function L2Field(arrs::Vector{<:AbstractArray{T, 3}}) where T
+        nfields = map(x -> size(x, 1), arrs)
+        @assert all(isequal(nfields[1]), nfields)
+        nfields = nfields[1]
+        nepes = map(x -> size(x, 2), arrs)
+        nelems = map(x -> size(x, 3), arrs)
+        offsets = Vector{eltype(nepes)}(undef, 0)
+        offset = 1
+        for b in 1:length(nepes)
+            push!(offsets, offset)
+            offset += nfields * nepes[b] * nelems[b]
+        end
+        data = mapreduce(vec, vcat, arrs)
+        return L2Field{T, typeof(data), nfields}(data, nepes, nelems, offsets)
     end
-    data = mapreduce(vec, vcat, arrs)
-    return L2Field(data, nfields, nepes, nelems, offsets)
-end
 
-function L2Field(::UndefInitializer, ::Type{T}, nfields::Int, qsizes::Vector{Tuple{Int, Int}}) where T
-    arrs = Array{T, 3}[]
-    for (nq, ne) in qsizes
-        push!(arrs, Array{T, 3}(undef, nfields, nq, ne))
+    function L2Field(::UndefInitializer, ::Type{T}, nfields::Int, qsizes::Vector{Tuple{Int, Int}}) where T
+        arrs = Array{T, 3}[]
+        for (nq, ne) in qsizes
+            push!(arrs, Array{T, 3}(undef, nfields, nq, ne))
+        end
+        return L2Field(arrs)
     end
-    return L2Field(arrs)
-end
 
-function L2Field(::UndefInitializer, ::Type{T}, nfields::Vector{Int}, qsizes::Vector{Tuple{Int, Int}}) where T
-    arrs = Array{T, 3}[]
-    for (nf, (nq, ne)) in zip(nfields, qsizes)
-        push!(arrs, Array{T, 3}(undef, nf, nq, ne))
+    function L2Field{T, D, NF}(::UndefInitializer, qsizes::Vector{Tuple{Int, Int}}) where {T, D, NF}
+        arrs = Array{T, 3}[]
+        for (nq, ne) in qsizes
+            push!(arrs, Array{T, 3}(undef, NF, nq, ne))
+        end
+        nepes = map(x -> size(x, 2), arrs)
+        nelems = map(x -> size(x, 3), arrs)
+        offsets = Vector{eltype(nepes)}(undef, 0)
+        offset = 1
+        for b in 1:length(nepes)
+            push!(offsets, offset)
+            offset += NF * nepes[b] * nelems[b]
+        end
+        data = mapreduce(vec, vcat, arrs)
+        return L2Field{T, typeof(data), NF}(data, nepes, nelems, offsets)
     end
-    return L2Field(arrs)
 end
 
-function Adapt.adapt_structure(to, field::L2Field{T, D}) where {T, D}
+function Adapt.adapt_structure(to, field::L2Field{T, D, NF}) where {T, D, NF}
     data = adapt(to, field.data)
-    return L2Field{T, typeof(data)}(
+    return L2Field{T, typeof(data), NF}(
         data,
-        field.nfields,
         field.nepes,
         field.nelems,
         field.offsets
     )
 end
 
-function Base.show(io::IO, field::L2Field)
-    println(io, "L2Field:")
-    for b in 1:num_blocks(field)
-        nf, nepe, ne = block_size(field, b)
-        println(io, "  Block $b:")
-        println(io, "    Number of fields               = $nf")
-        println(io, "    Number of entities per element = $nepe")
-        println(io, "    Number of elements             = $ne")
-    end
+function num_fields(::L2Field{T, D, NF}, b::Int) where {T, D, NF}
+    return NF
 end
 
-# just to have something to fall back to
-Base.getindex(field::L2Field, i::Int) = field.data[i]
-Base.IndexStyle(::Type{<:L2Field}) = IndexLinear()
-Base.size(field::L2Field) = size(field.data)
+######################################################################################################
+# StateVariableField
+######################################################################################################
+struct StateVariableField{
+    T, # Let it be anything to allow for structs
+    D <: AbstractVector{T}
+} <: AbstractDiscontinuousField{T, D}
+    data::D                    # flat storage (CPU or GPU)
+    nfields::Vector{Int}
+    nepes::Vector{Int} # num nodes, q points, etc.
+    nelems::Vector{Int}
+    offsets::Vector{Int}
 
-function Base.show(io::IO, ::MIME"text/plain", field::L2Field)
-    show(io, field)
-end
+    function StateVariableField{T, D}(data, nfields, nepes, nelems, offsets) where {T, D}
+        new{T, D}(data, nfields, nepes, nelems, offsets)
+    end
+
+    function StateVariableField(arrs::Vector{<:AbstractArray{T, 3}}) where T
+        nfields = map(x -> size(x, 1), arrs)
+        nepes = map(x -> size(x, 2), arrs)
+        nelems = map(x -> size(x, 3), arrs)
+        offsets = Vector{eltype(nepes)}(undef, 0)
+        offset = 1
+        for b in 1:length(nepes)
+            push!(offsets, offset)
+            offset += nfields[b] * nepes[b] * nelems[b]
+        end
+        data = mapreduce(vec, vcat, arrs)
+        return StateVariableField{T, typeof(data)}(data, nfields, nepes, nelems, offsets)
+    end
+
+    function StateVariableField(::UndefInitializer, ::Type{T}, nfields::Int, qsizes::Vector{Tuple{Int, Int}}) where T
+        arrs = Array{T, 3}[]
+        for (nq, ne) in qsizes
+            push!(arrs, Array{T, 3}(undef, nfields, nq, ne))
+        end
+        return StateVariableField(arrs)
+    end
 
-function block_size(field::L2Field, b::Int)
-    return (field.nfields[b], field.nepes[b], field.nelems[b])
+    function StateVariableField(::UndefInitializer, ::Type{T}, nfields::Vector{Int}, qsizes::Vector{Tuple{Int, Int}}) where T
+        arrs = Array{T, 3}[]
+        for (nf, (nq, ne)) in zip(nfields, qsizes)
+            push!(arrs, Array{T, 3}(undef, nf, nq, ne))
+        end
+        return StateVariableField(arrs)
+    end
 end
 
-function block_view(field::L2Field, b::Int)
-    nfield = field.nfields[b]
-    nepe = field.nepes[b]
-    nelem = field.nelems[b]
-    boffset = field.offsets[b]
-    bend = boffset + nfield * nepe * nelem - 1
-    return reshape(view(field.data, boffset:bend), nfield, nepe, nelem)
+function Adapt.adapt_structure(to, field::StateVariableField{T, D}) where {T, D}
+    data = adapt(to, field.data)
+    return StateVariableField{T, typeof(data)}(
+        data,
+        field.nfields,
+        field.nepes,
+        field.nelems,
+        field.offsets
+    )
 end
 
-function num_blocks(field::L2Field)
-    return length(field.nelems)
+function num_fields(field::StateVariableField, b::Int)
+    return field.nfields[b]
 end
diff --git a/src/FiniteElementContainers.jl b/src/FiniteElementContainers.jl
index 5504248..da33f4f 100644
--- a/src/FiniteElementContainers.jl
+++ b/src/FiniteElementContainers.jl
@@ -53,6 +53,7 @@ export HcurlField
 export HdivField
 export L2Field
 export Properties
+export StateVariableField
 export num_entities
 export num_fields
 
diff --git a/src/Parameters.jl b/src/Parameters.jl
index 938b91d..d53728c 100644
--- a/src/Parameters.jl
+++ b/src/Parameters.jl
@@ -17,8 +17,8 @@ function _setup_state_variables(fspace, physics)
     push!(state_old, state_old_temp)
     push!(state_new, state_new_temp)
   end
-  state_old = L2Field(state_old)
-  state_new = L2Field(state_new)
+  state_old = StateVariableField(state_old)
+  state_new = StateVariableField(state_new)
   return state_old, state_new
 end
 
@@ -63,8 +63,8 @@ struct Parameters{
   times::TimeStepper{RT}
   physics::Phys
   properties::Props
-  state_old::L2Field{RT, RV}
-  state_new::L2Field{RT, RV}
+  state_old::StateVariableField{RT, RV}
+  state_new::StateVariableField{RT, RV}
   coords::Coords
   field::Field
   field_old::Field
@@ -216,8 +216,8 @@ struct TypeStableParameters{
   times::TimeStepper{RT}
   physics::Phys
   properties::Props
-  state_old::L2Field{RT, RV}
-  state_new::L2Field{RT, RV}
+  state_old::StateVariableField{RT, RV}
+  state_new::StateVariableField{RT, RV}
   coords::Coords
   field::Field
   field_old::Field
diff --git a/src/assemblers/SparseMatrixAssembler.jl b/src/assemblers/SparseMatrixAssembler.jl
index 3c91336..74035b2 100644
--- a/src/assemblers/SparseMatrixAssembler.jl
+++ b/src/assemblers/SparseMatrixAssembler.jl
@@ -21,7 +21,7 @@ struct SparseMatrixAssembler{
   mass_storage::RV
   residual_storage::FieldStorage
   residual_unknowns::RV
-  scalar_quadrature_storage::L2Field{Float64, RV}
+  scalar_quadrature_storage::L2Field{Float64, RV, 1}
   stiffness_storage::RV
   stiffness_action_storage::FieldStorage
   stiffness_action_unknowns::RV
@@ -103,7 +103,8 @@ function SparseMatrixAssembler{SparseMatrixType, UseInPlaceMethods, UseSparseVec
 
   # setup quadrature scalar storage
   fspace = function_space(dof)
-  scalar_quadrature_storage = L2Field(undef, Float64, 1, block_quadrature_sizes(fspace))
+  # scalar_quadrature_storage = L2Field(undef, Float64, 1, block_quadrature_sizes(fspace))
+  scalar_quadrature_storage = L2Field{Float64, Vector{Float64}, 1}(undef, block_quadrature_sizes(fspace))
   fill!(scalar_quadrature_storage, 0.0)
 
   return SparseMatrixAssembler{
diff --git a/test/TestAssemblers.jl b/test/TestAssemblers.jl
index 90c7e0d..030a977 100644
--- a/test/TestAssemblers.jl
+++ b/test/TestAssemblers.jl
@@ -180,6 +180,7 @@ end
 end
 
 @testitem "Assembler - assembler_consistency_mechanics" setup=[AssemblerHelperMechanics] begin
+  # using KernelAbstractions
   include("TestUtils.jl")
   backends = _get_backends()
   for dev in backends
@@ -217,6 +218,7 @@ end
     
         # test vector consistency
         assemble_vector!(asm_1, residual, U_1, p_1)
+      
         assemble_vector!(asm_2, residual!, U_2, p_2)
         R_1 = residual(asm_1) |> copy
         R_2 = residual(asm_2) |> copy
diff --git a/test/TestFields.jl b/test/TestFields.jl
index 871cabe..68798ce 100644
--- a/test/TestFields.jl
+++ b/test/TestFields.jl
@@ -163,11 +163,11 @@ end
 
 @testitem "Fields - test_l2_field" begin
   a1 = rand(2, 3, 40)
-  a2 = rand(3, 4, 10)
+  a2 = rand(2, 4, 10)
   field = L2Field([a1, a2])
   @show field
   @test size(FiniteElementContainers.block_view(field, 1)) == (2, 3, 40)
-  @test size(FiniteElementContainers.block_view(field, 2)) == (3, 4, 10)
+  @test size(FiniteElementContainers.block_view(field, 2)) == (2, 4, 10)
 
   bview = FiniteElementContainers.block_view(field, 1)
   for k in axes(bview, 3)
@@ -187,3 +187,30 @@ end
     end
   end
 end
+
+@testitem "Fields - test_state_variable_field" begin
+  a1 = rand(2, 3, 40)
+  a2 = rand(3, 4, 10)
+  field = StateVariableField([a1, a2])
+  @show field
+  @test size(FiniteElementContainers.block_view(field, 1)) == (2, 3, 40)
+  @test size(FiniteElementContainers.block_view(field, 2)) == (3, 4, 10)
+
+  bview = FiniteElementContainers.block_view(field, 1)
+  for k in axes(bview, 3)
+    for j in axes(bview, 2)
+      for i in axes(bview, 1)
+        @test bview[i, j, k] ≈ a1[i, j, k]
+      end
+    end
+  end
+
+  bview = FiniteElementContainers.block_view(field, 2)
+  for k in axes(bview, 3)
+    for j in axes(bview, 2)
+      for i in axes(bview, 1)
+        @test bview[i, j, k] ≈ a2[i, j, k]
+      end
+    end
+  end
+end
\ No newline at end of file