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

@@ -2,49 +2,51 @@
 #                              Timestep                               #
 #######################################################################
 
+function compute_maxspeed(::Project{:spherical1d}, equation, mesh)
+  broadcast_volume!(maxspeed, equation, mesh)
+end
+function compute_maxspeed(::Project{:rescaled_spherical1d}, equation, mesh)
+  broadcast_volume!(maxspeed_rescaled_spherical1d, equation, mesh)
+end
 
-function timestep(env, P::Project, hrsc::Maybe{HRSC.AbstractHRSC})
-  @unpack cache, mesh = env
-  @unpack equation = P
-  @unpack max_v = get_static_variables(cache)
-
-  broadcast_volume!(cons2prim, equation, cache)
-  broadcast_volume!(speed, equation, cache)
-  vmax = dg1d.absolute_maximum(max_v)
-  vmax_limit = 0.9999
+function get_maxspeed(mesh, equation)
+  @unpack max_v = get_static_variables(mesh.cache)
+  vmax = maximum(abs, max_v)
+  vmax_limit = 0.99999999999
   if vmax > vmax_limit
     @warn "Limiting timestep due to maximum speed exceeding $vmax_limit"
     vmax = vmax_limit
   end
+  return vmax
+end
 
-  @unpack dl, x, element = mesh
-  @unpack N = element
-
-  dt = dl / (N^2 * vmax)
+function timestep(mesh, P::Project, hrsc::Maybe{HRSC.AbstractHRSC})
+  compute_maxspeed(P, P.equation, mesh)
+  vmax = get_maxspeed(mesh, P.equation)
+  dl = min_grid_spacing(mesh)
+  dt = dl / (vmax * dtfactor(mesh))
   return dt
 end
+dtfactor(mesh::Mesh1d{FVElement}) = 2
+dtfactor(mesh::Mesh1d{SpectralElement}) = mesh.element.N
 
 
-function timestep(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity)
-  @unpack cache, mesh = env
+function timestep(mesh, P::Project, hrsc::HRSC.AbstractArtificialViscosity)
+  @unpack cache = mesh
   @unpack equation = P
 
-  broadcast_volume!(cons2prim, equation, cache)
-  broadcast_volume!(speed, equation, cache)
-  max_v = get_variable(cache, :max_v)
-  vmax = dg1d.absolute_maximum(view(max_v,:))
-  vmax_limit = 0.9999
-  if vmax > vmax_limit
-    @warn "Limiting timestep due to maximum speed exceeding $vmax_limit"
-    vmax = vmax_limit
-  end
-  smoothed_mu = get_variable(cache, :smoothed_mu)
-  mumax = dg1d.absolute_maximum(view(smoothed_mu, :))
+  compute_maxspeed(P, P.equation, mesh)
+  vmax = get_maxspeed(mesh, P.equation)
+  @unpack mu = get_static_variables
+  mumax = maximum(abs, mu)
 
-  @unpack dl, x, element = mesh
+  @unpack element = mesh
   @unpack N = element
+  L, = dg1d.widths(mesh)
+  K, = mesh.tree.dims
+  dL = L/K
 
-  dt = 1 / (vmax * N / dl + mumax * N^3 / dl^2)
+  dt = 1 / (vmax * N / dL + mumax * N^3 / dL^2)
   dt_limit = 1e-7
   if dt < dt_limit
     @warn "Limiting timestep due to dt being smaller than $dt_limit"
@@ -59,38 +61,140 @@ end
 #######################################################################
 
 
-function rhs!(env, P::Project, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryconds)
+function rhs!(mesh::Mesh1d{FVElement}, P::Project{:spherical1d}, hrsc::Nothing)
 
-  @unpack cache, mesh = env
-  @unpack equation, rsolver = P
+  @unpack cache = mesh
+  @unpack equation = P
 
-  @unpack sD, sSr, stau             = get_dynamic_variables(cache)
-  @unpack flx_sD, flx_sSr, flx_stau = get_static_variables(cache)
-  @unpack rhs_sD, rhs_sSr, rhs_stau = get_rhs_variables(cache)
-  @unpack lhs_numflx_sD,   rhs_numflx_sD,
-          lhs_numflx_sSr,  rhs_numflx_sSr,
-          lhs_numflx_stau, rhs_numflx_stau = get_bdry_variables(cache)
+  @unpack D, Sr, τ                = get_dynamic_variables(cache)
+  @unpack flr_D, flr_Sr, flr_τ,
+          max_v, ρ, ϵ, vr, p,
+          src_D, src_Sr, src_τ    = get_static_variables(cache)
+  @unpack rhs_D, rhs_Sr, rhs_τ    = get_rhs_variables(cache)
+  @unpack nflr_D, nflr_Sr, nflr_τ,
+          bdry_D, bdry_Sr, bdry_τ = get_bdry_variables(cache)
 
