From 5a594e69195fa4675d2cdc7f5c1eaf878589c21c Mon Sep 17 00:00:00 2001
From: Florian Atteneder <florian.atteneder@uni-jena.de>
Date: Tue, 5 Dec 2023 21:17:41 +0000
Subject: [PATCH] Add a finite volume method for ScalarEq (dg/dg1d.jl!75)

* add reftest fv_advection_sine

* migrate all of ScalarEq to use broadcast_*_2! methods

* add isperiodic method for meshes

* fix grid spacing computation for SpectralElement

* fixup parametric types to subtype Mesh1d/2d

* make TODO error more comprehensive

* fixup cfl usage; add abstractions for min grid spacing

* enable add mesh kwarg to Evolution setup in order to enable FV evolution

* add timestep! and rhs! for FVElement

* do away with the State(...) syntax

* add docs to dg_rhs methods

* add FV stepper

* limit dg_rhs.jl methods to Meshes with SpectralElement

* add mesh kwarg to Evolution so that we can use a spearate stepper based on mesh type

* add a scheme parameter to Mesh

* add FVElement type
---
 src/ScalarEq/callbacks.jl                     |  14 ++--
 src/ScalarEq/equation.jl                      |  16 ++--
 src/ScalarEq/rhs.jl                           |  60 ++++++++-----
 src/ScalarEq/setup.jl                         |  12 +--
 src/TCI/threshold.jl                          |   4 +-
 src/dg1d.jl                                   |   9 +-
 src/dg_rhs.jl                                 |  35 ++++++--
 src/evolve.jl                                 |  11 ++-
 src/fv_rhs.jl                                 |  42 ++++++++++
 src/fvelement.jl                              |  36 ++++++++
 src/main.jl                                   |   5 +-
 src/mesh.jl                                   |  79 +++++++++++++-----
 src/parameters.jl                             |   8 ++
 src/utils.jl                                  |   2 +-
 .../fv_advection_sine/fv_advection_sine.toml  |  18 ++++
 .../refs/fv_advection_sine/output.h5          | Bin 0 -> 37736 bytes
 test/IntegrationTests/refs/testconfig.toml    |   3 +
 17 files changed, 270 insertions(+), 84 deletions(-)
 create mode 100644 src/fv_rhs.jl
 create mode 100644 src/fvelement.jl
 create mode 100644 test/IntegrationTests/refs/fv_advection_sine/fv_advection_sine.toml
 create mode 100644 test/IntegrationTests/refs/fv_advection_sine/output.h5

