GRHD: more specific boundry conditions [only incorporated for...

GRHD: more specific boundry conditions [only incorporated for FV+slope_limiter] (!179)

This adds new options under `GRHD.bc`:
- `tov_reflective_origin`: impose reflection symmetry at coordinate origin, impose vacuum values at outer bdry
- `tov_symmetric_domain`: impose vacuum values at both domain boundary ends

Also streamlined the setup of the bdry values for FV.
While doing so I realized that the inflow bdry in the Bondi test was wrongly treated.
There we cannot extrapolate interior data, nor just copy what we have in the interior, unless we are very far out (which in my test we aren't).
Instead, what we do now there is use the initial data, which is the analytical solution and extrapolate that.
This improves the behavior of the solution to no longer drift after the first time step. Nice.

src/GRHD/GRHD.jl

@@ -62,11 +62,15 @@ dg1d.@parameters GRHD begin
 
   """
   Available options:
-  - `periodic`
+  - `tov_reflective_origin`: for `id = tov` use reflection symmetry across coordinate origin;
+    requires `range = [ 0, R ]` where `R` shall be larger than the star's radius;
+    at the outer boundary we impose vacuum variables
+  - `tov_symmetric_domain`: for `id = tov` to impose vacuum values on both ends;
+    requires `range = [ -R, R ]` where `R` shall be larger than the star's radius
   - `from_id`
   """
-  bc = "periodic"
-  @check bc in [ "periodic", "from_id" ]
+  bc = "from_id"
+  @check bc in [ "tov_reflective_origin", "tov_symmetric_domain", "from_id" ]
 
   """
   Available options:

src/GRHD/initialdata.jl

@@ -77,7 +77,7 @@ function initialdata_equation!(::Val{:bondi_accretion}, env, P::Project{:valenci
 
   @unpack rD, rSr, rτ = get_dynamic_variables(mesh.cache)
   @unpack D, Sr, τ, ρ, p, ϵ, vr, r, max_v,
-          init_D, init_Sr, init_τ, init_vr, init_p,
+          init_D, init_Sr, init_τ, init_vr, init_p, init_ρ, init_ϵ,
           init_max_v = get_static_variables(mesh.cache)
   @unpack α, βru, ∂αr, ∂βrrdu, sqrt_detγ,
           γrr, γθθ, γϕϕ, ∂γrrr, ∂γrθθ, ∂γrϕϕ, ∂γθϕϕ,
@@ -175,6 +175,8 @@ function initialdata_equation!(::Val{:bondi_accretion}, env, P::Project{:valenci
   @. init_vr    = vr
   @. init_p     = p
   @. init_max_v = max_v
+  @. init_ρ     = ρ
+  @. init_ϵ     = ϵ
 
   return
 end
@@ -187,7 +189,7 @@ function initialdata_equation!(::Val{:bondi_accretion}, env, P::Project{:spheric
 
   @unpack D, Sr, τ = get_dynamic_variables(mesh.cache)
   @unpack r, max_v,
-          init_D, init_Sr, init_τ, init_vr, init_p,
+          init_D, init_Sr, init_τ, init_vr, init_p, init_ρ, init_ϵ,
           init_max_v = get_static_variables(mesh.cache)
   @unpack α, βru, ∂αr, ∂βrrdu,
           γrr, γθθ, γϕϕ, ∂γrrr, ∂γrθθ, ∂γrϕϕ, ∂γθϕϕ,
@@ -295,6 +297,8 @@ function initialdata_equation!(::Val{:bondi_accretion}, env, P::Project{:spheric
   @. init_vr    = vr
   @. init_p     = p
   @. init_max_v = max_v
+  @. init_ρ     = ρ
+  @. init_ϵ     = ϵ
 
   return
 end
@@ -307,7 +311,7 @@ function initialdata_equation!(::Val{:bondi_accretion_infall}, env, P::Project{:
 
   @unpack D, Sr, τ = get_dynamic_variables(mesh.cache)
   @unpack r, max_v,
-          init_D, init_Sr, init_τ, init_vr, init_p,
+          init_D, init_Sr, init_τ, init_vr, init_p, init_ρ, init_ϵ,
           init_max_v = get_static_variables(mesh.cache)
   @unpack α, βru, ∂αr, ∂βrrdu,
           γrr, γθθ, γϕϕ, ∂γrrr, ∂γrθθ, ∂γrϕϕ, ∂γθϕϕ,
@@ -418,6 +422,8 @@ function initialdata_equation!(::Val{:bondi_accretion_infall}, env, P::Project{:
   @. init_vr    = vr
   @. init_p     = p
   @. init_max_v = max_v
+  @. init_ρ     = ρ
+  @. init_ϵ     = ϵ
 
   let ref_ρ=ρ
     @unpack ρ = get_static_variables(mesh)
@@ -486,7 +492,7 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:valencia1d}, prms)
 
   @unpack rD, rSr, rτ   = get_dynamic_variables(mesh.cache)
   @unpack D, Sr, τ, ρ, p, ϵ, vr, r, max_v,
-          init_D, init_Sr, init_τ, init_vr, init_p,
+          init_D, init_Sr, init_τ, init_vr, init_p, init_ρ, init_ϵ,
           init_max_v    = get_static_variables(mesh.cache)
   @unpack α, βru, ∂αr, ∂βrrdu, sqrt_detγ,
           γrr, γθθ, γϕϕ, ∂γrrr, ∂γrθθ, ∂γrϕϕ, ∂γθϕϕ,
@@ -652,6 +658,8 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:valencia1d}, prms)
   @. init_vr    = vr
   @. init_p     = p
   @. init_max_v = max_v
+  @. init_ρ     = ρ
+  @. init_ϵ     = ϵ
 
 end
 
@@ -672,7 +680,7 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:spherical1d}, prms)
 
   @unpack D, Sr, τ = get_dynamic_variables(mesh.cache)
   @unpack ρ, p, ϵ, vr, r, max_v,
-          init_D, init_Sr, init_τ, init_vr, init_p,
+          init_D, init_Sr, init_τ, init_vr, init_p, init_ρ, init_ϵ,
           init_max_v = get_static_variables(mesh.cache)
   @unpack α, βru, ∂αr, ∂βrrdu,
           γrr, γθθ, γϕϕ, ∂γrrr, ∂γrθθ, ∂γrϕϕ, ∂γθϕϕ,
@@ -830,6 +838,8 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:spherical1d}, prms)
   @. init_vr    = vr
   @. init_p     = p
   @. init_max_v = max_v
+  @. init_ρ     = ρ
+  @. init_ϵ     = ϵ
 
 end
 
@@ -851,7 +861,7 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:rescaled_spherical1
   @unpack rD, rSr, rτ = get_dynamic_variables(mesh.cache)
   @unpack ρ, p, ϵ, vr, r, max_v,
           D, Sr, τ,
-          init_D, init_Sr, init_τ,
+          init_D, init_Sr, init_τ, init_ρ, init_ϵ,
           init_vr, init_p, init_max_v = get_static_variables(mesh.cache)
   @unpack A, ∂Ar, B, ∂Br, γrr = get_static_variables(mesh.cache)
 
@@ -980,6 +990,8 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:rescaled_spherical1
   @. init_vr    = vr
   @. init_p     = p
   @. init_max_v = max_v
+  @. init_ρ     = ρ
+  @. init_ϵ     = ϵ
 
 end
 

src/GRHD/rhs.jl

@@ -375,31 +375,65 @@ function step_fv_slope_limiter!(mesh, P::Project{:spherical1d}, state)
   dl = L/K
 
   id = P.prmsdb["GRHD"]["id"]
+  bc = P.prmsdb["GRHD"]["bc"]
   if id == "bondi_accretion"
-    @unpack r = get_static_variables(mesh)
-    dr = r[2]-r[1]
-    bdry_γrr_l, bdry_γrr_ll = extrapolate_quadratic(γrr, r, 1:3, r[1]-dr, r[1]-2*dr)
-    bdry_α_l,   bdry_α_ll   = extrapolate_quadratic(α,   r, 1:3, r[1]-dr, r[1]-2*dr)
-    bdry_βru_l, bdry_βru_ll = extrapolate_quadratic(βru, r, 1:3, r[1]-dr, r[1]-2*dr)
-    bdry_ρ_l,   bdry_ρ_ll   = extrapolate_quadratic(ρ,   r, 1:3, r[1]-dr, r[1]-2*dr)
-    bdry_vr_l,  bdry_vr_ll  = extrapolate_quadratic(vr,  r, 1:3, r[1]-dr, r[1]-2*dr)
-    bdry_ϵ_l,   bdry_ϵ_ll   = extrapolate_quadratic(ϵ,   r, 1:3, r[1]-dr, r[1]-2*dr)
+    if bc == "from_id"
+      @unpack init_ρ, init_vr, init_ϵ, r = get_static_variables(mesh)
+      dr = r[2]-r[1]
+      bdry_γrr_l, bdry_γrr_ll = extrapolate_quadratic(γrr,     r, 1:3,   r[1]-dr, r[1]-2*dr)
+      bdry_α_l,   bdry_α_ll   = extrapolate_quadratic(α,       r, 1:3,   r[1]-dr, r[1]-2*dr)
+      bdry_βru_l, bdry_βru_ll = extrapolate_quadratic(βru,     r, 1:3,   r[1]-dr, r[1]-2*dr)
+      bdry_ρ_l,   bdry_ρ_ll   = extrapolate_quadratic(ρ,       r, 1:3,   r[1]-dr, r[1]-2*dr)
+      bdry_vr_l,  bdry_vr_ll  = extrapolate_quadratic(vr,      r, 1:3,   r[1]-dr, r[1]-2*dr)
+      bdry_ϵ_l,   bdry_ϵ_ll   = extrapolate_quadratic(ϵ,       r, 1:3,   r[1]-dr, r[1]-2*dr)
+      bdry_γrr_r, bdry_γrr_rr = extrapolate_quadratic(γrr,     r, K-2:K, r[K]+dr, r[K]+2*dr)
+      bdry_α_r,   bdry_α_rr   = extrapolate_quadratic(α,       r, K-2:K, r[K]+dr, r[K]+2*dr)
+      bdry_βru_r, bdry_βru_rr = extrapolate_quadratic(βru,     r, K-2:K, r[K]+dr, r[K]+2*dr)
+      # only outer boundary is inflow bdry
+      bdry_ρ_r,   bdry_ρ_rr   = extrapolate_quadratic(init_ρ,  r, K-2:K, r[K]+dr, r[K]+2*dr)
+      bdry_vr_r,  bdry_vr_rr  = extrapolate_quadratic(init_vr, r, K-2:K, r[K]+dr, r[K]+2*dr)
+      bdry_ϵ_r,   bdry_ϵ_rr   = extrapolate_quadratic(init_ϵ,  r, K-2:K, r[K]+dr, r[K]+2*dr)
+    else
+      TODO("unknown boundary conditions for initial data $id: $bc")
+    end
   elseif id == "tov"
-    bdry_γrr_l, bdry_γrr_ll = γrr[1],  γrr[2]
-    bdry_α_l,   bdry_α_ll   = α[1],    α[2]
-    bdry_βru_l, bdry_βru_ll = -βru[1], -βru[2]
-    bdry_ρ_l,   bdry_ρ_ll   = ρ[1],    ρ[2]
-    bdry_vr_l,  bdry_vr_ll  = -vr[1],  -vr[2]
-    bdry_ϵ_l,   bdry_ϵ_ll   = ϵ[1],    ϵ[2]
+    if bc == "tov_reflective_origin"
+      bdry_γrr_l, bdry_γrr_ll = γrr[1],   γrr[2]
+      bdry_α_l,   bdry_α_ll   = α[1],     α[2]
+      bdry_βru_l, bdry_βru_ll = -βru[1],  -βru[2]
+      bdry_ρ_l,   bdry_ρ_ll   = ρ[1],     ρ[2]
+      bdry_vr_l,  bdry_vr_ll  = -vr[1],   -vr[2]
+      bdry_ϵ_l,   bdry_ϵ_ll   = ϵ[1],     ϵ[2]
+      bdry_γrr_r, bdry_γrr_rr = γrr[end], γrr[end]
+      bdry_α_r,   bdry_α_rr   = α[end],   α[end]
+      bdry_βru_r, bdry_βru_rr = βru[end], βru[end]
+      bdry_ρ_r,   bdry_ρ_rr   = ρ[end],   ρ[end]
+      bdry_vr_r,  bdry_vr_rr  = vr[end],  vr[end]
+      bdry_ϵ_r,   bdry_ϵ_rr   = ϵ[end],   ϵ[end]
+    elseif bc == "tov_symmetric_domain"
+      bdry_γrr_l, bdry_γrr_ll = γrr[1],   γrr[end]
+      bdry_α_l,   bdry_α_ll   = α[1],     α[end]
+      bdry_βru_l, bdry_βru_ll = βru[1],   βru[end]
+      bdry_ρ_l,   bdry_ρ_ll   = ρ[1],     ρ[end]
+      bdry_vr_l,  bdry_vr_ll  = vr[1],    vr[end]
+      bdry_ϵ_l,   bdry_ϵ_ll   = ϵ[1],     ϵ[end]
+      bdry_γrr_r, bdry_γrr_rr = γrr[end], γrr[end]
+      bdry_α_r,   bdry_α_rr   = α[end],   α[end]
+      bdry_βru_r, bdry_βru_rr = βru[end], βru[end]
+      bdry_ρ_r,   bdry_ρ_rr   = ρ[end],   ρ[end]
+      bdry_vr_r,  bdry_vr_rr  = vr[end],  vr[end]
+      bdry_ϵ_r,   bdry_ϵ_rr   = ϵ[end],   ϵ[end]
+    else
+      TODO("unknown boundary conditions for initial data $id: $bc")
+    end
   else
-    TODO("boundary conditions for initial data $id not implemented")
-  end
-  bγrr      = BdryBufferedData(bdry_γrr_ll, bdry_γrr_l, γrr, γrr[end], γrr[end-1])
-  bα        = BdryBufferedData(bdry_α_ll,   bdry_α_l,   α,   α[end],   α[end-1])
-  bβru      = BdryBufferedData(bdry_βru_ll, bdry_βru_l, βru, βru[end], βru[end-1])
-  bρ        = BdryBufferedData(bdry_ρ_ll,   bdry_ρ_l,   ρ,   ρ[end],   ρ[end-1])
-  bvr       = BdryBufferedData(bdry_vr_ll,  bdry_vr_l,  vr,  vr[end],  vr[end-1])
-  bϵ        = BdryBufferedData(bdry_ϵ_ll,   bdry_ϵ_l,   ϵ,   ϵ[end],   ϵ[end-1])
+  end
+  bγrr      = BdryBufferedData(bdry_γrr_ll, bdry_γrr_l, γrr, bdry_γrr_r, bdry_γrr_rr)
+  bα        = BdryBufferedData(bdry_α_ll,   bdry_α_l,   α,   bdry_α_r,   bdry_α_rr)
+  bβru      = BdryBufferedData(bdry_βru_ll, bdry_βru_l, βru, bdry_βru_r, bdry_βru_rr)
+  bρ        = BdryBufferedData(bdry_ρ_ll,   bdry_ρ_l,   ρ,   bdry_ρ_r,   bdry_ρ_rr)
+  bvr       = BdryBufferedData(bdry_vr_ll,  bdry_vr_l,  vr,  bdry_vr_r,  bdry_vr_rr)
+  bϵ        = BdryBufferedData(bdry_ϵ_ll,   bdry_ϵ_l,   ϵ,   bdry_ϵ_r,   bdry_ϵ_rr)
   bρL       = BdryBufferedData(0, 0, ρL,  0, 0)
   bvrL      = BdryBufferedData(0, 0, vrL, 0, 0)
   bϵL       = BdryBufferedData(0, 0, ϵL,  0, 0)
@@ -415,7 +449,7 @@ function step_fv_slope_limiter!(mesh, P::Project{:spherical1d}, state)
   bdϵ_m     = BdryBufferedData(0, 0, dϵ_m, 0, 0)
   bdϵ_p     = BdryBufferedData(0, 0, dϵ_p, 0, 0)
   @unpack max_v = get_static_variables(mesh)
-  bde_max_v = BdryBufferedData(max_v[1], max_v[1], max_v, max_v[end], max_v[end-1])
+  bde_max_v = BdryBufferedData(max_v[1], max_v[1], max_v, max_v[end], max_v[end])
 
   for vars in ( (bρL,bρR,bρ,bdρ_m,bdρ_p),
                 (bvrL,bvrR,bvr,bdvr_m,bdvr_p),
@@ -458,9 +492,11 @@ function step_fv_slope_limiter!(mesh, P::Project{:spherical1d}, state)
   end
 
   # for the outer bdry conditions for the TOV star we cheat here
-  rhs_D[end] = 0
-  rhs_Sr[end] = 0
-  rhs_τ[end] = 0
+  if bc == "tov_reflective_origin"
+    rhs_D[end] = rhs_Sr[end] = rhs_τ[end] = 0
+  elseif bc == "tov_symmetric_domain"
+    rhs_D[1] = rhs_Sr[1] = rhs_τ[1] = rhs_D[end] = rhs_Sr[end] = rhs_τ[end] = 0
+  end
 
   if !P.prmsdb["GRHD"]["atm_evolve"]
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
@@ -620,13 +656,47 @@ function step_fv_slope_limiter!(mesh, P::Project{:doublecartoon}, state)
   L, = dg1d.widths(mesh)
   dl = L/K
 
-  bγxx      = BdryBufferedData(γxx[2],      γxx[1],      γxx,      γxx[end],      γxx[end-1])
-  bα        = BdryBufferedData(α[2],        α[1],        α,        α[end],        α[end-1])
-  bβux      = BdryBufferedData(-βux[2],     -βux[1],     βux,      βux[end],      βux[end-1])
-  bsqrtdetγ = BdryBufferedData(sqrtdetγ[2], sqrtdetγ[1], sqrtdetγ, sqrtdetγ[end], sqrtdetγ[end-1])
-  bρ        = BdryBufferedData(ρ[2],   ρ[1],   ρ,  ρ[end],  ρ[end-1])
-  bvx       = BdryBufferedData(-vx[2], -vx[1], vx, vx[end], vx[end-1])
-  bϵ        = BdryBufferedData(ϵ[2],   ϵ[1],   ϵ,  ϵ[end],  ϵ[end-1])
+  bc = P.prmsdb["GRHD"]["bc"]
+  if bc == "tov_symmetric_domain"
+    bdry_γxx_l, bdry_γxx_ll           = γxx[1],        γxx[2]
+    bdry_α_l, bdry_α_ll               = α[1],          α[2]
+    bdry_βux_l, bdry_βux_ll           = βux[1],        βux[2]
+    bdry_sqrtdetγ_l, bdry_sqrtdetγ_ll = sqrtdetγ[1],   sqrtdetγ[2]
+    bdry_ρ_l, bdry_ρ_ll               = ρ[1],          ρ[2]
+    bdry_vx_l, bdry_vx_ll             = vx[1],         vx[2]
+    bdry_ϵ_l, bdry_ϵ_ll               = ϵ[1],          ϵ[2]
+    bdry_γxx_r,      bdry_γxx_rr      = γxx[end],      γxx[end]
+    bdry_α_r,        bdry_α_rr        = α[end],        α[end]
+    bdry_βux_r,      bdry_βux_rr      = βux[end],      βux[end]
+    bdry_sqrtdetγ_r, bdry_sqrtdetγ_rr = sqrtdetγ[end], sqrtdetγ[end]
+    bdry_ρ_r,        bdry_ρ_rr        = ρ[end],        ρ[end]
+    bdry_vx_r,       bdry_vx_rr       = vx[end],       vx[end]
+    bdry_ϵ_r,        bdry_ϵ_rr        = ϵ[end],        ϵ[end]
+  elseif bc == "tov_reflective_origin"
+    bdry_γxx_l, bdry_γxx_ll           = γxx[1],        γxx[2]
+    bdry_α_l, bdry_α_ll               = α[1],          α[2]
+    bdry_βux_l, bdry_βux_ll           = -βux[1],       -βux[2]
+    bdry_sqrtdetγ_l, bdry_sqrtdetγ_ll = sqrtdetγ[1],   sqrtdetγ[2]
+    bdry_ρ_l, bdry_ρ_ll               = ρ[1],          ρ[2]
+    bdry_vx_l, bdry_vx_ll             = -vx[1],        -vx[2]
+    bdry_ϵ_l, bdry_ϵ_ll               = ϵ[1],          ϵ[2]
+    bdry_γxx_r,      bdry_γxx_rr      = γxx[end],      γxx[end]
+    bdry_α_r,        bdry_α_rr        = α[end],        α[end]
+    bdry_βux_r,      bdry_βux_rr      = βux[end],      βux[end]
+    bdry_sqrtdetγ_r, bdry_sqrtdetγ_rr = sqrtdetγ[end], sqrtdetγ[end]
+    bdry_ρ_r,        bdry_ρ_rr        = ρ[end],        ρ[end]
+    bdry_vx_r,       bdry_vx_rr       = vx[end],       vx[end]
+    bdry_ϵ_r,        bdry_ϵ_rr        = ϵ[end],        ϵ[end]
+  else
+    TODO("implement boundary condition $bc")
+  end
+  bγxx      = BdryBufferedData(bdry_γxx_ll,      bdry_γxx_l,      γxx,      bdry_γxx_r,      bdry_γxx_rr)
+  bα        = BdryBufferedData(bdry_α_ll,        bdry_α_l,        α,        bdry_α_r,        bdry_α_rr)
+  bβux      = BdryBufferedData(bdry_βux_ll,      bdry_βux_l,      βux,      bdry_βux_r,      bdry_βux_rr)
+  bsqrtdetγ = BdryBufferedData(bdry_sqrtdetγ_ll, bdry_sqrtdetγ_l, sqrtdetγ, bdry_sqrtdetγ_r, bdry_sqrtdetγ_rr)
+  bρ        = BdryBufferedData(bdry_ρ_ll,        bdry_ρ_l,        ρ,        bdry_ρ_r,        bdry_ρ_rr)
+  bvx       = BdryBufferedData(bdry_vx_ll,       bdry_vx_l,       vx,       bdry_vx_r,       bdry_vx_rr)
+  bϵ        = BdryBufferedData(bdry_ϵ_ll,        bdry_ϵ_l,        ϵ,        bdry_ϵ_r,        bdry_ϵ_rr)
   bρL       = BdryBufferedData(0, 0, ρL,  0, 0)
   bvxL      = BdryBufferedData(0, 0, vxL, 0, 0)
   bϵL       = BdryBufferedData(0, 0, ϵL,  0, 0)
@@ -687,9 +757,11 @@ function step_fv_slope_limiter!(mesh, P::Project{:doublecartoon}, state)
   end
 
   # for the outer bdry conditions for the TOV star we cheat here
-  rhs_rD[end]  = 0
-  rhs_rSx[end] = 0
-  rhs_rτ[end]  = 0
+  if bc == "tov_reflective_origin"
+    rhs_rD[end] = rhs_rSx[end] = rhs_rτ[end] = 0
+  elseif bc == "tov_symmetric_domain"
+    rhs_rD[1] = rhs_rSx[1] = rhs_rτ[1] = rhs_rD[end] = rhs_rSx[end] = rhs_rτ[end] = 0
+  end
 
   if !P.prmsdb["GRHD"]["atm_evolve"]
     broadcast_volume!(determine_atmosphere, P.equation, mesh)
@@ -754,6 +826,8 @@ end
   @unpack bγxx, bα, bβux, bsqrtdetγ, bρ, bvx, bϵ, bρL, bvxL, bϵL, bpL,
           bρR, bvxR, bϵR, bpR, bdρ_m, bdρ_p, bdvx_m, bdvx_p, bdϵ_m, bdϵ_p, bde_max_v,
           nflx_rD_L, nflx_rSx_L, nflx_rτ_L, nflx_rD_R, nflx_rSx_R, nflx_rτ_R, eos = vars
+  R          = MMatrix{5,5,Float64,25}(0 for _ in 1:25)
+  Δu, Δω, up_u = [ MVector{5,Float64}(0 for _ in 1:5) for _ in 1:4 ]
   for k in 1:K
     # numerical fluxes
     # subscript + ... cell interior value
@@ -829,22 +903,21 @@ end
       # also see Font 2000, PRD 61 044011, eq 64, but only contains right eigenvectors
       #
       # here we assume vx=vz=0 and an isotropic metric
-      R = @SMatrix [ KK/(h₀*W₀)      0.0        0.0       1.0                                            1.0;
+      R .= @SMatrix [ KK/(h₀*W₀)      0.0        0.0       1.0                                            1.0;
                      vx₀             0.0        0.0       h₀*W₀*(vx₀-(vux₀-Λp)/(γuuxx₀-vux₀*Λp))         h₀*W₀*(vx₀-(vux₀-Λm)/(γuuxx₀-vux₀*Λm));
                      0.0             h₀*γxx₀    0.0       0.0                                            0.0;
                      0.0             0.0        h₀*γxx₀   0.0                                            0.0;
                      1.0-KK/(h₀*W₀)  0.0        0.0       h₀*W₀*(γuuxx₀-vux₀*vux₀)/(γuuxx₀-vux₀*Λp)-1.0  h₀*W₀*(γuuxx₀-vux₀*vux₀)/(γuuxx₀-vux₀*Λm)-1.0
                    ]
-      u₊ = @SVector [ D₊, Sx₊, 0.0, 0.0, τ₊ ]
-      u₋ = @SVector [ D₋, Sx₋, 0.0, 0.0, τ₋ ]
-      Δω = R \ (u₋-u₊)
+      Δu .= @SVector [ D₋-D₊, Sx₋-Sx₊, 0.0, 0.0, τ₋-τ₊ ]
+      Δω .= R \ Δu
 
       # roe upwinding term
       # cf. Font 2000, PRD 61 044011, eq 64
       # cf. Hesthaven 2018, Numerical methods for conservation laws,
       #   sec. 3.2 - Riemann solver
       #   sec. 6.2.1 - Roe flux
-      up_u = @. abs(λ0)*Δω[1]*R[:,1] + abs(λ0)*Δω[2]*R[:,2] + abs(λ0)*Δω[3]*R[:,3] + abs(λp)*Δω[4]*R[:,4] + abs(λm)*Δω[5]*R[:,5]
+      @views up_u .= @. abs(λ0)*Δω[1]*R[:,1] + abs(λ0)*Δω[2]*R[:,2] + abs(λ0)*Δω[3]*R[:,3] + abs(λp)*Δω[4]*R[:,4] + abs(λm)*Δω[5]*R[:,5]
       up_D, up_Sx, up_τ = up_u[1], up_u[2], up_u[5]
 
       f₊, f₋       = sdγ₀*α₀*D₊*(vxu₊-βux₀/α₀), sdγ₀*α₀*D₋*(vxu₋-βux₀/α₀)

src/GRHD/setup.jl

@@ -138,6 +138,7 @@ function register_evolution!(mesh, P::Project{:valenciad1d})
                              :ρ, :vr, :p, :ϵ,          # primitives
                              :max_v, :r,               # maximum charactersitic speed, radius
                              :init_D, :init_Sr, :init_τ, :init_max_v, :init_vr, :init_p,
+                             :init_ρ, :init_ϵ,
                              :E, :Em1, :flr_E, :flr_Em1, :EP,
                              :ρhW2),
     bdry_variablenames    = (:bdry_D, :bdry_Sr, :bdry_τ,
@@ -160,6 +161,7 @@ function register_evolution!(mesh, P::Project{:spherical1d})
                              :ρ, :vr, :p, :ϵ,          # primitives
                              :max_v, :r,               # maximum charactersitic speed, radius
                              :init_D, :init_Sr, :init_τ, :init_max_v, :init_vr, :init_p,
+                             :init_ρ, :init_ϵ,
                              :E, :Em1, :flr_E, :flr_Em1, :EP,
                              :ρhW2),
     bdry_variablenames    = (:bdry_D, :bdry_Sr, :bdry_τ,
@@ -183,6 +185,7 @@ function register_evolution!(mesh::Mesh1d{FVElement}, P::Project{:spherical1d})
                                :ρ, :vr, :p, :ϵ,          # primitives
                                :max_v, :r,               # maximum charactersitic speed, radius
                                :init_D, :init_Sr, :init_τ, :init_max_v, :init_vr, :init_p,
+                               :init_ρ, :init_ϵ,
                                :E, :Em1, :flr_E, :flr_Em1, :EP,
                                :ρhW2),
       bdry_variablenames    = (:bdry_D, :bdry_Sr, :bdry_τ,
@@ -214,6 +217,7 @@ function register_evolution!(mesh::Mesh1d{FVElement}, P::Project{:spherical1d})
                                :max_v, :r,               # maximum charactersitic speed, radius
                                :E, :Em1, :flr_E, :flr_Em1, :EP,
                                :init_D, :init_Sr, :init_τ, :init_max_v, :init_vr, :init_p,
+                               :init_ρ, :init_ϵ,
                                :ρhW2),
       global_variablenames  = (:t, :tm1) # time steps
     )
@@ -235,6 +239,7 @@ function register_evolution!(mesh, P::Project{:rescaled_spherical1d})
                              :ρ, :vr, :p, :ϵ,             # primitives
                              :max_v, :r,                  # maximum charactersitic speed, radius
                              :init_D, :init_Sr, :init_τ, :init_max_v, :init_vr, :init_p,
+                             :init_ρ, :init_ϵ,
                              :E, :Em1, :flr_E, :flr_Em1, :EP,
                              :ρhW2),
     bdry_variablenames    = (:bdry_D, :bdry_Sr, :bdry_τ,
@@ -291,6 +296,7 @@ function register_evolution!(mesh::Mesh1d{FVElement}, P::Project{:doublecartoon}
                                :D, :Sx, :τ,
                                :max_v, :r,                # maximum charactersitic speed, radius
                                :init_D, :init_Sx, :init_τ, :init_max_v, :init_vx, :init_p,
+                               :init_ρ, :init_ϵ,
                                :buf_cartoon, :ρhW2),
       bdry_variablenames    = (:bdry_rD, :bdry_rSx, :bdry_rτ,
                                :bdry_ρ, :bdry_ϵ, :bdry_vx, :bdry_p,

test/IntegrationTests/refs/fv_minmod_marquina_grhd_tov_doublecartoon/fv_minmod_marquina_grhd_tov_doublecartoon.toml

@@ -4,7 +4,7 @@ polytrope_k = 100.0
 polytrope_gamma = 2.0
 
 [GRHD]
-bc = "from_id"
+bc = "tov_reflective_origin"
 id = "tov"
 formulation = "doublecartoon"
 id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"

test/IntegrationTests/refs/fv_minmod_roe_grhd_tov_doublecartoon/fv_minmod_roe_grhd_tov_doublecartoon.toml

@@ -4,7 +4,7 @@ polytrope_k = 100.0
 polytrope_gamma = 2.0
 
 [GRHD]
-bc = "from_id"
+bc = "tov_reflective_origin"
 id = "tov"
 formulation = "doublecartoon"
 id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"