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/initialdata.jl

 function initialdata!(env, P::Project, prms)
   @unpack id = prms
-  let
-    @unpack r = get_static_variables(env.cache)
-    @. r = env.mesh.x[:]
-  end
   initialdata_equation!(Val(Symbol(id)), env, P, prms)
-  return
-  # initialdata_entropy_variables!(id, P, prms)
-  initialdata_hrsc!(Val(Symbol(id)), env, P, P.hrsc, prms)
-  initialdata_tci!(Val(Symbol(id)), env, P, P.tci, prms)
+  # # initialdata_entropy_variables!(id, P, prms)
+  # initialdata_hrsc!(Val(Symbol(id)), env, P, P.hrsc, prms)
+  # initialdata_tci!(Val(Symbol(id)), env, P, P.tci, prms)
 end
 
 
@@ -17,291 +12,382 @@ end
 #######################################################################
 
 
-function initialdata_equation!(::Val{T}, env, P, parameters) where T
-  TODO(initialdata_equation!)
-end
-
+function initialdata_equation!(::Val{:bondi_accretion}, env, P::Project{:spherical1d}, prms)
 
-function initialdata_equation!(::Val{:bondi_accretion}, env, P, prms)
-  TODO(initialdata_equation!)
-end
+  @unpack mesh = env
+  @unpack equation = P
 
+  id_filename = prms["id_filename"]
+  if !isfile(id_filename)
+    error("GRHD: Failed to locate bondi_accretion initial data at $id_filename")
+  end
 
-function readdata_tov(filename)
-  return h5open(filename, "r") do file
-    g = file["TOV"]
-    Γ = read_attribute(g, "Γ")
-    K = read_attribute(g, "K")
+  data = h5open(id_filename, "r") do file
+    g = file["bondi_accretion"]
+    Γ  = read_attribute(g, "Γ")
+    K  = read_attribute(g, "K")
     ρc = read_attribute(g, "ρc")
-    M = read_attribute(g, "M")
-    ϕ = read(g, "ϕ")
-    p = read(g, "p")
-    m = read(g, "m")
-    R = read(g, "R")
-    r = read(g, "r")
-    (; Γ, K, ρc, M, ϕ, p, m, R, r)
+    rc = read_attribute(g, "rc")
+    M  = read_attribute(g, "M")
+    r  = read(g, "r")
+    ρ  = read(g, "ρ")
+    p  = read(g, "p")
+    ϵ  = read(g, "ϵ")
+    vr = read(g, "vr")
+    (; Γ, K, ρc, rc, M, r, ρ, ϵ, p, vr)
   end
