Refactor GRHD project

Major rework of GRHD project to include new broadcasting interface. The goal is to get a TOV evolution with the AV method.

PR contains various refactors:

• initialdata/: We promote this to a separate env in order to add CairoMakie for plotting.
• refactored and improve initialdata/SphericalAccretion.jl solver so that it agrees with various references: solver from bugner's thesis (extract from bamps), plots from two papers from Papadoulos.
• refactored initialdata/TOV.jl to not depend on dg1d.jl anymore, instead use OrindaryDiffEq and CairoMakie for plots
• copy of cons2prim from SRHD. these two versions should be merged eventually, but will do in a separate PR
• implement spherical1d formulation: this is a rewrite of the evolution GRHD equations to work with non-zero radial shift vector. this is needed to evolve the spherical accretion test in cowling approximation where there should be no dynamics.
• implement rescaled_spherical1d formulation: this is a rewrite of the previous formulation based on the generalized Valencia formulation (cf. Montero, Baumgarte paper from 2017); the auxiliary metric is chosen such that coordinate singularities at r=0 in the source terms are cancelled; doesn't work atm
• allow the spherical1d to include r=0 in the domain by using symmetry conditions to evaluate source terms at r=0 (workaround for rescaled_spherical1d not working right now)
• new ref tests: bondi accretion and TOV star
• implement a differentiate! method on Mesh in order to compute derivatives by hand, if needed

src/GRHD/equation.jl

@@ -260,6 +260,176 @@ end
 end
 
 
