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

-### TODO Move include to GRHD.jl or merge with this file?
-include("cons2prim_implementations.jl")
+function make_Equation(eos, prms, m::Mesh)
+  @unpack atm_factor, atm_threshold_factor = prms
+  if m isa Mesh1d
+    return Equation(eos, -1.0, atm_factor, atm_threshold_factor)
+  elseif m isa Mesh2d
+    TODO()
+  else
+    error("unknown mesh type encountered")
+  end
+end
 
 
-module Formulae
-  ### Auxilary functions
+#######################################################################
+#                             "spherical1d"                           #
+#           Spherical symmetry (1d) + spherical coordinates           #
+#######################################################################
 
-  # Lorentz factor
-  Wprim(B, ρ, vr, ϵ, p) = 1.0 / sqrt(1.0 - B^2 * vr^2)
-  Wcons(B, D, Sr, τ, p) = (τ + p + D) / sqrt((τ + p + D)^2 - (Sr / B)^2)
 
-  ### Primitive variables
-  # in terms of conservative ones
+@with_signature function flux_source_spherical1d(eq::Equation)
+  @accepts D, Sr, τ, ϵ, ρ, p, vr
+  @accepts α, ∂αr, βru, ∂βrrdu,
+           γrr, γθθ, γϕϕ, ∂γrrr, ∂γrθθ, ∂γrϕϕ, ∂γθϕϕ,
+           Γrrr, Γrθθ, Γrϕϕ, Γθrθ, Γθϕϕ, Γϕθϕ, Γϕrϕ,
+           Krr, Kθθ, Kϕϕ
+  @accepts flr_D, flr_Sr, flr_τ
+
+  γuurr = 1/γrr
+  γuuθθ = 1/γθθ
+  γuuϕϕ = 1/γϕϕ
+  Γkrk = Γrrr + Γθrθ + Γϕrϕ
+
+  gtt = -α^2 + γrr*βru*βru
+  gtr = γrr*βru
+  grr = γrr
+  guutt = -1/α^2
+  guutr = βru/α^2
+  guurr = γuurr - βru*βru/α^2
+  guuθθ = γuuθθ
+  guuϕϕ = γuuϕϕ
+
+  vru = γuurr*vr
+  v2 = vru*vr
+  W = 1/sqrt(1-v2)
+  ut = W/α
+  ur = W*(vru - βru/α)
+
+  h = 1 + ϵ + p/ρ
+  Tuutt = ρ*h*ut*ut + p*guutt
+  Tuutr = ρ*h*ut*ur + p*guutr
+  Tuurr = ρ*h*ur*ur + p*guurr
+  Tuuθθ = p*guuθθ
+  Tuuϕϕ = p*guuϕϕ
+  Tudtr = Tuutt*gtr + Tuutr*grr
+
+  flr_D  = α*D*(vru - βru/α)
+  flr_Sr = α*( Sr*(vru - βru/α) + p )
+  flr_τ  = α*( τ*(vru - βru/α) + p*vru )
+
+  src_D = - flr_D*Γkrk
+
+  # tmp1 = 1/2*(Tuutt*βru*βru + 2*Tuutr*βru + Tuurr) * (∂γrrr - 2*Γrrr*γrr)
+  tmp1 = 0 # because Dj γik = 0 (Levi-Civita connection)
+  tmp2 = Tudtr*(∂βrrdu + Γrrr*βru)
+  tmp3 = -Tuutt*α*∂αr
+  flθ_Sθ = α*p
+  flϕ_Sϕ = α*p
+  src_Sr = α*(tmp1 + tmp2 + tmp3) - flr_Sr*Γkrk + flr_Sr*Γrrr + flθ_Sθ*Γθrθ + flϕ_Sϕ*Γϕrϕ
+
+  src_τ = Tuutt*(βru*βru*Krr - βru*∂αr) + Tuutr*(2*βru*Krr - ∂αr) + Tuurr*Krr + Tuuθθ*Kθθ + Tuuϕϕ*Kϕϕ
+  src_τ *= α
+  src_τ += - flr_τ*Γkrk
+
+  @returns src_D, src_Sr, src_τ, flr_D, flr_Sr, flr_τ
+end
 
