From c9e98bff1857510249a20efaa0c2f6de4c527b75 Mon Sep 17 00:00:00 2001
From: Florian Atteneder <florian.atteneder@uni-jena.de>
Date: Tue, 16 Jul 2024 17:21:28 +0200
Subject: [PATCH] teach spherical1d how deal with tov and bondi data

---
 src/GRHD/rhs.jl | 38 +++++++++++++++++++++++++++++++-------
 1 file changed, 31 insertions(+), 7 deletions(-)

diff --git a/src/GRHD/rhs.jl b/src/GRHD/rhs.jl
index 61f91d83..96c278bc 100644
--- a/src/GRHD/rhs.jl
+++ b/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
-- 
GitLab