diff --git a/src/ScalarEq/callbacks.jl b/src/ScalarEq/callbacks.jl
index 347e8c73ddd87d19cd59e68fd3290b07c7533c6e..439316b25ff3bb4e9c747c779054c36251d9aecd 100644
--- a/src/ScalarEq/callbacks.jl
+++ b/src/ScalarEq/callbacks.jl
@@ -1,11 +1,11 @@
-function make_callback(env, P::Project, bdryconds::BoundaryConditions)
-  cbfn_equation(u, t) = callback_equation(u, t, env, P, isperiodic(bdryconds), P.equation, env.mesh)
+function make_callback(env, P::Project, isperiodic::Bool)
+  cbfn_equation(u, t) = callback_equation(u, t, env, P, isperiodic, P.equation, env.mesh)
   cb_equation         = FunctionCallback(cbfn_equation,
                                          CallbackTiming(every_iteration=1,every_dt=0))
-  cbfn_tci(u, t)      = callback_tci(u, t, env, P, isperiodic(bdryconds), P.tci, env.mesh)
+  cbfn_tci(u, t)      = callback_tci(u, t, env, P, isperiodic, P.tci, env.mesh)
   cb_tci              = FunctionCallback(cbfn_tci,
                                          CallbackTiming(every_iteration=1, every_dt=0))
-  cbfn_hrsc(u, t)     = callback_hrsc(u, t, env, P, isperiodic(bdryconds), P.hrsc, env.mesh)
+  cbfn_hrsc(u, t)     = callback_hrsc(u, t, env, P, isperiodic, P.hrsc, env.mesh)
   cb_hrsc             = FunctionCallback(cbfn_hrsc,
                                          CallbackTiming(every_iteration=1,every_dt=0))
   callbackset = CallbackSet(cb_equation, cb_tci, cb_hrsc)
@@ -206,8 +206,8 @@ end
 function callback_tci(state_u, state_t, env, P, isperiodic,
     tci::Union{TCI.Diffusion, TCI.Minmod, TCI.Threshold,
                TCI.ModalDecayAverage, TCI.ModalDecayHighest}, mesh)
-  @unpack cache, mesh            = env
-  @unpack u                      = get_dynamic_variables(cache)
-  @unpack flag                   = get_cell_variables(cache)
+  @unpack cache, mesh = env
+  @unpack u           = get_dynamic_variables(cache)
+  @unpack flag        = get_cell_variables(cache)
   TCI.compute_indicator!(flag, u, tci)
 end
diff --git a/src/ScalarEq/equation.jl b/src/ScalarEq/equation.jl
index 652c4a2f9990ac174fa266790f3e393eb9504b9b..155e43abe6fa44be6cd745f0c8313c971371538b 100644
--- a/src/ScalarEq/equation.jl
+++ b/src/ScalarEq/equation.jl
@@ -84,13 +84,13 @@ end
 
 
 @with_signature function flux(eq::Advection)
-  @accepts State(u)
+  @accepts u
   flx_u = eq.v * u
   @returns flx_u
 end
 
 @with_signature function speed(eq::Advection)
-  @accepts State(u)
+  @accepts u
   v = eq.v
   @returns v
 end
@@ -113,13 +113,13 @@ end
 
 
 @with_signature function flux(eq::Advection2d)
-  @accepts State(u)
+  @accepts u
   flx_x, flx_y = eq.v_x * u, eq.v_y * u
   @returns flx_x, flx_y
 end
 
 @with_signature function speed(eq::Advection2d)
-  @accepts State(u)
+  @accepts u
   v = sqrt(eq.v_x^2 + eq.v_y^2)
   @returns v
 end
@@ -145,13 +145,13 @@ end
 
 
 @with_signature function flux(eq::Burgers)
-  @accepts State(u)
+  @accepts u
   flx_u = u^2 / 2
   @returns flx_u
 end
 
 @with_signature function speed(eq::Burgers)
-  @accepts State(u)
+  @accepts u
   v = u
   @returns v
 end
@@ -175,13 +175,13 @@ end
 
 
 @with_signature function flux(eq::Burgers2d)
-  @accepts State(u)
+  @accepts u
   flx_x, flx_y = u^2/2, u^2/2
   @returns flx_x, flx_y
 end
 
 @with_signature function speed(eq::Burgers2d)
-  @accepts State(u)
+  @accepts u
   v = sqrt(2) * u
   @returns v
 end
diff --git a/src/ScalarEq/rhs.jl b/src/ScalarEq/rhs.jl
index 63b0f03be5ace7330fe17feebe8ef4850dc27fde..0649c2d4667b9cde360ed05bfe2e95c5c38218d4 100644
--- a/src/ScalarEq/rhs.jl
+++ b/src/ScalarEq/rhs.jl
@@ -2,35 +2,35 @@
 #                              Timestep                               #
 #######################################################################
 
