diff --git a/src/EulerEq/equation.jl b/src/EulerEq/equation.jl
index e68070cd35012e69aaef0092779a16d196dfa363..7e8a3f3cbad7ce0cb0e51f2cace6800a9c9c0183 100644
--- a/src/EulerEq/equation.jl
+++ b/src/EulerEq/equation.jl
@@ -43,7 +43,7 @@ end
 
 
 @with_signature function entropy(equation::EulerEquation{EquationOfState.IdealGas})
-  @accepts rho, q, p, x
+  @accepts rho, p, x
   @unpack eos  = equation
   gamma = eos.gamma
   S = rho / (gamma - 1) * log(p / rho^gamma)
@@ -64,7 +64,7 @@ end
   @accepts Prefix(rhs), rho, flx_rho, q, flx_q, E, flx_E, v
   @accepts nx
 
-  vmax = max(lhs_v, rhs_v)
+  vmax = absmax(lhs_v, rhs_v)
 
   lhs_nflx = nx*lhs_flx_rho
   rhs_nflx = nx*rhs_flx_rho
@@ -108,7 +108,7 @@ end
   rhs_flx_q   = rhs_q^2 / rhs_rho + rhs_p
   rhs_flx_E   = rhs_q / rhs_rho * (rhs_E + rhs_p)
 
-  vmax = max(lhs_v, rhs_init_v)
+  vmax = absmax(lhs_v, rhs_init_v)
 
   lhs_nflx = nx*lhs_flx_rho
   rhs_nflx = nx*rhs_flx_rho
@@ -156,17 +156,41 @@ end
 end
 
 
+@with_signature function ldg_nflux(eq::EulerEquation)
+  @accepts Prefix(lhs), rho, q, E
+  @accepts Prefix(rhs), rho, q, E
+  @accepts nx
+
+  nflx_rho = 0.5 * nx * (lhs_rho + rhs_rho)
+  nflx_q   = 0.5 * nx * (lhs_q + rhs_q)
+  nflx_E   = 0.5 * nx * (lhs_E + rhs_E)
+
+  @returns nflx_rho, nflx_q, nflx_E
+end
+@with_signature function ldg_bdryllf(eq::EulerEquation)
+  @accepts Prefix(lhs), rho, q, E
+  @accepts Prefix(rhs), rho, q, E
+  @accepts nx
+
+  # set exterior state to 0 according to stability analysis
+  # see dg-notes/tex/evm.pdf for reasoning
+  nflx_rho = nx * lhs_rho
+  nflx_q   = nx * lhs_q
+  nflx_E   = nx * lhs_E
+
+  @returns nflx_rho, nflx_q, nflx_E
+end
+
+
 #######################################################################
 #                                 2D                                  #
 #######################################################################
 
 
-LLF(nfl,nfr,ul,ur,v) = 1/2 * (nfl + nfr - v * (ur - ul))
-
-
 @with_signature function internal_energy(::EulerEquation2d)
-  @accepts rho, q, E
-  eps = E / rho - q^2 / (2 * rho^2)
+  @accepts rho, qx, qy, E
+  q2 = qx^2 + qy^2
+  eps = E / rho - q2 / (2 * rho^2)
   @returns eps
 end
 
@@ -174,7 +198,8 @@ end
 @with_signature function pressure(equation::EulerEquation2d)
   @accepts rho, qx, qy, E
   q = sqrt(qx^2+qy^2)
-  eps, = internal_energy((rho,q,E), equation)
+  eps, = internal_energy((rho,qx,qy,E), equation)
+  @unpack eos = equation
   p = eos(Pressure, Density, rho, InternalEnergy, eps)
   @returns p
 end
@@ -184,7 +209,7 @@ end
   @accepts rho, qx, qy, E, p
   @unpack eos = equation
   q = sqrt(qx^2+qy^2)
-  eps, = internal_energy((rho,q,E), equation)
+  eps, = internal_energy((rho,qx,qy,E), equation)
   p = eos(Pressure, Density, rho, InternalEnergy, eps)
   flx_rho_x = qx
   flx_rho_y = qy
@@ -197,13 +222,14 @@ end
 end
 
 
+# TODO Replace with llf(::EulerEquation2d) below
 @with_signature function llf_rho(equation::EulerEquation2d)
   @accepts Prefix(lhs), rho, flx_rho_x, flx_rho_y, v
   @accepts Prefix(rhs), rho, flx_rho_x, flx_rho_y, v
   @accepts nx, ny
   lhs_nflx = nx*lhs_flx_rho_x + ny*lhs_flx_rho_y
   rhs_nflx = nx*rhs_flx_rho_x + ny*rhs_flx_rho_y
-  vmax = max(lhs_v, rhs_v)
+  vmax = absmax(lhs_v, rhs_v)
   nflx_rho = LLF(lhs_nflx, rhs_nflx, lhs_rho, rhs_rho, vmax)
   @returns nflx_rho
 end
@@ -213,7 +239,7 @@ end
   @accepts nx, ny
   lhs_nflx = nx*lhs_flx_qx_x + ny*lhs_flx_qx_y
   rhs_nflx = nx*rhs_flx_qx_x + ny*rhs_flx_qx_y
-  vmax = max(lhs_v, rhs_v)
+  vmax = absmax(lhs_v, rhs_v)
   nflx_qx = LLF(lhs_nflx, rhs_nflx, lhs_qx, rhs_qx, vmax)
   @returns nflx_qx
 end
@@ -223,7 +249,7 @@ end
   @accepts nx, ny
   lhs_nflx = nx*lhs_flx_qy_x + ny*lhs_flx_qy_y
   rhs_nflx = nx*rhs_flx_qy_x + ny*rhs_flx_qy_y
-  vmax = max(lhs_v, rhs_v)
+  vmax = absmax(lhs_v, rhs_v)
   nflx_qy = LLF(lhs_nflx, rhs_nflx, lhs_qy, rhs_qy, vmax)
   @returns nflx_qy
 end