-  broadcast_volume!(scaledcons2cons, equation, cache)
-  postprocess_scaledcons2cons!(cache, equation, mesh)
-  broadcast_volume!(cons2prim, equation, cache)
-  broadcast_volume!(flux, equation, cache)
-  broadcast_volume!(sources, equation, cache)
-  broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
-  broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
+  broadcast_volume!(cons2prim, equation, mesh)
+  broadcast_volume!(maxspeed, equation, mesh)
+  broadcast_volume!(flux_source_spherical1d, equation, mesh)
+  broadcast_bdry!(fv_bdry_flux, equation, mesh)
 
-  @unpack p, vr, rho, eps = get_static_variables(cache)
+  dt = timestep(mesh, P, hrsc)
+  fv_update_step!(rhs_D,  D,  flr_D,  src_D,  bdry_D,  nflr_D,  dt, mesh)
+  fv_update_step!(rhs_Sr, Sr, flr_Sr, src_Sr, bdry_Sr, nflr_Sr, dt, mesh)
+  fv_update_step!(rhs_τ,  τ,  flr_τ,  src_τ,  bdry_τ,  nflr_τ,  dt, mesh)
+
+  return
+end
+
+
+function rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
+
+  @unpack cache = mesh
+  @unpack equation = P
+
+  @unpack D, Sr, τ                    = get_dynamic_variables(cache)
+  @unpack rhs_D, rhs_Sr, rhs_τ        = get_rhs_variables(cache)
+  @unpack flr_D, flr_Sr, flr_τ,
+          max_v, vr, p,
+          src_D, src_Sr, src_τ        = get_static_variables(cache)
+  @unpack nflr_D, nflr_Sr, nflr_τ,
+          bdry_D, bdry_Sr, bdry_τ,
+          bdry_max_v, bdry_vr, bdry_p = get_bdry_variables(cache)
+
+  broadcast_volume!(cons2prim, equation, mesh)
+  broadcast_volume!(maxspeed, equation, mesh)
+
+  dg1d.interpolate_face_data!(mesh, D,     bdry_D)
+  dg1d.interpolate_face_data!(mesh, Sr,    bdry_Sr)
+  dg1d.interpolate_face_data!(mesh, τ,     bdry_τ)
+  dg1d.interpolate_face_data!(mesh, max_v, bdry_max_v)
+  dg1d.interpolate_face_data!(mesh, vr,    bdry_vr)
+  dg1d.interpolate_face_data!(mesh, p,     bdry_p)
+
+  broadcast_volume!(flux_source_spherical1d, equation, mesh)
+  impose_symmetry_sources!(P, mesh)
+  broadcast_faces!(llf, equation, mesh)
+  broadcast_bdry!(bdryllf, equation, mesh)
+
+  compute_rhs_weak_form!(rhs_D,  flr_D,  src_D,  nflr_D,  mesh)
+  compute_rhs_weak_form!(rhs_Sr, flr_Sr, src_Sr, nflr_Sr, mesh)
+  compute_rhs_weak_form!(rhs_τ,  flr_τ,  src_τ,  nflr_τ,  mesh)
+
+  return
+end
 