-function timestep(env, P::Project, ::Mesh1d, hrsc::Maybe{HRSC.AbstractHRSC})
-  @unpack cache, mesh = env
-  @unpack equation = P
-  @unpack v = get_static_variables(cache)
 
-  broadcast_volume!(speed, equation, cache)
+function maxspeed(mesh, equation, cache)
+  @unpack v = get_static_variables(cache)
+  broadcast_volume_2!(speed, equation, mesh)
   vmax = dg1d.absolute_maximum(v)
   vmax_limit = 1e4
   if vmax > vmax_limit
     @warn "Limiting timestep due to maximum speed exceeding $vmax_limit"
     vmax = vmax_limit
   end
+  return vmax
+end
 
-  @unpack element = mesh
-  @unpack N, z = element
-  dz = z[2]-z[1]
-  dx, = minimum(dg1d.widths(mesh.boxes[1]))
-  dl = dx * dz
 
-  dt = dl / (N * vmax)
+function timestep(env, P::Project, mesh::Mesh1d, hrsc::Maybe{HRSC.AbstractHRSC})
+  vmax = maxspeed(mesh, P.equation, mesh.cache)
+  dl = min_grid_spacing(mesh)
+  dt = dl / (vmax * dtfactor(mesh))
   return dt
 end
+dtfactor(mesh::Mesh1d{FVElement}) = 2
+dtfactor(mesh::Mesh1d{SpectralElement}) = mesh.element.N
 
 
-function timestep(env, P::Project, ::Mesh2d, hrsc::Maybe{HRSC.AbstractHRSC})
+function timestep(env, P::Project, ::Mesh2d{SpectralElement}, hrsc::Maybe{HRSC.AbstractHRSC})
   @unpack cache, mesh = env
   @unpack equation = P
 
-  broadcast_volume!(speed, equation, cache)
+  broadcast_volume_2!(speed, equation, mesh)
 
   @unpack v = get_static_variables(cache)
   vmax = dg1d.absolute_maximum(v)
@@ -66,11 +66,30 @@ function rhs!(env, P::Project, hrsc, bdryconds...)
 end
 
 
+function rhs!(env, mesh::Mesh1d{FVElement}, P::Project, hrsc::Nothing,
+    bdryconds, ldg_bdryconds::Nothing, av_bdryconds::Nothing)
+
+  @unpack cache, mesh = env
+  @unpack equation = P
+
+  broadcast_volume_2!(flux, equation, mesh)
+
+  @unpack flx_u  = get_static_variables(cache)
+  @unpack u      = get_dynamic_variables(cache)
+  @unpack rhs_u  = get_rhs_variables(cache)
+
+  dt = timestep(env, P, mesh, hrsc)
+  fv_update_step!(rhs_u, u, flx_u, dt, mesh)
+
+  return
+end
+
+
 function rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing,
     bdryconds, ldg_bdryconds::Nothing, av_bdryconds::Nothing)
 
   @unpack cache, mesh = env
-  @unpack equation, rsolver = P
+  @unpack equation = P
 
   broadcast_volume_2!(flux, equation, mesh)
   broadcast_faces_2!(llf_1d, equation, mesh)
@@ -90,13 +109,11 @@ function rhs!(env, mesh::Mesh2d, P::Project, hrsc::Nothing,
     bdryconds, ldg_bdryconds::Nothing, av_bdryconds::Nothing)
 
   @unpack cache = env
-  @unpack equation, rsolver = P
+  @unpack equation = P
 
   broadcast_volume_2!(flux, equation, mesh)
   broadcast_faces_2!(llf_2d, equation, mesh)
   dg1d.broadcast_bdry_2!(bdryllf_2d, equation, P.bdrycond, mesh)
-  # broadcast_boundaryconditions!(central_flux, bdryconds_x, cache, mesh, 0.0)
-  # broadcast_boundaryconditions!(central_flux, bdryconds_y, cache, mesh, 0.0)
 
   @unpack flx_x, flx_y = get_static_variables(cache)
   @unpack u            = get_dynamic_variables(cache)
