131 lines
3.7 KiB
Go
131 lines
3.7 KiB
Go
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
||
// SPDX-License-Identifier: MIT
|
||
|
||
package signal
|
||
|
||
import (
|
||
"math"
|
||
"testing"
|
||
|
||
"sourcedock.dev/petrbalvin/tensor/internal/core"
|
||
)
|
||
|
||
// poissonManufactured builds f = −Δu for u = sin(πx)·sin(πy) on the
|
||
// unit square over an (n × n) grid including the boundary, so the
|
||
// exact solution is known everywhere.
|
||
func poissonManufactured(t *testing.T, n int) (f, exact *core.Array) {
|
||
t.Helper()
|
||
flatF := make([]float64, n*n)
|
||
flatU := make([]float64, n*n)
|
||
for r := range n {
|
||
for c := range n {
|
||
x := float64(c) / float64(n-1)
|
||
y := float64(r) / float64(n-1)
|
||
i := r*n + c
|
||
flatU[i] = math.Sin(math.Pi*x) * math.Sin(math.Pi*y)
|
||
flatF[i] = 2 * math.Pi * math.Pi * flatU[i]
|
||
}
|
||
}
|
||
f = mustFloats(t, flatF, n, n)
|
||
exact = mustFloats(t, flatU, n, n)
|
||
return f, exact
|
||
}
|
||
|
||
// TestSolvePoissonDirichletManufactured checks the sine-solve against
|
||
// the manufactured solution: second-order convergence, with the
|
||
// coarse grid already inside 5e-3.
|
||
func TestSolvePoissonDirichletManufactured(t *testing.T) {
|
||
for _, n := range []int{9, 17, 33} {
|
||
f, exact := poissonManufactured(t, n)
|
||
u, err := SolvePoissonDirichlet(f, 1, 1)
|
||
if err != nil {
|
||
t.Fatalf("SolvePoissonDirichlet(%d): %v", n, err)
|
||
}
|
||
worst := 0.0
|
||
for i := range exact.Len() {
|
||
if e := math.Abs(u.FloatAt(i) - exact.FloatAt(i)); e > worst {
|
||
worst = e
|
||
}
|
||
}
|
||
h := 1 / float64(n-1)
|
||
if worst > 3*h*h {
|
||
t.Fatalf("grid %d: worst error %.3g above the O(h²) budget %.3g", n, worst, 3*h*h)
|
||
}
|
||
// The boundary is exactly zero by construction.
|
||
for c := range n {
|
||
if u.FloatAt(c) != 0 || u.FloatAt((n-1)*n+c) != 0 {
|
||
t.Fatalf("grid %d: boundary came back nonzero", n)
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestSolvePoissonNeumannManufactured checks the cosine solve with a
|
||
// zero-mean source whose solution is u = cos(πx)·cos(πy), the
|
||
// derivative-free data the Neumann solve exists for.
|
||
func TestSolvePoissonNeumannManufactured(t *testing.T) {
|
||
const n = 17
|
||
flatF := make([]float64, n*n)
|
||
flatU := make([]float64, n*n)
|
||
for r := range n {
|
||
for c := range n {
|
||
x := float64(c) / float64(n-1)
|
||
y := float64(r) / float64(n-1)
|
||
i := r*n + c
|
||
flatU[i] = math.Cos(math.Pi*x) * math.Cos(math.Pi*y)
|
||
flatF[i] = 2 * math.Pi * math.Pi * flatU[i]
|
||
}
|
||
}
|
||
// Zero-mean shift: subtract the mean, which also shifts u by a
|
||
// constant the Neumann problem cannot see.
|
||
meanF := 0.0
|
||
for _, v := range flatF {
|
||
meanF += v
|
||
}
|
||
meanF /= float64(n * n)
|
||
for i := range flatF {
|
||
flatF[i] -= meanF
|
||
}
|
||
f := mustFloats(t, flatF, n, n)
|
||
exact := mustFloats(t, flatU, n, n)
|
||
u, err := SolvePoissonNeumann(f, 1, 1)
|
||
if err != nil {
|
||
t.Fatalf("SolvePoissonNeumann: %v", err)
|
||
}
|
||
// Compare up to the free constant: recentre both to zero mean.
|
||
got := 0.0
|
||
for i := range u.Len() {
|
||
got += u.FloatAt(i)
|
||
}
|
||
got /= float64(u.Len())
|
||
worst := 0.0
|
||
for i := range u.Len() {
|
||
if e := math.Abs(u.FloatAt(i) - got - exact.FloatAt(i)); e > worst {
|
||
worst = e
|
||
}
|
||
}
|
||
h := 1 / float64(n-1)
|
||
if worst > 3*h*h {
|
||
t.Fatalf("worst error %.3g above the O(h²) budget %.3g", worst, 3*h*h)
|
||
}
|
||
}
|
||
|
||
// TestSolvePoissonNeumannMeanRefusal checks the compatibility
|
||
// condition: a nonzero-mean source has no Neumann solution.
|
||
func TestSolvePoissonNeumannMeanRefusal(t *testing.T) {
|
||
f := mustFloats(t, []float64{
|
||
1, 1, 1,
|
||
1, 1, 1,
|
||
1, 1, 1,
|
||
}, 3, 3)
|
||
if _, err := SolvePoissonNeumann(f, 1, 1); err == nil {
|
||
t.Fatal("nonzero-mean source accepted")
|
||
}
|
||
if _, err := SolvePoissonDirichlet(mustFloats(t, []float64{1, 1}, 1, 2), 1, 1); err == nil {
|
||
t.Fatal("grid under 3×3 accepted")
|
||
}
|
||
if _, err := SolvePoissonDirichlet(f, 0, 1); err == nil {
|
||
t.Fatal("non-positive length accepted")
|
||
}
|
||
}
|