-  compute_rhs_weak_form!(rhs_sD,   flx_sD,    lhs_numflx_sD,   rhs_numflx_sD,   mesh)
-  compute_rhs_weak_form!(rhs_sSr,  flx_sSr,   lhs_numflx_sSr,  rhs_numflx_sSr,  mesh)
-  compute_rhs_weak_form!(rhs_stau, flx_stau,  lhs_numflx_stau, rhs_numflx_stau, mesh)
+
+function rhs!(mesh::Mesh1d{FVElement}, P::Project{:rescaled_spherical1d}, hrsc::Nothing)
+
+  @unpack cache = mesh
+  @unpack equation = P
+
+  @unpack rD, rSr, rτ                = get_dynamic_variables(cache)
+  @unpack flr_rD, flr_rSr, flr_rτ,
+          D, Sr, τ,
+          src_rD, src_rSr, src_rτ    = get_static_variables(cache)
+  @unpack rhs_rD, rhs_rSr, rhs_rτ    = get_rhs_variables(cache)
+  @unpack nflr_rD, nflr_rSr, nflr_rτ,
+          bdry_D, bdry_Sr, bdry_τ = get_bdry_variables(cache)
+
+  broadcast_volume!(unscale_conservatives, equation, mesh)
+  impose_symmetry!(P, mesh)
+  broadcast_volume!(cons2prim_rescaled_spherical1d, equation, mesh)
+  broadcast_volume!(maxspeed_rescaled_spherical1d, equation, mesh)
+  broadcast_volume!(flux_source_rescaled_spherical1d, equation, mesh)
+  broadcast_bdry!(fv_bdry_flux_rescaled_spherical1d, equation, mesh)
+
+  dt = timestep(mesh, P, hrsc)
+  fv_update_step!(rhs_rD,  rD,  flr_rD,  src_rD,  bdry_D,  nflr_rD,  dt, mesh)
+  fv_update_step!(rhs_rSr, rSr, flr_rSr, src_rSr, bdry_Sr, nflr_rSr, dt, mesh)
+  fv_update_step!(rhs_rτ,  rτ,  flr_rτ,  src_rτ,  bdry_τ,  nflr_rτ,  dt, mesh)
+
+  return
+end
+
+
+function rhs!(mesh::Mesh1d, P::Project{:rescaled_spherical1d}, hrsc::Nothing)
+
+  @unpack cache = mesh
+  @unpack equation = P
+
+  @unpack rD, rSr, rτ                    = get_dynamic_variables(cache)
+  @unpack rhs_rD, rhs_rSr, rhs_rτ        = get_rhs_variables(cache)
+  @unpack flr_rD, flr_rSr, flr_rτ,
+          max_v, vr, p, D, Sr, τ,
+          src_rD, src_rSr, src_rτ        = get_static_variables(cache)
+  @unpack nflr_rD, nflr_rSr, nflr_rτ,
+          bdry_D, bdry_Sr, bdry_τ,
+          bdry_max_v, bdry_vr, bdry_p = get_bdry_variables(cache)
+
+  broadcast_volume!(unscale_conservatives, equation, mesh)
+  impose_symmetry!(P, mesh)
+  broadcast_volume!(cons2prim_rescaled_spherical1d, equation, mesh)
+  broadcast_volume!(maxspeed_rescaled_spherical1d, equation, mesh)
+
+  dg1d.interpolate_face_data!(mesh, D,     bdry_D)
+  dg1d.interpolate_face_data!(mesh, Sr,    bdry_Sr)
+  dg1d.interpolate_face_data!(mesh, τ,     bdry_τ)
+  dg1d.interpolate_face_data!(mesh, max_v, bdry_max_v)
+  dg1d.interpolate_face_data!(mesh, vr,    bdry_vr)
+  dg1d.interpolate_face_data!(mesh, p,     bdry_p)
+
+  broadcast_volume!(flux_source_rescaled_spherical1d, equation, mesh)
+  broadcast_faces!(llf_rescaled_spherical1d, equation, mesh)
+  broadcast_bdry!(bdryllf_rescaled_spherical1d, equation, mesh)
+
+  compute_rhs_weak_form!(rhs_rD,  flr_rD,  src_rD,  nflr_rD,  mesh)
+  compute_rhs_weak_form!(rhs_rSr, flr_rSr, src_rSr, nflr_rSr, mesh)
+  compute_rhs_weak_form!(rhs_rτ,  flr_rτ,  src_rτ,  nflr_rτ,  mesh)
 
   return
 end
 
 
-function rhs!(env, P::Project, hrsc::HRSC.AbstractReconstruction,
-    bdryconds, ldg_bdryconds, av_bdryconds)
+function rhs!(env, P::Project, hrsc::HRSC.AbstractReconstruction)
+
+  TODO()
 
   @unpack cache, mesh = env
   @unpack equation, rsolver = P
@@ -122,8 +226,9 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractReconstruction,
 end
 
 
-function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
-    bdryconds, ldg_bdryconds, av_bdryconds)
+function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity)
+
+  TODO()
 
   @unpack cache, mesh = env
   @unpack equation, rsolver, ldg_rsolver, av_rsolver = P