@@ -123,9 +140,9 @@ function rhs!(env, mesh::Mesh1d, P::Project, hrsc::Union{BernsteinReconstruction
 
   HRSC.reconstruct!(u, flag, hrsc, isperiodic=isperiodic(bdryconds))
 
-  broadcast_volume!(flux, equation, cache)
-  broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
-  broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
+  broadcast_volume_2!(flux, equation, mesh)
+  broadcast_faces_2!(llf, equation, mesh)
+  broadcast_bdry_2!(bdry_llf, bdryconds, mesh)
 
   compute_rhs_weak_form!(rhs_u, flx_u, lhs_numflx_u, rhs_numflx_u, mesh)
 
@@ -181,17 +198,14 @@ function rhs!(env, mesh::Mesh2d, P::Project, hrsc::HRSC.AbstractArtificialViscos
   broadcast_volume_2!(ldg_flux_2d, equation, mesh)
   broadcast_faces_2!(ldg_nflux_qx, equation, mesh)
   broadcast_faces_2!(ldg_nflux_qy, equation, mesh)
-  # broadcast_boundaryconditions!(central_flux, ldg_bdryconds, cache, mesh, 0.0)
   compute_rhs_weak_form!(qx, flx_qx_x, flx_qx_y, nflx_qx, mesh)
   compute_rhs_weak_form!(qy, flx_qy_x, flx_qy_y, nflx_qy, mesh)
 
   # ## compute rhs of regularized equation: ∂t u + ∂x f + ∂y f + ∂x μ qx + ∂y μ qy = 0
   broadcast_volume_2!(av_flux_2d, equation, mesh)
   broadcast_faces_2!(av_nflux_2d, equation, mesh)
-  # broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, 0.0)
   broadcast_volume_2!(flux, equation, mesh)
   broadcast_faces_2!(llf_2d, equation, mesh)
-  # broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
 
   # add up fluxs and numerical fluxes
   @. flx_x  += flx_g_x
diff --git a/src/ScalarEq/setup.jl b/src/ScalarEq/setup.jl
index c1edec36d868229a7008eafbabeec4561e213442..26230f212e70981c189dc04523df6f8afe42c132 100644
--- a/src/ScalarEq/setup.jl
+++ b/src/ScalarEq/setup.jl
@@ -45,7 +45,7 @@ function Project(env::Environment, mesh::Mesh2d, prms)
   equation = make_Equation(prms["ScalarEq"], dimension(mesh))
   tci  = TCI.make_TCI(mesh, prms["TCI"])
   hrsc = HRSC.make_HRSC(mesh, prms["HRSC"])
-  rsolver = ApproxRiemannSolver2d(flux, speed, equation)
+  rsolver = nothing
   ldg_rsolver, av_rsolver = nothing, nothing
 
   bdrycond = dg1d.DirichletBC2()
@@ -58,19 +58,18 @@ function Project(env::Environment, mesh::Mesh2d, prms)
   _register_variables!(mesh, tci)
   _register_variables!(mesh, hrsc)
   _register_variables!(cache, bdrycond)
-  register_variables!(cache, rsolver)
   display(cache)
 
   # setup initial data
   initialdata!(env, P, prms["ScalarEq"])
 
-  # setup boundary conditions
-  bdryconds_x = BoundaryConditions(PeriodicBC(), PeriodicBC(), rsolver)
-  bdryconds_y = BoundaryConditions(PeriodicBC(), PeriodicBC(), rsolver)
+  # # setup boundary conditions
+  bdryconds_x = nothing
+  bdryconds_y = nothing
   ldg_bdryconds, av_bdryconds = nothing, nothing
 
   # setup callbacks
-  projectcb = make_callback(env, P, bdryconds_x)
+  projectcb = make_callback(env, P, isperiodic(mesh))
   append!(env.callbacks, CallbackSet(projectcb.callbacks...))
 
   return P, ((bdryconds_x, bdryconds_y), ldg_bdryconds, av_bdryconds)
@@ -159,6 +158,7 @@ function _register_variables!(cache, eq::Union{Advection2d,Burgers2d})
     rhs_variablenames     = (:rhs_u,),
     static_variablenames  = (:flx, :flx_x, :flx_y, :v, :v_x, :v_y),
     cell_variablenames    = (:max_v,),
+    bdry_variablenames    = (:nflx_u,),
     global_variablenames  = (:t, :tm1))
 end
 
diff --git a/src/TCI/threshold.jl b/src/TCI/threshold.jl
index eb9cd04733b837d0347b17864939f97d1bf2873d..52cfd3f9b30a5cf28d58ccf1329d9b43c49501d5 100644
--- a/src/TCI/threshold.jl
+++ b/src/TCI/threshold.jl
@@ -8,7 +8,7 @@ end
 function compute_indicator!(
     flag, # mutated outputs
     u,    # inputs
-    tci::Threshold{Mesh1d})
+    tci::Threshold{<:Mesh1d})
 
   @unpack mesh, threshold = tci
   L       = layout(mesh)
@@ -34,7 +34,7 @@ end
 function compute_indicator!(
     flag, # mutated outputs
     u,    # inputs
-    tci::Threshold{Mesh2d})
+    tci::Threshold{<:Mesh2d})
 
   @unpack mesh, threshold = tci
   @unpack Npts = mesh.element
diff --git a/src/dg1d.jl b/src/dg1d.jl
index 16fe42fccf1a2cdf6f2913333f0cd9c31fc0cf96..451e650da949cc287a798c7cf94a558e90ec49cf 100644
--- a/src/dg1d.jl
+++ b/src/dg1d.jl
@@ -72,6 +72,9 @@ include("lgl.jl")
 export SpectralElement
 include("spectralelement.jl")
 