@@ -233,18 +259,18 @@ end
   @accepts nx, ny
   lhs_nflx = nx*lhs_flx_E_x + ny*lhs_flx_E_y
   rhs_nflx = nx*rhs_flx_E_x + ny*rhs_flx_E_y
-  vmax = max(lhs_v, rhs_v)
+  vmax = absmax(lhs_v, rhs_v)
   nflx_E = LLF(lhs_nflx, rhs_nflx, lhs_E, rhs_E, vmax)
   @returns nflx_E
 end
 
 
 @with_signature function llf(equation::EulerEquation2d)
-  @accepts Prefix(lhs), rho, flx_rho_x, flx_rho_y, qx, flx_qx_x, flx_qx_y, qy, flx_qy_x, flx_qy_y, E, flx_E_x, flx_E_y, v
-  @accepts Prefix(rhs), rho, flx_rho_x, flx_rho_y, qx, flx_qx_x, flx_qx_y, qy, flx_qy_x, flx_qy_y, E, flx_E_x, flx_E_y, v
+  @accepts Prefix(lhs), rho, flx_rho_x, flx_rho_y, qx, flx_qx_x, flx_qx_y, qy, flx_qy_y, E, flx_E_x, flx_E_y, v
+  @accepts Prefix(rhs), rho, flx_rho_x, flx_rho_y, qx, flx_qx_x, flx_qx_y, qy, flx_qy_y, E, flx_E_x, flx_E_y, v
   @accepts nx, ny
 
-  vmax = max(lhs_v, rhs_v)
+  vmax = absmax(lhs_v, rhs_v)
 
   lhs_nflx = nx*lhs_flx_rho_x + ny*lhs_flx_rho_y
   rhs_nflx = nx*rhs_flx_rho_x + ny*rhs_flx_rho_y
@@ -254,6 +280,8 @@ end
   rhs_nflx = nx*rhs_flx_qx_x + ny*rhs_flx_qx_y
   nflx_qx = LLF(lhs_nflx, rhs_nflx, lhs_qx, rhs_qx, vmax)
 
+  lhs_flx_qy_x = lhs_flx_qx_y
+  rhs_flx_qy_x = rhs_flx_qx_y
   lhs_nflx = nx*lhs_flx_qy_x + ny*lhs_flx_qy_y
   rhs_nflx = nx*rhs_flx_qy_x + ny*rhs_flx_qy_y
   nflx_qy = LLF(lhs_nflx, rhs_nflx, lhs_qy, rhs_qy, vmax)
@@ -266,11 +294,73 @@ end
 end
 
 
+@with_signature function bdryllf(equation::EulerEquation2d)
+  @accepts Prefix(lhs), rho, flx_rho_x, flx_rho_y,
+                        qx,  flx_qx_x, flx_qx_y,
+                        qy,  #=flx_qy_x,=# flx_qy_y,
+                        E,   flx_E_x, flx_E_y, v
+  @accepts init_rho, init_qx, init_qy, init_E, init_v
+  @accepts nx, ny
+
+  lhs_flx_qy_x = lhs_flx_qx_y
+
+  # enforce dirichlet/neumann conditions following stability analysis
+  # see dg-notes/tex/evm.pdf
+  # Hesthaven, Numerical Methods for Conservation Laws uses same logic
+
+  # Dirichlet conditions
+  rhs_rho = -lhs_rho + 2*init_rho
+  rhs_qx  = -lhs_qx  + 2*init_qx
+  rhs_qy  = -lhs_qy  + 2*init_qy
+  rhs_E   = -lhs_E   + 2*init_E
+
+  # Neumann conditions
+  # rhs_rho = lhs_rho
+  # rhs_qx  = lhs_qx
+  # rhs_qy  = lhs_qy
+  # rhs_E   = lhs_E
+
+  error()
+  @unpack eos = equation
+  rhs_q = sqrt(rhs_qx^2 + rhs_qy^2)
+  rhs_eps, = internal_energy((rhs_rho,rhs_qx,rhs_qy,rhs_E), equation)
+  rhs_p = eos(Pressure, Density, rhs_rho, InternalEnergy, rhs_eps)
+  rhs_flx_rho_x = rhs_qx
+  rhs_flx_rho_y = rhs_qy
+  rhs_flx_qx_x  = rhs_qx^2 / rhs_rho + rhs_p
+  rhs_flx_qx_y  = rhs_qx*rhs_qy / rhs_rho
+  rhs_flx_qy_x  = rhs_flx_qx_y
+  rhs_flx_qy_y  = rhs_qy^2 / rhs_rho + rhs_p
+  rhs_flx_E_x   = rhs_qx / rhs_rho * (rhs_E + rhs_p)
+  rhs_flx_E_y   = rhs_qy / rhs_rho * (rhs_E + rhs_p)
+
+  vmax = absmax(lhs_v, init_v)
+
+  lhs_nflx = nx*lhs_flx_rho_x + ny*lhs_flx_rho_y
+  rhs_nflx = nx*rhs_flx_rho_x + ny*rhs_flx_rho_y
+  nflx_rho = LLF(lhs_nflx, rhs_nflx, lhs_rho, rhs_rho, vmax)
+
+  lhs_nflx = nx*lhs_flx_qx_x + ny*lhs_flx_qx_y
+  rhs_nflx = nx*rhs_flx_qx_x + ny*rhs_flx_qx_y
+  nflx_qx  = LLF(lhs_nflx, rhs_nflx, lhs_qx, rhs_qx, vmax)
+
+  lhs_nflx = nx*lhs_flx_qy_x + ny*lhs_flx_qy_y
+  rhs_nflx = nx*rhs_flx_qy_x + ny*rhs_flx_qy_y
+  nflx_qy  = LLF(lhs_nflx, rhs_nflx, lhs_qy, rhs_qy, vmax)
+
+  lhs_nflx = nx*lhs_flx_E_x + ny*lhs_flx_E_y
+  rhs_nflx = nx*rhs_flx_E_x + ny*rhs_flx_E_y
+  nflx_E = LLF(lhs_nflx, rhs_nflx, lhs_E, rhs_E, vmax)
+
+  @returns nflx_rho, nflx_qx, nflx_qy, nflx_E
+end
+
+
 @with_signature function maxspeed(equation::EulerEquation2d)
   @accepts rho, qx, qy, E
   @unpack eos = equation
   q = sqrt(qx^2+qy^2)
