src/EulerEq/EulerEq.jl

@@ -12,7 +12,6 @@ include("equation.jl")
 include("rhs.jl")
 include("initialdata.jl")
 include("callbacks.jl")
-include("boundaryconditions.jl")
 include("setup.jl")
 
 

src/EulerEq/boundaryconditions.jl

-function make_BoundaryConditions(env, eq, rsolver, prms)
-
-  @unpack bc = prms
-
-  lhs_bc, rhs_bc = if bc == "periodic"
-    PeriodicBC(), PeriodicBC()
-  elseif bc == "from_id"
-    @unpack rho, q, E = get_dynamic_variables(env.cache)
-    lhs_state = (rho[1], q[1], E[1])
-    rhs_state = (rho[end], q[end], E[end])
-    (DirichletBC(lhs_state), DirichletBC(rhs_state))
-  else
-    error("Unknown boundary condition requested: '$bc'")
-  end
-
-  return BoundaryConditions(lhs_bc, rhs_bc, rsolver)
-end
-
-
-function make_BoundaryConditions_LDG(env, eq, ldg_rsolver, prms)
-
-  @unpack bc = prms
-
-  isnothing(ldg_rsolver) && return nothing
-
-  lhs_bc, rhs_bc = if type == "periodic"
-    PeriodicBC(), PeriodicBC()
-  elseif bc == "from_id"
-    OutflowBC(), OutflowBC()
-  else
-    error("Unknown boundary condition requested: '$bc'")
-  end
-
-  return BoundaryConditions(lhs_bc, rhs_bc, ldg_rsolver)
-end

src/EulerEq/equation.jl

@@ -15,6 +15,33 @@
   @returns flx_rho, flx_q, flx_E
 end
 
+@with_signature [legacy=false] function fv_bdry_flux(equation::EulerEquation)
+  @accepts rho, flx_rho, q, flx_q, E, flx_E, v
+  @accepts init_rho, init_q, init_E, init_v
+
+  # enforce dirichlet/neumann conditions following stability analysis
+  # see dg-notes/tex/evm.pdf
+  # Hesthaven, Numerical Methods for Conservation Laws uses same logic
+
+  # Dirichlet conditions
+  bdry_rho = -rho + 2*init_rho
+  bdry_q   = -q   + 2*init_q
+  bdry_E   = -E   + 2*init_E
+
+  # Neumann conditions
+  # bdry_rho = rho
+  # bdry_q   = q
+  # bdry_E   = E
+
+  @unpack eos  = equation
+  bdry_eps     = bdry_E / bdry_rho - bdry_q^2 / (2 * bdry_rho^2)
+  bdry_p       = eos(Pressure, Density, bdry_rho, InternalEnergy, bdry_eps)
+  nflx_rho = bdry_q
+  nflx_q   = bdry_q^2 / bdry_rho + bdry_p
+  nflx_E   = bdry_q / bdry_rho * (bdry_E + bdry_p)
+
+  @returns [bdry] bdry_rho, bdry_q, bdry_E, nflx_rho, nflx_q, nflx_E
+end
 
 @with_signature [legacy=false] function maxspeed(equation::EulerEquation)
   @accepts rho, q, E

src/EulerEq/initialdata.jl

@@ -179,8 +179,8 @@ end
 #   @unpack S, Sm1, flx_S, flx_Sm1 = get_static_variables(cache)
 #   tm1 .= 0.0
 #   t   .= 0.0
-#   broadcast_volume!(entropy, equation, cache)
-#   broadcast_volume!(entropy_flux, equation, cache)
+#   broadcast_volume_2!(entropy, equation, mesh)
+#   broadcast_volume_2!(entropy_flux, equation, mesh)
 #   @. Sm1 = S
 #   @. flx_Sm1 = flx_S
 # end

src/EulerEq/rhs.jl

@@ -2,26 +2,26 @@
 #                              Timestep                               #
 #######################################################################
 
-
-function timestep(env, P::Project, hrsc::Maybe{HRSC.AbstractHRSC}, ::Mesh1d)
-  @unpack cache, mesh = env
-  @unpack equation = P
-
+function compute_maxspeed(mesh, equation, cache)
+  @unpack v = get_static_variables(cache)
   dg1d.new_broadcast_volume!(maxspeed, equation, mesh)
-  v = get_variable(cache, :v)
-  vmax = dg1d.absolute_maximum(view(v,:))
+  vmax = dg1d.absolute_maximum(v)
   vmax_limit = 1e4
   if vmax > vmax_limit
     @warn "Limiting timestep due to maximum speed exceeding $vmax_limit"
     vmax = vmax_limit
   end
+  return vmax
+end
 
