diff --git a/src/EulerEq/equation.jl b/src/EulerEq/equation.jl
index e02e97032cf089da10489c5e445a0053a7c61e6a..c87d74bd8fd902b64041c2f2874956b8c4088983 100644
--- a/src/EulerEq/equation.jl
+++ b/src/EulerEq/equation.jl
@@ -82,6 +82,29 @@ end
 end
 
 
+@with_signature function bdryllf(equation::EulerEquation)
+  @accepts Prefix(lhs), rho, flx_rho, q, flx_q, E, flx_E, v
+  @accepts Prefix(rhs), init_rho, init_flx_rho, init_q, init_flx_q, init_E, init_flx_E, init_v
+  @accepts nx
+
+  vmax = max(lhs_v, rhs_init_v)
+
+  lhs_nflx = nx*lhs_flx_rho
+  rhs_init_nflx = nx*rhs_init_flx_rho
+  nflx_rho = LLF(lhs_nflx, rhs_init_nflx, lhs_rho, rhs_init_rho, vmax)
+
+  lhs_nflx = nx*lhs_flx_q
+  rhs_init_nflx = nx*rhs_init_flx_q
+  nflx_q = LLF(lhs_nflx, rhs_init_nflx, lhs_q, rhs_init_q, vmax)
+
+  lhs_nflx = nx*lhs_flx_E
+  rhs_init_nflx = nx*rhs_init_flx_E
+  nflx_E = LLF(lhs_nflx, rhs_init_nflx, lhs_E, rhs_init_E, vmax)
+
+  @returns nflx_rho, nflx_q, nflx_E
+end
+
+
 @with_signature function av_flux(eq::EulerEquation)
   @accepts flx_rho, ldg_rho, flx_q, ldg_q, flx_E, ldg_E, smoothed_mu
 
diff --git a/src/EulerEq/initialdata.jl b/src/EulerEq/initialdata.jl
index a84fdcf14ca48ce97c77d4eb5e606ee6120d3421..e6be57c962812120641ff377fe0455b7b6af7145 100644
--- a/src/EulerEq/initialdata.jl
+++ b/src/EulerEq/initialdata.jl
@@ -5,6 +5,7 @@ function initialdata!(env, P::Project, prms, mesh::Mesh1d)
   initialdata_entropy_variables!(env, P, prms)
   initialdata_hrsc!(env, P, P.hrsc, prms, mesh)
   initialdata_tci!(env, P, P.tci, prms)
+  initialdata_bdrycond!(env, P, P.bdrycond, mesh)
 end
 
 function initialdata!(env, P::Project2d, prms, mesh::Mesh2d)
@@ -220,3 +221,20 @@ function initialdata_tci!(env, P::Project, tci::TCI.AbstractTCI, prms)
   @unpack flag = get_cell_variables(env.cache)
   @. flag = 0.0
 end
+
+
+function initialdata_bdrycond!(env, P::Project, bdrycond::dg1d.DirichletBC2, mesh)
+  @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)
+  @unpack init_rho, init_q, init_E = get_static_variables(mesh.cache)
+  broadcast_volume_2!(speed, P.equation, mesh)
+  broadcast_volume_2!(flux, P.equation, mesh)
+  init_rho     .= rho
+  init_q       .= q
+  init_E       .= E
+  init_flx_rho .= flx_rho
+  init_flx_q   .= flx_q
+  init_flx_E   .= flx_E
+  init_v       .= v
+end
diff --git a/src/EulerEq/rhs.jl b/src/EulerEq/rhs.jl
index bc3e324af24b8e5bf98294bb1529ad23e18deb03..c9eb026e54898d355e02ff4f0f9de763c9db7593 100644
--- a/src/EulerEq/rhs.jl
+++ b/src/EulerEq/rhs.jl
@@ -203,6 +203,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)
 
   compute_rhs_weak_form!(rhs_rho, flx_rho,  nflx_rho, mesh)
   compute_rhs_weak_form!(rhs_q,   flx_q,    nflx_q,   mesh)
