debug EFE formulation

examples/efe_bh/dg.toml

@@ -4,7 +4,7 @@ id = "bh"
 formulation = "z3"
 
 [Mesh]
-range = [ 1.0, 40.0 ]
+range = [ 1.0, 100.0 ]
 n = 5
 k = 64
 basis = "lgl"
@@ -17,13 +17,13 @@ every_iteration = 1
 # every_dt  = 5.0
 # aligned_ts = "$(collect(range(5.0,5000.0,step=5.0)))"
 # aligned_ts      = "$(collect(range(0.2,120.0,step=0.2)))"
-# variables1d = """
-# grr, gθθ, Ar, Drrr, Drθθ, Zr, Krr, Kθθ, α,
-# rhs_grr, rhs_gθθ, rhs_Ar, rhs_Drrr, rhs_Drθθ, rhs_Zr, rhs_Krr, rhs_Kθθ, rhs_α
-# """
 variables1d = """
 grr, gθθ, Ar, Drrr, Drθθ, Zr, Krr, Kθθ, α,
+rhs_grr, rhs_gθθ, rhs_Ar, rhs_Drrr, rhs_Drθθ, rhs_Zr, rhs_Krr, rhs_Kθθ, rhs_α
 """
+# variables1d = """
+# rhs_grr, rhs_gθθ, rhs_Ar, rhs_Drrr, rhs_Drθθ, rhs_Zr, rhs_Krr, rhs_Kθθ, rhs_α,
+# """
 
 [Evolution]
 cfl = 0.5

src/EFE/equation.jl

@@ -49,7 +49,7 @@ end
 
 @with_signature function maxspeed(equation::Equation)
   @accepts α, grr
-  max_v = sqrt(α^2/grr)
+  max_v = α/sqrt(grr)
   @returns max_v
 end
 
@@ -129,18 +129,18 @@ end
 
   vmax     = absmax(max_v, init_max_v)
 
-  # Dirichlet conditions
-  if nx < 0.0 # bdry inside horizon should be purely outflow
-    bdry_grr  = init_grr
-    bdry_gθθ  = init_gθθ
-    bdry_Ar   = init_Ar
-    bdry_Drrr = init_Drrr
-    bdry_Drθθ = init_Drθθ
-    bdry_Zr   = init_Zr
-    bdry_Krr  = init_Krr
-    bdry_Kθθ  = init_Kθθ
-    bdry_α    = init_α
-  else
+  # # Dirichlet conditions
+  # if nx < 0.0 # bdry inside horizon should be purely outflow
+  #   bdry_grr  = init_grr
+  #   bdry_gθθ  = init_gθθ
+  #   bdry_Ar   = init_Ar
+  #   bdry_Drrr = init_Drrr
+  #   bdry_Drθθ = init_Drθθ
+  #   bdry_Zr   = init_Zr
+  #   bdry_Krr  = init_Krr
+  #   bdry_Kθθ  = init_Kθθ
+  #   bdry_α    = init_α
+  # else
     bdry_grr  = -grr  + 2*init_grr
     bdry_gθθ  = -gθθ  + 2*init_gθθ
     bdry_Ar   = -Ar   + 2*init_Ar
@@ -150,7 +150,7 @@ end
     bdry_Krr  = -Krr  + 2*init_Krr
     bdry_Kθθ  = -Kθθ  + 2*init_Kθθ
     bdry_α    = -α    + 2*init_α
-  end
+  # end
 
   bdry_flx_grr  = 0
   bdry_flx_gθθ  = 0

src/EFE/initialdata.jl

@@ -23,11 +23,11 @@ function initialdata_equation!(::Val{:bh}, env, P::Project{:z3}, prms)
   M = 1.0
   rL, rR = extrema(r)
   rL > 0 || error("bh initialdata can't contain coordinate origin")
-  rL/2 <= M || error("radial domain must penetrate the horizon")
+  # rL/2 <= M || error("radial domain must penetrate the horizon")
   K, = mesh.tree.dims
   Npts = mesh.element.Npts
 