-  @unpack dl, x, CFL, element = mesh
-  @unpack N = element
-
-  dt = CFL * dl / (N * vmax)
+function timestep(env, P::Project, hrsc::Maybe{HRSC.AbstractHRSC}, mesh::Mesh1d)
+  vmax = compute_maxspeed(mesh, P.equation, mesh.cache)
+  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, ::Mesh1d)
@@ -39,10 +39,10 @@ function timestep(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, ::Mes
   smoothed_mu = get_variable(cache, :smoothed_mu)
   mumax = dg1d.absolute_maximum(view(smoothed_mu, :))
 
-  @unpack dl, x, CFL, element = mesh
+  @unpack dl, x, element = mesh
   @unpack N = element
 
-  dt = CFL / (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"
@@ -65,7 +65,7 @@ function timestep(env, P::Project2d, hrsc::Maybe{HRSC.AbstractHRSC}, ::Mesh2d)
     vmax = vmax_limit
   end
 
-  @unpack CFL, element = mesh
+  @unpack element = mesh
   @unpack N, z = element
 
   dx, dy = dg1d.widths(mesh.boxes[1])
@@ -77,7 +77,7 @@ function timestep(env, P::Project2d, hrsc::Maybe{HRSC.AbstractHRSC}, ::Mesh2d)
   dz = z[2]-z[1]
   dl = min(dx, dy) * dz
 
-  dt = CFL * dl / (N * vmax)
+  dt = dl / (N * vmax)
 
   return dt
 end
@@ -91,7 +91,7 @@ end
 function rhs!(env, P::Project, hrsc::Nothing, ::Mesh1d, bdryconds, ldg_bdryconds, av_bdryconds)
 
   @unpack cache, mesh                   = env
-  @unpack equation, rsolver             = P
+  @unpack equation                      = P
   @unpack rho, q, E                     = get_dynamic_variables(cache)
   @unpack flx_rho, flx_q, flx_E, p, eps = get_static_variables(cache)
   @unpack rhs_rho, rhs_q, rhs_E         = get_rhs_variables(cache)
@@ -117,11 +117,33 @@ function rhs!(env, P::Project, hrsc::Nothing, ::Mesh1d, bdryconds, ldg_bdryconds
 end
 
 
+function rhs!(env, P::Project, hrsc::Nothing, ::Mesh1d{FVElement}, bdryconds, ldg_bdryconds, av_bdryconds)
+
+  @unpack cache, mesh              = env
+  @unpack equation                 = P
+  @unpack rho, q, E                = get_dynamic_variables(cache)
+  @unpack flx_rho, flx_q, flx_E    = get_static_variables(cache)
+  @unpack rhs_rho, rhs_q, rhs_E    = get_rhs_variables(cache)
+  @unpack nflx_rho, nflx_q, nflx_E = get_bdry_variables(cache)
+  @unpack bdry_rho, bdry_q, bdry_E = get_bdry_variables(cache)
+
+  dg1d.new_broadcast_volume!(flux, equation, mesh)
+  dg1d.new_broadcast_bdry!(fv_bdry_flux, equation, mesh)
+
+  dt = timestep(env, P, hrsc, mesh)
+  fv_update_step!(rhs_rho, rho, flx_rho, bdry_rho, nflx_rho, dt, mesh)
+  fv_update_step!(rhs_q,   q,   flx_q,   bdry_q,   nflx_q,   dt, mesh)
+  fv_update_step!(rhs_E,   E,   flx_E,   bdry_E,   nflx_E,   dt, mesh)
+
+  return
+end
+
+
+
 function rhs!(env, P::Project2d, hrsc::Nothing, ::Mesh2d, bdryconds, ldg_bdryconds, av_bdryconds)
 
   @unpack cache, mesh                        = env
-  @unpack equation, rsolver_rho, rsolver_qx,
-          rsolver_qy, rsolver_E              = P
+  @unpack equation                           = P
   @unpack rho, qx, qy, E                     = get_dynamic_variables(cache)
   @unpack flx_rho_x, flx_rho_y,
           flx_qx_x, flx_qx_y,
@@ -155,7 +177,7 @@ end
 function rhs!(env, P::Project, hrsc::AbstractReconstruction, ::Mesh1d, bdryconds, ldg_bdryconds, av_bdryconds)
 
   @unpack cache, mesh                   = env
-  @unpack equation, rsolver             = P
+  @unpack equation                      = P
   @unpack rho, q, E                     = get_dynamic_variables(cache)
   @unpack flx_rho, flx_q, flx_E, p, eps = get_static_variables(cache)
   @unpack flag                          = get_cell_variables(cache)
@@ -191,7 +213,7 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, ::Mesh1d,
     bdryconds, ldg_bdryconds, av_bdryconds)
 
   @unpack cache, mesh = env
-  @unpack equation, rsolver, ldg_rsolver, av_rsolver = P
+  @unpack equation    = P
 
   @unpack rho, q, E                            = get_dynamic_variables(cache)
   @unpack flx_rho, flx_q, flx_E, p, eps,

src/EulerEq/setup.jl

@@ -8,11 +8,9 @@ function Project(env::Environment, mesh::Mesh1d, prms)
   equation = EulerEquation(eos)
   tci  = TCI.make_TCI(env.mesh, prms["TCI"])
   hrsc = HRSC.make_HRSC(env.mesh, prms["HRSC"])
-  rsolver = nothing
-  ldg_rsolver, av_rsolver = nothing, nothing
 
   bdrycond = dg1d.DirichletBC2()
-  P = Project(equation, rsolver, hrsc, tci, ldg_rsolver, av_rsolver, bdrycond)
+  P = Project(equation, hrsc, tci, bdrycond)
 
   # register variables
   # TODO add a ::Nothing overload for register_variables!
@@ -50,11 +48,9 @@ function Project(env::Environment, mesh::Mesh2d, prms)
   eos = EquationOfState.make_EquationOfState(prms["EquationOfState"])
   hrsc = HRSC.make_HRSC(mesh, prms["HRSC"])
   equation = EulerEquation2d(eos)
-  ldg_rsolver, av_rsolver = nothing, nothing
 
   bdrycond = dg1d.DirichletBC2()
-  P = Project2d(equation, nothing, nothing,
-                nothing, nothing, hrsc, nothing, nothing, nothing)
+  P = Project2d(equation, hrsc, nothing)
 
   # register variables
   # TODO Somehow replace _register_variables! with register_variables!

src/EulerEq/types.jl

@@ -8,39 +8,21 @@ end
 
 
 mutable struct Project{T_Equation <:EulerEquation,
-                       T_RSolver  <:Maybe{AbstractRiemannSolver},
                        T_HRSC     <:Maybe{HRSC.AbstractHRSC},
                        T_TCI      <:Maybe{TCI.AbstractTCI},
-                       T_LDG_RSolver<:Maybe{AbstractRiemannSolver},
-                       T_AV_RSolver<:Maybe{AbstractRiemannSolver},
                        T_BC}
 
   equation::T_Equation
-  rsolver::T_RSolver
   hrsc::T_HRSC
   tci::T_TCI
-  ldg_rsolver::T_LDG_RSolver
-  av_rsolver::T_AV_RSolver
   bdrycond::T_BC
 end
 
 mutable struct Project2d{T_Equation <:EulerEquation2d,
-                         T_RSolver1  <:Maybe{AbstractRiemannSolver},
-                         T_RSolver2  <:Maybe{AbstractRiemannSolver},
-                         T_RSolver3  <:Maybe{AbstractRiemannSolver},
-                         T_RSolver4  <:Maybe{AbstractRiemannSolver},
                          T_HRSC     <:Maybe{HRSC.AbstractHRSC},
-                         T_TCI      <:Maybe{TCI.AbstractTCI},
-                         T_LDG_RSolver<:Maybe{AbstractRiemannSolver},
-                         T_AV_RSolver<:Maybe{AbstractRiemannSolver}}
+                         T_TCI      <:Maybe{TCI.AbstractTCI}}
 
   equation::T_Equation
-  rsolver_rho::T_RSolver1
-  rsolver_qx::T_RSolver2
-  rsolver_qy::T_RSolver3
-  rsolver_E::T_RSolver4
   hrsc::T_HRSC
   tci::T_TCI
-  ldg_rsolver::T_LDG_RSolver
-  av_rsolver::T_AV_RSolver
 end

src/GRHD/rhs.jl

@@ -17,10 +17,10 @@ function timestep(env, P::Project, hrsc::Maybe{HRSC.AbstractHRSC})
     vmax = vmax_limit
   end
 
-  @unpack dl, x, CFL, element = mesh
+  @unpack dl, x, element = mesh
   @unpack N = element
 
-  dt = CFL * dl / (N^2 * vmax)
+  dt = dl / (N^2 * vmax)
   return dt
 end
 
@@ -41,10 +41,10 @@ function timestep(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity)
   smoothed_mu = get_variable(cache, :smoothed_mu)
   mumax = dg1d.absolute_maximum(view(smoothed_mu, :))
 
-  @unpack dl, x, CFL, element = mesh
+  @unpack dl, x, element = mesh
   @unpack N = element
 
-  dt = CFL / (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"

src/HRSC/ArtificialViscosity.jl

@@ -7,7 +7,7 @@ Supertype for Artificial Viscosity (AV) interface.
 Artificial viscosity aims to resolve shocks and discontinuities
 by regularizing a first order hyperbolic conservation law with a second
 order spatial derivative term, e.g.
-  
+
   ∂t u + ∂x f(u) = 0   =>   ∂t u + ∂x f(u) = ∂x(ν(x) ∂x u)
 """
 abstract type AbstractArtificialViscosity <: AbstractHRSC end
@@ -20,7 +20,8 @@ Compute artificial viscosity `mu` inplace based on current solution `u`, local s
 time `t`.
 """
 function timestep(vmax, mumax, mesh::Mesh, av::AbstractArtificialViscosity)
-  dt = CFL / (vmax * N^2 / dl + mumax * N^4 / dl^2)
+  TODO("deprecated")
+  dt = 1 / (vmax * N^2 / dl + mumax * N^4 / dl^2)
   return dt
 end
 

src/HeatEq/callbacks.jl

-function make_callback(env, P::Project, bdryconds::BoundaryConditions) end
+function make_callback(env, P::Project, isperiodic) end

src/HeatEq/rhs.jl

@@ -15,13 +15,13 @@ function timestep(env, P::Project, ::Mesh1d)
     vmax = vmax_limit
   end
 
-  @unpack CFL, element = mesh
+  @unpack element = mesh
   @unpack N, z = element
   dz = z[2]-z[1]
   dx, = minimum(dg1d.widths(mesh.boxes[1]))
   dl = dx * dz
 
-  dt = CFL * dl^2 / (N^2 * vmax) # just a guess to use 1/N^2 here since we have 2nd derivatives
+  dt = dl^2 / (N^2 * vmax) # just a guess to use 1/N^2 here since we have 2nd derivatives
   return dt
 end
 

src/SRHD/SRHD.jl

@@ -9,7 +9,6 @@ include("equation.jl")
 include("rhs.jl")
 include("initialdata.jl")
 include("callbacks.jl")
-include("boundaryconditions.jl")
 include("setup.jl")
 
 

src/SRHD/boundaryconditions.jl

-function make_BoundaryConditions(env, eq, rsolver, prms)
-
-  @unpack mesh = env
-  @unpack bc = prms
-  @unpack xmin, xmax = mesh
-
-  lhs_bc, rhs_bc = if bc == "periodic"
-    PeriodicBC(), PeriodicBC()
-  elseif bc == "from_id"
-    @unpack D, S, tau = get_dynamic_variables(env.mesh.cache)
-    p0 = 0.0
-    DStau_lhs = (D[1], S[1], tau[1])
-    (p_lhs,) = find_pressure_insane((DStau_lhs..., p0, xmin), eq)
-    (rho_lhs, v_lhs, eps_lhs) = primitives((DStau_lhs..., p_lhs), eq)
-    lhs_state = (DStau_lhs..., p_lhs, rho_lhs, v_lhs, eps_lhs)
-    DStau_rhs = (D[end], S[end], tau[end])
-    (p_rhs,) = find_pressure_insane((DStau_rhs..., p0, xmax), eq)
-    (rho_rhs, v_rhs, eps_rhs) = primitives((DStau_rhs..., p_rhs), eq)
-    rhs_state = (DStau_rhs..., p_rhs, rho_rhs, v_rhs, eps_rhs)
-    (DirichletBC(lhs_state), DirichletBC(rhs_state))
-  else
-    error("Unknown boundary condition requested: '$bc'")
-  end
-
-  return BoundaryConditions(lhs_bc, rhs_bc, rsolver)
-end
-
-
-function make_BoundaryConditions_LDG(env, eq, ldg_rsolver, prms)
-
-  @unpack bc = prms
-
-  isnothing(ldg_rsolver) && return nothing
-
-  lhs_bc, rhs_bc = if bc == "periodic"
-    PeriodicBC(), PeriodicBC()
-  elseif bc == "no_inflow"
-    OutflowBC(), OutflowBC()
-  elseif bc == "from_id"
-    OutflowBC(), OutflowBC()
-  else
-    error("Unknown boundary condition requested: '$bc'")
-  end
-
-  return BoundaryConditions(lhs_bc, rhs_bc, ldg_rsolver)
-end

src/SRHD/callbacks.jl

@@ -13,14 +13,14 @@ function make_callback(P)
 end
 
 
-function make_callback(env, P::Project, bdryconds::BoundaryConditions)
-  cbfn_equation(u, t) = callback_equation(u, t, env, P, isperiodic(bdryconds), P.equation)
+function make_callback(env, P::Project, isperiodic)
+  cbfn_equation(u, t) = callback_equation(u, t, env, P, isperiodic, P.equation)
   cb_equation         = FunctionCallback(cbfn_equation,
                                          CallbackTiming(every_iteration=1,every_dt=0))
-  cbfn_tci(u, t)      = callback_tci(u, t, env, P, isperiodic(bdryconds), P.tci)
+  cbfn_tci(u, t)      = callback_tci(u, t, env, P, isperiodic, P.tci)
   cb_tci              = FunctionCallback(cbfn_tci,
                                          CallbackTiming(every_iteration=1, every_dt=0))
-  cbfn_hrsc(u, t)     = callback_hrsc(u, t, env, P, isperiodic(bdryconds), P.hrsc, env.mesh)
+  cbfn_hrsc(u, t)     = callback_hrsc(u, t, env, P, isperiodic, P.hrsc, env.mesh)
   cb_hrsc             = FunctionCallback(cbfn_hrsc,
                                          CallbackTiming(every_iteration=1,every_dt=0))
   callbackset = CallbackSet(cb_equation, cb_tci, cb_hrsc)
@@ -79,7 +79,7 @@ function callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc::HRSC.Abstract
   @unpack max_v                 = get_static_variables(cache)
   @unpack mu, cellmax_v         = get_cell_variables(cache)
   @unpack D                     = get_dynamic_variables(cache)
-  broadcast_volume!(speed, equation, cache)
+  broadcast_volume_2!(speed, equation, mesh)
   Npts, K = layout(mesh)
   mat_cellmax_v = reshape(view(max_v, :), (Npts, K))
   for k = 1:K

src/SRHD/equation.jl

@@ -179,6 +179,118 @@ end
 end
 
 
+@with_signature function llf(equation::Equation)
+  @accepts Prefix(lhs), D, flx_D, S, flx_S, tau, flx_tau, max_v
+  @accepts Prefix(rhs), D, flx_D, S, flx_S, tau, flx_tau, max_v
+  @accepts nx
+
+  vmax = max(lhs_max_v, rhs_max_v)
+
+  lhs_nflx = nx*lhs_flx_D
+  rhs_nflx = nx*rhs_flx_D
+  nflx_D = LLF(lhs_nflx, rhs_nflx, lhs_D, rhs_D, vmax)
+
+  lhs_nflx = nx*lhs_flx_S
+  rhs_nflx = nx*rhs_flx_S
+  nflx_S = LLF(lhs_nflx, rhs_nflx, lhs_S, rhs_S, vmax)
+
+  lhs_nflx = nx*lhs_flx_tau
+  rhs_nflx = nx*rhs_flx_tau
+  nflx_tau = LLF(lhs_nflx, rhs_nflx, lhs_tau, rhs_tau, vmax)
+
+  @returns nflx_D, nflx_S, nflx_tau
+end
+
+
+# @with_signature function bdry_llf(equation::Equation2d)
+#   @accepts D, flx_D_x, flx_D_y, Sx, flx_Sx_x, flx_Sx_y, Sy, flx_Sy_x, flx_Sy_y, tau, flx_tau_x, flx_tau_y, max_v
+#   @accepts init_D, init_flx_D_x, init_flx_D_y, init_Sx, init_flx_Sx_x, init_flx_Sx_y, init_Sy, init_flx_Sy_x,
+#            init_flx_Sy_y, init_tau, init_flx_tau_x, init_flx_tau_y, init_max_v
+#   @accepts nx, ny
+#
+#   vmax = max(max_v, init_max_v)
+#
+#   nflx      = nx*flx_D_x + ny*flx_D_y
+#   init_nflx = nx*init_flx_D_x + ny*init_flx_D_y
+#   nflx_D    = LLF(nflx, init_nflx, D, init_D, vmax)
+#
+#   nflx      = nx*flx_Sx_x + ny*flx_Sx_y
+#   init_nflx = nx*init_flx_Sx_x + ny*init_flx_Sx_y
+#   nflx_Sx   = LLF(nflx, init_nflx, Sx, init_Sx, vmax)
+#
+#   nflx      = nx*flx_Sy_x + ny*flx_Sy_y
+#   init_nflx = nx*init_flx_Sy_x + ny*init_flx_Sy_y
+#   nflx_Sy   = LLF(nflx, init_nflx, Sy, init_Sy, vmax)
+#
+#   nflx      = nx*flx_tau_x + ny*flx_tau_y
+#   init_nflx = nx*init_flx_tau_x + ny*init_flx_tau_y
+#   nflx_tau  = LLF(nflx, init_nflx, tau, init_tau, vmax)
+#
+#   @returns nflx_D, nflx_Sx, nflx_Sy, nflx_tau
+# end
+
+
+@with_signature function ldg_nflux(eq::Equation)
+  @accepts Prefix(lhs), D, S, tau
+  @accepts Prefix(rhs), D, S, tau
+  @accepts nx
+
+  nflx_D   = 0.5 * nx * (lhs_D + rhs_D)
+  nflx_S   = 0.5 * nx * (lhs_S + rhs_S)
+  nflx_tau = 0.5 * nx * (lhs_tau + rhs_tau)
+
+  @returns nflx_D, nflx_S, nflx_tau
+end
+@with_signature function bdry_ldg_nflux(eq::Equation)
+  @accepts Prefix(lhs), D, S, tau
+  @accepts Prefix(rhs), D, S, tau
+  @accepts nx, ny
+
+  nflx_D = nx * lhs_D
+  nflx_S = nx * lhs_S
+  nflx_tau = ny * lhs_tau
+
+  @returns nflx_D, nflx_S, nflx_tau
+end
+
+
+@with_signature function av_nflux(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
+
+  lhs_nflx = nx*lhs_ldg_D
+  rhs_nflx = nx*rhs_ldg_D
+  nflx_D += 0.5 * (lhs_smoothed_mu*lhs_nflx + rhs_smoothed_mu*rhs_nflx)
+
+  lhs_nflx = nx*lhs_ldg_S
+  rhs_nflx = nx*rhs_ldg_S
+  nflx_S += 0.5 * (lhs_smoothed_mu*lhs_nflx + rhs_smoothed_mu*rhs_nflx)
+
+  lhs_nflx = nx*lhs_ldg_tau
+  rhs_nflx = nx*rhs_ldg_tau
+  nflx_tau += 0.5 * (lhs_smoothed_mu*lhs_nflx + rhs_smoothed_mu*rhs_nflx)
+
+  @returns nflx_D, nflx_S, nflx_tau
+end
+@with_signature function bdry_av_nflux(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
+
+  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
+
+
 #######################################################################
 #                                 2d                                  #
 #######################################################################

src/SRHD/initialdata.jl

@@ -37,7 +37,7 @@ function initialdata_equation_of_state!(env, P::Project, eos::EquationOfState.Id
   @. S   = Formulae.S(ρ0, v0, ϵ0, p0)
   @. tau = Formulae.τ(ρ0, v0, ϵ0, p0)
 
-  broadcast_volume!(cons2prim, P.equation, cache)
+  broadcast_volume_2!(cons2prim, P.equation, mesh)
 
   return
 end

src/SRHD/rhs.jl

@@ -10,8 +10,8 @@ function timestep(env, P::Project, hrsc::Maybe{HRSC.AbstractHRSC}, mesh::Mesh1d)
   @unpack equation    = P
   @unpack max_v = get_static_variables(cache)
 
-  broadcast_volume!(cons2prim, equation, cache)
-  broadcast_volume!(speed, equation, cache)
+  broadcast_volume_2!(cons2prim, equation, mesh)
+  broadcast_volume_2!(speed, equation, mesh)
   vmax = dg1d.absolute_maximum(max_v)
   # vmax > 0.99 && error("Unreasonable maximum speed encountered: $vmax")
   vmax_limit = 0.9999
@@ -20,10 +20,10 @@ function timestep(env, P::Project, hrsc::Maybe{HRSC.AbstractHRSC}, mesh::Mesh1d)
     vmax = vmax_limit
   end
 
-  @unpack dl, x, CFL, element = mesh
+  @unpack dl, x, element = mesh
   @unpack N = element
 
-  dt = CFL * dl / (N^2 * vmax)
+  dt = dl / (N^2 * vmax)
   return dt
 end
 
@@ -32,8 +32,8 @@ function timestep(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, mesh:
   @unpack cache, mesh = env
   @unpack equation    = P
 
-  broadcast_volume!(cons2prim, equation, cache)
-  broadcast_volume!(speed, equation, cache)
+  broadcast_volume_2!(cons2prim, equation, mesh)
+  broadcast_volume_2!(speed, equation, mesh)
   max_v = get_variable(cache, :max_v)
   vmax = dg1d.absolute_maximum(view(max_v,:))
   vmax_limit = 0.99
@@ -44,10 +44,10 @@ function timestep(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, mesh:
   smoothed_mu = get_variable(cache, :smoothed_mu)
   mumax = dg1d.absolute_maximum(view(smoothed_mu, :))
 
-  @unpack dl, x, CFL, element = mesh
+  @unpack dl, x, element = mesh
   @unpack N = element
 
-  dt = CFL / (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"
@@ -78,7 +78,7 @@ function timestep(env, P, hrsc::Maybe{HRSC.AbstractHRSC}, mesh::Mesh2d)
     vmax = vmax_limit
   end
 
-  @unpack CFL, element = mesh
+  @unpack element = mesh
   @unpack N, z = element
 
   dx, dy = dg1d.widths(mesh.boxes[1])
@@ -90,7 +90,7 @@ function timestep(env, P, hrsc::Maybe{HRSC.AbstractHRSC}, mesh::Mesh2d)
   dz = z[2]-z[1]
   dl = min(dx, dy) * dz
 
-  dt = CFL * dl / (N * vmax)
+  dt = dl / (N * vmax)
 
   return dt
 end
@@ -107,16 +107,15 @@ rhs!(env, P, hrsc, bdryconds, ldg_bdryconds, av_bdryconds) =
 function rhs!(env, P::Project, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryconds, mesh::Mesh1d)
 
   @unpack cache, mesh              = env
-  @unpack equation, rsolver        = P
+  @unpack equation                 = P
   @unpack D, S, tau                = get_dynamic_variables(cache)
   @unpack flx_D, flx_S, flx_tau    = get_static_variables(cache)
   @unpack rhs_D, rhs_S, rhs_tau    = get_rhs_variables(cache)
   @unpack nflx_D, nflx_S, nflx_tau = get_bdry_variables(cache)
 
-  broadcast_volume!(cons2prim, equation, cache)
-  broadcast_volume!(flux, equation, cache)
-  broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
-  # broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
+  broadcast_volume_2!(cons2prim, equation, mesh)
+  broadcast_volume_2!(flux, equation, mesh)
+  broadcast_faces_2!(llf, equation, mesh)
 
   compute_rhs_weak_form!(rhs_D,   flx_D,    nflx_D,   mesh)
   compute_rhs_weak_form!(rhs_S,   flx_S,    nflx_S,   mesh)
@@ -130,7 +129,7 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractReconstruction,
     bdryconds, ldg_bdryconds, av_bdryconds, ::Mesh1d)
 
   @unpack cache, mesh              = env
-  @unpack equation, rsolver        = P
+  @unpack equation                 = P
   @unpack D, S, tau                = get_dynamic_variables(cache)
   @unpack flx_D, flx_S, flx_tau    = get_static_variables(cache)
   @unpack rhs_D, rhs_S, rhs_tau    = get_rhs_variables(cache)
@@ -141,9 +140,9 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractReconstruction,
   HRSC.reconstruct!(S,    flag, hrsc, isperiodic=isperiodic(bdryconds), limit_with_boundaries=false)
   HRSC.reconstruct!(tau,  flag, hrsc, isperiodic=isperiodic(bdryconds), m=1e-10, limit_with_boundaries=false)
 
-  broadcast_volume!(cons2prim, equation, cache)
-  broadcast_volume!(flux, equation, cache)
-  broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
+  broadcast_volume_2!(cons2prim, equation, mesh)
+  broadcast_volume_2!(flux, equation, mesh)
+  broadcast_faces_2!(llf, equation, mesh)
 
   compute_rhs_weak_form!(rhs_D,   flx_D,    nflx_D,   mesh)
   compute_rhs_weak_form!(rhs_S,   flx_S,    nflx_S,   mesh)
@@ -157,8 +156,7 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
     bdryconds, ldg_bdryconds, av_bdryconds, ::Mesh1d)
 
   @unpack cache, mesh                       = env
-  @unpack equation, rsolver, ldg_rsolver,
-          av_rsolver                        = P
+  @unpack equation                          = P
   @unpack D, S, tau                         = get_dynamic_variables(cache)
   @unpack flx_D, flx_S, flx_tau,
           ldg_D, ldg_S, ldg_tau,
@@ -168,22 +166,21 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
           nflx_ldg_S, nflx_ldg_tau          = get_bdry_variables(cache)
 
   ## solve auxiliary equation: q + ∂x u = 0
-  broadcast_volume!(ldg_flux, equation, cache)
-  broadcast_faces!(mean_flux, ldg_rsolver, cache, mesh)
-  broadcast_boundaryconditions!(mean_flux, ldg_bdryconds, cache, mesh, 0.0)
+  broadcast_volume_2!(ldg_flux, equation, mesh)
+  broadcast_faces_2!(ldg_nflux, equation, mesh)
   compute_rhs_weak_form!(ldg_D,   flx_ldg_D,   nflx_ldg_D,   mesh)
   compute_rhs_weak_form!(ldg_S,   flx_ldg_S,   nflx_ldg_S,   mesh)
   compute_rhs_weak_form!(ldg_tau, flx_ldg_tau, nflx_ldg_tau, mesh)
 
   ## compute rhs of regularized equation: ∂t u + ∂x f + ∂x μ q = 0
-  broadcast_volume!(av_flux, equation, cache)
-  broadcast_faces!(mean_flux, av_rsolver, cache, mesh)
-  broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, 0.0)
+  broadcast_volume_2!(av_flux, equation, mesh)
+  broadcast_faces_2!(av_nflux, equation, mesh)
+  # broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, 0.0)
 
-  broadcast_volume!(cons2prim, equation, cache)
-  broadcast_volume!(flux, equation, cache)
-  broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
-  broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
+  broadcast_volume_2!(cons2prim, equation, mesh)
+  broadcast_volume_2!(flux, equation, mesh)
+  broadcast_faces_2!(llf, equation, mesh)
+  # broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
 
   @. flx_D          += flx_ldg_D
   @. flx_S          += flx_ldg_S
@@ -284,13 +281,10 @@ function dg1d.rhs!(env, P::Project2d, hrsc::HRSC.AbstractArtificialViscosity, me
   broadcast_volume_2!(flux, equation, mesh)
   broadcast_faces_2!(llf, equation, mesh)
   dg1d.broadcast_bdry_2!(bdry_llf, equation, P.bdrycond, mesh)
-  # broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
 
   broadcast_volume_2!(av_flux_2d, equation, mesh)
   broadcast_faces_2!(av_nflux_2d, equation, mesh)
   dg1d.broadcast_bdry_2!(bdry_av_nflux_2d, equation, P.bdrycond, mesh)
-  # broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, 0.0)
-  # broadcast_faces_2!(mean_flux, equation, mesh, P.bdryconds)
 
   compute_rhs_weak_form!(rhs_D,   flx_D_x,   flx_D_y,   nflx_D,   mesh)
   compute_rhs_weak_form!(rhs_Sx,  flx_Sx_x,  flx_Sx_y,  nflx_Sx,  mesh)

src/SRHD/setup.jl

@@ -8,34 +8,17 @@ function Project(env::Environment, prms, mesh::Mesh1d)
   equation = make_Equation(eos, prms["SRHD"], mesh)
   tci  = TCI.make_TCI(env.mesh, prms["TCI"])
   hrsc = HRSC.make_HRSC(env.mesh, prms["HRSC"])
-  rsolver = ApproxRiemannSolver(flux, speed, equation)
-  ldg_rsolver, av_rsolver = if hrsc isa AbstractArtificialViscosity
-    ApproxRiemannSolver(ldg_flux, ldg_speed, equation),
-    ApproxRiemannSolver(av_flux, av_flux, equation)
-  else
-    nothing, nothing
-  end
 
-  P = Project(equation, rsolver, hrsc, tci, ldg_rsolver, av_rsolver)
+  P = Project(equation, hrsc, tci)
 
   # register variables
-  # TODO add a ::Nothing overload for register_variables!
   # TODO Somehow replace _register_variables! with register_variables!
   @unpack cache = env
   _register_variables!(mesh, P)
-  register_variables!(mesh, rsolver)
-  !isnothing(ldg_rsolver) && register_variables!(mesh, ldg_rsolver)
-  !isnothing(av_rsolver) && register_variables!(mesh, av_rsolver, dont_register=true)
   display(cache)
 
   # setup boundary conditions
-  bdryconds = make_BoundaryConditions(env, equation, rsolver, prms["SRHD"])
-  ldg_bdryconds, av_bdryconds = if hrsc isa AbstractArtificialViscosity
-    make_BoundaryConditions_LDG(env, equation, ldg_rsolver, prms["SRHD"]),
-    make_BoundaryConditions_LDG(env, equation, av_rsolver, prms["SRHD"])
-  else
-    nothing, nothing
-  end
+  bdryconds, ldg_bdryconds, av_bdryconds = nothing, nothing, nothing
 
   # TODO Bdry conditions should be moved to project and should only be just a 'lightweight'
   # descriptor. Furthermore, they should also register variables in cache themselves.
@@ -47,7 +30,7 @@ function Project(env::Environment, prms, mesh::Mesh1d)
   initialdata!(env, P, prms["SRHD"])
 
   # setup callbacks
-  projectcb = make_callback(env, P, bdryconds)
+  projectcb = make_callback(env, P, isperiodic(mesh))
 
   append!(env.callbacks, CallbackSet(projectcb.callbacks...))
 
@@ -114,6 +97,7 @@ function _register_variables!(mesh::Mesh, P::Project)
                              :max_v,                    # maximum charactersitic speed
                              :E, :Em1,                  # entropy
                              :flx_E, :flx_Em1),         # entropy flux
+    bdry_variablenames    = (:nflx_D, :nflx_S, :nflx_tau,),
     cell_variablenames    = (:cellmax_v,),
     global_variablenames  = (:t, :tm1) # time steps
   )

src/SRHD/types.jl

@@ -20,18 +20,12 @@ struct Equation2d{T_EoS<:EquationOfState.AbstractEquationOfState} <: AbstractEqu
 end
 
 
-struct Project{T_RSolver  <:AbstractRiemannSolver,
-               T_HRSC     <:Maybe{HRSC.AbstractHRSC},
-               T_TCI      <:Maybe{TCI.AbstractTCI},
-               T_LDG_RSolver<:Maybe{AbstractRiemannSolver},
-               T_AV_RSolver<:Maybe{AbstractRiemannSolver}}
+struct Project{T_HRSC     <:Maybe{HRSC.AbstractHRSC},
+               T_TCI      <:Maybe{TCI.AbstractTCI}}
 
   equation::Equation
-  rsolver::T_RSolver
   hrsc::T_HRSC
   tci::T_TCI
-  ldg_rsolver::T_LDG_RSolver
-  av_rsolver::T_AV_RSolver
 end
 
 

src/ScalarEq/ScalarEq.jl

@@ -10,7 +10,6 @@ include("equation.jl")
 include("rhs.jl")
 include("initialdata.jl")
 include("callbacks.jl")
-include("boundaryconditions.jl")
 include("setup.jl")