-  eps, = internal_energy((rho,q,E), equation)
+  eps, = internal_energy((rho,qx,qy,E), equation)
   p = eos(Pressure, Density, rho, InternalEnergy, eps)
   v = abs(q / p) + sqrt(abs(eos.gamma * p / rho))
   @returns v
@@ -294,82 +384,50 @@ end
   @returns flx_S_x, flx_S_y
 end
 
-#######################################################################
-#                         LDG implementation                          #
-#######################################################################
-
-
-@with_signature function bak_ldg_flux(equation::EulerEquation)
-  @accepts State(ldg_rho, ldg_q, ldg_E), rho, q, E
-  flx_ldg_rho, flx_ldg_q, flx_ldg_E = rho, q, E
-  @returns flx_ldg_rho, flx_ldg_q, flx_ldg_E
-end
-@with_signature function ldg_nflux(eq::EulerEquation)
-  @accepts Prefix(lhs), rho, q, E
-  @accepts Prefix(rhs), rho, q, E
-  @accepts nx
-
-  nflx_rho = 0.5 * nx * (lhs_rho + rhs_rho)
-  nflx_q   = 0.5 * nx * (lhs_q + rhs_q)
-  nflx_E   = 0.5 * nx * (lhs_E + rhs_E)
-
-  @returns nflx_rho, nflx_q, nflx_E
-end
-@with_signature function ldg_bdryllf(eq::EulerEquation)
-  @accepts Prefix(lhs), rho, q, E
-  @accepts Prefix(rhs), rho, q, E
-  @accepts nx
-
-  # set exterior state to 0 according to stability analysis
-  # see dg-notes/tex/evm.pdf for reasoning
-  nflx_rho = nx * lhs_rho
-  nflx_q   = nx * lhs_q
-  nflx_E   = nx * lhs_E
 
-  @returns nflx_rho, nflx_q, nflx_E
-end
+@with_signature function ldg_nflux(eq::EulerEquation2d)
+  @accepts Prefix(lhs), rho, qx, qy, E
+  @accepts Prefix(rhs), rho, qx, qy, E
+  @accepts nx, ny
 
+  avg_rho = 0.5 * (lhs_rho + rhs_rho)
+  avg_qx  = 0.5 * (lhs_qx + rhs_qx)
+  avg_qy  = 0.5 * (lhs_qy + rhs_qy)
+  avg_E   = 0.5 * (lhs_E + rhs_E)
 
-@with_signature function ldg_speed(equation::EulerEquation)
-  @accepts State(ldg_rho, ldg_q, ldg_E), rho, q, E
-  v_ldg_rho, v_ldg_q, v_ldg_E = 1, 1, 1
-  @returns v_ldg_rho, v_ldg_q, v_ldg_E
-end
+  nflx_rho_x = nx * avg_rho
+  nflx_rho_y = ny * avg_rho
 
+  nflx_qx_x  = nx * avg_qx
+  nflx_qx_y  = ny * avg_qx
 
-# TODO Remove this together with ldg_flux, ldg_speed, av_speed when
-# moving everything to the new broadcast interface
-@with_signature function bak_av_flux(equation::EulerEquation)
-  @accepts State(ldg_rho, ldg_q, ldg_E), smoothed_mu
-  flx_ldg_rho, flx_ldg_q, flx_ldg_E =
-    smoothed_mu * ldg_rho, smoothed_mu * ldg_q, smoothed_mu * ldg_E
-  @returns flx_ldg_rho, flx_ldg_q, flx_ldg_E
-end
+  nflx_qy_x  = nx * avg_qy
+  nflx_qy_y  = ny * avg_qy
 
+  nflx_E_x   = nx * avg_E
+  nflx_E_y   = ny * avg_E
 
-@with_signature function av_speed(equation::EulerEquation)
-  @accepts State(ldg_rho, ldg_q, ldg_E), smoothed_mu
-  v_av_rho, v_av_q, v_av_E = smoothed_mu, smoothed_mu, smoothed_mu
-  @returns v_av_rho, v_av_q, v_av_E
+  @returns nflx_rho_x, nflx_rho_y, nflx_qx_x, nflx_qx_y,
+           nflx_qy_x, nflx_qy_y, nflx_E_x, nflx_E_y
 end
-
-
-@with_signature function ldg_nflux(eq::EulerEquation2d)
-  @accepts Prefix(lhs), rho, qx, qy, E
-  @accepts Prefix(rhs), rho, qx, qy, E
+@with_signature function ldg_bdryllf(eq::EulerEquation2d)
+  @accepts Prefix(in), rho, qx, qy, E
+  @accepts Prefix(out), rho, qx, qy, E
   @accepts nx, ny
 
-  nflx_rho_x = 0.5 * nx * (lhs_rho + rhs_rho)
-  nflx_rho_y = 0.5 * ny * (lhs_rho + rhs_rho)
+  # set exterior state to 0 according to stability analysis
+  # see dg-notes/tex/evm.pdf for reasoning
+  nflx_rho_x = nx * in_rho
+  nflx_rho_y = ny * in_rho
 
-  nflx_qx_x = 0.5 * nx * (lhs_qx + rhs_qx)
-  nflx_qx_y = 0.5 * ny * (lhs_qx + rhs_qx)
+  nflx_qx_x  = nx * in_qx
+  nflx_qx_y  = ny * in_qx
 
-  nflx_qy_x = 0.5 * nx * (lhs_qy + rhs_qy)
-  nflx_qy_y = 0.5 * ny * (lhs_qy + rhs_qy)
+  nflx_qy_x  = nx * in_qy
+  nflx_qy_y  = ny * in_qy
 
