diff --git a/src/GRHD/rhs.jl b/src/GRHD/rhs.jl
index dd7ed4804fab56e678f0c073baff0dcd4c3db264..33984c2704e227ac24c56024599542a6db0b1494 100644
--- a/src/GRHD/rhs.jl
+++ b/src/GRHD/rhs.jl
@@ -32,49 +32,40 @@ end
 function timestep(mesh, P::Project, hrsc::Maybe{HRSC.AbstractHRSC})
   compute_maxspeed(P, P.equation, mesh)
   vmax = get_maxspeed(mesh, P.equation)
-  dl = min_grid_spacing(mesh)
-  dt = dl / (vmax * dtfactor(mesh))
+  L, = dg1d.widths(mesh)
+  K, = mesh.tree.dims
+  dL = L/K
+  dt = dL / (vmax * dtfactor(mesh))
   return dt
 end
 dtfactor(mesh::Mesh1d{FVElement}) = 2
 dtfactor(mesh::Mesh1d{SpectralElement}) = 1
 
 function timestep(mesh::Mesh2d{SpectralElement}, P::Project, hrsc::Maybe{HRSC.AbstractHRSC})
-  compute_maxspeed(P, P.equation, mesh)
-  vmax = get_maxspeed(mesh, P.equation)
   @unpack element = mesh
   @unpack N, z = element
-  dx, dy = dg1d.widths(mesh.boxes[1])
-  for box in view(mesh.boxes, 2:length(mesh.boxes))
-    ddx, ddy = dg1d.widths(box)
-    ddx < dx && (dx = ddx)
-    ddy < dy && (dy = ddy)
-  end
-  dz = z[2]-z[1]
-  dl = min(dx, dy) * dz
-  dt = dl / (N * vmax)
+  compute_maxspeed(P, P.equation, mesh)
+  vmax = get_maxspeed(mesh, P.equation)
+  Lx, Ly = dg1d.widths(mesh.boxes[1])
+  Kx, Ky = mesh.tree.dims
+  dL = min(Lx/Kx, Ly/Ky)
+  dt = dL / (vmax * dtfactor(mesh))
   return dt
 end
 
 
 function timestep(mesh, P::Project, hrsc::HRSC.AbstractArtificialViscosity)
   @unpack equation = P
-
-  compute_maxspeed(P, P.equation, mesh)
-  vmax = get_maxspeed(mesh, P.equation)
-
-  # @unpack mu = get_cell_variables(mesh)
-  # mumax = maximum(abs, mu)
   @unpack smoothed_mu = get_static_variables(mesh)
-  mumax = maximum(abs, smoothed_mu)
-
   @unpack element = mesh
   @unpack N = element
+  compute_maxspeed(P, P.equation, mesh)
+  vmax = get_maxspeed(mesh, P.equation)
+  mumax = maximum(abs, smoothed_mu)
   L, = dg1d.widths(mesh)
   K, = mesh.tree.dims
   dL = L/K
   dl = min_grid_spacing(mesh)
-
   dt1 = dl / (vmax * N)
   dt2 = dL^2 / (mumax * N^3)
   dt = min(dt1, dt2)
@@ -1260,28 +1251,60 @@ function rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
   broadcast_volume!(flux_source_spherical1d, equation, mesh)
   impose_symmetry_sources!(P, mesh)
   broadcast_faces!(llf_spherical1d, equation, mesh)
-  if P.prms.problem === :bondi
-    broadcast_bdry!(bdryllf_spherical1d_bondi, equation, mesh)
-  elseif P.prms.problem === :tov
-    if ismirrored(mesh)
+  id = P.prms.id
+  bc = P.prms.bc
+  if id == "bondi_accretion" || id == "bondi_accretion_infall"
+    if bc == "from_id"
+      # already accounts for inflow boundary
+      broadcast_bdry!(bdryllf_spherical1d_bondi, equation, mesh)
+    else
+      TODO("unknown boundary conditions for initial data $id: $bc")
+    end
+  elseif id == "tov"
+    if bc == "tov_symmetric_domain"
+      # nothing todo, atmosphering should take care of outer bdrys
+    elseif bc == "tov_reflective_origin"
       broadcast_bdry!(bdryllf_spherical1d_tov, equation, mesh)
     else
