teach spherical1d how deal with tov and bondi data

src/GRHD/rhs.jl

@@ -1278,17 +1278,24 @@ 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.prmsdb["GRHD"]["id"]
+  bc = P.prmsdb["GRHD"]["bc"]
+  if id == "bondi_accretion"
+    if bc == "from_id"
+      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
 
   if "D" ∉ P.prmsdb["GRHD"]["freeze_vars"]
@@ -1304,6 +1311,23 @@ function rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
   if !P.prmsdb["GRHD"]["atm_evolve"]
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_spherical1d, P.equation, mesh)
+
+  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