ScalarEq: Implement modal Bspline1 DG evolution...

ScalarEq: Implement modal Bspline1 DG evolution (!216)

Similar to !213

Same caveats apply, e.g. no HRSC working etc.

src/ScalarEq/rhs.jl

@@ -268,7 +268,7 @@ end
 
 
 function rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
-  if mesh.element.kind === :modal_bspline2
+  if mesh.element.kind ∈ (:modal_bspline2, :modal_bspline1)
     modal_rhs!(env, mesh, P, hrsc)
   else
     nodal_rhs!(env, mesh, P, hrsc)

src/ScalarEq/setup.jl

@@ -14,8 +14,12 @@ 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.bs2)
+  if mesh.element isa DGElement && mesh.element.kind ∈ (:modal_bspline2, :modal_bspline1)
+    if mesh.element.kind === :modal_bspline2
+      V_bspline = dg1d.Bspline.vandermonde_matrix(mesh.element.z, mesh.element.bs2)
+    elseif mesh.element.kind === :modal_bspline1
+      V_bspline = dg1d.Bspline.vandermonde_matrix(mesh.element.z, mesh.element.bs1)
+    end
     invV_bspline = pinv(V_bspline)
   else
     V_bspline, invV_bspline = nothing, nothing
@@ -132,7 +136,7 @@ function register_equation!(mesh::Mesh1d, P::Project)
       TODO(P.prms.hrsc)
     end
   end
-  if mesh.element isa DGElement && mesh.element.kind === :modal_bspline2
+  if mesh.element isa DGElement && mesh.element.kind ∈ (:modal_bspline2, :modal_bspline1)
     register_variables!(mesh.cache, static_variablenames=(:u_modal,:flx_u_modal,:src_u_modal))
   end
 end

src/bspline.jl

@@ -362,7 +362,7 @@ function integrate(f::Function, idxs::AbstractVector{<:Integer}, bs::Bspline1)
   @toggled_assert length(idxs) == bs.Npts
   sum = 0.0
   @inbounds for i in bs.Npts
-    fi, xl, xr = f[idxs[i]], bs.xs[i], bs.xs[i+2]
+    fi, xl, xr = f(idxs[i]), bs.xs[i], bs.xs[i+2]
     sum += fi * (xr-xl)/2
   end
   return sum

src/dgelement.jl

@@ -83,7 +83,7 @@ struct DGElement
       purge_zeros!(quadr_z)
       z = collect(range(-1.0,1.0,Npts))
       purge_zeros!(z)
-      D = [ Bernstein.derivative_bernstein_polynomial(z[i], j-1, N) for i in 1:Npts, j in 1:Npts ]
+      D = [ Bernstein.derivative_bernstein_polynomial(z[i], j, N) for i in 1:Npts, j in 1:Npts ]
       M = Bernstein.mass_matrix(Npts)
       S = Bernstein.stiffness_matrix(Npts)
       invM = inv(M)
@@ -110,6 +110,7 @@ struct DGElement
       # for which a 4th order accurate LG rule should suffice.
       # We use a Gauss rule to avoid evaluating the Bspline on subcell interfaces and bdrys.
       quadr_z, quadr_w = LG.rule(3)
+      # TODO Derivatives at bspline control points is ill-defined, because the splines are C^0 there
       D = [ Bspline.derivative_polynomial(z[i], j, bs1) for i in 1:Npts, j in 1:Npts ]
       M = Bspline.mass_matrix(bs1)
       S = Bspline.stiffness_matrix(bs1)
@@ -210,8 +211,10 @@ end
 
 
 @inline function integrate(u::AbstractVector, el::DGElement)
-  if el.kind ∈ (:modal_bspline2, :modal_bspline1)
-    return Bspline.integrate(u, el.bs)
+  if el.kind === :modal_bspline1
+    return Bspline.integrate(u, el.bs1)
+  elseif el.kind === :modal_bspline2
+    return Bspline.integrate(u, el.bs2)
   else
     if el.quadr === :lgl || el.quadr === :lgl_lumped
       return LGL.integrate(u, el.quadr_z, el.quadr_w)
@@ -223,8 +226,10 @@ end
 
 
 @inline function integrate(fn::Function, idxs::AbstractVector{<:Integer}, el::DGElement)
-  if el.kind ∈ (:modal_bspline2, :modal_bspline1)
-    return Bspline.integrate(fn, idxs, el.bs)
+  if el.kind === :modal_bspline1
+    return Bspline.integrate(fn, idxs, el.bs1)
+  elseif el.kind === :modal_bspline2
+    return Bspline.integrate(fn, idxs, el.bs2)
   else
     if el.quadr === :lgl || el.quadr === :lgl_lumped
       return LGL.integrate(fn, idxs, el.quadr_z, el.quadr_w)

test/IntegrationTests/refs/modal_bspline1_advection_sine/modal_bspline1_advection_sine.toml

+[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_bspline1"
+range = [0.0, 1.0]
+k = 30
+basis = "lgl"
+n = 5

test/IntegrationTests/refs/testconfig.toml

@@ -10,6 +10,9 @@ variables1d = [ "u", "u_err" ]
 [modal_bernstein_advection_sine]
 variables1d = [ "u" ]
 
+[modal_bspline1_advection_sine]
+variables1d = [ "u", "u_err" ]
+
 [modal_bspline2_advection_sine]
 variables1d = [ "u", "u_err" ]
 

test/IntegrationTests/src/IntegrationTests.jl

@@ -529,7 +529,9 @@ function impl_runtests(tests::Vector{String}=String[])
         (data, i, j) -> (j == 3 || j == 4) && !data[i,j],
         crayon"red"
     )
-    pretty_table(logs, hcat(1:length(tests), tests, isnothing.(runs_status), comparsion_success);
+    test_ran = isnothing.(runs_status)
+    test_passed = map(s -> ismissing(s) ? false : s, comparsion_success)
+    pretty_table(logs, hcat(1:length(tests), tests, test_ran, test_passed);
                  tf = tf_unicode_rounded,
                  alignment = [ :l, :l, :r, :r ],
                  highlighters = (hl_success,hl_fail),