fix missing lapse factor in source terms

src/GRHD/equation.jl

@@ -261,7 +261,7 @@ end
 
 
 @with_signature function flux_source_rescaled_spherical1d(eq::Equation)
-  @accepts D, Sr, τ, ϵ, ρ, p, vr
+  @accepts D, Sr, τ, p, vr
   @accepts A, ∂Ar, B, ∂Br, r
 
   γuurr = 1/B^2
@@ -273,8 +273,11 @@ end
   flr_rτ  = r*A*B3 * (τ + p) * vru
 
   src_rD  = - A*B3 * D * vru
-  src_rSr = -(τ+p)*r*B3*∂Ar + A*B^2*(r*∂Br+B/3)*(Sr*vru+3*p) - 4/3*A*B3*Sr*vru
-  src_rτ  = -r*B*Sr*∂Ar - A*B3*(τ+p)*vru
+  # src_rSr = -(τ+p)*r*B3*∂Ar + A*B^2*(r*∂Br+B/3)*(Sr*vru+3*p) - 4/3*A*B3*Sr*vru
+  # src_rSr = -(τ+p)*r*B3*∂Ar + A*B^2*(r*∂Br+B/3)*(Sr*vru/A^2+3*p) - 4/3*A*B3*Sr*vru
+  src_rSr = -(τ+p)*r*B3*∂Ar + A*B^2*(r*∂Br+B/3)*(Sr*vru/A^2+3*p) - 4/3*A*B3*Sr*vru
+  # src_rτ  = -r*B*Sr*∂Ar - A*B3*(τ+p)*vru
+  src_rτ  = -r*B*Sr*∂Ar/A - A*B3*(τ+p)*vru
 
   @returns src_rD, src_rSr, src_rτ, flr_rD, flr_rSr, flr_rτ
 end
@@ -305,8 +308,8 @@ end
   vmax     = absmax(max_v, bdry_max_v)
   γuurr    = 1/B^2
   bdry_vru = γuurr * bdry_vr
+  B3       = B^3
 
-  B3 = B^3
   rD       = r*B3*D
   rSr      = r*B3*Sr
   rτ       = r*B3*τ
@@ -343,33 +346,20 @@ end
   γuurr    = 1/B^2
   init_vru = γuurr * init_vr
   vru      = γuurr * vr
-
-  # enforce dirichlet/neumann conditions following stability analysis
-  # see dg-notes/tex/evm.pdf
-  # Hesthaven, Numerical Methods for Conservation Laws uses same logic
+  B3       = B^3
 
   # Dirichlet conditions
-  # ## if block for hand-tweaked outflow condition for bondi accretion at inner radius
-  # if nx > 0
-    bdry_D   = -D   + 2*init_D
-    bdry_Sr  = -Sr  + 2*init_Sr
-    bdry_τ   = -τ   + 2*init_τ
-    bdry_vru = -vru + 2*init_vru
-    bdry_p   = -p   + 2*init_p
-  # else
-  #   bdry_D   = D
-  #   bdry_Sr  = Sr
-  #   bdry_τ   = τ
-  #   bdry_vru = vru
-  #   bdry_p   = p
-  # end
+  bdry_D   = -D   + 2*init_D
+  bdry_Sr  = -Sr  + 2*init_Sr
+  bdry_τ   = -τ   + 2*init_τ
+  bdry_vru = -vru + 2*init_vru
+  bdry_p   = -p   + 2*init_p
 
   # Neumann conditions
   # bdry_D   = D
   # bdry_Sr   = Sr
   # bdry_τ = τ
 
-  B3 = B^3
   rD       = r*B3*D
   rSr      = r*B3*Sr
   rτ       = r*B3*τ
@@ -379,7 +369,7 @@ end
 
   bdry_flr_rD  = r*A*B3 * bdry_D * bdry_vru
   bdry_flr_rSr = r*A*B3 * ( bdry_Sr * bdry_vru + bdry_p )
-  bdry_flr_rτ  = r*A*B3 * ( (bdry_τ + bdry_p) * bdry_vru )
+  bdry_flr_rτ  = r*A*B3 * (bdry_τ + bdry_p) * bdry_vru
 
   nflr      = nx*flr_rD
   bdry_nflr = nx*bdry_flr_rD

src/GRHD/rhs.jl

@@ -2,6 +2,12 @@
 #                              Timestep                               #
 #######################################################################
 
+function compute_maxspeed(::Project{:spherical1d}, equation, mesh)
+  broadcast_volume!(maxspeed, equation, mesh)
+end
+function compute_maxspeed(::Project{:rescaled_spherical1d}, equation, mesh)
+  broadcast_volume!(maxspeed_rescaled_spherical1d, equation, mesh)
+end
 
 function get_maxspeed(mesh, equation)
   @unpack max_v = get_static_variables(mesh.cache)
@@ -14,13 +20,6 @@ function get_maxspeed(mesh, equation)
   return vmax
 end
 
-function compute_maxspeed(::Project{:spherical1d}, equation, mesh)
-  broadcast_volume!(maxspeed, equation, mesh)
-end
-function compute_maxspeed(::Project{:rescaled_spherical1d}, equation, mesh)
-  broadcast_volume!(maxspeed_rescaled_spherical1d, equation, mesh)
-end
-
 function timestep(mesh, P::Project, hrsc::Maybe{HRSC.AbstractHRSC})
   compute_maxspeed(P, P.equation, mesh)
   vmax = get_maxspeed(mesh, P.equation)