diff --git a/src/ScalarEq/rhs.jl b/src/ScalarEq/rhs.jl
index 75100c42b3a8369bd1470fe804dc7319af3a2d23..8f4f038abe838aee08e9f4da3d0be00f1815e35e 100644
--- a/src/ScalarEq/rhs.jl
+++ b/src/ScalarEq/rhs.jl
@@ -268,21 +268,54 @@ end
function rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
+ if mesh.element.kind === :modal_bspline2
+ modal_rhs!(env, mesh, P, hrsc)
+ else
+ nodal_rhs!(env, mesh, P, hrsc)
+ end
+end
- @unpack cache, mesh = env
- @unpack equation = P
- @unpack flx_u, src_u = get_static_variables(cache)
- @unpack u = get_dynamic_variables(cache)
- @unpack rhs_u = get_rhs_variables(cache)
- @unpack nflx_u, bdry_u = get_bdry_variables(cache)
+function nodal_rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
+ @unpack flx_u, src_u = get_static_variables(mesh)
+ @unpack u = get_dynamic_variables(mesh)
+ @unpack rhs_u = get_rhs_variables(mesh)
+ @unpack nflx_u, bdry_u = get_bdry_variables(mesh)
dg1d.interpolate_face_data!(mesh, u, bdry_u)
- broadcast_volume!(flux_source, equation, mesh)
- broadcast_faces!(llf_1d, equation, mesh)
- # broadcast_bdry!(bdryllf_1d, equation, P.bdrycond, mesh)
+ broadcast_volume!(flux_source, P.equation, mesh)
+ broadcast_faces!(llf_1d, P.equation, mesh)
compute_rhs_weak_form!(rhs_u, flx_u, src_u, nflx_u, mesh)
+ return
+end
+
+
+function modal_rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
+ @unpack flx_u, src_u,
+ u_modal, flx_u_modal,
+ src_u_modal = get_static_variables(mesh)
+ @unpack u = get_dynamic_variables(mesh)
+ @unpack rhs_u = get_rhs_variables(mesh)
+ @unpack nflx_u, bdry_u = get_bdry_variables(mesh)
+
+ 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)
+ end
+
+ dg1d.interpolate_face_data!(mesh, u, bdry_u)
+ broadcast_volume!(flux_source, P.equation, mesh)
+ broadcast_faces!(llf_1d, P.equation, mesh)
+
+ for (v_flx_u, v_flx_u_modal) in zip(eachcell(mesh,flx_u),eachcell(mesh,flx_u_modal))
+ mul!(v_flx_u_modal, P.prms.invV_bspline, v_flx_u)
+ end
+ for (v_src_u, v_src_u_modal) in zip(eachcell(mesh,src_u),eachcell(mesh,src_u_modal))
+ 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
end
diff --git a/src/ScalarEq/setup.jl b/src/ScalarEq/setup.jl
index 690f9b41c51a3d5d41f12405602125516820e6f9..7e85208e48a900fb8ecbf442ecc736491f85eab9 100644
--- a/src/ScalarEq/setup.jl
+++ b/src/ScalarEq/setup.jl
@@ -14,10 +14,16 @@ function Project(env::Environment, mesh::Mesh1d, prms)
analyze_error = prms["ScalarEq"]["analyze_error"]
analyze_error_norm = prms["ScalarEq"]["analyze_error_norm"]
V_bernstein = dg1d.Bernstein.vandermonde_matrix(mesh.element.z)
+ if mesh.element isa DGElement && mesh.element.kind === :modal_bspline2
+ V_bspline = dg1d.Bspline.vandermonde_matrix(mesh.element.z, mesh.element.bs)
+ invV_bspline = pinv(V_bspline)
+ else
+ V_bspline, invV_bspline = nothing, nothing
+ end
fixedprms = (; bernstein, av_derivative_scheme, av_drag, hrsc, muscl_omega,
slope_limiter_method, slope_limiter_tvb_M,
av_recompute_substeps,
- V_bernstein,
+ V_bernstein, V_bspline, invV_bspline,
analyze_error, analyze_error_norm)
# construct Project
@@ -126,6 +132,9 @@ function register_equation!(mesh::Mesh1d, P::Project)
TODO(P.prms.hrsc)
end
end
+ if mesh.element isa DGElement && mesh.element.kind === :modal_bspline2
+ register_variables!(mesh.cache, static_variablenames=(:u_modal,:flx_u_modal,:src_u_modal))
+ end
end
function register_equation!(mesh::Mesh2d, P::Project)
register_variables!(mesh.cache,
diff --git a/src/bspline.jl b/src/bspline.jl
index 2ac9014014547b404ef5d0131e9eb67b0b391d29..281336dc93e15b97bf7622234b1521fd436b2006 100644
--- a/src/bspline.jl
+++ b/src/bspline.jl
@@ -108,12 +108,13 @@ end
function unsafe_interpolate(x::Real, coeffs::AbstractVector, bs::Bspline2)
@unpack Npts, xs = bs
@inbounds xl, xr = xs[1], xs[end]
- xl ≤ x ≤ xr || return 0.0
+ xl == x && return coeffs[1]
+ xr == x && return coeffs[end]
+ xl < x < xr || return 0.0
o = 2 # offset duplicated bdry knots
v_xs = view(xs, o .+ (1:Npts))
ir = searchsortedfirst(v_xs, x)
il = ir-1
- il < 1 && return coeffs[1] # x == xl
is_supp = max(il-1,1):min(ir+1,Npts+1)
B = 0.0
@inbounds for i in is_supp
diff --git a/test/IntegrationTests/refs/modal_bspline2_advection_sine/modal_bspline2_advection_sine.toml b/test/IntegrationTests/refs/modal_bspline2_advection_sine/modal_bspline2_advection_sine.toml
new file mode 100644
index 0000000000000000000000000000000000000000..3d670bb74227fed93e1f0d23722e299d208b3a0c
--- /dev/null
+++ b/test/IntegrationTests/refs/modal_bspline2_advection_sine/modal_bspline2_advection_sine.toml
@@ -0,0 +1,21 @@
+[Evolution]
+cfl = 0.5
+tend = 2.2
+
+[ScalarEq]
+analyze_error = true
+equation = "advection"
+id = "sine"
+
+[Output]
+substep_every_iteration = 200
+substep_variables = ["u", "u_err"]
+aligned_ts = "$(collect(range(0.1,2.2,step=0.1)))"
+variables1d = ["u", "u_err"]
+
+[Mesh]
+dg_kind = "modal_bspline2"
+range = [0.0, 1.0]
+k = 30
+basis = "lgl"
+n = 5
diff --git a/test/IntegrationTests/refs/modal_bspline2_advection_sine/output1d.h5 b/test/IntegrationTests/refs/modal_bspline2_advection_sine/output1d.h5
new file mode 100644
index 0000000000000000000000000000000000000000..d0e6758edad6c361ad5196a967dd777fb917170e
Binary files /dev/null and b/test/IntegrationTests/refs/modal_bspline2_advection_sine/output1d.h5 differ
diff --git a/test/IntegrationTests/refs/modal_bspline2_advection_sine/substep_output.h5 b/test/IntegrationTests/refs/modal_bspline2_advection_sine/substep_output.h5
new file mode 100644
index 0000000000000000000000000000000000000000..14ee01b05bc913a144ae5617141275c5929794b4
Binary files /dev/null and b/test/IntegrationTests/refs/modal_bspline2_advection_sine/substep_output.h5 differ
diff --git a/test/IntegrationTests/refs/testconfig.toml b/test/IntegrationTests/refs/testconfig.toml
index b5cdbccf82d11f4257b10818dcd97077dd2bfe0b..f9687bb2450cee77da23cdc94519931fcb9c79f3 100644
--- a/test/IntegrationTests/refs/testconfig.toml
+++ b/test/IntegrationTests/refs/testconfig.toml
@@ -10,6 +10,9 @@ variables1d = [ "u", "u_err" ]
[modal_bernstein_advection_sine]
variables1d = [ "u" ]
+[modal_bspline2_advection_sine]
+variables1d = [ "u", "u_err" ]
+
[advection_sine_analysis]
variables0d = [ "u_err_norm", "u_norm2" ]
variables1d = [ "u" ]