-      impose_symmetry_bdry_values_spherical1d!(mesh)
-      broadcast_bdry!(bdryllf_spherical1d_tov_reflected, equation, mesh)
+      TODO("unknown boundary conditions for initial data $id: $bc")
     end
   else
-    TODO()
+    TODO("unknown boundary conditions for initial data $id: $bc")
   end
 
-  compute_rhs_weak_form!(rhs_D,  flr_D,  src_D,  nflr_D,  mesh)
-  compute_rhs_weak_form!(rhs_Sr, flr_Sr, src_Sr, nflr_Sr, mesh)
-  compute_rhs_weak_form!(rhs_Ï„,  flr_Ï„,  src_Ï„,  nflr_Ï„,  mesh)
+  if :D ∉ P.prms.freeze_vars
+    compute_rhs_weak_form!(rhs_D,  flr_D,  src_D,  nflr_D,  mesh)
+  end
+  if :Sr ∉ P.prms.freeze_vars
+    compute_rhs_weak_form!(rhs_Sr, flr_Sr, src_Sr, nflr_Sr, mesh)
+  end
+  if :τ ∉ P.prms.freeze_vars
+    compute_rhs_weak_form!(rhs_Ï„,  flr_Ï„,  src_Ï„,  nflr_Ï„,  mesh)
+  end
 
   if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_spherical1d, P.equation, mesh)
   end
 
+  if bc == "tov_symmetric_domain"
+    # Need to symmetrize data around origin to avoid asymmetries to blow up there.
+    # This is needed, because near the surface the solution becomes asymmetric,
+    # and those errors are amplified when propagating inwards.
+    # Not sure if this is caused by the LDG/AV method or by the DG method itself.
+    korigin = find_cell_index(0.0, mesh)
+    for (k,(v_rhs_D, v_rhs_Sr, v_rhs_Ï„)) in enumerate(zip(eachcell(mesh,rhs_D),
+                                                          eachcell(mesh,rhs_Sr),
+                                                          eachcell(mesh,rhs_Ï„)))
+      if korigin == k
+        @views v_rhs_D  .= (v_rhs_D .+ v_rhs_D[end:-1:1])./2
+        @views v_rhs_Sr .= (v_rhs_Sr .- v_rhs_Sr[end:-1:1])./2
+        @views v_rhs_Ï„  .= (v_rhs_Ï„ .+ v_rhs_Ï„[end:-1:1])./2
+        break
+      end
+    end
+  end
+
   return
 end
 
@@ -1448,23 +1471,60 @@ function rhs!(mesh, P::Project{:spherical1d}, hrsc::HRSC.AbstractReconstruction)
   broadcast_volume!(flux_source_spherical1d, equation, mesh)
   impose_symmetry_sources!(P, mesh)
   broadcast_faces!(llf_spherical1d, equation, mesh)
-  if P.prms.problem === :bondi
-    broadcast_bdry!(bdryllf_spherical1d_bondi, equation, mesh)
-  elseif P.prms.problem === :tov
-    broadcast_bdry!(bdryllf_spherical1d_tov, equation, mesh)
+  id = P.prms.id
+  bc = P.prms.bc
+  if id == "bondi_accretion" || id == "bondi_accretion_infall"
+    if bc == "from_id"
+      # already accounts for inflow boundary
+      broadcast_bdry!(bdryllf_spherical1d_bondi, equation, mesh)
+    else
+      TODO("unknown boundary conditions for initial data $id: $bc")
+    end
+  elseif id == "tov"
+    if bc == "tov_symmetric_domain"
+      # nothing todo, atmosphering should take care of outer bdrys
+    elseif bc == "tov_reflective_origin"
+      broadcast_bdry!(bdryllf_spherical1d_tov, equation, mesh)
+    else
+      TODO("unknown boundary conditions for initial data $id: $bc")
+    end
   else
