add test with non-equidistant nodes

ranocha · Jun 14, 2024 · 1a2a376 · 1a2a376
1 parent 914cab2
commit 1a2a376
Show file tree

Hide file tree

Showing 2 changed files with 29 additions and 8 deletions.
diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl
@@ -15,6 +15,7 @@ function SummationByPartsOperators.function_space_operator(basis_functions,
                                                            options = Options(g_tol = 1e-14, iterations = 10000)) where {T, SourceOfCoefficients}
 
     @assert derivative_order == 1 "Only first derivative operators are supported"
+    sort!(nodes)
     weights, D = construct_function_space_operator(basis_functions, x_min, x_max, nodes, source; options = options)
     return MatrixDerivativeOperator(x_min, x_max, nodes, weights, D, accuracy_order, source)
 end
@@ -131,7 +132,7 @@ function construct_function_space_operator(basis_functions, x_min, x_max, nodes,
     x0 = zeros(L + N)
     fg!(F, G, x) = optimization_function_and_grad!(F, G, x, p)
     result = optimize(Optim.only_fg!(fg!), x0, LBFGS(), options)
-    display(result)
+
     x = minimizer(result)
     sigma = x[1:L]
     rho = x[(L + 1):end]

diff --git a/test/matrix_operators_test.jl b/test/matrix_operators_test.jl
@@ -118,10 +118,9 @@ end
 if VERSION >= v"1.9"
   @testset "Function space operators" begin
     N = 5
-    nodes = collect(range(-1, 1, length=N))
     x_min = -1.0
     x_max = 1.0
-    x = [x_min, -0.5, 0.0, 0.5, x_max]
+    nodes = collect(range(x_min, x_max, length=N))
     source = GlaubitzNordströmÖffner2023()
     for compact in (true, false)
       show(IOContext(devnull, :compact=>compact), source)
@@ -132,20 +131,41 @@ if VERSION >= v"1.9"
     let basis_functions = [x -> x^i for i in 0:3]
       D = function_space_operator(basis_functions, x_min, x_max, nodes, source)
 
+      @test grid(D) ≈ nodes
       @test ≈(D * ones(N), zeros(N); atol = 1e-13)
-      @test D * x ≈ ones(N)
-      @test D * (x .^ 2) ≈ 2 * x
-      @test D * (x .^ 3) ≈ 3 * (x .^ 2)
+      @test D * nodes ≈ ones(N)
+      @test D * (nodes .^ 2) ≈ 2 * nodes
+      @test D * (nodes .^ 3) ≈ 3 * (nodes .^ 2)
       M = mass_matrix(D)
       @test M * D.D + D.D' * M ≈ B
     end
 
     let basis_functions = [one, identity, exp]
       D = function_space_operator(basis_functions, x_min, x_max, nodes, source)
 
+      @test grid(D) ≈ nodes
       @test ≈(D * ones(N), zeros(N); atol = 1e-13)
-      @test D * x ≈ ones(N)
-      @test D * exp.(x) ≈ exp.(x)
+      @test D * nodes ≈ ones(N)
+      @test D * exp.(nodes) ≈ exp.(nodes)
+      M = mass_matrix(D)
+      @test M * D.D + D.D' * M ≈ B
+    end
+
+    # test non-equidistant nodes generated by `nodes = [0.0, rand(8)..., 1.0]`
+    nodes = [0.0, 0.01585580467018155, 0.18010381213204507, 0.270467434432868,
+             0.37699483985320303, 0.5600831197666554, 0.5698824835924449, 0.623949064816263,
+             0.8574665549914025, 1.0]
+    N = length(nodes)
+    B = zeros(N, N)
+    B[1, 1] = -1.0
+    B[N, N] = 1.0
+    let basis_functions = [one, identity, exp]
+      D = function_space_operator(basis_functions, first(nodes), last(nodes), nodes, source)
+
+      @test grid(D) ≈ nodes
+      @test ≈(D * ones(N), zeros(N); atol = 5e-13)
+      @test D * nodes ≈ ones(N)
+      @test D * exp.(nodes) ≈ exp.(nodes)
       M = mass_matrix(D)
       @test M * D.D + D.D' * M ≈ B
     end