wip -- perhaps slopy merged

src/SRHD/callbacks.jl

@@ -93,9 +93,13 @@ function callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc::HRSC.Abstract
     max_v[k] = maximum(view(mat_max_v, :, k))
   end
 
+  # HRSC.compute_viscosity!(
+  #     mu, tmp_mu,
+  #     (D,S,tau), cellmax_v,
+  #     hrsc)
   HRSC.compute_viscosity!(
-      mu, tmp_mu,
-      (D,S,tau), cellmax_v,
+      mu,
+      D, cellmax_v,
       hrsc)
 
 end

src/SRHD/cons2prim.jl

@@ -57,15 +57,15 @@
   #   fixed = true
   # end
 
-  h0 = 1e-3
-  S2 = S^2 / B^2
-  r2 = S2 / D^2
-  z02 = r2 / h0^2
-  v02 = z02 / (1.0 + z02)
-  max_v2 = 0.999
-  if v02 > max_v2
-    S = sign(S) * sqrt(max_v2/(1-max_v2) * h0^2 * D^2)
-  end
+  # h0 = 1e-3
+  # S2 = S^2 / B^2
+  # r2 = S2 / D^2
+  # z02 = r2 / h0^2
+  # v02 = z02 / (1.0 + z02)
+  # max_v2 = 0.999
+  # if v02 > max_v2
+  #   S = sign(S) * sqrt(max_v2/(1-max_v2) * h0^2 * D^2)
+  # end
 
   flag_consfixed = Float64(fixed)
 
@@ -100,6 +100,7 @@ end
   q = tau / D
   S2 = (S/B)^2
   Si = S
+  # @assert sign(Si) > 0
   ri = Si / D
   # r = sqrt(S2) / D # =^= r in paper, r already denotes radius here
   z0i = ri / h0
@@ -194,6 +195,7 @@ end
 
         # compute thermodynamic vars
         phat = eos(Pressure, Density, rhohat, InternalEnergy, epshat)
+        # phat = eos(Pressure, Density, rhohat)
         # hhat = eos(Enthalpy, Density, rhohat, InternalEnergy, epshat)
         ahat = phat / (rhohat * (1.0 + epshat))
 
@@ -230,19 +232,20 @@ end
     # eq (15) [Kastaun] : rho*W = D
     rhohW2 = D / mu
 
-    # regularize velocity according to Dumbser et al arxiv:2307.06629
-    # eqs (82)-(84)
-    y, y0 = rhohW2, 1e-4
-    if y > y0
-      # v = S / mu
-      # compute v accoriding to eq (68) or (40) [Kastaun]
-      v = sign(Si) * min(mu * r, v0)
-    else
-      # @warn "limiting v"
-      # limit velocity ratio v/S
-      fy = 2.0 * (y/y0)^3 - 3.0 * (y/y0)^2 + 1.0
-      v = Si * y / (y^2 + fy * 5e-9)
-    end
+    v = sign(Si) * min(mu * r, v0)
+    # # regularize velocity according to Dumbser et al arxiv:2307.06629
+    # # eqs (82)-(84)
+    # y, y0 = rhohW2, 1e-4
+    # if y > y0
+    #   # v = S / mu
+    #   # compute v accoriding to eq (68) or (40) [Kastaun]
+    #   v = sign(Si) * min(mu * r, v0)
+    # else
+    #   # @warn "limiting v"
+    #   # limit velocity ratio v/S
+    #   fy = 2.0 * (y/y0)^3 - 3.0 * (y/y0)^2 + 1.0
+    #   v = Si * y / (y^2 + fy * 5e-9)
+    # end
 
     # eq (68) or (40)
     # v = min(mu * r, v0)
@@ -273,6 +276,7 @@ end
 
     # compute thermodynamic vars
     p = eos(Pressure, Density, rho, InternalEnergy, eps)
+    # p = eos(Pressure, Density, rho)
     if p <= 0
       error((rho,eps))
     end

src/SRHD/equation.jl

@@ -86,14 +86,20 @@ end
   nflx      = nx*flx_D
   bdry_nflx = nx*bdry_flx_D
   nflx_D    = LLF(nflx, bdry_nflx, D, bdry_D, vmax)
