diff --git a/src/GRHD/GRHD.jl b/src/GRHD/GRHD.jl
index 6891b3fce5c7c12606dacf3ead6a6b684906932a..ed28be39676899a17bab6da397a7aa031161116c 100644
--- a/src/GRHD/GRHD.jl
+++ b/src/GRHD/GRHD.jl
@@ -62,11 +62,15 @@ dg1d.@parameters GRHD begin
 
   """
   Available options:
-  - `periodic`
+  - `tov_reflective_origin`: for `id = tov` use reflection symmetry across coordinate origin;
+    requires `range = [ 0, R ]` where `R` shall be larger than the star's radius;
+    at the outer boundary we impose vacuum variables
+  - `tov_symmetric_domain`: for `id = tov` to impose vacuum values on both ends;
+    requires `range = [ -R, R ]` where `R` shall be larger than the star's radius
   - `from_id`
   """
-  bc = "periodic"
-  @check bc in [ "periodic", "from_id" ]
+  bc = "from_id"
+  @check bc in [ "tov_reflective_origin", "tov_symmetric_domain", "from_id" ]
 
   """
   Available options:
diff --git a/src/GRHD/initialdata.jl b/src/GRHD/initialdata.jl
index 87f7fd67dd717059f892101ad14dc55d1f50c6f2..5268e247b93031ee81b5fe5a7b38591c3e8f9a40 100644
--- a/src/GRHD/initialdata.jl
+++ b/src/GRHD/initialdata.jl
@@ -77,7 +77,7 @@ function initialdata_equation!(::Val{:bondi_accretion}, env, P::Project{:valenci
 
   @unpack rD, rSr, rτ = get_dynamic_variables(mesh.cache)
   @unpack D, Sr, τ, ρ, p, ϵ, vr, r, max_v,
-          init_D, init_Sr, init_Ï„, init_vr, init_p,
+          init_D, init_Sr, init_τ, init_vr, init_p, init_ρ, init_ϵ,
           init_max_v = get_static_variables(mesh.cache)
   @unpack α, βru, ∂αr, ∂βrrdu, sqrt_detγ,
           γrr, γθθ, γϕϕ, ∂γrrr, ∂γrθθ, ∂γrϕϕ, ∂γθϕϕ,
@@ -175,6 +175,8 @@ function initialdata_equation!(::Val{:bondi_accretion}, env, P::Project{:valenci
   @. init_vr    = vr
   @. init_p     = p
   @. init_max_v = max_v
+  @. init_ρ     = ρ
+  @. init_ϵ     = ϵ
 
   return
 end
@@ -187,7 +189,7 @@ function initialdata_equation!(::Val{:bondi_accretion}, env, P::Project{:spheric
 
   @unpack D, Sr, Ï„ = get_dynamic_variables(mesh.cache)
   @unpack r, max_v,
-          init_D, init_Sr, init_Ï„, init_vr, init_p,
+          init_D, init_Sr, init_τ, init_vr, init_p, init_ρ, init_ϵ,
           init_max_v = get_static_variables(mesh.cache)
   @unpack α, βru, ∂αr, ∂βrrdu,
           γrr, γθθ, γϕϕ, ∂γrrr, ∂γrθθ, ∂γrϕϕ, ∂γθϕϕ,
@@ -295,6 +297,8 @@ function initialdata_equation!(::Val{:bondi_accretion}, env, P::Project{:spheric
   @. init_vr    = vr
   @. init_p     = p
   @. init_max_v = max_v
+  @. init_ρ     = ρ
+  @. init_ϵ     = ϵ
 
   return
 end
@@ -307,7 +311,7 @@ function initialdata_equation!(::Val{:bondi_accretion_infall}, env, P::Project{:
 
   @unpack D, Sr, Ï„ = get_dynamic_variables(mesh.cache)
   @unpack r, max_v,
-          init_D, init_Sr, init_Ï„, init_vr, init_p,
+          init_D, init_Sr, init_τ, init_vr, init_p, init_ρ, init_ϵ,
           init_max_v = get_static_variables(mesh.cache)
   @unpack α, βru, ∂αr, ∂βrrdu,
           γrr, γθθ, γϕϕ, ∂γrrr, ∂γrθθ, ∂γrϕϕ, ∂γθϕϕ,
@@ -418,6 +422,8 @@ function initialdata_equation!(::Val{:bondi_accretion_infall}, env, P::Project{:
   @. init_vr    = vr
   @. init_p     = p
   @. init_max_v = max_v
+  @. init_ρ     = ρ
+  @. init_ϵ     = ϵ
 
   let ref_ρ=ρ
     @unpack ρ = get_static_variables(mesh)
@@ -486,7 +492,7 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:valencia1d}, prms)
 
   @unpack rD, rSr, rτ   = get_dynamic_variables(mesh.cache)
   @unpack D, Sr, τ, ρ, p, ϵ, vr, r, max_v,
-          init_D, init_Sr, init_Ï„, init_vr, init_p,
+          init_D, init_Sr, init_τ, init_vr, init_p, init_ρ, init_ϵ,
           init_max_v    = get_static_variables(mesh.cache)
   @unpack α, βru, ∂αr, ∂βrrdu, sqrt_detγ,
           γrr, γθθ, γϕϕ, ∂γrrr, ∂γrθθ, ∂γrϕϕ, ∂γθϕϕ,
@@ -652,6 +658,8 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:valencia1d}, prms)
   @. init_vr    = vr
   @. init_p     = p
   @. init_max_v = max_v
+  @. init_ρ     = ρ
+  @. init_ϵ     = ϵ
 
 end
 
@@ -672,7 +680,7 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:spherical1d}, prms)
 
   @unpack D, Sr, Ï„ = get_dynamic_variables(mesh.cache)
   @unpack ρ, p, ϵ, vr, r, max_v,
-          init_D, init_Sr, init_Ï„, init_vr, init_p,
+          init_D, init_Sr, init_τ, init_vr, init_p, init_ρ, init_ϵ,
           init_max_v = get_static_variables(mesh.cache)
   @unpack α, βru, ∂αr, ∂βrrdu,
           γrr, γθθ, γϕϕ, ∂γrrr, ∂γrθθ, ∂γrϕϕ, ∂γθϕϕ,
@@ -830,6 +838,8 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:spherical1d}, prms)
   @. init_vr    = vr
   @. init_p     = p
   @. init_max_v = max_v
+  @. init_ρ     = ρ
+  @. init_ϵ     = ϵ
 
 end
 
@@ -851,7 +861,7 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:rescaled_spherical1
   @unpack rD, rSr, rτ = get_dynamic_variables(mesh.cache)
   @unpack ρ, p, ϵ, vr, r, max_v,
           D, Sr, Ï„,
-          init_D, init_Sr, init_Ï„,
+          init_D, init_Sr, init_τ, init_ρ, init_ϵ,
           init_vr, init_p, init_max_v = get_static_variables(mesh.cache)
   @unpack A, ∂Ar, B, ∂Br, γrr = get_static_variables(mesh.cache)
 
@@ -980,6 +990,8 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:rescaled_spherical1
   @. init_vr    = vr
   @. init_p     = p
   @. init_max_v = max_v
+  @. init_ρ     = ρ
+  @. init_ϵ     = ϵ
 
 end
 
diff --git a/src/GRHD/rhs.jl b/src/GRHD/rhs.jl
index 73e7c818eb861d2042cbbe9c6a32ef5a68a8c766..357a65e8fdc469df4b4ac3d8b4cc2c61a8923916 100644
--- a/src/GRHD/rhs.jl
+++ b/src/GRHD/rhs.jl
@@ -375,31 +375,65 @@ function step_fv_slope_limiter!(mesh, P::Project{:spherical1d}, state)
   dl = L/K
 
   id = P.prmsdb["GRHD"]["id"]
+  bc = P.prmsdb["GRHD"]["bc"]
   if id == "bondi_accretion"
-    @unpack r = get_static_variables(mesh)
-    dr = r[2]-r[1]
-    bdry_γrr_l, bdry_γrr_ll = extrapolate_quadratic(γrr, r, 1:3, r[1]-dr, r[1]-2*dr)
-    bdry_α_l,   bdry_α_ll   = extrapolate_quadratic(α,   r, 1:3, r[1]-dr, r[1]-2*dr)
-    bdry_βru_l, bdry_βru_ll = extrapolate_quadratic(βru, r, 1:3, r[1]-dr, r[1]-2*dr)
-    bdry_ρ_l,   bdry_ρ_ll   = extrapolate_quadratic(ρ,   r, 1:3, r[1]-dr, r[1]-2*dr)
-    bdry_vr_l,  bdry_vr_ll  = extrapolate_quadratic(vr,  r, 1:3, r[1]-dr, r[1]-2*dr)
-    bdry_ϵ_l,   bdry_ϵ_ll   = extrapolate_quadratic(ϵ,   r, 1:3, r[1]-dr, r[1]-2*dr)
+    if bc == "from_id"
+      @unpack init_ρ, init_vr, init_ϵ, r = get_static_variables(mesh)
+      dr = r[2]-r[1]
+      bdry_γrr_l, bdry_γrr_ll = extrapolate_quadratic(γrr,     r, 1:3,   r[1]-dr, r[1]-2*dr)
+      bdry_α_l,   bdry_α_ll   = extrapolate_quadratic(α,       r, 1:3,   r[1]-dr, r[1]-2*dr)
+      bdry_βru_l, bdry_βru_ll = extrapolate_quadratic(βru,     r, 1:3,   r[1]-dr, r[1]-2*dr)
+      bdry_ρ_l,   bdry_ρ_ll   = extrapolate_quadratic(ρ,       r, 1:3,   r[1]-dr, r[1]-2*dr)
+      bdry_vr_l,  bdry_vr_ll  = extrapolate_quadratic(vr,      r, 1:3,   r[1]-dr, r[1]-2*dr)
+      bdry_ϵ_l,   bdry_ϵ_ll   = extrapolate_quadratic(ϵ,       r, 1:3,   r[1]-dr, r[1]-2*dr)
+      bdry_γrr_r, bdry_γrr_rr = extrapolate_quadratic(γrr,     r, K-2:K, r[K]+dr, r[K]+2*dr)
+      bdry_α_r,   bdry_α_rr   = extrapolate_quadratic(α,       r, K-2:K, r[K]+dr, r[K]+2*dr)
+      bdry_βru_r, bdry_βru_rr = extrapolate_quadratic(βru,     r, K-2:K, r[K]+dr, r[K]+2*dr)
+      # only outer boundary is inflow bdry
+      bdry_ρ_r,   bdry_ρ_rr   = extrapolate_quadratic(init_ρ,  r, K-2:K, r[K]+dr, r[K]+2*dr)
+      bdry_vr_r,  bdry_vr_rr  = extrapolate_quadratic(init_vr, r, K-2:K, r[K]+dr, r[K]+2*dr)
+      bdry_ϵ_r,   bdry_ϵ_rr   = extrapolate_quadratic(init_ϵ,  r, K-2:K, r[K]+dr, r[K]+2*dr)
+    else
+      TODO("unknown boundary conditions for initial data $id: $bc")
+    end
   elseif id == "tov"
-    bdry_γrr_l, bdry_γrr_ll = γrr[1],  γrr[2]
-    bdry_α_l,   bdry_α_ll   = α[1],    α[2]
-    bdry_βru_l, bdry_βru_ll = -βru[1], -βru[2]
-    bdry_ρ_l,   bdry_ρ_ll   = ρ[1],    ρ[2]
-    bdry_vr_l,  bdry_vr_ll  = -vr[1],  -vr[2]
-    bdry_ϵ_l,   bdry_ϵ_ll   = ϵ[1],    ϵ[2]
+    if bc == "tov_reflective_origin"
+      bdry_γrr_l, bdry_γrr_ll = γrr[1],   γrr[2]
+      bdry_α_l,   bdry_α_ll   = α[1],     α[2]
+      bdry_βru_l, bdry_βru_ll = -βru[1],  -βru[2]
+      bdry_ρ_l,   bdry_ρ_ll   = ρ[1],     ρ[2]
+      bdry_vr_l,  bdry_vr_ll  = -vr[1],   -vr[2]
+      bdry_ϵ_l,   bdry_ϵ_ll   = ϵ[1],     ϵ[2]
+      bdry_γrr_r, bdry_γrr_rr = γrr[end], γrr[end]
+      bdry_α_r,   bdry_α_rr   = α[end],   α[end]
+      bdry_βru_r, bdry_βru_rr = βru[end], βru[end]
+      bdry_ρ_r,   bdry_ρ_rr   = ρ[end],   ρ[end]
+      bdry_vr_r,  bdry_vr_rr  = vr[end],  vr[end]
+      bdry_ϵ_r,   bdry_ϵ_rr   = ϵ[end],   ϵ[end]
+    elseif bc == "tov_symmetric_domain"
+      bdry_γrr_l, bdry_γrr_ll = γrr[1],   γrr[end]
+      bdry_α_l,   bdry_α_ll   = α[1],     α[end]
+      bdry_βru_l, bdry_βru_ll = βru[1],   βru[end]
+      bdry_ρ_l,   bdry_ρ_ll   = ρ[1],     ρ[end]
+      bdry_vr_l,  bdry_vr_ll  = vr[1],    vr[end]
+      bdry_ϵ_l,   bdry_ϵ_ll   = ϵ[1],     ϵ[end]
+      bdry_γrr_r, bdry_γrr_rr = γrr[end], γrr[end]
+      bdry_α_r,   bdry_α_rr   = α[end],   α[end]
+      bdry_βru_r, bdry_βru_rr = βru[end], βru[end]
+      bdry_ρ_r,   bdry_ρ_rr   = ρ[end],   ρ[end]
+      bdry_vr_r,  bdry_vr_rr  = vr[end],  vr[end]
+      bdry_ϵ_r,   bdry_ϵ_rr   = ϵ[end],   ϵ[end]
+    else
+      TODO("unknown boundary conditions for initial data $id: $bc")
+    end
   else
-    TODO("boundary conditions for initial data $id not implemented")
-  end
-  bγrr      = BdryBufferedData(bdry_γrr_ll, bdry_γrr_l, γrr, γrr[end], γrr[end-1])
-  bα        = BdryBufferedData(bdry_α_ll,   bdry_α_l,   α,   α[end],   α[end-1])
-  bβru      = BdryBufferedData(bdry_βru_ll, bdry_βru_l, βru, βru[end], βru[end-1])
-  bρ        = BdryBufferedData(bdry_ρ_ll,   bdry_ρ_l,   ρ,   ρ[end],   ρ[end-1])
-  bvr       = BdryBufferedData(bdry_vr_ll,  bdry_vr_l,  vr,  vr[end],  vr[end-1])
-  bϵ        = BdryBufferedData(bdry_ϵ_ll,   bdry_ϵ_l,   ϵ,   ϵ[end],   ϵ[end-1])
+  end
+  bγrr      = BdryBufferedData(bdry_γrr_ll, bdry_γrr_l, γrr, bdry_γrr_r, bdry_γrr_rr)
+  bα        = BdryBufferedData(bdry_α_ll,   bdry_α_l,   α,   bdry_α_r,   bdry_α_rr)
+  bβru      = BdryBufferedData(bdry_βru_ll, bdry_βru_l, βru, bdry_βru_r, bdry_βru_rr)
+  bρ        = BdryBufferedData(bdry_ρ_ll,   bdry_ρ_l,   ρ,   bdry_ρ_r,   bdry_ρ_rr)
+  bvr       = BdryBufferedData(bdry_vr_ll,  bdry_vr_l,  vr,  bdry_vr_r,  bdry_vr_rr)
+  bϵ        = BdryBufferedData(bdry_ϵ_ll,   bdry_ϵ_l,   ϵ,   bdry_ϵ_r,   bdry_ϵ_rr)
   bρL       = BdryBufferedData(0, 0, ρL,  0, 0)
   bvrL      = BdryBufferedData(0, 0, vrL, 0, 0)
   bϵL       = BdryBufferedData(0, 0, ϵL,  0, 0)
@@ -415,7 +449,7 @@ function step_fv_slope_limiter!(mesh, P::Project{:spherical1d}, state)
   bdϵ_m     = BdryBufferedData(0, 0, dϵ_m, 0, 0)
   bdϵ_p     = BdryBufferedData(0, 0, dϵ_p, 0, 0)
   @unpack max_v = get_static_variables(mesh)
-  bde_max_v = BdryBufferedData(max_v[1], max_v[1], max_v, max_v[end], max_v[end-1])
+  bde_max_v = BdryBufferedData(max_v[1], max_v[1], max_v, max_v[end], max_v[end])
 
   for vars in ( (bρL,bρR,bρ,bdρ_m,bdρ_p),
                 (bvrL,bvrR,bvr,bdvr_m,bdvr_p),
@@ -458,9 +492,11 @@ function step_fv_slope_limiter!(mesh, P::Project{:spherical1d}, state)
   end
 
   # for the outer bdry conditions for the TOV star we cheat here
-  rhs_D[end] = 0
-  rhs_Sr[end] = 0
-  rhs_Ï„[end] = 0
+  if bc == "tov_reflective_origin"
+    rhs_D[end] = rhs_Sr[end] = rhs_Ï„[end] = 0
+  elseif bc == "tov_symmetric_domain"
+    rhs_D[1] = rhs_Sr[1] = rhs_Ï„[1] = rhs_D[end] = rhs_Sr[end] = rhs_Ï„[end] = 0
+  end
 
   if !P.prmsdb["GRHD"]["atm_evolve"]
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
@@ -620,13 +656,47 @@ function step_fv_slope_limiter!(mesh, P::Project{:doublecartoon}, state)
   L, = dg1d.widths(mesh)
   dl = L/K
 
-  bγxx      = BdryBufferedData(γxx[2],      γxx[1],      γxx,      γxx[end],      γxx[end-1])
-  bα        = BdryBufferedData(α[2],        α[1],        α,        α[end],        α[end-1])
-  bβux      = BdryBufferedData(-βux[2],     -βux[1],     βux,      βux[end],      βux[end-1])
-  bsqrtdetγ = BdryBufferedData(sqrtdetγ[2], sqrtdetγ[1], sqrtdetγ, sqrtdetγ[end], sqrtdetγ[end-1])
-  bρ        = BdryBufferedData(ρ[2],   ρ[1],   ρ,  ρ[end],  ρ[end-1])
-  bvx       = BdryBufferedData(-vx[2], -vx[1], vx, vx[end], vx[end-1])
-  bϵ        = BdryBufferedData(ϵ[2],   ϵ[1],   ϵ,  ϵ[end],  ϵ[end-1])
+  bc = P.prmsdb["GRHD"]["bc"]
+  if bc == "tov_symmetric_domain"
+    bdry_γxx_l, bdry_γxx_ll           = γxx[1],        γxx[2]
+    bdry_α_l, bdry_α_ll               = α[1],          α[2]
+    bdry_βux_l, bdry_βux_ll           = βux[1],        βux[2]
+    bdry_sqrtdetγ_l, bdry_sqrtdetγ_ll = sqrtdetγ[1],   sqrtdetγ[2]
+    bdry_ρ_l, bdry_ρ_ll               = ρ[1],          ρ[2]
+    bdry_vx_l, bdry_vx_ll             = vx[1],         vx[2]
+    bdry_ϵ_l, bdry_ϵ_ll               = ϵ[1],          ϵ[2]
+    bdry_γxx_r,      bdry_γxx_rr      = γxx[end],      γxx[end]
+    bdry_α_r,        bdry_α_rr        = α[end],        α[end]
+    bdry_βux_r,      bdry_βux_rr      = βux[end],      βux[end]
+    bdry_sqrtdetγ_r, bdry_sqrtdetγ_rr = sqrtdetγ[end], sqrtdetγ[end]
+    bdry_ρ_r,        bdry_ρ_rr        = ρ[end],        ρ[end]
+    bdry_vx_r,       bdry_vx_rr       = vx[end],       vx[end]
+    bdry_ϵ_r,        bdry_ϵ_rr        = ϵ[end],        ϵ[end]
+  elseif bc == "tov_reflective_origin"
+    bdry_γxx_l, bdry_γxx_ll           = γxx[1],        γxx[2]
+    bdry_α_l, bdry_α_ll               = α[1],          α[2]
+    bdry_βux_l, bdry_βux_ll           = -βux[1],       -βux[2]
+    bdry_sqrtdetγ_l, bdry_sqrtdetγ_ll = sqrtdetγ[1],   sqrtdetγ[2]
+    bdry_ρ_l, bdry_ρ_ll               = ρ[1],          ρ[2]
+    bdry_vx_l, bdry_vx_ll             = -vx[1],        -vx[2]
+    bdry_ϵ_l, bdry_ϵ_ll               = ϵ[1],          ϵ[2]
+    bdry_γxx_r,      bdry_γxx_rr      = γxx[end],      γxx[end]
+    bdry_α_r,        bdry_α_rr        = α[end],        α[end]
+    bdry_βux_r,      bdry_βux_rr      = βux[end],      βux[end]
+    bdry_sqrtdetγ_r, bdry_sqrtdetγ_rr = sqrtdetγ[end], sqrtdetγ[end]
+    bdry_ρ_r,        bdry_ρ_rr        = ρ[end],        ρ[end]
+    bdry_vx_r,       bdry_vx_rr       = vx[end],       vx[end]
+    bdry_ϵ_r,        bdry_ϵ_rr        = ϵ[end],        ϵ[end]
+  else
+    TODO("implement boundary condition $bc")
+  end
+  bγxx      = BdryBufferedData(bdry_γxx_ll,      bdry_γxx_l,      γxx,      bdry_γxx_r,      bdry_γxx_rr)
+  bα        = BdryBufferedData(bdry_α_ll,        bdry_α_l,        α,        bdry_α_r,        bdry_α_rr)
+  bβux      = BdryBufferedData(bdry_βux_ll,      bdry_βux_l,      βux,      bdry_βux_r,      bdry_βux_rr)
+  bsqrtdetγ = BdryBufferedData(bdry_sqrtdetγ_ll, bdry_sqrtdetγ_l, sqrtdetγ, bdry_sqrtdetγ_r, bdry_sqrtdetγ_rr)
+  bρ        = BdryBufferedData(bdry_ρ_ll,        bdry_ρ_l,        ρ,        bdry_ρ_r,        bdry_ρ_rr)
+  bvx       = BdryBufferedData(bdry_vx_ll,       bdry_vx_l,       vx,       bdry_vx_r,       bdry_vx_rr)
+  bϵ        = BdryBufferedData(bdry_ϵ_ll,        bdry_ϵ_l,        ϵ,        bdry_ϵ_r,        bdry_ϵ_rr)
   bρL       = BdryBufferedData(0, 0, ρL,  0, 0)
   bvxL      = BdryBufferedData(0, 0, vxL, 0, 0)
   bϵL       = BdryBufferedData(0, 0, ϵL,  0, 0)
@@ -687,9 +757,11 @@ function step_fv_slope_limiter!(mesh, P::Project{:doublecartoon}, state)
   end
 
   # for the outer bdry conditions for the TOV star we cheat here
-  rhs_rD[end]  = 0
-  rhs_rSx[end] = 0
-  rhs_rτ[end]  = 0
+  if bc == "tov_reflective_origin"
+    rhs_rD[end] = rhs_rSx[end] = rhs_rτ[end] = 0
+  elseif bc == "tov_symmetric_domain"
+    rhs_rD[1] = rhs_rSx[1] = rhs_rτ[1] = rhs_rD[end] = rhs_rSx[end] = rhs_rτ[end] = 0
+  end
 
   if !P.prmsdb["GRHD"]["atm_evolve"]
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
@@ -754,6 +826,8 @@ end
   @unpack bγxx, bα, bβux, bsqrtdetγ, bρ, bvx, bϵ, bρL, bvxL, bϵL, bpL,
           bρR, bvxR, bϵR, bpR, bdρ_m, bdρ_p, bdvx_m, bdvx_p, bdϵ_m, bdϵ_p, bde_max_v,
           nflx_rD_L, nflx_rSx_L, nflx_rτ_L, nflx_rD_R, nflx_rSx_R, nflx_rτ_R, eos = vars
+  R          = MMatrix{5,5,Float64,25}(0 for _ in 1:25)
+  Δu, Δω, up_u = [ MVector{5,Float64}(0 for _ in 1:5) for _ in 1:4 ]
   for k in 1:K
     # numerical fluxes
     # subscript + ... cell interior value
@@ -829,22 +903,21 @@ end
       # also see Font 2000, PRD 61 044011, eq 64, but only contains right eigenvectors
       #
       # here we assume vx=vz=0 and an isotropic metric
-      R = @SMatrix [ KK/(hâ‚€*Wâ‚€)      0.0        0.0       1.0                                            1.0;
+      R .= @SMatrix [ KK/(hâ‚€*Wâ‚€)      0.0        0.0       1.0                                            1.0;
                      vx₀             0.0        0.0       h₀*W₀*(vx₀-(vux₀-Λp)/(γuuxx₀-vux₀*Λp))         h₀*W₀*(vx₀-(vux₀-Λm)/(γuuxx₀-vux₀*Λm));
                      0.0             h₀*γxx₀    0.0       0.0                                            0.0;
                      0.0             0.0        h₀*γxx₀   0.0                                            0.0;
                      1.0-KK/(h₀*W₀)  0.0        0.0       h₀*W₀*(γuuxx₀-vux₀*vux₀)/(γuuxx₀-vux₀*Λp)-1.0  h₀*W₀*(γuuxx₀-vux₀*vux₀)/(γuuxx₀-vux₀*Λm)-1.0
                    ]
-      uâ‚Š = @SVector [ Dâ‚Š, Sxâ‚Š, 0.0, 0.0, Ï„â‚Š ]
-      uâ‚‹ = @SVector [ Dâ‚‹, Sxâ‚‹, 0.0, 0.0, Ï„â‚‹ ]
-      Δω = R \ (u₋-u₊)
+      Δu .= @SVector [ D₋-D₊, Sx₋-Sx₊, 0.0, 0.0, τ₋-τ₊ ]
+      Δω .= R \ Δu
 
       # roe upwinding term
       # cf. Font 2000, PRD 61 044011, eq 64
       # cf. Hesthaven 2018, Numerical methods for conservation laws,
       #   sec. 3.2 - Riemann solver
       #   sec. 6.2.1 - Roe flux
-      up_u = @. abs(λ0)*Δω[1]*R[:,1] + abs(λ0)*Δω[2]*R[:,2] + abs(λ0)*Δω[3]*R[:,3] + abs(λp)*Δω[4]*R[:,4] + abs(λm)*Δω[5]*R[:,5]
+      @views up_u .= @. abs(λ0)*Δω[1]*R[:,1] + abs(λ0)*Δω[2]*R[:,2] + abs(λ0)*Δω[3]*R[:,3] + abs(λp)*Δω[4]*R[:,4] + abs(λm)*Δω[5]*R[:,5]
       up_D, up_Sx, up_Ï„ = up_u[1], up_u[2], up_u[5]
 
       f₊, f₋       = sdγ₀*α₀*D₊*(vxu₊-βux₀/α₀), sdγ₀*α₀*D₋*(vxu₋-βux₀/α₀)
diff --git a/src/GRHD/setup.jl b/src/GRHD/setup.jl
index 775ab6ad0f0769922c29916f31ead1fa30d2caf4..02fbb31312f2164ae17a14f5e0ef8e29eeb8c1e6 100644
--- a/src/GRHD/setup.jl
+++ b/src/GRHD/setup.jl
@@ -138,6 +138,7 @@ function register_evolution!(mesh, P::Project{:valenciad1d})
                              :ρ, :vr, :p, :ϵ,          # primitives
                              :max_v, :r,               # maximum charactersitic speed, radius
                              :init_D, :init_Sr, :init_Ï„, :init_max_v, :init_vr, :init_p,
+                             :init_ρ, :init_ϵ,
                              :E, :Em1, :flr_E, :flr_Em1, :EP,
                              :ρhW2),
     bdry_variablenames    = (:bdry_D, :bdry_Sr, :bdry_Ï„,
@@ -160,6 +161,7 @@ function register_evolution!(mesh, P::Project{:spherical1d})
                              :ρ, :vr, :p, :ϵ,          # primitives
                              :max_v, :r,               # maximum charactersitic speed, radius
                              :init_D, :init_Sr, :init_Ï„, :init_max_v, :init_vr, :init_p,
+                             :init_ρ, :init_ϵ,
                              :E, :Em1, :flr_E, :flr_Em1, :EP,
                              :ρhW2),
     bdry_variablenames    = (:bdry_D, :bdry_Sr, :bdry_Ï„,
@@ -183,6 +185,7 @@ function register_evolution!(mesh::Mesh1d{FVElement}, P::Project{:spherical1d})
                                :ρ, :vr, :p, :ϵ,          # primitives
                                :max_v, :r,               # maximum charactersitic speed, radius
                                :init_D, :init_Sr, :init_Ï„, :init_max_v, :init_vr, :init_p,
+                               :init_ρ, :init_ϵ,
                                :E, :Em1, :flr_E, :flr_Em1, :EP,
                                :ρhW2),
       bdry_variablenames    = (:bdry_D, :bdry_Sr, :bdry_Ï„,
@@ -214,6 +217,7 @@ function register_evolution!(mesh::Mesh1d{FVElement}, P::Project{:spherical1d})
                                :max_v, :r,               # maximum charactersitic speed, radius
                                :E, :Em1, :flr_E, :flr_Em1, :EP,
                                :init_D, :init_Sr, :init_Ï„, :init_max_v, :init_vr, :init_p,
+                               :init_ρ, :init_ϵ,
                                :ρhW2),
       global_variablenames  = (:t, :tm1) # time steps
     )
@@ -235,6 +239,7 @@ function register_evolution!(mesh, P::Project{:rescaled_spherical1d})
                              :ρ, :vr, :p, :ϵ,             # primitives
                              :max_v, :r,                  # maximum charactersitic speed, radius
                              :init_D, :init_Sr, :init_Ï„, :init_max_v, :init_vr, :init_p,
+                             :init_ρ, :init_ϵ,
                              :E, :Em1, :flr_E, :flr_Em1, :EP,
                              :ρhW2),
     bdry_variablenames    = (:bdry_D, :bdry_Sr, :bdry_Ï„,
@@ -291,6 +296,7 @@ function register_evolution!(mesh::Mesh1d{FVElement}, P::Project{:doublecartoon}
                                :D, :Sx, :Ï„,
                                :max_v, :r,                # maximum charactersitic speed, radius
                                :init_D, :init_Sx, :init_Ï„, :init_max_v, :init_vx, :init_p,
+                               :init_ρ, :init_ϵ,
                                :buf_cartoon, :ρhW2),
       bdry_variablenames    = (:bdry_rD, :bdry_rSx, :bdry_rτ,
                                :bdry_ρ, :bdry_ϵ, :bdry_vx, :bdry_p,
diff --git a/test/IntegrationTests/refs/fv_minmod_marquina_grhd_tov_doublecartoon/fv_minmod_marquina_grhd_tov_doublecartoon.toml b/test/IntegrationTests/refs/fv_minmod_marquina_grhd_tov_doublecartoon/fv_minmod_marquina_grhd_tov_doublecartoon.toml
index aed5aeefad3d1288f3338e1ee3c09a1662c10306..7e5a15e1eecbfaeb82e8257139b346e595ab672e 100644
--- a/test/IntegrationTests/refs/fv_minmod_marquina_grhd_tov_doublecartoon/fv_minmod_marquina_grhd_tov_doublecartoon.toml
+++ b/test/IntegrationTests/refs/fv_minmod_marquina_grhd_tov_doublecartoon/fv_minmod_marquina_grhd_tov_doublecartoon.toml
@@ -4,7 +4,7 @@ polytrope_k = 100.0
 polytrope_gamma = 2.0
 
 [GRHD]
-bc = "from_id"
+bc = "tov_reflective_origin"
 id = "tov"
 formulation = "doublecartoon"
 id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
diff --git a/test/IntegrationTests/refs/fv_minmod_roe_grhd_tov_doublecartoon/fv_minmod_roe_grhd_tov_doublecartoon.toml b/test/IntegrationTests/refs/fv_minmod_roe_grhd_tov_doublecartoon/fv_minmod_roe_grhd_tov_doublecartoon.toml
index 72f68b90d9e98ce42ebb10375df751445abdf951..6580da921260896b21ad10d41a870b9ac03f84f2 100644
--- a/test/IntegrationTests/refs/fv_minmod_roe_grhd_tov_doublecartoon/fv_minmod_roe_grhd_tov_doublecartoon.toml
+++ b/test/IntegrationTests/refs/fv_minmod_roe_grhd_tov_doublecartoon/fv_minmod_roe_grhd_tov_doublecartoon.toml
@@ -4,7 +4,7 @@ polytrope_k = 100.0
 polytrope_gamma = 2.0
 
 [GRHD]
-bc = "from_id"
+bc = "tov_reflective_origin"
 id = "tov"
 formulation = "doublecartoon"
 id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
diff --git a/test/IntegrationTests/refs/fv_superbee_grhd_bondi_accretion/output1d.h5 b/test/IntegrationTests/refs/fv_superbee_grhd_bondi_accretion/output1d.h5
index d46e2cf78a0b786695539e95205fe36f35caeca7..9ef091fca03a0e6223b0b126cbae478c4357e44b 100644
Binary files a/test/IntegrationTests/refs/fv_superbee_grhd_bondi_accretion/output1d.h5 and b/test/IntegrationTests/refs/fv_superbee_grhd_bondi_accretion/output1d.h5 differ