-    TODO()
+    TODO("unknown boundary conditions for initial data $id: $bc")
   end
 
-  compute_rhs_weak_form!(rhs_D,  flr_D,  src_D,  nflr_D,  mesh)
-  compute_rhs_weak_form!(rhs_Sr, flr_Sr, src_Sr, nflr_Sr, mesh)
-  compute_rhs_weak_form!(rhs_Ï„,  flr_Ï„,  src_Ï„,  nflr_Ï„,  mesh)
+  if :D ∉ P.prms.freeze_vars
+    compute_rhs_weak_form!(rhs_D,  flr_D,  src_D,  nflr_D,  mesh)
+  end
+  if :Sr ∉ P.prms.freeze_vars
+    compute_rhs_weak_form!(rhs_Sr, flr_Sr, src_Sr, nflr_Sr, mesh)
+  end
+  if :τ ∉ P.prms.freeze_vars
+    compute_rhs_weak_form!(rhs_Ï„,  flr_Ï„,  src_Ï„,  nflr_Ï„,  mesh)
+  end
 
   if !P.prms.atm_evolve
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_spherical1d, P.equation, mesh)
   end
 
+  if bc == "tov_symmetric_domain"
+    # Need to symmetrize data around origin to avoid asymmetries to blow up there.
+    # This is needed, because near the surface the solution becomes asymmetric,
+    # and those errors are amplified when propagating inwards.
+    # Not sure if this is caused by the LDG/AV method or by the DG method itself.
+    korigin = find_cell_index(0.0, mesh)
+    for (k,(v_rhs_D, v_rhs_Sr, v_rhs_Ï„)) in enumerate(zip(eachcell(mesh,rhs_D),
+                                                          eachcell(mesh,rhs_Sr),
+                                                          eachcell(mesh,rhs_Ï„)))
+      if korigin == k
+        @views v_rhs_D  .= (v_rhs_D .+ v_rhs_D[end:-1:1])./2
+        @views v_rhs_Sr .= (v_rhs_Sr .- v_rhs_Sr[end:-1:1])./2
+        @views v_rhs_Ï„  .= (v_rhs_Ï„ .+ v_rhs_Ï„[end:-1:1])./2
+        break
+      end
+    end
+  end
+
   return
 end
 
diff --git a/test/IntegrationTests/refs/grhd_bondi_accretion/grhd_bondi_accretion.toml b/test/IntegrationTests/refs/grhd_bondi_accretion/grhd_bondi_accretion.toml
index b9e0c86d2bcc7fbf47b2505b70dbf241b0884efd..27ba7b96e6919694fc235c27c8a248220c600d92 100644
--- a/test/IntegrationTests/refs/grhd_bondi_accretion/grhd_bondi_accretion.toml
+++ b/test/IntegrationTests/refs/grhd_bondi_accretion/grhd_bondi_accretion.toml
@@ -25,5 +25,5 @@ variables1d = [ "D_err", "D", "Sr", "τ", "p", "ρ", "ϵ", "vr" ]
 aligned_ts  = "$(collect(range(0.5,5.0,step=0.5)))"
 
 [Evolution]
-cfl = 0.25
+cfl = 0.2
 tend = 5.0