-  nflx_E_x = 0.5 * nx * (lhs_E + rhs_E)
-  nflx_E_y = 0.5 * ny * (lhs_E + rhs_E)
+  nflx_E_x   = nx * in_E
+  nflx_E_y   = ny * in_E
 
   @returns nflx_rho_x, nflx_rho_y, nflx_qx_x, nflx_qx_y,
            nflx_qy_x, nflx_qy_y, nflx_E_x, nflx_E_y
@@ -379,25 +437,29 @@ end
 @with_signature function av_flux_2d(eq::EulerEquation2d)
   @accepts flx_rho_x, flx_rho_y, ldg_rho_x, ldg_rho_y,
            flx_qx_x, flx_qx_y, ldg_qx_x, ldg_qx_y,
-           flx_qy_x, flx_qy_y, ldg_qy_x, ldg_qy_y,
+           #=flx_qy_x,=# flx_qy_y, ldg_qy_x, ldg_qy_y,
            flx_E_x, flx_E_y, ldg_E_x, ldg_E_y,
            smoothed_mu
 
+  flx_qy_x = flx_qx_y
+
   flx_rho_x += smoothed_mu * ldg_rho_x
   flx_rho_y += smoothed_mu * ldg_rho_y
 
-  flx_qx_x += smoothed_mu * ldg_qx_x
-  flx_qx_y += smoothed_mu * ldg_qx_y
+  flx_qx_x  += smoothed_mu * ldg_qx_x
+  flx_qx_y  += smoothed_mu * ldg_qx_y
 
-  flx_qy_x += smoothed_mu * ldg_qy_x
-  flx_qy_y += smoothed_mu * ldg_qy_y
+  flx_qy_x  += smoothed_mu * ldg_qy_x
+  flx_qy_y  += smoothed_mu * ldg_qy_y
 
-  flx_E_x += smoothed_mu * ldg_E_x
-  flx_E_y += smoothed_mu * ldg_E_y
+  flx_E_x   += smoothed_mu * ldg_E_x
+  flx_E_y   += smoothed_mu * ldg_E_y
 
   @returns flx_rho_x, flx_rho_y, flx_qx_x, flx_qx_y,
            flx_qy_x, flx_qy_y, flx_E_x, flx_E_y
 end
+
+
 @with_signature function av_nflux_2d(eq::EulerEquation2d)
   @accepts Prefix(lhs), ldg_rho_x, ldg_rho_y, ldg_qx_x, ldg_qx_y, ldg_qy_x, ldg_qy_y, ldg_E_x, ldg_E_y, smoothed_mu
   @accepts Prefix(rhs), ldg_rho_x, ldg_rho_y, ldg_qx_x, ldg_qx_y, ldg_qy_x, ldg_qy_y, ldg_E_x, ldg_E_y, smoothed_mu
