diff --git a/src/bdryconditions.jl b/src/bdryconditions.jl
index 41b5d293f05939bfdc95b96afc574fc0adc552dd..7644897acf4f001062e30e172d4742918ed542b9 100644
--- a/src/bdryconditions.jl
+++ b/src/bdryconditions.jl
@@ -112,12 +112,12 @@ isperiodic(bdry::BoundaryConditions) = bdry.isperiodic
 
 
 """
-    broadcast_boundaryconditions!(f, bc::BoundaryConditions, cache::Cache, mesh::Mesh, t)
+    broadcast_boundaryconditions!(f, bc::BoundaryConditions, cache::Cache, mesh::AbstractMesh, t)
 
 Impose boundary conditions using the flux method `f` and boundary conditions `bc`
 for fields in `cache` on `mesh` at time `t`.
 """
-function broadcast_boundaryconditions!(f, bc::BoundaryConditions, cache::Cache, mesh::Mesh, t)
+function broadcast_boundaryconditions!(f, bc::BoundaryConditions, cache::Cache, mesh::AbstractMesh, t)
   @unpack lhs_bc, rhs_bc, rsolver = bc
   has_signature(f, rsolver) || error("broadcast_boundary!: Function '$f' has no signature!")
   bc.isperiodic && return
@@ -148,7 +148,7 @@ boundary_state(bc::OutflowBC, u_inner, t)  = u_inner
 
 # similar to _broadcast_volume!, but only iterating domain bdrys and broadcasting
 # over two returned data arrays