+
+  @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,
+          init_D, init_Sr, init_τ, init_vr, init_p,
+          init_max_v = get_static_variables(mesh.cache)
+  @unpack α, βru, ∂αr, ∂βrrdu,
+          γrr, γθθ, γϕϕ, ∂γrrr, ∂γrθθ, ∂γrϕϕ, ∂γθϕϕ,
+          Γrrr, Γrθθ, Γrϕϕ, Γθrθ, Γθϕϕ, Γϕrϕ, Γϕθϕ,
+          Krr, Kθθ, Kϕϕ = get_static_variables(mesh.cache)
+
+  p   .= interpolate(r, data.r, data.p)
+  ρ   .= interpolate(r, data.r, data.ρ)
+  ϵ   .= interpolate(r, data.r, data.ϵ)
+  vru  = interpolate(r, data.r, data.vr)
+
+  equation.atm_density = maximum(ρ) * equation.atm_factor
+  equation.atm_factor = prms["atm_factor"]
+  equation.atm_threshold = prms["atm_threshold_factor"]
+
+  # to avoid problems with 1/sin(θ) etc., this is also done in initialdata/SphericalAccretion.jl
+  θ = π/2
+  @unpack M = data
+
+  @. α   = sqrt(r/(r+2*M))
+  @. βru = 2*M/(r+2*M) # = β^r
+  @. ∂αr = M / sqrt(r*(2*M+r)^3)
+  @. ∂βrrdu = - 2*M/(2*M+r)^2
+
+  # γ_ij
+  @. γrr = 1 + 2*M/r
+  @. γθθ = r^2
+  @. γϕϕ = r^2 * sin(θ)^2
+
+  # # γijuu = γ^ij
+  # γrruu = @. 1/γrr
+  # γθθuu = @. 1/γθθ
+  # γϕϕuu = @. 1/γϕϕ
+
+  # ∂iγjk := ∂γijk
+  @. ∂γrrr = -2*M/r^2
+  @. ∂γrθθ = 2*r
+  @. ∂γrϕϕ = 2*r*sin(θ)^2
+  @. ∂γθϕϕ = r^2 * sin(θ)*cos(θ)
+
+  # Γijk := Γ^i_jk wrt γ_ij
+  @. Γrrr = -M/(r*(2*M+r))
+  @. Γrθθ = -r^2/(2*M+r)
+  @. Γrϕϕ = -r^2*sin(θ)^2/(2*M+r)
+  @. Γθrθ = 1/r
+  @. Γθϕϕ = -sin(θ)*cos(θ)
+  @. Γϕrϕ = 1/r
+  @. Γϕθϕ = cos(θ)/sin(θ)
+
+  # Kij := K_ij
+  @. Krr = -2*M*(M+r)/sqrt(r^5*(2*M+r))
+  @. Kθθ = 2*M*sqrt(r/(2*M+r))
+  @. Kϕϕ = 2*M*sqrt(r/(2*M+r))*sin(θ)^2
+
+  # # ∂Kijk = ∂_i K_jk
+  # ∂Krrr = @. 2*M*(5*M^2+2*r^2+6*M*r)/sqrt(r^7*(2*M+r)^3)
+  # ∂Krθθ = @. 2*M^2/sqrt(r*(2*M+r)^3)
+  # ∂Krϕϕ = @. ∂Krθθ * sin(θ)^2
+
+  # # verify trace of K
+  # trK_v1 = @. Krr*γrruu + Kθθ*γθθuu + Kϕϕ*γϕϕuu
+  # trK_v2 = @. 2*M*α^3/r^2*(1+3*M/r)
+  # @show(extrema(trK_v1.-trK_v2))
+
+
+  # from initialdata/SphericalAccretion.jl
+  # also see Papadopoulos and Font 1998, arXiv:gr-qc/9803087, eq (5)
+  v2  = @. γrr*vru^2
+  W   = @. 1/sqrt(1-v2)  # = αu^t, Lorentz factor
+  @. vr = γrr * vru # = γ_rr v^r
+
+  h = @. 1 + ϵ + p/ρ
+  @. D  = W * ρ
+  @. Sr = W^2 * ρ * h * vr
+  @. τ  = W^2 * ρ * h - p - D
+
+  broadcast_volume!(cons2prim, equation, mesh)
+  broadcast_volume!(maxspeed, equation, mesh)
+  broadcast_volume!(flux_source_spherical1d, equation, mesh)
+
+  @. init_D     = D
+  @. init_Sr    = Sr
+  @. init_τ     = τ
+  @. init_vr    = vr
+  @. init_p     = p
+  @. init_max_v = max_v
+
+  return
 end
 
 
-function initialdata_equation!(::Val{:tov}, env, P, prms)
+function initialdata_equation!(::Val{:tov}, env, P::Project{:spherical1d}, prms)
 
   filename = prms["id_filename"]
-  if length(filename) == 0
-    error("GRHD: Missing filename for TOV output, include with GRHD.id_filename = 'tov'")
-  end
-  filename = joinpath(get_output_dir(), filename)
-
-  # read TOV result from file
-  # should only contain solution in star interior
-  # output_dir = dg1d.get_output_dir()
-  # data  = TOV.read_data(joinpath(output_dir, "TOV.h5"))
-  data  = readdata_tov(filename)
-  ϕint  = data.ϕ
-  Rint  = data.R
-  rint  = data.r
-  pint  = data.p
-  K     = data.K
-  Γ     = data.Γ
-  M     = data.M
-  Rsurf = Rint[end]
-
-  # interpolate data onto our grid
-  # we don't assume a special grid layout for the loaded data
-
-  @unpack cache, mesh = env
-  @unpack Npts = mesh.element
-
-  intrp_R = vec(mesh.x) # DG coordinates
-  NN = length(intrp_R)
-
-  # determine if grid contains r = 0 point
-  R_ismirrored = intrp_R[1] ≈ -intrp_R[end]
-
-  if !R_ismirrored
-    idx_origin = findall(Ri -> isapprox(Ri, 0), intrp_R)
-    @assert length(idx_origin) > 0 "Non-symmetric grids must include coordinate origin"
+  if !isfile(filename)
+    error("GRHD: Failed to locate TOV initial data at id_filename = $filename")
   end
 