diff --git a/test/IntegrationTests/refs/grhd_bondi_accretion/output0d.h5 b/test/IntegrationTests/refs/grhd_bondi_accretion/output0d.h5
index 2372200c9d55488c1dd7572923b9d2a4700f6ac2..9931e08123ca1eec7aeddf3d3b8406eaa683877c 100644
Binary files a/test/IntegrationTests/refs/grhd_bondi_accretion/output0d.h5 and b/test/IntegrationTests/refs/grhd_bondi_accretion/output0d.h5 differ
diff --git a/test/IntegrationTests/refs/grhd_bondi_accretion/output1d.h5 b/test/IntegrationTests/refs/grhd_bondi_accretion/output1d.h5
index a0d6e090c3e1fc69cd25d3161a58ab8fea1f85c0..7219b39bc289421afdd21452654ec09cca851686 100644
Binary files a/test/IntegrationTests/refs/grhd_bondi_accretion/output1d.h5 and b/test/IntegrationTests/refs/grhd_bondi_accretion/output1d.h5 differ
diff --git a/test/IntegrationTests/refs/grhd_tov_doublecartoon/grhd_tov_doublecartoon.toml b/test/IntegrationTests/refs/grhd_tov_doublecartoon/grhd_tov_doublecartoon.toml
index abc308355908be55fadcdd73a1fab3a158064cb8..9ece736d26eb6b6372de3a489f4a3b632873bb55 100644
--- a/test/IntegrationTests/refs/grhd_tov_doublecartoon/grhd_tov_doublecartoon.toml
+++ b/test/IntegrationTests/refs/grhd_tov_doublecartoon/grhd_tov_doublecartoon.toml
@@ -1,26 +1,26 @@
 [EquationOfState]
-eos = "polytrope"
-polytrope_k = 100.0
 polytrope_gamma = 2.0
+polytrope_k = 100.0
+eos = "polytrope"
+
+[Evolution]
+cfl = 0.1
+tend = 10.0
+
+[Output]
+aligned_ts = "$(collect(range(2.5,10.0,step=2.5)))"
+variables1d = ["D", "Sx", "Ï„"]
 
 [GRHD]
-bc = "from_id"
-id = "tov"
 id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
-atm_factor = 1e-8
+atm_factor = 1.0e-8
+id = "tov"
 formulation = "doublecartoon"
+bc = "from_id"
 
 [Mesh]
-range = [ 0, 20.0 ]
-n = 3
+periodic = false
+range = [0, 20.0]
 k = 27
 basis = "lgl"
-periodic = false
-
-[Output]
-variables1d = [ "D", "Sx", "Ï„" ]
-aligned_ts = "$(collect(range(2.5,10.0,step=2.5)))"
-
-[Evolution]
-cfl = 0.3
-tend = 10.0
+n = 3
diff --git a/test/IntegrationTests/refs/grhd_tov_doublecartoon/output1d.h5 b/test/IntegrationTests/refs/grhd_tov_doublecartoon/output1d.h5
index 6771b41e68afe3ebac838ec7e3be40278b4f3f22..37c960b7ee85fa721f2a78af1e357397ca9f0e50 100644
Binary files a/test/IntegrationTests/refs/grhd_tov_doublecartoon/output1d.h5 and b/test/IntegrationTests/refs/grhd_tov_doublecartoon/output1d.h5 differ
diff --git a/test/IntegrationTests/refs/grhd_tov_rescaled_spherical1d/grhd_tov_rescaled_spherical1d.toml b/test/IntegrationTests/refs/grhd_tov_rescaled_spherical1d/grhd_tov_rescaled_spherical1d.toml
index e46f282512b7a189ec9eaa9a8cd631ec928f514e..4bd18776c100e5c21a9f72e5aa33967477a4665e 100644
--- a/test/IntegrationTests/refs/grhd_tov_rescaled_spherical1d/grhd_tov_rescaled_spherical1d.toml
+++ b/test/IntegrationTests/refs/grhd_tov_rescaled_spherical1d/grhd_tov_rescaled_spherical1d.toml
@@ -1,27 +1,27 @@
 [EquationOfState]
-eos = "polytrope"
-polytrope_k = 100.0
 polytrope_gamma = 2.0