+@with_signature function flux_source_rescaled_spherical1d(eq::Equation)
+  @accepts D, Sr, τ, ϵ, ρ, p, vr
+  @accepts A, ∂Ar, B, ∂Br, r
+
+  γuurr = 1/B^2
+  vru   = γuurr * vr
+
+  B3 = B^3
+  flr_rD  = r*A*B3 * D * vru
+  flr_rSr = r*A*B3 * (Sr*vru + p)
+  flr_rτ  = r*A*B3 * (τ + p) * vru
+
+  src_rD  = - A*B3 * D * vru
+  src_rSr = -(τ+p)*r*B3*∂Ar + A*B^2*(r*∂Br+B/3)*(Sr*vru+3*p) - 4/3*A*B3*Sr*vru
+  src_rτ  = -r*B*Sr*∂Ar - A*B3*(τ+p)*vru
+
+  @returns src_rD, src_rSr, src_rτ, flr_rD, flr_rSr, flr_rτ
+end
+
+@with_signature function maxspeed_rescaled_spherical1d(equation::Equation)
+  @accepts vr, ρ, ϵ
+  @accepts B, A
+  @unpack eos = equation
+  cs2 = eos(SquareSpeedOfSound, Density, ρ, InternalEnergy, ϵ)
+  cs = sqrt(abs(cs2))
+  γuurr = 1/B^2
+  vru = γuurr * vr
+  v2 = vru*vr
+  λ0 = A*vru
+  W = 1/sqrt(1-v2)
+  λp = A/(1-v2)*(vru*(1-cs2)+cs/W*sqrt((1-v2*cs2)*γuurr-vru^2*(1-cs2)))
+  λm = A/(1-v2)*(vru*(1-cs2)-cs/W*sqrt((1-v2*cs2)*γuurr-vru^2*(1-cs2)))
+  max_v = absmax(λm, λ0, λp)
+  @returns max_v
+end
+
+@with_signature function llf_rescaled_spherical1d(equation::Equation)
+  @accepts        D, flr_rD, Sr, flr_rSr, τ, flr_rτ, max_v, p
+  @accepts        A, B, r
+  @accepts [bdry] bdry_D, bdry_Sr, bdry_τ, bdry_max_v, bdry_vr, bdry_p
+  @accepts [bdry] nx
+
+  vmax     = absmax(max_v, bdry_max_v)
+  γuurr    = 1/B^2
+  bdry_vru = γuurr * bdry_vr
+
+  B3 = B^3
+  rD       = r*B3*D
+  rSr      = r*B3*Sr
+  rτ       = r*B3*τ
+  bdry_rD  = r*B3*bdry_D
+  bdry_rSr = r*B3*bdry_Sr
+  bdry_rτ  = r*B3*bdry_τ
+
+  bdry_flr_rD  = r*A*B3 * bdry_D * bdry_vru
+  bdry_flr_rSr = r*A*B3 * ( bdry_Sr * bdry_vru + bdry_p )
+  bdry_flr_rτ  = r*A*B3 * ( (bdry_τ + bdry_p) * bdry_vru )
+
+  nflr      = nx*flr_rD
+  bdry_nflr = nx*bdry_flr_rD
+  nflr_rD   = LLF(nflr, bdry_nflr, rD, bdry_rD, vmax)
+
+  nflr      = nx*flr_rSr
+  bdry_nflr = nx*bdry_flr_rSr
+  nflr_rSr  = LLF(nflr, bdry_nflr, rSr, bdry_rSr, vmax)
+
+  nflr      = nx*flr_rτ
+  bdry_nflr = nx*bdry_flr_rτ
+  nflr_rτ   = LLF(nflr, bdry_nflr, rτ, bdry_rτ, vmax)
+
+  @returns [bdry] nflr_rD, nflr_rSr, nflr_rτ
+end
+
+@with_signature function bdryllf_rescaled_spherical1d(equation::Equation)
+  @accepts        D, flr_rD, Sr, flr_rSr, τ, flr_rτ, max_v, vr, p
+  @accepts        init_D, init_Sr, init_τ, init_vr, init_max_v, init_p
+  @accepts        A, B, r
+  @accepts [bdry] nx
+
+  vmax     = absmax(max_v, init_max_v)
+  γuurr    = 1/B^2
+  init_vru = γuurr * init_vr
+  vru      = γuurr * vr
+
+  # enforce dirichlet/neumann conditions following stability analysis
+  # see dg-notes/tex/evm.pdf
+  # Hesthaven, Numerical Methods for Conservation Laws uses same logic
+
+  # Dirichlet conditions
+  # ## if block for hand-tweaked outflow condition for bondi accretion at inner radius
+  # if nx > 0
+    bdry_D   = -D   + 2*init_D
+    bdry_Sr  = -Sr  + 2*init_Sr
+    bdry_τ   = -τ   + 2*init_τ
+    bdry_vru = -vru + 2*init_vru
+    bdry_p   = -p   + 2*init_p
+  # else
+  #   bdry_D   = D
+  #   bdry_Sr  = Sr
+  #   bdry_τ   = τ
+  #   bdry_vru = vru
+  #   bdry_p   = p
+  # end
+
+  # Neumann conditions
+  # bdry_D   = D
+  # bdry_Sr   = Sr
+  # bdry_τ = τ
+
+  B3 = B^3
+  rD       = r*B3*D
+  rSr      = r*B3*Sr
+  rτ       = r*B3*τ
+  bdry_rD  = r*B3*bdry_D
+  bdry_rSr = r*B3*bdry_Sr
+  bdry_rτ  = r*B3*bdry_τ
+
+  bdry_flr_rD  = r*A*B3 * bdry_D * bdry_vru
+  bdry_flr_rSr = r*A*B3 * ( bdry_Sr * bdry_vru + bdry_p )
+  bdry_flr_rτ  = r*A*B3 * ( (bdry_τ + bdry_p) * bdry_vru )
+
+  nflr      = nx*flr_rD
+  bdry_nflr = nx*bdry_flr_rD
+  nflr_rD   = LLF(nflr, bdry_nflr, rD, bdry_rD, vmax)
+
+  nflr      = nx*flr_rSr
+  bdry_nflr = nx*bdry_flr_rSr
+  nflr_rSr  = LLF(nflr, bdry_nflr, rSr, bdry_rSr, vmax)
+
+  nflr      = nx*flr_rτ
+  bdry_nflr = nx*bdry_flr_rτ
+  nflr_rτ   = LLF(nflr, bdry_nflr, rτ, bdry_rτ, vmax)
+
+  @returns [bdry] nflr_rD, nflr_rSr, nflr_rτ
+end
+
+@with_signature function unscale_conservatives(equation::Equation)
+  @accepts rD, rSr, rτ
+  @accepts r, B
+  B3 = B^3
+  D  = rD/(r*B3)
+  Sr = rSr/(r*B3)
+  τ  = rτ/(r*B3)
+  @returns D, Sr, τ
+end
+
+function impose_symmetry!(equations::Equation, mesh::Mesh)
+
+  @unpack cache = mesh
+  @unpack D, Sr, τ = get_static_variables(cache)
+  L = layout(mesh)
+  Npts, _ = L
+
+  vD  = dg1d.vreshape(D, L)
+  vSr = dg1d.vreshape(Sr, L)
+  vτ  = dg1d.vreshape(τ, L)
+
+  # symmetric variables are interpolated such that derivative vanishes at r=0
+  for var in (vD, vτ)
+    tmp = 0
+    for j = 2:Npts
+      tmp += mesh.element.D[1,j] * var[j,1]
+    end
+    var[1,1] = -tmp/mesh.element.D[1,1]
+  end
+  # asymmetric variables are set to zero at r=0
+  Sr[1] = 0
+
+  return
+end
 
 # @with_signature function flux(equation::Equation)
 #   @accepts sD, sSr, stau, p, rho, vr, eps, A, B, r