-  # determine our grid setup
-  R_includes_origin = (
-    !R_ismirrored # LGL grids always contain bdry points
-    || (
-      iseven(mesh.K) ||  # origin falls at interface between two cells
-      (isodd(mesh.K) && isodd(mesh.element.Npts)) # origin sits in the center
-                                                  # of a cell with a point in the center
-    )
-  )
-  R_origin_isdoubled = (
-    R_ismirrored && # only happens for mirrored domains
-    iseven(mesh.K)  # and origin falls at interface between two cells
-  )
-  # if r = 0 only covered by one point we have to exclude the boundary
-  # point when mirroring/refelcting back
-  exclude_origin = !R_origin_isdoubled && R_includes_origin
-
-  # on a mirrored grid we only do interpolation on one half and then
-  # mirror/reflect onto the other
-  @show NN
-  @show R_ismirrored
-  @show R_origin_isdoubled
-  @show R_includes_origin
-  if R_ismirrored
-    idx_innermost_R = if R_origin_isdoubled
-      Int64(NN/2) + 1
-    elseif R_includes_origin
-      Int64((NN+1)/2)
-    else
-      Int64(NN/2) + 1
-    end
-    intrp_R = intrp_R[idx_innermost_R:end]
-  end
-  Rmax = maximum(intrp_R) # star's radius
-
-  @info "TOV star surface located at R = $Rsurf"
-
-  #### setup metric functions A, B and derivatives in isotropic coordinates
-
-  # star interior solution
-  Aint = exp.(ϕint)
-  # B := r / R, needs special treatment of coordinate singularity.
-  Bint = [ (ri == 0) ? 1.0 : ri / Ri for (ri, Ri) in zip(rint, Rint) ]
-  # use linear interpolation to obtain value at r = R = 0.
-  Bint[1] = (Bint[3] - Bint[2]) / (Rint[3] - Rint[2]) * (0.0 - Rint[2]) + Bint[2]
-
-  # derivatives
-  dAint = dg1d.FD.central_dx(Aint, Rint)
-  dBint = dg1d.FD.central_dx(Bint, Rint)
-
-  # star exterior solution =^= Schwarzschild solution
-  Npts_ext = 1000 # number of points used for interpolation of exterior solution
-  Rext = collect(LinRange(Rsurf, Rmax, Npts_ext))
-  Aext, Bext = schwarzschild_solution_isotropic_coordinates(Rext, M)
-
-  # merge interior and exterior solutions
-  # note: exclude surface point in one of the arrays when merging
-  R = cat(Rint, Rext[2:end], dims=1)
-  A = cat(Aint, Aext[2:end], dims=1)
-  B = cat(Bint, Bext[2:end], dims=1)
-
-  # derivatives
-  dA = dg1d.FD.central_dx(A, R)
-  dB = dg1d.FD.central_dx(B, R)
-
-  intrp_A   = interpolate(intrp_R[:], R, A)
-  intrp_B   = interpolate(intrp_R[:], R, B)
-  intrp_A_r  = interpolate(intrp_R[:], R, dA)
-  intrp_B_r  = interpolate(intrp_R[:], R, dB)
-
-  if R_ismirrored
-    intrp_A    = mirror(intrp_A, exclude_origin)
-    intrp_B    = mirror(intrp_B, exclude_origin)
-    intrp_A_r  = reflect(intrp_A_r, exclude_origin)
-    intrp_B_r  = reflect(intrp_B_r, exclude_origin)
+  data = h5open(filename, "r") do file
+    g = file["TOV"]
+    Γ  = read_attribute(g, "Γ")
+    K  = read_attribute(g, "K")
+    ρc = read_attribute(g, "ρc")
+    M  = read_attribute(g, "M")
+    rsch = read(g, "r") # schwarzschild radius coordinate
+    riso = read(g, "R") # isotropic radius coordinate
+    ρ   = read(g, "ρ")
+    p   = read(g, "p")
+    ϵ   = read(g, "ϵ")
+    ϕ   = read(g, "ϕ")
+    m   = read(g, "m")
+
+    # metric potentials: ds^2 = -A^2 dt^2 + B^2 γij dx^i dx^j
+    A = @. exp(ϕ)
+    B = @. rsch / riso
+    # B[1] is nan, because rsch = riso = 0 there
+    # use a parabola to interpolate data point at riso = 0
+    # this also ensures that B is symmetric arcoss riso = 0
+    B2, B3, r2, r3 = B[2], B[3], riso[2], riso[3]
+    B[1] = (r3^2*B2-r2^2*B3)/(r3^2-r2^2)
+
+    (; Γ, K, ρc, M, riso, ρ, ϵ, p, A, B)
   end
 