-  @. α    = (1 - M/(2*r))/(1+M/(2*r))
+  @. α    = (1-M/(2*r))/(1+M/(2*r))
   @. grr  = (1 + M/(2*r))^4
   @. gθθ  = grr
   Ar     .= dg1d.differentiate(α, mesh) ./ α
@@ -56,6 +56,18 @@ function initialdata_equation!(::Val{:bh}, env, P::Project{:z3}, prms)
   @. init_Sθθ   = Sθθ
   @. init_max_v = max_v
 
+  # This is a simplification of the RHS of the evolution equation for Krr.
+  # This should be zero, or at least converge to zero, but it does not.
+  # So this indicates that something is wrong at least with the ∂tKrr equation,
+  # but there is likely also a problem with ∂tKθθ due to it also being not zero initially.
+  flx = @. α/grr * (Ar + Drrr*2/3)
+  ∂rflx = dg1d.differentiate(flx, mesh)
+  src = @. α*( Drrr*Ar/(3*grr) - Drrr^2*2/(3*grr) )
+  check = @. -flx + src
+  @show extrema(check)
+
+  TODO()
+
 end
 
 

src/EFE/rhs.jl

@@ -70,7 +70,7 @@ function rhs!(mesh::Mesh1d, P::Project{:z3})
   dg1d.interpolate_face_data!(mesh, max_v, bdry_max_v)
 
   broadcast_volume!(flux_source_z3, equation, mesh)
-  impose_symmetry_sources!(P, mesh)
+  # impose_symmetry_sources!(P, mesh)
   broadcast_faces!(llf_z3, equation, mesh)
   if P.prms.problem === :bh
     broadcast_bdry!(bdryllf_bh, equation, mesh)
@@ -78,15 +78,17 @@ function rhs!(mesh::Mesh1d, P::Project{:z3})
     TODO()
   end
 
-  compute_rhs_weak_form!(rhs_grr,  flx_grr,  src_grr,  nflx_grr,  mesh)
-  compute_rhs_weak_form!(rhs_gθθ,  flx_gθθ,  src_gθθ,  nflx_gθθ,  mesh)
+  # compute_rhs_weak_form!(rhs_grr,  flx_grr,  src_grr,  nflx_grr,  mesh)
+  # compute_rhs_weak_form!(rhs_gθθ,  flx_gθθ,  src_gθθ,  nflx_gθθ,  mesh)
+  # TODO Problematic
   compute_rhs_weak_form!(rhs_Krr,  flx_Krr,  src_Krr,  nflx_Krr,  mesh)
-  compute_rhs_weak_form!(rhs_Kθθ,  flx_Kθθ,  src_Kθθ,  nflx_Kθθ,  mesh)
-  compute_rhs_weak_form!(rhs_Ar,   flx_Ar,   src_Ar,   nflx_Ar,   mesh)
-  compute_rhs_weak_form!(rhs_Zr,   flx_Zr,   src_Zr,   nflx_Zr,   mesh)
-  compute_rhs_weak_form!(rhs_α,    flx_α,    src_α,    nflx_α,    mesh)
-  compute_rhs_weak_form!(rhs_Drrr, flx_Drrr, src_Drrr, nflx_Drrr, mesh)
-  compute_rhs_weak_form!(rhs_Drθθ, flx_Drθθ, src_Drθθ, nflx_Drθθ, mesh)
+  # compute_rhs_weak_form!(rhs_Kθθ,  flx_Kθθ,  src_Kθθ,  nflx_Kθθ,  mesh)
+  #
+  # compute_rhs_weak_form!(rhs_Ar,   flx_Ar,   src_Ar,   nflx_Ar,   mesh)
+  # compute_rhs_weak_form!(rhs_Zr,   flx_Zr,   src_Zr,   nflx_Zr,   mesh)
+  # compute_rhs_weak_form!(rhs_α,    flx_α,    src_α,    nflx_α,    mesh)
+  # compute_rhs_weak_form!(rhs_Drrr, flx_Drrr, src_Drrr, nflx_Drrr, mesh)
+  # compute_rhs_weak_form!(rhs_Drθθ, flx_Drθθ, src_Drθθ, nflx_Drθθ, mesh)
 
   return
 end