-  ρ(B, D, Sr, τ, p)  = D / Wcons(B, D, Sr, τ, p)
-  vr(B, D, Sr, τ, p) = Sr / (B^2 * (τ + D + p))
-  ϵ(B, D, Sr, τ, p) = (sqrt((τ + p + D)^2 - (Sr/B)^2) - Wcons(B, D, Sr, τ, p) * p) / D - 1
 
-  ### Conservative variables
-  # in terms of primitive variables
+@with_signature function maxspeed(equation::Equation)
+  @accepts vr, ρ, ϵ
+  @accepts α, βru, γrr
+  @unpack eos = equation
+  cs2 = eos(SquareSpeedOfSound, Density, ρ, InternalEnergy, ϵ)
+  cs = sqrt(abs(cs2))
+  γuurr = 1/γrr
+  vru = γuurr * vr
+  v2 = vru*vr
+  λ0 = α*vru - βru
+  W = 1/sqrt(1-v2)
+  λp = α/(1-v2)*(vru*(1-cs2)+cs/W*sqrt((1-v2*cs2)*γuurr-vru^2*(1-cs2)))
+  λm = α/(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
 
-  D(B, ρ, vr, ϵ, p) = ρ * Wprim(B, ρ, vr, ϵ, p)
-  Sr(B, ρ, vr, ϵ, p) = ρ * (1 + ϵ + p / ρ) * Wprim(B, ρ, vr, ϵ, p)^2 * B^2 * vr
-  τ(B, ρ, vr, ϵ, p) = ρ * (1 + ϵ + p / ρ) * Wprim(B, ρ, vr, ϵ, p)^2 - p - ρ * Wprim(B, ρ, vr, ϵ, p)
 
-  ### Fluxes
-  # in terms of primitive variables
+@with_signature function llf(equation::Equation)
+  @accepts        D, flr_D, Sr, flr_Sr, τ, flr_τ, max_v, p
+  @accepts        α, βru, γrr
+  @accepts [bdry] bdry_D, bdry_Sr, bdry_τ, bdry_max_v, bdry_ρ, bdry_vr, bdry_ϵ, bdry_p
+  @accepts [bdry] nx
 
-  fD(B, ρ, vr, ϵ, p) = D(B, ρ, vr, ϵ, p) * vr
-  fSr(B, ρ, vr, ϵ, p) = Sr(B, ρ, vr, ϵ, p) * vr + p
-  fτ(B, ρ, vr, ϵ, p) = (τ(B, ρ, vr, ϵ, p) + p) * vr
+  vmax     = absmax(max_v, bdry_max_v)
+  γuurr    = 1/γrr
+  bdry_vru = γuurr * bdry_vr
 
-  #### Sources
+  bdry_flr_D  = α * bdry_D * (bdry_vru - βru/α)
+  bdry_flr_Sr = α * ( bdry_Sr * (bdry_vru - βru/α) + bdry_p )
+  bdry_flr_τ  = α * ( bdry_τ * (bdry_vru - βru/α) + bdry_p * bdry_vru )
 
-  function sD(r, A, A_r, B, B_r, ρ, vr, ϵ, p)
-    val = - A * B^3 * D(B, ρ, vr, ϵ, p) * vr
-    return val
-  end
+  nflr      = nx*flr_D
+  bdry_nflr = nx*bdry_flr_D
+  nflr_D    = LLF(nflr, bdry_nflr, D, bdry_D, vmax)
 
-  function sSr(r, A, A_r, B, B_r, ρ, vr, ϵ, p)
-    h = 1 + ϵ + p / ρ
-    Sr = Formulae.Sr(B, ρ, vr, ϵ, p)
-    W = Wprim(B, ρ, vr, ϵ, p)
-    v2 = B^2 * vr^2
-    ρhW2 = ρ * h * W^2
-    val = - (ρhW2 - p) * r * B^3 * A_r
-    val += A * B^4 * (r * B_r + B / 3) * (ρhW2 * vr^2 + 3 * p  / B^2)
-    val += - 4 / 3 * A * B^3 * Sr * vr
-    return val
-  end
+  nflr      = nx*flr_Sr
+  bdry_nflr = nx*bdry_flr_Sr
+  nflr_Sr   = LLF(nflr, bdry_nflr, Sr, bdry_Sr, vmax)
 
-  function sτ(r, A, A_r, B, B_r, ρ, vr, ϵ, p)
-    τ = Formulae.τ(B, ρ, vr, ϵ, p)
-    h = 1 + ϵ + p / ρ
-    W = Wprim(B, ρ, vr, ϵ, p)
-    ρhW2 = ρ * h * W^2
-    val = - r * B^3 * ρhW2 * vr * A_r - A * B^3 * (τ + p) * vr
-    return val
-  end
+  nflr      = nx*flr_τ
+  bdry_nflr = nx*bdry_flr_τ
+  nflr_τ    = LLF(nflr, bdry_nflr, τ, bdry_τ, vmax)
+
+  @returns [bdry] nflr_D, nflr_Sr, nflr_τ
 end
 
 
-@with_signature function flux(equation::Equation)
-  @accepts State(sD, sSr, stau), p, rho, vr, eps, A, B, r
-  rAB3     = r * A * B^3
-  flx_sD, flx_sSr, flx_stau = 0, 0, 0
-  try
-    flx_sD   = rAB3 * Formulae.fD(B, rho, vr, eps, p)
-    flx_sSr  = rAB3 * Formulae.fSr(B, rho, vr, eps, p)
-    flx_stau = rAB3 * Formulae.fτ(B, rho, vr, eps, p)
-  catch e
-    println("failed @ r = $r")
-    rethrow(e)
-  end
-  @returns flx_sD, flx_sSr, flx_stau
+@with_signature function bdryllf(equation::Equation)
+  @accepts        D, flr_D, Sr, flr_Sr, τ, flr_τ, max_v, vr, p
+  @accepts        init_D, init_Sr, init_τ, init_vr, init_max_v, init_p
+  @accepts        α, βru, γrr
+  @accepts [bdry] nx
+
+  vmax     = absmax(max_v, init_max_v)
+  γuurr    = 1/γrr
+  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_τ = τ
+
+  bdry_flr_D  = α * bdry_D * (bdry_vru - βru/α)
+  bdry_flr_Sr = α * ( bdry_Sr * (bdry_vru - βru/α) + bdry_p )
+  bdry_flr_τ  = α * ( bdry_τ * (bdry_vru - βru/α) + bdry_p * bdry_vru )
+
+  nflr      = nx*flr_D
+  bdry_nflr = nx*bdry_flr_D
+  nflr_D    = LLF(nflr, bdry_nflr, D, bdry_D, vmax)
+
+  nflr      = nx*flr_Sr
+  bdry_nflr = nx*bdry_flr_Sr
+  nflr_Sr   = LLF(nflr, bdry_nflr, Sr, bdry_Sr, vmax)
+
+  nflr      = nx*flr_τ
+  bdry_nflr = nx*bdry_flr_τ
+  nflr_τ    = LLF(nflr, bdry_nflr, τ, bdry_τ, vmax)
+
+  @returns [bdry] nflr_D, nflr_Sr, nflr_τ
 end
 
-
-@with_signature function speed(equation::Equation)
-  @accepts State(sD, sSr, stau), p, rho, vr, eps, A, B, r
-  # TODO Implement characteristics and derive maxiumum wave speeds.
-  # Might be as simple as adding a missing factor for the auxiliary metric to the
-  # eigenvalues, since ts is the only thing the generalized Valencia formulation adds
-  # to the principal part.
-  # for now we use speed of light
-  max_v = 0.999
-  @returns max_v
+@with_signature function fv_bdry_flux(equation::Equation)
+  @accepts D, Sr, τ, vr, p
+  @accepts init_D, init_Sr, init_τ, init_vr, init_p
+  @accepts α, βru, γrr
+  @accepts [bdry] nx
+
+  γuurr    = 1/γrr
+  init_vru = γuurr * init_vr
+  vru      = γuurr * vr
+
+  # Dirichlet conditions
+  # 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
+
+  nflr_D  = α * bdry_D * (bdry_vru - βru/α)
+  nflr_Sr = α * ( bdry_Sr * (bdry_vru - βru/α) + bdry_p )
+  nflr_τ  = α * ( bdry_τ * (bdry_vru - βru/α) + bdry_p * bdry_vru )
+
+  @returns [bdry] bdry_D, bdry_Sr, bdry_τ, nflr_D, nflr_Sr, nflr_τ
 end
 
 
-# @with_signature function impose_atmosphere(equation::Equation)
-#   @accepts rho, p, eps, v, D, S, tau
-#   @unpack atmosphere = equation
-#
-#   ρatm = atmosphere.rho
-#   patm = 0.0
-#   ϵatm = tau / ρatm
-#   sm = 1.0
-#   Watm = 1.0
-#   vatm = 1.0
-#
-#   if rho < atmosphere.rho * atmosphere.threshold
-#     rho = rhoatm
-#     p   = patm
-#     eps = epsatm
-#     v   = vatm
-#     D   = Formulae.D(rhoatm, vatm, epsatm, patm)
-#     S   = Formulae.S(rhoatm, vatm, epsatm, patm)
-#     tau = Formulae.τ(rhoatm, vatm, epsatm, patm)
-#   end
-#
-#   @returns rho, p, eps, v, D, S, tau
-# end
-
-
-@with_signature function primitives(equation::Equation)
-  @accepts D, Sr, tau, p, B
-  rho = Formulae.ρ(B, D, Sr, tau, p)
-  vr  = Formulae.vr(B, D, Sr, tau, p)
-  eps = Formulae.ϵ(B, D, Sr, tau, p)
-  @returns rho, vr, eps
-end
+#######################################################################
+#                       "rescaled_spherical1d"                        #
+#           Spherical symmetry (1d) + spherical coordinates           #
+#######################################################################
 
 
-@with_signature function prim2cons(equation::Equation)
-  @accepts rho, vr, eps, p, B
-  D   = Formulae.D(B, rho, vr, eps, p)
-  Sr  = Formulae.Sr(B, rho, vr, eps, p)
-  tau = Formulae.τ(B, rho, vr, eps, p)
-  @returns D, Sr, tau
-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
 
-@with_signature function entropy_variables(equation::Equation{EquationOfState.IdealGas})
-  @accepts rho, p, E, flx_E
-  TODO(entropy_variables)
-  # @unpack eos = equation
-  # γ     = eos.gamma
-  # E     = rho / (γ - 1) * log(p / rho^γ)
-  # flx_E = q * E / rho
-  @returns E, flx_E
+  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_rSr = -(τ+p)*r*B3*∂Ar + A*B^2*(r*∂Br+B/3)*(Sr*vru/A^2+3*p) - 4/3*A*B3*Sr*vru
+  src_rSr = -(τ+p)*r*B3*∂Ar + A*B^2*(r*∂Br+B/3)*(Sr*vru/A^2+3*p) - 4/3*A*B3*Sr*vru
+  # src_rτ  = -r*B*Sr*∂Ar - A*B3*(τ+p)*vru
+  src_rτ  = -r*B*Sr*∂Ar/A - 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 sources(equation::Equation)
-  @accepts rho, vr, p, eps, r, A, A_r, B, B_r
-  src_sD   = Formulae.sD(r, A, A_r, B, B_r, rho, vr, eps, p)
-  src_sSr  = Formulae.sSr(r, A, A_r, B, B_r, rho, vr, eps, p)
-  src_stau = Formulae.sτ(r, A, A_r, B, B_r, rho, vr, eps, p)
-  @returns src_sD, src_sSr, src_stau
+@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
+  B3       = B^3
+
+  # Dirichlet conditions
+  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
+
+  # Neumann conditions
+  # bdry_D   = D
+  # bdry_Sr   = Sr
+  # bdry_τ = τ
+
+  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 cons2scaledcons(equations::Equation)
-  @accepts D, Sr, tau, B, r
-  sD   = r * B^3 * D
-  sSr  = r * B^3 * Sr
-  stau = r * B^3 * tau
-  @returns sD, sSr, stau
+@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
 
 
-@with_signature function scaledcons2cons(equations::Equation)
-  @accepts sD, sSr, stau, B, r
-  D   = sD / (r * B^3)
-  Sr  = sSr / (r * B^3)
-  tau = stau / (r * B^3)
-  @returns D, Sr, tau
+@with_signature function fv_bdry_flux_rescaled_spherical1d(equation::Equation)
+  @accepts D, Sr, τ, vr, p
+  @accepts init_D, init_Sr, init_τ, init_vr, init_p
+  @accepts A, B, ∂Br, r
+  @accepts [bdry] nx
+
+  γuurr    = 1/B^2
+  init_vru = γuurr * init_vr
+  vru      = γuurr * vr
+  B3       = B^3
+
+  # # Dirichlet conditions
+  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
+  # bdry_D   = init_D
+  # bdry_Sr  = init_Sr
+  # bdry_τ   = init_τ
+  # bdry_vru = init_vru
+  # bdry_p   = init_p
+
+  @show init_D, D
+
+  nflr_rD  = r*A*B3* bdry_D * bdry_vru
+  nflr_rSr = r*A*B3* (bdry_Sr * bdry_vru + bdry_p)
+  nflr_rτ  = r*A*B3* (bdry_τ + bdry_p) * bdry_vru
+
+  # FIXME We are abusing the bdry_D,Sr,τ buffers to store the rescaled bdry values
+  bdry_D  *= r*B3
+  bdry_Sr *= r*B3
+  bdry_τ  *= r*B3
+
+  @returns [bdry] bdry_D, bdry_Sr, bdry_τ, nflr_rD, nflr_rSr, nflr_rτ
 end
 
 
-"""
-Postprocess points where we divided by zero in `scaledcons2cons`.
-Do so by interpolating value where x = 0 using symmetry/anti-symmetry
-"""
-function postprocess_scaledcons2cons!(cache::Cache, equations::Equation, mesh::Mesh)
+#######################################################################
+#            helpers to resolve coordinate singularities              #
+#            by imposing symmetry conditions at r=0                   #
+#######################################################################
 
-  @unpack x = mesh
-  Npts, K = layout(mesh)
 
-  if isodd(Npts) && isodd(K)
-    error("Cannot yet handle case where r = 0 falls into symmetric cell with even polynomial!")
-  end
+function impose_symmetry!(P::Project{:rescaled_spherical1d}, mesh::Mesh)
 
-  # find points where x coordinate is zero
-  idxs = findall(xi->abs(xi) < 1e-14, mesh.x)
+  @unpack cache = mesh
+  @unpack D, Sr, τ = get_static_variables(cache)
+  L = layout(mesh)
+  Npts, _ = L
 
-  @unpack D, Sr, tau    = get_static_variables(cache)
+  vD  = dg1d.vreshape(D, L)
+  vSr = dg1d.vreshape(Sr, L)
+  vτ  = dg1d.vreshape(τ, L)
 
-  for con in (D, tau) # symmetric variables need to be interpolated
-    mat_con  = reshape(view(con, :), (Npts, K))
-    for idx in idxs
-      i, j = Tuple(idx)
-      mat_con[i,j] = 0 # remove nan, because we are computing inplace
-      Di = mesh.element.D[i,:] # derivative weights for ith point
-      Dii = Di[i]
-      mat_con[i,j] = - dot(Di, mat_con[:,j]) / Dii
-    end
-  end
-  for con in (Sr,) # asymmetric variables need to be set to zero
-    for idx in idxs
-      i, j = Tuple(idx)
-      con[i+(j-1)*Npts] = 0
+  # even 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
+  # odd variables are set to zero at r=0
+  Sr[1] = 0
+
+  return
 end
 
 
-#######################################################################
-#                                 LDG                                 #
-#######################################################################
+function impose_symmetry!(P::Project{:rescaled_spherical1d}, mesh::Mesh{FVElement})
 
+  @unpack cache = mesh
+  @unpack rD = get_dynamic_variables(cache)
+  @unpack D, Sr, τ, r = get_static_variables(cache)
+  L = layout(mesh)
+  Npts, _ = L
 
-@with_signature function ldg_flux(equations::Equation)
-  @accepts State(sD, sSr, stau)
-  flx_ldg_sD, flx_ldg_sSr, flx_ldg_stau = sD, sSr, stau
-  @returns flx_ldg_sD, flx_ldg_sSr, flx_ldg_stau
-end
-
+  @show rD[1], rD[2]
 
-@with_signature function ldg_speed(equations::Equation)
-  @accepts State(sD, sSr, stau)
-  v_ldg_sD, v_ldg_sSr, v_ldg_stau = 1, 1, 1
-  @returns v_ldg_sD, v_ldg_sSr, v_ldg_stau
-end
+  # symmetric variables are interpolated such that derivative vanishes at r=0
+  for var in (D, τ)
+    tmp = 0
+    # use quadratic interpolation around origin
+    var2, var3, r2, r3 = var[2], var[3], r[2], r[3]
+    var[1] = (r3^2*var2-r2^2*var3)/(r3^2-r2^2)
+  end
+  # asymmetric variables are set to zero at r=0
+  Sr[1] = 0
 
+  @show D[1], D[2]
 
-@with_signature function av_flux(equations::Equation)
-  @accepts State(ldg_sD, ldg_sSr, ldg_stau), smoothed_mu
-  flx_ldg_sD, flx_ldg_sSr, flx_ldg_stau =
-    smoothed_mu * ldg_sD, smoothed_mu * ldg_sSr, smoothed_mu * ldg_stau
-  @returns flx_ldg_sD, flx_ldg_sSr, flx_ldg_stau
+  return
 end
 
 
-@with_signature function av_speed(equations::Equation)
-  @accepts State(ldg_sD, ldg_sSr, ldg_stau), smoothed_mu
-  flx_ldg_sD, flx_ldg_sSr, flx_ldg_stau = smoothed_mu, smoothed_mu, smoothed_mu
-  @returns v_flx_ldg_sD, v_flx_ldg_sSr, v_flx_ldg_stau
+function impose_symmetry_sources!(P::Project{:spherical1d}, mesh::Mesh)
+  @unpack r = get_static_variables(mesh.cache)
+  r[1] ≈ 0.0 || return
+  @unpack src_D, src_Sr, src_τ = get_static_variables(mesh.cache)
+  L = layout(mesh)
+  Npts, _ = L
+  v_src_D  = dg1d.vreshape(src_D, L)
+  v_src_τ  = dg1d.vreshape(src_τ, L)
+  # even variables are interpolated such that derivative vanishes at r=0
+  for var in (v_src_D, v_src_τ)
+    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
+  # odd variables are set to zero at r=0
+  src_Sr[1] = 0
+  return
 end