-  @assert all(map(X->all(isfinite.(X)), [intrp_A, intrp_B, intrp_A_r, intrp_B_r])) == true ""
-    "Initial data interpolation failed!"
-
-  @unpack A, B, A_r, B_r = get_static_variables(cache)
-  A   .= intrp_A
-  B   .= intrp_B
-  A_r .= intrp_A_r
-  B_r .= intrp_B_r
-
-  ### setup hydrodynamic variables
   @unpack equation = P
-  # @unpack atmosphere = equation
-  eos = Polytrope(K=K, Gamma=Γ)
-
-  # apply atmosphere treatment
-  pmax = maximum(pint)
-  ρmax = eos(Density, Pressure, pmax)
-  ρatm = equation.atm_factor * ρmax
-  # initialize atmosphere
-  equation.atm_density = ρatm
-  patm = eos(Pressure, Density, ρatm)
-  pint[@.(pint < 0)] .= 0.0 # just to be sure, they will be replace with patm below anyways
-  ρint = @. eos(Density, Pressure, pint)
-  idxatm_treatment = @. ρint < ρatm * equation.atm_threshold
-  pint[idxatm_treatment] .= patm
-  ρint[idxatm_treatment] .= ρatm
-  vint = zeros(size(Rint))
-  ϵint = @. eos(InternalEnergy, Density, ρint)
-
-  # star exterior data
-  ρext = ρatm * ones(size(Rext))
-  pext = @. eos(Pressure, Density, ρext)
-  vext = zeros(size(Rext))
-  ϵext = @. eos(InternalEnergy, Density, ρext)
-
-  # merge interior and exterior solutions
-  p  = cat(pint, pext[2:end], dims=1)
-  ρ  = cat(ρint, ρext[2:end], dims=1)
-  vr = cat(vint, vext[2:end], dims=1)
-  ϵ  = cat(ϵint, ϵext[2:end], dims=1)
-
-  intrp_p  = interpolate(intrp_R[:], R, p)
-  intrp_ρ  = interpolate(intrp_R[:], R, ρ)
-  intrp_vr = interpolate(intrp_R[:], R, vr)
-  intrp_ϵ  = interpolate(intrp_R[:], R, ϵ)
-
-  if R_ismirrored
-    intrp_p  = mirror(intrp_p, exclude_origin)
-    intrp_ρ  = mirror(intrp_ρ, exclude_origin)
-    intrp_vr = reflect(intrp_vr, exclude_origin)
-    intrp_ϵ  = mirror(intrp_ϵ, exclude_origin)
-  end
-
-  @unpack p, vr, rho, eps = get_static_variables(cache)
-  @unpack sD, sSr, stau   = get_dynamic_variables(cache)
-  @unpack D, Sr, tau,     = get_static_variables(cache)
-
-  p .= intrp_p
-  rho .= intrp_ρ
-  vr .= intrp_vr
-  eps .= intrp_ϵ
-
-  broadcast_volume!(prim2cons, equation, cache)
-  broadcast_volume!(cons2scaledcons, equation, cache)
-
-  @assert all(map(X->all(isfinite.(X)),
-    [p, vr, rho, eps, D, Sr, tau, sD, sSr, stau])) == true "Initial data interpolation failed!"
+  @unpack mesh = env
+  @unpack eos = P.equation
+  @assert data.K ≈ eos.K
+  @assert data.Γ ≈ eos.Gamma
+
+  equation.atm_density = data.ρc * equation.atm_factor
+  equation.atm_factor = prms["atm_factor"]
+  equation.atm_threshold = prms["atm_threshold_factor"]
+
+  @unpack D, Sr, τ = get_dynamic_variables(mesh.cache)
+  @unpack ρ, p, ϵ, vr, r, max_v,
+          init_D, init_Sr, init_τ, init_vr, init_p,
+          init_max_v = get_static_variables(mesh.cache)
+  @unpack α, βru, ∂αr, ∂βrrdu,
+          γrr, γθθ, γϕϕ, ∂γrrr, ∂γrθθ, ∂γrϕϕ, ∂γθϕϕ,
+          Γrrr, Γrθθ, Γrϕϕ, Γθrθ, Γθϕϕ, Γϕrϕ, Γϕθϕ,
+          Krr, Kθθ, Kϕϕ = get_static_variables(mesh.cache)
+
+  # need to split initial data into portion that lies inside and outside of star,
+  # because derivatives of the hydro variables are discontinuous at the surface
+  # although metric functions are supposed to be better conditioned there,
+  # we also split them for the interpolation, because of potential contamination
+  # with numerical error from patching
+  rmax = data.riso[end]
+  idx_exterior = findfirst(ri -> ri > rmax, r) |> something
+  @assert idx_exterior > 0
+  interior = 1:idx_exterior-1
+  exterior = idx_exterior:length(mesh)
+
+  rint = view(r, interior)
+  rext = view(r, exterior)
+  A, B, ∂Ar, ∂Br = [ similar(r) for _ in 1:4 ]
+
+  ρ[interior] .= interpolate(rint, data.riso, data.ρ)
+  ϵ[interior] .= interpolate(rint, data.riso, data.ϵ)
+  p[interior] .= interpolate(rint, data.riso, data.p)
+  A[interior] .= interpolate(rint, data.riso, data.A)
+  B[interior] .= interpolate(rint, data.riso, data.B)
+
+  ρ[exterior] .= equation.atm_density
+  p[exterior] .= @. eos(Pressure, Density, ρ[exterior])
+  ϵ[exterior] .= @. eos(InternalEnergy, Density, ρ[exterior])
+  # metric potentials of Schwarzschild solution in isotropic coordinates:
+  #   ds^2 = -A^2 dt^2 + B^2 γij dx^i dx^j
+  @unpack M = data
+  A[exterior] .= @. (1-M/(2*rext))/(1+M/(2*rext))
+  B[exterior] .= @. (1+M/(2*rext))^2
+
+  # BasicInterpolators.jl cannot differentiate ...  so we diff numerically
+  ∂Ar .= dg1d.differentiate(A, mesh)
+  ∂Br .= dg1d.differentiate(B, mesh)
+
+  θ = π/2
+
+  @. α = A
+  @. βru = 0.0
+  @. ∂αr = ∂Ar
+  @. ∂βrrdu = 0
+
+  @. γrr = B^2
+  @. γθθ = B^2*r^2
+  @. γϕϕ = B^2*r^2*sin(θ)^2
+
+  @. ∂γrrr = 2*B*∂Br
+  @. ∂γrθθ = 2*(r*B^2+B*∂Br*r^2)
+  @. ∂γrϕϕ = ∂γrθθ * sin(θ)^2
+  @. ∂γθϕϕ = 0.0
+
+  @. Γrrr = ∂Br/B
+  @. Γrθθ = (-r^2*B*∂Br-r*B^2)/B^2
+  @. Γrϕϕ = Γrθθ*sin(θ)^2
+  @. Γθrθ = (r^2*B*∂Br+r*B^2)/(r^2*B^2)
+  @. Γθϕϕ = -sin(θ)*cos(θ)
+  @. Γϕrϕ = (r^2*B*∂Br+r*B^2)/(r^2*B^2)
+  @. Γϕθϕ = cos(θ)/sin(θ)
+
+  @. Krr = 0.0
+  @. Kθθ = 0.0
+  @. Kϕϕ = 0.0
+
+  @. vr = 0
+  W = 1.0
+  h = @. 1 + ϵ + p/ρ
+  @. D  = W * ρ
+  @. Sr = W^2 * ρ * h * vr
+  @. τ  = W^2 * ρ * h - p - D
+
+  broadcast_volume!(cons2prim, equation, mesh)
+  broadcast_volume!(maxspeed, equation, mesh)
+  broadcast_volume!(flux_source_spherical1d, equation, mesh)
+  impose_symmetry_sources!(P, mesh)
+
+  @. init_D     = D
+  @. init_Sr    = Sr
+  @. init_τ     = τ
+  @. init_vr    = vr
+  @. init_p     = p
+  @. init_max_v = max_v
 
