src/GRHD/initialdata.jl

@@ -55,8 +55,8 @@ function initialdata_equation!(::Val{:bondi_accretion}, env, P::Project{:spheric
     (; Γ, K, ρc, rc, M, r, ρ, ϵ, p, vr)
   end
 
-  # @show data.K - equation.eos.K
-  # @show data.Γ - equation.eos.Gamma
+  @assert data.K ≈ equation.eos.K
+  @assert data.Γ ≈ equation.eos.Gamma
 
   @unpack D, Sr, τ = get_dynamic_variables(mesh.cache)
   @unpack ρ, p, ϵ, vr, r, max_v,

src/mesh.jl

@@ -404,6 +404,18 @@ norm_L2(u, mesh::Mesh) = sqrt(broken_inner_product(u, u, mesh))
 norm_L1(u, mesh::Mesh) = broken_inner_product(abs.(u), ones(Float64, size(u)), mesh)
 
 
+function differentiate!(du, u, mesh::Mesh1d)
+  for (duk,uk) in zip(eachcell(mesh,du),eachcell(mesh,u))
+    differentiate!(duk, uk, mesh.element)
+  end
+end
+function differentiate(u, mesh)
+  du = similar(u)
+  differentiate!(du, u, mesh)
+  return du
+end
+
+
 """
     find_nearest(xi, mesh::Mesh)
 

src/spectralelement.jl

@@ -93,3 +93,12 @@ end
     return GLGL.integrate(el.w, el.v, el.D, u)
   end
 end
+
+
+@inline function differentiate(u, el::SpectralElement)
+  return el.D * u
+end
+
+@inline function differentiate!(du, u, el::SpectralElement)
+  mul!(du, el.D, u)
+end