compare residuals from another paper

src/EFE/initialdata.jl

@@ -28,8 +28,9 @@ function initialdata_equation!(::Val{:bh}, env, P::Project{:z3}, prms)
   Npts = mesh.element.Npts
 
   @. α    = (1-M/(2*r))/(1+M/(2*r))
-  @. grr  = (1 + M/(2*r))^4
-  @. gθθ  = grr
+  ψ4      = @. (1+M/(2*r))^4
+  @. grr  = ψ4
+  @. gθθ  = ψ4
   Ar     .= dg1d.differentiate(α, mesh) ./ α
   Drrr   .= dg1d.differentiate(grr, mesh) ./ (2*grr)
   Drθθ   .= dg1d.differentiate(gθθ, mesh) ./ (2*gθθ)
@@ -63,9 +64,29 @@ function initialdata_equation!(::Val{:bh}, env, P::Project{:z3}, prms)
   flx = @. α/grr * (Ar + Drrr*2/3)
   ∂rflx = dg1d.differentiate(flx, mesh)
   src = @. α*( Drrr*Ar/(3*grr) - Drrr^2*2/(3*grr) )
-  check = @. -flx + src
+  check = @. -∂rflx + src
   @show extrema(check)
 
+  # lets check the version from arXiv:gr-qc/9805084v1 21 May 1998, eq (19)
+  Drrr   .= dg1d.differentiate(grr, mesh) ./ grr
+  Drθθ   .= dg1d.differentiate(gθθ, mesh) ./ gθθ
+  DDrθθ = @. Drθθ-1/r
+  Lr = Ar
+  C = @. α/sqrt(grr)
+  Vr = @. 2*DDrθθ
+  flx = @. -2*C * (3/4*Vr - DDrθθ + Lr)
+  ∂rflx = dg1d.differentiate(flx, mesh)
+  src = @. C * ( Vr^2/16 + DDrθθ*(Vr/2-DDrθθ) + 2*Lr*DDrθθ
+                -3/2*Lr*Vr + 2/r*(4*DDrθθ-Vr/2-(Drrr+2*DDrθθ)) )
+  residual = @. ∂rflx - src
+  @show extrema(residual)
+
+  flx = @. - C*(DDrθθ-Vr/2)
+  ∂rflx = dg1d.differentiate(flx, mesh)
+  src = @. C * ( Lr*Vr/2 + DDrθθ^2/2 - Vr^2/16 + 1/r*(2*Lr+4*DDrθθ-(Drrr+2*DDrθθ)) )
+  residual = @. ∂rflx - src
+  @show extrema(residual)
+
   TODO()
 
 end