diff --git a/src/EulerEq/setup.jl b/src/EulerEq/setup.jl
index 2f52d1448a4c0e0b27c60e47c13d6b8ac8bcc4a2..594183ec1a93e13314c4cabcb887419f3ba164a3 100644
--- a/src/EulerEq/setup.jl
+++ b/src/EulerEq/setup.jl
@@ -17,18 +17,21 @@ function Project(env::Environment, mesh::Mesh1d, prms)
     nothing, nothing
   end
 
-  P = Project(equation, rsolver, hrsc, tci, ldg_rsolver, av_rsolver)
+  bdrycond = dg1d.DirichletBC2()
+  P = Project(equation, rsolver, hrsc, tci, ldg_rsolver, av_rsolver, bdrycond)
 
   # register variables
   # TODO add a ::Nothing overload for register_variables!
   # TODO Somehow replace _register_variables! with register_variables!
   @unpack cache = env
   _register_variables!(mesh, P)
+  _register_variables!(mesh.cache, bdrycond)
   register_variables!(mesh, rsolver, dont_register=true)
-  display(cache)
 
   # setup initial data
   initialdata!(env, P, prms["EulerEq"])
+  # TODO Need to setup some initial viscosity in order to have a chance to even
+  # start evolving with the EV method and higher orders
 
   # setup boundary conditions
   bdryconds = make_BoundaryConditions(env, equation, rsolver, prms["EulerEq"])
@@ -179,3 +182,13 @@ function _register_variables!(mesh, P::Project2d, ::Mesh2d)
     global_variablenames  = (:t,) # recent time stamps
   )
 end
+
+
+# TODO Already overwritten above
+# function _register_variables!(cache, bdryconds::Nothing) end
+function _register_variables!(cache, bdryconds::dg1d.AbstractBC)
+  register_variables!(cache,
+    static_variablenames  = (:init_rho, :init_q, :init_E, :init_v,
+                             :init_flx_rho, :init_flx_q, :init_flx_E)
+  )
+end
diff --git a/src/EulerEq/types.jl b/src/EulerEq/types.jl
index f34fad55990b67638b402a05451041387bcb5c78..9e8a4f087e788e2f746e03c42cc86081d84f33f9 100644
--- a/src/EulerEq/types.jl
+++ b/src/EulerEq/types.jl
@@ -12,7 +12,8 @@ mutable struct Project{T_Equation <:EulerEquation,
                        T_HRSC     <:Maybe{HRSC.AbstractHRSC},
                        T_TCI      <:Maybe{TCI.AbstractTCI},
                        T_LDG_RSolver<:Maybe{AbstractRiemannSolver},
-                       T_AV_RSolver<:Maybe{AbstractRiemannSolver}}
+                       T_AV_RSolver<:Maybe{AbstractRiemannSolver},
+                       T_BC}
 
   equation::T_Equation
   rsolver::T_RSolver
diff --git a/test/IntegrationTests/refs/euler_sod_shock_tube_entropyav.toml b/test/IntegrationTests/refs/euler_sod_shock_tube_entropyav.toml
index e2e48127a93a9734714d7c020dc91d1958089c8d..d215c02511843b8c15db9a6ad525caaf0e6d4be0 100644
--- a/test/IntegrationTests/refs/euler_sod_shock_tube_entropyav.toml
+++ b/test/IntegrationTests/refs/euler_sod_shock_tube_entropyav.toml
@@ -10,7 +10,9 @@ id = "sod_shock_tube"
 range = [ -1.0, 1.0 ]
 cfl = 0.5
 n = 1
-k = 1000
+k = 500
+# n = 1
+# k = 1000
 basis = "lgl"
 
 [Output]
@@ -22,7 +24,7 @@ variables       = [ "rho", "q", "E", "smoothed_mu", "S", "v" ]
 tend = 0.6
 
 [Log]
-progress_stdout = false
+progress_stdout = true
 
 [HRSC]
 method = "av"