+export FVElement
+include("fvelement.jl")
+
 include("tensorbasis.jl")
 include("box.jl")
 include("tree.jl")
@@ -94,7 +97,8 @@ include("cache.jl")
 
 export Mesh, Mesh1d, Mesh2d, MeshInterpolator, layout, n_points, n_cells, dimension,
        grid_average, cellwise_average, cellwise_inner_product, broken_inner_product,
-       locate_point, find_cell_index, eachcell, cellindices
+       locate_point, find_cell_index, eachcell, cellindices,
+       min_grid_spacing
 include("mesh.jl")
 
 export AbstractEquation
@@ -126,6 +130,9 @@ include("with_signature.jl")
 export compute_rhs_weak_form!
 include("dg_rhs.jl")
 
+export fv_update_step!
+include("fv_rhs.jl")
+
 # TODO Move this to top; also requires gather all type defintions
 # like SpectralElement, AbstractMesh, Tree etc. into types.jl
 export Maybe, Environment
diff --git a/src/dg_rhs.jl b/src/dg_rhs.jl
index d668df23a6f1db8b22157f57521a5293aa75aeaf..02e50bbee11c108ae1a0611fa8a3194822272656 100644
--- a/src/dg_rhs.jl
+++ b/src/dg_rhs.jl
@@ -1,8 +1,11 @@
 """
-    compute_rhs_weak_form!(rhs, f, nf_lhs, nf_rhs, mesh)
-    compute_rhs_weak_form!(rhs, f, s, nf_lhs, nf_rhs, mesh)
+    compute_rhs_weak_form!(rhs, f, nf, mesh::Mesh1d{SpectralElement})
+    compute_rhs_weak_form!(rhs, f, s, nf, mesh::Mesh1d{SpectralElement})
+
+Compute the `rhs` of the weak DG formulation using
+the (bulk) flux `f`, the numerical flux `nf` and sources `s` for a `mesh`.
 """
-function compute_rhs_weak_form!(rhs, f, nf, mesh::Mesh1d)
+function compute_rhs_weak_form!(rhs, f, nf, mesh::Mesh1d{SpectralElement})
   @unpack invjac, element = mesh
   @unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
   @unpack K = mesh
@@ -18,13 +21,21 @@ function compute_rhs_weak_form!(rhs, f, nf, mesh::Mesh1d)
   end
   return
 end
-function compute_rhs_weak_form!(rhs, f, s, nf, mesh::Mesh1d)
+function compute_rhs_weak_form!(rhs, f, s, nf, mesh::Mesh1d{SpectralElement})
   compute_rhs_weak_form!(rhs, f, nf, mesh)
   rhs .+= s
   return
 end
 
-function compute_rhs_weak_form!(rhs, fx, fy, nf, mesh::Mesh2d)
+
+"""
+    compute_rhs_weak_form!(rhs, fx, fy, nf, mesh::Mesh2d{SpectralElement})
+    compute_rhs_weak_form!(rhs, fx, fy, s, nf, mesh::Mesh2d{SpectralElement})
+
+Compute the `rhs` of the weak DG formulation using
+the (bulk) flux `fx, fy`, the numerical flux `nf` and sources `s` for a `mesh`.
+"""
+function compute_rhs_weak_form!(rhs, fx, fy, nf, mesh::Mesh2d{SpectralElement})
   @unpack element = mesh
   @unpack dxdX, dydY = get_static_variables(mesh.cache)
   @unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
@@ -93,7 +104,15 @@ function compute_rhs_weak_form!(rhs, fx, fy, nf, mesh::Mesh2d)
   end
   return
 end
-function compute_rhs_weak_form!(rhs, fx, fy::Real, nf, mesh::Mesh2d)
+
+"""
+    compute_rhs_weak_form!(rhs, fx, _::Real, nf, mesh::Mesh2d{SpectralElement})
+    compute_rhs_weak_form!(rhs, _::Real, fy, nf, mesh::Mesh2d{SpectralElement})
+
+Specialized version of `compute_rhs_weak_form!` where only the fluxs
+`fx` or `fy` are applied, respectively.
+"""
+function compute_rhs_weak_form!(rhs, fx, _::Real, nf, mesh::Mesh2d{SpectralElement})
   @unpack element = mesh
   @unpack dxdX, dydY = get_static_variables(mesh.cache)
   @unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
@@ -149,7 +168,7 @@ function compute_rhs_weak_form!(rhs, fx, fy::Real, nf, mesh::Mesh2d)
   end
   return
 end
-function compute_rhs_weak_form!(rhs, fx::Real, fy, nf, mesh::Mesh2d)
+function compute_rhs_weak_form!(rhs, _::Real, fy, nf, mesh::Mesh2d{SpectralElement})
   @unpack element = mesh
   @unpack dxdX, dydY = get_static_variables(mesh.cache)
   @unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
@@ -205,7 +224,7 @@ function compute_rhs_weak_form!(rhs, fx::Real, fy, nf, mesh::Mesh2d)
   end
   return
 end