+  # nflx_D    = 0.5 * (nflx + bdry_nflx)
+  # display((flx_D, bdry_flx_D, nx, nflx, bdry_nflx, nflx_D))
 
   nflx      = nx*flx_S
   bdry_nflx = nx*bdry_flx_S
   nflx_S    = LLF(nflx, bdry_nflx, S, bdry_S, vmax)
+  # nflx_S    = 0.5 * (nflx + bdry_nflx)
+  # display((flx_S, bdry_flx_S, nx, nflx, bdry_nflx, nflx_S))
 
   nflx      = nx*flx_tau
   bdry_nflx = nx*bdry_flx_tau
   nflx_tau  = LLF(nflx, bdry_nflx, tau, bdry_tau, vmax)
+  # nflx_tau    = 0.5 * (nflx + bdry_nflx)
+  # display((flx_D, bdry_flx_D, nx, nflx, bdry_nflx, nflx_D))
 
   @returns [bdry] nflx_D, nflx_S, nflx_tau
 end
@@ -180,8 +186,8 @@ end
 
 @with_signature [legacy=false] function av_flux(equations::Equation)
   @accepts flx_D, flx_S, flx_tau, ldg_D, ldg_S, ldg_tau, smoothed_mu
-  flx_D += smoothed_mu * ldg_D
-  flx_S += smoothed_mu * ldg_S
+  flx_D   += smoothed_mu * ldg_D
+  flx_S   += smoothed_mu * ldg_S
   flx_tau += smoothed_mu * ldg_tau
   @returns flx_D, flx_S, flx_tau
 end

src/SRHD/initialdata.jl

@@ -168,6 +168,13 @@ function initialdata_equation_of_state(
 
   xl, xr = -0.5, 0.5
 
+  ρ0 = [ (xi < xl || xi > xr) ? 1.0  :  0.125  for xi in x ]
+  v0 = [ (xi < xl || xi > xr) ? 0.0  :  0.0    for xi in x ]
+  p0 = [ (xi < xl || xi > xr) ? 1.0  :  0.1    for xi in x ]
+
+  (ρ0, v0, p0)
+end
+
 
 function initialdata_equation_of_state(
     ::Val{:simple_wave}, eos::EquationOfState.Polytrope, x)

src/SRHD/rhs.jl

@@ -245,8 +245,8 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
   # HRSC.reconstruct!(S,    flag, bernstein, isperiodic=true, limit_with_boundaries=false)
   # HRSC.reconstruct!(tau,  flag, bernstein, isperiodic=true, m=1e-10, limit_with_boundaries=false)
 
-  dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
-  dg1d.new_broadcast_volume!(new_speed, equation, mesh)
+  # dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
+  # dg1d.new_broadcast_volume!(new_speed, equation, mesh)
 
   dg1d.interpolate_face_data!(mesh, D, bdry_D)
   dg1d.interpolate_face_data!(mesh, S, bdry_S)
@@ -258,8 +258,8 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
   dg1d.interpolate_face_data!(mesh, p, bdry_p)
 
   # # # dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
-  # dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
-  # dg1d.new_broadcast_volume!(new_speed, equation, mesh)
+  dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
+  dg1d.new_broadcast_volume!(new_speed, equation, mesh)
 
   dg1d.new_broadcast_volume!(flux, equation, mesh)
   dg1d.new_broadcast_volume!(av_flux, equation, mesh)

src/mesh.jl

@@ -12,6 +12,11 @@ struct Mesh{N_Dim,N_Sides} <: AbstractMesh
   extends::NTuple{N_Dim,Tuple{Float64,Float64}} # total extend of mesh
   element::SpectralElement
   cache::Cache
+  # TODO Shall we replace these with start:stride:stride. Will require change in @with_signature
+  # to loop over ranges. start should include also cell offsets, because then these
+  # vectors don't have to be as long as the number of cells we have (which we require otherwise
+  # to find the right offset from mesh.offsets)
+  # TODO This should be benchmarked!
   bulkfaceindices::Vector{Int64}
   bulkbdryindices::Vector{Int64}
   faceindices::Vector{Int64}