-  return
-end
-
-
-function schwarzschild_solution_isotropic_coordinates(R, M)
-  Rmin = minimum(R)
-  @assert 2 * Rmin >= M
-  A = @. (1 - 0.5 * M / R) / (1 + 0.5 * M / R)
-  B = @. (1 + 0.5 * M / R)^2
-  return A, B
 end
 
 
-function interpolate(x, x_ref, y_ref)
-  x_intrp = LinearInterpolation(x_ref, y_ref, extrapolation_bc=Line())
-  y = [ x_intrp(xi) for xi in x ]
-  return y
-end
+function initialdata_equation!(::Val{:tov}, env, P::Project{:rescaled_spherical1d}, prms)
 
-
-function mirror(A, exclude_bdry::Bool)
-  if exclude_bdry
-    return cat(A[end:-1:2], A, dims=1)
-  else
-    return cat(A[end:-1:1], A, dims=1)
+  filename = prms["id_filename"]
+  if !isfile(filename)
+    error("GRHD: Failed to locate TOV initial data at id_filename = $filename")
   end
-end
 
-
-function reflect(A, exclude_bdry::Bool)
-  if exclude_bdry
-    return cat(-A[end:-1:2], A, dims=1)
-  else
-    return cat(-A[end:-1:1], A, dims=1)
+  data = h5open(filename, "r") do file
+    g = file["TOV"]
+    Γ  = read_attribute(g, "Γ")
+    K  = read_attribute(g, "K")
+    ρc = read_attribute(g, "ρc")
+    M  = read_attribute(g, "M")
+    rsch = read(g, "r") # schwarzschild radius coordinate
+    riso = read(g, "R") # isotropic radius coordinate
+    ρ   = read(g, "ρ")
+    p   = read(g, "p")
+    ϵ   = read(g, "ϵ")
+    ϕ   = read(g, "ϕ")
+    m   = read(g, "m")
+
+    # metric potentials: ds^2 = -A^2 dt^2 + B^2 γij dx^i dx^j
+    A = @. exp(ϕ)
+    B = @. rsch / riso
+    # B[1] is nan, because rsch = riso = 0 there
+    # use a parabola to interpolate data point at riso = 0
+    # this also ensures that B is symmetric arcoss riso = 0
+    B2, B3, r2, r3 = B[2], B[3], riso[2], riso[3]
+    B[1] = (r3^2*B2-r2^2*B3)/(r3^2-r2^2)
+
+    (; Γ, K, ρc, M, riso, ρ, ϵ, p, A, B)
   end