diff --git a/src/ScalarEq/callbacks.jl b/src/ScalarEq/callbacks.jl
index 347e8c73..439316b2 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 652c4a2f..155e43ab 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 63b0f03b..0649c2d4 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 c1edec36..26230f21 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 eb9cd047..52cfd3f9 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 16fe42fc..451e650d 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 d668df23..02e50bbe 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 9e7520d8..900a0f8f 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 00000000..18c6dcfc
--- /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 00000000..f6b3fb50
--- /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 8ef12534..9bec1992 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 030a7e4f..65326e6c 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 a7e6c333..9d05954f 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 8d9731f7..29648728 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 00000000..331a64e5
--- /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
GIT binary patch
literal 37736
zcmeI5cU+Wb*8T_V8WekpF=%YDVv9=9gNTaJ=!hj)gV;MNDk@Q8j4c-IU9oFo!GZ;m
zV;2z#iUkpI=xqj25s4D=8#vcJYIc3!&u7=&U-n)7!)u;)mvgS~x#xKphQ*E!oh?dM
zE2+u<SG>5Uu*S0BrTXs={2RSm_NZ-j+yyVppo^gU&E<uevBmsjRj}U@?dyuj4oB?o
z<J75B{!7)qQXKvu|K$;I?9@S<9|C8@p*S=MKEq~#*YaN%oUC^K6Z$Uywkal!e(>)<
z5#0ZxOW)<+MW}-+6snHCFaMgGDzDWT>Ejqwwn}P>o0yS3Z`c3V<0>lWVSr|6*UZec
zLr%qJA>;A<mu0QU*hBt5J59?{nwEbQD%j8zLBHna<^@+2M|*|LtV+olxC|TRu$g~y
z{rQtdXfz%;^8ZI($5pLVN@YBMj%!h@duEnqW*+s6ng(rVj70vN{=7cA$%?FT!I>Y{
zM^nf?y!bm2cvpQ%F7fus1+W=+y|4Q6rmx)_j5T%k=Kg6u(@gvS8E)LI%A4AGt9|PV
zru~X1JM8&0hxzY63wWz4C$zzNZ?&hzK^)U)Tu5GeK<7i3LRUf8K{r9SL3coRLH9uS
zK@UKGfgXY$5%uzg9)%u*9*3TQ`aw@ZPeD&Z&p`d5XQ2VmbI?F2i`tzRJrA|JfNiMV
zMQlUuE@2yLcNyDIyDQj++Fivq)b1L#p>{#ohT2`nHq<T{+fch3*oNALU>j<86WdU`
zP;5i(e#JJ_DonOrpgvF?)an+FLw%q+sMT#8hx$NuP^)kphx$NuP^&vQ4)uZRpjLNr
z9O?trL9OoLIMfHKgIe9kai|Ye2epd8ai|Ye2eo>D<4_-{4r&#N<4_-{4r=uf$Duw@
z9n>lc$Duw@9n|U(jzfK*I;d4NjzfK*I;hoS9EbWqbx=#49JhvQp>9xbs4p}OnhCXx
z!Szrr)D7wl^@WB(GohBTxE`v7x<S36zR)meCe$(x*F&{XH>fw%7a9i5gj&YqdZ-ra
z2K9#eLc^e$P|E~d57k25px#hlXc#mTYMF@Zp<1XL)Envx4TEMvEt7CPR10;3dP9Ao
zVbDydWiqaZYN2jWZ>TRc44MhGOu_X~Ez}L_4fTbFK{KJ2skk1hg}On#p}x>CXeQJ$
z4c9}pP&cSI)E62C&4gO&aXnNEb%T0CeW78{OsJ&+*F&{XH>fw%7a9i5gj%NKdZ-ra
z2K9#eLc^e$P|FNl57k25px#hlXc#mzgUmF4zVDH{8TCHupSaKYSl&~pCmFK!xA!@V
zOz%Se3NQXP9&+dM{_g|pBR+ra`@j$PWq$=%{&mOjuKO~t%5RTZ)sy#iU#2+yo97Di
zrC&#MZw*!F%R?TGQQlXbtN)pF@t;-a{<r(oH}4nSi>oN-$L-JS#lNaPS|hH1?fR(V
zp(3Py{^Q-`Z6Cgill|3TP?JyP#rxu<bMM}rG!b|tYBthLvnjv*vAcg;hz=q1uSMtH
zs=1hM_FXe=M8r<gS>?=}IN}IZ-My(~ShgSCy}0e#kF5h}`t!)p%(fS3K>N>+PVak}
znz1(e!GvqntFUh!pLxMlwqvJ+o~v)t;on=AKHW2n+CA(}4&k@R>HC<AmVLr$Y?6n4
z?5#U=wH`hCuKPVo4EdzlflK%4{_22X4O%~->olQZ(B4Q&s$mw<wOkY}j$1Hr#>_`_
z*CWHJYE(3xYTtcfZgm~4SUmV^x9d9U`Y?Lf*HdH2IjzOUnbl$`_1Vih4+CTAS^uI(
zon7MSPQOmQALPYRn-ZHhJ8X)lb$xfOn%6vmmYR(SoP8~U+IwWysoFo04)0Ix<B*g{
zL#Gt_HfeGaxmErqPM?!R+{C|(oo6yxtjX9F@+z4++`pDJXGRLGIr4c%$BYz;eQEx2
ziD9W^dAH`afZM55?wjZBziyvKx4m19`gu<pb?rN)?&aco+B|hp%b!N+X;5zNv&;*6
zdf2MHds1Zs{gilp>d-L;+U(QD)#;>x8U{9YJ)CEtITsps8*86VAKmKy%k1yd$-UI$
zpop{Sba86@E{`Ya^uT8N3#VEc)VOJ>FI{?M&<v{`yEabCpv@(&*ZysNhH;#0C(rcD
zpxJvT{9Z9MgGvlMmm3_DL5C(UAF<Z>zHji0yo@O?Gw9)}&6QmLc&%!`I!|4v+EMMQ
zepJ6|9FZ888s{{O+XC~Li+QPes(Gt8JVIOsBTg!ADvsfZtBP|u#61RiP<gROo=PKc
z=aEO1SC!`@$h-2va`<8ce4>1#e58Dp2A>Ut?;PPn<;&CXX#{+`3_ebPua(au;QLf?
z5DYHdz)3G~L*U34TvY{U1HoNWa5xZL&H|^g;I<Js4hPqY^W)&Yt|%W6{@GW~<B#P!
zz9sGO6ls^8q#v#;{j!scgDquTd{D;8HZpGhTIRvyWL^yZIYH*lcO(vcMfm3!i4!l8
zxbYo{BcGJG@-pF{J4@WTzvO|zKR*@zIb8VX*1|u7e+K{DPWb0T;sbUO{&}MC&)}aw
z5+5=6=hebL>%@0lMtsO+gntebpYm72KW`EK8T@mU_?(vr|GY`~=LW(*e<z%9DdC2v
z3P*fN_-Am&;Eur|Zxa5E<1FBx!7YPh2G=}D_-F8MQt)p&PhF?lQSE|%_K>>6O{5O-
zL*bt<W1Og44E`B)jln;Ie_o0>fPWq${PR}Ctvlj)3vn%rIHT?|>M*ZGUOEc@Tn2f&
zApG-G<aIdmd<=QFg%6Z3luz2iH+SKqEci+<{PQ6At~-3V8omsGPs8EcMes4!HNF9#
zN5c1da1bK=^AvD05ZpL|qhG+41vt}!yB6T^TW~oSoTh=>R^Yf!&4P8EYxFBv=lQpw
zf_0zeZ9o{Ep1fo2Q|k=QO4_jSCv65-shR1h9hbrJ7xkxYzt7+)S1V8ddQ}FWG+d0H
zyDfvC?{O=8drt;Ccd22p-=D!b*!IV951gNh>#fjUF4_-3f3-b^J1zOh*v{GN1D)g2
z*?MTPlEeJc`O~IjrX8J^&ONS7pV+NiIyWh5<<+)KIzN37XqX>r;2QG|9BsYWz=Qnu
z&mKeuUU<IOk-SVjhZr1O8*J0_fEn5W`p$Z;(bJ=pQ(783mWW$B=BG4X^X0)y9qOgA
z-M%eFQv*}^#KNFyQ=C#+GrYo$w3rmmzU<>NYC;O%Y4u5y7g@<%?l<$Rhdh$`Wb`Sg
zZ8=G7duZe8DEA~jz4GY;PeUSmmbqVKjcX!DrW)q8zLUUK%WHM|&OU*gjI+No`%pah
zd;H_?T`R_O_gY0p{xLI-pWLYW<oLr_&YOFvn3qi~H|aI~N6&RJyw9zs-^C0aUq3&-
z`}oc}&I|rFtZu8voE`PE<(`N~e9e7v<;63ic>L&D&Lzq{<lO3AvIgycz||do@I7l2
z!A-26r4PP-k9U@f-0j)tE+3rT%&AvYIM<(kBd6lf+q`P4OZg&kVSKtq<*LUohjP*v
zBgbtGyumYDU7u}tB8V#-pk5XBUg4;KhKpCOxyW~xPhRY*Vcu4HV57*}r@8c_{!!oU
zIK}}#?D%Yc<bJN;(mZlv)CTzH=YT`s7HrcZIZoM81?SU3xt;>0J-Q|BQ%UKMCQAPl
zCF7xkGCo==<E5H1e%dGVp{+7MY9RBaN-}@SmUz%8i4R?rcu}Urk8&iQbWP$*dnDe}
zK;lm?B_EV7`Jt~QU(`qPM-3#OR73Jh5t486ko?mL@dLFHf6xi>3t5VP=(hNYvcz8m
zzY+XL@FT4gf6`atSF#fS(jf6O!QZq<{7yr}{{$Wgd=Pjc@Iy6(Cu%Hwk(2O7HH1Gp
zB|Oq|;ghNguaqGC(p2G@{Dg0MD!h}E@Ne2af#awbg!*A>Pqlv%{i*)dc+~jRc-8pT
ze4YvarunM*FGf5<g@04LRQyytReV*vRs2;xz(0Y10{;a5$p-m7iG0sO{*@n;KN`U=
z%0J3a8St0#TO;^yH~e@N{(K3)X2HK#;ODXM_e=P_8vNe~Jotc*pTSFA@Us#;odjQr
z;B7kivjmS1z~@cyx&iza7yii?eCL7pg7w48$JSrAc?If2j_}V5P(S=pPku#xDUNzG
z9`)w|>d}7GrzNOYpQ3*4K|R}o`c@zH4*WCtXYkM9pTR$ae+K^y{u%r;^1<MrkuOI6
z82mHx%g8r_e+K^y{u%r;_-F9X;NP^~sP#v!NAM@ZuMGb(_&2R@;Gf}t(|V}ZN3~w6
z^%FcX_~JfTZ>_NYet6#d_H&UBTW?3;t>=O+KbQ`_xqsSeJcKnKRnN!X|8rOC7Sj3u
zt>><P@AEwQ^Gt3kvwipT?0?mB*fxmsU;8=ioBnh6YAnGr{{QCwX%ToRYuvAV`RC`h
zj<#~C`uV|!dT0vehZlb*0`GdR@74V6li&Bxzuxe5;Rz$Aj<zZ!{k^f@tYsBdE#B?-
z{T<{w<#*ZA(UR?b%SPY6OW&TH>DAli4%w~PyKHyrZSowuG@|5`TU0%~NV1pVR~l0C
zbN@e{gixi<&*F+by-q`hl(=0c{VJK0*^1?fmuT|pqDNmnJWtUP;mfyeJxjq8ZM+N3
z^`nSI-)A-TKTJnHn}2Qd<ee09O;aLp!y;-e%1dNFACvR=j9kZI(hiT6b~(3wK|kDI
z`ejQQ2RqBSxT=hkXUn)5^WZ@;FJ2_`<bW1ys+Y3A&j%ZgDYWxY1RtZVC&pHc<o#1y
z_%)vSkZrrasDC0NimT6j8hNcnG(YtZ>eOP%W3E3S_oi-wjyqdxMt@o%hL5~x7JBwn
z3|EQ|zOuzRmOa})o_ioWmeT{8ZN0uSjwkMT-Yntscn&Nz_j-I_Jdf4u#vIiqaQl)T
z>svlf-~r2Bot?)fa^|wQ*2SJCviV2Hx;C1g#8WbxH$I-7#M$m!Ef-BqX3M1|POR4_
zv(5EU<9{8R!nOLk*8L<rg~xi<9W%);l>^6C>hE_TmGkzw`bU>e<9nG)o;{qD#>vI5
zmf04P#(S-_7b?}!bKs)w!`jT$bJ4bQe}5jN=cg}AckEcfz%#=4udguBz*|QDQs&#O
z1}^gP=ToXYG;pf-Y@5~<(%IGh?8ak`=^Wpny2m_^bY7!vIDPvr<Mx$j2J_%_j;K{Q
zYj0XQcOK>Jyr6gnH@;kV-xjM3_8hX}^TZYzY@OL8@mu>0Zr^Opj$++2xODdV=!eGV
zO-kEge-oTH3)eM4JF4CL=qDEauEaRvFfKJtHSWTghnklM=9z_g7e^e<BQ6IJr&PpE
z#ZkpI0C66LxPOW~v_@XaAWvDyn;Y_|L0<L9^Eb%5EqpK<zED1y0N-T7M|0q-D)8Ac
z_^vB_*bu(VflsaATRZspAbfoqK0gNEUjzrs!G$Y0aS;A_E;y<Mu3W)cb#ONl9F`RR
z*#n%Wg4;XbI0;<)gY!&qZ-qLr2zB8k>cm;pjd0YF@u(|#!aof}-6@MY1pW#96ZkjP
zG4M~|pTIwXe**sm{t5gO_$TmB;Ge)h`AQxL{1f;mgXE3INgl~x_@|MQXF}e|S$sh0
z;tPUL2)-fsh)RpEsEzoHti*RTT6{?GCFKhLbVvB7MDa207GG0<_?!xh?<rF_plQMd
zZ4&;ekZ?obh=vOP6f2yOw{S-~;gIeKmt-&e)91o1%@B^Mk#J343g=W(_&2Ey_N((e
z)OBdb0qx#HKj5EMVH|2)E*NLB)Gac{Jk-3VW1gscggQv5i?k1MN<rKr5J%1ag1Snm
zvor>AuZcXUyr53gbK##*#|d@alxLOqcEZ2O7s@BfH+uL;`AYdr`OZ=Jr_Y6d$`<~q
zhVV}${F{8Oe11syCl>zyg`W#-(^1}I{CC{9`LvVSs&``QzT@@%@W1N5%?@$?54vw_
z^Q8pK`2U;lW2brJrp&ni!~K{kL?2%Kod~?^e$3MD?UVn7_hV+w!|=<E#xCFVajW9$
zQOfK4vDbA#>Y*CfKmNIEyKXYm|Bm{qL6!UW*H=fx`9G+>Zg>sBo4%ZuiH87hBVyc2
zrV#_RE9$LIrSG+_)fzufrS|i-PA%Cxjp`+N*F1VKjRtwt>{Z@ePcfOl-2be%o_eqP
zsMWpgdMdYl>Dn5}dRpGjFZ@YE17&&^o7HKYfl9XCc)84816A-|x3f!>fqWzOor<?e
zr{6byIry<nI>l~3vCDc`I++K~?tOcKaetTRhQK}PG$gD4;U}C<t3H{!WaQm+%C@y>
zQ7ti@MlKpyvE<Wqs<CZz|CYa{leVc<!unUQx7BfVzPetumy7mQf2#ji7|#)mPmQ-Y
z#-EP)EXDlFV7^Y6e>CDT0P(4fc-bI+a}ZBg#J6wPXBovCr_$pCTVCldrqIHk*_P+~
zq|k`Uo()gMCew@WAL;IoPo~ZDy9{ahEQvz$){S^PEr||xaQ{l5l}L@NJ7-LEOC;x*
zSr=y~Cs0((A^&=Vjh}D4`)Qv=q48vU?4xDJZR6>Snf+sVUmVSwUcLL}GI6x($8t6i
z?y<B}XXai0ehl3xRchz*FJkE2_+)F>)jHbJ*(RjN!^c!@YCWHfCDD|)Xh1+@t49=4
z^M^AIQ4h)cdCgeYIgzxbdanlKsz0ES4X4ywbMihlnsljH#jo#C%MFhf*D?Mc-~5^%
zXZ=1YoJL<yEp<Ql7WJ51+9Gv97}<TdqC-seP|7}<bf9a!8#J$L)4KB-2T|6^vrRPC
zSIF{)<*#j<UZm#Z`q%6BO(4a0^2%`UcABoO_<7IcFON}R;LfVgo%U1ZDt-g+_Sqo6
z2cy|wn^?JEoA=3a4v_OXLayif(jI>#?Q^E|$6nGuXUTZDgN%>0GF~1h<L3o3AO1q-
z$4z9uJX+?@hzBD+94zr-#E%hA-X`(o0Est;Nc?%X<b&5ret3c8i)|!-Y$y3-d&w{B
zB;Pzx^3P4h4}4hs!SD-j68~_v_=$&$zj%T8jdkKbb{0SKNbx7n6Th;Z_?PR8pLvM*
zn`6cA+(-P+WrYU@|NK;V;g5tLt}OiXap8+Kn+y2oMB$%{3jb{8SinEK3je%L_~)0x
zKaUmu`7`04`wIW2?a?@%gY(zmdKa{pC;W3^^cRo*QIGg2##agBMg3ybGt>OkeAWDQ
zhzIH;*G0T+5kD`)b2#Gr4dNYx_%B92Mj=1HAzxdOzfj@dlwXza0my$e{LlgZXa~Qv
zfq#^r*1%uNZ_0n*pSQ!G)rEgf75+I4eqJN|vybr4W8r^m@X!`~ECDYQ!H*Moas*$m
zz}q75mjfOvfKMCnst3Pe;5h_*Hxd5jb4Xr$0-}U}+Ku`E{wWgmqao@^Y1Ee|s5f&_
zf1aTpwMTtY^=b_27x*{TH&yRc{ZsXD1?nUCH`UK7sHdAyU;R;UgHeA~Jy!Mk2h{6U
zsNbrd6Y9IF_p1J@^+2r;!B{WAKY@P&{{;Ri8tV`Ih<chy$EJ+heD3RbZt`8nW2tB3
zcyBNF3m?^s<J;aR&)KYqW&7;GmP@l@xJ8{atE+qy!#zq2Uy*ZC$IAvSKH*S7$2KAL
zS~iM~<}6E}$setF#MLL|h7SHRiqFR#I@BjFlAE6}`?&7H2fW+r*|lQzBDnj5IF0qC
zdmQd*yXWq}yL^6Zb`9$%;e4UB<$xJK+~&3w)6d40yTwJ#dbS!wzp`_!g`+?16T(`r
z;BNDWUuR8;+gI}3uJV_TsoxZub%`ss$SqaO>pZ`_Y<4BD=vh8p@~YS65`JvkujAD9
z<qosYv?2lK7CSlj&hLS1%PnFTQSy-eR9?=bU2+}8OFL9j+9hl0hq8(n^h+IO9LCR&
z6^x5sY6`|l&1Kx=FY}-?GA{~|d6I+7o3=?DsDZ?VPD`AqqQs5TC63fw;!2Ao&eTQX
zPK_lG)L!yJp^_)sEqNp4k&a7VDJ$MVGbt*b<}ZElIC)e8MFd-&u4w%IWQWtk7wg;;
zDY@$G2Q9M`$%hRm^QI<|$4!fX2I)yucgB;ldxs~}v{M784!Li9pW?l&seR`Z3L33x
zY<Dz;4wfmQ?@}d|+{zAH+;w&;)k@rTzExx@joZC8{r8q>R5rq+{F<NAXvJrZ@1{OV
zqeb7?egC|no<a`CmEJ#1PX}YJTD|bsQ@3joMI&D6DaiO4PrEh-YE#$oRXsNYd01Ff
z`*D|ntXHhI*4{SI(wyzZ=H(h_SZ%LQD%47+m;FB4r|FnZKR^4)+y2{h>YVc3oV(N0
zsqIOJP`}0Lq)YKzr(19QJqwFv`+M(5r^Lq>mRsz7t!lqI&j#12c2v8npONTyImV&J
zrN(K4aR*}_Jut5$n5UYzii0=eQUP(Ah`7~79H$_zhY{!Nh<gt5;Euf1MxJIMZ#$63
zG~~4l@*Ii0FMtoS;S2CjjfH=ckHX+9JNT>_eAf*=TmoNy0iP=0R)UZH;p;~5x$=Ey
zaBve`_=6M0O$%`J3S4yqXPv=aHaKh!E=z+`1Gqg1j_ZT#<>34xxL<=hFduc{Q`Cw5
zs2i!MBWRb=55JIp8RIa1eoX4p3mGR{qi&r>9XpM>b^~>;6YAb})Ir3B15hWcp>94!
z9c_cUx)gP`JL+z8)M4bs`1=S_r&Zl<iaM_9I`YiOJHrRuReZq)@d@t~-|#H)5!V-A
z@pSPSgMSVdA2NK&;GaJe-|__UF_#x#^GfkKFA@GZRQTt^!UcaVoESehU%)@x3P(I!
z_~$Q#e_kcr@mS%I^MrpsEc|mf;h)<H#~d&Gv!(FQslvZW)&4k~r>^@7?WlHDKdN6f
zj#e0#C&t+h<GzV`XfZEe%u~%<#bGt#QW<fYjJVZD9NmR~J|g^^;+`%1Gx%rJX;a=*
z9`(qpBl4{B{-f~ESQl75xBi>Y%X}jb4b<lcK0KfOn};YLcDxaRcRio=?fb?_#+T|j
zsP^sM|M1qE&%2s=;ma+Jt-ZPb3&-B3_Lbk&K;GZ)i#raKPX0Td&-ysaZ2$ewXSIm)
zf6(*UHvLTgd()TGyisJ0N0l%CT!$lu%cURcuqlKeUi=Rcc++1kTwZ?(!JGSAhXGzw
zJ6Ys3<m)XV`(oL>W#x81*~W2pNw--}p>bSPyZAw=!SP%@td8@J<al1QIB(MLZV4Re
zTA^C0>;#Ud^jm&261nz^iyKm&Cvxik<vFX|l6dDI<%(@jN@8v7h0(7DB(p=^RlU~*
zC$n|fxs^k{Na0Z)W9E0<mckRqw5;38Je5PsAGN6Mn#%6eu1qpB{@z~U%Fk+Y#Wap}
zSg@t%m^8lK_+ajy(`kHW&7`6y3+cJVwas;nPsoj*tNYRY?ovGu4?lDK*I+&G^IO@s
ze_;b33IAo*n6?J)u<Wa>7h{a?EsL-3^!r)^@4Xtdd&)TjoBK4&?WHqt-RF(d*5(;F
zGbiV2+bZci-!Oc0+s5hqecKrgSG7y$8tsyw*mX_k<>3ozI%w0m{^YrPHaWfCR>!??
z{ybctj`mdhsz24gJ;q~?@g2r^_hS5NK5Bl=F<*boU&SLH@lo+o@e4ydCnLUN5bu_V
zza#R|2l?5Hd>N3xI>=`t@~iT#@?Q;pQ2tPUnGXLbKlQVCwXm_J@i`K09C~G13?G}Y
zF6WmB9XH+b^Mc>n=s2nH7tYHXKjy_D5AQq2KjKT%YG@a&i{eSPH*I=4Jmk<VfgEFa
z!0Sd|Iq1GRg2#1SI-y;w``kCkxzLpGyX>&_RIcO1JM3|;+=SnXgmZf57V-X@Zt?S~
zz9X(Q3gZ(tYtQfWxXFeW)4poDI+z!Xj11g);2Q7U=u*)w;4-&uckFJ-kPAHh{HRi)
zN6+!$6&;#By?TnrecZ@*T*OiC-j>|DhVJ7=Qy)zKG;BS7A98O%_>_Wesv*azkDN~%
z<a%<I_Nb|}PurwFnj`&Fyo`sA$oMEr#!DDKVLoIh^P^reUm7LzCvS-d1xb8pl*Eht
zC4Mwd;z_X*Uuq)pCSQp^*-1XAx8#R<2>&!i@<-2ve=02e6Y@=~gnzm({F6@nK~~}y
zN)-RlHt`en5PuQ)r)J_m0{=8g{7E6=S8@~nNhADIwD_Asg@5uB|C6=wK>dV&auNP%
zx9~#;g(tETzNo$MM!SW7(g^=lTlgov@Je37FL?^jR9*O|rNTee75+`zpWyghoUg7|
z?Wy)vf2x1Ncu=1R{1f;m@K4~Mz(0Y10{;a53H%fIC-6_;pTIwXe**sm{%H^LrAPiy
z&nW@<eTaOoMgFlKP%-!e{1es>(!x)D;4kI3qwt^dV;ubXBmBzn?>P8582(-bzbpSM
z9uyz$!at=7|MUbr{Q|yX!J8-e+X5b~!RJWux*hxu0ne?#cOl^)&&}uAW(xn@81-QQ
z>ctk+k71}MO;BIZAEST9co^g3Ow_An)Gy43F+X-fy&H-82mTrS^Jvt|0MyU<sHcc8
zBi;=D8TsJ8sL$`Ze{&mVdXrG_;(g!0srS9_|9(t6PDcLUaX)4~QfB+_zaMi(od1LF
z$AVs)`px&n)2xhN&yRn~mw&z=bH-PKf4Co0k@}$jQv}{sSB*bMd+Xo(uB)A$9JQMK
z$6@C90}chfslWAL{Jyw?f9eb5jAzL|rPzG^c}MJchZN552mhpLD`P?Tov!w^OPbTn
zzK^4W@sCCr&nsozJ-oj2R5@Bwqk5^bve}c)y^UvS)c90Xl&ba$8~+Y~hsMgdi6i*V
zRdv6sw3hMXV~S*T>ic`*=R8}+^PjcCs2Kh|&7q<f%)>VLq>^sZ!yxXGN?N-*JszZ{
z(Ckxm{Mq6ZO7-iR@<|=z_Le7}qXLrYBaiUax}M48wW@XM#>Yvt=*G8UyT>I_OpVxj
z<&B?TTXrGQ<%LHgMcHd-41bwGvucKBOz=pcj~yfL)O!|BtCDvd-#a#*9GgaN_#--w
zdUpTb!ciMXC7qU*UuXQjzn$af)zUPGrItBGYHwW?L(}uRxy{ef(Xf5<f^QDgQB>a^
z{<UU4rYHN?kN>q*G<EE^bz!(6ih}o>Zy&njA?;lGWA;I9Buz`HG%Mv*1eL!S;n!`~
zeX3QsaA@%^_b7eC$?5%)jNe;0?3c8$^TX+4uSZFqHEvVK;7VIB`GrxZO^^B+)`Zf;
zPsTRea^?oj3fi)waAXkm`K5?%efAX^@BCf;pfZ<e%b^un#p<4?+IhPcu9@mjIsQE|
z3$H&xwL<TeSvm0#jgN0TbotLas9DfM>#i#nQc+QMmHk{z&f{lt9oLX{c!RXdqop6#
zO20f^#=&D{Tx>7nWIGu*pO<;?WtkUao{V`j;=t|_7xs}jagfA~?IezzA#vqMi8CYa
zj6862$qOS-jJz@O$Xz9`yjt?ipGe*rK4AEQmy1t0R(!)t#78_&e8raHGxieSak}`B
z=ZG)4kNA}1#kX8Se9R%<U-<TWp2{J=eeSI9l*U26zYO?#RT>Ab^3I$PnZ`vo&VOlN
zL(gZ^E|^~&tmkc^huX(%&~sw_f!^zX)$@??PZvg+8MyWCyy5$s8o!UqV|(Q+X9Me>
zZ|`vCJL7Ytr;}Sh-D=!^Y+<(PjDd&doE&g4%)rA|@A9h?XW+z5i`x`>V&HX}GUx4I
z8u+U^H6F(rpEoI0`_*|vah(g=NkzMVJdu7<(eD6^;|nz|j58kN-i&#i#=Ht+p8lBk
zRm34!#RYK!{~U`rIw7v#A<lZl{X68LH1eYIr1GZnsPd}vJX`qZ8t}my_@Wzp(o6U^
z`AGSy1$;J5_-9w)pTR#rg-<)dw;}NHZuoizd|pxb=i0(QXA1uu1Wqo48wYUY2(H|~
znJ2jG2M#BK%WdGaA-L6?E#RN6!L<iC9{}#-u1OsjfVxltb>b=NMs?H?@J}OAXL_LS
zxT6k%f3icJ0{;|<I(89tO^-VFDe7KX)Iso1;GeFdZnj4q)uXOHK%GV0$sKhV{L>E9
z=@O{hs*Zn!x~}Scb<}+itOF5P7Zzcih{n3H0PDyctScq4&Zu=K8S9V-)}`)Pr=qZK
zX|Rr67GKj<@i|=--;<wkKpNqKP6{WqNw}d{;fMwZSJYBCBXCFHkTk+Sfl~sv<SHB!
zxF&E;;NGNazd8^6o7z$Bf`3!LY8>F724S4wpTIwXe**sm{t5gO_@`jR$-Ynb&YB|E
z_dKeg>2G|8V4+F;qhf7M_kx!%=Cx{~@hW)v-0DnAO-18NOU=|R#@|R=R<PG|$Er=5
zIt4H5_ww@7<bPL!g~n`f&oi3*2PYPqCZQHFn*2Le3(YrMf?sIz-vU}_PAmyoZG7?J
z-(MgC|J(P6zw3K!?C>W3-~ar<1@FTD+Rq<4J382F)SoiU&+~|9U;7#VsNP#&@{#{X
z_DmmQ$loyi###P`>Bo8b8>a67$=@*j0$Tot>C+DR8>S!i<!`ic!9$|_jkXbg`X>Ls
zfcgFp58=b-r~VHo_dQJQUC)U!XT9CHdJ^)!K0oz#>uM{PzArWMp`VTM@4)SyRk!uF
z+th#8Add%m;q+>2Pq&eS?oy!_J?rg1d5=Cj>Dp$DMFi2HTS4wF59r)5t@pYkk@Rb$
zGBKN8J|vgQzDIg=ctn-G>usL<Lo{7{;Fb02{9|h1+i}=+PaSPPejt5ur5Kv%9+kN7
zPz*($@0M4}K9<IRR;2g#@K_3I_2raqu5na0?CijohB)&7sn3t=rW>CdRHzeOGB=)@
zT#mkFKPQ2n7JJgI!S4yws=`kHL&oO~<HnXR=WqNwS{7%U{c1ZRiAEh+@$$2KNwl#+
z!-hjTCeycTsxOW_m`rCbH4R%-E`@%rHPWfSdkTGiVCYNR+bMLk`{}RiG)kp2%bHFq
zJTH}Et~|E78J<d06WU*$QZtP<^(Z$oWo#Pxf8XU`9p5xc>2*A<`ja#|Gw9T(R~qZ7
z_57QL3M2Jo-Q>5ug*F*K_g<z(*TvWMv_@;|wJl3e9qz2$XZx{%*39t<9^B4AI>)Me
za|alxNRfm4E4Ulz#G$j#b}TZ`p!zQc?$}`5KlI)h%iRWgdGJ+#!_NjPH)7@|tq;B4
zE{x*?aK5@;wWr!w{n?}cbr{c0jISBStH$q-`J`fg=P=*xnEwpKBLVR_hj^Vt{F)%1
zcM)GkyoVwF+mMg($WK}1YZCHjgM4Nqzs-^Fb;y4p{7?@5xCy@~|2%}BCc$4m@LMDJ
z?=t*Y9sV2+zutp?FT>BB;qM0UdvW+b4?GM2A7<d?2k>J&|8|uFztT<HofYHnhfw`S
zO;6dH2UF;09k;Hke~p6YjO*U7+hy9gpy|Dwu@`9ntqv>hmpMn-bgZd;r&E+y_-3O@
zGmg@%Q-?=R?!J$@H})CRuE%=(9CD>?ec9$Qa-6f}eAdhLyhYmMmC`<6mi{<M`sdm*
z9uAiAaeWyt&y(>p_-F9XxiViaEc53&5)bYm@!|6lFAkIVG5F^y5??+q@n-PP;Ge-i
zKahMe^2gwxEhWEvR`ShjCI8$){J>e_4{jrVVfcqF#ZNp|{KcEZZ#-1|$Nu6+zAgUb
z#^P7@6aO;&%<wnE@4QX?&w;`N8-x!&C%o{N!auhZo;XbS;tj$ZI|zTgPWb0o;h&2N
z|9n#TXYkKOg>TLk{<(tiZ`ua`4E`DXGx%rl&)}cIKZAcped6XAZySt10P{)5{LW*(
zKA69Xhl-Dimx^B_#Pbf~3;r4WGx%rl&*0ybKk(1spTR$ae+K^y{y9|m=LxcYF!*Qi
z&pU*FZY=!sRrnG7^GM;J!9QPtpS!@{4dM3^@V{nT!FtGpgnu>{{u%r;_-F9XpMy8>
z&o{s$_~)O&>r3!!UZjA3o(H}Mf%k%y+RG=P3hKd7)Q6|27s;p}8&FTYQC}{i-dslg
zu|hq%g8Ebo_3B^uQxlh;<Xy&p)=y3T(feX8zOwY+e_uQrasF%H7pr%iAM`^6K1AR{
z1U^LIFBJj1A%B$k-lB<t%kLSG=TX_fw$6E_M(67JeCtW~qNDWuMert@*wcDmTDaDg
z8Qyw6+-A<vJ;U^zVo@q(N^?DD?|*KVl$*xEjmx$#dCB;9daipfPMDv@w$t{_oM)fL
z=MP$ROM8{dtGZ3lH8`5ebMGBL-`)6o0FINoRUBL<mFxT4>>qtPg_H9d#(EA;;UB^W
z2JOpB=9nFmUcT6z%;rN{`g*o9{{5aSolbSSk;E<wwz&-(oW#`*lnqTxOXS%7hkC|L
zPvrQGCCo<TCGheP&*}5$C2+$59fF+-C2+lu=537gis#Pd`j~Zl701?pY}veZN*v$+
zF=2LQVl2OMuiB`;@%w&TPHuI;a5aW^X2yPdxp@p1%`^WzYqO4HyVdW!#$3k^oy-F~
zZbtLTiNl^0U-^h9X1C8Np^f6PX#=%~N<QR3&xdvDTzJ6M`g?ZWJvD+8f4N+LNUi(q
zWGH`i!Ogqu6u8ybY1$obW6?7=yJ|Q`T-?_7z}Z{e<;!jpwhsv7V4qT!R_QnS<aZwj
z46Aj6?OkWI^Kc5{(VlthL*`uJv~RW#3)^>*r&RAWJ!v!ZjFlJr_J}*fMZE5`KUMbx
z-(LN#wnXSbKDum7(H39r;7vL6t_77|$TdW1mF%bEavnXOTW}qXly+#Cv`ewl4-Jui
z$wtOOhh<!JS;k2YGH$vf^PrY8FX|!lq#iPFij_D}C5a39OPpw=#Er&E9O;C_mFy+X
zgt${-$piUFUZ}a`iC#+H2zeyrl^i6`R9EtD{65VBAJ9_q1?j{m6fM4?D&iv=B)+0!
z;xoz{T;My379WzA_>$tqr&L>fOP$2W<SD);Z}B<7_XG~;E8&7dgcEux+)%1;MBs|R
z8MPMf2pkf)qzb|*xd^wkbn4GlYS|em>(K01$G<j^L#dl(FZ44$cWR-{?lSPTs{PM#
z-cPvhTePFvRsCEP{;3Vdp~j`g3H}NE6Zj|aPhSfErZ}uYT)YscuMjuzPvD=xKj{#6
zjqpz=ke4CIlggV8d5l3`e?p!cAn(cto#2aR!aqd_|C9?Kfq#?Fl<$-e!9Rh20{;a5
z3H%fIC-6_;pTIwXe**sm{t5gO_$Tmh{9fP!{t5gO_$Lo=TS)k)P;h+^oG%3TU!V?b
zL|r(II*}v%^ElLzQK&2Ns57ppJMB=1PM|IYp-y!|-2(p%{u%r;_-F9X;Ge-igMSA9
z4F1^;bykPEYmPd+4Rsm(Gx#^vaaGsBKO^r9{u%r;_-F9X@DYQ52LJrK@Xt|Lhvr~i
zip4rr3+q-#tYfpVt}Vwp7l?Jw2J4_&7u7nMgLM=9voF@w<ydFiV%@!jb+|p&<w{tm
ZU9oO26OI}DGdO4PZ&I}%{4@CH{|EHSeNg}a

literal 0
HcmV?d00001

diff --git a/test/IntegrationTests/refs/testconfig.toml b/test/IntegrationTests/refs/testconfig.toml
index ee94aa62..35fc1305 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" ]
 
-- 
GitLab