diff --git a/src/SRHD/SRHD.jl b/src/SRHD/SRHD.jl
index f049253e70378003b348b0bd780d9362b376f02d..5bcc785e1015605672d98c9605fd9a40835010b7 100644
--- a/src/SRHD/SRHD.jl
+++ b/src/SRHD/SRHD.jl
@@ -19,16 +19,17 @@ include("setup.jl")
# dg1d project interface definitions #
#######################################################################
-
dg1d.project_type(::Val{:SRHD}) = Project
# TODO Remove bdryconds
dg1d.rhs!(env::Environment, P::Project, bdryconds) = rhs!(env, P, P.hrsc, bdryconds...)
-dg1d.rhs!(env::Environment, P::Project2d, bdryconds) = rhs!(env, P, env.mesh)
+# dg1d.rhs!(env::Environment, P::Project2d, bdryconds) = rhs!(env, P, env.mesh)
dg1d.timestep!(env::Environment, P::Project) = timestep(env, P, P.hrsc)
dg1d.timestep!(env::Environment, P::Project2d) = timestep(env, P, nothing)
+# TODO Can we combine dg1d.project_type definition and dg1d.@parameters
+# dg1d.project_type(::Val{:SRHD}) = Project
dg1d.@parameters SRHD begin
"""
diff --git a/src/SRHD/equation.jl b/src/SRHD/equation.jl
index f8b794bb0d989371826c5fa2f9dc0c59f262cbd3..6e0b8457634486a93d5f3c2a45cc0efe74663eed 100644
--- a/src/SRHD/equation.jl
+++ b/src/SRHD/equation.jl
@@ -359,11 +359,11 @@ end
@with_signature function llf(equation::Equation2d)
- @accepts Prefix(lhs), 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, v
- @accepts Prefix(rhs), 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, v
+ @accepts Prefix(lhs), 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 Prefix(rhs), 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 nx, ny
- vmax = max(lhs_v, rhs_v)
+ vmax = max(lhs_max_v, rhs_max_v)
lhs_nflx = nx*lhs_flx_D_x + ny*lhs_flx_D_y
rhs_nflx = nx*rhs_flx_D_x + ny*rhs_flx_D_y
@@ -385,6 +385,34 @@ end
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 entropy(equation::Equation2d{EquationOfState.IdealGas})
@accepts rho, p, vx, vy, E, flx_E_x, flx_E_y
@unpack eos = equation
@@ -453,6 +481,26 @@ end
@returns nflx_D_x, nflx_D_y, nflx_Sx_x, nflx_Sx_y,
nflx_Sy_x, nflx_Sy_y, nflx_tau_x, nflx_tau_y
end
+@with_signature function bdry_ldg_nflux(eq::Equation2d)
+ @accepts Prefix(lhs), D, Sx, Sy, tau
+ @accepts Prefix(rhs), D, Sx, Sy, tau
+ @accepts nx, ny
+
+ nflx_D_x = nx * lhs_D
+ nflx_D_y = ny * lhs_D
+
+ nflx_Sx_x = nx * lhs_Sx
+ nflx_Sx_y = ny * lhs_Sx
+
+ nflx_Sy_x = nx * lhs_Sy
+ nflx_Sy_y = ny * lhs_Sy
+
+ nflx_tau_x = nx * lhs_tau
+ nflx_tau_y = ny * lhs_tau
+
+ @returns nflx_D_x, nflx_D_y, nflx_Sx_x, nflx_Sx_y,
+ nflx_Sy_x, nflx_Sy_y, nflx_tau_x, nflx_tau_y
+end
@with_signature function av_flux_2d(eq::Equation2d)
@@ -500,3 +548,77 @@ end
@returns nflx_D, nflx_Sx, nflx_Sy, nflx_tau
end
+@with_signature function bdry_av_nflux_2d(eq::Equation2d)
+ @accepts Prefix(lhs), ldg_D_x, ldg_D_y, ldg_Sx_x, ldg_Sx_y, ldg_Sy_x, ldg_Sy_y, ldg_tau_x, ldg_tau_y, smoothed_mu
+ @accepts Prefix(rhs), ldg_D_x, ldg_D_y, ldg_Sx_x, ldg_Sx_y, ldg_Sy_x, ldg_Sy_y, ldg_tau_x, ldg_tau_y, smoothed_mu
+ @accepts nflx_D, nflx_Sx, nflx_Sy, nflx_tau, nx, ny
+
+ lhs_nflx = nx*lhs_ldg_D_x + ny*lhs_ldg_D_y
+ nflx_D += lhs_smoothed_mu*lhs_nflx
+
+ lhs_nflx = nx*lhs_ldg_Sx_x + ny*lhs_ldg_Sx_y
+ nflx_Sx += lhs_smoothed_mu*lhs_nflx
+
+ lhs_nflx = nx*lhs_ldg_Sy_x + ny*lhs_ldg_Sy_y
+ nflx_Sy += lhs_smoothed_mu*lhs_nflx
+
+ lhs_nflx = nx*lhs_ldg_tau_x + ny*lhs_ldg_tau_y
+ nflx_tau += lhs_smoothed_mu*lhs_nflx
+
+ @returns nflx_D, nflx_Sx, nflx_Sy, nflx_tau
+end
+
+
+# @meshloop function av_nflux_2d(out, in1, in2, in3, mesh, eq::Equation2d)
+# @accepts Prefix(lhs), ldg_D_x, ldg_D_y, ldg_Sx_x, ldg_Sx_y, ldg_Sy_x, ldg_Sy_y, ldg_tau_x, ldg_tau_y, smoothed_mu = in1
+# @accepts Prefix(rhs), ldg_D_x, ldg_D_y, ldg_Sx_x, ldg_Sx_y, ldg_Sy_x, ldg_Sy_y, ldg_tau_x, ldg_tau_y, smoothed_mu = ini2
+# @accepts nflx_D, nflx_Sx, nflx_Sy, nflx_tau, nx, ny = in3
+#
+# @batch begin
+# lhs_nflx = nx*lhs_ldg_D_x + ny*lhs_ldg_D_y
+# rhs_nflx = nx*rhs_ldg_D_x + ny*rhs_ldg_D_y
+# nflx_D += 0.5 * (lhs_smoothed_mu*lhs_nflx + rhs_smoothed_mu*rhs_nflx)
+# end
+#
+# @batch begin
+# lhs_nflx = nx*lhs_ldg_Sx_x + ny*lhs_ldg_Sx_y
+# rhs_nflx = nx*rhs_ldg_Sx_x + ny*rhs_ldg_Sx_y
+# nflx_Sx += 0.5 * (lhs_smoothed_mu*lhs_nflx + rhs_smoothed_mu*rhs_nflx)
+# end
+#
+# @batch begin
+# lhs_nflx = nx*lhs_ldg_Sy_x + ny*lhs_ldg_Sy_y
+# rhs_nflx = nx*rhs_ldg_Sy_x + ny*rhs_ldg_Sy_y
+# nflx_Sy += 0.5 * (lhs_smoothed_mu*lhs_nflx + rhs_smoothed_mu*rhs_nflx)
+# end
+#
+# @batch begin
+# lhs_nflx = nx*lhs_ldg_tau_x + ny*lhs_ldg_tau_y
+# rhs_nflx = nx*rhs_ldg_tau_x + ny*rhs_ldg_tau_y
+# nflx_tau += 0.5 * (lhs_smoothed_mu*lhs_nflx + rhs_smoothed_mu*rhs_nflx)
+# end
+#
+# @returns nflx_D, nflx_Sx, nflx_Sy, nflx_tau = out
+# end
+#
+# ->
+#
+# # How to mix bulk and bdry variables using @turbo and indices?
+# function av_flux_2d(out, in1, in2, in3, mesh::Mesh2d, eq::Equation2d)
+# # extract variables ...
+# @turbo for i in 1:npoints(mesh)
+# # batch1 ...
+# end
+# @turbo for i in 1:npoints(mesh)
+# # batch2 ...
+# end
+# @turbo for i in 1:npoints(mesh)
+# # batch3 ...
+# end
+# @turbo for i in 1:npoints(mesh)
+# # batch4 ...
+# end
+# @turbo for i in 1:npoints(mesh)
+# # batch5 ...
+# end
+# end
diff --git a/src/SRHD/initialdata.jl b/src/SRHD/initialdata.jl
index 40caf99f7291d5747d8aa48fe3157775ab425c84..38ee023992b4ee87c48972447c6c27a8d4cc14b4 100644
--- a/src/SRHD/initialdata.jl
+++ b/src/SRHD/initialdata.jl
@@ -11,6 +11,7 @@ end
function initialdata!(env, P::Project2d, prms, mesh::Mesh2d)
initialdata_equation_of_state!(env, P, P.equation.eos, prms, mesh)
initialdata_hrsc!(env, P, P.hrsc, prms, mesh)
+ initialdata_bdrycond!(env, P, P.bdrycond, mesh)
end
@@ -197,7 +198,7 @@ function initialdata_equation_of_state!(env, P, eos::EquationOfState.IdealGas, p
p0 = @. 1.0 * z^0
Ï0, vx0, vy0, p0
- elseif id == "id_shocktube_2d"
+ elseif id == "shocktube_2d"
@toggled_assert xmin == -1.0 && xmax == 1.0 && ymin == -1.0 && ymax == 1.0
# TODO What EoS is this test for?
Ï0, vx0, vy0, p0 = [ similar(x) for _ = 1:4 ]
@@ -310,3 +311,39 @@ function initialdata_tci!(env, P::Project, tci::TCI.AbstractTCI, prms, mesh)
@unpack flag = get_cell_variables(env.cache)
@. flag = 0.0
end
+
+
+#######################################################################
+# Bdryconds #
+#######################################################################
+
+
+function initialdata_bdrycond!(env, P::Project2d, bdrycond, mesh)
+ TODO()
+end
+
+
+function initialdata_bdrycond!(env, P::Project2d, bdrycond::dg1d.DirichletBC2, mesh)
+ @unpack D, Sx, Sy, tau = get_dynamic_variables(mesh.cache)
+ @unpack flx_D_x, flx_D_y, flx_Sx_x, flx_Sx_y, flx_Sy_x, flx_Sy_y,
+ flx_tau_x, flx_tau_y, max_v = get_static_variables(mesh.cache)
+ @unpack init_flx_D_x, init_flx_D_y, init_flx_Sx_x, init_flx_Sx_y, init_flx_Sy_x, init_flx_Sy_y,
+ init_flx_tau_x, init_flx_tau_y, init_max_v = get_static_variables(mesh.cache)
+ @unpack init_D, init_Sx, init_Sy, init_tau = get_static_variables(mesh.cache)
+ broadcast_volume_2!(cons2prim, P.equation, mesh)
+ broadcast_volume_2!(maxspeed, P.equation, mesh)
+ broadcast_volume_2!(flux, P.equation, mesh)
+ init_D .= D
+ init_Sx .= Sx
+ init_Sy .= Sy
+ init_tau .= tau
+ init_flx_D_x .= flx_D_x
+ init_flx_D_y .= flx_D_y
+ init_flx_Sx_x .= flx_Sx_x
+ init_flx_Sx_y .= flx_Sx_y
+ init_flx_Sy_x .= flx_Sy_x
+ init_flx_Sy_y .= flx_Sy_y
+ init_flx_tau_x .= flx_tau_x
+ init_flx_tau_y .= flx_tau_y
+ init_max_v .= max_v
+end
diff --git a/src/SRHD/plots.jl b/src/SRHD/plots.jl
index de7922ad78d92e20a51adff078a10ec35496aa93..718dcaf0feb90b6a65a086da1a4997f489363f02 100644
--- a/src/SRHD/plots.jl
+++ b/src/SRHD/plots.jl
@@ -230,6 +230,6 @@ function plot_rhs(env, P)
end
plt = Plots.plot(plts..., layout=(2,2))
-
+
return plt
end
diff --git a/src/SRHD/rhs.jl b/src/SRHD/rhs.jl
index 9dd7f5fb8b1eb7c1a6eb67f012c33c2a783bd9a0..1e30e996a9af1b029861f3eac6d93952e72edbcc 100644
--- a/src/SRHD/rhs.jl
+++ b/src/SRHD/rhs.jl
@@ -69,10 +69,13 @@ function timestep(env, P, hrsc::Maybe{HRSC.AbstractHRSC}, mesh::Mesh2d)
broadcast_volume_2!(cons2prim, equation, mesh)
broadcast_volume_2!(maxspeed, equation, mesh)
- vmax = dg1d.absolute_maximum(max_v)
+ # vmax = dg1d.absolute_maximum(max_v)
+ i_vmax, vmax = dg1d.absolute_maximum_index(max_v)
vmax_limit = 0.9999
if vmax > vmax_limit
- @warn "Limiting timestep due to maximum speed exceeding $vmax_limit"
+ i_cell = find_cell_index(i_vmax, mesh)
+ cell = mesh.tree.cells[i_cell]
+ @warn "Limiting timestep due to maximum speed exceeding $vmax_limit in cell '$cell'"
vmax = vmax_limit
end
@@ -202,8 +205,11 @@ end
# 2d #
#######################################################################
+function dg1d.rhs!(env::Environment, P::Project2d, bdryconds)
+ dg1d.rhs!(env, P, P.hrsc, env.mesh)
+end
-function rhs!(env, P, ::Mesh2d)
+function dg1d.rhs!(env, P, ::Nothing, ::Mesh2d)
@unpack cache, mesh = env
@unpack equation = P
@@ -229,8 +235,7 @@ function rhs!(env, P, ::Mesh2d)
end
-function rhs!(env, P::Project2d, hrsc::HRSC.AbstractArtificialViscosity, mesh::Mesh2d,
- bdryconds, ldg_bdryconds, av_bdryconds)
+function dg1d.rhs!(env, P::Project2d, hrsc::HRSC.AbstractArtificialViscosity, mesh::Mesh2d)
@unpack equation = P
@unpack D, Sx, Sy, tau = get_dynamic_variables(mesh.cache)
@@ -263,6 +268,8 @@ function rhs!(env, P::Project2d, hrsc::HRSC.AbstractArtificialViscosity, mesh::M
## solve auxiliary equation: q + ∂x u = 0
broadcast_faces_2!(ldg_nflux, equation, mesh)
+ dg1d.broadcast_bdry_2!(bdry_ldg_nflux, equation, P.bdrycond, mesh)
+
# broadcast_boundaryconditions!(central_flux, ldg_bdryconds, cache, mesh, 0.0)
compute_rhs_weak_form!(ldg_D_x, D, 0, nflx_D_x, mesh)
compute_rhs_weak_form!(ldg_D_y, 0, D, nflx_D_y, mesh)
@@ -276,11 +283,14 @@ function rhs!(env, P::Project2d, hrsc::HRSC.AbstractArtificialViscosity, mesh::M
## compute rhs of regularized equation: ∂t u + ∂x f + ∂x μ q = 0
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)
diff --git a/src/SRHD/setup.jl b/src/SRHD/setup.jl
index 488f0bc394662c09885683154ecf426ac453fde4..6f4e02f12f70ac4ff6d356aa01e29026ff0b0fda 100644
--- a/src/SRHD/setup.jl
+++ b/src/SRHD/setup.jl
@@ -61,8 +61,9 @@ function Project(env::Environment, prms, mesh::Mesh2d)
eos = EquationOfState.make_EquationOfState(prms["EquationOfState"])
hrsc = HRSC.make_HRSC(mesh, prms["HRSC"])
equation = make_Equation(eos, prms["SRHD"], mesh)
+ bdrycond = dg1d.DirichletBC2()
- P = Project2d(equation, hrsc)
+ P = Project2d(equation, hrsc, bdrycond)
# register variables
# TODO add a ::Nothing overload for register_variables!
@@ -70,6 +71,7 @@ function Project(env::Environment, prms, mesh::Mesh2d)
@unpack cache = env
_register_variables!(mesh, P)
_register_variables!(mesh, hrsc)
+ _register_variables!(mesh, bdrycond)
display(cache)
# TODO Bdry conditions should be moved to project and should only be just a 'lightweight'
@@ -191,3 +193,15 @@ function _register_variables!(mesh, P::Project2d)
cell_variablenames = (:cellmax_v,),
)
end
+
+
+function _register_variables!(mesh, bdrycond::dg1d.DirichletBC2)
+ register_variables!(mesh.cache,
+ static_variablenames = (:init_D, :init_Sx, :init_Sy, :init_tau,
+ :init_flx_D_x, :init_flx_D_y,
+ :init_flx_Sx_x, :init_flx_Sx_y,
+ :init_flx_Sy_x, :init_flx_Sy_y,
+ :init_flx_tau_x, :init_flx_tau_y,
+ :init_max_v)
+ )
+end
diff --git a/src/SRHD/types.jl b/src/SRHD/types.jl
index a8d0f404acbd93b02146aff49cb027c38b14ab73..2ca4d228d5edb1897e5368289ad098527db1a3da 100644
--- a/src/SRHD/types.jl
+++ b/src/SRHD/types.jl
@@ -36,7 +36,9 @@ end
struct Project2d{Eq<:Equation2d,
- T_HRSC<:Maybe{HRSC.AbstractHRSC}}
+ T_HRSC<:Maybe{HRSC.AbstractHRSC},
+ T_BC<:dg1d.AbstractBC}
equation::Eq
hrsc::T_HRSC
+ bdrycond::T_BC
end
diff --git a/src/bdryconditions.jl b/src/bdryconditions.jl
index fbbb17d5395a74563f7bd9ea60c57e467e56b80a..da58abbe06d6c93c428492a6f062aa59349afd18 100644
--- a/src/bdryconditions.jl
+++ b/src/bdryconditions.jl
@@ -16,12 +16,15 @@ abstract type AbstractBC end
Dirichlet boundary condition, `bdryval` might be time dependent.
"""
-struct DirichletBC <: AbstractBC
+struct DirichletBC <: AbstractBC
bdryval::Function
end
(bc::DirichletBC)(t) = bc.bdryval(t)
+struct DirichletBC2<:AbstractBC end
+
+
"""
Constructor for Dirichlet BC with constant boundary value.
"""
diff --git a/src/mesh.jl b/src/mesh.jl
index c9d0f660b66199ee8372c4807be22a429e64045a..d40b9d8710a18064b34ebe0d43d6748d91660878 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -84,7 +84,7 @@ function Mesh2d(; N=[5,5], K=[4,4], range=[-1.0,1.0, -1.0,1.0],
# TODO Remove this after updating parameter names in SpectralElement
quadrature_method = Dict("lgl" => :LGL, "glgl" => :GLGL, "lgl_lumped" => :LGL_lumped)[basis]
- tree = Tree2d(Kx,Ky,periodic=(true,true))
+ tree = Tree2d(Kx,Ky,periodic=(false,false))
boxes = place_boxes(tree, Float64.(xrange), Float64.(yrange))
element_x = SpectralElement(Nx, Symbol(quadrature_method))
element_y = SpectralElement(Ny, Symbol(quadrature_method))
@@ -514,7 +514,7 @@ end
return rv, state+1
end
-eachcell(mesh::AbstractMesh) = CellIterator(mesh)
+# eachcell(mesh::AbstractMesh) = CellIterator(mesh)
eachcell(mesh::AbstractMesh, data) = CellDataIterator(mesh, data)
diff --git a/src/with_signature.jl b/src/with_signature.jl
index 8fbe5c5c381bf5cab378522e77273717e7ecf880..35193f5d64d74b62b90971ab1daa41b54aa01b61 100644
--- a/src/with_signature.jl
+++ b/src/with_signature.jl
@@ -533,49 +533,125 @@ function _broadcast_face_2!(f, mesh::Mesh2d)
Kx, Ky = mesh.tree.dims
Nx, Ny = Npts, Npts
- @unpack nx, ny = get_bdry_variables(mesh.cache)
+ # @unpack nx, ny = get_bdry_variables(mesh.cache)
- for cell = 1:ncells(mesh)
+ for (icell,cell) = enumerate(mesh.tree.cells)
- bulk_in = cellindices(mesh, cell)
- rng_lhs, rng_rhs, rng_down, rng_up = faceindices(mesh, cell)
+ bulk_in = cellindices(mesh, icell)
+ rng_lhs, rng_rhs, rng_down, rng_up = faceindices(mesh, icell)
# xmin face
- cell_out = mesh.tree.cells[cell].neighbors[Cart2d.Xmin]
- bulk_out = cellindices(mesh, cell_out)
- for (fidx, face_idx) in enumerate(rng_lhs)
- idx_out = bulk_out[end,fidx]
- idx_in = bulk_in[1,fidx]
- @inline f(face_idx, idx_in, idx_out, face_idx)
+ if has_neighbor(cell, Cart2d.Xmin)
+ cell_out = mesh.tree.cells[icell].neighbors[Cart2d.Xmin]
+ bulk_out = cellindices(mesh, cell_out)
+ for (fidx, face_idx) in enumerate(rng_lhs)
+ idx_out = bulk_out[end,fidx]
+ idx_in = bulk_in[1,fidx]
+ @inline f(face_idx, idx_in, idx_out, face_idx)
+ end
end
# xmax face
- cell_out = mesh.tree.cells[cell].neighbors[Cart2d.Xmax]
- bulk_out = cellindices(mesh, cell_out)
- for (fidx, face_idx) in enumerate(rng_rhs)
- idx_out = bulk_out[1,fidx]
- idx_in = bulk_in[end,fidx]
- @inline f(face_idx, idx_in, idx_out, face_idx)
+ if has_neighbor(cell, Cart2d.Xmax)
+ cell_out = mesh.tree.cells[icell].neighbors[Cart2d.Xmax]
+ bulk_out = cellindices(mesh, cell_out)
+ for (fidx, face_idx) in enumerate(rng_rhs)
+ idx_out = bulk_out[1,fidx]
+ idx_in = bulk_in[end,fidx]
+ @inline f(face_idx, idx_in, idx_out, face_idx)
+ end
end
# ymin face
- cell_out = mesh.tree.cells[cell].neighbors[Cart2d.Ymin]
- bulk_out = cellindices(mesh, cell_out)
- for (fidx, face_idx) in enumerate(rng_down)
- idx_out = bulk_out[fidx,end]
- idx_in = bulk_in[fidx,1]
- @inline f(face_idx, idx_in, idx_out, face_idx)
+ if has_neighbor(cell, Cart2d.Ymin)
+ cell_out = mesh.tree.cells[icell].neighbors[Cart2d.Ymin]
+ bulk_out = cellindices(mesh, cell_out)
+ for (fidx, face_idx) in enumerate(rng_down)
+ idx_out = bulk_out[fidx,end]
+ idx_in = bulk_in[fidx,1]
+ @inline f(face_idx, idx_in, idx_out, face_idx)
+ end
end
# ymax face
- cell_out = mesh.tree.cells[cell].neighbors[Cart2d.Ymax]
- bulk_out = cellindices(mesh, cell_out)
- for (fidx, face_idx) in enumerate(rng_up)
- idx_out = bulk_out[fidx,1]
- idx_in = bulk_in[fidx,end]
- @inline f(face_idx, idx_in, idx_out, face_idx)
+ if has_neighbor(cell, Cart2d.Ymax)
+ cell_out = mesh.tree.cells[icell].neighbors[Cart2d.Ymax]
+ bulk_out = cellindices(mesh, cell_out)
+ for (fidx, face_idx) in enumerate(rng_up)
+ idx_out = bulk_out[fidx,1]
+ idx_in = bulk_in[fidx,end]
+ @inline f(face_idx, idx_in, idx_out, face_idx)
+ end
+ end
+
+ end
+ return
+end
+
+
+function broadcast_bdry_2!(f, equation, bdrycond::AbstractBC, mesh::Mesh)
+ @toggled_assert has_signature(f, equation) "broadcast_bdry!: '$f' was not defined using @with_signature"
+ vars1, vars2, vars3 = accepts(f, equation)
+ returns_vars = returns(f, equation)
+ args1, args2, args3 = wrap_variables_in_soa(mesh.cache, vars1, vars2, vars3)
+ rets = wrap_variables_in_soa(mesh.cache, returns_vars)
+ wrapped_f = (i,j1,j2,j3) -> rets[i] = f(args1[j1],args2[j2],args3[j3],equation)
+ _broadcast_bdry_2!(wrapped_f, bdrycond, mesh)
+ return
+end
+
+
+function _broadcast_bdry_2!(f, bdrycond::DirichletBC2, mesh::Mesh2d)
+
+ @unpack Npts = mesh.element
+ Kx, Ky = mesh.tree.dims
+ Nx, Ny = Npts, Npts
+
+ # @unpack nx, ny = get_bdry_variables(mesh.cache)
+
+ for bdrycell in mesh.tree.bdrycells
+
+ rng_lhs, rng_rhs, rng_down, rng_up = faceindices(mesh, bdrycell.index)
+ idxs = cellindices(mesh, bdrycell.index)
+
+ # @show rng_lhs, rng_rhs, rng_down, rng_up
+ # TODO("need to use bulk and face indices here!")
+
+ if !has_neighbor(bdrycell, Cart2d.Xmin)
+
+ for (fidx, face_idx) in enumerate(rng_lhs)
+ bidx = idxs[fidx,1]
+ @inline f(face_idx, bidx, bidx, bidx)
+ end
+
+ elseif !has_neighbor(bdrycell, Cart2d.Xmax)
+
+ for (fidx, face_idx) in enumerate(rng_rhs)
+ bidx = idxs[fidx,end]
+ @inline f(face_idx, bidx, bidx, bidx)
+ end
+
+ elseif !has_neighbor(bdrycell, Cart2d.Ymin)
+
+ for (fidx, face_idx) in enumerate(rng_down)
+ bidx = idxs[1,fidx]
+ @inline f(face_idx, bidx, bidx, bidx)
+ end
+
+ elseif !has_neighbor(bdrycell, Cart2d.Ymax)
+
+ for (fidx, face_idx) in enumerate(rng_up)
+ bidx = idxs[end,fidx]
+ @inline f(face_idx, bidx, bidx, bidx)
+ end
+
+ else
+
+ error("Something went wrong ...")
+
end
end
+ # TODO()
return
end
diff --git a/test/IntegrationTests/refs/srhd2d_shocktube_entropyav.toml b/test/IntegrationTests/refs/srhd2d_shocktube_entropyav.toml
new file mode 100644
index 0000000000000000000000000000000000000000..fe523d2fbb2ca1d542a612168f9b86d75c7c8d1e
--- /dev/null
+++ b/test/IntegrationTests/refs/srhd2d_shocktube_entropyav.toml
@@ -0,0 +1,31 @@
+[EquationOfState]
+eos = "idealgas"
+idealgas_gamma = 1.6666666666666666666667
+
+[SRHD]
+id = "shocktube_2d"
+
+[Mesh]
+range = [ -1.0,1.0, -1.0,1.0 ]
+cfl = 0.01
+n = [1,1]
+k = [88,88]
+basis = "lgl"
+
+[Output]
+every_iteration = 1
+every_dt = 0.05
+variables = [ "D", "Sx", "Sy", "tau", "smoothed_mu" ]
+enable1d = false
+enable2d = true
+format2d = "vtk"
+
+[Log]
+progress_stdout = false
+
+[Evolution]
+tend = 4.4
+
+[HRSC]
+method = "av"
+av_method = "entropy"