examples/simple_wave/aventropy.toml

+[EquationOfState]
+eos = "polytrope"
+polytrope_gamma = 1.6666666666666666667
+polytrope_k = 100.0
+
+[SRHD]
+id = "simple_wave"
+bc = "from_id"
+
+[Mesh]
+range = [ -1.5, 1.5 ]
+n = 5
+k = 125
+basis = "lgl"
+periodic = true
+
+[Output]
+variables       = [ "D", "S", "tau", "p", "rho", "eps", "v", "EP", "smoothed_mu", "E", "flx_E" ]
+# variables       = [ "D", "S", "tau" ]
+aligned_ts      = "$(collect(range(0.01,6.4,step=0.01)))"
+enable1d        = true
+
+[Evolution]
+cfl = 0.125
+tend = 6.4
+# tend = 0.2
+
+[Log]
+progress_stdout = true
+
+[HRSC]
+method = "av"
+av_method = "entropy"
+entropy_cmax = 1.0
+entropy_ce   = 1.0
+
+[TCI]
+method = "entropy_production"

examples/simple_wave/avmda.toml

+[EquationOfState]
+eos = "polytrope"
+polytrope_gamma = 1.6666666666666666667
+polytrope_k = 100.0
+
+[SRHD]
+id = "simple_wave"
+bc = "from_id"
+
+[Mesh]
+range = [ -1.5, 1.5 ]
+n = 5
+k = 200
+basis = "lgl"
+periodic = true
+
+[Output]
+variables       = [ "D", "S", "tau", "p", "rho", "eps", "v", "EP", "smoothed_mu", "E", "flx_E" ]
+aligned_ts      = "$(collect(range(0.01,6.4,step=0.01)))"
+enable1d        = true
+
+[Evolution]
+cfl = 0.125
+tend = 6.4
+# tend = 0.2
+
+[Log]
+progress_stdout = true
+
+[HRSC]
+method = "av"
+av_method = "mda"
+
+[TCI]
+method = "entropy_production"

examples/simple_wave/dg.toml

+[EquationOfState]
+eos = "polytrope"
+polytrope_gamma = 1.6666666666666666667
+polytrope_k = 100.0
+
+[SRHD]
+id = "simple_wave"
+bc = "from_id"
+
+[Mesh]
+range = [ -1.5, 1.5 ]
+n = 5
+k = 200
+basis = "lgl"
+periodic = false
+
+[Output]
+# variables       = [ "D", "S", "tau", "p", "rho", "eps", "E", "Em1", "flx_E", "flx_Em1", "v", "EP", "smoothed_mu" ]
+variables       = [ "D", "S", "tau" ]
+aligned_ts      = "$(collect(range(0.01,0.6,step=0.01)))"
+enable1d        = true
+
+[Evolution]
+cfl = 0.2
+tend = 0.6
+# tend = 0.2
+
+[Log]
+progress_stdout = true
+
+# [HRSC]
+# method = "av"
+# av_method = "entropy"
+# entropy_cmax = 0.1
+# entropy_ce   = 1.0
+
+[TCI]
+method = "entropy_production"

examples/simple_wave/fv.toml

+[EquationOfState]
+eos = "polytrope"
+polytrope_gamma = 1.6666666666666666667
+polytrope_k = 100.0
+
+[SRHD]
+id = "simple_wave"
+bc = "from_id"
+
+[Mesh]
+range = [ -1.5, 1.5 ]
+k = 2056
+periodic = false
+scheme = "FV"
+
+[Output]
+variables       = [ "D", "S", "tau" ]
+aligned_ts      = "$(collect(range(0.01,0.6,step=0.01)))"
+enable1d        = true
+
+[Evolution]
+cfl = 1.0
+tend = 0.6

examples/srhd_sod_shocktube/avmda.toml

+[EquationOfState]
+eos = "idealgas"
+idealgas_gamma = 1.4
+
+[SRHD]
+id = "sod_shocktube"
+bc = "from_id"
+
+[Mesh]
+range = [ -1.0, 1.0 ]
+n = 4
+k = 205
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# variables       = [ "D", "S", "tau", "p", "rho", "eps", "v", "smoothed_mu",
+# variables       = [ "D", "S", "tau", "smoothed_mu",
+#                     "ldg_D", "ldg_S", "ldg_tau", "rhs_D", "rhs_S", "rhs_tau" ]
+variables       = [ "D", "S", "tau", "p", "rho", "eps", "v" ]
+aligned_ts      = "$(collect(range(0.01,0.6,step=0.01)))"
+enable1d        = true
+
+[Evolution]
+cfl = 0.4
+tend = 0.6
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 20.0
+
+[TCI]
+method = "entropy_production"