-function _broadcast_boundaries!(f, cache::Cache, mesh::Mesh,
+function _broadcast_boundaries!(f, cache::Cache, mesh::AbstractMesh,
     accepts_vars::NTuple{Na,Vector{Float64}}, 
     returns_lhs_data::NTuple{Nr,Vector{Float64}},
     returns_rhs_data::NTuple{Nr,Vector{Float64}},
diff --git a/src/cache.jl b/src/cache.jl
index 4b6e19a3b0570d7ef7bc8faa9c210ca91f9d238f..54cd2858ea2a5811ca4771b48c9756962d65d132 100644
--- a/src/cache.jl
+++ b/src/cache.jl
@@ -44,7 +44,7 @@ struct Cache
 end
 
 
-function Cache(mesh::Mesh)
+function Cache(mesh::AbstractMesh)
   verbose = false
   variables_to_fields = Dict{Symbol, Symbol}()
   fields_to_variables = Dict{Symbol, Vector{Symbol}}()
diff --git a/src/callbacks.jl b/src/callbacks.jl
index a08054925c72324792b176cee1808c8e6dba2fc8..739aa288f9d781fab77b3c70eabd5b155dc6efaf 100644
--- a/src/callbacks.jl
+++ b/src/callbacks.jl
@@ -90,7 +90,8 @@ Base.show(io::IO, cbt::CallbackTiming) = print(io, state(cbt))
     istimed!(timing::CallbackTiming, t)
 
 Check if callback `timing` should trigger at time `t`.
-If `true`, `timing` is updated and calling `istimed!(timing, t)` again will yield `false`.
+If `true`, `timing` is updated and calling `istimed!(timing, t)` again will yield `false`,
+unless `timing.force = true`.
 """
 function istimed!(timing::CallbackTiming, t)
   @unpack every_iteration, every_dt = timing
@@ -178,9 +179,9 @@ struct FunctionCallback{T_Function<:Function} <: AbstractCallback
 end
 
 
-function (cb::FunctionCallback)(u, t, iter; force=false)
+function (cb::FunctionCallback)(u, t, iter)
   @unpack func, timing = cb
-  if force || istimed!(timing, t)
+  if istimed!(timing, t)
     func(u, t)
     return true
   end
@@ -206,7 +207,7 @@ Base.show(io::IO, fcb::FunctionCallback) = print(io, state(fcb))
 
 
 """
-    SaveCallback(variables::Symbol..., C::Cache, mesh::Mesh, filename, timing::CallbackTiming)
+    SaveCallback(variables::Symbol..., C::Cache, mesh::AbstractMesh, filename, timing::CallbackTiming)
 
 Make callback for saving `variables...` from `C` to file `filename`.
 """
@@ -216,7 +217,7 @@ mutable struct SaveCallback <: AbstractCallback
   filename::String
   variablenames::Vector{Tuple{Symbol,Symbol}}
 
-  function SaveCallback(variables::Vector{Symbol}, C::Cache, mesh::Mesh, filename,
+  function SaveCallback(variables::Vector{Symbol}, C::Cache, mesh::AbstractMesh, filename,
       timing::CallbackTiming)
 
     # normalize and sanitize filename
@@ -402,7 +403,10 @@ struct TimeAlignedSaveCallback <: AbstractCallback
     # time step is computed before the time step is taken after which saving should happen
     function save_on_time(u, t)
       if schedule_save
+        prev_force = scb.fcb.timing.force
+        scb.fcb.timing.force = true # force evaluation
         scb(u, t, 0)
+        scb.fcb.timing.force = prev_force # reset force flag
         schedule_save = false
       end
     end
diff --git a/src/dg1d.jl b/src/dg1d.jl
index 8aead18e4ce59dc52f3249b24cf11888b2bc1750..c263cdaa55fadb6e749ccdde024bd6d224cd5175 100644
--- a/src/dg1d.jl
+++ b/src/dg1d.jl
@@ -56,7 +56,6 @@ module dg1d
   include("cache.jl")
   include("callbacks.jl")
 
-  include("numericalflux.jl")
   include("riemannsolver.jl")
   include("bdryconditions.jl")
   include("dg_rhs.jl")
diff --git a/src/dg_rhs.jl b/src/dg_rhs.jl
index 3f4fbf1cc597feb44f1c6adf39dd80f9ecc4afe7..c980efa1f9f5120475c966dffe64a7afff014801 100644
--- a/src/dg_rhs.jl
+++ b/src/dg_rhs.jl
@@ -6,9 +6,9 @@ export compute_rhs_weak_form!,
     compute_rhs_weak_form!(rhs, f, nf_lhs, nf_rhs, mesh)
     compute_rhs_weak_form!(rhs, f, s, nf_lhs, nf_rhs, mesh)
 """
-function compute_rhs_weak_form!(rhs, f, nf_lhs, nf_rhs, mesh::Mesh)
+function compute_rhs_weak_form!(rhs, f, nf_lhs, nf_rhs, mesh::AbstractMesh)
   @unpack invjac, element = mesh
-  @unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
+  @unpack MDM, Ml_lhs, Ml_rhs, Npts = element
   @unpack K = mesh
   shape        = layout(mesh)
   mat_rhs      = reshape(view(rhs,:), shape)
@@ -23,7 +23,7 @@ function compute_rhs_weak_form!(rhs, f, nf_lhs, nf_rhs, mesh::Mesh)
   mat_rhs    .*= invjac
   return
 end
-function compute_rhs_weak_form!(rhs, f, s, nf_lhs, nf_rhs, mesh::Mesh)
+function compute_rhs_weak_form!(rhs, f, s, nf_lhs, nf_rhs, mesh::AbstractMesh)
   compute_rhs_weak_form!(rhs, f, nf_lhs, nf_rhs, mesh)
   rhs .+= s
   return
@@ -34,15 +34,15 @@ end
     compute_rhs_strong_form!(rhs, f, nf_lhs, nf_rhs, mesh)
     compute_rhs_strong_form!(rhs, f, s, nf_lhs, nf_rhs, mesh)
 """
-function compute_rhs_strong_form!(rhs, f, nf_lhs, nf_rhs, mesh::Mesh)
+function compute_rhs_strong_form!(rhs, f, nf_lhs, nf_rhs, mesh::AbstractMesh)
   @unpack invjac, element = mesh
-  @unpack Npts, invM, D, Ml_lhs, Ml_rhs = element
+  @unpack Npts, D, Ml_lhs, Ml_rhs = element
   @unpack K = mesh
   shape        = layout(mesh)
   mat_rhs      = reshape(view(rhs,:), shape)
   mat_f        = reshape(view(f,:), shape)
   # compute rhs = - D * f inplace to rhs
-  # benchmarking showed that splitting of the factor of -1 is faster 
+  # benchmarking showed that splitting of the factor of -1 is faster
   # because a temporary is not created then
   mul!(mat_rhs, D, mat_f)
   mat_rhs     .*= -1
@@ -50,7 +50,7 @@ function compute_rhs_strong_form!(rhs, f, nf_lhs, nf_rhs, mesh::Mesh)
   mat_rhs    .*= invjac
   return
 end
-function compute_rhs_strong_form!(rhs, f, s, nf_lhs, nf_rhs, mesh::Mesh)
+function compute_rhs_strong_form!(rhs, f, s, nf_lhs, nf_rhs, mesh::AbstractMesh)
   compute_rhs_strong_form!(rhs, f, nf_lhs, nf_rhs, mesh)
   rhs .+= s
   return
diff --git a/src/flux.jl b/src/flux.jl
deleted file mode 100644
index 53c32cd110dc0162b2a98a15f1d3017b4de85a43..0000000000000000000000000000000000000000
--- a/src/flux.jl
+++ /dev/null
@@ -1,340 +0,0 @@
-export NumericalFluxes
-
-
-#######################################################################
-#                                Types                                #
-#######################################################################
-
-
-@doc """
-  NumericalFluxes{TVar<:AbstractVariables}
-
-Object for storing numerical fluxes appearing on `lhs, rhs` element faces.
-"""
-struct NumericalFluxes
-  lhs::Variables
-  rhs::Variables
-  mesh::Mesh
-end
-
-
-#######################################################################
-#                            Constructors                             #
-#######################################################################
-
-
-function NumericalFluxes(mesh::Mesh, names::Symbol...)
-  _, ncells = layout(mesh)
-  lhs, rhs = [ Variables(ncells, names...) for _ in 1:2 ]
-  return NumericalFluxes(lhs, rhs, mesh)
-end
-
-
-#######################################################################
-#                              Constants                              #
-#######################################################################
-
-
-@doc """
-Outward pointing unit normal vectors for lhs/rhs edges of 1d cells, e.g.
-
-              <--x nv_lhs           x--> nv_rhs
------|-----|-----|-----|-----|------|-----|-----|-----|-----> x
-"""
-const nv_lhs = -1.0;
-const nv_rhs =  1.0;
-
-
-#######################################################################
-#                              Deprecacted                              #
-#######################################################################
-
-
-# @doc """
-#   flux_lax_friedrich!(numfluxes::NumericalFluxes{T}, fluxes::T, conservatives::T, speeds::T)
-#
-# Compute Lax Friedrich numerical flux inplace.
-# """
-# function numericalflux_lax_friedrich!(numfluxes::NumericalFluxes{T},
-#     fluxes::T, conservatives::T, speeds::T) where T<:AbstractVariables
-#   @unpack rhs, lhs, mesh = numfluxes
-#   shape = layout(mesh)
-#   npts, ncells = shape
-#   vars = fieldnames(T)
-#   for var in vars
-#     mat_flux = reshape(getfield(fluxes, var), shape)
-#     mat_cons = reshape(getfield(conservatives, var), shape)
-#     mat_speed = reshape(getfield(speeds, var), shape)
-#     numflux_lhs, numflux_rhs = getfield(lhs, var), getfield(rhs, var)
-#     for idx in 1:ncells
-#       idx_lhs, idx_rhs = [ dg1d.periodic_index(i, ncells) for i in [ idx-1, idx ] ]
-#       f_lhs, f_rhs = mat_flux[end,idx_lhs], mat_flux[1,idx_rhs]
-#       u_lhs, u_rhs = mat_cons[end,idx_lhs], mat_cons[1,idx_rhs]
-#       C = maximum(mat_speed[:,idx])
-#       avg_f = (f_lhs + f_rhs) / 2
-#       dif_u = u_lhs - u_rhs
-#       nf_lhs, nf_rhs = avg_f + C * dif_u * nv_lhs, avg_f + C * dif_u * nv_rhs
-#       numflux_lhs[idx], numflux_rhs[idx] = nf_lhs, nf_rhs
-#     end
-#   end
-#   return
-# end
-
-
-#######################################################################
-#                             Deprecated                              #
-#######################################################################
-
-
-@doc """
-Collect interface values, lhs/rhs refer to the lhs/rhs of the interface.
-"""
-struct Interface
-  lhs
-  rhs
-end
-
-
-@doc """
-Evaluate fun with arguments intf and return result in a new Interface struct.
-"""
-function Interface(fun, intf::Interface...)
-  intf_lhs = fun.(map(i->i.lhs, collect(intf))...)
-  intf_rhs = fun.(map(i->i.rhs, collect(intf))...)
-  return Interface(intf_lhs, intf_rhs)
-end
-
-
-@doc """
-Project u onto all interfaces of mesh.
-Periodic bdry conditions are assumed to set the ghost points.
-"""
-function to_interface(u, mesh::Mesh)
-  @unpack element, K = mesh
-  @unpack l_lhs, l_rhs, quadr = element
-  F = K + 1   # number of interfaces
-
-  lhs = similar(u, F)
-  rhs = similar(u, F)
-  if (quadr == :LGL || quadr == :LGL_lumped || quadr == :GLGL)
-    # no interpolation needed
-    lhs[2:end]    .= u[end,:]
-    rhs[1:end-1]  .= u[1,:]
-  else
-    error("We don't handle grids without boundary points yet!")
-    # lhs[2:end]    .= transpose(transpose(l_rhs) * u)
-    # rhs[1:end-1]  .= transpose(transpose(l_lhs) * u)
-  end # if
-  # assume periodic grid
-  lhs[1]        = lhs[end]
-  rhs[end]      = rhs[1]
-
-  intf = Interface(lhs, rhs)
-  return intf
-end
-
-
-@doc """
-Compute average between lhs and rhs interface values.
-"""
-function average(intf::Interface)
-  return 0.5 * (intf.lhs + intf.rhs)
-end
-
-
-@doc """
-Compute difference between lhs and rhs interface values.
-"""
-function difference(intf::Interface)
-  return intf.lhs - intf.rhs
-end
-
-
-@doc """
-Lax-Friedrich flux, see Hesthaven 2007, section 2.2.
-C ... parameter to control upwinding contribution
-      C = 0 turns this into central flux
-      Common choices include absolute maxmium wave speed on global or
-      local grid.
-"""
-function lax_friedrich_flux(avg_f, dif_u, normal, C)
-  return avg_f + 0.5 * C * dif_u * normal
-end
-
-
-@doc """
-Compute Lax-Friedrich flux times normal vector evaluated at lhs/rhs ends
-of the subdomains.
-C ... parameter to control upwinding contribution
-      C = 0 turns this into central flux
-      Common choices include absolute maxmium wave speed on global or
-      local grid.
-Flux is such that the conservation identity holds, 
-cf. PhD thesis Zingan 2012, p. 77
-     F(u_m, u_p, nv_m) = F(u_p, u_m, nv_p)
-"""
-function lax_friedrich_flux(intf_f::Interface, intf_u::Interface, C)
-  avg_f = average(intf_f)
-  dif_u = difference(intf_u)
-  nf = lax_friedrich_flux.(avg_f, dif_u, nv_rhs, C)
-  nf_lhs = nf[1:end-1]
-  nf_rhs = nf[2:end]
-  return nf_lhs, nf_rhs
-end
-
-
-@doc """
-Special case of lax_friedrich_flux with C = 0.
-Flux is such that the conservation identity holds, 
-cf. PhD thesis Zingan 2012, p. 77
-     F(u_m, u_p, nv_m) = F(u_p, u_m, nv_p)
-"""
-function central_flux(intf_u::Interface)
-  nu = @. 0.5 * ( intf_u.lhs + intf_u.rhs )
-  nu_lhs = nu[1:end-1]
-  nu_rhs = nu[2:end]
-  return nu_lhs, nu_rhs
-end
-
-
-@doc """
-Flux based on interior values.
-Useful for local DG method.
-"""
-function inner_flux(intf::Interface)
-  nf_lhs = intf.rhs[1:end-1]
-  nf_rhs = intf.lhs[2:end]
-  return nf_lhs, nf_rhs
-end
-
-
-@doc """
-Flux based on exterior values.
-Useful for local DG method.
-"""
-function outer_flux(intf::Interface)
-  nf_lhs = intf.lhs[1:end-1]
-  nf_rhs = intf.rhs[2:end]
-  return nf_lhs, nf_rhs
-end
-
-
-@doc """
-Local DG flux according to Hesthaven, Warburton 2007, p. 252.
-"""
-function ldg_flux_asym_lhs(intf::Interface)
-  nf_lhs = intf.lhs[1:end-1]
-  nf_rhs = intf.lhs[2:end]
-  return nf_lhs, nf_rhs
-end
-
-
-@doc """
-Local DG flux according to Hesthaven, Warburton 2007, p. 252.
-"""
-function ldg_flux_asym_rhs(intf::Interface)
-  nf_lhs = intf.rhs[1:end-1]
-  nf_rhs = intf.rhs[2:end]
-  return nf_lhs, nf_rhs
-end
-
-
-@doc """
-  um ... cell interior value of u
-  up ... cell exterior value of u
-  nm ... cell outward pointing unit vector (at this bdry)
-  β  ... 'upwinding direction'
-"""
-function ldg_flux(um, up, nm, β)
-  avg_u = 0.5 * (um + up)
-  dif_u = um - up
-  return avg_u + nm * β * dif_u
-end
-
-
-@doc """
-Symmetrized local DG flux according to Hesthaven, Warburton 2007, p. 252
-Reduces to local DG rhs/lhs versions from above if β = ±1/2.
-"""
-function ldg_flux_sym(intf::Interface, β)
-  nf_lhs =  ldg_flux.(intf.rhs[1:end-1], intf.lhs[1:end-1], nv_lhs, β)
-  nf_rhs =  ldg_flux.(intf.lhs[2:end], intf.rhs[2:end], nv_rhs, β)
-  return nf_lhs, nf_rhs
-end
-
-
-@doc """
-  qm ... cell interior value of q
-  qp ... cell exterior value of q
-  um ... cell interior value of u
-  up ... cell exterior value of u
-  nm ... cell outward pointing unit vector (at this bdry)
-  β  ... 'upwinding direction'
-"""
-function ldg_flux_penalty(qm, qp, um, up, nm, β, τ)
-  avg_q = 0.5 * (qm + qp)
-  dif_q = qm - qp
-  dif_u = um - up
-  return avg_q - nm * β * dif_q - nm * τ * dif_u
-end
-
-
-@doc """
-Symmetrized and penalized local DG flux according to Hesthaven, 
-Warburton 2007, p. 267
-"""
-function ldg_flux_sym_penalty( intf_q::Interface, intf_u::Interface, β, τ)
-  nf_lhs =  ldg_flux_penalty.(
-    intf_q.rhs[1:end-1], intf_q.lhs[1:end-1], 
-    intf_u.rhs[1:end-1], intf_u.lhs[1:end-1],
-    nv_lhs, β, τ)
-  nf_rhs =  ldg_flux_penalty.(
-    intf_q.lhs[2:end], intf_q.rhs[2:end], 
-    intf_u.lhs[2:end], intf_u.rhs[2:end],
-    nv_rhs, β, τ)
-  return nf_lhs, nf_rhs
-end
-
-
-@doc """
-Compute boundary flux from dirichlet boundary condition at given bdry.
-"""
-function bdry_flux(u, t, bdry::BDRY, bc::dg1d.Dirichlet_BC)
-  return bc(t)
-end
-
-
-@doc """
-Compute boundary flux from outflow boundary condition at given bdry.
-"""
-function bdry_flux(u, t, bdry::BDRY, bc::dg1d.Outflow_BC)
-  bdryflux = 0.0
-  if bdry == LHS
-    bdryflux = u[1]
-  else
-    bdryflux = u[end]
-  end
-  return bdryflux
-end
-
-
-@doc """
-Compute fluxes for given boundary conditions at both boundaries.
-"""
-function bdry_flux(u, t, bdryconds::dg1d.BdryConditions)
-  bflux_lhs = bdry_flux(u, t, LHS, bdryconds.lhs)
-  bflux_rhs = bdry_flux(u, t, RHS, bdryconds.rhs)
-  return bflux_lhs, bflux_rhs
-end
-
-
-@doc """
-Compute fluxes for LDG variable for given boundary conditions.
-
-TODO Dispatch this on bc if we ever happen to need vonNeumann conditions.
-"""
-function bdry_flux_LDG_variable(q, t, bc::dg1d.Abstract_BdryConditions)
-  bdryflux_lhs = 0.0
-  bdryflux_rhs = 0.0
-  return bdryflux_lhs, bdryflux_rhs
-end
diff --git a/src/mesh.jl b/src/mesh.jl
index 4533fc16fb5ab268a52e9648550468a789be587d..6ba7451ecd1dbfd5a583ac0eca6d6197a845bc6f 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -17,13 +17,16 @@ export Mesh, MeshInterpolator,
 #######################################################################
 
 
+abstract type AbstractMesh end
+
+
 @doc """
-    Mesh
+    Mesh <: AbstractMesh
 
 Object that bundles grid and operators needed for modeling a Discontinuous Galerkin 
 scheme in 1D.
 """
-struct Mesh
+struct Mesh <: AbstractMesh
   element::SpectralElement  # collection of spectral element operators
   K::Int64                  # number of cells
   xmin::Float64             # left most grid coordinate
@@ -106,26 +109,26 @@ end
 #######################################################################
 
 
-Base.broadcastable(mesh::Mesh) = Ref(mesh)
+Base.broadcastable(mesh::AbstractMesh) = Ref(mesh)
 
 
 """
-    extract_cell_bdrys(mesh::Mesh)
+    extract_cell_bdrys(mesh::AbstractMesh)
 
 Return coordinates of cell boundaries from a `mesh`.
 """
-function extract_cell_bdrys(mesh::Mesh)
+function extract_cell_bdrys(mesh::AbstractMesh)
   @unpack x = mesh
   return vcat(x[1,:], [x[end,end]])
 end
 
 
 @doc """
-  `grid_average(u, mesh::Mesh)`
+    grid_average(u, mesh::AbstractMesh)
 
-Return the grid average of all cell averages 
+Return the grid average of all cell averages.
 """
-function grid_average(u, mesh::Mesh)
+function grid_average(u, mesh::AbstractMesh)
   @unpack element, dl, L = mesh
   cellavg = cellwise_average(u, mesh) .* dl
   avg = sum(cellavg) / L
@@ -134,11 +137,11 @@ end
 
 
 @doc """
-  `cellwise_average(u, mesh::Mesh)`
+    cellwise_average(u, mesh::AbstractMesh)
 
 Return the cellwise average of `u` over `mesh.x`.
 """
-function cellwise_average(u, mesh::Mesh)
+function cellwise_average(u, mesh::AbstractMesh)
   @unpack element, invjac, K, dl = mesh
   @unpack Npts = element
   v = ones(Npts)
@@ -152,11 +155,11 @@ end
 
 
 @doc """
-  `cellwise_inner_product(u, v, mesh::Mesh)`
+    cellwise_inner_product(u, v, mesh::AbstractMesh)
 
 Return vector of inner products between `u, v` on every cell in `mesh`.
 """
-function cellwise_inner_product(u, v, mesh::Mesh)
+function cellwise_inner_product(u, v, mesh::AbstractMesh)
   @unpack element, K, invjac = mesh
   @unpack Npts = element
   L = layout(mesh)
@@ -169,51 +172,55 @@ end
 
 
 @doc """
-  `broken_inner_product(u, v, mesh::Mesh)`
-  
+    broken_inner_product(u, v, mesh::AbstractMesh)
+
 Return sum of all inner products between `u,v` on every cell in `mesh`.
 """
-function broken_inner_product(u, v, mesh::Mesh)
+function broken_inner_product(u, v, mesh::AbstractMesh)
   inner_products = cellwise_inner_product(u, v, mesh)
   total = sum(inner_products[:])
   return total
 end
 
 
-function LinearAlgebra.norm(u, mesh::Mesh, order)
+function LinearAlgebra.norm(u, mesh::AbstractMesh, order)
   @match order begin
     1 => norm_L1(u, mesh)
     2 => norm_L2(u, mesh)
     _ => error("Norm for order = $order not implemented.")
   end
 end
-norm_L2(u, mesh::Mesh) = sqrt(broken_inner_product(u, u, mesh))
-norm_L1(u, mesh::Mesh) = broken_inner_product(abs.(u), ones(Float64, size(u)), mesh)
+norm_L2(u, mesh::AbstractMesh) = sqrt(broken_inner_product(u, u, mesh))
+norm_L1(u, mesh::AbstractMesh) = broken_inner_product(abs.(u), ones(Float64, size(u)), mesh)
 
 
 @doc """
+    local_maxabs(v::AbstractMatrix, mesh::AbstractMesh)
+
 Compute local maximum absolute value.
 """
-function local_maxabs(v::AbstractMatrix, mesh::Mesh) 
+function local_maxabs(v::AbstractMatrix, mesh::AbstractMesh)
   return vec(maximum(abs.(v), dims=1))
 end
 
 
 @doc """
+    local_maxabs(v::AbstractVector, mesh::AbstractMesh)
+
 Compute local minimum absolute value.
 """
-function local_maxabs(v::AbstractVector, mesh::Mesh) 
+function local_maxabs(v::AbstractVector, mesh::AbstractMesh)
   return local_maxabs(reshape(v, (mesh.element.Npts, mesh.K)))
 end
 
 
 """
-    find_nearest(xi, mesh::Mesh)
+    find_nearest(xi, mesh::AbstractMesh)
 
 Find index of node in `mesh` that lies closest to `xi`.
 Returns the index of the outermost point if `xi` falls outside `mesh`.
 """
-function find_nearest(xi, mesh::Mesh)
+function find_nearest(xi, mesh::AbstractMesh)
   @unpack x = mesh
   min_index = 1
   min_distance = abs(x[min_index] - xi)
@@ -228,28 +235,28 @@ function find_nearest(xi, mesh::Mesh)
 end
 
 
-@deprecate locate_point(xi, mesh::Mesh) find_nearest(xi, mesh::Mesh)
+@deprecate locate_point(xi, mesh::AbstractMesh) find_nearest(xi, mesh::AbstractMesh)
 
 
 """
-    find_cell_index(xi, mesh::Mesh)
+    find_cell_index(xi, mesh::AbstractMesh)
 
 Find the cell index into which coordinate `xi` falls into the `mesh`.
 Note that it returns the first index found in case `xi` falls onto an interface boundary.
 Returns nothing if `xi` lies outside the `mesh`.
 """
-function find_cell_index(xi, mesh::Mesh)
+function find_cell_index(xi, mesh::AbstractMesh)
   cell_bdrys = extract_cell_bdrys(mesh)
   cell_index = findfirst(idx -> cell_bdrys[idx] <= xi <= cell_bdrys[idx+1], 1:length(cell_bdrys)-1)
   return cell_index
 end
 
 
-Base.length(mesh::Mesh)::Int64 = length(mesh.x)
-Base.extrema(mesh::Mesh) = extrema(mesh.x)
-n_points(mesh::Mesh)::Int = mesh.element.Npts
-n_cells(mesh::Mesh) = mesh.K
-layout(mesh::Mesh)::Tuple{Int,Int} = (n_points(mesh), n_cells(mesh))
+Base.length(mesh::AbstractMesh)::Int64 = length(mesh.x)
+Base.extrema(mesh::AbstractMesh) = extrema(mesh.x)
+n_points(mesh::AbstractMesh)::Int = mesh.element.Npts
+n_cells(mesh::AbstractMesh) = mesh.K
+layout(mesh::AbstractMesh)::Tuple{Int,Int} = (n_points(mesh), n_cells(mesh))
 
 
 """
@@ -274,15 +281,15 @@ julia> meshintrp(intrp_data, reference_data) # result is stored inside intrp_dat
 
 Note that for inplace interpolation `intrp_data` and `reference_data` cannot be aliased.
 """
-struct MeshInterpolator
-  reference_mesh::Mesh
-  sample_mesh::Mesh
+struct MeshInterpolator{T_Mesh<:AbstractMesh}
+  reference_mesh::T_Mesh
+  sample_mesh::T_Mesh
   interpolation_matrices::Vector{Tuple{Int64,Matrix{Float64}}}
   reference_cell_indices::Vector{Int64}
 end
 
 
-function MeshInterpolator(smpl_mesh::Mesh, ref_mesh::Mesh)
+function MeshInterpolator(smpl_mesh::T_Mesh, ref_mesh::T_Mesh) where T_Mesh<:AbstractMesh
 
   smpl_xmin, smpl_xmax = extrema(smpl_mesh)
   ref_xmin, ref_xmax   = extrema(ref_mesh)
@@ -307,7 +314,7 @@ function MeshInterpolator(smpl_mesh::Mesh, ref_mesh::Mesh)
     # extract the sample nodes
     smpl_xs   = smpl_x[smpl_idxs]
     # generate interpolation matrix
-    push!(interpolation_matrices, 
+    push!(interpolation_matrices,
           (cell_idx, dg1d.lagrange_interp_matrix(smpl_xs, ref_x[:,cell_idx])))
 
     current_idx += length(smpl_idxs)
@@ -316,7 +323,7 @@ function MeshInterpolator(smpl_mesh::Mesh, ref_mesh::Mesh)
     end
   end
 
-  return MeshInterpolator(ref_mesh, smpl_mesh, interpolation_matrices, ref_cell_indices)
+  return MeshInterpolator{T_Mesh}(ref_mesh, smpl_mesh, interpolation_matrices, ref_cell_indices)
 end
 
 
diff --git a/src/numericalflux.jl b/src/numericalflux.jl
deleted file mode 100644
index e3eb2c89f3ebee92faeeac7690c5b2570692b532..0000000000000000000000000000000000000000
--- a/src/numericalflux.jl
+++ /dev/null
@@ -1,69 +0,0 @@
-export numericalflux_lax_friedrich_local!,
-       numericalflux_central!,
-       numericalflux_strongform!
-
-
-function numericalflux_lax_friedrich_local!(nf_lhs, nf_rhs, f, u, local_v, mesh::Mesh)
-  shape = layout(mesh)
-  _, K = shape
-  Idx = LinearIndices(shape)
-  # domain bdry interfaces are assumed to be periodically connected
-  for k = 1:K # interface index
-    k_lhs, k_rhs = periodic_index(k-1, K), periodic_index(k, K)
-    # values **to** the lhs/rhs of kth interface
-    fbdry_lhs, fbdry_rhs = f[Idx[end,k_lhs]], f[Idx[1,k_rhs]]
-    ubdry_lhs, ubdry_rhs = u[Idx[end,k_lhs]], u[Idx[1,k_rhs]]
-    v = max(abs(local_v[k_lhs]), abs(local_v[k_rhs])) # take maximum in surrounding cells
-    # Convention:
-    # nf_lhs[k] corresponds to the numerical flux at the lhs end of the kth cell
-    # Sketch
-    # numerical flux:           nf_lhs[k]  nf_rhs[k]
-    # cell idx:               k-1   ↓    k    ↓    k+1
-    #                         ↓          ↓         ↓
-    # |---------|---------|---------|---------|---------|---------|---------|
-    #                       (k_lhs)   (k_rhs)
-    #                     ↑         ↑         ↑
-    # interface idx:     k-1        k        k+1
-    #
-    avg_fbdry = 0.5 * (fbdry_lhs + fbdry_rhs)
-    dif_ubdry = ubdry_lhs - ubdry_rhs
-    nf_lhs[k_rhs] = avg_fbdry + 0.5 * v * dif_ubdry
-    nf_rhs[k_lhs] = avg_fbdry + 0.5 * v * dif_ubdry
-  end
-  return
-end
-
-
-function numericalflux_central!(nf_lhs, nf_rhs, u, mesh::Mesh)
-  shape = layout(mesh)
-  Idx = LinearIndices(shape)
-  _, K = shape
-  # domain bdry interfaces are assumed to be periodically connected
-  for k = 1:K # interface index
-    k_lhs, k_rhs = periodic_index(k-1, K), periodic_index(k, K)
-    # values to the lhs/rhs of kth interface
-    ubdry_lhs, ubdry_rhs = u[Idx[end,k_lhs]], u[Idx[1,k_rhs]]
-    avg_ubdry = 0.5 * (ubdry_lhs + ubdry_rhs)
-    # for index convention see numericalflux_lax_friedrich_local!
-    nf_lhs[k_rhs] = avg_ubdry
-    nf_rhs[k_lhs] = avg_ubdry
-  end
-  return
-end
-
-
-function numericalflux_strongform!(nf_lhs, nf_rhs, f, mesh::Mesh)
-  shape = layout(mesh)
-  Idx = LinearIndices(shape)
-  _, K = shape
-  # domain bdry interfaces are assumed to be periodically connected
-  for k = 1:K # interface index
-    k_lhs, k_rhs = periodic_index(k-1, K), periodic_index(k, K)
-    # values to the lhs/rhs of kth interface
-    fbdry_lhs, fbdry_rhs = f[Idx[end,k_lhs]], f[Idx[1,k_rhs]]
-    # for index convention see numericalflux_lax_friedrich_local!
-    nf_lhs[k_rhs] = fbdry_rhs - nf_lhs[k_rhs]
-    nf_rhs[k_lhs] = fbdry_lhs - nf_rhs[k_lhs]
-  end
-  return
-end
diff --git a/src/rhs.jl b/src/rhs.jl
deleted file mode 100644
index a39e072bcba52e19b3866625b9541f32e0f4f371..0000000000000000000000000000000000000000
--- a/src/rhs.jl
+++ /dev/null
@@ -1,76 +0,0 @@
-@doc raw"""
-Compute RHS of weak DG approximation of a conservation law:
-```math
-  \partial_t u = - \partial_x f(u)
-  \rightarrow
-  u_i = M^{-1}_{ij} S^{jk} f_k 
-        - M^{-1}_{ij} ( l_j(LHS) (nf)^{*}(LHS) + l_j(RHS) (nf)^{*}(RHS) )
-```
-"""
-function rhs_weak!(du, f, nf_lhs, nf_rhs, mesh::Mesh)
-  # extract operators
-  @unpack invjac, element = mesh
-  @unpack MDM, Ml_lhs, Ml_rhs = element
-  # evaluate rhs
-  mul!(du, MDM, f)
-  du .-= (transpose(nf_rhs) .* Ml_rhs .- transpose(nf_lhs) .* Ml_lhs)
-  # account for coordinate transformation onto physical grid
-  du .*= invjac
-end # function
-
-
-@doc raw"""
-Compute RHS of strong DG approximation of a conservation law:
-```math
-  \partial_t u = - \partial_x f(u)
-  \rightarrow
-  u_i = - D^{jk} u_k 
-        + M^{-1}_{ij} ( l_j(LHS) (nf(LHS) - (nf)^{*}(LHS)) 
-                      + l_j(RHS) (nf(RHS) - (nf)^{*}(RHS)) )
-```
-"""
-function rhs_strong!(du, f, f_min_nf_lhs, f_min_nf_rhs, mesh::Mesh)
-  # extract operators
-  @unpack invjac, element = mesh
-  @unpack D, Ml_lhs, Ml_rhs = element
-  # evaluate rhs
-  mul!(du, D, - f)
-  du .+= (transpose(f_min_nf_rhs) .* Ml_rhs .-
-          transpose(f_min_nf_lhs) .* Ml_lhs)
-  # account for coordinate transformation onto physical grid
-  du .*= invjac
-end # function
-
-
-
-function new_rhs_weak!(vec_du, vec_f, nf_lhs, nf_rhs, mesh::Mesh)
-  # extract operators
-  @unpack invjac, element = mesh
-  @unpack Npts, MDM, Ml_lhs, Ml_rhs = element
-  @unpack K = mesh
-  du = reshape(vec_du, (Npts, K))
-  f = reshape(vec_f, (Npts, K))
-  # evaluate rhs
-  mul!(du, MDM, f)
-  du .-= (transpose(nf_rhs) .* Ml_rhs .- transpose(nf_lhs) .* Ml_lhs)
-  # account for coordinate transformation onto physical grid
-  du .*= invjac
-  return
-end # function
-
-
-
-function new_rhs_weak!(vec_du, vec_f, vec_s, nf_lhs, nf_rhs, mesh::Mesh)
-  @unpack invjac, element = mesh
-  @unpack Npts, invM, MDM, Ml_lhs, Ml_rhs = element
-  @unpack K = mesh
-  du = reshape(vec_du, (Npts, K))
-  f = reshape(vec_f, (Npts, K))
-  s = reshape(vec_s, (Npts, K))
-  mul!(du, MDM, f)
-  @. du += - (nf_rhs' * Ml_rhs - nf_lhs' * Ml_lhs)
-  # account for coordinate transformation onto physical grid
-  du .*= invjac
-  du .+= s
-  return
-end # function
diff --git a/src/riemannsolver.jl b/src/riemannsolver.jl
index 1d857c01ce021a0aa8e0c3aff0137d52ed0a763b..d993d17efd3001dbae65e472fde676fdff0fb9cf 100644
--- a/src/riemannsolver.jl
+++ b/src/riemannsolver.jl
@@ -8,14 +8,14 @@ export AbstractRiemannSolver,
 
 
 """
-    abstract AbstractRiemannSolver 
-    
+    abstract AbstractRiemannSolver
+
 Supertype for `RiemannSolver` interface.
 Any Subtype `Solver <: RiemannSolver` must implement the following methods
 - `register_variables!(cache, rsolver::Solver)`,
-- `broadcast_faces!(f, rsolver::Solver, cache::Cache, mesh::Mesh)` where `f` is a 
+- `broadcast_faces!(f, rsolver::Solver, cache::Cache, mesh::AbstractMesh)` where `f` is a
   `@with_signature` function.
-- `broadcast_boundaries!(f, rsolver::Solver, cache::Cache, mesh::Mesh)` where `f` is a 
+- `broadcast_boundaries!(f, rsolver::Solver, cache::Cache, mesh::AbstractMesh)` where `f` is a
   `@with_signature` function.
 
 See: `ApproxRiemannSolver`
@@ -111,12 +111,12 @@ end
 
 
 """
-    broadcast_faces!(f, rsolver::ApproxRiemannSolver, cache::Cache, mesh::Mesh)
+    broadcast_faces!(f, rsolver::ApproxRiemannSolver, cache::Cache, mesh::AbstractMesh)
 
 Broadcast `@with_signature` function `f` with dispatch `rsolver` over `cache` and
 faces of `mesh`.
 """
-function broadcast_faces!(f, rsolver::ApproxRiemannSolver, cache::Cache, mesh::Mesh)
+function broadcast_faces!(f, rsolver::ApproxRiemannSolver, cache::Cache, mesh::AbstractMesh)
   has_signature(f, rsolver) || error("broadcast_volume!: Function '$f' has no signature!")
   accepts_vars = accepts(f, rsolver)
   returns_vars_lhs, returns_vars_rhs = returns(f, rsolver)
@@ -137,7 +137,7 @@ end
 
 # similar to _broadcast_volume!, but only iterating cell bdrys and broadcasting
 # over two returned data arrays
-function _broadcast_face!(f, cache::Cache, mesh::Mesh,
+function _broadcast_face!(f, cache::Cache, mesh::AbstractMesh,
     accepts_data::NTuple{Na,Vector{Float64}},
     returns_lhs_data::NTuple{Nr,Vector{Float64}},
     returns_rhs_data::NTuple{Nr,Vector{Float64}}) where {Na,Nr}
diff --git a/src/troubled_cells.jl b/src/troubled_cells.jl
deleted file mode 100644
index 3d69776042a683bcbc7ade9890d2906dc1034d7e..0000000000000000000000000000000000000000
--- a/src/troubled_cells.jl
+++ /dev/null
@@ -1,100 +0,0 @@
-module TroubledCells
-
-
-import dg1d: @with_kw, @unpack, Mesh
-
-
-abstract type AbstractTroubledCellIndicator end
-
-
-@with_kw struct MinmodIndicator <: AbstractTroubledCellIndicator
-  "Threshold for troubled cell indication, the lower the more sensitive"
-  M::Float64 = 1e-5
-  @assert M > 0
-end
-
-
-
-function minmod(args::T...)::T where T <: Real
-  signs = sign.(args)
-  if sum(abs.(diff(collect(signs)))) != T(0)
-    return signs[1] * minimum(abs.(args))
-  else
-    return T(0)
-  end
-end
-
-
-function minmod_M(dx::T, M::T, args::T...)::T where T <: Real
-  if abs(args[1]) <= M * dx^2
-    return args[1]
-  else
-    return minmod(args...)
-  end
-end
-
-
-@doc """
-  `find_troubled_cells(indicator::MinmodIndicator, u::AbstractArray, u_averages::AbstractArray, mesh::Mesh)`
-
-Use modified `minmod` function to scan data `u` with averages `u_averages` to detect irregularities
-in gradients. Returns an array of indices to troubled cells.
-"""
-function find_troubled_cells(indicator::MinmodIndicator, 
-    u::AbstractArray, u_averages::AbstractArray, mesh::Mesh)
-
-  @unpack K, dl = mesh
-  @unpack Npts = mesh.element
-  @unpack M = indicator
-  list_trbld_cells = Int64[]
-  shape = (Npts, K)
-  mat_u = reshape(u, shape)
-
-  # skip bydr cells for now
-  for k in 2:K-1 # loop over cells
-    for i in 1:Npts # loop over pts in cell
-      ui = u[i]
-      uavg = u_averages[k]
-      ubdry_lhs = uavg - mat_u[1,k]
-      ubdry_rhs = mat_u[end,k] - uavg
-      delta_lhs = uavg - u_averages[k-1]
-      delta_rhs = u_averages[k+1] - uavg
-
-      # troubled check, cf. (3.64) in Bugner 2017
-      if !(minmod_M(dl, M, ubdry_lhs, delta_lhs, delta_rhs) ≈ ubdry_lhs) ||
-         !(minmod_M(dl, M, ubdry_rhs, delta_lhs, delta_rhs) ≈ ubdry_rhs)
-        push!(list_trbld_cells, k)
-        break
-      end
-    end
-  end
-
-  return list_trbld_cells
-end
-
-
-@with_kw struct ArtificialViscosityIndicator <: AbstractTroubledCellIndicator
-  "Threshold for troubled cell indication, the lower the more sensitive"
-  threshold::Float64 = 0.1
-  @assert threshold >= 0
-end
-
-
-@doc """
-  `find_troubled_cells(indicator::ArtificialViscosityIndicator, u::AbstractArray, mesh::Mesh)`
-
-Flag cell troubled if `u` exceeds indicator.threshold value.
-Returns list of troubled cell indices.
-"""
-function find_troubled_cells(indicator::ArtificialViscosityIndicator,
-  u::AbstractArray, mesh::Mesh)
-
-  idxs = findall(x -> x >= indicator.threshold, u)
-  cells = [ Tuple(idx)[2] for idx in idxs ]
-  unique!(cells)
-
-  return cells
-end
-
-
-end
diff --git a/src/variablesarray.jl b/src/variablesarray.jl
index 522b4f7911f2cec3ced1e8724934f6c7f1793420..a8fcbb0fc21203c39e49dac2916d50c30b705403 100644
--- a/src/variablesarray.jl
+++ b/src/variablesarray.jl
@@ -1,6 +1,3 @@
-import dg1d: Mesh
-
-
 #######################################################################
 #                                Type                                 #
 #######################################################################
diff --git a/test/dg1dTests/src/test_callbacks.jl b/test/dg1dTests/src/test_callbacks.jl
index 35534591dabc5cb533988c6a5b71e4ee180ff9c5..d99042a9bd1fe6278836c158671ad81b0e66cdfe 100644
--- a/test/dg1dTests/src/test_callbacks.jl
+++ b/test/dg1dTests/src/test_callbacks.jl
@@ -78,20 +78,6 @@ end
   fcb(0, 2, 0)
   @test counter == 3
 
-  # force callback evaluation
-  counter = 0
-  fcb2 = FunctionCallback(increment_counter, CallbackTiming(every_iteration=100))
-  fcb2(0, 0, 0)
-  fcb2(0, 1, 0)
-  @test counter == 0
-  fcb2(0, 0, 0, force=true)
-  fcb2(0, 1, 0, force=true)
-  # forcing ignores any previous calls
-  @test counter == 2
-  fcb2(0, 0, 0, force=true)
-  fcb2(0, 1, 0, force=true)
-  @test counter == 4
-
   # force evaluation by setting force flag in CallbackTiming
   counter = 0
   fcb3 = FunctionCallback(increment_counter, CallbackTiming(force=true))
@@ -103,7 +89,6 @@ end
   fcb3(0, 1, 0)
   @test counter == 4
 
-
 end
 
 
@@ -218,7 +203,7 @@ end
 @testset "TimeAlignedSaveCallback" begin
 
   thisdir = Filesystem.abspath(@__DIR__)
-  mockdir = Filesystem.abspath(joinpath(thisdir, "../mockdir"))
+  mockdir = Filesystem.normpath(joinpath(thisdir, "..", "mockdir"))
   ispath(mockdir) && rm(mockdir, recursive=true)
   mkdir(mockdir)
 
@@ -229,7 +214,7 @@ end
   # fixed time step
   timestep_fn(u, t, n) = 1.0
 
-  # artificial rhs method to unit
+  # artificial rhs method
   u0 = zeros(length(mesh)) # initial data
   rhs!(du, u, t) = u .= randn(length(mesh))
 
@@ -239,17 +224,16 @@ end
 
   # request time aligned save callback
   save_times = [ 0.5, 1.3, 3.1 ]
-  atscb = TimeAlignedSaveCallback(save_times, timestep_fn, scb)
-  aligned_timestep_fn = atscb.aligned_timestep_fn
+  tascb = TimeAlignedSaveCallback(save_times, timestep_fn, scb)
+  aligned_timestep_fn = tascb.aligned_timestep_fn
 
-  cbs = CallbackSet(scb, atscb)
+  cbs = CallbackSet(scb, tascb)
 
   # evolve
   dg1d.evolve_ERK(rhs!, u0, aligned_timestep_fn, (0.0, 5.0), :RK4, callback_fullstep=cbs)
 
   # test output
   @test isfile(filename)
-
   ts = h5open(filename) do file
     groups = keys(file)
     filter!(g->g!="mesh", groups)
@@ -264,11 +248,11 @@ end
 
   # explanation for test results below:
   # we scheduled a SaveCallback with every_dt=2.0; if we evolve till t=5.0 then we should
-  # see at least 2 outputs
+  # see at least 3 outputs, e.g. t ∈ [0, 2, 4]
   # furthermore, we scheduled a TimeAlignedSaveCallback with three time stamps before t = 5.0,
-  # hence, we expect to 3 outputs for this one
+  # hence, we expect 3 more outputs at those times
   # in total we should find 5 outputs
-  @test length(ts) == 5
+  @test length(ts) == 6
   @test all(save_times .∈ Ref(ts))
 
   # test for wrong argument type of timestep_fn in constructor