-function compute_rhs_weak_form!(rhs, fx, fy, s, nf, mesh::Mesh2d)
+function compute_rhs_weak_form!(rhs, fx, fy, s, nf, mesh::Mesh2d{SpectralElement})
   compute_rhs_weak_form!(rhs, fx, fy, nf, mesh)
   rhs .+= s
   return
diff --git a/src/evolve.jl b/src/evolve.jl
index 9e7520d8558addbd61eb37eb9ca6aa5dd18eb149..900a0f8f1806211ce398556fafeef1ccad9f2128 100644
--- a/src/evolve.jl
+++ b/src/evolve.jl
@@ -27,7 +27,8 @@ Optional values
 """
 function Evolution(rhs_fn!, u0::Vector{Float64}, timestep, tend, alg::TimeStepAlgorithm;
     callback_fullstep::Union{Function,CallbackSet,Nothing}=nothing,
-    callback_substep::Union{Function,CallbackSet,Nothing}=nothing)
+    callback_substep::Union{Function,CallbackSet,Nothing}=nothing,
+    mesh::Union{Mesh,Nothing}=nothing)
 
   # enforce interface
   timestep_fn = if timestep isa Number
@@ -58,7 +59,11 @@ function Evolution(rhs_fn!, u0::Vector{Float64}, timestep, tend, alg::TimeStepAl
   end
 
   stages = make_RK_stages(size(u0), alg.nstages)
-  stepper!(up1, u0, t, dt) = step!(up1, rhs_fn!, u0, t, dt, stages, alg, cb_substep)
+  stepper! = if mesh isa Mesh1d{FVElement}
+    (up1, u0, t, dt) -> rhs_fn!(up1, u0, t)
+  else
+    (up1, u0, t, dt) -> step!(up1, rhs_fn!, u0, t, dt, stages, alg, cb_substep)
+  end
   up1 = deepcopy(u0) # we need to initialize memory anyways
                      # and the stepper will swap up1, u before he steps
 
@@ -106,7 +111,7 @@ function step!(evo::Evolution)
   evo.it += 1
 
   if cb_fullstep(up1, evo.t, evo.it) == false
-    @info """Callback after step failed at t = $(t)!"""
+    @info """Callback failed after step at t = $(t)!"""
     return :callback_failed
   end
 
diff --git a/src/fv_rhs.jl b/src/fv_rhs.jl
new file mode 100644
index 0000000000000000000000000000000000000000..18c6dcfc7b5ba1d0258b26e70136e77ee8278f7b
--- /dev/null
+++ b/src/fv_rhs.jl
@@ -0,0 +1,42 @@
+"""
+    fv_update_step!(up1, u, f, mesh::Mesh1d{FVElement})
+    fv_update_step!(up1, u, f, s, mesh::Mesh1d{FVElement})
+
+Update the state `up1` of the FV formulation using
+the (bulk) flux `f` and sources `s` for a `mesh`.
+"""
+function fv_update_step!(up1, u, f, dt, mesh::Mesh1d{FVElement})
+  @unpack invjac = mesh
+  @unpack K = mesh
+  dl = widths(mesh)[1] / K
+  dtdl = dt/dl
+  @turbo for j = 2:K-1
+    up1[j] = (u[j+1] + 2*u[j] + u[j-1])/4 - dtdl/2 * (f[j+1] - f[j-1])
+  end
+  if mesh.tree.periodic[1]
+    up1[1]   = (u[2] + 2*u[1]   + u[end])/4   - dtdl/2 * (f[2] - f[end])
+    up1[end] = (u[1] + 2*u[end] + u[end-1])/4 - dtdl/2 * (f[1] - f[end-1])
+  else
+    TODO()
+  end
+  return
+end
+
+
+function fv_update_step!(up1, u, f, s, dt, mesh::Mesh1d{FVElement})
+  TODO()
+  @unpack invjac = mesh
+  @unpack K = mesh
+  dl = widths(mesh)[1] / K
+  dtdl = dt/dl
+  @turbo for j = 2:K-1
+    up1[j] = (u[j+1] + 2*u[j] + u[j-1])/4 - dtdl/2 * (f[j+1] - f[j-1]) + s[j]
+  end
+  if mesh.tree.periodic[1]
+    up1[1]   = (u[2] + 2*u[1]   + u[end])/4   - dtdl/2 * (f[2] - f[end]) + s[1]
+    up1[end] = (u[1] + 2*u[end] + u[end-1])/4 - dtdl/2 * (f[1] - f[end-1]) + s[end]
+  else
+    TODO()
+  end
+  return
+end
diff --git a/src/fvelement.jl b/src/fvelement.jl
new file mode 100644
index 0000000000000000000000000000000000000000..f6b3fb50ec744464eaa8a922172356c96006bd7e
--- /dev/null
+++ b/src/fvelement.jl
@@ -0,0 +1,36 @@
+struct FVElement
+  N::Int64                # polynomial order
+  Npts::Int64             # number of quadrature points = N + 1
+  z::Vector{Float64}      # quadrature points
+
+  function FVElement()
+    N = 0
+    Npts = N + 1
+    z = [0.0]
+    return new(N, Npts, z)
+  end
+end
+
+
+@inline function inner_product(u1, u2, el::FVElement)
+  TODO()
+  @unpack quadr = el
+  @toggled_assert quadr in (:LGL, :LGL_lumped, :GLGL)
+  if quadr === :LGL || quadr === :LGL_lumped
+    return LGL.integrate(el.w, u1, u2)
+  else # quadr === :GLGL
+    return GLGL.integrate(el.w, el.v, el.D, u1, u2)
+  end
+end
+
+
+@inline function integrate(u, el::FVElement)
+  TODO()
+  @unpack quadr = el
+  @toggled_assert quadr in (:LGL, :LGL_lumped, :GLGL)
+  if quadr === :LGL || quadr === :LGL_lumped
+    return LGL.integrate(el.w, u)
+  else # quadr === :GLGL
+    return GLGL.integrate(el.w, el.v, el.D, u)
+  end
+end
diff --git a/src/main.jl b/src/main.jl
index 8ef125347cb8c03277c1528fb715c211d9e33109..9bec199229a10f1271fa8d587b6ad4fdc349385f 100644
--- a/src/main.jl
+++ b/src/main.jl
@@ -174,7 +174,8 @@ function setup(project_name, prms, outputdir)
 
   mprms = prms["Mesh"]
   mesh = Mesh(; N=mprms["n"], K=mprms["k"], range=mprms["range"],
-              basis=mprms["basis"], periodic=mprms["periodic"])
+              basis=mprms["basis"], periodic=mprms["periodic"],
+              scheme=mprms["scheme"])
   callbacks = CallbackSet()
   callbacks_substeps = CallbackSet()
   env = Environment(mesh, mesh.cache, callbacks, callbacks_substeps)
@@ -285,7 +286,7 @@ function evolve!(env, rhs_fn, timestep_fn, prms)
   # TODO This should be hidden behind an interface
   u0 = reduce(vcat, fields(cache.dynamic_variables))
   alg = algorithm(Symbol(method))
-  evolution = Evolution(rhs_fn, u0, timestep_fn, tend, alg;
+  evolution = Evolution(rhs_fn, u0, timestep_fn, tend, alg; mesh=mesh,
                         callback_fullstep=callbacks)#, callback_substep=callbacks_substep)
 
   # run once on initial data
diff --git a/src/mesh.jl b/src/mesh.jl
index 030a7e4fd58b6073c21ed29c714f32844432f4b8..65326e6c56b49bae624b5f35d2fbe8e8a1358a8a 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -6,11 +6,11 @@
 abstract type AbstractMesh end
 
 
-struct Mesh{N_Dim,N_Sides} <: AbstractMesh
+struct Mesh{Element,N_Dim,N_Sides} <: AbstractMesh
   tree::Tree{N_Dim,N_Sides} # abstract representation of mesh
   boxes::Vector{Box{N_Dim}} # physical extends of each cell
   extends::NTuple{N_Dim,Tuple{Float64,Float64}} # total extend of mesh
-  element::SpectralElement
+  element::Element
   cache::Cache
   bulkfaceindices::Vector{Int64}
   bulkbdryindices::Vector{Int64}
@@ -20,8 +20,8 @@ struct Mesh{N_Dim,N_Sides} <: AbstractMesh
 end
 
 
-const Mesh1d = Mesh{1,2}
-const Mesh2d = Mesh{2,4}
+const Mesh1d{Element} = Mesh{Element,1,2}
+const Mesh2d{Element} = Mesh{Element,2,4}
 
 
 #######################################################################
@@ -29,7 +29,7 @@ const Mesh2d = Mesh{2,4}
 #######################################################################
 
 
-function Mesh1d(; N=5, K=10, range=[-1.0,1.0], basis="lgl_lumped", periodic=(true,))
+function Mesh1d(; N=5, K=10, range=[-1.0,1.0], basis="lgl_lumped", periodic=(true,), scheme="DG")
 
   @toggled_assert N > 0
   @toggled_assert K > 0
@@ -37,12 +37,17 @@ function Mesh1d(; N=5, K=10, range=[-1.0,1.0], basis="lgl_lumped", periodic=(tru
   @toggled_assert range[1] < range[2]
   @toggled_assert basis in ["lgl", "glgl", "lgl_lumped"]
 
-  # TODO Remove this after updating parameter names in SpectralElement
-  quadrature_method = Dict("lgl" => :LGL, "glgl" => :GLGL, "lgl_lumped" => :LGL_lumped)[basis]
-
   tree = Tree1d(K,periodic=tuple(periodic...))
   boxes = place_boxes(tree, Float64.(range))
-  element = SpectralElement(N, Symbol(quadrature_method))
+  element = if scheme == "DG"
+    # TODO Remove this after updating parameter names in SpectralElement
+    quadrature_method = Dict("lgl" => :LGL, "glgl" => :GLGL, "lgl_lumped" => :LGL_lumped)[basis]
+    SpectralElement(N, Symbol(quadrature_method))
+  elseif scheme == "FV"
+    FVElement()
+  else
+    error("unknown approximation scheme $scheme")
+  end
   @unpack Npts, z = element
   offsets = [ (i-1)*Npts for i = 1:length(tree) ]
   extends = (Tuple(range),)
@@ -80,13 +85,14 @@ function Mesh1d(; N=5, K=10, range=[-1.0,1.0], basis="lgl_lumped", periodic=(tru
     Int64[1, Npts*K], Int64[1, 2*K]
   end
 
-  return Mesh1d(tree, boxes, extends, element, cache,
+  return Mesh1d{typeof(element)}(tree, boxes, extends, element, cache,
                 bulkfaceindices, bulkbdryindices, faceindices, bdryindices, offsets)
 end
 
 
 function Mesh2d(; N=[5,5], K=[4,4], range=[-1.0,1.0, -1.0,1.0],
-                  basis="lgl_lumped", periodic=(true,true))
+                  basis="lgl_lumped", periodic=(true,true),
+                  scheme="DG")
 
   @toggled_assert length(N) == 2
   @toggled_assert length(K) == 2
@@ -103,13 +109,20 @@ function Mesh2d(; N=[5,5], K=[4,4], range=[-1.0,1.0, -1.0,1.0],
   @toggled_assert yrange[1] < yrange[2]
   @toggled_assert basis in ["lgl", "glgl", "lgl_lumped"]
 
-  # TODO Remove this after updating parameter names in SpectralElement
-  quadrature_method = Dict("lgl" => :LGL, "glgl" => :GLGL, "lgl_lumped" => :LGL_lumped)[basis]
-
   tree = Tree2d(Kx,Ky,periodic=tuple(periodic...))
   boxes = place_boxes(tree, Float64.(xrange), Float64.(yrange))
-  element_x = SpectralElement(Nx, Symbol(quadrature_method))
-  element_y = SpectralElement(Ny, Symbol(quadrature_method))
+  if scheme == "DG"
+    # TODO Remove this after updating parameter names in SpectralElement
+    quadrature_method = Dict("lgl" => :LGL, "glgl" => :GLGL, "lgl_lumped" => :LGL_lumped)[basis]
+    element_x = SpectralElement(Nx, Symbol(quadrature_method))
+    element_y = SpectralElement(Ny, Symbol(quadrature_method))
+  elseif scheme == "FV"
+    TODO()
+    element_x = FVElement()
+    element_y = FVElement()
+  else
+    error("unknown approximation scheme $scheme")
+  end
   Nptsx, Nptsy = element_x.Npts, element_y.Npts
   offsets = [ (i-1)*Nptsx*Nptsy for i = 1:length(tree) ]
   extends = (Tuple(xrange),Tuple(yrange))
@@ -199,16 +212,16 @@ function Mesh2d(; N=[5,5], K=[4,4], range=[-1.0,1.0, -1.0,1.0],
     bo += Nptsx
   end
 
-  return Mesh2d(tree, boxes, extends, element_x, cache,
+  return Mesh2d{typeof(element_x)}(tree, boxes, extends, element_x, cache,
                 bulkfaceindices, bulkbdryindices, faceindices, bdryindices, offsets)
 end
 
 
 # TODO Deprecate and remove this
-function Mesh(; N, K, range, basis, periodic)
+function Mesh(; N, K, range, basis, periodic, scheme)
   dims = Int(length(range)/2)
   MType = dims == 1 ? Mesh1d : Mesh2d
-  return MType(; N, K, range, basis, periodic)
+  return MType(; N, K, range, basis, periodic, scheme)
 end
 
 
@@ -265,11 +278,31 @@ end
 #######################################################################
 
 
-function widths(m::Mesh{N}) where N
+function widths(m::Mesh{Element,N}) where {Element,N}
   return NTuple{N,Float64}(abs(r-l) for (l,r) in m.extends)
 end
 
 
+function min_grid_spacing(m::Mesh{SpectralElement})
+  @unpack z = m.element
+  ws = widths(m.boxes[1])
+  minw = minimum(ws)
+  dz = z[2]-z[1]
+  return minw .* dz
+end
+function min_grid_spacing(m::Mesh{FVElement})
+  return minimum(widths(m)) / m.K
+end
+
+
+function isperiodic(m::Mesh)
+  if !all(p -> p == first(m.tree.periodic), m.tree.periodic)
+    TODO("mixed periodic domain")
+  end
+  return all(m.tree.periodic)
+end
+
+
 Base.broadcastable(mesh::Mesh) = Ref(mesh)
 
 
@@ -450,7 +483,7 @@ Base.extrema(mesh::Mesh) = mesh.extends
 npoints(mesh::Mesh)::Int = mesh.element.Npts
 ncells(mesh::Mesh) = prod(mesh.tree.dims)
 layout(mesh::Mesh)::Tuple{Int,Int} = (n_points(mesh), n_cells(mesh))
-dimension(mesh::Mesh{N}) where N = N
+dimension(mesh::Mesh{Element,N}) where {Element,N} = N
 cache(mesh::Mesh) = mesh.cache
 
 @deprecate n_points(mesh::Mesh) npoints(mesh)
@@ -568,14 +601,14 @@ end
 Base.length(it::CellDataIterator) = ncells(it.mesh)
 Base.size(it::CellDataIterator) = size(it.mesh)
 
-@inline function Base.iterate(it::CellDataIterator{Mesh1d}, state=1)
+@inline function Base.iterate(it::CellDataIterator{<:Mesh1d}, state=1)
   state > length(it) && return nothing
   idx = it.mesh.offsets[state]
   Npts = it.mesh.element.Npts
   v = view(it.data, idx+1:idx+Npts)
   return v, state+1
 end
-@inline function Base.iterate(it::CellDataIterator{Mesh2d}, state=1)
+@inline function Base.iterate(it::CellDataIterator{<:Mesh2d}, state=1)
   state > length(it) && return nothing
   idx = it.mesh.offsets[state]
   Npts = it.mesh.element.Npts
diff --git a/src/parameters.jl b/src/parameters.jl
index a7e6c333f8fc642992538754a787af9f482165d3..9d05954faff7ce657428db0ef43ba742a2d8ed75 100644
--- a/src/parameters.jl
+++ b/src/parameters.jl
@@ -495,6 +495,14 @@ end
   basis = "lgl_lumped"
   @check basis in [ "lgl", "glgl", "lgl_lumped" ]
 
+  """
+  Cellwise approximation scheme of solution. Available options
+  - `DG` ... Discontinuous Galerkin
+  - `FV` ... Finite Volume (central scheme)
+  """
+  scheme = "DG"
+  @check scheme in [ "DG", "FV" ]
+
   """
   Periodicity in cartesian directions `x,y`
   - 1d: [ `x_dir` ]