examples/srhd_sod_shocktube/fv.toml

+[EquationOfState]
+eos = "idealgas"
+idealgas_gamma = 1.4
+
+[SRHD]
+id = "sod_shocktube"
+bc = "from_id"
+
+[Mesh]
+range = [ -1.0, 1.0 ]
+k = 1024
+scheme = "FV"
+periodic = false
+
+[Output]
+# every_iteration = 1
+variables       = [ "D", "S", "tau", "p", "rho", "eps", "v" ]
+aligned_ts      = "$(collect(range(0.01,0.6,step=0.01)))"
+enable1d        = true
+
+[Evolution]
+cfl = 1.0
+tend = 0.6

src/SRHD/callbacks.jl

@@ -49,6 +49,7 @@ function callback_equation(state_t, isperiodic, eq, mesh)
   Em1     .= E
   flx_Em1 .= flx_E
   dg1d.new_broadcast_volume!(entropy_variables, eq, mesh)
+  dg1d.new_broadcast_volume!(maxspeed, eq, mesh)
   dt <= 0.0 && return
 
   LL      = layout(mesh)
@@ -58,7 +59,6 @@ function callback_equation(state_t, isperiodic, eq, mesh)
   Em1     = dg1d.vreshape(Em1, LL)
   flx_Em1 = dg1d.vreshape(flx_Em1, LL)
   invdetJ = dg1d.vreshape(invdetJ, LL)
-  v       = dg1d.vreshape(v, LL)
 
   # compute entropy residual
   @turbo @. flx_Em1 += flx_E