-end
-
 
-function initialdata_equation!(::Val{:atmosphere}, env, P, prms)
-
-  @unpack cache, mesh = env
   @unpack equation = P
-  @unpack Npts = mesh.element
-  @unpack eos = equation
-  @unpack A, B, A_r, B_r = get_static_variables(P.cache)
-  R = vec(mesh.x)
-
-  # minkowski metric
-  A   .= ones(size(R))
-  B   .= ones(size(R))
-  A_r .= zeros(size(R))
-  B_r .= zeros(size(R))
+  @unpack mesh = env
+  @unpack eos = P.equation
+  @assert data.K ≈ eos.K
+  @assert data.Γ ≈ eos.Gamma
+
+  equation.atm_density = data.ρc * equation.atm_factor
+  equation.atm_factor = prms["atm_factor"]
+  equation.atm_threshold = prms["atm_threshold_factor"]
+
+  @unpack rD, rSr, rτ = get_dynamic_variables(mesh.cache)
+  @unpack ρ, p, ϵ, vr, r, max_v,
+          D, Sr, τ,
+          init_D, init_Sr, init_τ,
+          init_vr, init_p, init_max_v = get_static_variables(mesh.cache)
+  @unpack A, ∂Ar, B, ∂Br = get_static_variables(mesh.cache)
+
+  # need to split initial data into portion that lies inside and outside of star,
+  # because derivatives of the hydro variables are discontinuous at the surface
+  # although metric functions are supposed to be better conditioned there,
+  # we also split them for the interpolation, because of potential contamination
+  # with numerical error from patching
+  rmax = data.riso[end]
+  idx_exterior = findfirst(ri -> ri > rmax, r) |> something
+  @assert idx_exterior > 0
+  interior = 1:idx_exterior-1
+  exterior = idx_exterior:length(mesh)
+
+  rint = view(r, interior)
+  rext = view(r, exterior)
+
+  ρ[interior] .= interpolate(rint, data.riso, data.ρ)
+  ϵ[interior] .= interpolate(rint, data.riso, data.ϵ)
+  p[interior] .= interpolate(rint, data.riso, data.p)
+  A[interior] .= interpolate(rint, data.riso, data.A)
+  B[interior] .= interpolate(rint, data.riso, data.B)
+
+  ρ[exterior] .= equation.atm_density
+  p[exterior] .= @. eos(Pressure, Density, ρ[exterior])
+  ϵ[exterior] .= @. eos(InternalEnergy, Density, ρ[exterior])
+  # metric potentials of Schwarzschild solution in isotropic coordinates:
+  #   ds^2 = -A^2 dt^2 + B^2 γij dx^i dx^j
+  @unpack M = data
+  A[exterior] .= @. (1-M/(2*rext))/(1+M/(2*rext))
+  B[exterior] .= @. (1+M/(2*rext))^2
+
+  # BasicInterpolators.jl cannot differentiate ...  so we diff numerically
+  ∂Ar .= dg1d.differentiate(A, mesh)
+  ∂Br .= dg1d.differentiate(B, mesh)
+
+  @. vr = 0
+  W = 1.0
+  h = @. 1 + ϵ + p/ρ
+  @. D  = W * ρ
+  @. Sr = W^2 * ρ * h * vr
+  @. τ  = W^2 * ρ * h - p - D
+  @. rD  = r*B^3*D
+  @. rSr = r*B^3*Sr
+  @. rτ  = r*B^3*τ
+
+  broadcast_volume!(unscale_conservatives, equation, mesh)
+  impose_symmetry!(P, mesh)
+  broadcast_volume!(cons2prim_rescaled_spherical1d, equation, mesh)
+  broadcast_volume!(flux_source_rescaled_spherical1d, equation, mesh)
+  broadcast_volume!(maxspeed_rescaled_spherical1d, equation, mesh)
+
+  @. init_D     = D
+  @. init_Sr    = Sr
+  @. init_τ     = τ
+  @. init_vr    = vr
+  @. init_p     = p
+  @. init_max_v = max_v
 
-  # low density atmosphere
-  # rho_atm = query(prms, "density", default=1e-11)
-  rho_atm = prms["id_atmosphere_density"]
-  @assert rho_atm > 0.0
-  @assert eos isa Polytrope "Require model for cold atmosphere"
-
-  @unpack rho, vr, p, eps = get_static_variables(P.cache)
+end
 
-  NN = length(mesh)
-  @. rho = rho_atm * $ones(NN)
-  @. p   = eos(Pressure, Density, rho, InternalEnergy, eps)
-  @. eps = eos(InternalEnergy, Density, rho, Pressure, p)
-  @. vr  = $zeros(NN)
 
-  @assert all(map(X->all(isfinite.(X)), [rho, p, eps, vr])) == true ""
-    "Invalid initial data encountered!"
+function interpolate(x, x_ref, y_ref)
+  intrp = CubicSplineInterpolator(x_ref, y_ref)
+  return intrp.(x)
+end
 
-  broadcast_volume!(prim2cons, equation, P.cache)
-  broadcast_volume!(cons2scaledcons, equation, P.cache)
 
-  return
+function interpolate_derivative(x, x_ref, y_ref)
+  intrp = CubicSplineInterpolator(x_ref, y_ref)
+  return intrp.(x)
 end