diff --git a/src/EulerEq/rhs.jl b/src/EulerEq/rhs.jl
index 2d78894584c1da9cec7a6b89709dca14098cfde0..f8891df8245ad9eb121329d1c1bcfc10bac7e60b 100644
--- a/src/EulerEq/rhs.jl
+++ b/src/EulerEq/rhs.jl
@@ -203,7 +203,7 @@ function step_fv_muscl!(P::Project, mesh::Mesh1d{FVElement}, state)
# limit slope and compute left (L) and (R) interface values
for k in 1:K
- dvL = limit_slope(dvâ‚Š[k], dvâ‚‹[k], dl, prms_slope_limiter)
+ dvL = dg1d.limit_slope(dvâ‚Š[k], dvâ‚‹[k], dl, prms_slope_limiter)
vL[k] = v[k] - dvL/2
vR[k] = v[k] + dvL/2
end
@@ -317,28 +317,6 @@ function step_fv_muscl!(P::Project, mesh::Mesh1d{FVElement}, state)
end
-@inline function limit_slope(a, b, dx, prms)
- if prms.method === :none
- return 0.0
- elseif prms.method === :minmod
- return TCI.minmod(a, b)
- elseif prms.method === :muscl
- return TCI.minmod((a+b)/2, 2*a, 2*b)
- elseif prms.method === :superbee
- return TCI.minmod(TCI.maxmod(a,b),TCI.minmod(2*a,2*b))
- elseif prms.method === :van_albada
- c = dx^3
- return TCI.minmod( ((a^2+c^2)*b+(b^2+c^2)*a)/(a^2+b^2+c^2), 2*a, 2*b )
- elseif prms.method === :van_leer
- return TCI.minmod(2*a*b/(a+b), 2*a, 2*b)
- elseif prms.method === :tvb
- return TCI.minmod_M(dx, prms.tvb_M, (a+b)/2, 2*a, 2*b)
- else
- TODO(prms.method)
- end
-end
-
-
#######################################################################
# DG rhs #
#######################################################################
diff --git a/src/GRHD/rhs.jl b/src/GRHD/rhs.jl
index 09187952e0de799664d0e55680354dab50626cad..6a1b86368d6ddb77521157338e0fc236ea76ca01 100644
--- a/src/GRHD/rhs.jl
+++ b/src/GRHD/rhs.jl
@@ -452,7 +452,7 @@ function step_fv_slope_limiter!(mesh, P::Project{:spherical1d}, state)
kâ‚‹, kâ‚Š = k-1, k+1
wâ‚‹, wâ‚€, wâ‚Š = w[kâ‚‹], w[k], w[kâ‚Š]
dwâ‚Š[k], dwâ‚‹[k] = wâ‚Š-wâ‚€, wâ‚€-wâ‚‹
- dwL = limit_slope(dwâ‚Š[k], dwâ‚‹[k], dl, prms_slope_limiter)
+ dwL = dg1d.limit_slope(dwâ‚Š[k], dwâ‚‹[k], dl, prms_slope_limiter)
wL[k] = w[k] - dwL/2
wR[k] = w[k] + dwL/2
end
@@ -725,7 +725,7 @@ function step_fv_slope_limiter!(mesh, P::Project{:doublecartoon}, state)
kâ‚‹, kâ‚Š = k-1, k+1
wâ‚‹, wâ‚€, wâ‚Š = w[kâ‚‹], w[k], w[kâ‚Š]
dwâ‚Š[k], dwâ‚‹[k] = wâ‚Š-wâ‚€, wâ‚€-wâ‚‹
- dwL = limit_slope(dwâ‚Š[k], dwâ‚‹[k], dl, prms_slope_limiter)
+ dwL = dg1d.limit_slope(dwâ‚Š[k], dwâ‚‹[k], dl, prms_slope_limiter)
wL[k] = w[k] - dwL/2
wR[k] = w[k] + dwL/2
end
@@ -1129,30 +1129,6 @@ end
end
-@inline function limit_slope(a, b, dx, prms)
- if prms.method === :none
- return 0.0
- elseif prms.method === :minmod
- return TCI.minmod(a, b)
- elseif prms.method === :muscl
- return TCI.minmod((a+b)/2, 2*a, 2*b)
- elseif prms.method === :superbee
- return TCI.minmod(TCI.maxmod(a,b),TCI.minmod(2*a,2*b))
- elseif prms.method === :van_albada
- c = dx^3
- return TCI.minmod( ((a^2+c^2)*b+(b^2+c^2)*a)/(a^2+b^2+c^2), 2*a, 2*b )
- elseif prms.method === :van_leer
- return TCI.minmod(2*a*b/(a+b), 2*a, 2*b)
- elseif prms.method === :tvb
- return TCI.minmod_M(dx, prms.tvb_M, (a+b)/2, 2*a, 2*b)
- elseif prms.method === :mc
- return dg1d.MC(b, a)
- else
- TODO(prms.method)
- end
-end
-
-
#######################################################################
# DG RHS #
#######################################################################
diff --git a/src/HRSC/mda.jl b/src/HRSC/mda.jl
index a85e51f24205b331c4fc10df19d2be9f1b346ba2..18779814d62b32c9ec48a2c7489534f75732ed7c 100644
--- a/src/HRSC/mda.jl
+++ b/src/HRSC/mda.jl
@@ -88,7 +88,7 @@ function compute_viscosity!(
@. data_y = log(skylined_mod_uh)
@. data_y = max(data_y, -20.0)
_, _, tau = dg1d.least_sqr_linear(data_x, data_y)
- isnan(tau) && error("undesired nan encountered")
+ !isfinite(tau) && error("undesired nan encountered")
mu[j] *= v[j]
if tau > n_high
diff --git a/src/SRHD/rhs.jl b/src/SRHD/rhs.jl
index 271cf2ff5765e0d6662cd3fb97e5f43de830731a..58ac2239d3762d7f13ff200d1f19cee62fdbbc6b 100644
--- a/src/SRHD/rhs.jl
+++ b/src/SRHD/rhs.jl
@@ -221,7 +221,7 @@ function step_fv_muscl!(P::Project, mesh::Mesh1d{FVElement}, state)
# limit slope and compute left (L) and (R) interface values
for k in 1:K
- dwL = limit_slope(dwâ‚Š[k], dwâ‚‹[k], dl, prms_slope_limiter)
+ dwL = dg1d.limit_slope(dwâ‚Š[k], dwâ‚‹[k], dl, prms_slope_limiter)
wL[k] = w[k] - dwL/2
wR[k] = w[k] + dwL/2
end
@@ -334,28 +334,6 @@ function step_fv_muscl!(P::Project, mesh::Mesh1d{FVElement}, state)
end
-@inline function limit_slope(a, b, dx, prms)
- if prms.method === :none
- return 0.0
- elseif prms.method === :minmod
- return TCI.minmod(a, b)
- elseif prms.method === :muscl
- return TCI.minmod((a+b)/2, 2*a, 2*b)
- elseif prms.method === :superbee
- return TCI.minmod(TCI.maxmod(a,b),TCI.minmod(2*a,2*b))
- elseif prms.method === :van_albada
- c = dx^3
- return TCI.minmod( ((a^2+c^2)*b+(b^2+c^2)*a)/(a^2+b^2+c^2), 2*a, 2*b )
- elseif prms.method === :van_leer
- return TCI.minmod(2*a*b/(a+b), 2*a, 2*b)
- elseif prms.method === :tvb
- return TCI.minmod_M(dx, prms.tvb_M, (a+b)/2, 2*a, 2*b)
- else
- TODO(prms.method)
- end
-end
-
-
#######################################################################
# DG rhs #
#######################################################################
diff --git a/src/ScalarEq/callbacks.jl b/src/ScalarEq/callbacks.jl
index 49af083706b0bc070390b166eb5130e62800148a..48a1f43d5e5c923c69f2d992b29284d3c04f08e4 100644
--- a/src/ScalarEq/callbacks.jl
+++ b/src/ScalarEq/callbacks.jl
@@ -246,7 +246,6 @@ end
#######################################################################
-# TODO
callback_analyze(state_u, state_t, env, P, mesh::Mesh) = nothing
@@ -255,6 +254,8 @@ function callback_analyze(state_u, state_t, env, P, mesh::Mesh1d)
@unpack u_norm2 = get_global_variables(mesh)
u_norm2 .= dg1d.norm_L2(u, mesh)^2
end
+# TODO Make this use dg1d.limit_slope!
+# after moving src/ScalarEq/rhs.jl:limit_slope! to src/limters.jl
function callback_analyze(state_u, state_t, env, P, mesh::Mesh1d{<:FVElement})
@unpack u = get_dynamic_variables(mesh)
@unpack u_norm2 = get_global_variables(mesh)
diff --git a/src/ScalarEq/rhs.jl b/src/ScalarEq/rhs.jl
index 9cd9b0e5b2e4c512451566262cad9f68f9442d42..e1f3cc61b6b9b2ab8265a5bc55d1a6bb10319a96 100644
--- a/src/ScalarEq/rhs.jl
+++ b/src/ScalarEq/rhs.jl
@@ -132,6 +132,7 @@ function rhs_fv_muscl!(env, mesh::Mesh1d{FVElement}, P::Project)
end
+# TODO Move this to src/limiters.jl
function limit_slope!(P::Project, mesh::Mesh1d{FVElement}, u)
@unpack equation = P
@unpack uL, uR,
@@ -158,7 +159,7 @@ function limit_slope!(P::Project, mesh::Mesh1d{FVElement}, u)
# limit slope and compute left (L) and (R) interface values
for k in 1:K
- duL = limit_slope(duâ‚Š[k], duâ‚‹[k], dl, prms_slope_limiter)
+ duL = dg1d.limit_slope(duâ‚Š[k], duâ‚‹[k], dl, prms_slope_limiter)
uL[k] = u[k] - duL/2
uR[k] = u[k] + duL/2
end
@@ -240,28 +241,6 @@ function step_fv_muscl!(P::Project, mesh::Mesh1d{FVElement}, u)
end
-@inline function limit_slope(a, b, dx, prms)
- if prms.method === :none
- return 0.0
- elseif prms.method === :minmod
- return TCI.minmod(a, b)
- elseif prms.method === :muscl
- return TCI.minmod((a+b)/2, 2*a, 2*b)
- elseif prms.method === :superbee
- return TCI.minmod(TCI.maxmod(a,b),TCI.minmod(2*a,2*b))
- elseif prms.method === :van_albada
- c = dx^3
- return TCI.minmod( ((a^2+c^2)*b+(b^2+c^2)*a)/(a^2+b^2+c^2), 2*a, 2*b )
- elseif prms.method === :van_leer
- return TCI.minmod(2*a*b/(a+b), 2*a, 2*b)
- elseif prms.method === :tvb
- return TCI.minmod_M(dx, prms.tvb_M, (a+b)/2, 2*a, 2*b)
- else
- TODO(prms.method)
- end
-end
-
-
#######################################################################
# DG RHS #
#######################################################################
@@ -269,14 +248,14 @@ end
function rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
if mesh.element.kind ∈ (:modal_bspline2, :modal_bspline1)
- modal_rhs!(env, mesh, P, hrsc)
+ modal_rhs!(env, mesh, P)
else
- nodal_rhs!(env, mesh, P, hrsc)
+ nodal_rhs!(env, mesh, P)
end
end
-function nodal_rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
+function nodal_rhs!(env, mesh::Mesh1d, P::Project)
@unpack flx_u, src_u = get_static_variables(mesh)
@unpack u = get_dynamic_variables(mesh)
@unpack rhs_u = get_rhs_variables(mesh)
@@ -291,7 +270,7 @@ function nodal_rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
end
-function modal_rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
+function modal_rhs!(env, mesh::Mesh1d, P::Project)
@unpack flx_u, src_u,
u_modal, flx_u_modal,
src_u_modal = get_static_variables(mesh)
@@ -299,6 +278,11 @@ function modal_rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
@unpack rhs_u = get_rhs_variables(mesh)
@unpack nflx_u, bdry_u = get_bdry_variables(mesh)
+ if P.prms.hrsc === :slope_limiter
+ @unpack flag = get_cell_variables(mesh)
+ dg1d.limit_slope_convex_hull!(u, flag, P.prms.workspace, mesh)
+ end
+
u_modal .= u
for (v_u, v_u_modal) in zip(eachcell(mesh,u),eachcell(mesh,u_modal))
mul!(v_u, P.prms.V_bspline, v_u_modal)
@@ -315,6 +299,7 @@ function modal_rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
mul!(v_src_u_modal, P.prms.invV_bspline, v_src_u)
end
compute_rhs_weak_form!(rhs_u, flx_u_modal, src_u_modal, nflx_u, mesh)
+
u .= u_modal
return
diff --git a/src/ScalarEq/setup.jl b/src/ScalarEq/setup.jl
index 727b4242262ea5d758db63360ae1e5230e1031c5..0d1f7926d6b7d25fa67487c9bcdf4fe1e1a78d8c 100644
--- a/src/ScalarEq/setup.jl
+++ b/src/ScalarEq/setup.jl
@@ -24,10 +24,11 @@ function Project(env::Environment, mesh::Mesh1d, prms)
else
V_bspline, invV_bspline = nothing, nothing
end
+ workspace = dg1d.Workspace()
fixedprms = (; bernstein, av_derivative_scheme, av_drag, hrsc, muscl_omega,
slope_limiter_method, slope_limiter_tvb_M,
av_recompute_substeps,
- V_bernstein, V_bspline, invV_bspline,
+ V_bernstein, V_bspline, invV_bspline, workspace,
analyze_error, analyze_error_norm)
# construct Project
diff --git a/src/TCI/minmod.jl b/src/TCI/minmod.jl
index c87c0a0113b85445ee63d8f7bf764873baf1228c..2efdea49619b5206ac43ac9be374dc1e6aa75acf 100644
--- a/src/TCI/minmod.jl
+++ b/src/TCI/minmod.jl
@@ -11,8 +11,6 @@ function compute_indicator!(
tci::Minmod)
@unpack mesh = tci
- @unpack invdetJ = get_static_variables(mesh.cache)
-
Npts, K = layout(mesh)
L, = dg1d.widths(mesh)
dL = L/K
@@ -23,7 +21,7 @@ function compute_indicator!(
@unpack threshold = tci
# fill flag by applying minmod criterion
- for k in 1:K # loop over cells
+ for k in 1:K
# assuming periodic grid
k_lhs, k_rhs = periodic_index(k-1, K), periodic_index(k+1, K)
uavg, uavg_lhs, uavg_rhs = u_averages[k], u_averages[k_lhs], u_averages[k_rhs]
diff --git a/src/TCI/modaldecayaverage.jl b/src/TCI/modaldecayaverage.jl
index 12b20c653a35555b0c1e4c9c344e63365b40b926..e8bacc4ed8a81763b80f93cfd95bf2f610c4e09f 100644
--- a/src/TCI/modaldecayaverage.jl
+++ b/src/TCI/modaldecayaverage.jl
@@ -81,7 +81,7 @@ function compute_indicator!(
@views for m in 1:N
uh_skylined = maximum(uh_mod[min(m,N-1):end])
log_uh = log(uh_skylined)
- lsq_y[m] = isinf(log_uh) ? - 20.0 : log_uh
+ lsq_y[m] = !isfinite(log_uh) ? - 20.0 : log_uh
end
# estimate effective decay rate by least square fit of modal coefficients, (3.9)
diff --git a/src/bspline.jl b/src/bspline.jl
index 3919fd54d332437a2c6b3404d004c79358246df6..fe51f5678d9da134e05631745aa8e8199ddd8fac 100644
--- a/src/bspline.jl
+++ b/src/bspline.jl
@@ -352,7 +352,7 @@ end
function integrate(f::AbstractVector, bs::Bspline1)
@toggled_assert length(f) == bs.Npts
sum = 0.0
- @inbounds for i in bs.Npts
+ @inbounds for i in 1:bs.Npts
fi, xl, xr = f[i], bs.xs[i], bs.xs[i+2]
sum += fi * (xr-xl)/2
end
@@ -361,7 +361,7 @@ end
function integrate(f::Function, idxs::AbstractVector{<:Integer}, bs::Bspline1)
@toggled_assert length(idxs) == bs.Npts
sum = 0.0
- @inbounds for i in bs.Npts
+ @inbounds for i in 1:bs.Npts
fi, xl, xr = f(idxs[i]), bs.xs[i], bs.xs[i+2]
sum += fi * (xr-xl)/2
end
@@ -374,7 +374,7 @@ end
function integrate(f::AbstractVector, bs::Bspline2)
@toggled_assert length(f) == bs.Npts+1
sum = 0.0
- @inbounds for i in bs.Npts+1
+ @inbounds for i in 1:bs.Npts+1
fi, xl, xr = f[i], bs.xs[i], bs.xs[i+3]
sum += fi * (xr-xl)/3
end
@@ -383,7 +383,7 @@ end
function integrate(f::Function, idxs::AbstractVector{<:Integer}, bs::Bspline2)
@toggled_assert length(idxs) == bs.Npts+1
sum = 0.0
- @inbounds for i in bs.Npts+1
+ @inbounds for i in 1:bs.Npts+1
fi, xl, xr = f(idxs[i]), bs.xs[i], bs.xs[i+3]
sum += fi * (xr-xl)/3
end
diff --git a/src/dg1d.jl b/src/dg1d.jl
index 48be5debad46941b2bbe6866f9dae41d58828395..d52d1f48f9965a21c6414e569cf4e22d5f370d94 100644
--- a/src/dg1d.jl
+++ b/src/dg1d.jl
@@ -88,6 +88,8 @@ export Mesh, Mesh1d, Mesh2d, MeshInterpolator, layout, n_points, n_cells, dimens
min_grid_spacing
include("mesh.jl")
+include("limiters.jl")
+
export AbstractEquation
include("equations.jl")
diff --git a/src/limiters.jl b/src/limiters.jl
new file mode 100644
index 0000000000000000000000000000000000000000..c6c388ec659c02efc73af3e375b87b989fb842fc
--- /dev/null
+++ b/src/limiters.jl
@@ -0,0 +1,222 @@
+@inline function limit_slope(a, b, dx, prms)
+ if prms.method === :none
+ return 0.0
+ elseif prms.method === :minmod
+ return TCI.minmod(a, b)
+ elseif prms.method === :muscl
+ return TCI.minmod((a+b)/2, 2*a, 2*b)
+ elseif prms.method === :superbee
+ return TCI.minmod(TCI.maxmod(a,b),TCI.minmod(2*a,2*b))
+ elseif prms.method === :van_albada
+ c = dx^3
+ return TCI.minmod( ((a^2+c^2)*b+(b^2+c^2)*a)/(a^2+b^2+c^2), 2*a, 2*b )
+ elseif prms.method === :van_leer
+ return TCI.minmod(2*a*b/(a+b), 2*a, 2*b)
+ elseif prms.method === :tvb
+ return TCI.minmod_M(dx, prms.tvb_M, (a+b)/2, 2*a, 2*b)
+ elseif prms.method === :mc
+ return dg1d.MC(b, a)
+ else
+ TODO(prms.method)
+ end
+end
+
+
+#######################################################################
+# slope limiters for FVElement mesh #
+#######################################################################
+
+
+# TODO Merge this with ScalarEq/rhs.jl:limit_slope!
+function limit_slope!(δu, u, mesh::Mesh1d{<:FVElement}, @nospecialize(prms::NamedTuple))
+ @unpack hrsc = prms
+ if hrsc === :muscl
+ limit_slope_muscl!(δu, u, mesh, prms)
+ elseif hrsc === :fv_central
+ # TODO Not sure if correct. The discussion in Hesthaven book 2018, sec 10.3 does not
+ # specify how to compute the slope. The only requirement seems to be that
+ # the slope obeys δu_j/h = du/dx(x_j) + O(h). And in the codes they provide they
+ # seem to reuse the same slope limiters as for the MUSCL scheme etc.
+ prms = (hrsc = :muscl, muscl_omega=0)
+ limit_slope_muscl!(δu, u, mesh, prms)
+ else
+ TODO(hrsc)
+ end
+end
+
+
+function limit_slope_muscl!(δu, u, mesh::Mesh1d{<:FVElement}, prms::NamedTuple)
+ # assume that FV solution in jth cell [x_(j-1/2),x_(j+1/2)] can be approximated by
+ # u_j ≈ ̄u_j + δu_j (x-x_j)/h
+ # where δu_j is the slope at x_j and
+ # x_j = x_j-1/2 + h/2
+ # h = x_j+1/2 - x_j-1/2
+ # Compare discussion in Hesthaven book 2018, sec 10.2
+ @unpack muscl_omega = prms
+ ω = muscl_omega
+ # TODO Enabling upwinding in fv_muscl_burgers_sine (ω = 1.0) gives asymmetric results.
+ # Do we need to compute two slopes like in the slope limiter implementations in EulerEq etc?
+ K, = mesh.tree.dims
+ @turbo for k = 2:K-1
+ kâ‚‹, kâ‚Š = k-1, k+1
+ uâ‚‹, uâ‚€, uâ‚Š = u[kâ‚‹], u[k], u[kâ‚Š]
+ Δâ»u = uâ‚€-uâ‚‹
+ Δâºu = uâ‚Š-uâ‚€
+ # (10.12) in Hesthaven book 2018, sec 10.2
+ δu[k] = (1+ω)/2*Δâ»u + (1-ω)/2*Δâºu
+ end
+ if isperiodic(mesh)
+ for k in (1,K)
+ kâ‚‹, kâ‚Š = periodic_index(k-1,K), periodic_index(k+1,K)
+ uâ‚‹, uâ‚€, uâ‚Š = u[kâ‚‹], u[k], u[kâ‚Š]
+ Δâ»u = uâ‚€-uâ‚‹
+ Δâºu = uâ‚Š-uâ‚€
+ δu[k] = (1+ω)/2*Δâ»u + (1-ω)/2*Δâºu
+ end
+ else
+ TODO("bdry conditions")
+ end
+end
+
+
+#######################################################################
+# slope limiter that preserves C^0 constraints #
+#######################################################################
+
+
+function limit_slope_preserve_continuity!(var::AbstractArray, mask::AbstractArray, ws::Workspace, mesh::Mesh1d{<:DGElement})
+ @toggled_assert mesh.element.kind === :modal_bspline1
+ enter(ws) do
+
+ @unpack x = get_static_variables(mesh)
+ Npts = mesh.element.Npts
+ K = mesh.tree.dims[1]
+ tmp_var = borrow(ws, Npts*K)
+ avgs = borrow(ws, Npts+1)
+ δs = borrow(ws, Npts+1)
+ lim_avgs = borrow(ws, Npts-1)
+ lim_δs = borrow(ws, Npts-1)
+ mat_A = borrow(ws, (Npts-1)*(Npts-1))
+ b = borrow(ws, Npts-1)
+ A = reshape(mat_A,(Npts-1,Npts-1))
+
+ for k in 1:K
+ if k == 1 || k == K || mask[k] ≈ 0.0
+ # skip bdry cells, but backup data because we restore the whole vector at once later
+ idxs = cellindices(mesh, k)
+ view(tmp_var, idxs) .= view(var, idxs)
+ continue
+ end
+
+ idxs_lhs, idxs, idxs_rhs = cellindices.(Ref(mesh), (k-1,k,k+1))
+ vvar_lhs, vvar, vvar_rhs = view.(Ref(var), (idxs_lhs,idxs,idxs_rhs))
+ vx_lhs, vx, vx_rhs = view.(Ref(x), (idxs_lhs,idxs,idxs_rhs))
+ vtmp_var = view(tmp_var, idxs)
+
+ # determine averages and slopes, including neighbors
+ v1, v2 = vvar_lhs[end-1], vvar_lhs[end]
+ x1, x2 = vx_lhs[end-1], vx_lhs[end]
+ avgs[1] = (v1+v2)/2; δs[1] = (v2-v1)/(x2-x1)
+ v1, v2 = vvar_rhs[1], vvar_rhs[2]
+ x1, x2 = vx_rhs[1], vx_rhs[2]
+ avgs[end] = (v1+v2)/2; δs[end] = (v2-v1)/(x2-x1)
+ for i in 1:Npts-1
+ v1, v2 = vvar[i], vvar[i+1]
+ x1, x2 = vx[i], vx[i+1]
+ avgs[1+i] = (v1+v2)/2; δs[1+i] = (v2-v1)/(x2-x1)
+ end
+
+ # limit slopes
+ for i in 2:Npts
+ # lim_δs[i-1] = TCI.minmod(δs[i-1],δs[i],δs[i+1])
+ x1, x2 = vx[i-1], vx[i]
+ dx = x2-x1
+ M = 0.1
+ δl, δc, δr = δs[i-1], δs[i], δs[i+1]
+ lim_δs[i-1] = TCI.minmod_M(dx, M, δc, δl, δr)
+ end
+
+ # recompute avgs so that continuity and total avg is preserved
+ fill!(A, 0.0) # because we use lu! below
+ for i in 1:Npts-2
+ x1, x2, x3 = vx[i], vx[i+1], vx[i+2]
+ dx1, dx2 = x2-x1, x3-x2
+ A[i,i] = 1.0
+ A[i,i+1] = -1.0
+ b[i] = - 0.5*dx1*lim_δs[i] - 0.5*dx2*lim_δs[i+1]
+ end
+ A[end,:] .= 1.0
+ v_avgs = view(avgs, 2:Npts)
+ b[end] = sum(v_avgs)
+
+ # solve the equations
+ F = lu!(A)
+ ldiv!(lim_avgs, F, b)
+
+ # resample data
+ x1, x2 = vx[1], vx[2]
+ dx = x2-x1
+ vtmp_var[1] = lim_avgs[1] - 0.5*dx*lim_δs[1]
+ x1, x2 = vx[end-1], vx[end]
+ dx = x2-x1
+ vtmp_var[end] = lim_avgs[end] + 0.5*dx*lim_δs[end]
+ for i in 2:Npts-1
+ x1, x2 = vx[i], vx[i+1]
+ dx = x2-x1
+ vtmp_var[i] = lim_avgs[i] - 0.5*dx*lim_δs[i]
+ end
+ end
+
+ copyto!(var, tmp_var)
+ end
+end
+
+
+function limit_slope_convex_hull!(var::AbstractArray, mask::AbstractArray,
+ ws::Workspace, mesh::Mesh1d{<:DGElement})
+ enter(ws) do
+
+ @unpack x = get_static_variables(mesh)
+ Npts = mesh.element.Npts
+ K = mesh.tree.dims[1]
+ tmp_var = borrow(ws, Npts*K)
+ dvar = borrow(ws, Npts-1)
+ w = 0.5
+
+ for k in 1:K
+
+ idxs = cellindices(mesh, k)
+ vvar = view(var, idxs)
+ vx = view(x, idxs)
+ vtmp_var = view(tmp_var, idxs)
+
+ if k == 1 || k == K || mask[k] ≈ 0.0
+ # skip bdry cells, but backup data because we restore the whole vector at once later
+ vtmp_var .= vvar
+ continue
+ end
+
+ for i in 1:Npts-1
+ dvar[i] = vvar[i+1]-vvar[i]
+ end
+ if abs(sum(sign, dvar)) ≈ Float64(Npts-1)
+ vtmp_var .= vvar
+ continue
+ end
+
+ idxs_lhs, idxs_rhs = cellindices.(Ref(mesh), (k-1,k+1))
+ vvar_lhs, vvar_rhs = view.(Ref(var), (idxs_lhs,idxs_rhs))
+
+ vtmp_var[1] = vvar[1] + (vvar_lhs[end]-vvar[1])*w
+ vtmp_var[end] = vvar[end] + (vvar_rhs[1]-vvar[end])*w
+ for i in 2:Npts-1
+ ul, uc, ur = vvar[i-1], vvar[i], vvar[i+1]
+ xl, xc, xr = vx[i-1], vx[i], vx[i+1]
+ um = ul + (ur-ul)/(xr-xl)*(xc-xl)
+ vtmp_var[i] = uc + (um-uc)*w
+ end
+ end
+
+ copyto!(var, tmp_var)
+ end
+end
diff --git a/src/mesh.jl b/src/mesh.jl
index d400bfd9d55e5c2623452e58883fce25543a30bb..9ef12398ba69c115f2e8a97e3eed6f0147496a70 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -920,59 +920,3 @@ function interpolate_face_data!(mesh::Mesh2d, var::AbstractVector, bdrybuf::Abst
bo += Nptsx
end
end
-
-
-#######################################################################
-# slope limiters for FVElement mesh #
-#######################################################################
-
-
-function limit_slope!(δu, u, mesh::Mesh1d{<:FVElement}, @nospecialize(prms::NamedTuple))
- @unpack hrsc = prms
- if hrsc === :muscl
- limit_slope_muscl!(δu, u, mesh, prms)
- elseif hrsc === :fv_central
- # TODO Not sure if correct. The discussion in Hesthaven book 2018, sec 10.3 does not
- # specify how to compute the slope. The only requirement seems to be that
- # the slope obeys δu_j/h = du/dx(x_j) + O(h). And in the codes they provide they
- # seem to reuse the same slope limiters as for the MUSCL scheme etc.
- prms = (hrsc = :muscl, muscl_omega=0)
- limit_slope_muscl!(δu, u, mesh, prms)
- else
- TODO(hrsc)
- end
-end
-
-
-function limit_slope_muscl!(δu, u, mesh::Mesh1d{<:FVElement}, prms::NamedTuple)
- # assume that FV solution in jth cell [x_(j-1/2),x_(j+1/2)] can be approximated by
- # u_j ≈ ̄u_j + δu_j (x-x_j)/h
- # where δu_j is the slope at x_j and
- # x_j = x_j-1/2 + h/2
- # h = x_j+1/2 - x_j-1/2
- # Compare discussion in Hesthaven book 2018, sec 10.2
- @unpack muscl_omega = prms
- ω = muscl_omega
- # TODO Enabling upwinding in fv_muscl_burgers_sine (ω = 1.0) gives asymmetric results.
- # Do we need to compute two slopes like in the slope limiter implementations in EulerEq etc?
- K, = mesh.tree.dims
- @turbo for k = 2:K-1
- kâ‚‹, kâ‚Š = k-1, k+1
- uâ‚‹, uâ‚€, uâ‚Š = u[kâ‚‹], u[k], u[kâ‚Š]
- Δâ»u = uâ‚€-uâ‚‹
- Δâºu = uâ‚Š-uâ‚€
- # (10.12) in Hesthaven book 2018, sec 10.2
- δu[k] = (1+ω)/2*Δâ»u + (1-ω)/2*Δâºu
- end
- if isperiodic(mesh)
- for k in (1,K)
- kâ‚‹, kâ‚Š = periodic_index(k-1,K), periodic_index(k+1,K)
- uâ‚‹, uâ‚€, uâ‚Š = u[kâ‚‹], u[k], u[kâ‚Š]
- Δâ»u = uâ‚€-uâ‚‹
- Δâºu = uâ‚Š-uâ‚€
- δu[k] = (1+ω)/2*Δâ»u + (1-ω)/2*Δâºu
- end
- else
- TODO("bdry conditions")
- end
-end
diff --git a/test/IntegrationTests/refs/burgers_sine_bspline1_limit_slope_convex_hull/burgers_sine_bspline1_limit_slope_convex_hull.toml b/test/IntegrationTests/refs/burgers_sine_bspline1_limit_slope_convex_hull/burgers_sine_bspline1_limit_slope_convex_hull.toml
new file mode 100644
index 0000000000000000000000000000000000000000..fc3b95a3f54f06caa9228445f72d07db4f2f2f49
--- /dev/null
+++ b/test/IntegrationTests/refs/burgers_sine_bspline1_limit_slope_convex_hull/burgers_sine_bspline1_limit_slope_convex_hull.toml
@@ -0,0 +1,23 @@
+[Evolution]
+cfl = 0.8
+tend = 0.24
+
+[ScalarEq]
+equation = "burgers"
+id = "sine"
+hrsc = "slope_limiter"
+
+[TCI]
+variables = ["u"]
+method = "mda"
+
+[Output]
+aligned_ts = "$(collect(range(0.02,0.24,step=0.02)))"
+variables1d = ["u"]
+
+[Mesh]
+dg_kind = "modal_bspline1"
+range = [0, 1]
+k = 13
+basis = "lgl"
+n = 7
diff --git a/test/IntegrationTests/refs/burgers_sine_bspline1_limit_slope_convex_hull/output1d.h5 b/test/IntegrationTests/refs/burgers_sine_bspline1_limit_slope_convex_hull/output1d.h5
new file mode 100644
index 0000000000000000000000000000000000000000..79ea9422160d52f88919d291a01f34d97a283563
Binary files /dev/null and b/test/IntegrationTests/refs/burgers_sine_bspline1_limit_slope_convex_hull/output1d.h5 differ
diff --git a/test/IntegrationTests/refs/burgers_sine_bspline2_limit_slope_convex_hull/burgers_sine_bspline2_limit_slope_convex_hull.toml b/test/IntegrationTests/refs/burgers_sine_bspline2_limit_slope_convex_hull/burgers_sine_bspline2_limit_slope_convex_hull.toml
new file mode 100644
index 0000000000000000000000000000000000000000..2849e536f24ede9677b05603a630cd600090bce4
--- /dev/null
+++ b/test/IntegrationTests/refs/burgers_sine_bspline2_limit_slope_convex_hull/burgers_sine_bspline2_limit_slope_convex_hull.toml
@@ -0,0 +1,23 @@
+[Evolution]
+cfl = 0.8
+tend = 0.24
+
+[ScalarEq]
+equation = "burgers"
+id = "sine"
+hrsc = "slope_limiter"
+
+[TCI]
+variables = ["u"]
+method = "mda"
+
+[Output]
+aligned_ts = "$(collect(range(0.02,0.24,step=0.02)))"
+variables1d = ["u"]
+
+[Mesh]
+dg_kind = "modal_bspline2"
+range = [0, 1]
+k = 14
+basis = "lgl"
+n = 7
diff --git a/test/IntegrationTests/refs/burgers_sine_bspline2_limit_slope_convex_hull/output1d.h5 b/test/IntegrationTests/refs/burgers_sine_bspline2_limit_slope_convex_hull/output1d.h5
new file mode 100644
index 0000000000000000000000000000000000000000..5e6c34bba2e750a0ed685d6e25bcb2ab679188c1
Binary files /dev/null and b/test/IntegrationTests/refs/burgers_sine_bspline2_limit_slope_convex_hull/output1d.h5 differ
diff --git a/test/IntegrationTests/refs/testconfig.toml b/test/IntegrationTests/refs/testconfig.toml
index fd03eac817b42f2d7a39325bf49540e2cf9e6b91..e9e6bbd4f0d74f502bb7652715d404d7bae1396d 100644
--- a/test/IntegrationTests/refs/testconfig.toml
+++ b/test/IntegrationTests/refs/testconfig.toml
@@ -61,6 +61,11 @@ variables1d = [ "u" ]
# this test now runs into the discontinuity
[burgers_sine_avmda]
variables1d = [ "u" ]
+# bspline tests
+[burgers_sine_bspline1_limit_slope_convex_hull]
+variables1d = [ "u" ]
+[burgers_sine_bspline2_limit_slope_convex_hull]
+variables1d = [ "u" ]
[balancedburgers_kink_avmda]
variables1d = [ "u" ]