+polytrope_k = 100.0
+eos = "polytrope"
+
+[Evolution]
+cfl = 0.05
+tend = 10.0
+
+[Output]
+aligned_ts = "$(collect(range(2.5,10.0,step=2.5)))"
+variables1d = ["D", "Sr", "Ï„"]
 
 [GRHD]
-bc = "from_id"
-id = "tov"
-id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
-atm_factor = 1e-8
 atm_threshold_factor = 100.0
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1.0e-8
+id = "tov"
 formulation = "rescaled_spherical1d"
+bc = "from_id"
 
 [Mesh]
-range = [ 0.0, 20.0 ]
-n = 7
+periodic = false
+range = [0.0, 20.0]
 k = 47
 basis = "lgl"
-periodic = false
-
-[Output]
-variables1d = [ "D", "Sr", "Ï„" ]
-aligned_ts = "$(collect(range(2.5,10.0,step=2.5)))"
-
-[Evolution]
-cfl = 0.3
-tend = 10.0
+n = 5
diff --git a/test/IntegrationTests/refs/grhd_tov_rescaled_spherical1d/output1d.h5 b/test/IntegrationTests/refs/grhd_tov_rescaled_spherical1d/output1d.h5
index f981f13a38a55aa8170bfe6313e554effc89c87c..06cef9798b88e4f004a0ef0debf924eb618d2ecb 100644
Binary files a/test/IntegrationTests/refs/grhd_tov_rescaled_spherical1d/output1d.h5 and b/test/IntegrationTests/refs/grhd_tov_rescaled_spherical1d/output1d.h5 differ
diff --git a/test/IntegrationTests/refs/grhd_tov_spherical1d/grhd_tov_spherical1d.toml b/test/IntegrationTests/refs/grhd_tov_spherical1d/grhd_tov_spherical1d.toml
index b2dfe766a6d17bfa010c026be759dfe3b1036716..a37582e4f2c92427242b06092bc3e8c46ee15be0 100644
--- a/test/IntegrationTests/refs/grhd_tov_spherical1d/grhd_tov_spherical1d.toml
+++ b/test/IntegrationTests/refs/grhd_tov_spherical1d/grhd_tov_spherical1d.toml
@@ -1,28 +1,28 @@
 [EquationOfState]
-eos = "polytrope"
-polytrope_k = 100.0
 polytrope_gamma = 2.0
+polytrope_k = 100.0
+eos = "polytrope"
+
+[Evolution]
+cfl = 0.1
+tend = 10.0
+
+[Output]
+aligned_ts = "$(collect(range(2.5,10.0,step=2.5)))"
+variables1d = ["D", "Sr", "Ï„"]
 
 [GRHD]
-bc = "from_id"
-id = "tov"
-id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
-atm_factor = 1e-8
 atm_threshold_factor = 100.0
-formulation = "spherical1d"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1.0e-8
+id = "tov"
 c2p_enforce_causal_atm = false
+formulation = "spherical1d"
+bc = "tov_symmetric_domain"
 
 [Mesh]
-range = [ 0.0, 20.0 ]
-n = 4
-k = 124
-basis = "lgl"
 periodic = false
-
-[Output]
-variables1d = [ "D", "Sr", "Ï„" ]
-aligned_ts  = "$(collect(range(2.5,10.0,step=2.5)))"
-
-[Evolution]
-cfl = 0.4
-tend = 10.0
+range = [-20.0, 20.0]
+k = 203
+basis = "lgl"
+n = 5
diff --git a/test/IntegrationTests/refs/grhd_tov_spherical1d/output1d.h5 b/test/IntegrationTests/refs/grhd_tov_spherical1d/output1d.h5
index 3dc52a15f02e0ba2ae99337fa71a500e2db9347e..0b49f432e61657eb59d28216f8e5aab9841df6f2 100644
Binary files a/test/IntegrationTests/refs/grhd_tov_spherical1d/output1d.h5 and b/test/IntegrationTests/refs/grhd_tov_spherical1d/output1d.h5 differ