diff --git a/src/utils.jl b/src/utils.jl
index 8d9731f70d0b9edbb6fb8f495b93c1808c7b82fb..29648728fe86352cfe6a3ad9742d80171ba85132 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -270,7 +270,7 @@ end
 
 
 TODO() = error("Not implemented yet!")
-TODO(msg::String) = error(msg)
+TODO(msg::String) = error("Not implemented: $msg")
 TODO(fn::Function) = error("'$fn': Not implemented yet!")
 TODO(fn::Function, msg) = error("""
 '$fn'' Not implemented yet!
diff --git a/test/IntegrationTests/refs/fv_advection_sine/fv_advection_sine.toml b/test/IntegrationTests/refs/fv_advection_sine/fv_advection_sine.toml
new file mode 100644
index 0000000000000000000000000000000000000000..331a64e5d303fbf3edf03f1d68505586e2eeae18
--- /dev/null
+++ b/test/IntegrationTests/refs/fv_advection_sine/fv_advection_sine.toml
@@ -0,0 +1,18 @@
+[ScalarEq]
+equation = "advection"
+id = "sine"
+
+[Mesh]
+range = [ 0.0, 1.0 ]
+k = 256
+basis = "lgl_lumped"
+scheme = "FV"
+
+[Output]
+every_iteration = 0
+variables  = [ "u" ]
+aligned_ts = [ 0.25, 0.50, 0.75, 1.0, 1.25, 1.50, 1.75, 2.0 ]
+
+[Evolution]
+tend = 2.2
+cfl = 0.9
diff --git a/test/IntegrationTests/refs/fv_advection_sine/output.h5 b/test/IntegrationTests/refs/fv_advection_sine/output.h5
new file mode 100644
index 0000000000000000000000000000000000000000..1813b5e200253908e35815ff827dbcc768cc0233
Binary files /dev/null and b/test/IntegrationTests/refs/fv_advection_sine/output.h5 differ
diff --git a/test/IntegrationTests/refs/testconfig.toml b/test/IntegrationTests/refs/testconfig.toml
index ee94aa62253d0f601a3a91c08d135da02307cc7c..35fc1305fc69d5bf0dfb737f19092d7ddf64b2db 100644
--- a/test/IntegrationTests/refs/testconfig.toml
+++ b/test/IntegrationTests/refs/testconfig.toml
@@ -7,6 +7,9 @@ default_abstol = 1e-8
 [advection_sine]
 variables = [ "u" ]
 
+[fv_advection_sine]
+variables = [ "u" ]
+
 [advection_sine_glgl]
 variables = [ "u" ]