diff --git a/WENO.F90 b/WENO.F90
index 200d3a935ea354321c32066657855d5342a5d19a..712b8b773cbce21fb5f00be9b3fc365330225ebb 100644
--- a/WENO.F90
+++ b/WENO.F90
@@ -23,34 +23,34 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, &
   implicit none
 
   integer, intent(in) :: nvars, ndim, nx, ny, nz
-  real, intent(in) :: states(nvars, -2:nx+3, -2:ny+3, -2:nz+3)
-  real, intent(in) :: fluxes(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
-  real, intent(in) :: metrics(ndim, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
-  real, intent(in) :: metric_jacobians(-2:nx+3, -2:ny+3, -2:nz+3)
-  real, intent(out) :: flux_derivatives(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
-
-  real flux_derivatives_generalized(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
-  real weno_flux_tmp(nvars, 0:nx+1, 0:ny+1, 0:nz+1)
-
-  real generalized_fluxes(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
-  real generalized_states_frozen(nvars, -2:3)
-  real generalized_fluxes_frozen(nvars, ndim, -2:3)
-  real metric_solve_tmp(ndim)
-  real R(nvars, nvars), R_inv(nvars, nvars), lambda_roe(nvars)
-  real lambda_pointwise(nvars, -2:3)
-  real wavespeeds(nvars)
+  real*8, intent(in) :: states(nvars, -2:nx+3, -2:ny+3, -2:nz+3)
+  real*8, intent(in) :: fluxes(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
+  real*8, intent(in) :: metrics(ndim, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
+  real*8, intent(in) :: metric_jacobians(-2:nx+3, -2:ny+3, -2:nz+3)
+  real*8, intent(out) :: flux_derivatives(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
+
+  real*8 flux_derivatives_generalized(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
+  real*8 weno_flux_tmp(nvars, 0:nx+1, 0:ny+1, 0:nz+1)
+
+  real*8 generalized_fluxes(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
+  real*8 generalized_states_frozen(nvars, -2:3)
+  real*8 generalized_fluxes_frozen(nvars, ndim, -2:3)
+  real*8 metric_solve_tmp(ndim)
+  real*8 R(nvars, nvars), R_inv(nvars, nvars), lambda_roe(nvars)
+  real*8 lambda_pointwise(nvars, -2:3)
+  real*8 wavespeeds(nvars)
   integer i, j, k
-  real delta_xi, delta_eta, delta_zeta
-  real characteristic_fluxes_pos(nvars, -2:3)
-  real characteristic_fluxes_neg(nvars, -2:3)
-  real metrics_frozen(ndim, ndim)
-  real combined_frozen_metrics(ndim)
+  real*8 delta_xi, delta_eta, delta_zeta
+  real*8 characteristic_fluxes_pos(nvars, -2:3)
+  real*8 characteristic_fluxes_neg(nvars, -2:3)
+  real*8 metrics_frozen(ndim, ndim)
+  real*8 combined_frozen_metrics(ndim)
   integer v, d
 
   ! grid spacing in generalized coordinates
-  delta_xi = 1.0
-  delta_eta = 1.0
-  delta_zeta = 1.0
+  delta_xi = 1.0d0
+  delta_eta = 1.0d0
+  delta_zeta = 1.0d0
 
   !$loopy begin tagged: to_generalized
   !call convert_to_generalized(nvars, ndim, nx, ny, nz, &
@@ -61,7 +61,7 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, &
     do j=-2,ny+3
       do i=-2,nx+3
         do v=1,nvars
-          call mult_mat_vec(ndim, ndim, 1.0/metric_jacobians(i,j,k), metrics(:,:,i,j,k), &
+          call mult_mat_vec(ndim, ndim, 1.0d0/metric_jacobians(i,j,k), metrics(:,:,i,j,k), &
             fluxes(v,:,i,j,k), generalized_fluxes(v,:,i,j,k))
         end do
       end do
@@ -250,21 +250,21 @@ subroutine pointwise_eigenvalues(nvars, d, states, lambda_pointwise)
 
   integer, intent(in) :: nvars
   integer, intent(in) :: d
-  real, intent(in) :: states(nvars, -2:3)
-  real, intent(out) :: lambda_pointwise(nvars, -2:3)
+  real*8, intent(in) :: states(nvars, -2:3)
+  real*8, intent(out) :: lambda_pointwise(nvars, -2:3)
 
-  real u(3), c, p, rho, ke
+  real*8 u(3), c, p, rho, ke
   integer v, k, i
 
   do k=-2,3
     rho = states(1,k)
-    ke = 0.0
+    ke = 0.0d0
     do i=1,3
       u(i) = states(i+1,k)/rho
       ke = ke + u(i)**2
     end do
-    p = (1.4 - 1)*(states(nvars,k) - 0.5*rho*ke)
-    c = sqrt(1.4*p/rho)
+    p = (1.4d0 - 1)*(states(nvars,k) - 0.5d0*rho*ke)
+    c = sqrt(1.4d0*p/rho)
 
     do v=1,nvars-2
       lambda_pointwise(v,k) = u(d)
@@ -280,21 +280,21 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam
   integer, intent(in) :: nvars
   integer, intent(in) :: ndim
   integer, intent(in) :: d
-  real, intent(in) :: states(nvars, 2)
-  real, intent(in) :: metrics_frozen(ndim, ndim)
-  real, intent(out) :: R(nvars, nvars)
-  real, intent(out) :: R_inv(nvars, nvars)
-  real, intent(out) :: lambda_roe(nvars)
+  real*8, intent(in) :: states(nvars, 2)
+  real*8, intent(in) :: metrics_frozen(ndim, ndim)
+  real*8, intent(out) :: R(nvars, nvars)
+  real*8, intent(out) :: R_inv(nvars, nvars)
+  real*8, intent(out) :: lambda_roe(nvars)
 
   integer ik, il, im
-  real metric_norm(ndim)
-  real p, ke
-  real u_orig(ndim, 2), H_orig(2)
-  real sqrt_rho(2), sqrt_rho_sum
-  real u(ndim), rho, c, H, q
-  real u_tilde(ndim), k(ndim), l(ndim), m(ndim)
+  real*8 metric_norm(ndim)
+  real*8 p, ke
+  real*8 u_orig(ndim, 2), H_orig(2)
+  real*8 sqrt_rho(2), sqrt_rho_sum
+  real*8 u(ndim), rho, c, H, q
+  real*8 u_tilde(ndim), k(ndim), l(ndim), m(ndim)
   integer i, j
-  real b1, b2, alpha, beta
+  real*8 b1, b2, alpha, beta
 
   ik = d
   il = d+1
@@ -309,17 +309,17 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam
   end do
 
   do j=1,2
-    ke = 0.0
+    ke = 0.0d0
     do i=1,ndim
       ke = ke + u_orig(i,j)**2
     end do
 
-    p = (1.4 - 1.0)*(states(nvars,j) - 0.5*states(1,j)*ke)
+    p = (1.4d0 - 1)*(states(nvars,j) - 0.5d0*states(1,j)*ke)
     H_orig(j) = (states(nvars,j) + p)/states(1,j)
   end do
 
   do i=1,ndim
-    metric_norm(i) = 0.0
+    metric_norm(i) = 0.0d0
     do j=1,ndim
       metric_norm(i) = metric_norm(i) + metrics_frozen(i,j)**2
     end do
@@ -342,33 +342,33 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam
   rho = sqrt_rho(1)*sqrt_rho(2)
   H = (H_orig(1)*sqrt_rho(1) + H_orig(2)*sqrt_rho(2))/sqrt_rho_sum
 
-  q = 0.0
+  q = 0.0d0
   do i=1,ndim
     q = q + u(i)**2
   end do
 
-  c = sqrt((1.4 - 1.0)*(H - 0.5*q))
+  c = sqrt((1.4d0 - 1.0d0)*(H - 0.5d0*q))
 
-  b1 = (1.4 - 1.0)/(c**2)
-  b2 = 1.0 + b1*q - b1*H
+  b1 = (1.4d0 - 1.0d0)/(c**2)
+  b2 = 1.0d0 + b1*q - b1*H
 
-  u_tilde(1) = 0.0
+  u_tilde(1) = 0.0d0
   do i=1,ndim
     u_tilde(1) = u_tilde(1) + k(i)*u(i)
   end do
 
-  u_tilde(2) = 0.0
+  u_tilde(2) = 0.0d0
   do i=1,ndim
     u_tilde(2) = u_tilde(2) + l(i)*u(i)
   end do
 
-  u_tilde(3) = 0.0
+  u_tilde(3) = 0.0d0
   do i=1,ndim
     u_tilde(3) = u_tilde(3) + m(i)*u(i)
   end do
 
-  alpha = rho/(2.0*c)
-  beta = 1.0/(2.0*alpha)
+  alpha = rho/(2.0d0*c)
+  beta = 1.0d0/(2.0d0*alpha)
 
   lambda_roe(1) = u_tilde(1)*metric_norm(ik)
   lambda_roe(2) = u_tilde(1)*metric_norm(ik)
@@ -376,19 +376,19 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam
   lambda_roe(4) = u_tilde(1)*metric_norm(ik) + c*metric_norm(ik)
   lambda_roe(5) = u_tilde(1)*metric_norm(ik) - c*metric_norm(ik)
 
-  R(1,1) = 1.0
+  R(1,1) = 1.0d0
   R(2,1) = u(1)
   R(3,1) = u(2)
   R(4,1) = u(3)
-  R(5,1) = H - 1.0/b1
+  R(5,1) = H - 1.0d0/b1
 
-  R(1,2) = 0.0
+  R(1,2) = 0.0d0
   R(2,2) = rho*l(1)
   R(3,2) = rho*l(2)
   R(4,2) = rho*l(3)
   R(5,2) = rho*u_tilde(2)
 
-  R(1,3) = 0.0
+  R(1,3) = 0.0d0
   R(2,3) = rho*m(1)
   R(3,3) = rho*m(2)
   R(4,3) = rho*m(3)
@@ -406,7 +406,7 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam
   R(4,5) = alpha*(u(3) - c*k(3))
   R(5,5) = alpha*(H - c*u_tilde(1))
 
-  R_inv(1,1) = 1.0 - b2
+  R_inv(1,1) = 1.0d0 - b2
   R_inv(1,2) = b1*u(1)
   R_inv(1,3) = b1*u(2)
   R_inv(1,4) = b1*u(3)
@@ -416,13 +416,13 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam
   R_inv(2,2) = l(1)/rho
   R_inv(2,3) = l(2)/rho
   R_inv(2,4) = l(3)/rho
-  R_inv(2,5) = 0.0
+  R_inv(2,5) = 0.0d0
 
   R_inv(3,1) = -u_tilde(3)/rho
   R_inv(3,2) = m(1)/rho
   R_inv(3,3) = m(2)/rho
   R_inv(3,4) = m(3)/rho
-  R_inv(3,5) = 0.0
+  R_inv(3,5) = 0.0d0
 
   R_inv(4,1) = beta*(b2 - u_tilde(1)/c)
   R_inv(4,2) = -beta*(b1*u(1) - k(1)/c)
@@ -441,15 +441,15 @@ subroutine lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds)
   implicit none
 
   integer, intent(in) :: nvars
-  real, intent(in) :: lambda_pointwise(nvars, -2:3)
-  real, intent(in) :: lambda_roe(nvars)
-  real, intent(out) :: wavespeeds(nvars)
+  real*8, intent(in) :: lambda_pointwise(nvars, -2:3)
+  real*8, intent(in) :: lambda_roe(nvars)
+  real*8, intent(out) :: wavespeeds(nvars)
 
-  real kappa
+  real*8 kappa
   integer v
   integer k
 
-  kappa = 1.1
+  kappa = 1.1d0
 
   do v=1,nvars
     wavespeeds(v) = abs(lambda_roe(v))
@@ -471,32 +471,32 @@ subroutine convert_to_generalized_frozen( &
             generalized_states_frozen, &
             generalized_fluxes_frozen)
   integer, intent(in) :: nvars, ndim
-  real, intent(in) :: states(nvars, -2:3)
-  real, intent(in) :: fluxes(nvars, ndim, -2:3)
-  real, intent(in) :: metrics(ndim, ndim, 2)
-  real, intent(in) :: metric_jacobians(2)
-  real, intent(out) :: metrics_frozen(ndim, ndim)
-  real, intent(out) :: combined_frozen_metrics(ndim)
-  real, intent(out) :: generalized_states_frozen(nvars, -2:3)
-  real, intent(out) :: generalized_fluxes_frozen(nvars, ndim, -2:3)
+  real*8, intent(in) :: states(nvars, -2:3)
+  real*8, intent(in) :: fluxes(nvars, ndim, -2:3)
+  real*8, intent(in) :: metrics(ndim, ndim, 2)
+  real*8, intent(in) :: metric_jacobians(2)
+  real*8, intent(out) :: metrics_frozen(ndim, ndim)
+  real*8, intent(out) :: combined_frozen_metrics(ndim)
+  real*8, intent(out) :: generalized_states_frozen(nvars, -2:3)
+  real*8, intent(out) :: generalized_fluxes_frozen(nvars, ndim, -2:3)
 
-  real jacobian_frozen
+  real*8 jacobian_frozen
 
   integer c, k, v
   integer i, j
 
   do j=1,ndim
     do i=1,ndim
-      metrics_frozen(i,j) = 0.0
+      metrics_frozen(i,j) = 0.0d0
       do c=1,2
-        metrics_frozen(i,j) = metrics_frozen(i,j) + 0.5*metrics(i,j,c)/metric_jacobians(c)
+        metrics_frozen(i,j) = metrics_frozen(i,j) + 0.5d0*metrics(i,j,c)/metric_jacobians(c)
       end do
     end do
   end do
-  jacobian_frozen = 0.5*(metric_jacobians(1) + metric_jacobians(2))
+  jacobian_frozen = 0.5d0*(metric_jacobians(1) + metric_jacobians(2))
 
   do i=1,ndim
-    combined_frozen_metrics(i) = 0.0
+    combined_frozen_metrics(i) = 0.0d0
     do j=1,ndim
       combined_frozen_metrics(i) = combined_frozen_metrics(i) + metrics_frozen(i,j)**2
     end do
@@ -510,7 +510,7 @@ subroutine convert_to_generalized_frozen( &
 
   do k=-2,3
     do v=1,nvars
-      call mult_mat_vec(ndim, ndim, 1.0, metrics_frozen, &
+      call mult_mat_vec(ndim, ndim, 1.0d0, metrics_frozen, &
         fluxes(v,:,k), generalized_fluxes_frozen(v,:,k))
     end do
   end do
@@ -525,10 +525,10 @@ subroutine convert_to_generalized(nvars, ndim, nx, ny, nz, &
   implicit none
 
   integer, intent(in) :: nvars, ndim, nx, ny, nz
-  real, intent(in) :: fluxes(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
-  real, intent(in) :: metrics(ndim, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
-  real, intent(in) :: metric_jacobians(-2:nx+3, -2:ny+3, -2:nz+3)
-  real, intent(out) :: generalized_fluxes(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
+  real*8, intent(in) :: fluxes(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
+  real*8, intent(in) :: metrics(ndim, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
+  real*8, intent(in) :: metric_jacobians(-2:nx+3, -2:ny+3, -2:nz+3)
+  real*8, intent(out) :: generalized_fluxes(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
 
   integer i, j, k, v
 
@@ -536,7 +536,7 @@ subroutine convert_to_generalized(nvars, ndim, nx, ny, nz, &
     do j=-2,ny+3
       do i=-2,nx+3
         do v=1,nvars
-          call mult_mat_vec(ndim, ndim, 1.0/metric_jacobians(i,j,k), metrics(:,:,i,j,k), &
+          call mult_mat_vec(ndim, ndim, 1.0d0/metric_jacobians(i,j,k), metrics(:,:,i,j,k), &
             fluxes(v,:,i,j,k), generalized_fluxes(v,:,i,j,k))
         end do
       end do
@@ -552,10 +552,10 @@ subroutine convert_from_generalized(nvars, ndim, nx, ny, nz, &
   implicit none
 
   integer, intent(in) :: nvars, ndim, nx, ny, nz
-  real, intent(in) :: flux_derivatives_generalized(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
-  real, intent(in) :: metrics(ndim, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
-  real, intent(in) :: metric_jacobians(-2:nx+3, -2:ny+3, -2:nz+3)
-  real, intent(out) :: flux_derivatives(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
+  real*8, intent(in) :: flux_derivatives_generalized(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
+  real*8, intent(in) :: metrics(ndim, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
+  real*8, intent(in) :: metric_jacobians(-2:nx+3, -2:ny+3, -2:nz+3)
+  real*8, intent(out) :: flux_derivatives(nvars, ndim, -2:nx+3, -2:ny+3, -2:nz+3)
 
   integer i, j, k, v
   integer d
@@ -579,11 +579,11 @@ end subroutine
 subroutine solve3x3(a, b, x)
   implicit none
 
-  real, intent(in) :: a(3,3), b(3)
-  real, intent(out) :: x(3)
+  real*8, intent(in) :: a(3,3), b(3)
+  real*8, intent(out) :: x(3)
 
-  real temp(3,3)
-  real det_a(1), det_temp(1)
+  real*8 temp(3,3)
+  real*8 det_a(1), det_temp(1)
   integer k
   integer i,j
 
@@ -603,8 +603,8 @@ end subroutine
 subroutine determinant3x3(a, det)
   implicit none
 
-  real, intent(in) :: a(3,3)
-  real, intent(out) :: det(1)
+  real*8, intent(in) :: a(3,3)
+  real*8, intent(out) :: det(1)
 
   det(1) = a(1,1)*(a(2,2)*a(3,3) - a(2,3)*a(3,2)) &
         + a(1,2)*(a(2,3)*a(3,1) - a(2,1)*a(3,3)) &
@@ -621,26 +621,26 @@ subroutine split_characteristic_fluxes(nvars, &
   implicit none
 
   integer, intent(in) :: nvars
-  real, intent(in) :: generalized_states_frozen(nvars, -2:3)
-  real, intent(in) :: generalized_fluxes_frozen(nvars, -2:3)
-  real, intent(in) :: R_inv(nvars, nvars)
-  real, intent(in) :: wavespeeds(nvars)
-  real, intent(out) :: characteristic_fluxes_pos(nvars, -2:3)
-  real, intent(out) :: characteristic_fluxes_neg(nvars, -2:3)
+  real*8, intent(in) :: generalized_states_frozen(nvars, -2:3)
+  real*8, intent(in) :: generalized_fluxes_frozen(nvars, -2:3)
+  real*8, intent(in) :: R_inv(nvars, nvars)
+  real*8, intent(in) :: wavespeeds(nvars)
+  real*8, intent(out) :: characteristic_fluxes_pos(nvars, -2:3)
+  real*8, intent(out) :: characteristic_fluxes_neg(nvars, -2:3)
 
   integer k, m
   integer l
 
   do k=-2,3
     do m=1,nvars
-      characteristic_fluxes_pos(m,k) = 0.0
-      characteristic_fluxes_neg(m,k) = 0.0
+      characteristic_fluxes_pos(m,k) = 0.0d0
+      characteristic_fluxes_neg(m,k) = 0.0d0
       do l=1,nvars
         characteristic_fluxes_pos(m,k) = characteristic_fluxes_pos(m,k) &
-          + 0.5*R_inv(m,l) &
+          + 0.5d0*R_inv(m,l) &
           *(generalized_fluxes_frozen(l,k) + wavespeeds(m)*generalized_states_frozen(l,k))
         characteristic_fluxes_neg(m,k) = characteristic_fluxes_neg(m,k) &
-          + 0.5*R_inv(m,l) &
+          + 0.5d0*R_inv(m,l) &
           *(generalized_fluxes_frozen(l,k) - wavespeeds(m)*generalized_states_frozen(l,k))
       end do
     end do
@@ -652,16 +652,16 @@ subroutine weno_flux(nvars, generalized_fluxes, characteristic_fluxes_pos, &
   implicit none
 
   integer, intent(in) :: nvars
-  real, intent(in) :: generalized_fluxes(nvars, -2:3)
-  real, intent(in) :: characteristic_fluxes_pos(nvars, -2:3)
-  real, intent(in) :: characteristic_fluxes_neg(nvars, -2:3)
-  real, intent(in) :: combined_frozen_metrics
-  real, intent(in) :: R(nvars, nvars)
-  real, intent(out) :: flux(nvars)
-
-  real consistent(nvars)
-  real dissipation_pos(nvars)
-  real dissipation_neg(nvars)
+  real*8, intent(in) :: generalized_fluxes(nvars, -2:3)
+  real*8, intent(in) :: characteristic_fluxes_pos(nvars, -2:3)
+  real*8, intent(in) :: characteristic_fluxes_neg(nvars, -2:3)
+  real*8, intent(in) :: combined_frozen_metrics
+  real*8, intent(in) :: R(nvars, nvars)
+  real*8, intent(out) :: flux(nvars)
+
+  real*8 consistent(nvars)
+  real*8 dissipation_pos(nvars)
+  real*8 dissipation_neg(nvars)
   integer v
 
   call consistent_part(nvars, generalized_fluxes, consistent)
@@ -679,8 +679,8 @@ subroutine consistent_part(nvars, generalized_fluxes, consistent)
   implicit none
 
   integer, intent(in) :: nvars
-  real, intent(in) :: generalized_fluxes(nvars, -2:3)
-  real, intent(out) :: consistent(nvars)
+  real*8, intent(in) :: generalized_fluxes(nvars, -2:3)
+  real*8, intent(out) :: consistent(nvars)
 
   integer v
 
@@ -688,7 +688,7 @@ subroutine consistent_part(nvars, generalized_fluxes, consistent)
     consistent(v) = generalized_fluxes(v,-2) - 8*generalized_fluxes(v,-1) &
       + 37*generalized_fluxes(v,0) + 37*generalized_fluxes(v,1) &
       - 8*generalized_fluxes(v,2) + generalized_fluxes(v,3)
-    consistent(v) = consistent(v)/60.0
+    consistent(v) = consistent(v)/60.0d0
   end do
 end subroutine
 
@@ -697,16 +697,16 @@ subroutine dissipation_part_pos(nvars, characteristic_fluxes, &
   implicit none
 
   integer, intent(in) :: nvars
-  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
-  real, intent(in) :: combined_frozen_metrics
-  real, intent(in) :: R(nvars, nvars)
-  real, intent(out) :: dissipation_pos(nvars)
+  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+  real*8, intent(in) :: combined_frozen_metrics
+  real*8, intent(in) :: R(nvars, nvars)
+  real*8, intent(out) :: dissipation_pos(nvars)
 
-  real combined_fluxes(nvars)
+  real*8 combined_fluxes(nvars)
 
   call weno_combination_pos(nvars, characteristic_fluxes, combined_frozen_metrics, combined_fluxes)
 
-  call mult_mat_vec(nvars, nvars, -1.0/60, R, combined_fluxes, dissipation_pos)
+  call mult_mat_vec(nvars, nvars, -1.0d0/60, R, combined_fluxes, dissipation_pos)
 end subroutine
 
 subroutine weno_combination_pos(nvars, characteristic_fluxes, &
@@ -714,12 +714,12 @@ subroutine weno_combination_pos(nvars, characteristic_fluxes, &
   implicit none
 
   integer, intent(in) :: nvars
-  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
-  real, intent(in) :: combined_frozen_metrics
-  real, intent(out) :: combined_fluxes(nvars)
+  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+  real*8, intent(in) :: combined_frozen_metrics
+  real*8, intent(out) :: combined_fluxes(nvars)
 
-  real w(nvars, 3)
-  real flux_differences(nvars, 3)
+  real*8 w(nvars, 3)
+  real*8 flux_differences(nvars, 3)
   integer v
 
   call weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metrics, w)
@@ -736,26 +736,48 @@ subroutine weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metric
   implicit none
 
   integer, intent(in) :: nvars
-  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
-  real, intent(in) :: combined_frozen_metrics
-  real, intent(out) :: w(nvars, 3)
+  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+  real*8, intent(in) :: combined_frozen_metrics
+  real*8, intent(out) :: w(nvars, 3)
 
-  real :: C(3) = (/0.1, 0.6, 0.3/)
-  real :: weights1(0:2) = (/1, -4, 3/)
-  real :: weights2(0:2) = (/-1, 0, 1/)
-  real :: weights3(0:2) = (/-3, 4, 1/)
-  real :: weightsc(0:2) = (/1, -2, 1/)
-
-  real IS(3), alpha(3), eps
-
-  real w1sum(1), w2sum(1), w3sum(1), c1sum(1), c2sum(1), c3sum(1)
+  real*8 :: C(3) = (/0.1d0, 0.6d0, 0.3d0/)
+  real*8 IS(nvars, 3), alpha(3), eps
 
   integer p, i, j
-  real sum_alpha(1)
+  real*8 sum_alpha(1)
 
-  eps = 1e-6*combined_frozen_metrics
+  eps = 1.0d-6!*combined_frozen_metrics
   p = 2
 
+  call oscillation_pos(nvars, characteristic_fluxes, IS)
+
+  do i=1,nvars
+    sum_alpha(1) = 0.0d0
+    do j=1,3
+      alpha(j) = C(j)/(IS(i, j) + eps)**p
+      sum_alpha(1) = sum_alpha(1) + alpha(j)
+    end do
+    !call sum_array(3, alpha, sum_alpha)
+    do j=1,3
+      w(i,j) = alpha(j)/sum_alpha(1)
+    end do
+  end do
+end subroutine
+
+subroutine oscillation_pos(nvars, characteristic_fluxes, oscillation)
+  integer, intent(in) :: nvars
+  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+  real*8, intent(out) :: oscillation(nvars, 3)
+
+  real*8 :: weights1(0:2) = (/ 1.0d0, -4.0d0,  3.0d0/)
+  real*8 :: weights2(0:2) = (/-1.0d0,  0.0d0,  1.0d0/)
+  real*8 :: weights3(0:2) = (/-3.0d0,  4.0d0, -1.0d0/)
+  real*8 :: weightsc(0:2) = (/ 1.0d0, -2.0d0,  1.0d0/)
+
+  real*8 w1sum(1), w2sum(1), w3sum(1), c1sum(1), c2sum(1), c3sum(1)
+
+  integer i
+
   do i=1,nvars
     call weighted_sum(3, weights1, characteristic_fluxes(i,-2:0), w1sum)
     call weighted_sum(3, weightsc, characteristic_fluxes(i,-2:0), c1sum)
@@ -764,17 +786,9 @@ subroutine weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metric
     call weighted_sum(3, weights3, characteristic_fluxes(i,0:2), w3sum)
     call weighted_sum(3, weightsc, characteristic_fluxes(i,0:2), c3sum)
 
-    IS(1) = (1.0/4)*w1sum(1)**2 + (13.0/12)*c1sum(1)**2
-    IS(2) = (1.0/4)*w2sum(1)**2 + (13.0/12)*c2sum(1)**2
-    IS(3) = (1.0/4)*w3sum(1)**2 + (13.0/12)*c3sum(1)**2
-
-    do j=1,3
-      alpha(j) = C(j)/(IS(j) + eps)**p
-    end do
-    call sum_array(3, alpha, sum_alpha)
-    do j=1,3
-      w(i,j) = alpha(j)/sum_alpha(1)
-    end do
+    oscillation(i, 1) = (1.0d0/4)*w1sum(1)**2 + (13.0d0/12)*c1sum(1)**2
+    oscillation(i, 2) = (1.0d0/4)*w2sum(1)**2 + (13.0d0/12)*c2sum(1)**2
+    oscillation(i, 3) = (1.0d0/4)*w3sum(1)**2 + (13.0d0/12)*c3sum(1)**2
   end do
 end subroutine
 
@@ -782,8 +796,8 @@ subroutine flux_differences_pos(nvars, characteristic_fluxes, flux_differences)
   implicit none
 
   integer, intent(in) :: nvars
-  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
-  real, intent(out) :: flux_differences(nvars, 3)
+  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+  real*8, intent(out) :: flux_differences(nvars, 3)
 
   integer i, v
 
@@ -800,16 +814,16 @@ subroutine dissipation_part_neg(nvars, characteristic_fluxes, &
   implicit none
 
   integer, intent(in) :: nvars
-  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
-  real, intent(in) :: combined_frozen_metrics
-  real, intent(in) :: R(nvars, nvars)
-  real, intent(out) :: dissipation_neg(nvars)
+  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+  real*8, intent(in) :: combined_frozen_metrics
+  real*8, intent(in) :: R(nvars, nvars)
+  real*8, intent(out) :: dissipation_neg(nvars)
 
-  real combined_fluxes(nvars)
+  real*8 combined_fluxes(nvars)
 
   call weno_combination_neg(nvars, characteristic_fluxes, combined_frozen_metrics, combined_fluxes)
 
-  call mult_mat_vec(nvars, nvars, 1.0/60, R, combined_fluxes, dissipation_neg)
+  call mult_mat_vec(nvars, nvars, 1.0d0/60, R, combined_fluxes, dissipation_neg)
 end subroutine
 
 subroutine weno_combination_neg(nvars, characteristic_fluxes, &
@@ -817,12 +831,12 @@ subroutine weno_combination_neg(nvars, characteristic_fluxes, &
   implicit none
 
   integer, intent(in) :: nvars
-  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
-  real, intent(in) :: combined_frozen_metrics
-  real, intent(out) :: combined_fluxes(nvars)
+  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+  real*8, intent(in) :: combined_frozen_metrics
+  real*8, intent(out) :: combined_fluxes(nvars)
 
-  real w(nvars, 3)
-  real flux_differences(nvars, 3)
+  real*8 w(nvars, 3)
+  real*8 flux_differences(nvars, 3)
   integer v
 
   call weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metrics, w)
@@ -839,25 +853,48 @@ subroutine weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metric
   implicit none
 
   integer, intent(in) :: nvars
-  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
-  real, intent(in) :: combined_frozen_metrics
-  real, intent(out) :: w(nvars, 3)
-
-  real :: C(3) = (/0.1, 0.6, 0.3/)
-  real :: weights1(0:2) = (/1, -4, 3/)
-  real :: weights2(0:2) = (/-1, 0, 1/)
-  real :: weights3(0:2) = (/-3, 4, 1/)
-  real :: weightsc(0:2) = (/1, -2, 1/)
+  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+  real*8, intent(in) :: combined_frozen_metrics
+  real*8, intent(out) :: w(nvars, 3)
 
-  real IS(3), alpha(3), eps
-  real w1sum(1), w2sum(1), w3sum(1), c1sum(1), c2sum(1), c3sum(1)
+  real*8 :: C(3) = (/0.1d0, 0.6d0, 0.3d0/)
+  real*8 IS(nvars, 3), alpha(3), eps
 
   integer p, i, j
-  real sum_alpha(1)
+  real*8 sum_alpha(1)
 
-  eps = 1e-6*combined_frozen_metrics
+  eps = 1.0d-6!*combined_frozen_metrics
   p = 2
 
+  call oscillation_neg(nvars, characteristic_fluxes, IS)
+
+  do i=1,nvars
+    do j=1,3
+      alpha(j) = C(j)/(IS(i, j) + eps)**p
+    end do
+    call sum_array(3, alpha, sum_alpha)
+    do j=1,3
+      w(i,j) = alpha(j)/sum_alpha(1)
+    end do
+  end do
+end subroutine
+
+subroutine oscillation_neg(nvars, characteristic_fluxes, oscillation)
+  implicit none
+
+  integer, intent(in) :: nvars
+  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+  real*8, intent(out) :: oscillation(nvars, 3)
+
+  real*8 :: weights1(0:2) = (/ 1.0d0, -4.0d0,  3.0d0/)
+  real*8 :: weights2(0:2) = (/-1.0d0,  0.0d0,  1.0d0/)
+  real*8 :: weights3(0:2) = (/-3.0d0,  4.0d0, -1.0d0/)
+  real*8 :: weightsc(0:2) = (/ 1.0d0, -2.0d0,  1.0d0/)
+
+  real*8 w1sum(1), w2sum(1), w3sum(1), c1sum(1), c2sum(1), c3sum(1)
+
+  integer i
+
   do i=1,nvars
     call weighted_sum(3, weights1, characteristic_fluxes(i,3:1:-1), w1sum)
     call weighted_sum(3, weightsc, characteristic_fluxes(i,3:1:-1), c1sum)
@@ -866,17 +903,9 @@ subroutine weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metric
     call weighted_sum(3, weights3, characteristic_fluxes(i,1:-1:-1), w3sum)
     call weighted_sum(3, weightsc, characteristic_fluxes(i,1:-1:-1), c3sum)
 
-    IS(1) = (1.0/4)*w1sum(1)**2 + (13.0/12)*c1sum(1)**2
-    IS(2) = (1.0/4)*w2sum(1)**2 + (13.0/12)*c2sum(1)**2
-    IS(3) = (1.0/4)*w3sum(1)**2 + (13.0/12)*c3sum(1)**2
-
-    do j=1,3
-      alpha(j) = C(j)/(IS(j) + eps)**p
-    end do
-    call sum_array(3, alpha, sum_alpha)
-    do j=1,3
-      w(i,j) = alpha(j)/sum_alpha(1)
-    end do
+    oscillation(i, 1) = (1.0d0/4)*w1sum(1)**2 + (13.0d0/12)*c1sum(1)**2
+    oscillation(i, 2) = (1.0d0/4)*w2sum(1)**2 + (13.0d0/12)*c2sum(1)**2
+    oscillation(i, 3) = (1.0d0/4)*w3sum(1)**2 + (13.0d0/12)*c3sum(1)**2
   end do
 end subroutine
 
@@ -884,8 +913,8 @@ subroutine flux_differences_neg(nvars, characteristic_fluxes, flux_differences)
   implicit none
 
   integer, intent(in) :: nvars
-  real, intent(in) :: characteristic_fluxes(nvars, -2:3)
-  real, intent(out) :: flux_differences(nvars, 3)
+  real*8, intent(in) :: characteristic_fluxes(nvars, -2:3)
+  real*8, intent(out) :: flux_differences(nvars, 3)
 
   integer i, v
 
@@ -901,16 +930,16 @@ subroutine mult_mat_vec(m, n, alpha, a, b, c)
   implicit none
 
   integer, intent(in) :: m, n
-  real, intent(in) :: alpha
-  real, intent(in) :: a(m, n)
-  real, intent(in) :: b(n)
-  real, intent(out) :: c(m)
+  real*8, intent(in) :: alpha
+  real*8, intent(in) :: a(m, n)
+  real*8, intent(in) :: b(n)
+  real*8, intent(out) :: c(m)
 
   integer i, j
-  real accumulator
+  real*8 accumulator
 
   do i=1,m
-    accumulator = 0.0
+    accumulator = 0.0d0
     do j=1,n
       accumulator = accumulator + alpha*a(i,j)*b(j)
     end do
@@ -922,12 +951,12 @@ subroutine sum_array(n, a, a_sum)
   implicit none
 
   integer, intent(in) :: n
-  real, intent(in) :: a(n)
-  real, intent(out) :: a_sum(1)
+  real*8, intent(in) :: a(n)
+  real*8, intent(out) :: a_sum(1)
 
   integer i
 
-  a_sum(1) = 0.0
+  a_sum(1) = 0.0d0
   do i=1,n
     a_sum(1) = a_sum(1) + a(i)
   end do
@@ -937,12 +966,12 @@ subroutine weighted_sum(n, w, a, a_sum)
   implicit none
 
   integer, intent(in) :: n
-  real, intent(in) :: a(n), w(n)
-  real, intent(out) :: a_sum(1)
+  real*8, intent(in) :: a(n), w(n)
+  real*8, intent(out) :: a_sum(1)
 
   integer i
 
-  a_sum(1) = 0.0
+  a_sum(1) = 0.0d0
   do i=1,n
     a_sum(1) = a_sum(1) + w(i)*a(i)
   end do
diff --git a/data_for_test.py b/data_for_test.py
new file mode 100644
index 0000000000000000000000000000000000000000..bde2cd974fb711349107e23ebf51e308b3c05316
--- /dev/null
+++ b/data_for_test.py
@@ -0,0 +1,938 @@
+import pytest
+import utilities as u
+
+
+# {{{ FluxDataSingle
+
+class FluxDataSingle:
+    nvars = 5
+    ndim = 3
+    dirs = {"x": 1, "y": 2, "z": 3}
+
+    def __init__(self, states_str, fluxes_str,
+            lam_str, R_str, R_inv_str, wavespeeds_str,
+            char_fluxes_pos_str, char_fluxes_neg_str,
+            oscillation_pos_str, oscillation_neg_str,
+            weno_weights_pos_str, weno_weights_neg_str,
+            consistent_str, dissipation_pos_str, dissipation_neg_str,
+            weno_flux_str, direction):
+
+        self.direction = self.dirs[direction]
+
+        self.state_pair = self.swap_array_rows(
+                u.transposed_array_from_string(states_str), self.direction)
+        self.states = u.expand_to_6(self.state_pair)
+
+        self.flux_pair = self.swap_array_rows(
+                u.transposed_array_from_string(fluxes_str), self.direction)
+        self.fluxes = u.expand_to_6(self.flux_pair)
+
+        self.lam_pointwise = u.expand_to_6(
+                u.transposed_array_from_string(lam_str))
+
+        self.R = self.swap_array_rows(
+                u.array_from_string(R_str), self.direction)
+        self.R_inv = self.swap_array_cols(
+                u.array_from_string(R_inv_str), self.direction)
+        self.wavespeeds = u.array_from_string(wavespeeds_str)
+
+        self.char_fluxes_pos = u.expand_to_6(
+                u.transposed_array_from_string(char_fluxes_pos_str))
+        self.char_fluxes_neg = u.expand_to_6(
+                u.transposed_array_from_string(char_fluxes_neg_str))
+
+        self.oscillation_pos = u.transposed_array_from_string(oscillation_pos_str)
+        self.oscillation_neg = u.transposed_array_from_string(oscillation_neg_str)
+
+        self.weno_weights_pos = u.transposed_array_from_string(weno_weights_pos_str)
+        self.weno_weights_neg = u.transposed_array_from_string(weno_weights_neg_str)
+
+        self.consistent = self.swap_array(
+                u.array_from_string(consistent_str), self.direction)
+        self.dissipation_pos = self.swap_array(
+                u.array_from_string(dissipation_pos_str), self.direction)
+        self.dissipation_neg = self.swap_array(
+                u.array_from_string(dissipation_neg_str), self.direction)
+
+        self.weno_flux = self.swap_array(
+                u.array_from_string(weno_flux_str), self.direction)
+
+    def swap_array(self, arr, d):
+        p = self.permutation(d)
+        arr[p] = arr[[1, 2, 3]]
+        return arr
+
+    def swap_array_rows(self, arr, d):
+        p = self.permutation(d)
+        arr[p, :] = arr[[1, 2, 3], :]
+        return arr
+
+    def swap_array_cols(self, arr, d):
+        p = self.permutation(d)
+        arr[:, p] = arr[:, [1, 2, 3]]
+        return arr
+
+    def permutation(self, d):
+        return [(d-1+i) % 3 + 1 for i in range(3)]
+
+# }}}
+
+
+single_data = {}
+
+# {{{ Input data: "flat" case
+
+single_data["Case flat:x"] = FluxDataSingle(
+        states_str="1 1 1 1 5.5,1 1 1 1 5.5",
+        fluxes_str="1 2.6 1 1 7.1,1 2.6 1 1 7.1",
+        lam_str=("1. 1. 1. 2.4966629547095764 -0.49666295470957644,"
+            "1. 1. 1. 2.4966629547095764 -0.49666295470957644"),
+        R_str=("1 0 0 0.3340765523905305 0.3340765523905305,"
+            "1 0 0 0.8340765523905306 -0.1659234476094695,"
+            "1 1 0 0.3340765523905305 0.3340765523905305,"
+            "1 0 1 0.3340765523905305 0.3340765523905305,"
+            "1.5000000000000009 1 1 2.8719435219727663 1.8719435219727667"),
+        R_inv_str=("0.7321428571428572 0.1785714285714286 0.1785714285714286 "
+            "0.1785714285714286 -0.1785714285714286,"
+            "-1. 0 1. 0 0,"
+            "-1. 0 0 1. 0,"
+            "-0.5991081371313636 0.7327387580875756 -0.2672612419124244 "
+            "-0.2672612419124244 0.2672612419124244,"
+            "1.4008918628686364 -1.2672612419124245 -0.2672612419124244 "
+            "-0.2672612419124244 0.2672612419124244"),
+        wavespeeds_str="1.1 1.1 1.1 2.7463292501805343 0.5463292501805341",
+        char_fluxes_pos_str=(
+            "0.30000000000000004 0. 0. 2.8024972160321826 0.026547751617514326,"
+            "0.30000000000000004 0. 0. 2.8024972160321826 0.026547751617514326"),
+        char_fluxes_neg_str=(
+            "-0.014285714285714235 0. 0. -0.13345224838248493 -0.5575027839678178,"
+            "-0.014285714285714235 0. 0. -0.13345224838248493 -0.5575027839678178"),
+        oscillation_pos_str=("0 0 0 0 0,0 0 0 0 0,0 0 0 0 0"),
+        oscillation_neg_str=("0 0 0 0 0,0 0 0 0 0,0 0 0 0 0"),
+        weno_weights_pos_str=("0.1 0.1 0.1 0.1 0.1,"
+            "0.6 0.6 0.6 0.6 0.6,0.3 0.3 0.3 0.3 0.3"),
+        weno_weights_neg_str=("0.1 0.1 0.1 0.1 0.1,"
+            "0.6 0.6 0.6 0.6 0.6,0.3 0.3 0.3 0.3 0.3"),
+        consistent_str="1 2.6 1 1 7.1",
+        dissipation_pos_str=("0 0 0 0 0"),
+        dissipation_neg_str=("0 0 0 0 0"),
+        weno_flux_str=("1 2.6 1 1 7.1"),
+        direction="x")
+single_data["Case flat:y"] = FluxDataSingle(
+        states_str="1 1 1 1 5.5,1 1 1 1 5.5",
+        fluxes_str="1 2.6 1 1 7.1,1 2.6 1 1 7.1",
+        lam_str=("1. 1. 1. 2.4966629547095764 -0.49666295470957644,"
+            "1. 1. 1. 2.4966629547095764 -0.49666295470957644"),
+        R_str=("1 0 0 0.3340765523905305 0.3340765523905305,"
+            "1 0 0 0.8340765523905306 -0.1659234476094695,"
+            "1 1 0 0.3340765523905305 0.3340765523905305,"
+            "1 0 1 0.3340765523905305 0.3340765523905305,"
+            "1.5000000000000009 1 1 2.8719435219727663 1.8719435219727667"),
+        R_inv_str=("0.7321428571428572 0.1785714285714286 0.1785714285714286 "
+            "0.1785714285714286 -0.1785714285714286,"
+            "-1. 0 1. 0 0,"
+            "-1. 0 0 1. 0,"
+            "-0.5991081371313636 0.7327387580875756 -0.2672612419124244 "
+            "-0.2672612419124244 0.2672612419124244,"
+            "1.4008918628686364 -1.2672612419124245 -0.2672612419124244 "
+            "-0.2672612419124244 0.2672612419124244"),
+        wavespeeds_str="1.1 1.1 1.1 2.7463292501805343 0.5463292501805341",
+        char_fluxes_pos_str=(
+            "0.30000000000000004 0. 0. 2.8024972160321826 0.026547751617514326,"
+            "0.30000000000000004 0. 0. 2.8024972160321826 0.026547751617514326"),
+        char_fluxes_neg_str=(
+            "-0.014285714285714235 0. 0. -0.13345224838248493 -0.5575027839678178,"
+            "-0.014285714285714235 0. 0. -0.13345224838248493 -0.5575027839678178"),
+        oscillation_pos_str=("0 0 0 0 0,0 0 0 0 0,0 0 0 0 0"),
+        oscillation_neg_str=("0 0 0 0 0,0 0 0 0 0,0 0 0 0 0"),
+        weno_weights_pos_str=("0.1 0.1 0.1 0.1 0.1,"
+            "0.6 0.6 0.6 0.6 0.6,0.3 0.3 0.3 0.3 0.3"),
+        weno_weights_neg_str=("0.1 0.1 0.1 0.1 0.1,"
+            "0.6 0.6 0.6 0.6 0.6,0.3 0.3 0.3 0.3 0.3"),
+        consistent_str="1 2.6 1 1 7.1",
+        dissipation_pos_str=("0 0 0 0 0"),
+        dissipation_neg_str=("0 0 0 0 0"),
+        weno_flux_str=("1 2.6 1 1 7.1"),
+        direction="y")
+single_data["Case flat:z"] = FluxDataSingle(
+        states_str="1 1 1 1 5.5,1 1 1 1 5.5",
+        fluxes_str="1 2.6 1 1 7.1,1 2.6 1 1 7.1",
+        lam_str=("1. 1. 1. 2.4966629547095764 -0.49666295470957644,"
+            "1. 1. 1. 2.4966629547095764 -0.49666295470957644"),
+        R_str=("1 0 0 0.3340765523905305 0.3340765523905305,"
+            "1 0 0 0.8340765523905306 -0.1659234476094695,"
+            "1 1 0 0.3340765523905305 0.3340765523905305,"
+            "1 0 1 0.3340765523905305 0.3340765523905305,"
+            "1.5000000000000009 1 1 2.8719435219727663 1.8719435219727667"),
+        R_inv_str=("0.7321428571428572 0.1785714285714286 0.1785714285714286 "
+            "0.1785714285714286 -0.1785714285714286,"
+            "-1. 0 1. 0 0,"
+            "-1. 0 0 1. 0,"
+            "-0.5991081371313636 0.7327387580875756 -0.2672612419124244 "
+            "-0.2672612419124244 0.2672612419124244,"
+            "1.4008918628686364 -1.2672612419124245 -0.2672612419124244 "
+            "-0.2672612419124244 0.2672612419124244"),
+        wavespeeds_str="1.1 1.1 1.1 2.7463292501805343 0.5463292501805341",
+        char_fluxes_pos_str=(
+            "0.30000000000000004 0. 0. 2.8024972160321826 0.026547751617514326,"
+            "0.30000000000000004 0. 0. 2.8024972160321826 0.026547751617514326"),
+        char_fluxes_neg_str=(
+            "-0.014285714285714235 0. 0. -0.13345224838248493 -0.5575027839678178,"
+            "-0.014285714285714235 0. 0. -0.13345224838248493 -0.5575027839678178"),
+        oscillation_pos_str=("0 0 0 0 0,0 0 0 0 0,0 0 0 0 0"),
+        oscillation_neg_str=("0 0 0 0 0,0 0 0 0 0,0 0 0 0 0"),
+        weno_weights_pos_str=("0.1 0.1 0.1 0.1 0.1,"
+            "0.6 0.6 0.6 0.6 0.6,0.3 0.3 0.3 0.3 0.3"),
+        weno_weights_neg_str=("0.1 0.1 0.1 0.1 0.1,"
+            "0.6 0.6 0.6 0.6 0.6,0.3 0.3 0.3 0.3 0.3"),
+        consistent_str="1 2.6 1 1 7.1",
+        dissipation_pos_str=("0 0 0 0 0"),
+        dissipation_neg_str=("0 0 0 0 0"),
+        weno_flux_str=("1 2.6 1 1 7.1"),
+        direction="z")
+
+# }}}
+
+# {{{ Input data: Case (a)
+
+single_data["Case a:x"] = FluxDataSingle(
+        states_str="2 4 4 4 20,1 1 1 1 5.5",
+        fluxes_str="4 11.2 8 8 46.4,1 2.6 1 1 7.1",
+        lam_str=("2. 2. 2. 3.4966629547095764 0.5033370452904236,"
+            "1. 1. 1. 2.4966629547095764 -0.49666295470957644"),
+        R_str=("1 0 0 0.45781245952886396 0.45781245952886396,"
+            "1.5857864376269053 0 0 1.4330995704840366 0.018886008110941373,"
+            "1.5857864376269053 1.4142135623730951 0 0.7259927892974889 "
+            "0.7259927892974889,"
+            "1.5857864376269053 0 1.4142135623730951 0.7259927892974889 "
+            "0.7259927892974889,"
+            "3.7720779386421466 2.2426406871192857 2.2426406871192857 "
+            "5.578600290173388 3.335959603054103"),
+        R_inv_str=("0.36752136386566203 0.265894835574803 0.265894835574803 "
+            "0.265894835574803 -0.16767379847989416,"
+            "-1.1213203435596428 0 0.7071067811865475 0 0,"
+            "-1.1213203435596428 0 0 0.7071067811865475 0,"
+            "-0.43055863210991147 0.41670966546755195 -0.2903971157189956 "
+            "-0.2903971157189956 0.1831249838115456,"
+            "1.8120820550093746 -0.9975038969055432 -0.2903971157189956 "
+            "-0.2903971157189956 0.1831249838115456"),
+        wavespeeds_str="2.2 2.2 2.2 3.8463292501805344 0.553670749819466",
+        char_fluxes_pos_str=("1.0907156303492833 1.2301515190164993 "
+            "1.2301515190164993 7.523052586013577 0.23295627081394965,"
+            "0.4673767959969717 -0.6627416997969526 "
+            "-0.6627416997969526 1.4795302623674356 0.3125273038443238"),
+        char_fluxes_neg_str=("-0.16835489671908455 -0.05857864376269045 "
+            "-0.05857864376269045 -0.7274934638648571 -0.30602629876675014,"
+            "-0.06722315769446469 0.24852813742385726 "
+            "0.24852813742385726 -0.10725061063772967 -0.37456222716537935"),
+        oscillation_pos_str=("0. 0. 0. 0. 0.,"
+            "0.5180684032155981 4.777392983773268 4.777392983773268 "
+            "48.69888276854567 0.00844206573002786,"
+            "1.2951710080389953 11.94348245943317 11.94348245943317 "
+            "121.7472069213642 0.02110516432506965"),
+        oscillation_neg_str=("0. 0. 0. 0. 0.,"
+            "0.01363683818419176 0.12575276673434946 0.12575276673434946 "
+            "0.5129349293057707 0.006262897975282706,"
+            "0.0340920954604794 0.31438191683587363 0.31438191683587363 "
+            "1.2823373232644266 0.015657244938206753"),
+        weno_weights_pos_str=("0.999999997585736 0.9999999999716082 "
+            "0.9999999999716082 0.9999999999997268 0.9999909282503908,"
+            "2.2354258683870077e-9 2.6288602366075745e-11 "
+            "2.6288602366075745e-11 2.5299566294114733e-13 8.39888391917283e-6,"
+            "1.788382117901561e-10 2.1030934718853924e-12 "
+            "2.1030934718853924e-12 2.0239658022589593e-14 6.728656898188601e-7"),
+        weno_weights_neg_str=("0.9999965203327451 0.999999959029252 "
+            "0.999999959029252 0.9999999975371708 0.9999835300215946,"
+            "3.2217041402391125e-6 3.793560951911702e-8 "
+            "3.793560951911702e-8 2.2803933579018906e-9 0.000015247816238787948,"
+            "2.579631146567429e-7 3.0351383606886035e-9 "
+            "3.0351383606886035e-9 1.824357365679662e-10 1.2221621667041045e-6"),
+        consistent_str="2.5 6.9 4.5 4.5 26.75",
+        dissipation_pos_str=("1.6768551393982627 4.8239743723334865 "
+            "3.997611768980551 3.997611768980551 22.145196359962693"),
+        dissipation_neg_str=("0.1768550804014161 0.523974175351619 "
+            "0.4976116690805634 0.4976116690805634 2.495196349669178"),
+        weno_flux_str=("4.353710219799678 12.247948547685105 8.995223438061114 "
+            "8.995223438061114 51.39039270963187"),
+        direction="x")
+single_data["Case a:y"] = FluxDataSingle(
+        states_str="2 4 4 4 20,1 1 1 1 5.5",
+        fluxes_str="4 11.2 8 8 46.4,1 2.6 1 1 7.1",
+        lam_str=("2. 2. 2. 3.4966629547095764 0.5033370452904236,"
+            "1. 1. 1. 2.4966629547095764 -0.49666295470957644"),
+        R_str=("1 0 0 0.45781245952886396 0.45781245952886396,"
+            "1.5857864376269053 0 0 1.4330995704840366 0.018886008110941373,"
+            "1.5857864376269053 1.4142135623730951 0 0.7259927892974889 "
+            "0.7259927892974889,"
+            "1.5857864376269053 0 1.4142135623730951 0.7259927892974889 "
+            "0.7259927892974889,"
+            "3.7720779386421466 2.2426406871192857 2.2426406871192857 "
+            "5.578600290173388 3.335959603054103"),
+        R_inv_str=("0.36752136386566203 0.265894835574803 0.265894835574803 "
+            "0.265894835574803 -0.16767379847989416,"
+            "-1.1213203435596428 0 0.7071067811865475 0 0,"
+            "-1.1213203435596428 0 0 0.7071067811865475 0,"
+            "-0.43055863210991147 0.41670966546755195 -0.2903971157189956 "
+            "-0.2903971157189956 0.1831249838115456,"
+            "1.8120820550093746 -0.9975038969055432 -0.2903971157189956 "
+            "-0.2903971157189956 0.1831249838115456"),
+        wavespeeds_str="2.2 2.2 2.2 3.8463292501805344 0.553670749819466",
+        char_fluxes_pos_str=("1.0907156303492833 1.2301515190164993 "
+            "1.2301515190164993 7.523052586013577 0.23295627081394965,"
+            "0.4673767959969717 -0.6627416997969526 "
+            "-0.6627416997969526 1.4795302623674356 0.3125273038443238"),
+        char_fluxes_neg_str=("-0.16835489671908455 -0.05857864376269045 "
+            "-0.05857864376269045 -0.7274934638648571 -0.30602629876675014,"
+            "-0.06722315769446469 0.24852813742385726 "
+            "0.24852813742385726 -0.10725061063772967 -0.37456222716537935"),
+        oscillation_pos_str=("0. 0. 0. 0. 0.,"
+            "0.5180684032155981 4.777392983773268 4.777392983773268 "
+            "48.69888276854567 0.00844206573002786,"
+            "1.2951710080389953 11.94348245943317 11.94348245943317 "
+            "121.7472069213642 0.02110516432506965"),
+        oscillation_neg_str=("0. 0. 0. 0. 0.,"
+            "0.01363683818419176 0.12575276673434946 0.12575276673434946 "
+            "0.5129349293057707 0.006262897975282706,"
+            "0.0340920954604794 0.31438191683587363 0.31438191683587363 "
+            "1.2823373232644266 0.015657244938206753"),
+        weno_weights_pos_str=("0.999999997585736 0.9999999999716082 "
+            "0.9999999999716082 0.9999999999997268 0.9999909282503908,"
+            "2.2354258683870077e-9 2.6288602366075745e-11 "
+            "2.6288602366075745e-11 2.5299566294114733e-13 8.39888391917283e-6,"
+            "1.788382117901561e-10 2.1030934718853924e-12 "
+            "2.1030934718853924e-12 2.0239658022589593e-14 6.728656898188601e-7"),
+        weno_weights_neg_str=("0.9999965203327451 0.999999959029252 "
+            "0.999999959029252 0.9999999975371708 0.9999835300215946,"
+            "3.2217041402391125e-6 3.793560951911702e-8 "
+            "3.793560951911702e-8 2.2803933579018906e-9 0.000015247816238787948,"
+            "2.579631146567429e-7 3.0351383606886035e-9 "
+            "3.0351383606886035e-9 1.824357365679662e-10 1.2221621667041045e-6"),
+        consistent_str="2.5 6.9 4.5 4.5 26.75",
+        dissipation_pos_str=("1.6768551393982627 4.8239743723334865 "
+            "3.997611768980551 3.997611768980551 22.145196359962693"),
+        dissipation_neg_str=("0.1768550804014161 0.523974175351619 "
+            "0.4976116690805634 0.4976116690805634 2.495196349669178"),
+        weno_flux_str=("4.353710219799678 12.247948547685105 8.995223438061114 "
+            "8.995223438061114 51.39039270963187"),
+        direction="y")
+single_data["Case a:z"] = FluxDataSingle(
+        states_str="2 4 4 4 20,1 1 1 1 5.5",
+        fluxes_str="4 11.2 8 8 46.4,1 2.6 1 1 7.1",
+        lam_str=("2. 2. 2. 3.4966629547095764 0.5033370452904236,"
+            "1. 1. 1. 2.4966629547095764 -0.49666295470957644"),
+        R_str=("1 0 0 0.45781245952886396 0.45781245952886396,"
+            "1.5857864376269053 0 0 1.4330995704840366 0.018886008110941373,"
+            "1.5857864376269053 1.4142135623730951 0 0.7259927892974889 "
+            "0.7259927892974889,"
+            "1.5857864376269053 0 1.4142135623730951 0.7259927892974889 "
+            "0.7259927892974889,"
+            "3.7720779386421466 2.2426406871192857 2.2426406871192857 "
+            "5.578600290173388 3.335959603054103"),
+        R_inv_str=("0.36752136386566203 0.265894835574803 0.265894835574803 "
+            "0.265894835574803 -0.16767379847989416,"
+            "-1.1213203435596428 0 0.7071067811865475 0 0,"
+            "-1.1213203435596428 0 0 0.7071067811865475 0,"
+            "-0.43055863210991147 0.41670966546755195 -0.2903971157189956 "
+            "-0.2903971157189956 0.1831249838115456,"
+            "1.8120820550093746 -0.9975038969055432 -0.2903971157189956 "
+            "-0.2903971157189956 0.1831249838115456"),
+        wavespeeds_str="2.2 2.2 2.2 3.8463292501805344 0.553670749819466",
+        char_fluxes_pos_str=("1.0907156303492833 1.2301515190164993 "
+            "1.2301515190164993 7.523052586013577 0.23295627081394965,"
+            "0.4673767959969717 -0.6627416997969526 "
+            "-0.6627416997969526 1.4795302623674356 0.3125273038443238"),
+        char_fluxes_neg_str=("-0.16835489671908455 -0.05857864376269045 "
+            "-0.05857864376269045 -0.7274934638648571 -0.30602629876675014,"
+            "-0.06722315769446469 0.24852813742385726 "
+            "0.24852813742385726 -0.10725061063772967 -0.37456222716537935"),
+        oscillation_pos_str=("0. 0. 0. 0. 0.,"
+            "0.5180684032155981 4.777392983773268 4.777392983773268 "
+            "48.69888276854567 0.00844206573002786,"
+            "1.2951710080389953 11.94348245943317 11.94348245943317 "
+            "121.7472069213642 0.02110516432506965"),
+        oscillation_neg_str=("0. 0. 0. 0. 0.,"
+            "0.01363683818419176 0.12575276673434946 0.12575276673434946 "
+            "0.5129349293057707 0.006262897975282706,"
+            "0.0340920954604794 0.31438191683587363 0.31438191683587363 "
+            "1.2823373232644266 0.015657244938206753"),
+        weno_weights_pos_str=("0.999999997585736 0.9999999999716082 "
+            "0.9999999999716082 0.9999999999997268 0.9999909282503908,"
+            "2.2354258683870077e-9 2.6288602366075745e-11 "
+            "2.6288602366075745e-11 2.5299566294114733e-13 8.39888391917283e-6,"
+            "1.788382117901561e-10 2.1030934718853924e-12 "
+            "2.1030934718853924e-12 2.0239658022589593e-14 6.728656898188601e-7"),
+        weno_weights_neg_str=("0.9999965203327451 0.999999959029252 "
+            "0.999999959029252 0.9999999975371708 0.9999835300215946,"
+            "3.2217041402391125e-6 3.793560951911702e-8 "
+            "3.793560951911702e-8 2.2803933579018906e-9 0.000015247816238787948,"
+            "2.579631146567429e-7 3.0351383606886035e-9 "
+            "3.0351383606886035e-9 1.824357365679662e-10 1.2221621667041045e-6"),
+        consistent_str="2.5 6.9 4.5 4.5 26.75",
+        dissipation_pos_str=("1.6768551393982627 4.8239743723334865 "
+            "3.997611768980551 3.997611768980551 22.145196359962693"),
+        dissipation_neg_str=("0.1768550804014161 0.523974175351619 "
+            "0.4976116690805634 0.4976116690805634 2.495196349669178"),
+        weno_flux_str=("4.353710219799678 12.247948547685105 8.995223438061114 "
+            "8.995223438061114 51.39039270963187"),
+        direction="z")
+
+# }}}
+
+# {{{ Input data: Case (b)
+
+single_data["Case b:x"] = FluxDataSingle(
+        states_str="1 -1 -1 -1 5.5,2 -4 -4 -4 20",
+        fluxes_str="-1 2.6 1 1 -7.1,-4 11.2 8 8 -46.4",
+        lam_str=("-1. -1. -1. 0.49666295470957644 -2.4966629547095764,"
+            "-2. -2. -2. -0.5033370452904236 -3.4966629547095764"),
+        R_str=("1 0 0 0.45781245952886396 0.45781245952886396,"
+            "-1.5857864376269053 0 0 -0.018886008110941373 "
+            "-1.4330995704840366,"
+            "-1.5857864376269053 1.4142135623730951 0 -0.7259927892974889 "
+            "-0.7259927892974889,"
+            "-1.5857864376269053 0 1.4142135623730951 -0.7259927892974889 "
+            "-0.7259927892974889,"
+            "3.7720779386421466 -2.2426406871192857 -2.2426406871192857 "
+            "3.335959603054103 5.578600290173388"),
+        R_inv_str=("0.36752136386566203 -0.265894835574803 -0.265894835574803 "
+            "-0.265894835574803 -0.16767379847989416,"
+            "1.1213203435596428 0 0.7071067811865475 0 0,"
+            "1.1213203435596428 0 0 0.7071067811865475 0,"
+            "1.8120820550093746 0.9975038969055432 0.2903971157189956 "
+            "0.2903971157189956 0.1831249838115456,"
+            "-0.43055863210991147 -0.41670966546755195 0.2903971157189956 "
+            "0.2903971157189956 0.1831249838115456"),
+        wavespeeds_str="2.2 2.2 2.2 0.553670749819466 3.8463292501805344",
+        char_fluxes_pos_str=("0.06722315769446469 0.24852813742385726 "
+            "0.24852813742385726 0.37456222716537935 0.10725061063772967,"
+            "0.16835489671908455 -0.05857864376269045 "
+            "-0.05857864376269045 0.30602629876675014 0.7274934638648571"),
+        char_fluxes_neg_str=("-0.4673767959969717 -0.6627416997969526 "
+            "-0.6627416997969526 -0.3125273038443238 -1.4795302623674356,"
+            "-1.0907156303492833 1.2301515190164993 "
+            "1.2301515190164993 -0.23295627081394965 -7.523052586013577"),
+        oscillation_pos_str=("0. 0. 0. 0. 0.,"
+            "0.01363683818419176 0.12575276673434946 0.12575276673434946 "
+            "0.006262897975282706 0.5129349293057707,"
+            "0.0340920954604794 0.31438191683587363 0.31438191683587363 "
+            "0.015657244938206753 1.2823373232644266"),
+        oscillation_neg_str=("0. 0. 0. 0. 0.,"
+            "0.5180684032155981 4.777392983773268 4.777392983773268 "
+            "0.00844206573002786 48.69888276854567,"
+            "1.2951710080389953 11.94348245943317 11.94348245943317 "
+            "0.02110516432506965 121.7472069213642"),
+        weno_weights_pos_str=("0.9999965203327451 0.999999959029252 "
+            "0.999999959029252 0.9999835300215946 0.9999999975371708,"
+            "3.2217041402391125e-6 3.793560951911702e-8 "
+            "3.793560951911702e-8 0.000015247816238787948 2.2803933579018906e-9,"
+            "2.579631146567429e-7 3.0351383606886035e-9 "
+            "3.0351383606886035e-9 1.2221621667041045e-6 1.824357365679662e-10"),
+        weno_weights_neg_str=("0.999999997585736 0.9999999999716082 "
+            "0.9999999999716082 0.9999909282503908 0.9999999999997268,"
+            "2.2354258683870077e-9 2.6288602366075745e-11 "
+            "2.6288602366075745e-11 8.39888391917283e-6 2.5299566294114733e-13,"
+            "1.788382117901561e-10 2.1030934718853924e-12 "
+            "2.1030934718853924e-12 6.728656898188601e-7 2.0239658022589593e-14"),
+        consistent_str="-2.5 6.899999999999999 4.5 4.5 -26.75",
+        dissipation_pos_str=("-0.17685508040141606 0.523974175351619 "
+            "0.49761166908056337 0.49761166908056337 -2.495196349669178"),
+        dissipation_neg_str=("-1.6768551393982627 4.823974372333487 "
+            "3.9976117689805504 3.9976117689805504 -22.145196359962693"),
+        weno_flux_str=("-4.353710219799678 12.247948547685105 "
+            "8.995223438061114 8.995223438061114 -51.390392709631875"),
+        direction="x")
+single_data["Case b:y"] = FluxDataSingle(
+        states_str="1 -1 -1 -1 5.5,2 -4 -4 -4 20",
+        fluxes_str="-1 2.6 1 1 -7.1,-4 11.2 8 8 -46.4",
+        lam_str=("-1. -1. -1. 0.49666295470957644 -2.4966629547095764,"
+            "-2. -2. -2. -0.5033370452904236 -3.4966629547095764"),
+        R_str=("1 0 0 0.45781245952886396 0.45781245952886396,"
+            "-1.5857864376269053 0 0 -0.018886008110941373 "
+            "-1.4330995704840366,"
+            "-1.5857864376269053 1.4142135623730951 0 -0.7259927892974889 "
+            "-0.7259927892974889,"
+            "-1.5857864376269053 0 1.4142135623730951 -0.7259927892974889 "
+            "-0.7259927892974889,"
+            "3.7720779386421466 -2.2426406871192857 -2.2426406871192857 "
+            "3.335959603054103 5.578600290173388"),
+        R_inv_str=("0.36752136386566203 -0.265894835574803 -0.265894835574803 "
+            "-0.265894835574803 -0.16767379847989416,"
+            "1.1213203435596428 0 0.7071067811865475 0 0,"
+            "1.1213203435596428 0 0 0.7071067811865475 0,"
+            "1.8120820550093746 0.9975038969055432 0.2903971157189956 "
+            "0.2903971157189956 0.1831249838115456,"
+            "-0.43055863210991147 -0.41670966546755195 0.2903971157189956 "
+            "0.2903971157189956 0.1831249838115456"),
+        wavespeeds_str="2.2 2.2 2.2 0.553670749819466 3.8463292501805344",
+        char_fluxes_pos_str=("0.06722315769446469 0.24852813742385726 "
+            "0.24852813742385726 0.37456222716537935 0.10725061063772967,"
+            "0.16835489671908455 -0.05857864376269045 "
+            "-0.05857864376269045 0.30602629876675014 0.7274934638648571"),
+        char_fluxes_neg_str=("-0.4673767959969717 -0.6627416997969526 "
+            "-0.6627416997969526 -0.3125273038443238 -1.4795302623674356,"
+            "-1.0907156303492833 1.2301515190164993 "
+            "1.2301515190164993 -0.23295627081394965 -7.523052586013577"),
+        oscillation_pos_str=("0. 0. 0. 0. 0.,"
+            "0.01363683818419176 0.12575276673434946 0.12575276673434946 "
+            "0.006262897975282706 0.5129349293057707,"
+            "0.0340920954604794 0.31438191683587363 0.31438191683587363 "
+            "0.015657244938206753 1.2823373232644266"),
+        oscillation_neg_str=("0. 0. 0. 0. 0.,"
+            "0.5180684032155981 4.777392983773268 4.777392983773268 "
+            "0.00844206573002786 48.69888276854567,"
+            "1.2951710080389953 11.94348245943317 11.94348245943317 "
+            "0.02110516432506965 121.7472069213642"),
+        weno_weights_pos_str=("0.9999965203327451 0.999999959029252 "
+            "0.999999959029252 0.9999835300215946 0.9999999975371708,"
+            "3.2217041402391125e-6 3.793560951911702e-8 "
+            "3.793560951911702e-8 0.000015247816238787948 2.2803933579018906e-9,"
+            "2.579631146567429e-7 3.0351383606886035e-9 "
+            "3.0351383606886035e-9 1.2221621667041045e-6 1.824357365679662e-10"),
+        weno_weights_neg_str=("0.999999997585736 0.9999999999716082 "
+            "0.9999999999716082 0.9999909282503908 0.9999999999997268,"
+            "2.2354258683870077e-9 2.6288602366075745e-11 "
+            "2.6288602366075745e-11 8.39888391917283e-6 2.5299566294114733e-13,"
+            "1.788382117901561e-10 2.1030934718853924e-12 "
+            "2.1030934718853924e-12 6.728656898188601e-7 2.0239658022589593e-14"),
+        consistent_str="-2.5 6.899999999999999 4.5 4.5 -26.75",
+        dissipation_pos_str=("-0.17685508040141606 0.523974175351619 "
+            "0.49761166908056337 0.49761166908056337 -2.495196349669178"),
+        dissipation_neg_str=("-1.6768551393982627 4.823974372333487 "
+            "3.9976117689805504 3.9976117689805504 -22.145196359962693"),
+        weno_flux_str=("-4.353710219799678 12.247948547685105 "
+            "8.995223438061114 8.995223438061114 -51.390392709631875"),
+        direction="y")
+single_data["Case b:z"] = FluxDataSingle(
+        states_str="1 -1 -1 -1 5.5,2 -4 -4 -4 20",
+        fluxes_str="-1 2.6 1 1 -7.1,-4 11.2 8 8 -46.4",
+        lam_str=("-1. -1. -1. 0.49666295470957644 -2.4966629547095764,"
+            "-2. -2. -2. -0.5033370452904236 -3.4966629547095764"),
+        R_str=("1 0 0 0.45781245952886396 0.45781245952886396,"
+            "-1.5857864376269053 0 0 -0.018886008110941373 "
+            "-1.4330995704840366,"
+            "-1.5857864376269053 1.4142135623730951 0 -0.7259927892974889 "
+            "-0.7259927892974889,"
+            "-1.5857864376269053 0 1.4142135623730951 -0.7259927892974889 "
+            "-0.7259927892974889,"
+            "3.7720779386421466 -2.2426406871192857 -2.2426406871192857 "
+            "3.335959603054103 5.578600290173388"),
+        R_inv_str=("0.36752136386566203 -0.265894835574803 -0.265894835574803 "
+            "-0.265894835574803 -0.16767379847989416,"
+            "1.1213203435596428 0 0.7071067811865475 0 0,"
+            "1.1213203435596428 0 0 0.7071067811865475 0,"
+            "1.8120820550093746 0.9975038969055432 0.2903971157189956 "
+            "0.2903971157189956 0.1831249838115456,"
+            "-0.43055863210991147 -0.41670966546755195 0.2903971157189956 "
+            "0.2903971157189956 0.1831249838115456"),
+        wavespeeds_str="2.2 2.2 2.2 0.553670749819466 3.8463292501805344",
+        char_fluxes_pos_str=("0.06722315769446469 0.24852813742385726 "
+            "0.24852813742385726 0.37456222716537935 0.10725061063772967,"
+            "0.16835489671908455 -0.05857864376269045 "
+            "-0.05857864376269045 0.30602629876675014 0.7274934638648571"),
+        char_fluxes_neg_str=("-0.4673767959969717 -0.6627416997969526 "
+            "-0.6627416997969526 -0.3125273038443238 -1.4795302623674356,"
+            "-1.0907156303492833 1.2301515190164993 "
+            "1.2301515190164993 -0.23295627081394965 -7.523052586013577"),
+        oscillation_pos_str=("0. 0. 0. 0. 0.,"
+            "0.01363683818419176 0.12575276673434946 0.12575276673434946 "
+            "0.006262897975282706 0.5129349293057707,"
+            "0.0340920954604794 0.31438191683587363 0.31438191683587363 "
+            "0.015657244938206753 1.2823373232644266"),
+        oscillation_neg_str=("0. 0. 0. 0. 0.,"
+            "0.5180684032155981 4.777392983773268 4.777392983773268 "
+            "0.00844206573002786 48.69888276854567,"
+            "1.2951710080389953 11.94348245943317 11.94348245943317 "
+            "0.02110516432506965 121.7472069213642"),
+        weno_weights_pos_str=("0.9999965203327451 0.999999959029252 "
+            "0.999999959029252 0.9999835300215946 0.9999999975371708,"
+            "3.2217041402391125e-6 3.793560951911702e-8 "
+            "3.793560951911702e-8 0.000015247816238787948 2.2803933579018906e-9,"
+            "2.579631146567429e-7 3.0351383606886035e-9 "
+            "3.0351383606886035e-9 1.2221621667041045e-6 1.824357365679662e-10"),
+        weno_weights_neg_str=("0.999999997585736 0.9999999999716082 "
+            "0.9999999999716082 0.9999909282503908 0.9999999999997268,"
+            "2.2354258683870077e-9 2.6288602366075745e-11 "
+            "2.6288602366075745e-11 8.39888391917283e-6 2.5299566294114733e-13,"
+            "1.788382117901561e-10 2.1030934718853924e-12 "
+            "2.1030934718853924e-12 6.728656898188601e-7 2.0239658022589593e-14"),
+        consistent_str="-2.5 6.899999999999999 4.5 4.5 -26.75",
+        dissipation_pos_str=("-0.17685508040141606 0.523974175351619 "
+            "0.49761166908056337 0.49761166908056337 -2.495196349669178"),
+        dissipation_neg_str=("-1.6768551393982627 4.823974372333487 "
+            "3.9976117689805504 3.9976117689805504 -22.145196359962693"),
+        weno_flux_str=("-4.353710219799678 12.247948547685105 "
+            "8.995223438061114 8.995223438061114 -51.390392709631875"),
+        direction="z")
+
+# }}}
+
+# {{{ Input data: Case (c)
+
+single_data["Case c:x"] = FluxDataSingle(
+        states_str="2 4 8 12 64,1 1 2 3 11",
+        fluxes_str="4 11.2 16 24 134.4,1 2.6 2 3 12.6",
+        lam_str=("2. 2. 2. 3.4966629547095764 0.5033370452904236,"
+            "1. 1. 1. 2.4966629547095764 -0.49666295470957644"),
+        R_str=("1 0 0 0.41384589551844686 0.41384589551844686,"
+            "1.5857864376269053 0 0 1.363377989567262 -0.05083557280583326,"
+            "3.1715728752538106 1.4142135623730951 0 1.3125424167614286 "
+            "1.3125424167614286,"
+            "4.757359312880716 0 1.4142135623730951 1.968813625142143 "
+            "1.968813625142143,"
+            "17.60303038033002 4.485281374238571 6.727922061357858 "
+            "11.42671019719969 9.184069510080404"),
+        R_inv_str=("-1.4118746341171051 0.21727611674823194 0.4345522334964639 "
+            "0.6518283502446959 -0.13701474018997217,"
+            "-2.2426406871192857 0 0.7071067811865475 0 0,"
+            "-3.3639610306789285 0 0 0.7071067811865475 0,"
+            "1.7926564050747988 0.4445982978342618 -0.5250169667045714 "
+            "-0.7875254500568571 0.16553835820737872,"
+            "4.035297092194084 -0.9696152645388333 -0.5250169667045714 "
+            "-0.7875254500568571 0.16553835820737872"),
+        wavespeeds_str="2.2 2.2 2.2 3.8463292501805344 0.553670749819466",
+        char_fluxes_pos_str=("1.1162114029572359 2.4603030380329987 "
+            "3.690454557049497 7.986924884679773 0.4290108979594782,"
+            "0.26073527447612554 -1.3254833995939053 "
+            "-1.9882250993908572 2.052422890589809 0.5017887559695788"),
+        char_fluxes_neg_str=("-0.1482823715615963 -0.1171572875253809 "
+            "-0.1757359312880714 -0.8893252542028551 -0.20004041758408686,"
+            "-0.009488213714461069 0.49705627484771453 "
+            "0.7455844122715716 -0.4306378813126712 -0.31431832174320373"),
+        oscillation_pos_str=("0. 0. 0. 0. 0.,0.9757858752013722 "
+            "19.109571935093072 42.996536853959384 46.95775189047701 "
+            "0.007062155488717821,2.4394646880034303 47.77392983773268 "
+            "107.49134213489849 117.39437972619254 0.017655388721794552"),
+        oscillation_neg_str=("0. 0. 0. 0. 0.,0.025685091003327314 "
+            "0.5030110669373978 1.1317749006091449 0.28052547473186484 "
+            "0.017412585838667068,0.06421272750831829 1.2575276673434945 "
+            "2.8294372515228616 0.7013136868296621 0.04353146459666767"),
+        weno_weights_pos_str=("0.999999999319454 0.9999999999982254 "
+            "0.9999999999996494 0.9999999999997061 0.9999870425235151,"
+            "6.301345524776077e-10 1.6430428067143077e-12 "
+            "3.245518542322869e-13 2.721049684463335e-13 0.000011996153712066484,"
+            "5.041138413805998e-11 1.3144350707803127e-13 "
+            "2.5964155584975223e-14 2.1768403038595738e-14 9.613227729623656e-7"),
+        weno_weights_neg_str=("0.9999990185022956 0.9999999974390363 "
+            "0.9999999994941201 0.9999999917661908 0.9999978651320311,"
+            "9.08762722250494e-7 2.3712585027633853e-9 "
+            "4.684070887233427e-10 7.62387333908512e-9 1.976628702372801e-6,"
+            "7.27349821618269e-8 1.897052057748541e-10 "
+            "3.747296441221644e-11 6.099359570650648e-10 1.5823926647977797e-7"),
+        consistent_str="2.5 6.9 9. 13.5 73.5",
+        dissipation_pos_str=("1.6406634409601506 4.725635754575325 "
+            "7.88043892893644 11.820658393408754 68.69422377699841"),
+        dissipation_neg_str=("0.14066328824586258 0.4256356884430299 "
+            "0.8804384437989111 1.3206576666570364 7.7942206041633595"),
+        weno_flux_str=("4.281326729206013 12.051271443018354 "
+            "17.76087737273535 26.64131606006579 149.98844438116177"),
+        direction="x")
+single_data["Case c:y"] = FluxDataSingle(
+        states_str="2 8 12 4 64,1 2 3 1 11",
+        fluxes_str="8 35.2 48 16 268.8,2 5.6 6 2 25.2",
+        lam_str=("4. 4. 4. 5.496662954709576 2.5033370452904236,"
+            "2. 2. 2. 3.4966629547095764 0.5033370452904236"),
+        R_str=("1 0 0 0.41384589551844686 0.41384589551844686,"
+            "3.1715728752538106 0 0 2.019649197947976 0.605435635574881,"
+            "4.757359312880716 1.4142135623730951 0 1.968813625142143 "
+            "1.968813625142143,"
+            "1.5857864376269053 0 1.4142135623730951 0.6562712083807143 "
+            "0.6562712083807143,"
+            "17.60303038033002 6.727922061357858 2.2426406871192857 "
+            "12.548030540759331 8.062749166520762"),
+        R_inv_str=("-1.4118746341171051 0.4345522334964639 0.6518283502446959 "
+            "0.21727611674823194 -0.13701474018997217,"
+            "-3.3639610306789285 0 0.7071067811865475 0 0,"
+            "-1.1213203435596428 0 0 0.7071067811865475 0,"
+            "0.671336061515156 0.18208981448197611 -0.7875254500568571 "
+            "-0.2625084833522857 0.16553835820737872,"
+            "5.156617435753728 -1.2321237478911191 -0.7875254500568571 "
+            "-0.2625084833522857 0.16553835820737872"),
+        wavespeeds_str="4.4 4.4 4.4 6.046329250180534 2.753670749819466",
+        char_fluxes_pos_str=("2.2324228059144673 7.380909114098994 "
+            "2.4603030380329987 15.885347333048884 0.9465242322295992,"
+            "0.5214705489522515 -3.9764501987817145 "
+            "-1.3254833995939053 1.3413036144786457 3.767119678640133"),
+        char_fluxes_neg_str=("-0.2965647431231918 -0.3514718625761428 "
+            "-0.1171572875253809 -1.6097440213843948 -0.5689873221894901,"
+            "-0.018976427428922138 1.4911688245431431 "
+            "0.49705627484771453 -0.05753157056903602 -1.432380835542713"),
+        oscillation_pos_str=("0. 0. 0. 0. 0.,3.9031435008054665 "
+            "171.98614741583754 19.109571935093072 282.0389435835765 "
+            "10.607678229749117,9.757858752013666 429.96536853959395 "
+            "47.77392983773268 705.0973589589414 26.519195574372795"),
+        oscillation_neg_str=("0. 0. 0. 0. 0.,0.10274036401330869 "
+            "4.5270996024365795 0.5030110669373978 3.21248465662163 "
+            "0.9939311452005626,0.25685091003327176 11.317749006091447 "
+            "1.2575276673434945 8.031211641554076 2.484827863001407"),
+        weno_weights_pos_str=("0.9999999999574652 0.999999999999978 "
+            "0.9999999999982254 0.9999999999999919 0.9999999999942413,"
+            "3.9384014966385437e-11 2.0284497966079935e-14 "
+            "1.6430428067143077e-12 7.542808138536806e-15 5.332240836330361e-12,"
+            "3.1507308840273685e-12 1.622759950511315e-15 "
+            "1.3144350707803127e-13 6.034246767570432e-16 4.265797494767529e-13"),
+        weno_weights_neg_str=("0.9999999386221037 0.9999999999683822 "
+            "0.9999999974390363 0.9999999999372101 0.9999999993440752,"
+            "5.68308936788955e-8 2.9275831062158654e-11 "
+            "2.3712585027633853e-9 5.8138848037057226e-11 6.073372407373299e-10,"
+            "4.547002513715645e-9 2.3420726930957915e-12 "
+            "1.897052057748541e-10 4.6511252168302e-12 4.8587565862160355e-11"),
+        consistent_str="5. 20.4 27. 9. 147.",
+        dissipation_pos_str=("3.281326602835503 16.54629350160538 "
+            "23.641315459112736 7.8804384863675505 137.3884413587739"),
+        dissipation_neg_str=("0.2813265968286516 1.7462934823903076 "
+            "2.641315430506613 0.8804384760489334 15.588441251607525"),
+        weno_flux_str=("8.562653199664155 38.69258698399568 "
+            "53.28263088961935 17.760876962416482 299.9768826103814"),
+        direction="y")
+single_data["Case c:z"] = FluxDataSingle(
+        states_str="2 12 4 8 64,1 3 1 2 11",
+        fluxes_str="12 75.2 24 48 403.2,3 10.6 3 6 37.8",
+        lam_str=("6. 6. 6. 7.496662954709576 4.503337045290424,"
+            "3. 3. 3. 4.496662954709576 1.5033370452904236"),
+        R_str=("1 0 0 0.41384589551844686 0.41384589551844686,"
+            "4.757359312880716 0 0 2.6759204063286903 1.2617068439555954,"
+            "1.5857864376269053 1.4142135623730951 0 0.6562712083807143 "
+            "0.6562712083807143,"
+            "3.1715728752538106 0 1.4142135623730951 1.3125424167614286 "
+            "1.3125424167614286,"
+            "17.60303038033002 2.2426406871192857 4.485281374238571 "
+            "13.669350884318975 6.941428822961119"),
+        R_inv_str=("-1.4118746341171051 0.6518283502446959 "
+            "0.21727611674823194 0.4345522334964639 -0.13701474018997217,"
+            "-1.1213203435596428 0 0.7071067811865475 0 0,"
+            "-2.2426406871192857 0 0 0.7071067811865475 0,"
+            "-0.4499842820444868 -0.08041866887030964 -0.2625084833522857 "
+            "-0.5250169667045714 0.16553835820737872,"
+            "6.27793777931337 -1.4946322312434048 -0.2625084833522857 "
+            "-0.5250169667045714 0.16553835820737872"),
+        wavespeeds_str="6.6 6.6 6.6 8.246329250180533 4.953670749819467",
+        char_fluxes_pos_str=("3.348634208871708 3.690454557049497 "
+            "7.380909114098994 26.24407281945101 -0.9962654715332988,"
+            "0.7822058234283737 -1.9882250993908581 "
+            "-3.9764501987817162 -0.6952990612264269 8.357934000904585"),
+        char_fluxes_neg_str=("-0.44484711468478866 -0.1757359312880714 "
+            "-0.3514718625761428 -2.447320076091316 -0.8207769392695052,"
+            "-0.02846464114338254 0.7455844122715716 "
+            "1.4911688245431431 0.8126310150223119 -3.0474996241899346"),
+        oscillation_pos_str=("0. 0. 0. 0. 0.,8.782072876812371 "
+            "42.9965368539594 171.9861474158376 967.6396764339121 "
+            "116.6680636935429,21.955182192030932 107.49134213489847 "
+            "429.9653685395939 2419.099191084781 291.67015923385725"),
+        oscillation_neg_str=("0. 0. 0. 0. 0.,0.23116581902994632 "
+            "1.1317749006091449 4.5270996024365795 14.169708155270577 "
+            "6.6110585540523275,0.5779145475748658 2.8294372515228616 "
+            "11.317749006091447 35.42427038817644 16.527646385130815"),
+        weno_weights_pos_str=("0.999999999991598 0.9999999999996494 "
+            "0.999999999999978 0.9999999999999993 0.9999999999999524,"
+            "7.779580658268568e-12 3.2455185423228674e-13 "
+            "2.028449796607992e-14 6.408020704124379e-16 4.408056940300416e-14,"
+            "6.223673030753551e-13 2.596415558497523e-14 "
+            "1.6227599505113154e-15 5.126416626873783e-17 3.526445914956101e-15"),
+        weno_weights_neg_str=("0.9999999878747177 0.9999999994941201 "
+            "0.9999999999683822 0.9999999999967727 0.9999999999851737,"
+            "1.1227070122726888e-8 4.684070887233427e-10 "
+            "2.9275831062158654e-11 2.988331869942141e-12 1.372802081812638e-11,"
+            "8.982122341016933e-10 3.747296441221644e-11 "
+            "2.3420726930957915e-12 2.390667520553204e-13 1.0982436589127136e-12"),
+        consistent_str="7.5 42.9 13.5 27. 220.5",
+        dissipation_pos_str=("4.921989904281101 36.24738971766925 "
+            "11.82065772959958 23.641315459201046 206.08266203865145"),
+        dissipation_neg_str=("0.4219899024845117 3.9473897091113366 "
+            "1.3206577265155963 2.6413154534736667 23.382662006532225"),
+        weno_flux_str=("12.843979806765613 83.09477942678058 "
+            "26.641315456115176 53.28263091267471 449.9653240451837"),
+        direction="z")
+
+# }}}
+
+# {{{ Input data: Case (d)
+
+single_data["Case d:x"] = FluxDataSingle(
+        states_str="1 -1 -2 -3 11,2 -4 -8 -12 64",
+        fluxes_str="-1 2.6 2 3 -12.6,-4 11.2 16 24 -134.4",
+        lam_str=("-1. -1. -1. 0.49666295470957644 -2.4966629547095764,"
+            "-2. -2. -2. -0.5033370452904236 -3.4966629547095764"),
+        R_str=("1 0 0 0.41384589551844686 0.41384589551844686,"
+            "-1.5857864376269053 0 0 0.05083557280583326 -1.363377989567262,"
+            "-3.1715728752538106 1.4142135623730951 0 -1.3125424167614286 "
+            "-1.3125424167614286,"
+            "-4.757359312880716 0 1.4142135623730951 -1.968813625142143 "
+            "-1.968813625142143,"
+            "17.60303038033002 -4.485281374238571 -6.727922061357858 "
+            "9.184069510080404 11.42671019719969"),
+        R_inv_str=("-1.4118746341171051 -0.21727611674823194 "
+            "-0.4345522334964639 -0.6518283502446959 -0.13701474018997217,"
+            "2.2426406871192857 0 0.7071067811865475 0 0,"
+            "3.3639610306789285 0 0 0.7071067811865475 0,"
+            "4.035297092194084 0.9696152645388333 0.5250169667045714 "
+            "0.7875254500568571 0.16553835820737872,"
+            "1.7926564050747988 -0.4445982978342618 0.5250169667045714 "
+            "0.7875254500568571 0.16553835820737872"),
+        wavespeeds_str="2.2 2.2 2.2 0.553670749819466 3.8463292501805344",
+        char_fluxes_pos_str=("0.009488213714461069 0.49705627484771453 "
+            "0.7455844122715716 0.31431832174320373 0.4306378813126712,"
+            "0.1482823715615963 -0.1171572875253809 "
+            "-0.1757359312880714 0.20004041758408686 0.8893252542028551"),
+        char_fluxes_neg_str=("-0.26073527447612554 -1.3254833995939053 "
+            "-1.9882250993908572 -0.5017887559695788 -2.052422890589809,"
+            "-1.1162114029572359 2.4603030380329987 "
+            "3.690454557049497 -0.4290108979594782 -7.986924884679773"),
+        oscillation_pos_str=("0. 0. 0. 0. 0.,0.025685091003327314 "
+            "0.5030110669373978 1.1317749006091449 0.017412585838667068 "
+            "0.28052547473186484,0.06421272750831829 1.2575276673434945 "
+            "2.8294372515228616 0.04353146459666767 0.7013136868296621"),
+        oscillation_neg_str=("0. 0. 0. 0. 0.,0.9757858752013722 "
+            "19.109571935093072 42.996536853959384 0.007062155488717821 "
+            "46.95775189047701,2.4394646880034303 47.77392983773268 "
+            "107.49134213489849 0.017655388721794552 117.39437972619254"),
+        weno_weights_pos_str=("0.9999990185022956 0.9999999974390363 "
+            "0.9999999994941201 0.9999978651320311 0.9999999917661908,"
+            "9.08762722250494e-7 2.3712585027633853e-9 "
+            "4.684070887233427e-10 1.976628702372801e-6 7.62387333908512e-9,"
+            "7.27349821618269e-8 1.897052057748541e-10 "
+            "3.747296441221644e-11 1.5823926647977797e-7 6.099359570650648e-10"),
+        weno_weights_neg_str=("0.999999999319454 0.9999999999982254 "
+            "0.9999999999996494 0.9999870425235151 0.9999999999997061,"
+            "6.301345524776077e-10 1.6430428067143077e-12 "
+            "3.245518542322869e-13 0.000011996153712066484 2.721049684463335e-13,"
+            "5.041138413805998e-11 1.3144350707803127e-13 "
+            "2.5964155584975223e-14 9.613227729623656e-7 2.1768403038595738e-14"),
+        consistent_str="-2.5 6.9 9. 13.5 -73.5",
+        dissipation_pos_str=("-0.14066328824586258 0.42563568844302985 "
+            "0.8804384437989112 1.3206576666570364 -7.7942206041633595"),
+        dissipation_neg_str=("-1.6406634409601506 4.725635754575325 "
+            "7.88043892893644 11.820658393408754 -68.69422377699841"),
+        weno_flux_str=("-4.281326729206013 12.051271443018354 "
+            "17.76087737273535 26.641316060065794 -149.98844438116177"),
+        direction="x")
+single_data["Case d:y"] = FluxDataSingle(
+        states_str="1 -2 -3 -1 11,2 -8 -12 -4 64",
+        fluxes_str="-2 5.6 6 2 -25.2,-8 35.2 48 16 -268.8",
+        lam_str=("-2. -2. -2. -0.5033370452904236 -3.4966629547095764,"
+            "-4. -4. -4. -2.5033370452904236 -5.496662954709576"),
+        R_str=("1 0 0 0.41384589551844686 0.41384589551844686,"
+            "-3.1715728752538106 0 0 -0.605435635574881 -2.019649197947976,"
+            "-4.757359312880716 1.4142135623730951 0 -1.968813625142143 "
+            "-1.968813625142143,"
+            "-1.5857864376269053 0 1.4142135623730951 -0.6562712083807143 "
+            "-0.6562712083807143,"
+            "17.60303038033002 -6.727922061357858 -2.2426406871192857 "
+            "8.062749166520762 12.548030540759331"),
+        R_inv_str=("-1.4118746341171051 -0.4345522334964639 "
+            "-0.6518283502446959 -0.21727611674823194 -0.13701474018997217,"
+            "3.3639610306789285 0 0.7071067811865475 0 0,"
+            "1.1213203435596428 0 0 0.7071067811865475 0,"
+            "5.156617435753728 1.2321237478911191 0.7875254500568571 "
+            "0.2625084833522857 0.16553835820737872,"
+            "0.671336061515156 -0.18208981448197611 0.7875254500568571 "
+            "0.2625084833522857 0.16553835820737872"),
+        wavespeeds_str="4.4 4.4 4.4 2.753670749819466 6.046329250180534",
+        char_fluxes_pos_str=("0.018976427428922138 1.4911688245431431 "
+            "0.49705627484771453 1.432380835542713 0.05753157056903602,"
+            "0.2965647431231918 -0.3514718625761428 "
+            "-0.1171572875253809 0.5689873221894901 1.6097440213843948"),
+        char_fluxes_neg_str=("-0.5214705489522515 -3.9764501987817145 "
+            "-1.3254833995939053 -3.767119678640133 -1.3413036144786457,"
+            "-2.2324228059144673 7.380909114098994 "
+            "2.4603030380329987 -0.9465242322295992 -15.885347333048884"),
+        oscillation_pos_str=("0. 0. 0. 0. 0.,0.10274036401330869 "
+            "4.5270996024365795 0.5030110669373978 0.9939311452005626 "
+            "3.21248465662163,0.25685091003327176 11.317749006091447 "
+            "1.2575276673434945 2.484827863001407 8.031211641554076"),
+        oscillation_neg_str=("0. 0. 0. 0. 0.,3.9031435008054665 "
+            "171.98614741583754 19.109571935093072 10.607678229749117 "
+            "282.0389435835765,9.757858752013666 429.96536853959395 "
+            "47.77392983773268 26.519195574372795 705.0973589589414"),
+        weno_weights_pos_str=("0.9999999386221037 0.9999999999683822 "
+            "0.9999999974390363 0.9999999993440752 0.9999999999372101,"
+            "5.68308936788955e-8 2.9275831062158654e-11 "
+            "2.3712585027633853e-9 6.073372407373299e-10 5.8138848037057226e-11,"
+            "4.547002513715645e-9 2.3420726930957915e-12 "
+            "1.897052057748541e-10 4.8587565862160355e-11 4.6511252168302e-12"),
+        weno_weights_neg_str=("0.9999999999574652 0.999999999999978 "
+            "0.9999999999982254 0.9999999999942413 0.9999999999999919,"
+            "3.9384014966385437e-11 2.0284497966079935e-14 "
+            "1.6430428067143077e-12 5.332240836330361e-12 7.542808138536806e-15,"
+            "3.1507308840273685e-12 1.622759950511315e-15 "
+            "1.3144350707803127e-13 4.265797494767529e-13 6.034246767570432e-16"),
+        consistent_str="-5. 20.4 27. 9. -147.",
+        dissipation_pos_str=("-0.2813265968286516 1.7462934823903074 "
+            "2.641315430506612 0.8804384760489334 -15.588441251607525"),
+        dissipation_neg_str=("-3.281326602835503 16.54629350160538 "
+            "23.641315459112736 7.880438486367551 -137.3884413587739"),
+        weno_flux_str=("-8.562653199664155 38.69258698399569 "
+            "53.28263088961934 17.760876962416486 -299.97688261038144"),
+        direction="y")
+single_data["Case d:z"] = FluxDataSingle(
+        states_str="1 -3 -1 -2 11,2 -12 -4 -8 64",
+        fluxes_str="-3 10.6 3 6 -37.8,-12 75.2 24 48 -403.2",
+        lam_str=("-3. -3. -3. -1.5033370452904236 -4.496662954709576,"
+            "-6. -6. -6. -4.503337045290424 -7.496662954709576"),
+        R_str=("1 0 0 0.41384589551844686 0.41384589551844686,"
+            "-4.757359312880716 0 0 -1.2617068439555954 -2.6759204063286903,"
+            "-1.5857864376269053 1.4142135623730951 0 -0.6562712083807143 "
+            "-0.6562712083807143,"
+            "-3.1715728752538106 0 1.4142135623730951 -1.3125424167614286 "
+            "-1.3125424167614286,"
+            "17.60303038033002 -2.2426406871192857 -4.485281374238571 "
+            "6.941428822961119 13.669350884318975"),
+        R_inv_str=("-1.4118746341171051 -0.6518283502446959 "
+            "-0.21727611674823194 -0.4345522334964639 -0.13701474018997217,"
+            "1.1213203435596428 0 0.7071067811865475 0 0,"
+            "2.2426406871192857 0 0 0.7071067811865475 0,"
+            "6.27793777931337 1.4946322312434048 0.2625084833522857 "
+            "0.5250169667045714 0.16553835820737872,"
+            "-0.4499842820444868 0.08041866887030964 0.2625084833522857 "
+            "0.5250169667045714 0.16553835820737872"),
+        wavespeeds_str="6.6 6.6 6.6 4.953670749819467 8.246329250180533",
+        char_fluxes_pos_str=("0.02846464114338254 0.7455844122715716 "
+            "1.4911688245431431 3.0474996241899346 -0.8126310150223119,"
+            "0.44484711468478866 -0.1757359312880714 "
+            "-0.3514718625761428 0.8207769392695052 2.447320076091316"),
+        char_fluxes_neg_str=("-0.7822058234283737 -1.9882250993908581 "
+            "-3.9764501987817162 -8.357934000904585 0.6952990612264269,"
+            "-3.348634208871708 3.690454557049497 "
+            "7.380909114098994 0.9962654715332988 -26.24407281945101"),
+        oscillation_pos_str=("0. 0. 0. 0. 0.,0.23116581902994632 "
+            "1.1317749006091449 4.5270996024365795 6.6110585540523275 "
+            "14.169708155270577,0.5779145475748658 2.8294372515228616 "
+            "11.317749006091447 16.527646385130815 35.42427038817644"),
+        oscillation_neg_str=("0. 0. 0. 0. 0.,8.782072876812371 "
+            "42.9965368539594 171.9861474158376 116.6680636935429 "
+            "967.6396764339121,21.955182192030932 107.49134213489847 "
+            "429.9653685395939 291.67015923385725 2419.099191084781"),
+        weno_weights_pos_str=("0.9999999878747177 0.9999999994941201 "
+            "0.9999999999683822 0.9999999999851737 0.9999999999967727,"
+            "1.1227070122726888e-8 4.684070887233427e-10 "
+            "2.9275831062158654e-11 1.372802081812638e-11 2.988331869942141e-12,"
+            "8.982122341016933e-10 3.747296441221644e-11 "
+            "2.3420726930957915e-12 1.0982436589127136e-12 2.390667520553204e-13"),
+        weno_weights_neg_str=("0.999999999991598 0.9999999999996494 "
+            "0.999999999999978 0.9999999999999524 0.9999999999999993,"
+            "7.779580658268568e-12 3.2455185423228674e-13 "
+            "2.028449796607992e-14 4.408056940300416e-14 6.408020704124379e-16,"
+            "6.223673030753551e-13 2.596415558497523e-14 "
+            "1.6227599505113154e-15 3.526445914956101e-15 5.126416626873783e-17"),
+        consistent_str="-7.5 42.9 13.5 27. -220.5",
+        dissipation_pos_str=("-0.42198990248451174 3.947389709111337 "
+            "1.3206577265155963 2.6413154534736667 -23.382662006532225"),
+        dissipation_neg_str=("-4.9219899042811 36.247389717669265 "
+            "11.820657729599578 23.641315459201046 -206.08266203865145"),
+        weno_flux_str=("-12.843979806765612 83.09477942678059 "
+            "26.641315456115173 53.28263091267471 -449.9653240451837"),
+        direction="z")
+
+# }}}
+
+
+@pytest.fixture(scope="session", params=[
+    "Case flat:x", "Case flat:y", "Case flat:z",
+    "Case a:x", "Case a:y", "Case a:z",
+    "Case b:x", "Case b:y", "Case b:z",
+    "Case c:x", "Case c:y", "Case c:z",
+    "Case d:x", "Case d:y", "Case d:z"])
+def flux_test_data_fixture(request):
+    return single_data[request.param]
+
+
+# vim: foldmethod=marker
diff --git a/requirements.txt b/requirements.txt
index 36d95d70721dcf19677fc5fd2fad9cee4d514ab4..b6462ed615a84e3f73737d7eb6572f7819eb5380 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -5,7 +5,7 @@ git+https://github.com/inducer/pyopencl.git
 git+https://github.com/inducer/pymbolic.git
 git+https://github.com/inducer/genpy.git
 git+https://github.com/inducer/codepy.git
-git+https://gitlab.tiker.net/inducer/loopy.git@kernel_callables_v3
+git+https://gitlab.tiker.net/inducer/loopy.git@fix-slice-processing-at-base
 
 git+https://github.com/inducer/f2py
 
diff --git a/setup.cfg b/setup.cfg
index 6fdc175dd5bae21563adf96162fd2a2ad06e992e..7963e687c9e69eeeeafbee1326357c0b5b7f6d04 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,3 +1,4 @@
 [flake8]
 ignore = E126,E127,E128,E123,E226,E241,E242,E265,N802,N803,N806,W503,E402,N814,W504
 max-line-length=85
+per-file-ignores = test.py:F811
diff --git a/test.py b/test.py
index 878e48c8519491c5462f3aab80328283fb844643..671b7fe606f27d44b93dfd6f761ad7ab6a5878c6 100644
--- a/test.py
+++ b/test.py
@@ -15,339 +15,188 @@ from pyopencl.tools import (  # noqa
         as pytest_generate_tests)
 
 import utilities as u
+from data_for_test import (  # noqa: F401
+        flux_test_data_fixture,
+        single_data as std  # "single_test_data", sorry
+        )
 
 
-@pytest.mark.xfail
-@pytest.mark.parametrize(("gen_fluxes_str,char_fluxes_pos_str,char_fluxes_neg_str,"
-    "combined_metric,R_str,flux_expected_str"), [
-        ("4 11.2 8 8 46.4,1 2.6 1 1 7.1",
-            ("1.09071563 1.23015152 1.23015152 7.52305259 0.232956271,"
-                "0.467376796 -0.6627417 -0.6627417 1.47953026 0.312527304"),
-            ("-0.168354897 -0.0585786438 -0.0585786438 -0.727493464 -0.306026299,"
-                "-0.0672231577 0.248528137 0.248528137 -0.107250611 -0.374562227"),
-            1.0,
-            ("1 0 0 0.45781246 0.45781246,"
-                "1.58578644 0 0 1.43309957 0.0188860081,"
-                "1.58578644 1.41421356 0 0.725992789 0.725992789,"
-                "1.58578644 0 1.41421356 0.725992789 0.725992789,"
-                "3.77207794 2.24264069 2.24264069 5.57860029 3.3359596"),
-            "4.35371022 12.2479485 8.99522344 8.99522344 51.3903927"),
-        ("-1 2.6 1 1 -7.1,-4 11.2 8 8 -46.4",
-            ("0.0672231577 0.248528137 0.248528137 0.374562227 0.107250611,"
-                "0.168354897 -0.0585786438 -0.0585786438 0.306026299 0.727493464"),
-            ("-0.467376796 -0.6627417 -0.6627417 -0.312527304 -1.47953026,"
-                "-1.09071563 1.23015152 1.23015152 -0.232956271 -7.52305259"),
-            1.0,
-            ("1 0 0 0.45781246 0.45781246,"
-                "-1.58578644 0 0 -0.0188860081 -1.43309957,"
-                "-1.58578644 1.41421356 0 -0.725992789 -0.725992789,"
-                "-1.58578644 0 1.41421356 -0.725992789 -0.725992789,"
-                "3.77207794 -2.24264069 -2.24264069 3.3359596 5.57860029"),
-            "-4.35371022 12.2479485 8.99522344 8.99522344 -51.3903927")
-        ])
-def test_weno_flux_uniform_grid(
-        ctx_factory, gen_fluxes_str, char_fluxes_pos_str, char_fluxes_neg_str,
-        combined_metric, R_str, flux_expected_str):
+def test_weno_flux_uniform_grid(ctx_factory, flux_test_data_fixture):
+    data = flux_test_data_fixture
+
     prg = u.get_weno_program_with_root_kernel("weno_flux")
     queue = u.get_queue(ctx_factory)
 
-    nvars = 5
+    flux_dev = u.empty_array_on_device(queue, data.nvars)
 
-    gen_fluxes = u.expand_to_6(
-            u.transposed_array_from_string(gen_fluxes_str))
-    char_fluxes_pos = u.expand_to_6(
-            u.transposed_array_from_string(char_fluxes_pos_str))
-    char_fluxes_neg = u.expand_to_6(
-            u.transposed_array_from_string(char_fluxes_neg_str))
-    R = u.array_from_string(R_str)
-    flux_dev = u.empty_array_on_device(queue, nvars)
-
-    prg(queue, nvars=nvars,
-            generalized_fluxes=gen_fluxes,
-            characteristic_fluxes_pos=char_fluxes_pos,
-            characteristic_fluxes_neg=char_fluxes_neg,
-            combined_frozen_metrics=combined_metric,
-            R=R,
+    prg(queue, nvars=data.nvars,
+            generalized_fluxes=data.fluxes,
+            characteristic_fluxes_pos=data.char_fluxes_pos,
+            characteristic_fluxes_neg=data.char_fluxes_neg,
+            combined_frozen_metrics=1.0,
+            R=data.R,
             flux=flux_dev)
 
-    flux_expected = u.array_from_string(flux_expected_str)
-    u.compare_arrays(flux_dev.get(), flux_expected)
+    u.compare_arrays(flux_dev.get(), data.weno_flux)
 
 
-@pytest.mark.xfail
-@pytest.mark.parametrize(("gen_fluxes_str,consistent_expected_str"), [
-    ("4 11.2 8 8 46.4,1 2.6 1 1 7.1",
-        "2.5 6.9 4.5 4.5 26.75"),
-    ("-1 2.6 1 1 -7.1,-4 11.2 8 8 -46.4",
-        "-2.5 6.9 4.5 4.5 -26.75")
-    ])
-def test_consistent_part_uniform_grid(
-        ctx_factory, gen_fluxes_str, consistent_expected_str):
+def test_consistent_part_uniform_grid(ctx_factory, flux_test_data_fixture):
+    data = flux_test_data_fixture
+
     prg = u.get_weno_program_with_root_kernel("consistent_part")
     queue = u.get_queue(ctx_factory)
 
-    nvars = 5
-
-    gen_fluxes = u.expand_to_6(
-            u.transposed_array_from_string(gen_fluxes_str))
-    consistent_dev = u.empty_array_on_device(queue, nvars)
+    consistent_dev = u.empty_array_on_device(queue, data.nvars)
 
-    prg(queue, nvars=nvars,
-            generalized_fluxes=gen_fluxes,
+    prg(queue, nvars=data.nvars,
+            generalized_fluxes=data.fluxes,
             consistent=consistent_dev)
 
-    consistent_expected = u.array_from_string(consistent_expected_str)
-    u.compare_arrays(consistent_dev.get(), consistent_expected)
-
-
-@pytest.mark.xfail
-@pytest.mark.parametrize(("char_fluxes_str,combined_metric,"
-    "R_str,dissipation_expected_str"), [
-        (("1.09071563 1.23015152 1.23015152 7.52305259 0.232956271,"
-            "0.467376796 -0.6627417 -0.6627417 1.47953026 0.312527304"),
-            1.0,
-            ("1 0 0 0.45781246 0.45781246,"
-                "1.58578644 0 0 1.43309957 0.0188860081,"
-                "1.58578644 1.41421356 0 0.725992789 0.725992789,"
-                "1.58578644 0 1.41421356 0.725992789 0.725992789,"
-                "3.77207794 2.24264069 2.24264069 5.57860029 3.3359596"),
-            "1.67685514 4.82397437 3.99761177 3.99761177 22.1451964"),
-        (("0.0672231577 0.248528137 0.248528137 0.374562227 0.107250611,"
-            "0.168354897 -0.0585786438 -0.0585786438 0.306026299 0.727493464"),
-            1.0,
-            ("1 0 0 0.45781246 0.45781246,"
-                "-1.58578644 0 0 -0.0188860081 -1.43309957,"
-                "-1.58578644 1.41421356 0 -0.725992789 -0.725992789,"
-                "-1.58578644 0 1.41421356 -0.725992789 -0.725992789,"
-                "3.77207794 -2.24264069 -2.24264069 3.3359596 5.57860029"),
-            "-0.17685508 0.523974175 0.497611669 0.497611669 -2.49519635")
-        ])
-def test_dissipation_part_pos_uniform_grid(
-        ctx_factory, char_fluxes_str, combined_metric, R_str,
-        dissipation_expected_str):
+    u.compare_arrays(consistent_dev.get(), data.consistent)
+
+
+def test_dissipation_part_pos_uniform_grid(ctx_factory, flux_test_data_fixture):
+    data = flux_test_data_fixture
+
     prg = u.get_weno_program_with_root_kernel("dissipation_part_pos")
     queue = u.get_queue(ctx_factory)
 
-    nvars = 5
-
-    char_fluxes = u.expand_to_6(
-            u.transposed_array_from_string(char_fluxes_str))
-    R = u.array_from_string(R_str)
-    dissipation_dev = u.empty_array_on_device(queue, nvars)
+    dissipation_dev = u.empty_array_on_device(queue, data.nvars)
 
-    prg(queue, nvars=nvars,
-            characteristic_fluxes=char_fluxes,
-            combined_frozen_metrics=combined_metric,
-            R=R,
+    prg(queue, nvars=data.nvars,
+            characteristic_fluxes=data.char_fluxes_pos,
+            combined_frozen_metrics=1.0,
+            R=data.R,
             dissipation_pos=dissipation_dev)
 
-    dissipation_expected = u.array_from_string(dissipation_expected_str)
-    u.compare_arrays(dissipation_dev.get(), dissipation_expected)
-
-
-@pytest.mark.xfail
-@pytest.mark.parametrize(("char_fluxes_str,combined_metric,"
-    "R_str,dissipation_expected_str"), [
-        (("-0.168354897 -0.0585786438 -0.0585786438 -0.727493464 -0.306026299,"
-            "-0.0672231577 0.248528137 0.248528137 -0.107250611 -0.374562227"),
-            1.0,
-            ("1 0 0 0.45781246 0.45781246,"
-                "1.58578644 0 0 1.43309957 0.0188860081,"
-                "1.58578644 1.41421356 0 0.725992789 0.725992789,"
-                "1.58578644 0 1.41421356 0.725992789 0.725992789,"
-                "3.77207794 2.24264069 2.24264069 5.57860029 3.3359596"),
-            "0.17685508 0.523974175 0.497611669 0.497611669 2.49519635"),
-        (("-0.467376796 -0.6627417 -0.6627417 -0.312527304 -1.47953026,"
-            "-1.09071563 1.23015152 1.23015152 -0.232956271 -7.52305259"),
-            1.0,
-            ("1 0 0 0.45781246 0.45781246,"
-                "-1.58578644 0 0 -0.0188860081 -1.43309957,"
-                "-1.58578644 1.41421356 0 -0.725992789 -0.725992789,"
-                "-1.58578644 0 1.41421356 -0.725992789 -0.725992789,"
-                "3.77207794 -2.24264069 -2.24264069 3.3359596 5.57860029"),
-            "-1.67685514 4.82397437 3.99761177 3.99761177 -22.1451964")
-        ])
-def test_dissipation_part_neg_uniform_grid(
-        ctx_factory, char_fluxes_str, combined_metric, R_str,
-        dissipation_expected_str):
+    u.compare_arrays(dissipation_dev.get(), data.dissipation_pos)
+
+
+def test_dissipation_part_neg_uniform_grid(ctx_factory, flux_test_data_fixture):
+    data = flux_test_data_fixture
+
     prg = u.get_weno_program_with_root_kernel("dissipation_part_neg")
     queue = u.get_queue(ctx_factory)
 
-    nvars = 5
-
-    char_fluxes = u.expand_to_6(
-            u.transposed_array_from_string(char_fluxes_str))
-    R = u.array_from_string(R_str)
-    dissipation_dev = u.empty_array_on_device(queue, nvars)
+    dissipation_dev = u.empty_array_on_device(queue, data.nvars)
 
-    prg(queue, nvars=nvars,
-            characteristic_fluxes=char_fluxes,
-            combined_frozen_metrics=combined_metric,
-            R=R,
+    prg(queue, nvars=data.nvars,
+            characteristic_fluxes=data.char_fluxes_neg,
+            combined_frozen_metrics=1.0,
+            R=data.R,
             dissipation_neg=dissipation_dev)
 
-    dissipation_expected = u.array_from_string(dissipation_expected_str)
-    u.compare_arrays(dissipation_dev.get(), dissipation_expected)
+    u.compare_arrays(dissipation_dev.get(), data.dissipation_neg)
 
 
-@pytest.mark.slow
-@pytest.mark.parametrize(("states_str,fluxes_str,R_inv_str,wavespeeds_str,"
-        "fluxes_pos_expected_str,fluxes_neg_expected_str"), [
-    ("2 4 4 4 20,1 1 1 1 5.5", "4 11.2 8 8 46.4,1 2.6 1 1 7.1",
-        ("0.367521364 0.265894836 0.265894836 0.265894836 -0.167673798,"
-        "-1.12132034 0 0.707106781 0 0,"
-        "-1.12132034 0 0 0.707106781 0,"
-        "-0.430558632 0.416709665 -0.290397116 -0.290397116 0.183124984,"
-        "1.81208206 -0.997503897 -0.290397116 -0.290397116 0.183124984"),
-        "2.2 2.2 2.2 3.84632925 0.55367075",
-        ("1.09071563 1.23015152 1.23015152 7.52305259 0.232956271,"
-        "0.467376796 -0.6627417 -0.6627417 1.47953026 0.312527304"),
-        ("-0.168354897 -0.0585786438 -0.0585786438 -0.727493464 -0.306026299,"
-        "-0.0672231577 0.248528137 0.248528137 -0.107250611 -0.374562227")),
-    ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "-1 2.6 1 1 -7.1,-4 11.2 8 8 -46.4",
-        ("0.367521364 -0.265894836 -0.265894836 -0.265894836 -0.167673798,"
-        "1.12132034 0 0.707106781 0 0,"
-        "1.12132034 0 0 0.707106781 0,"
-        "1.81208206 0.997503897 0.290397116 0.290397116 0.183124984,"
-        "-0.430558632 -0.416709665 0.290397116 0.290397116 0.183124984"),
-        "2.2 2.2 2.2 0.55367075 3.84632925",
-        ("0.0672231577 0.248528137 0.248528137 0.374562227 0.107250611,"
-        "0.168354897 -0.0585786438 -0.0585786438 0.306026299 0.727493464"),
-        ("-0.467376796 -0.6627417 -0.6627417 -0.312527304 -1.47953026,"
-        "-1.09071563 1.23015152 1.23015152 -0.232956271 -7.52305259")),
-    ("2 4 8 12 64,1 1 2 3 11", "4 11.2 16 24 134.4,1 2.6 2 3 12.6",
-        ("-1.41187463 0.217276117 0.434552233 0.65182835 -0.13701474,"
-        "-2.24264069 0 0.707106781 0 0,"
-        "-3.36396103 0 0 0.707106781 0,"
-        "1.79265641 0.444598298 -0.525016967 -0.78752545 0.165538358,"
-        "4.03529709 -0.969615265 -0.525016967 -0.78752545 0.165538358"),
-        "2.2 2.2 2.2 3.84632925 0.55367075",
-        ("1.1162114 2.46030304 3.69045456 7.98692488 0.429010898,"
-        "0.260735274 -1.3254834 -1.9882251 2.05242289 0.501788756"),
-        ("-0.148282372 -0.117157288 -0.175735931 -0.889325254 -0.200040418,"
-        "-0.00948821371 0.497056275 0.745584412 -0.430637881 -0.314318322")),
-    ])
-def test_flux_splitting_uniform_grid(
-        ctx_factory, states_str, fluxes_str, R_inv_str, wavespeeds_str,
-        fluxes_pos_expected_str, fluxes_neg_expected_str):
-    prg = u.get_weno_program_with_root_kernel("split_characteristic_fluxes")
+def test_weno_weights_pos_uniform_grid(ctx_factory, flux_test_data_fixture):
+    data = flux_test_data_fixture
+
+    prg = u.get_weno_program_with_root_kernel("weno_weights_pos")
     queue = u.get_queue(ctx_factory)
 
-    nvars = 5
+    weights_dev = u.empty_array_on_device(queue, data.nvars, 3)
 
-    states = u.expand_to_6(u.transposed_array_from_string(states_str))
-    fluxes = u.expand_to_6(u.transposed_array_from_string(fluxes_str))
-    R_inv = u.array_from_string(R_inv_str)
-    wavespeeds = u.array_from_string(wavespeeds_str)
-    fluxes_pos_dev = u.empty_array_on_device(queue, nvars, 6)
-    fluxes_neg_dev = u.empty_array_on_device(queue, nvars, 6)
-
-    prg(queue, nvars=nvars,
-            generalized_states_frozen=states,
-            generalized_fluxes_frozen=fluxes,
-            R_inv=R_inv,
-            wavespeeds=wavespeeds,
-            characteristic_fluxes_pos=fluxes_pos_dev,
-            characteristic_fluxes_neg=fluxes_neg_dev)
+    prg(queue, nvars=data.nvars,
+            characteristic_fluxes=data.char_fluxes_pos,
+            combined_frozen_metrics=1.0,
+            w=weights_dev)
 
-    fluxes_pos_expected = u.expand_to_6(
-        u.transposed_array_from_string(fluxes_pos_expected_str))
-    u.compare_arrays(fluxes_pos_dev.get(), fluxes_pos_expected)
+    sum_weights = np.sum(weights_dev.get(), axis=1)
+    u.compare_arrays(sum_weights, np.ones(data.nvars))
 
-    fluxes_neg_expected = u.expand_to_6(
-        u.transposed_array_from_string(fluxes_neg_expected_str))
-    u.compare_arrays(fluxes_neg_dev.get(), fluxes_neg_expected)
+    u.compare_arrays(weights_dev.get(), data.weno_weights_pos)
 
 
-@pytest.mark.slow
-@pytest.mark.parametrize("lam_pointwise_str,lam_roe_str,lam_expected_str", [
-    ("1 2 3 4 5,2 4 6 8 10", "1.5 3 4.5 6 7.5", "2.2 4.4 6.6 8.8 11"),
-    ("1 2 3 4 5,-2 -4 -6 -8 -10", "1.5 3 4.5 6 7.5", "2.2 4.4 6.6 8.8 11"),
-    ("1 2 3 4 5,-2 -4 -6 -8 -10", "3 6 9 12 15", "3.3 6.6 9.9 13.2 16.5"),
-    ("1 2 3 4 5,2 4 6 8 10", "-3 -6 -9 -12 -15", "3.3 6.6 9.9 13.2 16.5"),
-    ("3 2 9 4 5,2 6 6 12 10", "-1 -4 -3 -8 -15", "3.3 6.6 9.9 13.2 16.5")
-    ])
-def test_lax_wavespeeds(
-        ctx_factory, lam_pointwise_str, lam_roe_str, lam_expected_str):
-    prg = u.get_weno_program_with_root_kernel("lax_wavespeeds")
+def test_weno_weights_neg_uniform_grid(ctx_factory, flux_test_data_fixture):
+    data = flux_test_data_fixture
+
+    prg = u.get_weno_program_with_root_kernel("weno_weights_neg")
     queue = u.get_queue(ctx_factory)
 
-    nvars = 5
+    weights_dev = u.empty_array_on_device(queue, data.nvars, 3)
 
-    lam_pointwise = u.expand_to_6(u.transposed_array_from_string(lam_pointwise_str))
-    lam_roe = u.array_from_string(lam_roe_str)
-    lam_dev = u.empty_array_on_device(queue, nvars)
+    prg(queue, nvars=data.nvars,
+            characteristic_fluxes=data.char_fluxes_neg,
+            combined_frozen_metrics=1.0,
+            w=weights_dev)
 
-    prg(queue, nvars=nvars, lambda_pointwise=lam_pointwise,
-            lambda_roe=lam_roe, wavespeeds=lam_dev)
+    sum_weights = np.sum(weights_dev.get(), axis=1)
+    u.compare_arrays(sum_weights, np.ones(data.nvars))
 
-    lam_expected = u.array_from_string(lam_expected_str)
-    u.compare_arrays(lam_dev.get(), lam_expected)
+    u.compare_arrays(weights_dev.get(), data.weno_weights_neg)
 
 
-@pytest.mark.slow
-@pytest.mark.parametrize("states_str,direction,lam_expected_str", [
-    ("2 4 4 4 20,1 1 1 1 5.5", "x",
-        "2 2 2 3.49666295 0.503337045,1 1 1 2.49666295 -0.496662955"),
-    ("2 4 4 4 20,1 1 1 1 5.5", "y",
-        "2 2 2 3.49666295 0.503337045,1 1 1 2.49666295 -0.496662955"),
-    ("2 4 4 4 20,1 1 1 1 5.5", "z",
-        "2 2 2 3.49666295 0.503337045,1 1 1 2.49666295 -0.496662955"),
-    ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "x",
-        "-1 -1 -1 0.496662955 -2.49666295,-2 -2 -2 -0.503337045 -3.49666295"),
-    ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "y",
-        "-1 -1 -1 0.496662955 -2.49666295,-2 -2 -2 -0.503337045 -3.49666295"),
-    ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "z",
-        "-1 -1 -1 0.496662955 -2.49666295,-2 -2 -2 -0.503337045 -3.49666295"),
-    ("2 4 8 12 64,1 1 2 3 11", "x",
-        "2 2 2 3.49666295 0.503337045,1 1 1 2.49666295 -0.496662955"),
-    ("2 4 8 12 64,1 1 2 3 11", "y",
-        "4 4 4 5.49666295 2.50333705,2 2 2 3.49666295 0.503337045"),
-    ("2 4 8 12 64,1 1 2 3 11", "z",
-        "6 6 6 7.49666295 4.50333705,3 3 3 4.49666295 1.503337045"),
-    ("1 -1 -2 -3 11,2 -4 -8 -12 64", "x",
-        "-1 -1 -1 0.496662955 -2.49666295,-2 -2 -2 -0.503337045 -3.49666295"),
-    ("1 -1 -2 -3 11,2 -4 -8 -12 64", "y",
-        "-2 -2 -2 -0.503337045 -3.49666295,-4 -4 -4 -2.50333705 -5.49666295"),
-    ("1 -1 -2 -3 11,2 -4 -8 -12 64", "z",
-        "-3 -3 -3 -1.50333705 -4.49666295,-6 -6 -6 -4.50333705 -7.49666295")
-    ])
-def test_pointwise_eigenvalues_ideal_gas(
-        ctx_factory, states_str, direction, lam_expected_str):
+def test_oscillation_pos_uniform_grid(ctx_factory, flux_test_data_fixture):
+    data = flux_test_data_fixture
+
+    prg = u.get_weno_program_with_root_kernel("oscillation_pos")
+    queue = u.get_queue(ctx_factory)
+
+    oscillation_dev = u.empty_array_on_device(queue, data.nvars, 3)
+
+    prg(queue, nvars=data.nvars,
+            characteristic_fluxes=data.char_fluxes_pos,
+            oscillation=oscillation_dev)
+
+    u.compare_arrays(oscillation_dev.get(), data.oscillation_pos)
+
+
+def test_oscillation_neg_uniform_grid(ctx_factory, flux_test_data_fixture):
+    data = flux_test_data_fixture
+
+    prg = u.get_weno_program_with_root_kernel("oscillation_neg")
+    queue = u.get_queue(ctx_factory)
+
+    oscillation_dev = u.empty_array_on_device(queue, data.nvars, 3)
+
+    prg(queue, nvars=data.nvars,
+            characteristic_fluxes=data.char_fluxes_neg,
+            oscillation=oscillation_dev)
+
+    u.compare_arrays(oscillation_dev.get(), data.oscillation_neg)
+
+
+def test_flux_splitting_uniform_grid(ctx_factory, flux_test_data_fixture):
+    data = flux_test_data_fixture
+
+    prg = u.get_weno_program_with_root_kernel("split_characteristic_fluxes")
+    queue = u.get_queue(ctx_factory)
+
+    fluxes_pos_dev = u.empty_array_on_device(queue, data.nvars, 6)
+    fluxes_neg_dev = u.empty_array_on_device(queue, data.nvars, 6)
+
+    prg(queue, nvars=data.nvars,
+            generalized_states_frozen=data.states,
+            generalized_fluxes_frozen=data.fluxes,
+            R_inv=data.R_inv,
+            wavespeeds=data.wavespeeds,
+            characteristic_fluxes_pos=fluxes_pos_dev,
+            characteristic_fluxes_neg=fluxes_neg_dev)
+
+    u.compare_arrays(fluxes_pos_dev.get(), data.char_fluxes_pos)
+    u.compare_arrays(fluxes_neg_dev.get(), data.char_fluxes_neg)
+
+
+def test_pointwise_eigenvalues_ideal_gas(ctx_factory, flux_test_data_fixture):
+    data = flux_test_data_fixture
+
     prg = u.get_weno_program_with_root_kernel("pointwise_eigenvalues")
     queue = u.get_queue(ctx_factory)
 
-    nvars = 5
-    dirs = {"x": 1, "y": 2, "z": 3}
+    lam_dev = u.empty_array_on_device(queue, data.nvars, 6)
 
-    states = u.expand_to_6(u.transposed_array_from_string(states_str))
-    lam_dev = u.empty_array_on_device(queue, nvars, 6)
+    prg(queue, nvars=data.nvars, d=data.direction,
+            states=data.states, lambda_pointwise=lam_dev)
 
-    prg(queue, nvars=nvars, d=dirs[direction],
-            states=states, lambda_pointwise=lam_dev)
+    u.compare_arrays(lam_dev.get(), data.lam_pointwise)
 
-    lam_expected = u.expand_to_6(u.transposed_array_from_string(lam_expected_str))
-    u.compare_arrays(lam_dev.get(), lam_expected)
 
+def test_roe_uniform_grid_ideal_gas(ctx_factory, flux_test_data_fixture):
+    data = flux_test_data_fixture
 
-@pytest.mark.slow
-@pytest.mark.parametrize("states_str,fluxes_str,direction", [
-    ("2 4 4 4 20,1 1 1 1 5.5", "4 11.2 8 8 46.4,1 2.6 1 1 7.1", "x"),
-    ("2 4 4 4 20,1 1 1 1 5.5", "4 8 11.2 8 46.4,1 1 2.6 1 7.1", "y"),
-    ("2 4 4 4 20,1 1 1 1 5.5", "4 8 8 11.2 46.4,1 1 1 2.6 7.1", "z"),
-    ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "-1 2.6 1 1 -7.1,-4 11.2 8 8 -46.4", "x"),
-    ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "-1 1 2.6 1 -7.1,-4 8 11.2 8 -46.4", "y"),
-    ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "-1 1 1 2.6 -7.1,-4 8 8 11.2 -46.4", "z"),
-    ("2 4 8 12 64,1 1 2 3 11", "4 11.2 16 24 134.4,1 2.6 2 3 12.6", "x"),
-    ("2 4 8 12 64,1 1 2 3 11", "8 16 35.2 48 268.8,2 2 5.6 6 25.2", "y"),
-    ("2 4 8 12 64,1 1 2 3 11", "12 24 48 75.2 403.2,3 3 6 10.6 37.8", "z"),
-    ("1 -1 -2 -3 11,2 -4 -8 -12 64", "-1 2.6 2 3 -12.6,-4 11.2 16 24 -134.4", "x"),
-    ("1 -1 -2 -3 11,2 -4 -8 -12 64", "-2 2 5.6 6 -25.2,-8 16 35.2 48 -268.8", "y"),
-    ("1 -1 -2 -3 11,2 -4 -8 -12 64", "-3 3 6 10.6 -37.8,-12 24 48 75.2 -403.2", "z")
-    ])
-def test_roe_uniform_grid_ideal_gas(ctx_factory, states_str, fluxes_str, direction):
     def identity_matrix(n):
-        return np.identity(n).astype(np.float32).copy(order="F")
+        return np.identity(n).astype(np.float64).copy(order="F")
 
     def check_roe_identity(states, R, R_inv):
         d_state = states[:, 1] - states[:, 0]
@@ -364,33 +213,49 @@ def test_roe_uniform_grid_ideal_gas(ctx_factory, states_str, fluxes_str, directi
     prg = u.get_weno_program_with_root_kernel("roe_eigensystem")
     queue = u.get_queue(ctx_factory)
 
-    nvars = 5
-    ndim = 3
-    dirs = {"x": 1, "y": 2, "z": 3}
-
-    states = u.transposed_array_from_string(states_str)
-    fluxes = u.transposed_array_from_string(fluxes_str)
-    metrics_frozen = identity_matrix(ndim)
+    metrics_frozen = identity_matrix(data.ndim)
 
-    R_dev = u.empty_array_on_device(queue, nvars, nvars)
-    R_inv_dev = u.empty_array_on_device(queue, nvars, nvars)
-    lam_dev = u.empty_array_on_device(queue, nvars)
+    R_dev = u.empty_array_on_device(queue, data.nvars, data.nvars)
+    R_inv_dev = u.empty_array_on_device(queue, data.nvars, data.nvars)
+    lam_dev = u.empty_array_on_device(queue, data.nvars)
 
-    prg(queue, nvars=nvars, ndim=ndim, d=dirs[direction],
-            states=states, metrics_frozen=metrics_frozen,
+    prg(queue, nvars=data.nvars, ndim=data.ndim, d=data.direction,
+            states=data.state_pair, metrics_frozen=metrics_frozen,
             R=R_dev, R_inv=R_inv_dev, lambda_roe=lam_dev)
 
     R = R_dev.get()
     R_inv = R_inv_dev.get()
     lam = lam_dev.get()
 
-    print(lam)
+    check_roe_identity(data.state_pair, R, R_inv)
+    check_roe_property(data.state_pair, data.flux_pair, R, R_inv, lam)
 
-    check_roe_identity(states, R, R_inv)
-    check_roe_property(states, fluxes, R, R_inv, lam)
+
+@pytest.mark.parametrize("lam_pointwise_str,lam_roe_str,lam_expected_str", [
+    ("1 2 3 4 5,2 4 6 8 10", "1.5 3 4.5 6 7.5", "2.2 4.4 6.6 8.8 11"),
+    ("1 2 3 4 5,-2 -4 -6 -8 -10", "1.5 3 4.5 6 7.5", "2.2 4.4 6.6 8.8 11"),
+    ("1 2 3 4 5,-2 -4 -6 -8 -10", "3 6 9 12 15", "3.3 6.6 9.9 13.2 16.5"),
+    ("1 2 3 4 5,2 4 6 8 10", "-3 -6 -9 -12 -15", "3.3 6.6 9.9 13.2 16.5"),
+    ("3 2 9 4 5,2 6 6 12 10", "-1 -4 -3 -8 -15", "3.3 6.6 9.9 13.2 16.5")
+    ])
+def test_lax_wavespeeds(
+        ctx_factory, lam_pointwise_str, lam_roe_str, lam_expected_str):
+    prg = u.get_weno_program_with_root_kernel("lax_wavespeeds")
+    queue = u.get_queue(ctx_factory)
+
+    nvars = 5
+
+    lam_pointwise = u.expand_to_6(u.transposed_array_from_string(lam_pointwise_str))
+    lam_roe = u.array_from_string(lam_roe_str)
+    lam_dev = u.empty_array_on_device(queue, nvars)
+
+    prg(queue, nvars=nvars, lambda_pointwise=lam_pointwise,
+            lambda_roe=lam_roe, wavespeeds=lam_dev)
+
+    lam_expected = u.array_from_string(lam_expected_str)
+    u.compare_arrays(lam_dev.get(), lam_expected)
 
 
-@pytest.mark.slow
 def test_matvec(ctx_factory):
     prg = u.get_weno_program_with_root_kernel("mult_mat_vec")
     queue = u.get_queue(ctx_factory)
diff --git a/utilities.py b/utilities.py
index cce38044b20a7323e342ce9e55900cff18f68761..c8bde8a3b7baf9e034541188c01b6e73e3bfc00e 100644
--- a/utilities.py
+++ b/utilities.py
@@ -11,7 +11,7 @@ from pytest import approx
 # {{{ arrays
 
 def compare_arrays(a, b):
-    assert a == approx(b, rel=1e-5)
+    assert a == approx(b, rel=1e-5, abs=2e-5)
 
 
 def random_array_on_device(queue, *shape):
@@ -21,7 +21,7 @@ def random_array_on_device(queue, *shape):
 
 
 def empty_array_on_device(queue, *shape):
-    return cl.array.empty(queue, shape, dtype=np.float32, order="F")
+    return cl.array.empty(queue, shape, dtype=np.float64, order="F")
 
 
 def arrays_from_string(string_arrays):
@@ -38,7 +38,7 @@ def array_from_string(string_array):
             return np.array(split_map_to_list(string_array[1:], int, " "))
         else:
             return np.array(
-                    split_map_to_list(string_array, float, " "), dtype=np.float32)
+                    split_map_to_list(string_array, float, " "), dtype=np.float64)
 
     def array_from_string_2d(string_array):
         if string_array[0] == ",":