@@ -78,6 +78,7 @@ function callback_equation(state_t, isperiodic, eq::Equation, mesh::Mesh1d{FVEle
   Em1     .= E
   flx_Em1 .= flx_E
   dg1d.new_broadcast_volume!(entropy_variables, eq, mesh)
+  dg1d.new_broadcast_volume!(maxspeed, eq, mesh)
 end
 
 
@@ -103,12 +104,13 @@ function callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc::HRSC.Abstract
   @unpack D, S, tau             = get_dynamic_variables(cache)
   @unpack t, tm1                = get_global_variables(cache)
 
-  max_v = cellmax_v
-  t[1] <= 0.0 && tm1[1] <= 0 && return
-  Npts, K = layout(mesh)
-  mat_max_v = reshape(view(v, :), (Npts, K))
-  for k = 1:K
-    max_v[k] = maximum(view(mat_max_v, :, k))
+  # prevent AV recomputation when callbacks are forced on initial data
+  # need this for shock tube tests so that we have non-zero viscosity in troublesome places
+  # although we start with smoothed initial data
+  t[1] <= tm1[1] && return
+
+  for (k,vi) in enumerate(eachcell(mesh, max_v))
+    cellmax_v[k] = maximum(vi)
   end
 
   HRSC.compute_viscosity!(

src/SRHD/cons2prim.jl

@@ -109,7 +109,7 @@ end
   r = abs(ri)
   v0 = abs(v0i)
 
-  if v0 >= 1
+  if v0 > 1
     error((v0,z0i,Si,D))
   end
 
@@ -134,10 +134,31 @@ end
   if D < rhomin
     @warn "atmosphering I"
 
+    error((D,S,tau))
+
+    @assert eos isa EquationOfState.IdealGas
+    @unpack Gamma = eos
+
+    # employ atmosphere values by taking the zero temperature limit
+    # for an ideal gas eos, we have that eos.Gamma becomes the Gamma for
+    # a polytrope, see grhd notes for derivation
+
     rho = rhoatm
-    eps = eos.cold_eos(InternalEnergy, Density, rho)
     v = 0.0
-    p = eos.cold_eos(Pressure, Density, rho)
+    W = 1.0
+    # p = K * rho^Gamma
+    # eps = Gamma * rho^
+    # D = W * rho
+    # S = rho * (1+eps+p/rho) * W^2 * v
+    # tau = rho * (1+eps+p/rho) * W^2 - p - D
+
+    # if eos isa ColdEquationOfState
+    #   eps = eos.cold_eos(InternalEnergy, Density, rho)
+    #   p = eos.cold_eos(Pressure, Density, rho)
+    # else
+    #   eps = eos(InternalEnergy, Density, rho)
+    #   p = eos(Pressure, Density, rho)
+    # end
 
   else
 

src/SRHD/equation.jl

@@ -189,6 +189,15 @@ end
 end
 
 
+@with_signature [legacy=false] function av_flux_covariant(equations::Equation)
+  @accepts rhs_D, rhs_S, rhs_tau, flx_D, flx_S, flx_tau, ldg_D, ldg_S, ldg_tau, smoothed_mu
+  flx_D   += -smoothed_mu * (rhs_D - ldg_D)
+  flx_S   += -smoothed_mu * (rhs_S - ldg_S)
+  flx_tau += -smoothed_mu * (rhs_tau - ldg_tau)
+  @returns flx_D, flx_S, flx_tau
+end
+
+
 # @with_signature [legacy=false] function av_nflux(eq::Equation)
 #   @accepts        ldg_D, ldg_S, ldg_tau, smoothed_mu
 #   @accepts [bdry] bdry_ldg_D, bdry_ldg_S, bdry_ldg_tau, bdry_smoothed_mu
@@ -277,6 +286,45 @@ end
 end
 
 
+@with_signature [legacy=false] function av_nflux_covariant(eq::Equation)
+  @accepts        rhs_D, rhs_S, rhs_tau, ldg_D, ldg_S, ldg_tau, smoothed_mu
+  @accepts [bdry] bdry_rhs_D, bdry_rhs_S, bdry_rhs_tau, bdry_ldg_D, bdry_ldg_S, bdry_ldg_tau, bdry_smoothed_mu
+  @accepts [bdry] nflx_D, nflx_S, nflx_tau, nx
+
+  nflx       = -nx*(rhs_D - ldg_D)
+  bdry_nflx  = -nx*(bdry_rhs_D - bdry_ldg_D)
+  nflx_D    += 0.5 * (smoothed_mu*nflx + bdry_smoothed_mu*bdry_nflx)
+
+  nflx       = -nx*(rhs_S - ldg_S)
+  bdry_nflx  = -nx*(bdry_rhs_S - bdry_ldg_S)
+  nflx_S    += 0.5 * (smoothed_mu*nflx + bdry_smoothed_mu*bdry_nflx)
+
+  nflx       = -nx*(rhs_tau - ldg_tau)
+  bdry_nflx  = -nx*(bdry_rhs_tau - bdry_ldg_tau)
+  nflx_tau  += 0.5 * (smoothed_mu*nflx + bdry_smoothed_mu*bdry_nflx)
+
+  @returns [bdry] nflx_D, nflx_S, nflx_tau
+end
+@with_signature function bdry_av_nflux_covariant(eq::Equation)
+  @accepts Prefix(lhs), ldg_D, ldg_S, ldg_tau, smoothed_mu
+  @accepts Prefix(rhs), ldg_D, ldg_S, ldg_tau, smoothed_mu
+  @accepts nflx_D, nflx_S, nflx_tau, nx
+
+  TODO()
+
+  lhs_nflx = nx*lhs_ldg_D
+  nflx_D += lhs_smoothed_mu*lhs_nflx
+
+  lhs_nflx = nx*lhs_ldg_S
+  nflx_S += lhs_smoothed_mu*lhs_nflx
+
+  lhs_nflx = nx*lhs_ldg_tau
+  nflx_tau += lhs_smoothed_mu*lhs_nflx
+
+  @returns nflx_D, nflx_S, nflx_tau
+end
+
+
 @with_signature [legacy=false] function fv_bdry_flux(equation::Equation)
   @accepts D, S, tau, v, p
   @accepts init_D, init_S, init_tau, init_v, init_p

src/SRHD/initialdata.jl

@@ -6,6 +6,9 @@ function initialdata!(env, P::Project, prms, mesh::Mesh1d)
   initialdata_hrsc!(env, P, P.hrsc, prms, mesh)
   initialdata_tci!(env, P, P.tci, prms, mesh)
   initialdata_bdrycond!(env, P, dg1d.DirichletBC2(), mesh)
+  if get(P.prms, :id_smooth, false)
+    initialdata_smooth!(env, P, mesh)
+  end
   return
 end
 
@@ -163,19 +166,16 @@ end
 function initialdata_equation_of_state(
     ::Val{:sod_shocktube}, eos::EquationOfState.IdealGas, mesh)
 
-  TODO()
-
+  # from EFL paper: arxiv 2202.08839
   @unpack x = get_static_variables(mesh.cache)
   (xmin, xmax), = mesh.extends
 
   @toggled_assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
-  @toggled_assert eos.gamma == 5/3
+  @toggled_assert eos.gamma ≈ 1.4
   x0 = 0.0
-  # ρ0 = [ (xi < x0) ? 1.0  :  0.125  for xi in x ]
-  # v0 = [ (xi < x0) ? 0.0  :  0.0    for xi in x ]
-  # p0 = [ (xi < x0) ? 1.0  :  0.1    for xi in x ]
-
-  xl, xr = -0.5, 0.5
+  ρ0 = [ (xi < x0) ? 1.0  :  0.125  for xi in x ]
+  v0 = [ (xi < x0) ? 0.0  :  0.0    for xi in x ]
+  p0 = [ (xi < x0) ? 1.0  :  0.1    for xi in x ]
 
   (ρ0, v0, p0)
 end
@@ -225,7 +225,7 @@ function initialdata_equation_of_state(
   @unpack x = get_static_variables(mesh.cache)
   (xmin, xmax), = mesh.extends
 
-  @toggled_assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
+  @toggled_assert isapprox(xmin, -0.5) && isapprox(xmax, 0.5)
   @toggled_assert eos.gamma == 5/3
   x0 = 0.0
   ρ0 = [ (xi < x0) ? 10.0  : 1.0 for xi in x ]
@@ -516,3 +516,24 @@ function initialdata_bdrycond!(env, P::Project2d, bdrycond::dg1d.DirichletBC2, m
   init_flx_tau_y .= flx_tau_y
   init_max_v     .= max_v
 end
+
+
+#######################################################################
+#                       Initial data smoothing                        #
+#######################################################################
+
+
+function initialdata_smooth!(env, P::Project, mesh)
+  @unpack D, S, tau = get_dynamic_variables(mesh.cache)
+  @unpack mu = get_cell_variables(mesh.cache)
+  bernstein = BernsteinReconstruction(mesh)
+  flag = zeros(mesh.K)
+  for k = 1:mesh.K
+    flag[k] = Float64(mu[k] > 0.0)
+  end
+  HRSC.reconstruct!(D,    flag, bernstein, isperiodic=true)
+  HRSC.reconstruct!(S,    flag, bernstein, isperiodic=true)
+  HRSC.reconstruct!(tau,  flag, bernstein, isperiodic=true)
+  dg1d.new_broadcast_volume!(cons2prim_kastaun, P.equation, mesh)
+  dg1d.new_broadcast_volume!(maxspeed, P.equation, mesh)
+end

src/SRHD/rhs.jl

@@ -6,9 +6,9 @@ timestep(env, P, hrsc) = timestep(env, P, hrsc, env.mesh)
 
 
 function compute_maxspeed(mesh, equation, cache)
-  @unpack v = get_static_variables(cache)
+  @unpack max_v = get_static_variables(cache)
   dg1d.new_broadcast_volume!(maxspeed, equation, mesh)
-  vmax = dg1d.absolute_maximum(v)
+  vmax = dg1d.absolute_maximum(max_v)
   vmax_limit = 1e4
   if vmax > vmax_limit
     @warn "Limiting timestep due to maximum speed exceeding $vmax_limit"
@@ -206,11 +206,53 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractReconstruction,
 end
 
 
+function fd_diff_1d!(_df, _f, _x, mesh)
+  @unpack Npts = mesh.element
+  for (df, f, x) in zip(eachcell(mesh,_df),eachcell(mesh,_f),eachcell(mesh,_x))
+    for n in 2:Npts-1
+      fl,fm,fr = f[n-1],f[n],f[n+1]
+      xl,xm,xr = x[n-1],x[n],x[n+1]
+      dll = (xm-xr)/((xl-xm)*(xl-xr))
+      dlm = (2*xm-xl-xr)/((xm-xl)*(xm-xr))
+      dlr = (xm-xl)/((xr-xl)*(xr-xm))
+      df[n] = fl*dll + fm*dlm + fr*dlr
+    end
+    fl,fm,fr = f[1],f[2],f[3]
+    xl,xm,xr = x[1],x[2],x[3]
+    dll = (2*xl-xr-xm)/((xl-xm)*(xl-xr))
+    dlm = (xl-xr)/((xm-xl)*(xm-xr))
+    dlr = (xl-xm)/((xr-xl)*(xr-xm))
+    df[1] = fl*dll + fm*dlm + fr*dlr
+    fl,fm,fr = f[Npts-2],f[Npts-1],f[Npts]
+    xl,xm,xr = x[Npts-2],x[Npts-1],x[Npts]
+    dll = (xr-xm)/((xl-xm)*(xl-xr))
+    dlm = (xr-xl)/((xm-xl)*(xm-xr))
+    dlr = (2*xr-xl-xm)/((xr-xl)*(xr-xm))
+    df[Npts] = fl*dll + fm*dlm + fr*dlr
+  end
+end
+
+
 function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
     bdryconds, ldg_bdryconds, av_bdryconds, ::Mesh1d)
 
-  @unpack cache, mesh           = env
-  @unpack equation              = P
+  @unpack cache, mesh = env
+  @unpack equation, prms = P
+  reg = prms.av_regularization
+
+  if reg === :mono
+    rhs_mono!(cache, mesh, equation, P)
+  elseif reg === :covariant
+    rhs_covariant!(cache, mesh, equation, P)
+  else
+    TODO(reg)
+  end
+
+end
+
+
+function rhs_mono!(cache, mesh, equation, P)
+
   @unpack D, S, tau             = get_dynamic_variables(cache)
   @unpack flx_D, flx_S, flx_tau,
           ldg_D, ldg_S, ldg_tau,
@@ -224,6 +266,7 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
           bdry_ldg_D, bdry_ldg_S,
           bdry_ldg_tau,
           bdry_smoothed_mu      = get_bdry_variables(cache)
+  @unpack x = get_static_variables(cache)
 
   ## solve auxiliary equation: q + ∂x u = 0
   dg1d.new_broadcast_faces!(ldg_nflux, equation, mesh)
@@ -232,7 +275,12 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
   compute_rhs_weak_form!(ldg_D,   D,   nflx_D,   mesh)
   compute_rhs_weak_form!(ldg_S,   S,   nflx_S,   mesh)
   compute_rhs_weak_form!(ldg_tau, tau, nflx_tau, mesh)
+  # fd_diff_1d!(ldg_D,   D,   x, mesh)
+  # fd_diff_1d!(ldg_S,   S,   x, mesh)
+  # fd_diff_1d!(ldg_tau, tau, x, mesh)
+
 
+  # dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
   dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
   dg1d.new_broadcast_volume!(maxspeed, equation, mesh)
 
@@ -270,6 +318,85 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
 end
 
 
+function rhs_covariant!(cache, mesh, equation, P)
+
+  @unpack D, S, tau             = get_dynamic_variables(cache)
+  @unpack flx_D, flx_S, flx_tau,
+          ldg_D, ldg_S, ldg_tau,
+          smoothed_mu, max_v,
+          rho, v, eps, p        = get_static_variables(cache)
+  @unpack rhs_D, rhs_S, rhs_tau = get_rhs_variables(cache)
+  @unpack nflx_D, nflx_S, nflx_tau,
+          bdry_D, bdry_S, bdry_tau,
+          bdry_rho, bdry_v, bdry_eps,
+          bdry_p, bdry_max_v,
+          bdry_ldg_D, bdry_ldg_S, bdry_ldg_tau,
+          bdry_rhs_D, bdry_rhs_S, bdry_rhs_tau,
+          bdry_smoothed_mu      = get_bdry_variables(cache)
+  @unpack x = get_static_variables(cache)
+
+  # compute rhs for ∂t u + ∂x f(x) = 0
+
+  # dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
+  dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
+  dg1d.new_broadcast_volume!(maxspeed, equation, mesh)
+
+  dg1d.interpolate_face_data!(mesh, D,     bdry_D)
+  dg1d.interpolate_face_data!(mesh, S,     bdry_S)
+  dg1d.interpolate_face_data!(mesh, tau,   bdry_tau)
+  dg1d.interpolate_face_data!(mesh, max_v, bdry_max_v)
+  dg1d.interpolate_face_data!(mesh, rho,   bdry_rho)
+  dg1d.interpolate_face_data!(mesh, v,     bdry_v)
+  dg1d.interpolate_face_data!(mesh, eps,   bdry_eps)
+  dg1d.interpolate_face_data!(mesh, p,     bdry_p)
+
+  dg1d.new_broadcast_volume!(flux, equation, mesh)
+
+  # numerical fluxes
+  dg1d.new_broadcast_faces!(llf,    equation, mesh)
+  dg1d.new_broadcast_bdry!(bdryllf, equation, mesh)
+
+  compute_rhs_weak_form!(rhs_D,   flx_D,   nflx_D,   mesh)
+  compute_rhs_weak_form!(rhs_S,   flx_S,   nflx_S,   mesh)
+  compute_rhs_weak_form!(rhs_tau, flx_tau, nflx_tau, mesh)
+
+  # now compute rhs for ∂t u + ∂x f(x) = ∂x( μ (-∂t+∂x) u )
+
+  ## solve auxiliary equation: q + ∂x u = 0
+  dg1d.new_broadcast_faces!(ldg_nflux,  equation, mesh)
+  dg1d.new_broadcast_bdry!(ldg_bdryllf, equation, mesh)
+  compute_rhs_weak_form!(ldg_D,   D,   nflx_D,   mesh)
+  compute_rhs_weak_form!(ldg_S,   S,   nflx_S,   mesh)
+  compute_rhs_weak_form!(ldg_tau, tau, nflx_tau, mesh)
+  # fd_diff_1d!(ldg_D,   D,   x, mesh)
+  # fd_diff_1d!(ldg_S,   S,   x, mesh)
+  # fd_diff_1d!(ldg_tau, tau, x, mesh)
+
+  # need to recompute, because ldg computation overwrites nflx_{D,S,tau}
+  # numerical fluxes
+  dg1d.new_broadcast_faces!(llf,    equation, mesh)
+  dg1d.new_broadcast_bdry!(bdryllf, equation, mesh)
+
+  dg1d.new_broadcast_volume!(av_flux_covariant, equation, mesh)
+
+  dg1d.interpolate_face_data!(mesh, ldg_D,       bdry_ldg_D)
+  dg1d.interpolate_face_data!(mesh, ldg_S,       bdry_ldg_S)
+  dg1d.interpolate_face_data!(mesh, ldg_tau,     bdry_ldg_tau)
+  dg1d.interpolate_face_data!(mesh, rhs_D,       bdry_rhs_D)
+  dg1d.interpolate_face_data!(mesh, rhs_S,       bdry_rhs_S)
+  dg1d.interpolate_face_data!(mesh, rhs_tau,     bdry_rhs_tau)
+  dg1d.interpolate_face_data!(mesh, smoothed_mu, bdry_smoothed_mu)
+
+  dg1d.new_broadcast_faces!(av_nflux_covariant, equation, mesh)
+
+  compute_rhs_weak_form!(rhs_D,   flx_D,   nflx_D,   mesh)
+  compute_rhs_weak_form!(rhs_S,   flx_S,   nflx_S,   mesh)
+  compute_rhs_weak_form!(rhs_tau, flx_tau, nflx_tau, mesh)
+
+  return
+end
+
+
 #######################################################################
 #                                 2d                                  #
 #######################################################################

src/SRHD/setup.jl

@@ -9,7 +9,8 @@ function Project(env::Environment, prms, mesh::Mesh1d)
   tci  = TCI.make_TCI(env.mesh, prms["TCI"])
   hrsc = HRSC.make_HRSC(env.mesh, prms["HRSC"])
 
-  P = Project(equation, hrsc, tci)
+  fixedprms = (; av_regularization=:covariant, id_smooth=true)
+  P = Project(equation, hrsc, tci, fixedprms)
 
   # register variables
   # TODO Somehow replace _register_variables! with register_variables!
@@ -55,7 +56,7 @@ function Project(env::Environment, prms, mesh::Mesh2d)
   # TODO Somehow replace _register_variables! with register_variables!
   @unpack cache = env
   _register_variables!(mesh, P)
-  _register_variables!(mesh, hrsc)
+  _register_variables!(mesh, hrsc, nothing)
   _register_variables!(mesh, bdrycond)
   display(cache)
 
@@ -108,23 +109,23 @@ function _register_variables!(mesh::Mesh, P::Project)
     cell_variablenames    = (:cellmax_v,),
     global_variablenames  = (:t, :tm1) # time steps
   )
-  _register_variables!(mesh, P.hrsc)
+  _register_variables!(mesh, P.hrsc, P.prms)
   _register_variables!(mesh, P.tci)
   @unpack x = get_static_variables(mesh.cache)
   # @. x = env.mesh.x[:]
 end
 
 
-function _register_variables!(mesh, tci_or_hrsc::Nothing) end
+function _register_variables!(mesh, tci_or_hrsc::Nothing, args...) end
 
 
-_register_variables!(mesh, hrsc::AbstractHRSC) =
+_register_variables!(mesh, hrsc::AbstractHRSC, prms) =
   error("_register_variables: Not implemented for hrsc of type '$(typeof(hrsc))'")
 
-_register_variables!(mesh, hrsc::HRSC.AbstractReconstruction) = nothing
+_register_variables!(mesh, hrsc::HRSC.AbstractReconstruction, prms) = nothing
 
 
-function _register_variables!(mesh, hrsc::HRSC.AbstractArtificialViscosity)
+function _register_variables!(mesh, hrsc::HRSC.AbstractArtificialViscosity, prms)
   register_variables!(mesh.cache,
     static_variablenames  = (:ldg_D, :ldg_S, :ldg_tau,
                              :flx_ldg_D, :flx_ldg_S, :flx_ldg_tau),
@@ -133,17 +134,23 @@ function _register_variables!(mesh, hrsc::HRSC.AbstractArtificialViscosity)
 end
 
 
-function _register_variables!(mesh, hrsc::HRSC.SmoothedArtificialViscosity)
-  _register_variables!(mesh, hrsc.av)
+function _register_variables!(mesh, hrsc::HRSC.SmoothedArtificialViscosity, prms)
+  _register_variables!(mesh, hrsc.av, prms)
   register_variables!(mesh.cache,
     static_variablenames  = (:smoothed_mu,),
     bdry_variablenames    = (:bdry_ldg_D, :bdry_ldg_S, :bdry_ldg_tau,
                              :bdry_smoothed_mu),
   )
+  reg = prms.av_regularization
+  if reg == :covariant
+    register_variables!(mesh.cache,
+      bdry_variablenames    = (:bdry_rhs_D, :bdry_rhs_S, :bdry_rhs_tau)
+    )
+  end
 end
 
 
-function _register_variables!(mesh::Mesh2d, hrsc::HRSC.EntropyViscosity)
+function _register_variables!(mesh::Mesh2d, hrsc::HRSC.EntropyViscosity, prms)
   register_variables!(mesh.cache,
     static_variablenames = (:E, :Em1, :flx_E_x, :flx_E_y, :flxm1_E_x, :flxm1_E_y,
                             :ldg_D_x, :ldg_D_y, :ldg_Sx_x, :ldg_Sx_y,
@@ -156,7 +163,7 @@ function _register_variables!(mesh::Mesh2d, hrsc::HRSC.EntropyViscosity)
 end
 
 
-function _register_variables!(mesh, tci::TCI.AbstractTCI)
+function _register_variables!(mesh, tci::TCI.AbstractTCI, prms)
   register_variables!(mesh.cache, cell_variablenames = (:flag,:D_flag,:S_flag,:tau_flag))
 end
 

src/SRHD/types.jl

@@ -21,11 +21,13 @@ end
 
 
 struct Project{T_HRSC     <:Maybe{HRSC.AbstractHRSC},
-               T_TCI      <:Maybe{TCI.AbstractTCI}}
+               T_TCI      <:Maybe{TCI.AbstractTCI},
+               T_Prms}
 
   equation::Equation
   hrsc::T_HRSC
   tci::T_TCI
+  prms::T_Prms
 end