fix simple wave initial data and disable conservative fixing as it

kills the simple wave propagation

src/SRHD/initialdata.jl

@@ -21,7 +21,7 @@ end
 #######################################################################
 
 
-function initialdata_equation_of_state!(env, P::Project, eos::EquationOfState.IdealGas, prms, mesh::Mesh1d)
+function initialdata_equation_of_state!(env, P::Project, eos, prms, mesh::Mesh1d)
 
   @unpack id = prms
   @unpack cache, mesh = env
@@ -31,14 +31,17 @@ function initialdata_equation_of_state!(env, P::Project, eos::EquationOfState.Id
 
   ϵ0 = @. eos(InternalEnergy, Density, ρ0, Pressure, p0)
 
-  @unpack D, S, tau = get_dynamic_variables(cache)
-  @unpack p         = get_static_variables(cache)
+  @unpack D, S, tau      = get_dynamic_variables(cache)
+  @unpack p, rho, v, eps = get_static_variables(cache)
   @. p   = p0
   @. D   = Formulae.D(ρ0, v0, ϵ0, p0)
   @. S   = Formulae.S(ρ0, v0, ϵ0, p0)
   @. tau = Formulae.τ(ρ0, v0, ϵ0, p0)
+  @. v   = v0
+  @. eps = ϵ0
+  @. rho = ρ0
 
-  dg1d.new_broadcast_volume!(conservative_fixing, P.equation, mesh)
+  # dg1d.new_broadcast_volume!(conservative_fixing, P.equation, mesh)
   dg1d.new_broadcast_volume!(cons2prim_kastaun, P.equation, mesh)
 
   return
@@ -165,9 +168,39 @@ 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 ]
+
+function initialdata_equation_of_state(
+    ::Val{:simple_wave}, eos::EquationOfState.Polytrope, x)
+
+  # from EFL paper which refers to Sebastiano's WENO paper
+  # also these refers from the WENO paper
+  # - [55]: Laim, 1977, in particular eq. II.15 for the speed of sound
+  # - [56]: Anile, 1990, Relativistic Fluids and Magneto-fluids, chapter 5 and sec 5.5.
+
+  xmin, xmax = extrema(x)
+  @toggled_assert isapprox(xmin, -1.5) && isapprox(xmax, 1.5)
+  @toggled_assert eos.K ≈ 100
+  @toggled_assert eos.Gamma ≈ 5/3
+  X, a, Γ, K = 0.3, 0.5, eos.Gamma, eos.K
+
+  # reference state
+  ρref = @. 1.0 * x^0
+  vref = @. 0.0 * x^0
+  ϵref = @. eos(InternalEnergy, Density, ρref)
+  cs0  = @. eos(SpeedOfSound, Density, ρref, InternalEnergy, ϵref)
+
+  # velocity perturbation
+  v0 = @. (abs(x) <= X) * a * sin( pi/2 * (x/X - 1) )^6
+
+  # Liang, 1977, eq. II.15
+  aa = sqrt(Γ-1)
+  bb = aa/2 # = sqrt(Γ-1)/2
+  # speed of sound such that cs(v=0) = cs0 and cs(v=1) = sqrt(Γ-1)
+  cs = @. aa * ( (aa+cs0)/(aa-cs0) * ((1+v0)/(1-v0))^bb - 1 ) / ( (aa+cs0)/(aa-cs0) * ((1+v0)/(1-v0))^bb + 1 )
+
+  # solving cs^2 = ∂p/∂ρ / h, where h is the enthalpy
+  ρ0 = @. (K * Γ * (1/cs^2-1/(Γ-1)))^(1/(1-Γ))
+  p0 = @. eos(Pressure, Density, ρ0)
 
   (ρ0, v0, p0)
 end

src/SRHD/rhs.jl

@@ -10,7 +10,7 @@ function timestep(env, P::Project, hrsc::Maybe{HRSC.AbstractHRSC}, mesh::Mesh1d)
   @unpack equation    = P
   @unpack max_v = get_static_variables(cache)
 
-  dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
+  # dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
   dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
   dg1d.new_broadcast_volume!(new_speed, equation, mesh)
   vmax = dg1d.absolute_maximum(max_v)
@@ -32,7 +32,7 @@ function timestep(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, mesh:
   @unpack cache, mesh = env
   @unpack equation    = P
 
-  dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
+  # dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
   dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
   dg1d.new_broadcast_volume!(new_speed, equation, mesh)
   max_v = get_variable(cache, :max_v)
@@ -119,7 +119,7 @@ function rhs!(env, P::Project, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryc
           bdry_rho, bdry_v, bdry_eps,
           bdry_p, bdry_max_v       = get_bdry_variables(cache)
 
-  dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
+  # 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!(flux, equation, mesh)
@@ -134,8 +134,7 @@ function rhs!(env, P::Project, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryc
   dg1d.interpolate_face_data!(mesh, p, bdry_p)
 
   dg1d.new_broadcast_faces!(llf, equation, mesh)
-  # TODO Implement bdry conditions
-  # broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
+  dg1d.new_broadcast_bdry!(bdryllf, equation, mesh)
 
   compute_rhs_weak_form!(rhs_D,   flx_D,    nflx_D,   mesh)
   compute_rhs_weak_form!(rhs_S,   flx_S,    nflx_S,   mesh)
@@ -166,7 +165,7 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractReconstruction,
   HRSC.reconstruct!(S,    flag, hrsc, isperiodic=isperiodic(bdryconds), limit_with_boundaries=false)
   HRSC.reconstruct!(tau,  flag, hrsc, isperiodic=isperiodic(bdryconds), m=1e-10, limit_with_boundaries=false)
 
-  dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
+  # 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!(flux, equation, mesh)
@@ -216,14 +215,13 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
   ## solve auxiliary equation: q + ∂x u = 0
   dg1d.new_broadcast_faces!(ldg_nflux, equation, mesh)
   dg1d.new_broadcast_bdry!(ldg_bdryllf, equation, mesh)
-  # dg1d.new_broadcast_faces!(ldg_bdryllf, equation, mesh)
   # broadcast_boundaryconditions!(mean_flux, ldg_bdryconds, cache, mesh, 0.0)
   compute_rhs_weak_form!(ldg_D,   D,   nflx_D,   mesh)
   compute_rhs_weak_form!(ldg_S,   S,   nflx_S,   mesh)
   compute_rhs_weak_form!(ldg_tau, tau, nflx_tau, mesh)
 
 
-  dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
+  # dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
 
   ## applying bernsteining where viscosity appeared cures many artificial oscillations
   ## occuring from the AV, but then the AV method appears to be useless...
@@ -274,7 +272,7 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
   # numerical fluxes
   dg1d.new_broadcast_faces!(llf, equation, mesh)
   dg1d.new_broadcast_faces!(av_nflux, equation, mesh)
-  # dg1d.new_broadcast_bdry!(llf, equation, mesh)
+  dg1d.new_broadcast_bdry!(bdryllf, equation, mesh)
 
   compute_rhs_weak_form!(rhs_D,   flx_D,   nflx_D,   mesh)
   compute_rhs_weak_form!(rhs_S,   flx_S,   nflx_S,   mesh)