diff --git a/.githooks/pre-push b/.githooks/pre-push new file mode 100755 index 0000000..239da8f --- /dev/null +++ b/.githooks/pre-push @@ -0,0 +1,57 @@ +#!/bin/sh +# +# Run the test suite before a push, and refuse the push if it fails. +# +# Enable with +# +# git config core.hooksPath .githooks +# +# and skip a single push with `SIMPLESPLINES_SKIP_TESTS=1 git push` — for a +# documentation-only change, say, or when the failure is already known and being +# pushed to a branch on purpose. + +set -eu + +# git feeds the hook one ` ` line per ref +# being pushed, and an all-zero local oid means that ref is being deleted. Read them: +# a push that only deletes has nothing to test, and a hook that exits without draining +# leaves git writing into a closed pipe. +refs=0 +content=0 +while read -r _local_ref local_oid _remote_ref _remote_oid; do + refs=$((refs + 1)) + case "$local_oid" in + *[!0]*) content=1 ;; + esac +done + +if [ "$refs" -gt 0 ] && [ "$content" -eq 0 ]; then + printf 'pre-push: nothing but ref deletions, skipping the test suite.\n' >&2 + exit 0 +fi + +if [ "${SIMPLESPLINES_SKIP_TESTS:-0}" != "0" ]; then + printf 'pre-push: SIMPLESPLINES_SKIP_TESTS is set, skipping the test suite.\n' >&2 + exit 0 +fi + +if ! command -v julia >/dev/null 2>&1; then + printf 'pre-push: julia is not on PATH, cannot run the test suite; refusing the push.\n' >&2 + printf ' Set SIMPLESPLINES_SKIP_TESTS=1 to push anyway.\n' >&2 + exit 1 +fi + +# git runs its hooks from the top level of the work tree, but asking rather than assuming +# also gives the right answer in a linked worktree. +root=$(git rev-parse --show-toplevel) + +printf 'pre-push: running the SimpleSplines test suite...\n' >&2 + +if julia --project="$root" -e 'using Pkg; Pkg.test()'; then + printf 'pre-push: tests passed.\n' >&2 + exit 0 +fi + +printf 'pre-push: tests FAILED; refusing the push.\n' >&2 +printf ' Set SIMPLESPLINES_SKIP_TESTS=1 to push anyway.\n' >&2 +exit 1 diff --git a/.github/dependabot.yml b/.github/dependabot.yml index 700707c..bebf22d 100644 --- a/.github/dependabot.yml +++ b/.github/dependabot.yml @@ -5,3 +5,8 @@ updates: directory: "/" # Location of package manifests schedule: interval: "weekly" + # One pull request for all action bumps rather than one per action. + groups: + github-actions: + patterns: + - "*" diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 923306b..3d10615 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -16,65 +16,55 @@ jobs: name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} runs-on: ${{ matrix.os }} timeout-minutes: 60 + # A pre-release or nightly failure is information, not a broken build. + continue-on-error: ${{ matrix.experimental }} permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created actions: write contents: read strategy: fail-fast: false matrix: + # 'min' resolves the lower bound of the julia compat entry in Project.toml, so the + # matrix follows the declared support window instead of having to be edited alongside + # it; '1' is the current stable release. version: - - '1.10' - - '1.12' - - 'nightly' + - 'min' + - '1' os: - ubuntu-latest - macOS-latest - windows-latest arch: - - x64 + # 'default' is the runner's native architecture, which is aarch64 on macOS-latest; + # a hardcoded 'x64' would test Julia under Rosetta there. + - default + experimental: + - false + include: + - version: 'lts' + os: ubuntu-latest + arch: default + experimental: false + - version: 'pre' + os: ubuntu-latest + arch: default + experimental: true + - version: 'nightly' + os: ubuntu-latest + arch: default + experimental: true steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 + - uses: actions/checkout@v7 + - uses: julia-actions/setup-julia@v3 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v2 + - uses: julia-actions/cache@v3 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v4 + - uses: codecov/codecov-action@v7 with: files: lcov.info token: ${{ secrets.CODECOV_TOKEN }} fail_ci_if_error: false - docs: - name: Documentation - runs-on: ubuntu-latest - permissions: - actions: write # needed to allow julia-actions/cache to proactively delete old caches that it has created - contents: write - statuses: write - steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 - with: - version: '1' - - uses: julia-actions/cache@v2 - - name: Configure doc environment - shell: julia --project=docs --color=yes {0} - run: | - using Pkg - Pkg.develop(PackageSpec(path=pwd())) - Pkg.instantiate() - - uses: julia-actions/julia-buildpkg@v1 - - uses: julia-actions/julia-docdeploy@v1 - env: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} - - name: Run doctests - shell: julia --project=docs --color=yes {0} - run: | - using Documenter: DocMeta, doctest - using SimpleSplines - DocMeta.setdocmeta!(SimpleSplines, :DocTestSetup, :(using SimpleSplines); recursive=true) - doctest(SimpleSplines) diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml index cba9134..18352eb 100644 --- a/.github/workflows/CompatHelper.yml +++ b/.github/workflows/CompatHelper.yml @@ -3,10 +3,21 @@ on: schedule: - cron: 0 0 * * * workflow_dispatch: +permissions: + actions: write # needed to allow julia-actions/cache to proactively delete old caches that it has created + contents: write + pull-requests: write jobs: CompatHelper: runs-on: ubuntu-latest + timeout-minutes: 30 steps: + # The runner images no longer ship a Julia, so one has to be installed rather than + # assumed; without this the `julia -e` steps below fail before CompatHelper starts. + - uses: julia-actions/setup-julia@v3 + with: + version: '1' + - uses: julia-actions/cache@v3 - name: Pkg.add("CompatHelper") run: julia -e 'using Pkg; Pkg.add("CompatHelper")' - name: CompatHelper.main() diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml new file mode 100644 index 0000000..cbb6e17 --- /dev/null +++ b/.github/workflows/Documenter.yml @@ -0,0 +1,43 @@ +name: Documentation +on: + push: + branches: + - main + tags: ['*'] + pull_request: + workflow_dispatch: +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} +jobs: + build: + name: Documentation + runs-on: ubuntu-latest + timeout-minutes: 60 + permissions: + actions: write # needed to allow julia-actions/cache to proactively delete old caches that it has created + contents: write + pull-requests: read + statuses: write + steps: + - uses: actions/checkout@v7 + - uses: julia-actions/setup-julia@v3 + with: + version: '1' + - uses: julia-actions/cache@v3 + - uses: julia-actions/julia-buildpkg@v1 + # julia-docdeploy resolves the docs environment against the checkout and instantiates + # it itself, so there is no separate step for that here. + - uses: julia-actions/julia-docdeploy@v1 + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} + - name: Run doctests + shell: julia --project=docs --color=yes {0} + run: | + using Documenter: DocMeta, doctest + using SimpleSplines + DocMeta.setdocmeta!(SimpleSplines, :DocTestSetup, :(using SimpleSplines); recursive=true) + doctest(SimpleSplines) diff --git a/.github/workflows/Register.yml b/.github/workflows/Register.yml index 5b7cd3b..0367532 100644 --- a/.github/workflows/Register.yml +++ b/.github/workflows/Register.yml @@ -8,9 +8,10 @@ on: jobs: register: runs-on: ubuntu-latest + timeout-minutes: 30 permissions: - contents: write + contents: write steps: - - uses: julia-actions/RegisterAction@latest + - uses: julia-actions/RegisterAction@v0.3.2 with: token: ${{ secrets.GITHUB_TOKEN }} diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml index 0cd3114..234915d 100644 --- a/.github/workflows/TagBot.yml +++ b/.github/workflows/TagBot.yml @@ -24,6 +24,7 @@ jobs: TagBot: if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' runs-on: ubuntu-latest + timeout-minutes: 30 steps: - uses: JuliaRegistries/TagBot@v1 with: diff --git a/AUTHORS.md b/AUTHORS.md new file mode 100644 index 0000000..81447b3 --- /dev/null +++ b/AUTHORS.md @@ -0,0 +1,19 @@ +# Authors + +SimpleSplines' development is coordinated by a group of *principal developers*, who are also +its main contributors and who can be contacted in case of questions about SimpleSplines. In +addition, there are *contributors* who have provided substantial additions or modifications. +Together, these two groups form "The SimpleSplines Authors" as mentioned in the +[LICENSE](LICENSE.md) file. + +## Principal Developers + +* [Michael Kraus](https://www.michael-kraus.org/), + Max Planck Institute for Plasma Physics, Garching, Germany + +## Contributors + +Everyone who has contributed to SimpleSplines, the principal developers above included, in +alphabetical order: + +* Michael Kraus diff --git a/CHANGELOG.md b/CHANGELOG.md index 33cf2f3..b5a017f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -86,7 +86,22 @@ was the difference between 6.0 ms and 0.4 ms for assembling one Hessian. The constant assemblies — `mass_matrix`, `stiffness_matrix`, `derivative_matrix` and any `mixed_matrix` — are memoised on first use. They do not depend on the field, but a downstream -time integrator asks for them inside every Newton iteration of every step. +time integrator asks for them inside every Newton iteration of every step. `basis_integrals` +is one of these constants too, and is now assembled with the quadrature and returned by +reference rather than recomputed per call. + +The paths a time integrator actually runs per step allocate nothing. Evaluation is +allocation-free at any degree and derivative order on any mesh; a `CirculantMass` solve is +allocation-free; and `l2_projection!` is now allocation-free end to end on a uniform mesh, +the `f ⊙ w` product going into a buffer held by the quadrature and the load vector being +formed with `mul!` straight into the output. The buffer is taken only when the product lands +in the quadrature's element type, so a wider sample — a complex `f` — is still projected, +through a product of its own, rather than narrowed into it. On a non-uniform mesh a CHOLMOD +temporary remains, CHOLMOD having no in-place `ldiv!`. + +That buffer is mutable state on a struct that reads as immutable, as the memoising `cache` +behind `mixed_matrix` already was, so the `SplineQuadrature` docstring now warns that one +quadrature must not be shared between threads. ### Dependencies @@ -94,9 +109,74 @@ time integrator asks for them inside every Newton iteration of every step. assemblies do not map onto that package's Galerkin interface. `FFTW` and `SparseArrays` were added for the above, and `CompactBasisFunctions` for the shared `Basis` hierarchy. +### Repository and CI + +The documentation build moved out of `CI.yml` into its own `Documenter.yml` workflow, so that +a docs failure and a test failure are separate signals and the docs job is not queued behind +the test matrix. The README gained a badge for it. + +The workflows were brought up to current action versions — `actions/checkout@v7`, +`julia-actions/setup-julia@v3`, `julia-actions/cache@v3`, `codecov/codecov-action@v7` — and +the test matrix now names Julia versions by alias rather than by number: + +- `min` resolves the lower bound of the `julia` compat entry, so the matrix tracks the + declared support window instead of having to be edited alongside it; +- `lts` and `1` cover the long-term-support and current stable releases; +- `pre` and `nightly` run on Linux only and are `continue-on-error`, since an upcoming-release + failure is information rather than a broken build. + +`arch` is now `default` rather than `x64`, which is what tests Julia natively on the ARM64 +macOS runners instead of under Rosetta. + +`AUTHORS.md` was added and `LICENSE` renamed to `LICENSE.md`, whose copyright line now names +"The SimpleSplines Authors" and points at it. + +The `pre-push` hook the README asks the reader to enable now exists in `.githooks/`. It runs +the test suite and refuses the push if it fails; `SIMPLESPLINES_SKIP_TESTS=1` overrides it. + ### Fixed - `QuadratureRules` compat was `"0.1"`, which could not co-resolve with `CompactBasisFunctions`; it is now `"0.2"`. - `LinearAlgebra` compat was `"1.12.0"` alongside `julia = "1.10"`, which contradicted the 1.10 row of the CI matrix; it is now `"1"`. +- `docs/Project.toml` pinned `CompactBasisFunctions` to an absolute path on a developer's + machine through a `[sources]` entry. That path does not exist on a CI runner, so the + documentation build could only ever have succeeded locally; the entry is removed and the + dependency now resolves from the registry. +- `CompatHelper.yml` invoked `julia` without installing it. The runner images no longer ship + a Julia, so the workflow failed before CompatHelper started; it now sets Julia up first. + +## Open Issues + +### `weighted_matrix` allocates a fresh matrix per call + +`weighted_matrix` is the one assembly that cannot be memoised — it depends on the field, so a +downstream time integrator asks for a *different* one inside every Newton iteration of every +step, which is exactly the call pattern under which allocation matters most. Measured at +`N = 128`, `p = 3`, `nq = 5`: + +| | bytes | +|:--|--:| +| `f .* q.w` temporary | 5 kB | +| `Φₐ * Diagonal(f ⊙ w)` | 47 kB | +| `(Φₐ D) * Φᵦᵀ` sparse-sparse product | 208 kB | +| **total per call** | **260 kB** | + +The `f ⊙ w` temporary is the same one `l2_projection!` no longer pays and could be removed the +same way, with the buffer the quadrature already holds. The other 255 kB are not a temporary +at all: they are the result, a freshly built `SparseMatrixCSC` with its `colptr`, `rowval` and +`nzval` allocated and its structure recomputed from scratch. + +That structure does not depend on `f`. For a fixed `(a, b)` the sparsity pattern of +`Φₐ diag(f ⊙ w) Φᵦᵀ` is the same for every coefficient — a basis function overlaps only the +`2p+1` others whose supports meet its own — so the pattern could be assembled once per +`(a, b)`, cached beside the `mixed_matrix` results, and only `nzval` refilled per +call. That turns 260 kB into zero. + +What it needs is an in-place entry point, `weighted_matrix!(A, q, f, a, b)`, since the present +signature has nowhere to write. Callers holding a matrix across steps would use it and callers +wanting a value would keep the allocating form. Deferred rather than done because it widens +the API, and because the sparse triple product would have to be written out by hand against +the cached pattern instead of delegating to `SparseArrays`, which is the part that needs to be +got right rather than merely written. diff --git a/LICENSE b/LICENSE.md similarity index 92% rename from LICENSE rename to LICENSE.md index b472609..70b5dd4 100644 --- a/LICENSE +++ b/LICENSE.md @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2025 Michael Kraus +Copyright (c) 2025-present The SimpleSplines Authors (see [AUTHORS.md](AUTHORS.md)) Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/README.md b/README.md index 0bae0c0..b8f18ed 100644 --- a/README.md +++ b/README.md @@ -3,6 +3,7 @@ [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaDEC.github.io/SimpleSplines.jl/stable/) [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaDEC.github.io/SimpleSplines.jl/dev/) [![Build Status](https://github.com/JuliaDEC/SimpleSplines.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/JuliaDEC/SimpleSplines.jl/actions/workflows/CI.yml?query=branch%3Amain) +[![Documentation](https://github.com/JuliaDEC/SimpleSplines.jl/actions/workflows/Documenter.yml/badge.svg?branch=main)](https://github.com/JuliaDEC/SimpleSplines.jl/actions/workflows/Documenter.yml?query=branch%3Amain) [![Coverage](https://codecov.io/gh/JuliaDEC/SimpleSplines.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/JuliaDEC/SimpleSplines.jl) [![PkgEval](https://JuliaCI.github.io/NanosoldierReports/pkgeval_badges/S/SimpleSplines.svg)](https://JuliaCI.github.io/NanosoldierReports/pkgeval_badges/S/SimpleSplines.html) @@ -19,3 +20,15 @@ To run the test suite before every push, enable the repository's git hooks: ```sh git config core.hooksPath .githooks ``` + +`.githooks/pre-push` refuses the push if the suite fails. To push regardless — a +documentation-only change, or a failure that is already known — set + +```sh +SIMPLESPLINES_SKIP_TESTS=1 git push +``` + +## License + +SimpleSplines is licensed under the [MIT License](LICENSE.md). See [AUTHORS.md](AUTHORS.md) +for the list of authors the copyright refers to. diff --git a/docs/Project.toml b/docs/Project.toml index ef2915a..7b31ff2 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -4,5 +4,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" SimpleSplines = "ac480f6c-9a94-4d53-a283-9db69f109154" [sources] -CompactBasisFunctions = {path = "/Users/mkraus/Datashare/Julia/CompactBasisFunctions"} SimpleSplines = {path = ".."} + +[compat] +Documenter = "1" diff --git a/src/quadrature.jl b/src/quadrature.jl index b17427a..184dec0 100644 --- a/src/quadrature.jl +++ b/src/quadrature.jl @@ -70,6 +70,13 @@ true support of `p+1` cells are ever nonzero and only those are computed, but keeping the array dense lets the contractions above run as one BLAS call, which is the faster arrangement at the sizes these discretisations are used at. + +!!! warning "One quadrature per thread" + A `SplineQuadrature` carries mutable state behind an otherwise read-only interface: + [`mixed_matrix`](@ref) memoises its results into a `Dict`, and [`l2_projection!`](@ref) + forms its ``f \odot w`` product in a shared buffer. Neither is synchronised, so one + quadrature must not be used from two threads at once. Give each thread its own, or stay + with the allocating [`l2_projection`](@ref), which forms its own product. """ struct SplineQuadrature{T, BT <: PeriodicBSplineBasis{T}, MO <: MassOperator{T}} basis::BT @@ -78,6 +85,8 @@ struct SplineQuadrature{T, BT <: PeriodicBSplineBasis{T}, MO <: MassOperator{T}} w::Vector{T} Φ::Vector{SparseMatrixCSC{T, Int}} mass::MO + integrals::Vector{T} + scratch::Vector{T} cache::Dict{Tuple{Int,Int}, SparseMatrixCSC{T, Int}} function SplineQuadrature(basis::BT; @@ -149,7 +158,13 @@ struct SplineQuadrature{T, BT <: PeriodicBSplineBasis{T}, MO <: MassOperator{T}} rethrow() end - new{T, BT, typeof(mass)}(basis, Int(nq), x, w, Φ, mass, + # ∫ φ_i dx is a constant of the discretisation like the mass matrix, so it is + # assembled here rather than recomputed on every call. `scratch` is the f ⊙ w + # buffer that keeps `l2_projection!` from allocating one per call. + integrals = Φ[1] * w + scratch = Vector{T}(undef, n * nq) + + new{T, BT, typeof(mass)}(basis, Int(nq), x, w, Φ, mass, integrals, scratch, Dict{Tuple{Int,Int}, SparseMatrixCSC{T, Int}}()) end end @@ -306,8 +321,11 @@ This is the gradient of the total mass ``C_0 = \int_\Omega u \, dx`` with respec degrees of freedom, and it equals ``\mathbb{M} \mathbf{1}`` because the basis is a partition of unity. It spans the kernel of the first discrete bracket, which is why the mass is a Casimir there. + +Assembled once when the [`SplineQuadrature`](@ref) is built and returned by reference, as +[`mass_matrix`](@ref) is — do not mutate the result. """ -basis_integrals(q::SplineQuadrature) = basis_values(q, 0) * q.w +basis_integrals(q::SplineQuadrature) = q.integrals @doc raw""" l2_projection(q::SplineQuadrature, f) @@ -339,11 +357,33 @@ end l2_projection!(û, q::SplineQuadrature, f) In-place [`l2_projection`](@ref), writing the coefficients into `û`. + +On a [`UniformMesh`](@ref) this allocates nothing: the ``f \odot w`` product goes into a +buffer held by the quadrature, the load vector is formed with `mul!` straight into `û`, and +the [`CirculantMass`](@ref) solve is itself allocation-free. On a non-uniform mesh the +CHOLMOD solve still allocates a temporary, as [`mass_solve!`](@ref) notes. + +A sample whose element type is wider than the quadrature's — a complex `f` — gets its own +product instead of being narrowed into that buffer, so this method accepts exactly what +[`l2_projection`](@ref) accepts and merely stops being allocation-free there. + +The shared buffer is one of the two reasons a [`SplineQuadrature`](@ref) may not be used from +two threads at once; see the warning there. """ function l2_projection!(û::AbstractVector, q::SplineQuadrature, f::AbstractVector) length(û) == nbasis(q) || throw(DimensionMismatch( "the coefficient vector has $(length(û)) entries but the basis has $(nbasis(q))")) - û .= basis_values(q, 0) * (q.w .* f) + length(f) == length(q.x) || throw(DimensionMismatch( + "the function was sampled at $(length(f)) points but the quadrature has " * + "$(length(q.x))")) + # The buffer has the quadrature's element type, so it can only take the product when the + # product lands in that type; a complex or extended-precision sample gets its own array + # rather than an `InexactError` or a silent narrowing. The test is on types alone, hence + # resolved when the method is compiled, so the ordinary path still allocates nothing. + S = promote_type(eltype(q.w), eltype(f)) + fw = S === eltype(q.scratch) ? q.scratch : similar(f, S) + fw .= q.w .* f + mul!(û, basis_values(q, 0), fw) mass_solve!(û, q.mass, û) return û end diff --git a/test/quadrature_tests.jl b/test/quadrature_tests.jl index 7ef5adc..5fe58e0 100644 --- a/test/quadrature_tests.jl +++ b/test/quadrature_tests.jl @@ -59,6 +59,12 @@ using Test @test mass_matrix(q) * ones(16) ≈ Iv @test all(>(0), Iv) end + # a constant of the discretisation, assembled once and returned by reference rather + # than recomputed on every call + q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(16), 3)) + @test basis_integrals(q) === basis_integrals(q) + basis_integrals(q) + @test (@allocated basis_integrals(q)) == 0 end @testset "$(rpad("S is already antisymmetric, given nq >= p",76))" begin @@ -156,6 +162,35 @@ using Test @test l2_projection!(ŵ, q, sin) ≈ û @test_throws DimensionMismatch l2_projection(q, [1.0, 2.0]) @test_throws DimensionMismatch l2_projection!(zeros(3), q, sin) + @test_throws DimensionMismatch l2_projection!(similar(û), q, [1.0, 2.0]) + + # on a uniform mesh the whole projection is allocation-free: the f ⊙ w buffer lives + # in the quadrature, the load vector is formed with mul! straight into û, and the + # CirculantMass solve allocates nothing + fq = sin.(quadrature_nodes(q)) + l2_projection!(ŵ, q, fq) + @test (@allocated l2_projection!(ŵ, q, fq)) == 0 + @test ŵ ≈ û + + # The in-place and allocating methods are separate implementations -- the buffer is + # the whole point of the split -- so they are checked against each other on every + # mesh family, and not only on the uniform mesh the allocation test above needs. + for (nm, mk) in MESHES, p in 1:4 + qm = SplineQuadrature(PeriodicBSplineBasis(mk(16), p)) + fm = sin.(quadrature_nodes(qm)) + ûm = Vector{Float64}(undef, nbasis(basis(qm))) + @test l2_projection!(ûm, qm, fm) ≈ l2_projection(qm, fm) + end + + # A sample wider than the quadrature's element type gets its own f ⊙ w product + # instead of being narrowed into the shared buffer, which would be an InexactError + # here and a silent loss of precision for a BigFloat. Graded rather than uniform: + # the rfft plan of a CirculantMass takes a real argument only. + qc = SplineQuadrature(PeriodicBSplineBasis(GradedMesh(16, 2π), 3)) + fc = cis.(quadrature_nodes(qc)) + ûc = Vector{ComplexF64}(undef, nbasis(basis(qc))) + @test l2_projection!(ûc, qc, fc) ≈ l2_projection(qc, fc) + @test eltype(l2_projection(qc, fc)) == ComplexF64 end @testset "$(rpad("L2 projection converges at order p+1",76))" begin