diff --git a/src/EulerEq/initialdata.jl b/src/EulerEq/initialdata.jl
index 1b34bccced150bd75d6a4748d87a4b99a3638448..7ba9dd2aa44c3b7715e6ad21ee2feb081ed05b33 100644
--- a/src/EulerEq/initialdata.jl
+++ b/src/EulerEq/initialdata.jl
@@ -136,10 +136,8 @@ function initialdata_equation_of_state!(env, P, eos::EquationOfState.IdealGas, p
     ρ0, vx0, vy0, p0 = [ similar(x) for _ = 1:4 ]
     CI = CartesianIndices((Npts*Kx,Npts*Ky))
     for (i,(xi,yi)) in enumerate(zip(x,y))
-      ix, _ = Tuple(CI[i])
       ylow = 0.25 + Y(xi)
       yup  = 0.75 + Y(xi)
-      # if yi >= ylow[ix] && yi <= yup[ix]
       if yi >= ylow && yi <= yup
         ρ0[i]  = 2.0
         vx0[i] = -0.5
@@ -190,13 +188,20 @@ end
 
 function initialdata_hrsc!(env, P, hrsc::Maybe{HRSC.AbstractHRSC}, prms, mesh) end
 
-function initialdata_hrsc!(env, P, hrsc::HRSC.SmoothedArtificialViscosity, prms, mesh)
+function initialdata_hrsc!(env, P, hrsc::HRSC.SmoothedArtificialViscosity, prms, mesh::Mesh1d)
   @unpack cache       = env
   @unpack mu          = get_cell_variables(mesh.cache)
   @unpack smoothed_mu = get_static_variables(mesh.cache)
   initialdata_hrsc!(env, P, hrsc.av, prms, mesh)
   hrsc.smoother(smoothed_mu, mu, env.mesh, false)
 end
+function initialdata_hrsc!(env, P, hrsc::HRSC.SmoothedArtificialViscosity, prms, mesh::Mesh2d)
+  @unpack cache       = env
+  @unpack mu          = get_cell_variables(mesh.cache)
+  @unpack smoothed_mu = get_static_variables(mesh.cache)
+  initialdata_hrsc!(env, P, hrsc.av, prms, mesh)
+  hrsc.smoother(smoothed_mu, mu, env.mesh, true, true)
+end
 
 function initialdata_hrsc!(env, P::Project, hrsc::HRSC.EntropyViscosity, prms, ::Mesh1d)
   @unpack mu                     = get_cell_variables(env.cache)
@@ -265,7 +270,7 @@ function initialdata_tci!(env, P::Project, tci::TCI.AbstractTCI, prms)
 end
 
 
-function initialdata_bdrycond!(env, P::Project, bdrycond::dg1d.DirichletBC2, mesh)
+function initialdata_bdrycond!(env, P::Project, bdrycond::dg1d.DirichletBC2, mesh::Mesh1d)
   @unpack rho, q, E = get_dynamic_variables(mesh.cache)
   @unpack flx_rho, flx_q, flx_E, v = get_static_variables(mesh.cache)
   @unpack init_flx_rho, init_flx_q, init_flx_E, init_v = get_static_variables(mesh.cache)
@@ -280,6 +285,17 @@ function initialdata_bdrycond!(env, P::Project, bdrycond::dg1d.DirichletBC2, mes
   init_flx_E   .= flx_E
   init_v       .= v
 end
+function initialdata_bdrycond!(env, P::Project, bdrycond::dg1d.DirichletBC2, mesh::Mesh2d)
+  @unpack rho, qx, qy, E                                = get_dynamic_variables(mesh.cache)
+  @unpack init_rho, init_qx, init_qy, init_E, v, init_v = get_static_variables(mesh.cache)
+  broadcast_volume_2!(speed, P.equation, mesh)
+  broadcast_volume_2!(flux, P.equation, mesh)
+  init_rho .= rho
+  init_qx  .= qx
+  init_qy  .= qy
+  init_E   .= E
+  init_v   .= v
+end
 
 
 function initialdata_smooth!(env, P, mesh)
diff --git a/src/EulerEq/rhs.jl b/src/EulerEq/rhs.jl
index a46e85c7f262ba70b6ebde91cd0d895d1bfa4157..d2c5808b311d3ea951a53a67fc25fa412ce51cac 100644
--- a/src/EulerEq/rhs.jl
+++ b/src/EulerEq/rhs.jl
@@ -179,7 +179,7 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, ::Mesh1d,
 
   ## solve auxiliary equation: q + ∂x u = 0
   broadcast_faces_2!(ldg_nflux, equation, mesh)
-  dg1d.broadcast_bdry_2!(ldg_bdryllf, equation, P.bdrycond, mesh)
+  broadcast_bdry_2!(ldg_bdryllf, equation, P.bdrycond, mesh)
   compute_rhs_weak_form!(ldg_rho, rho, nflx_rho, mesh)
   compute_rhs_weak_form!(ldg_q,   q,   nflx_q,   mesh)
   compute_rhs_weak_form!(ldg_E,   E,   nflx_E,   mesh)
@@ -189,7 +189,7 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, ::Mesh1d,
   broadcast_faces_2!(llf, equation, mesh)
   broadcast_volume_2!(av_flux, equation, mesh)
   broadcast_faces_2!(av_nflux, equation, mesh)
-  dg1d.broadcast_bdry_2!(bdryllf, equation, P.bdrycond, mesh)
+  broadcast_bdry_2!(bdryllf, equation, P.bdrycond, mesh)
   # TODO don't we also need av_bdrylff here?
 
   compute_rhs_weak_form!(rhs_rho, flx_rho,  nflx_rho, mesh)
@@ -219,7 +219,8 @@ function rhs!(env, P::Project2d, hrsc::HRSC.AbstractArtificialViscosity, mesh::M
           nflx_qx_x, nflx_qx_y,
           nflx_qy_x, nflx_qy_y,
           nflx_E_x, nflx_E_y             = get_bdry_variables(mesh.cache)
-  @unpack t = get_global_variables(mesh.cache)
+
+  broadcast_volume_2!(pressure, equation, mesh)
 
   ## TODO Maybe implement a macro with this syntax
   # @weakform(mesh) ∂t ldg_rho_x + ∂x rho = [ nflx_rho_x ]
@@ -234,7 +235,8 @@ function rhs!(env, P::Project2d, hrsc::HRSC.AbstractArtificialViscosity, mesh::M
 
   ## solve auxiliary equation: q + ∂x u = 0
   broadcast_faces_2!(ldg_nflux, equation, mesh)
-  # broadcast_boundaryconditions!(central_flux, ldg_bdryconds, cache, mesh, 0.0)
+  # TODO Need to revise Bdry conditions
+  # broadcast_bdry_2!(ldg_bdryllf, equation, P.bdrycond, mesh)
   compute_rhs_weak_form!(ldg_rho_x, rho, 0,   nflx_rho_x, mesh)
   compute_rhs_weak_form!(ldg_rho_y, 0,   rho, nflx_rho_y, mesh)
   compute_rhs_weak_form!(ldg_qx_x,  qx,  0,   nflx_qx_x,  mesh)
@@ -247,11 +249,12 @@ function rhs!(env, P::Project2d, hrsc::HRSC.AbstractArtificialViscosity, mesh::M
   ## compute rhs of regularized equation: ∂t u + ∂x f + ∂x μ q = 0
   broadcast_volume_2!(flux, equation, mesh)
   broadcast_faces_2!(llf, equation, mesh)
+  # TODO Need to revise Bdry conditions
+  # broadcast_bdry_2!(bdryllf, equation, mesh)
 
-  # broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
   broadcast_volume_2!(av_flux_2d, equation, mesh)
   broadcast_faces_2!(av_nflux_2d, equation, mesh)
-  # broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, 0.0)
+  # TODO don't we also need av_bdrylff here?
 
   compute_rhs_weak_form!(rhs_rho, flx_rho_x, flx_rho_y, nflx_rho, mesh)
   compute_rhs_weak_form!(rhs_qx,  flx_qx_x,  flx_qx_y,  nflx_qx,  mesh)
diff --git a/src/EulerEq/setup.jl b/src/EulerEq/setup.jl
index 85da1d33eb53edac63a018b1fbad252266adb226..297ba3c55fa67b42d79d950df90281e944021694 100644
--- a/src/EulerEq/setup.jl
+++ b/src/EulerEq/setup.jl
@@ -25,7 +25,7 @@ function Project(env::Environment, mesh::Mesh1d, prms)
   # TODO Somehow replace _register_variables! with register_variables!
   @unpack cache = env
   _register_variables!(mesh, P)
-  _register_variables!(mesh.cache, bdrycond)
+  _register_variables!(mesh, bdrycond)
   register_variables!(mesh, rsolver, dont_register=true)
   display(cache)
 
@@ -65,6 +65,7 @@ function Project(env::Environment, mesh::Mesh2d, prms)
   equation = EulerEquation2d(eos)
   ldg_rsolver, av_rsolver = nothing, nothing
 
+  bdrycond = dg1d.DirichletBC2()
   P = Project2d(equation, nothing, nothing,
                 nothing, nothing, hrsc, nothing, nothing, nothing)
 
@@ -72,6 +73,7 @@ function Project(env::Environment, mesh::Mesh2d, prms)
   # TODO Somehow replace _register_variables! with register_variables!
   @unpack cache = env
   _register_variables!(mesh, P, mesh)
+  _register_variables!(mesh, bdrycond)
   _register_variables!(mesh, hrsc)
   display(cache)
 
@@ -187,9 +189,14 @@ end
 
 # TODO Already overwritten above
 # function _register_variables!(cache, bdryconds::Nothing) end
-function _register_variables!(cache, bdryconds::dg1d.AbstractBC)
-  register_variables!(cache,
+function _register_variables!(mesh::Mesh1d, bdryconds::dg1d.AbstractBC)
+  register_variables!(mesh.cache,
     static_variablenames  = (:init_rho, :init_q, :init_E, :init_v,
                              :init_flx_rho, :init_flx_q, :init_flx_E)
   )
 end
+function _register_variables!(mesh::Mesh2d, bdryconds::dg1d.AbstractBC)
+  register_variables!(mesh.cache,
+        static_variablenames  = (:init_rho, :init_qx, :init_qy, :init_E, :init_v)
+  )
+end
diff --git a/src/SRHD/equation.jl b/src/SRHD/equation.jl
index 6c7740c7ac08e145ef185a619cdbc85c49db26a6..400d54c38aaef6a440071b20072d279688f52986 100644
--- a/src/SRHD/equation.jl
+++ b/src/SRHD/equation.jl
@@ -312,9 +312,6 @@ end
 end
 
 
-@inline LLF(nfl,nfr,ul,ur,v) = 1/2 * (nfl + nfr - v * (ur - ul))
-
-
 # TODO Remove llf_D, llf_Sx, llf_Sy, llf_tau in favor of lld below
 @with_signature function llf_D(equation::Equation2d)
   @accepts Prefix(lhs), D, flx_D_x, flx_D_y, max_v
diff --git a/src/ScalarEq/equation.jl b/src/ScalarEq/equation.jl
index dd56db2784d4e1fec4abcc93bc51aed45f639ba2..652c4a2f9990ac174fa266790f3e393eb9504b9b 100644
--- a/src/ScalarEq/equation.jl
+++ b/src/ScalarEq/equation.jl
@@ -112,9 +112,6 @@ end
 #######################################################################
 
 
-LLF(nfl,nfr,ul,ur,v) = 1/2 * (nfl + nfr - v * (ur - ul))
-
-
 @with_signature function flux(eq::Advection2d)
   @accepts State(u)
   flx_x, flx_y = eq.v_x * u, eq.v_y * u
diff --git a/src/dg1d.jl b/src/dg1d.jl
index 05e4851b73c4a16d944f45a4077ac72b3d86c15a..ed481f15e15d546c4581ed42dadd11d666841d59 100644
--- a/src/dg1d.jl
+++ b/src/dg1d.jl
@@ -34,9 +34,6 @@ export # parameter handling
        query, splitallext, trim_keys, match_keys, prepend_keys,
        # type utils
        concretetypes,
-       # numeric utilities
-       periodize, periodic_index, # deprecated
-       estimate_convergence_rate,
        # parameter files and file system operations
        mk_outputdir, dirname_path, flatten, get_parameter_file, get_parameter_file_index,
        select_parameter_files, gather_output_directories, gather_parameter_files,
@@ -44,10 +41,18 @@ export # parameter handling
        find_duplicated_parameter_files,
        run_projects, analyze_projects,
        TODO, @TODO, @stringify_names,
-       indent,
-       absmax
+       indent
 include("utils.jl")
 
+
+export LLF,
+       absmax,
+       periodize,
+       estimate_convergence_rate,
+       # deprecated
+       periodic_index
+include("numutils.jl")
+
 export Output, read_variable, read_variable!, read_mesh
 include("output.jl")
 
@@ -86,7 +91,7 @@ export Cache, register_variables!,
        get_bdry_variables, get_rhs_variables, get_static_variables,
        get_variable, wrap_dynamic_variables!, wrap_rhs_variables!,
        broadcast_volume!, variablenames, variablegroups,
-       broadcast_volume_2!, broadcast_faces_2!
+       broadcast_volume_2!, broadcast_faces_2!, broadcast_bdry_2!
 include("cache.jl")
 
 export Mesh, Mesh1d, Mesh2d, MeshInterpolator, layout, n_points, n_cells, dimension,
diff --git a/src/dg_rhs.jl b/src/dg_rhs.jl
index 727330c04be335f7c612eabb8b5e22c5f0026b45..26841c1ee7295a1af98ca9720aace1cd8cfb2977 100644
--- a/src/dg_rhs.jl
+++ b/src/dg_rhs.jl
@@ -73,29 +73,37 @@ function compute_rhs_weak_form!(rhs, fx, fy::Real, nf, mesh::Mesh2d)
   @unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
   Ml_down, Ml_up = Ml_lhs, Ml_rhs
   Nx, Ny = Npts, Npts
-  @inbounds for k = 1:n_cells(mesh.tree)
+  @inbounds for kk = 1:n_cells(mesh.tree)
 
     # bulk
-    idxs_bulk = cellindices(mesh, k)
+    idxs_bulk = cellindices(mesh, kk)
     v_rhs     = view(rhs,  idxs_bulk)
     v_fx      = view(fx,   idxs_bulk)
     v_dxdX    = view(dxdX, idxs_bulk)
+    v_dydY    = view(dydY, idxs_bulk)
 
     # faces
-    idxs_lhs, idxs_rhs, _, _ = faceindices(mesh, k)
+    idxs_lhs, idxs_rhs, idxs_down, idxs_up = faceindices(mesh, kk)
 
     v_nf_lhs  = view(nf, idxs_lhs)
     v_nf_rhs  = view(nf, idxs_rhs)
+    v_nf_down = view(nf, idxs_down)
+    v_nf_up   = view(nf, idxs_up)
 
     # TODO Benchmark and see whether LoopVectorization can used here.
     @inbounds for j=1:Ny, i=1:Nx
       rhsij = 0
       dxdXij = v_dxdX[i,j]
+      dydYij = v_dydY[i,j]
       for k = 1:Nx
         rhsij += MDM[i,k] * v_fx[k,j]
       end
       v_rhs[i,j] = rhsij * dxdXij
       v_rhs[i,j] -= (v_nf_rhs[j] * Ml_rhs[i] + v_nf_lhs[j]  * Ml_lhs[i])  * dxdXij
+      # TODO Can we really neglect num fluxes in orthogonal direction in general?
+      #      Or we need to adjust for triangular meshes?
+      # v_rhs[i,j] -= (v_nf_rhs[j] * Ml_rhs[i] + v_nf_lhs[j]  * Ml_lhs[i])  * dxdXij +
+      #               (v_nf_up[i]  * Ml_up[j]  + v_nf_down[i] * Ml_down[j]) * dydYij
     end
 
   end
@@ -107,29 +115,37 @@ function compute_rhs_weak_form!(rhs, fx::Real, fy, nf, mesh::Mesh2d)
   @unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
   Ml_down, Ml_up = Ml_lhs, Ml_rhs
   Nx, Ny = Npts, Npts
-  @inbounds for k = 1:n_cells(mesh.tree)
+  @inbounds for kk = 1:n_cells(mesh.tree)
 
     # bulk
-    idxs_bulk = cellindices(mesh, k)
+    idxs_bulk = cellindices(mesh, kk)
     v_rhs     = view(rhs,  idxs_bulk)
     v_fy      = view(fy,   idxs_bulk)
+    v_dxdX    = view(dxdX, idxs_bulk)
     v_dydY    = view(dydY, idxs_bulk)
 
     # faces
-    _, _, idxs_down, idxs_up = faceindices(mesh, k)
+    idxs_lhs, idxs_rhs, idxs_down, idxs_up = faceindices(mesh, kk)
 
+    v_nf_lhs  = view(nf, idxs_lhs)
+    v_nf_rhs  = view(nf, idxs_rhs)
     v_nf_down = view(nf, idxs_down)
     v_nf_up   = view(nf, idxs_up)
 
     # TODO Benchmark and see whether LoopVectorization can used here.
     @inbounds for j=1:Ny, i=1:Nx
       rhsij = 0
+      dxdXij = v_dxdX[i,j]
       dydYij = v_dydY[i,j]
       for k = 1:Ny
         rhsij += MDM[j,k] * v_fy[i,k]
       end
-      v_rhs[i,j] += rhsij * dydYij
+      v_rhs[i,j] = rhsij * dydYij
       v_rhs[i,j] -= (v_nf_up[i]  * Ml_up[j]  + v_nf_down[i] * Ml_down[j]) * dydYij
+      # TODO Can we really neglect num fluxes in orthogonal direction in general?
+      #      Or we need to adjust for triangular meshes?
+      # v_rhs[i,j] -= (v_nf_rhs[j] * Ml_rhs[i] + v_nf_lhs[j]  * Ml_lhs[i])  * dxdXij +
+      #               (v_nf_up[i]  * Ml_up[j]  + v_nf_down[i] * Ml_down[j]) * dydYij
     end
 
   end
diff --git a/src/numutils.jl b/src/numutils.jl
new file mode 100644
index 0000000000000000000000000000000000000000..d803d82d4ddd4dbb306719d5610a5ab76852e8da
--- /dev/null
+++ b/src/numutils.jl
@@ -0,0 +1,106 @@
+"""
+    LLF(nfl,nfr,ul,ur,v)
+
+Local Lax-Friedrich flux:
+
+  nf* = 1/2 ( f(ul) + f(ur) - v * (ur - ul) )
+"""
+@inline LLF(nfl,nfr,ul,ur,v) = 1/2 * (nfl + nfr - v * (ur - ul))
+
+
+"""
+    estimate_convergence_rate(n1, n2, f1, f2)
+
+Experimental rate of convergence based on resolutions `n1, n2` and data `f1, f2`.
+"""
+function estimate_convergence_rate(n1, n2, f1, f2)
+  return log10(f1/f2) / log10(n2/n1)
+end
+
+
+"""
+    purge_zeros!(M, thrs=1e-13)
+
+Scan `M` for elements smaller than `thrs` and set them to 0.
+"""
+purge_zeros!(M, thrs=1e-13) = M[abs.(M).<thrs] .= 0.0
+
+
+"""
+    least_sqr_linear(x, y)
+
+Compute a linear least squares fit for data `x, y`.
+Returns fit function `f, a, b`, where `f(x) = a + b * x` and `a,b` are the fitting coefficients.
+Alogrithm taken from https://mathworld.wolfram.com/LeastSquaresFitting.html
+"""
+function least_sqr_linear(x, y)
+
+  n = length(x)
+  @assert n == length(y) "Mismatch in number of data points x, y!"
+
+  yx = dot(y, x)
+  xx = dot(x, x)
+  mean_y = sum(y) / n
+  mean_x = sum(x) / n
+
+  ss_xx = xx - mean_x^2 * n
+  ss_xy = yx - mean_x * mean_y * n
+
+  b = ss_xy / ss_xx
+  a = mean_y - b * mean_x
+
+  f(x) = a + b * x
+  return f, a, b
+end
+
+
+stepdata(x, x0, fL, fR) = [ (xi < x0) ? fL : fR for xi in x ]
+
+
+"""
+    periodize(i, N)
+
+Wrap `i` around period `N` if `i > N`, otherwise return `i`.
+Behaves similar when called with `Int`s or `Real`s.
+"""
+@inline function periodize(i::T, N::S) where {T<:Int,S<:Int}
+  mod = i % N
+  return mod <= 0 ? N+mod : mod
+end
+@inline function periodize(x::T, L::S) where {T<:Real,S<:Real} 
+  xround = floor(x/L) * L
+  return x <= xround ? x : x - xround
+end
+
+
+@deprecate periodic_index(i, N) periodize(i, N)
+
+
+@inline absolute_maximum(x) = norm(x, Inf)
+
+
+"""
+    absolute_maximum_index(x::AbstractArray)
+
+Return index and value of the absolute maximum of array `x`.
+"""
+function absolute_maximum_index(x::AbstractArray)
+  max_i = 1
+  max_x = abs(x[1])
+  for (i,xi) in enumerate(x)
+    abs_x = abs(xi)
+    if abs_x > max_x
+      max_i = i
+      max_x = abs_x
+    end
+  end
+  return max_i, max_x
+end
+
+
+"""
+    absmax(vs...)
+
+Return absolute maximum of `vs...`.
+"""
+@inline absmax(vs...) = mapreduce(abs, max, vs)
diff --git a/src/utils.jl b/src/utils.jl
index 647929159a9b50357cac8131418f472d12898875..088b1485773f907068c1935ab33eb3fc9648c978 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,21 +1,3 @@
-"""
-    estimate_convergence_rate(n1, n2, f1, f2)
-
-Experimental rate of convergence based on resolutions `n1, n2` and data `f1, f2`.
-"""
-function estimate_convergence_rate(n1, n2, f1, f2)
-  return log10(f1/f2) / log10(n2/n1)
-end
-
-
-"""
-    purge_zeros!(M, thrs=1e-13)
-
-Scan `M` for elements smaller than `thrs` and set them to 0.
-"""
-purge_zeros!(M, thrs=1e-13) = M[abs.(M).<thrs] .= 0.0
-
-
 function get_outputdir()
   thisdir = dirname(@__FILE__)
   outputdir = "$thisdir/../outputs/"
@@ -30,37 +12,6 @@ end
 @deprecate get_output_dir() get_outputdir()
 
 
-"""
-    least_sqr_linear(x, y)
-
-Compute a linear least squares fit for data `x, y`.
-Returns fit function `f, a, b`, where `f(x) = a + b * x` and `a,b` are the fitting coefficients.
-Alogrithm taken from https://mathworld.wolfram.com/LeastSquaresFitting.html
-"""
-function least_sqr_linear(x, y)
-
-  n = length(x)
-  @assert n == length(y) "Mismatch in number of data points x, y!"
-
-  yx = dot(y, x)
-  xx = dot(x, x)
-  mean_y = sum(y) / n
-  mean_x = sum(x) / n
-
-  ss_xx = xx - mean_x^2 * n
-  ss_xy = yx - mean_x * mean_y * n
-
-  b = ss_xy / ss_xx
-  a = mean_y - b * mean_x
-
-  f(x) = a + b * x
-  return f, a, b
-end
-
-
-stepdata(x, x0, fL, fR) = [ (xi < x0) ? fL : fR for xi in x ]
-
-
 # stolen from https://discourse.julialang.org/t/convert-type-to-dictionary/6982/7
 typedict(x::T) where T = Dict(fn=>getfield(x,fn) for fn ∈ fieldnames(T))
 
@@ -78,45 +29,6 @@ end
 
 
 
-"""
-    periodize(i, N)
-
-Wrap `i` around period `N` if `i > N`, otherwise return `i`.
-Behaves similar when called with `Int`s or `Real`s.
-"""
-@inline function periodize(i::T, N::S) where {T<:Int,S<:Int}
-  mod = i % N
-  return mod <= 0 ? N+mod : mod
-end
-@inline function periodize(x::T, L::S) where {T<:Real,S<:Real} 
-  xround = floor(x/L) * L
-  return x <= xround ? x : x - xround
-end
-
-
-@deprecate periodic_index(i, N) periodize(i, N)
-
-
-@inline absolute_maximum(x) = norm(x, Inf)
-
-
-function absolute_maximum_index(x::AbstractArray)
-  max_i = 1
-  max_x = abs(x[1])
-  for (i,xi) in enumerate(x)
-    abs_x = abs(xi)
-    if abs_x > max_x
-      max_i = i
-      max_x = abs_x
-    end
-  end
-  return max_i, max_x
-end
-
-
-@inline absmax(vs...) = mapreduce(abs, max, vs)
-
-
 function query_key(prms, key, default)
   value = get(prms, key, missing)
   if ismissing(value)
diff --git a/test/IntegrationTests/refs/euler2d_kelvin_helmholtz_entropyav.toml b/test/IntegrationTests/refs/euler2d_kelvin_helmholtz_entropyav.toml
index c844f0867104f34a91074f5812ceb4bb71487c1f..860254bc7047b6a313e14f1fe975fa8345bfa62a 100644
--- a/test/IntegrationTests/refs/euler2d_kelvin_helmholtz_entropyav.toml
+++ b/test/IntegrationTests/refs/euler2d_kelvin_helmholtz_entropyav.toml
@@ -1,5 +1,6 @@
 [EulerEq]
 id = "kelvin_helmholtz"
+bc = "periodic"
 
 [EquationOfState]
 eos = "idealgas"
@@ -8,12 +9,13 @@ idealgas_gamma = 1.4
 [Mesh]
 range = [ 0.0,1.0, 0.0,1.0 ]
 cfl = 0.1
+# n = [4,4]
 n = [1,1]
-k = [256,256]
+k = [64,64]
 basis = "lgl"
+periodic = [true, true]
 
 [Output]
-every_iteration  = 1
 every_dt  = 0.02
 variables = [ "rho", "qx", "qy", "E", "p", "smoothed_mu" ]
 enable1d  = false
@@ -26,3 +28,5 @@ tend = 2.0
 [HRSC]
 method = "av"
 av_method = "entropy"
+entropy_cmax = 0.1
+entropy_ce = 0.1