GRHD: add tov symmetrization for all DG rhs' that use the spherical1d...

GRHD: add tov symmetrization for all DG rhs' that use the spherical1d (!195)

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
 

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

test/IntegrationTests/refs/grhd_tov_doublecartoon/grhd_tov_doublecartoon.toml

 [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

test/IntegrationTests/refs/grhd_tov_rescaled_spherical1d/grhd_tov_rescaled_spherical1d.toml

 [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

test/IntegrationTests/refs/grhd_tov_spherical1d/grhd_tov_spherical1d.toml

 [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