From f820bccbd9039b2253181599ec29ef38d738f03b Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 18 Oct 2026 19:28:02 -0700 Subject: [PATCH 01/17] Add an individual Octave testGradient --- tests/matlab/CMakeLists.txt | 25 +++ tests/matlab/testAddScalarBC.m | 310 +++++++++++++++++++++++++++ tests/matlab/testBCConsistency.m | 78 +++++++ tests/matlab/testDivPolynomial.m | 145 +++++++++++++ tests/matlab/testDivergence.m | 26 +++ tests/matlab/testEnergy.m | 43 ++++ tests/matlab/testGI13Curv.m | 159 ++++++++++++++ tests/matlab/testGradPolynomial.m | 143 ++++++++++++ tests/matlab/testGradient.m | 24 +++ tests/matlab/testGridOrientation.m | 127 +++++++++++ tests/matlab/testGridOrientation3D.m | 95 ++++++++ tests/matlab/testLapPolynomial.m | 134 ++++++++++++ tests/matlab/testLaplacian.m | 24 +++ tests/matlab/testPoissonAccuracy.m | 56 +++++ tests/matlab/testTTMGrids.m | 116 ++++++++++ tests/octave/Makefile | 2 + tests/octave/testGradient.m | 46 ++-- 17 files changed, 1530 insertions(+), 23 deletions(-) create mode 100644 tests/matlab/CMakeLists.txt create mode 100644 tests/matlab/testAddScalarBC.m create mode 100644 tests/matlab/testBCConsistency.m create mode 100644 tests/matlab/testDivPolynomial.m create mode 100644 tests/matlab/testDivergence.m create mode 100644 tests/matlab/testEnergy.m create mode 100644 tests/matlab/testGI13Curv.m create mode 100644 tests/matlab/testGradPolynomial.m create mode 100644 tests/matlab/testGradient.m create mode 100644 tests/matlab/testGridOrientation.m create mode 100644 tests/matlab/testGridOrientation3D.m create mode 100644 tests/matlab/testLapPolynomial.m create mode 100644 tests/matlab/testLaplacian.m create mode 100644 tests/matlab/testPoissonAccuracy.m create mode 100644 tests/matlab/testTTMGrids.m create mode 100644 tests/octave/Makefile diff --git a/tests/matlab/CMakeLists.txt b/tests/matlab/CMakeLists.txt new file mode 100644 index 000000000..bee31a4d0 --- /dev/null +++ b/tests/matlab/CMakeLists.txt @@ -0,0 +1,25 @@ +# Find all GNU Octave test files (test*.m) +file(GLOB TEST_SOURCES test*.m) + +# Enable testing +enable_testing() + +# Loop through all GNU Octave test files and register them as tests +foreach(TEST_SOURCE ${TEST_SOURCES}) + get_filename_component(TEST_NAME ${TEST_SOURCE} NAME_WE) + + add_test( + NAME ${TEST_NAME} + COMMAND matlab + -nojvm + -batch + "addpath('${CMAKE_CURRENT_SOURCE_DIR}'); addpath('${CMAKE_SOURCE_DIR}/src/octave'); assertSuccess(runtests('${TEST_NAME}'))" + ) +endforeach() + +# Custom target to run all GNU Octave tests using CTest +add_custom_target(run_octave_tests + COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure + COMMENT "Running all GNU Octave unit tests..." +) + diff --git a/tests/matlab/testAddScalarBC.m b/tests/matlab/testAddScalarBC.m new file mode 100644 index 000000000..5feeba89d --- /dev/null +++ b/tests/matlab/testAddScalarBC.m @@ -0,0 +1,310 @@ +classdef testAddScalarBC < matlab.unittest.TestCase + methods(Test) + function testBC1D(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')); + + k = 2; + dx = 0.1; + m = 8; + tol = 1e-10; + + A = rand(m + 2); + b = rand(m + 2, 1); + + v = [7; 8]; + + % First test case: no boundary conditions + + dc0 = [0; 0]; + nc0 = [0; 0]; + + [A0, b0] = addScalarBC1D(A, b, k, m, dx, dc0, nc0, v); + + testCase.verifyEqual(A0, A, ... + "Matrix A with no boundary conditions doesn't match expected result") + testCase.verifyEqual(b0, b, ... + "Vector b with no boundary conditions doesn't match expected result") + + % Second test case: with boundary conditions + + dc1 = [2; 3]; + nc1 = [4; 5]; + + [A1, b1] = addScalarBC1D(A, b, k, m, dx, dc1, nc1, v); + + A_ref = A; + b_ref = b; + + A_ref(1,:) = 0; + A_ref(end,:) = 0; + b_ref(1) = 0; + b_ref(end) = 0; + + G = grad(k, m, dx); + + Al = zeros(m + 2); + Al(1,1) = dc1(1); + Bl = zeros(m + 2, m + 1); + Bl(1,1) = -nc1(1); + + Ar = zeros(m + 2); + Ar(end,end) = dc1(2); + Br = zeros(m + 2, m + 1); + Br(end,end) = nc1(2); + + A_ref = A_ref + (Al + Bl*G) + (Ar + Br*G); + + b_ref(1) = v(1); + b_ref(end) = v(2); + + testCase.verifyEqual(A1, A_ref, ... + "Matrix A with boundary conditions doesn't match expected result", ... + 'AbsTol', tol); + testCase.verifyEqual(b1, b_ref, ... + "Vector b with boundary conditions doesn't match expected result", ... + 'AbsTol', tol); + end + + function testBC2D(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')); + + k = 2; + dx = 0.1; + dy = 0.1; + m = 8; + n = m + 2; + tol = 1e-10; + + A = rand((m+2)*(n+2)); + b = rand((m+2)*(n+2),1); + + bclr = 4*ones(n,1); + bcbt = 5*ones(m+2,1); + + v = {bclr; bclr; bcbt; bcbt}; + + % First test case: no boundary conditions + + dc0 = [0; 0; 0; 0]; + nc0 = [0; 0; 0; 0]; + + [A0, b0] = addScalarBC2D(A, b, k, m, dx, n, dy, dc0, nc0, v); + + testCase.verifyEqual(A0, A, ... + "Matrix A with no boundary conditions doesn't match expected result") + testCase.verifyEqual(b0, b, ... + "Vector b with no boundary conditions doesn't match expected result") + + % Second test case: with boundary conditions + + A = rand((m+2)*(n+2)); + b = rand((m+2)*(n+2), 1); + + dc1 = [2; 3; 4; 5]; + nc1 = [6; 7; 8; 9]; + + [A1, b1] = addScalarBC2D(A, b, k, m, dx, n, dy, dc1, nc1, v); + + A_ref = A; + b_ref = b; + + Gx = grad(k, m, dx); + In = eye(n + 2); + In(1,1) = 0; + In(end,end) = 0; + + Al = zeros(m + 2, m + 2); + Al(1,1) = dc1(1); + Bl = zeros(m + 2, m + 1); + Bl(1,1) = -nc1(1); + + Ar = zeros(m + 2, m + 2); + Ar(end,end) = dc1(2); + Br = zeros(m + 2, m + 1); + Br(end,end) = nc1(2); + + Al = Al + Bl*Gx; + Ar = Ar + Br*Gx; + Abclr = kron(In, Al + Ar); + [bcrowslr,~,~] = find(Abclr); + + Gy = grad(k, n, dy); + Im = eye(m + 2); + + Ab = zeros(n + 2, n + 2); + Ab(1,1) = dc1(3); + Bb = zeros(n + 2, n + 1); + Bb(1,1) = -nc1(3); + + At = zeros(n + 2, n + 2); + At(end,end) = dc1(4); + Bt = zeros(n + 2, n + 1); + Bt(end,end) = nc1(4); + + Ab = Ab + Bb*Gy; + At = At + Bt*Gy; + Abcbt = kron(Ab + At, Im); + [bcrowsbt,~,~] = find(Abcbt); + + bcrows = unique([bcrowslr; bcrowsbt]); + A_ref(bcrows, :) = zeros(numel(bcrows), size(A_ref, 2)); + + A_ref = A_ref + Abclr + Abcbt; + + b_ref = reshape(b_ref, m+2, n+2); + b_ref(1,2:n+1) = bclr; + b_ref(end,2:n+1) = bclr; + b_ref(1:m+2,1) = bcbt; + b_ref(1:m+2,end) = bcbt; + b_ref = reshape(b_ref, [], 1); + + testCase.verifyEqual(A1, A_ref, ... + "Matrix A with boundary conditions doesn't match expected result", ... + 'AbsTol', tol); + testCase.verifyEqual(b1, b_ref, ... + "Vector b with boundary conditions doesn't match expected result", ... + 'AbsTol', tol); + end + + function testBC3D(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')); + + k = 2; + dx = 0.1; + dy = 0.1; + dz = 0.1; + m = 8; + n = m + 2; + o = n + 2; + tol = 1e-10; + + A = rand((m+2)*(n+2)*(o+2)); + b = rand((m+2)*(n+2)*(o+2),1); + + bclr = 4*ones(n*o,1); + bcbt = 5*ones((m+2)*o,1); + bcfz = 6*ones((m+2)*(n+2),1); + + v = {bclr; bclr; bcbt; bcbt; bcfz; bcfz}; + + % First test case: no boundary conditions + + dc0 = [0; 0; 0; 0; 0; 0]; + nc0 = [0; 0; 0; 0; 0; 0]; + + [A0, b0] = addScalarBC3D(A, b, k, m, dx, n, dy, o, dz, dc0, nc0, v); + + testCase.verifyEqual(A0, A, ... + "Matrix A with no boundary conditions doesn't match expected result") + testCase.verifyEqual(b0, b, ... + "Vector b with no boundary conditions doesn't match expected result") + + % Second test case: with boundary conditions + + A = rand((m+2)*(n+2)*(o+2)); + b = rand((m+2)*(n+2)*(o+2), 1); + + dc1 = [2; 3; 4; 5; 6; 7]; + nc1 = [8; 9; 10; 11; 12; 13]; + + [A1, b1] = addScalarBC3D(A, b, k, m, dx, n, dy, o, dz, dc1, nc1, v); + + A_ref = A; + b_ref = b; + + Gx = grad(k, m, dx); + Gy = grad(k, n, dy); + Gz = grad(k, o, dz); + + Im = eye(m + 2); + + In = eye(n + 2); + In(1,1) = 0; + In(end,end) = 0; + + Io = eye(o + 2); + Io(1,1) = 0; + Io(end,end) = 0; + + Al = zeros(m+2, m+2); + Al(1,1) = dc1(1); + Bl = zeros(m + 2, m + 1); + Bl(1,1) = -nc1(1); + + Ar = zeros(m+2, m+2); + Ar(end,end) = dc1(2); + Br = zeros(m + 2, m + 1); + Br(end,end) = nc1(2); + + Al = Al + Bl*Gx; + Ar = Ar + Br*Gx; + + Abclr = kron(kron(Io, In), Al+Ar); + [bcrowslr, ~, ~] = find(Abclr); + + Ab = zeros(n + 2, n + 2); + Ab(1,1) = dc1(3); + Bb = zeros(n + 2, n + 1); + Bb(1,1) = -nc1(3); + + At = zeros(n + 2, n + 2); + At(end,end) = dc1(4); + Bt = zeros(n + 2, n + 1); + Bt(end,end) = nc1(4); + + Ab = Ab + Bb*Gy; + At = At + Bt*Gy; + + Abcbt = kron(kron(Io, Ab+At), Im); + [bcrowsbt, ~, ~] = find(Abcbt); + + Af = zeros(o + 2, o + 2); + Af(1,1) = dc1(5); + Bf = zeros(o + 2, o + 1); + Bf(1,1) = -nc1(5); + + Az = zeros(o + 2, o + 2); + Az(end,end) = dc1(6); + Bz = zeros(o + 2, o + 1); + Bz(end,end) = nc1(6); + + Af = Af + Bf*Gz; + Az = Az + Bz*Gz; + In(1,1) = 1; + In(end,end) = 1; + Abcfz = kron(kron(Af+Az, In), Im); + [bcrowsfz, ~, ~] = find(Abcfz); + + bcrows = unique([bcrowslr; bcrowsbt; bcrowsfz]); + A_ref(bcrows, :) = zeros(numel(bcrows), size(A_ref, 2)); + + A_ref = A_ref + Abclr + Abcbt + Abcfz; + + b_ref(1:(m+2)*(n+2)) = bcfz; + b_ref((end-(m+2)*(n+2)+1):end) = bcfz; + + b_ref = reshape(b_ref, m+2, n+2, o+2); + for i=1:o + b_ref(1,2:n+1,i+1) = bclr(((i-1)*n+1):i*n); + b_ref(end,2:n+1,i+1) = bclr(((i-1)*n+1):i*n); + + b_ref(1:m+2,1,i+1) = bcbt(((i-1)*(m+2)+1):i*(m+2)); + b_ref(1:m+2,end,i+1) = bcbt(((i-1)*(m+2)+1):i*(m+2)); + end + b_ref = reshape(b_ref, [], 1); + + testCase.verifyEqual(A1, A_ref, ... + "Matrix A with boundary conditions doesn't match expected result", ... + 'AbsTol', tol); + testCase.verifyEqual(b1, b_ref, ... + "Vector b with boundary conditions doesn't match expected result", ... + 'AbsTol', tol); + end + end +end \ No newline at end of file diff --git a/tests/matlab/testBCConsistency.m b/tests/matlab/testBCConsistency.m new file mode 100644 index 000000000..75fbb8600 --- /dev/null +++ b/tests/matlab/testBCConsistency.m @@ -0,0 +1,78 @@ +classdef testBCConsistency < matlab.unittest.TestCase + methods(Test) + function test1DCase(testCase) + % 1D Boundary Conditions Consistency Test + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')); + + k = 4; % Order of accuracy + m = 50; % Number of cells + dx = 0.314; % Grid spacing + + % Dirichlet-Dirichlet Boundary Conditions + RBC = robinBC(k, m, dx, 1, 0); + MBC = mixedBC(k, m, dx, 'Dirichlet', 1, 'Dirichlet', 1); + testCase.verifyEqual(RBC(1), MBC(1), "1D Dirichlet BCs do not match at west boundary"); + testCase.verifyEqual(RBC(end), MBC(end), "1D Dirichlet BCs do not match at east boundary"); + + % Neumann-Neumann Boundary Conditions + RBC = robinBC(k, m, dx, 0, 2); + MBC = mixedBC(k, m, dx, 'Neumann', 2, 'Neumann', 2); + testCase.verifyEqual(RBC(1, :), MBC(1, :), "1D Neumann BCs do not match at west boundary"); + testCase.verifyEqual(RBC(end, :), MBC(end, :), "1D Neumann BCs do not match at east boundary"); + + % Robin-Robin Boundary Conditions + RBC = robinBC(k, m, dx, 3, 4); + MBC = mixedBC(k, m, dx, 'Robin', [3, 4], 'Robin', [3, 4]); + testCase.verifyEqual(RBC(1, :), MBC(1, :), "1D Robin BCs do not match at west boundary"); + testCase.verifyEqual(RBC(end, :), MBC(end, :), "1D Robin BCs do not match at east boundary"); + + % Mixed Robin BCs + RBC2 = robinBC(k, m, dx, 5, 6); + MBC = mixedBC(k, m, dx, 'Robin', [3, 4], 'Robin', [5, 6]); + testCase.verifyEqual(RBC(1, :), MBC(1, :), "1D Mixed Robin BCs do not match at west boundary"); + testCase.verifyEqual(RBC2(end, :), MBC(end, :), "1D Mixed Robin BCs do not match at east boundary"); + end + + function test2DCase(testCase) + % 2D Boundary Conditions Consistency Test + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')); + + k = 4; % Order of accuracy + m = 50; % Number of cells (x-direction) + n = 73; % Number of cells (y-direction) + dx = 0.314; % x-direction grid spacing + dy = 0.123; % y-direction grid spacing + + % Robin-Robin Boundary Conditions + RBC = robinBC2D(k, m, dx, n, dy, 3, 4); + MBC = mixedBC2D(k, m, dx, n, dy, 'Robin', [3, 4], 'Robin', [3, 4], 'Robin', [3, 4], 'Robin', [3, 4]); + testCase.verifyEqual(RBC, MBC, "2D Robin BCs are not consistent"); + end + + function test3DCase(testCase) + % 3D Boundary Conditions Consistency Test + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')); + + k = 4; % Order of accuracy + m = 50; % Number of cells (x-direction) + n = 73; % Number of cells (y-direction) + o = 18; % Number of cells (z-direction) + dx = 0.314; % x-direction grid spacing + dy = 0.123; % y-direction grid spacing + dz = 0.198; % z-direction grid spacing + + % Robin-Robin Boundary Conditions + RBC = robinBC3D(k, m, dx, n, dy, o, dz, 6, 7); + MBC = mixedBC3D(k, m, dx, n, dy, o, dz, 'Robin', [6, 7], 'Robin', [6, 7], ... + 'Robin', [6, 7], 'Robin', [6, 7], 'Robin', [6, 7], 'Robin', [6, 7]); + testCase.verifyEqual(RBC, MBC, "3D Robin BCs are not consistent"); + end + end +end + diff --git a/tests/matlab/testDivPolynomial.m b/tests/matlab/testDivPolynomial.m new file mode 100644 index 000000000..0264ba886 --- /dev/null +++ b/tests/matlab/testDivPolynomial.m @@ -0,0 +1,145 @@ +classdef testDivPolynomial < matlab.unittest.TestCase + methods (Test) + function testDivOfPolynomial1D(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + ks = [2, 4, 6, 8]; % Different orders of accuracy + tol = 1e-11; + x1 = 0; + x2 = 1; + + for k = ks + m = 2 * k + 1; + dx = (x2 - x1) / m; + + D = div(k, m, dx); + + sfield = [x1, x1+dx/2:dx:x2-dx/2, x2]'; + vfield = (x1:dx:x2)'; + + approx = D * vfield.^k; + analytic = k * sfield.^(k-1); + + % Note: ignore boundary points for divergence + err = approx(2:end-1) - analytic(2:end-1); + + testCase.verifyLessThan(norm(err), tol, ... + sprintf("Nullity test failed for k = %d", k)); + end + end + + function testDivOfPolynomial2D(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + ks = [2, 4, 6, 8]; % Different orders of accuracy + tol = 1e-11; + x1 = 0; + x2 = 1; + y1 = 0; + y2 = 1; + + for k = ks + m = 2 * k + 1; + dx = (x2 - x1) / m; + n = m + 2; + dy = (y2 - y1) / n; + + D = div2D(k, m, dx, n, dy); + + % staggered grid components + nodes_x = (x1:dx:x2)'; + nodes_y = (y1:dy:y2)'; + centers_x = (x1+dx/2:dx:x2-dx/2)'; + centers_y = (y1+dy/2:dy:y2-dy/2)'; + + % set up vector field + [~, vfield_x] = meshgrid(centers_y, nodes_x); + vfield_x = reshape(vfield_x, [], 1); + [vfield_y, ~] = meshgrid(nodes_y, centers_x); + vfield_y = reshape(vfield_y, [], 1); + + % set up scalar field + [sfield_y, sfield_x] = meshgrid([y1,centers_y',y2]', [x1,centers_x',x2]'); + + % compute approximate and analytical solutions + approx = D * [vfield_x; vfield_y].^k; + analytic = k * sfield_x.^(k-1) + k * sfield_y.^(k-1); + + % remove boundary for divergence; int(erior) only + analytic_int = zeros(m+2, n+2); + analytic_int(2:end-1, 2:end-1) = analytic(2:end-1, 2:end-1); + analytic_int = reshape(analytic_int, [], 1); + + err = approx - analytic_int; + + testCase.verifyLessThan(norm(err), tol, ... + sprintf("Nullity test failed for k = %d", k)); + end + end + + function testDivOfPolynomial3D(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + ks = [2, 4, 6, 8]; % Different orders of accuracy + tol = 1e-11; + x1 = 0; + x2 = 1; + y1 = 0; + y2 = 1; + z1 = 0; + z2 = 1; + + for k = ks + m = 2 * k + 1; + dx = (x2 - x1) / m; + n = m + 2; + dy = (y2 - y1) / n; + o = n + 2; + dz = (z2 - z1) / o; + + D = div3D(k, m, dx, n, dy, o, dz); + + % staggered grid components + nodes_x = (x1:dx:x2)'; + nodes_y = (y1:dy:y2)'; + nodes_z = (z1:dz:z2)'; + centers_x = (x1+dx/2:dx:x2-dx/2)'; + centers_y = (y1+dy/2:dy:y2-dy/2)'; + centers_z = (z1+dz/2:dz:z2-dz/2)'; + + % set up vector field + [~, vfield_x, ~] = meshgrid(centers_y, nodes_x, centers_z); + vfield_x = reshape(vfield_x, [], 1); + [vfield_y, ~, ~] = meshgrid(nodes_y, centers_x, centers_z); + vfield_y = reshape(vfield_y, [], 1); + [~, ~, vfield_z] = meshgrid(centers_y, centers_x, nodes_z); + vfield_z = reshape(vfield_z, [], 1); + + % set up scalar field + [sfield_y, sfield_x, sfield_z] = meshgrid([y1,centers_y',y2]', ... + [x1,centers_x',x2]', [z1,centers_z',z2]'); + + % compute approximate and analytical solutions + approx = D * [vfield_x; vfield_y; vfield_z].^k; + analytic = k * sfield_x.^(k-1) + k * sfield_y.^(k-1) + ... + k * sfield_z.^(k-1); + + % remove boundary for divergence; int(erior) only + analytic_int = zeros(m+2, n+2, o+2); + analytic_int(2:end-1, 2:end-1, 2:end-1) = analytic(2:end-1, 2:end-1, 2:end-1); + analytic_int = reshape(analytic_int, [], 1); + + err = approx - analytic_int; + + testCase.verifyLessThan(norm(err), tol, ... + sprintf("Nullity test failed for k = %d", k)); + end + end + end +end \ No newline at end of file diff --git a/tests/matlab/testDivergence.m b/tests/matlab/testDivergence.m new file mode 100644 index 000000000..85f1cb3e6 --- /dev/null +++ b/tests/matlab/testDivergence.m @@ -0,0 +1,26 @@ +classdef testDivergence < matlab.unittest.TestCase + methods (Test) + function testNullityOfDivergence(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + ks = [2, 4, 6, 8]; % Different orders of accuracy + tol = 1e-10; + + for k = ks + m = 2 * k + 1; + dx = 1 / m; + + D = div(k, m, dx); + + field = ones(m + 1, 1); + + sol = D * field; + + testCase.verifyLessThan(norm(sol), tol, ... + sprintf("Nullity test failed for k = %d", k)); + end + end + end +end diff --git a/tests/matlab/testEnergy.m b/tests/matlab/testEnergy.m new file mode 100644 index 000000000..25b04d3f8 --- /dev/null +++ b/tests/matlab/testEnergy.m @@ -0,0 +1,43 @@ +classdef testEnergy < matlab.unittest.TestCase + methods(Test) + function testforEnergy(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + %parameteres + k = 4; + a = -5; + b = 5; + m = 500; + grid = linspace(a, b, m); + dx = grid(2) - grid(1); + tol = 1e-6; + + % Laplacian + L = lap(k, m - 2, dx); + + % Transform grid using an anonymous function + grid = arrayfun(@(x) x^2, grid); + + % Potential matrix V + V = spdiags(grid', 0, m, m); % Transpose grid to match dimension + + % Hamiltonian H + H = -0.5 * L + V; + + % Eigenvalue computation + eigval = eig(full(H)); % Convert sparse matrix to full for eig function + eigval = sort(eigval); + + % Expected eigenvalues + expected = [1, 3, 5, 7, 9]; + + for i= 1:length(expected) + testCase.verifyLessThan(abs(real(eigval(i) / eigval(1)) - expected(i)),tol, ... + sprintf("Nullity test failed for k = %d", k)); + end + end + end +end + diff --git a/tests/matlab/testGI13Curv.m b/tests/matlab/testGI13Curv.m new file mode 100644 index 000000000..de3c59c32 --- /dev/null +++ b/tests/matlab/testGI13Curv.m @@ -0,0 +1,159 @@ +classdef testGI13Curv < matlab.unittest.TestCase +% Regression tests for GI13, the face-to-face interpolation used by +% grad3DCurvLegacy (which the four-argument grad3DCurv dispatches to). +% +% The defect these guard against was an index-map error: GI13 built +% kron(speye(n*o), I1), walking the source in n*o blocks when the source has +% (n+1)*o blocks, so the map drifted by one eta-plane per zeta-level. +% +% NOTE ON TEST CHOICE. The broken operator had the RIGHT SHAPE and mapped a +% constant to a constant, so shape checks and row-sum checks both passed on it. +% Only an explicit index map (testIndexMap) or a convergence study +% (testGradConvergesOnCurvilinearGrid) detects it. Please keep both. + + methods(Static) + function [X, Y, Z, xc, yc, zc] = sineGrid(N, amp) + % Nodes, and the cell centres with a ghost layer, for a smoothly + % distorted grid. Matches examples/matlab_octave conventions. + [a, b, c] = meshgrid(linspace(0,1,N), linspace(0,1,N), linspace(0,1,N)); + X = a + amp*sin(2*pi*b); + Y = b + amp*sin(2*pi*a); + Z = c + 0.5*amp*sin(2*pi*a).*sin(2*pi*b); + s = linspace(0,1,N); + s = [0, s(1:end-1)+0.5/(N-1), 1]; + [p, q, t] = meshgrid(s, s, s); + xc = p + amp*sin(2*pi*q); + yc = q + amp*sin(2*pi*p); + zc = t + 0.5*amp*sin(2*pi*p).*sin(2*pi*q); + end + end + + methods(Test) + + function testIndexMap(testCase) + % Feed GI13 a field whose values encode their own (i,j,k) and + % check every output reads from the correct zeta-plane. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + m = 4; + n = 3; + o = 3; + [I, J, K] = ndgrid(1:m, 1:n+1, 1:o); + v = I(:) + 100*J(:) + 10000*K(:); + W = reshape(GI13(speye(m*(n+1)*o), m, n, o, 'Gn')*v, m+1, n, o); + + bad = 0; + for k = 1:o + for j = 1:n + for i = 1:m+1 + val = W(i,j,k); + if val == 0 + continue; + end + kk = floor(val/10000); + if kk ~= k + bad = bad + 1; + end + end + end + end + testCase.verifyEqual(bad, 0, ... + sprintf('%d of %d GI13 outputs read from the wrong zeta-plane', ... + bad, (m+1)*n*o)); + end + + function testShapesAndConstants(testCase) + % All six shifts must have the right shape and map a constant + % field to a constant field. Weak, but cheap. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + m = 5; + n = 4; + o = 3; + ty = {'Gn','Gc','Ge','Gcy','Gee','Gnn'}; + src = [m*(n+1)*o, m*n*(o+1), (m+1)*n*o, m*n*(o+1), (m+1)*n*o, m*(n+1)*o]; + tgt = [(m+1)*n*o, (m+1)*n*o, m*(n+1)*o, m*(n+1)*o, m*n*(o+1), m*n*(o+1)]; + + for i = 1:numel(ty) + A = GI13(speye(src(i)), m, n, o, ty{i}); + testCase.verifyEqual(size(A), [tgt(i) src(i)], ... + sprintf('GI13 %s has the wrong shape', ty{i})); + testCase.verifyLessThan(max(abs(A*ones(src(i),1) - 1)), 1e-12, ... + sprintf('GI13 %s does not preserve a constant', ty{i})); + end + end + + function testGradExactOnCartesianGrid(testCase) + % On an undistorted grid the composition must be exact to roundoff. + % Guards against a "fix" that trades curvilinear accuracy for + % Cartesian accuracy. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + for N = [11 17 25] + [X, Y, Z, xc, yc, zc] = testGI13Curv.sineGrid(N, 0.0); + f = reshape(permute(xc.^2 + yc.^2 + zc.^2, [2 1 3]), [], 1); + L = div3DCurv(2, X, Y, Z)*grad3DCurv(2, X, Y, Z); + v = reshape(L*f, N+1, N+1, N+1); + d = v(3:end-2, 3:end-2, 3:end-2); + d = d(:) - 6; + testCase.verifyLessThan(sqrt(mean(d.^2)), 1e-9, ... + sprintf('Cartesian L=D*G not exact at n=%d', N)); + end + end + + function testGradExactOnLinearField(testCase) + % The gradient of f = x + 2y + 3z is a constant. The stock + % operator gave rms 0.39 at distortion 0.10; anything near that + % means the face shift is wrong again. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + N = 21; + for amp = [0.05 0.10] + [X, Y, Z, xc, yc, zc] = testGI13Curv.sineGrid(N, amp); + f = reshape(permute(xc + 2*yc + 3*zc, [2 1 3]), [], 1); + T = grad3DCurv(2, X, Y, Z)*f; + gx = T(1:N*(N-1)*(N-1)) - 1; + testCase.verifyLessThan(sqrt(mean(gx.^2)), 0.05, ... + sprintf('grad of a linear field is badly wrong at amp=%g', amp)); + end + end + + function testGradConvergesOnCurvilinearGrid(testCase) + % The headline check. The stock operator gave order 0.03-0.09 -- + % no convergence at all. Require better than first order, which + % the stock operator misses by a wide margin and the corrected one + % clears at ~2.1. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + amp = 0.10; + Ns = [13 21 33]; + errs = zeros(size(Ns)); + for i = 1:numel(Ns) + N = Ns(i); + [X, Y, Z, xc, yc, zc] = testGI13Curv.sineGrid(N, amp); + f = reshape(permute(xc.^2 + yc.^2 + zc.^2, [2 1 3]), [], 1); + T = grad3DCurv(2, X, Y, Z)*f; + cx = linspace(0,1,N); + cx = cx(1:end-1) + 0.5/(N-1); + [aa, bb, ~] = meshgrid(linspace(0,1,N), cx, cx); + xf = aa + amp*sin(2*pi*bb); + gx = T(1:N*(N-1)*(N-1)); + errs(i) = sqrt(mean((gx - reshape(permute(2*xf,[2 1 3]),[],1)).^2)); + end + p = log(errs(1)/errs(end)) / log((Ns(end)-1)/(Ns(1)-1)); + testCase.verifyGreaterThan(p, 1.5, ... + sprintf('grad3DCurv observed order %.2f on a curvilinear grid', p)); + end + + end +end diff --git a/tests/matlab/testGradPolynomial.m b/tests/matlab/testGradPolynomial.m new file mode 100644 index 000000000..a7068058b --- /dev/null +++ b/tests/matlab/testGradPolynomial.m @@ -0,0 +1,143 @@ +classdef testGradPolynomial < matlab.unittest.TestCase + methods (Test) + function testGradOfPolynomial1D(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + ks = [2, 4, 6, 8]; % Different orders of accuracy + tol = 1e-11; + x1 = 0; + x2 = 1; + + for k = ks + m = 2 * k + 1; + dx = (x2-x1) / m; + + G = grad(k, m, dx); + + sfield = [x1, x1+dx/2:dx:x2-dx/2, x2]'; + vfield = (x1:dx:x2)'; + + approx = G * sfield.^k; + analytic = k * vfield.^(k-1); + + err = approx - analytic; + + testCase.verifyLessThan(norm(err), tol, ... + sprintf("Nullity test failed for k = %d", k)); + end + end + + function testGradOfPolynomial2D(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + ks = [2, 4, 6, 8]; % Different orders of accuracy + tol = 1e-11; + x1 = 0; + x2 = 1; + y1 = 0; + y2 = 1; + + for k = ks + m = 2 * k + 1; + dx = (x2 - x1) / m; + n = m + 2; + dy = (y2 - y1) / n; + + G = grad2D(k, m, dx, n, dy); + + % staggered grid components + nodes_x = (x1:dx:x2)'; + nodes_y = (y1:dy:y2)'; + centers_x = (x1+dx/2:dx:x2-dx/2)'; + centers_y = (y1+dy/2:dy:y2-dy/2)'; + + % set up vector field + [~, vfield_x] = meshgrid(centers_y, nodes_x); + vfield_x = reshape(vfield_x, [], 1); + [vfield_y, ~] = meshgrid(nodes_y, centers_x); + vfield_y = reshape(vfield_y, [], 1); + + % set up scalar field + [sfield_y, sfield_x] = meshgrid([y1,centers_y',y2]', [x1,centers_x',x2]'); + sfield_x = reshape(sfield_x, [], 1); + sfield_y = reshape(sfield_y, [], 1); + + % compute approximate and analytical solutions + approx = G * (sfield_x.^k + sfield_y.^k); + analytic_x = k * vfield_x.^(k-1); + analytic_y = k * vfield_y.^(k-1); + analytic_combined = [analytic_x; analytic_y]; + + err = approx - analytic_combined; + + testCase.verifyLessThan(norm(err), tol, ... + sprintf("Nullity test failed for k = %d", k)); + end + end + + function testGradOfPolynomial3D(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + ks = [2, 4, 6, 8]; % Different orders of accuracy + tol = 1e-11; + x1 = 0; + x2 = 1; + y1 = 0; + y2 = 1; + z1 = 0; + z2 = 1; + + for k = ks + m = 2 * k + 1; + dx = (x2 - x1) / m; + n = m + 2; + dy = (y2 - y1) / n; + o = n + 2; + dz = (z2 - z1) / o; + + G = grad3D(k, m, dx, n, dy, o, dz); + + % staggered grid components + nodes_x = (x1:dx:x2)'; + nodes_y = (y1:dy:y2)'; + nodes_z = (z1:dz:z2)'; + centers_x = (x1+dx/2:dx:x2-dx/2)'; + centers_y = (y1+dy/2:dy:y2-dy/2)'; + centers_z = (z1+dz/2:dz:z2-dz/2)'; + + % set up vector field + [~, vfield_x, ~] = meshgrid(centers_y, nodes_x, centers_z); + vfield_x = reshape(vfield_x, [], 1); + [vfield_y, ~, ~] = meshgrid(nodes_y, centers_x, centers_z); + vfield_y = reshape(vfield_y, [], 1); + [~, ~, vfield_z] = meshgrid(centers_y, centers_x, nodes_z); + vfield_z = reshape(vfield_z, [], 1); + + % set up scalar field + [sfield_y, sfield_x, sfield_z] = meshgrid([y1,centers_y',y2]', ... + [x1,centers_x',x2]', [z1,centers_z',z2]'); + sfield_x = reshape(sfield_x, [], 1); + sfield_y = reshape(sfield_y, [], 1); + sfield_z = reshape(sfield_z, [], 1); + + % compute approximate and analytical solutions + approx = G * (sfield_x.^k + sfield_y.^k + sfield_z.^k); + analytic_x = k * vfield_x.^(k-1); + analytic_y = k * vfield_y.^(k-1); + analytic_z = k * vfield_z.^(k-1); + analytic_combined = [analytic_x; analytic_y; analytic_z]; + + err = approx - analytic_combined; + + testCase.verifyLessThan(norm(err), tol, ... + sprintf("Nullity test failed for k = %d", k)); + end + end + end +end \ No newline at end of file diff --git a/tests/matlab/testGradient.m b/tests/matlab/testGradient.m new file mode 100644 index 000000000..f80561823 --- /dev/null +++ b/tests/matlab/testGradient.m @@ -0,0 +1,24 @@ +classdef testGradient < matlab.unittest.TestCase + methods(Test) + function testNullityofGradient(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + ks=[2,4,6,8]; + tol = 1e-10; + for k = ks + m = 2*k+1; + dx=1/m; + + G= grad(k,m,dx); + field = ones(m+2,1); + sol = G* field; + + testCase.verifyLessThan(norm(sol),tol, ... + sprintf("Nullity test failed for k = %d", k)); + end + end + end +end + diff --git a/tests/matlab/testGridOrientation.m b/tests/matlab/testGridOrientation.m new file mode 100644 index 000000000..2ccf787a7 --- /dev/null +++ b/tests/matlab/testGridOrientation.m @@ -0,0 +1,127 @@ +classdef testGridOrientation < matlab.unittest.TestCase +% Tests for the left-handed / tangled grid detection added to jacobian2D. +% +% The curvilinear operators read [n, m] = size(X): rows are eta, columns are +% xi. A grid supplied the other way round is still self-consistent for grad +% and div -- J flips sign, the cofactors flip with it, gradients still come +% out exact -- so nothing fails and nothing used to complain. These tests +% check that it complains now, and, just as importantly, that it stays quiet +% on correctly oriented grids. + + methods(Static) + function [X, Y] = sineGrid(N, amp) + [X, Y] = meshgrid(linspace(0,1,N), linspace(0,1,N)); + X = X + amp*sin(2*pi*Y); + Y = Y + amp*sin(2*pi*X); + end + end + + methods(Test) + + function testRightHandedIsSilent(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + [X, Y] = testGridOrientation.sineGrid(21, 0.10); + testCase.verifyWarningFree(@() jacobian2D(2, X, Y)); + end + + function testTransposedWarns(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + [X, Y] = testGridOrientation.sineGrid(21, 0.10); + testCase.verifyWarning(@() jacobian2D(2, X.', Y.'), ... + 'jacobian2D:leftHandedGrid'); + end + + function testGridGenOutputWarns(testCase) + % The reported case: gridGen returns xi-nodes as ROWS, which is the + % opposite of what the operators read, so passing its output + % straight in yields a globally left-handed grid. + % + % tfi does addpath(['grids/' grid_name]) RELATIVE TO THE CURRENT + % DIRECTORY, so the working directory has to be src/octave, not + % src/octave/grids. And the path is added non-recursively on + % purpose: genpath would put every grid folder on the path at once, + % and the curve functions (right.m, top.m, ...) share names across + % grids, so another grid would silently answer for this one. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + origDir = pwd; + cleanupDir = onCleanup(@() cd(origDir)); + srcDir = fullfile(pwd, '..', '..', 'src', 'octave'); + addpath(srcDir) + cd(srcDir) + + % chevron, not swan: swan is itself tangled (mixed-sign Jacobian + % either way round), which is the other warning. See + % testTangledGridWarns. + [X, Y] = gridGen('TFI', 'chevron', 39, 99, false); + testCase.verifyNotEmpty(strfind(which('right'), 'chevron'), ... + 'the chevron curves are not the ones on the path'); + testCase.verifyWarning(@() jacobian2D(2, X, Y), ... + 'jacobian2D:leftHandedGrid'); + testCase.verifyWarningFree(@() jacobian2D(2, X.', Y.')); + end + + function testMagnitudesAreIdentical(testCase) + % Documents the symptom that makes this hard to notice: same grid, + % same |J| to the last bit, opposite sign. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + [X, Y] = testGridOrientation.sineGrid(21, 0.10); + warning('off', 'jacobian2D:leftHandedGrid'); + restoreWarn = onCleanup(@() warning('on', 'jacobian2D:leftHandedGrid')); + a = sort(abs(nonzeros(jacobian2D(2, X, Y)))); + b = sort(abs(nonzeros(jacobian2D(2, X.', Y.')))); + testCase.verifyEqual(a, b, 'AbsTol', 1e-12); + end + + function testTangledGridWarns(testCase) + % A mesh that folds over itself has mixed-sign J. Different + % problem, different warning. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + [X, Y] = meshgrid(linspace(0,1,15), linspace(0,1,15)); + X(6:10, 6:10) = X(6:10, 6:10) - 0.45; + testCase.verifyWarning(@() jacobian2D(2, X, Y), ... + 'jacobian2D:tangledGrid'); + end + + function testWarningReachesOperatorUsers(testCase) + % grad2DCurv and div2DCurv both route through jacobian2D, so users + % who never call it directly still get told. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + [X, Y] = testGridOrientation.sineGrid(21, 0.10); + testCase.verifyWarning(@() grad2DCurv(2, X.', Y.'), ... + 'jacobian2D:leftHandedGrid'); + testCase.verifyWarning(@() div2DCurv(2, X.', Y.'), ... + 'jacobian2D:leftHandedGrid'); + end + + function testFiveArgPath(testCase) + % jacobian2D has two entry paths: 3 arguments dispatches to + % jacobian2DLegacy, 5 arguments computes inline. Both are checked. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + N = 21; + dc = [1;1;1;1]; + nc = [0;0;0;0]; + s = linspace(0,1,N); + s = [0, s(1:end-1)+0.5/(N-1), 1]; + [X, Y] = meshgrid(s, s); + X = X + 0.08*sin(2*pi*Y); + Y = Y + 0.08*sin(2*pi*X); + testCase.verifyWarningFree(@() jacobian2D(2, X, Y, dc, nc)); + testCase.verifyWarning(@() jacobian2D(2, X.', Y.', dc, nc), ... + 'jacobian2D:leftHandedGrid'); + end + + end +end diff --git a/tests/matlab/testGridOrientation3D.m b/tests/matlab/testGridOrientation3D.m new file mode 100644 index 000000000..5cd4c1de9 --- /dev/null +++ b/tests/matlab/testGridOrientation3D.m @@ -0,0 +1,95 @@ +classdef testGridOrientation3D < matlab.unittest.TestCase +% Extends the 2D orientation guard to jacobian3D, and adds detection of +% vanishing Jacobians in both dimensions. +% +% In 3D, swapping two axes is an odd permutation, so handedness inverts and J +% goes uniformly negative -- the same failure as the 2D transpose, with one +% more way to trigger it (mirroring a single axis does it too). +% +% A vanishing Jacobian is a different problem: the mapping is degenerate at +% that cell and is not invertible there, and grad/div divide by J. It is not +% caught by the sign checks, because zero is neither positive nor negative. + + methods(Static) + function [X, Y, Z] = sineGrid3(N, amp) + [a, b, c] = meshgrid(linspace(0,1,N), linspace(0,1,N), ... + linspace(0,1,N)); + X = a + amp*sin(2*pi*b); + Y = b + amp*sin(2*pi*a); + Z = c + 0.5*amp*sin(2*pi*a).*sin(2*pi*b); + end + end + + methods(Test) + + function testRightHanded3DIsSilent(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + [X, Y, Z] = testGridOrientation3D.sineGrid3(13, 0.08); + testCase.verifyWarningFree(@() jacobian3D(2, X, Y, Z)); + end + + function testSwappedAxesWarn(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + [X, Y, Z] = testGridOrientation3D.sineGrid3(13, 0.08); + p = @(A) permute(A, [2 1 3]); + testCase.verifyWarning(@() jacobian3D(2, p(X), p(Y), p(Z)), ... + 'jacobian3D:leftHandedGrid'); + end + + function testMirroredAxisWarns(testCase) + % Reversing one axis is also an odd permutation. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + [X, Y, Z] = testGridOrientation3D.sineGrid3(13, 0.08); + f = @(A) A(end:-1:1, :, :); + testCase.verifyWarning(@() jacobian3D(2, f(X), f(Y), f(Z)), ... + 'jacobian3D:leftHandedGrid'); + end + + function testWarningReaches3DOperatorUsers(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + [X, Y, Z] = testGridOrientation3D.sineGrid3(13, 0.08); + p = @(A) permute(A, [2 1 3]); + testCase.verifyWarning(@() grad3DCurv(2, p(X), p(Y), p(Z)), ... + 'jacobian3D:leftHandedGrid'); + testCase.verifyWarning(@() div3DCurv(2, p(X), p(Y), p(Z)), ... + 'jacobian3D:leftHandedGrid'); + end + + function testVanishingJacobianWarns(testCase) + % Three identical columns make the central difference across the + % middle one exactly zero, so those cells are degenerate. The sign + % checks cannot see this: zero is neither positive nor negative. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + N = 21; + [X, Y] = meshgrid(linspace(0,1,N), linspace(0,1,N)); + X = X + 0.10*sin(2*pi*Y); + Y = Y + 0.10*sin(2*pi*X); + X(:, 9:11) = repmat(X(:, 10), 1, 3); + Y(:, 9:11) = repmat(Y(:, 10), 1, 3); + testCase.verifyWarning(@() jacobian2D(2, X, Y), ... + 'jacobian2D:vanishingJacobian'); + end + + function testHealthyGridHasNoVanishingWarning(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + N = 21; + [X, Y] = meshgrid(linspace(0,1,N), linspace(0,1,N)); + X = X + 0.10*sin(2*pi*Y); + Y = Y + 0.10*sin(2*pi*X); + testCase.verifyWarningFree(@() jacobian2D(2, X, Y)); + end + + end +end diff --git a/tests/matlab/testLapPolynomial.m b/tests/matlab/testLapPolynomial.m new file mode 100644 index 000000000..68151bfbf --- /dev/null +++ b/tests/matlab/testLapPolynomial.m @@ -0,0 +1,134 @@ +classdef testLapPolynomial < matlab.unittest.TestCase + methods (Test) + function testLapOfPolynomial1D(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + ks = [2, 4, 6, 8]; % Different orders of accuracy + tol = 1e-10; + x1 = 0; + x2 = 1; + + for k = ks + m = 2 * k + 1; + dx = (x2-x1) / m; + + L = lap(k, m, dx); + + sfield = [x1, x1+dx/2:dx:x2-dx/2, x2]'; + + approx = L * sfield.^k; + analytic = k * (k-1) * sfield.^(k-2); + + % Note: ignore boundary points (because of divergence) + err = approx(2:end-1) - analytic(2:end-1); + + testCase.verifyLessThan(norm(err), tol, ... + sprintf("Nullity test failed for k = %d", k)); + end + end + + function testLapOfPolynomial2D(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + ks = [2, 4, 6, 8]; % Different orders of accuracy + tol = 1e-10; + x1 = 0; + x2 = 1; + y1 = 0; + y2 = 1; + + for k = ks + m = 2 * k + 1; + dx = (x2 - x1) / m; + n = m + 2; + dy = (y2 - y1) / n; + + L = lap2D(k, m, dx, n, dy); + + % staggered grid components (only centers needed) + centers_x = (x1+dx/2:dx:x2-dx/2)'; + centers_y = (y1+dy/2:dy:y2-dy/2)'; + + % set up scalar field only (vector field not needed) + [sfield_y, sfield_x] = meshgrid([y1,centers_y',y2]', [x1,centers_x',x2]'); + sfield_x = reshape(sfield_x, [], 1); + sfield_y = reshape(sfield_y, [], 1); + + % compute approximate and analytical solutions + approx = L * (sfield_x.^k + sfield_y.^k); + analytic = k * (k-1) * sfield_x.^(k-2) + ... + k * (k-1) * sfield_y.^(k-2); + + % remove boundary (because of divergence); int(erior) only + analytic = reshape(analytic, m+2, n+2); + analytic_int = zeros(m+2, n+2); + analytic_int(2:end-1, 2:end-1) = analytic(2:end-1, 2:end-1); + analytic_int = reshape(analytic_int, [], 1); + + err = approx - analytic_int; + + testCase.verifyLessThan(norm(err), tol, ... + sprintf("Nullity test failed for k = %d", k)); + end + end + + function testLapOfPolynomial3D(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + ks = [2, 4, 6, 8]; % Different orders of accuracy + tol = 1e-10; + x1 = 0; + x2 = 1; + y1 = 0; + y2 = 1; + z1 = 0; + z2 = 1; + + for k = ks + m = 2 * k + 1; + dx = (x2 - x1) / m; + n = m + 2; + dy = (y2 - y1) / n; + o = n + 2; + dz = (z2 - z1) / o; + + L = lap3D(k, m, dx, n, dy, o, dz); + + % staggered grid components (only centers needed) + centers_x = (x1+dx/2:dx:x2-dx/2)'; + centers_y = (y1+dy/2:dy:y2-dy/2)'; + centers_z = (z1+dz/2:dz:z2-dz/2)'; + + % set up scalar field only (vector field not needed) + [sfield_y, sfield_x, sfield_z] = meshgrid([y1,centers_y',y2]', ... + [x1,centers_x',x2]', [z1,centers_z',z2]'); + sfield_x = reshape(sfield_x, [], 1); + sfield_y = reshape(sfield_y, [], 1); + sfield_z = reshape(sfield_z, [], 1); + + % compute approximate and analytical solutions + approx = L * (sfield_x.^k + sfield_y.^k + sfield_z.^k); + analytic = k * (k-1) * sfield_x.^(k-2) + ... + k * (k-1) * sfield_y.^(k-2) + ... + k * (k-1) * sfield_z.^(k-2); + + % remove boundary (because of divergence); int(erior) only + analytic = reshape(analytic, m+2, n+2, o+2); + analytic_int = zeros(m+2, n+2, o+2); + analytic_int(2:end-1, 2:end-1, 2:end-1) = analytic(2:end-1, 2:end-1, 2:end-1); + analytic_int = reshape(analytic_int, [], 1); + + err = approx - analytic_int; + + testCase.verifyLessThan(norm(err), tol, ... + sprintf("Nullity test failed for k = %d", k)); + end + end + end +end \ No newline at end of file diff --git a/tests/matlab/testLaplacian.m b/tests/matlab/testLaplacian.m new file mode 100644 index 000000000..65b1a07c7 --- /dev/null +++ b/tests/matlab/testLaplacian.m @@ -0,0 +1,24 @@ +classdef testLaplacian < matlab.unittest.TestCase + methods(Test) + function testNullityofLaplacian(testCase) + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + addpath(genpath('../../src/octave')) + + ks=[2,4,6,8]; + tol = 1e-10; + for k = ks + m = 2*k+1; + dx=1/m; + + L= lap(k,m,dx); + field = ones(m+2,1); + sol = L * field; + + testCase.verifyLessThan(norm(sol),tol, ... + sprintf("Nullity test failed for k = %d", k)); + end + end + end +end + diff --git a/tests/matlab/testPoissonAccuracy.m b/tests/matlab/testPoissonAccuracy.m new file mode 100644 index 000000000..5026236b2 --- /dev/null +++ b/tests/matlab/testPoissonAccuracy.m @@ -0,0 +1,56 @@ +classdef testPoissonAccuracy < matlab.unittest.TestCase + methods(Test) + function testforEnergy(testCase) + + west = 0; % Domain's limits + east = 1; + + ks = [2, 4, 6]; % Different orders of accuracy + % tests for ks = 8 removed for issues + grid_sizes = [20, 40]; % Grid sizes to test + + for k = ks + + errors = zeros(size(grid_sizes)); + + for i = 1:numel(grid_sizes) + m = grid_sizes(i); % Number of cells + dx = (east - west) / m; % Step length + + L = lap(k, m, dx); % 1D Mimetic Laplacian operator + + % Impose Robin BC on Laplacian operator + a = 1; + b = 1; + L = L + robinBC(k, m, dx, a, b); + + % 1D Staggered grid + grid = [west west+dx/2 : dx : east-dx/2 east]; + + % RHS + U = exp(grid)'; + U(1) = 0; % West BC + U(end) = 2*exp(1); % East BC + + % Solve a linear system of equations + computed_solution = L\U; + + % Compute error using L2 norm + analytical_solution = exp(grid); + errors(i) = max(abs(computed_solution' - analytical_solution)); + end + + % Compute order of accuracy + order = zeros(numel(errors) - 1, 1); + for i = 1:numel(errors) - 1 + order(i) = log2(errors(i) / errors(i + 1)); + + testCase.verifyGreaterThan(order(i),k-0.5, ... + sprintf("Test failed for k = %d", k)); + end + end + end + end +end + + diff --git a/tests/matlab/testTTMGrids.m b/tests/matlab/testTTMGrids.m new file mode 100644 index 000000000..f4a3bdc69 --- /dev/null +++ b/tests/matlab/testTTMGrids.m @@ -0,0 +1,116 @@ +classdef testTTMGrids < matlab.unittest.TestCase +% Tests that ttm produces unfolded grids. +% +% Before the TFI initial guess was added, ttm started with the interior at +% zeros(m,n) -- every interior node at the origin -- so the first sweeps +% computed alpha, beta and gamma from a maximally degenerate grid. The +% iteration then had to untangle a fold it had created itself, and on every +% grid shipped with MOLE it failed to: swan, chevron and horseshoe all came +% out with mixed-sign Jacobians at 500 iterations. +% +% A folded grid is locally left-handed. jacobian2D returns the signed +% determinant and grad2DCurv divides by it, so a run on one produces finite, +% plausible output with the metric tensor inverted in a patch. Nothing +% raises. +% +% NOTE ON ORIENTATION. ttm returns X with xi-nodes as rows, which is the +% opposite of what jacobian2D reads, so these tests transpose. Without the +% transpose a correct grid reports ALL NEGATIVE rather than ALL POSITIVE. + + methods(Static) + function [npos, nneg] = jacobianSigns(X, Y) + origState = warning('off', 'jacobian2D:leftHandedGrid'); + restore1 = onCleanup(@() warning(origState)); + origState2 = warning('off', 'jacobian2D:tangledGrid'); + restore2 = onCleanup(@() warning(origState2)); + d = nonzeros(jacobian2D(2, X.', Y.')); + npos = sum(d > 0); + nneg = sum(d < 0); + end + end + + methods(Test) + + function testShippedGridsAreNotFolded(testCase) + % chevron and horseshoe are clean under TFI, so there is no excuse + % for the elliptic generator to fold them. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + origDir = pwd; + cleanupDir = onCleanup(@() cd(origDir)); + addpath(genpath('../../src/octave')) + cd('../../src/octave') + + for g = {'chevron', 'horseshoe'} + [X, Y] = ttm(g{1}, 39, 99, 500, false); + [~, nneg] = testTTMGrids.jacobianSigns(X, Y); + testCase.verifyEqual(nneg, 0, ... + sprintf('ttm folded the %s grid: %d negative Jacobian entries', ... + g{1}, nneg)); + end + end + + function testEllipticSmoothingRepairsAFoldedTFI(testCase) + % swan is folded under TFI itself (113 negative entries). Repairing + % that is exactly what an elliptic generator is for, and with a + % sensible initial guess ttm does it. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + origDir = pwd; + cleanupDir = onCleanup(@() cd(origDir)); + addpath(genpath('../../src/octave')) + cd('../../src/octave') + + [Xt, Yt] = tfi('swan', 39, 99, false); + [~, nnegTFI] = testTTMGrids.jacobianSigns(Xt, Yt); + testCase.verifyGreaterThan(nnegTFI, 0, ... + 'swan is expected to be folded under TFI; test premise changed'); + + [X, Y] = ttm('swan', 39, 99, 500, false); + [~, nneg] = testTTMGrids.jacobianSigns(X, Y); + testCase.verifyEqual(nneg, 0, ... + sprintf('ttm left swan folded: %d negative entries', nneg)); + end + + function testFewIterationsDoNotFold(testCase) + % The failure was not slow degradation from over-iteration -- the + % stock version was folded within 5 sweeps. Check the low end. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + origDir = pwd; + cleanupDir = onCleanup(@() cd(origDir)); + addpath(genpath('../../src/octave')) + cd('../../src/octave') + + for it = [5 25 100] + [X, Y] = ttm('chevron', 39, 99, it, false); + [~, nneg] = testTTMGrids.jacobianSigns(X, Y); + testCase.verifyEqual(nneg, 0, ... + sprintf('ttm folded chevron after %d iterations (%d negative)', ... + it, nneg)); + end + end + + function testBoundariesAreUnchanged(testCase) + % The initial guess must only touch the interior. Boundary nodes + % come from the curve definitions and must survive untouched. + origPath = path; + cleanupObj = onCleanup(@() path(origPath)); + origDir = pwd; + cleanupDir = onCleanup(@() cd(origDir)); + addpath(genpath('../../src/octave')) + cd('../../src/octave') + + m = 39; + n = 99; + [X, Y] = ttm('chevron', m, n, 50, false); + [Xt, Yt] = tfi('chevron', m, n, false); + + testCase.verifyEqual(X(1, :), Xt(1, :), 'AbsTol', 1e-12); + testCase.verifyEqual(X(end, :), Xt(end, :), 'AbsTol', 1e-12); + testCase.verifyEqual(Y(:, 1), Yt(:, 1), 'AbsTol', 1e-12); + testCase.verifyEqual(Y(:, end), Yt(:, end), 'AbsTol', 1e-12); + end + + end +end diff --git a/tests/octave/Makefile b/tests/octave/Makefile new file mode 100644 index 000000000..7d9d41f9f --- /dev/null +++ b/tests/octave/Makefile @@ -0,0 +1,2 @@ +run_octave_tests: + octave --quiet --eval "addpath(pwd()); test('testGradient')" diff --git a/tests/octave/testGradient.m b/tests/octave/testGradient.m index f80561823..8774f0778 100644 --- a/tests/octave/testGradient.m +++ b/tests/octave/testGradient.m @@ -1,24 +1,24 @@ -classdef testGradient < matlab.unittest.TestCase - methods(Test) - function testNullityofGradient(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - ks=[2,4,6,8]; - tol = 1e-10; - for k = ks - m = 2*k+1; - dx=1/m; - - G= grad(k,m,dx); - field = ones(m+2,1); - sol = G* field; - - testCase.verifyLessThan(norm(sol),tol, ... - sprintf("Nullity test failed for k = %d", k)); - end - end - end -end +function testGradient () +endfunction + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! +%! addpath (genpath ("../../src/octave")); +%! +%! ks = [2, 4, 6, 8]; +%! tol = 1e-10; +%! +%! for k = ks +%! m = 2*k + 1; +%! dx = 1/m; +%! +%! G = grad (k, m, dx); +%! field = ones (m + 2, 1); +%! sol = G * field; +%! +%! assert (norm (sol) < tol, ... +%! sprintf ("Nullity test failed for k = %d", k)); +%! endfor From 832883ee95fef7989ec96ec5d407b2b9137e9868 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 18 Oct 2026 19:54:56 -0700 Subject: [PATCH 02/17] Add an individual Octave testDivergence --- tests/octave/Makefile | 2 ++ tests/octave/testDivergence.m | 46 +++++++++++++++++------------------ 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/tests/octave/Makefile b/tests/octave/Makefile index 7d9d41f9f..57fd75a7a 100644 --- a/tests/octave/Makefile +++ b/tests/octave/Makefile @@ -1,2 +1,4 @@ run_octave_tests: octave --quiet --eval "addpath(pwd()); test('testGradient')" + octave --quiet --eval "addpath(pwd()); test('testDivergence')" + diff --git a/tests/octave/testDivergence.m b/tests/octave/testDivergence.m index 85f1cb3e6..ca21eeb1c 100644 --- a/tests/octave/testDivergence.m +++ b/tests/octave/testDivergence.m @@ -1,26 +1,24 @@ -classdef testDivergence < matlab.unittest.TestCase - methods (Test) - function testNullityOfDivergence(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - ks = [2, 4, 6, 8]; % Different orders of accuracy - tol = 1e-10; +function testDivergence () +endfunction - for k = ks - m = 2 * k + 1; - dx = 1 / m; - D = div(k, m, dx); - - field = ones(m + 1, 1); - - sol = D * field; - - testCase.verifyLessThan(norm(sol), tol, ... - sprintf("Nullity test failed for k = %d", k)); - end - end - end -end +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! +%! addpath (genpath ("../../src/octave")); +%! +%! ks = [2, 4, 6, 8]; % Different orders of accuracy +%! tol = 1e-10; +%! +%! for k = ks +%! m = 2 * k + 1; +%! dx = 1 / m; +%! +%! D = div (k, m, dx); +%! field = ones (m + 1, 1); +%! sol = D * field; +%! +%! assert (norm (sol) < tol, ... +%! sprintf ("Nullity test failed for k = %d", k)); +%! endfor From 63deba04aae178317db732a0e8097b659ac88b89 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 18 Oct 2026 21:05:06 -0700 Subject: [PATCH 03/17] Add an individual Octave testLaplacian --- tests/octave/Makefile | 1 + tests/octave/testLaplacian.m | 46 ++++++++++++++++++------------------ 2 files changed, 24 insertions(+), 23 deletions(-) diff --git a/tests/octave/Makefile b/tests/octave/Makefile index 57fd75a7a..bf1be4f42 100644 --- a/tests/octave/Makefile +++ b/tests/octave/Makefile @@ -1,4 +1,5 @@ run_octave_tests: octave --quiet --eval "addpath(pwd()); test('testGradient')" octave --quiet --eval "addpath(pwd()); test('testDivergence')" + octave --quiet --eval "addpath(pwd()); test('testLaplacian')" diff --git a/tests/octave/testLaplacian.m b/tests/octave/testLaplacian.m index 65b1a07c7..933a5671b 100644 --- a/tests/octave/testLaplacian.m +++ b/tests/octave/testLaplacian.m @@ -1,24 +1,24 @@ -classdef testLaplacian < matlab.unittest.TestCase - methods(Test) - function testNullityofLaplacian(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - ks=[2,4,6,8]; - tol = 1e-10; - for k = ks - m = 2*k+1; - dx=1/m; - - L= lap(k,m,dx); - field = ones(m+2,1); - sol = L * field; - - testCase.verifyLessThan(norm(sol),tol, ... - sprintf("Nullity test failed for k = %d", k)); - end - end - end -end +function testLaplacian () +endfunction + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! +%! addpath (genpath ("../../src/octave")); +%! +%! ks = [2, 4, 6, 8]; +%! tol = 1e-10; +%! +%! for k = ks +%! m = 2 * k + 1; +%! dx = 1 / m; +%! +%! L = lap (k, m, dx); +%! field = ones (m + 2, 1); +%! sol = L * field; +%! +%! assert (norm (sol) < tol, ... +%! sprintf ("Nullity test failed for k = %d", k)); +%! endfor From d40d111de296ec84850443ec7e1f14abc76e39d6 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 18 Oct 2026 21:13:28 -0700 Subject: [PATCH 04/17] Add an individual Octave testAddScalarBC --- tests/octave/Makefile | 1 + tests/octave/testAddScalarBC.m | 609 ++++++++++++++++----------------- 2 files changed, 300 insertions(+), 310 deletions(-) diff --git a/tests/octave/Makefile b/tests/octave/Makefile index bf1be4f42..639b45ba6 100644 --- a/tests/octave/Makefile +++ b/tests/octave/Makefile @@ -2,4 +2,5 @@ run_octave_tests: octave --quiet --eval "addpath(pwd()); test('testGradient')" octave --quiet --eval "addpath(pwd()); test('testDivergence')" octave --quiet --eval "addpath(pwd()); test('testLaplacian')" + octave --quiet --eval "addpath(pwd()); test('testAddScalarBC')" diff --git a/tests/octave/testAddScalarBC.m b/tests/octave/testAddScalarBC.m index 5feeba89d..41f246f81 100644 --- a/tests/octave/testAddScalarBC.m +++ b/tests/octave/testAddScalarBC.m @@ -1,310 +1,299 @@ -classdef testAddScalarBC < matlab.unittest.TestCase - methods(Test) - function testBC1D(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')); - - k = 2; - dx = 0.1; - m = 8; - tol = 1e-10; - - A = rand(m + 2); - b = rand(m + 2, 1); - - v = [7; 8]; - - % First test case: no boundary conditions - - dc0 = [0; 0]; - nc0 = [0; 0]; - - [A0, b0] = addScalarBC1D(A, b, k, m, dx, dc0, nc0, v); - - testCase.verifyEqual(A0, A, ... - "Matrix A with no boundary conditions doesn't match expected result") - testCase.verifyEqual(b0, b, ... - "Vector b with no boundary conditions doesn't match expected result") - - % Second test case: with boundary conditions - - dc1 = [2; 3]; - nc1 = [4; 5]; - - [A1, b1] = addScalarBC1D(A, b, k, m, dx, dc1, nc1, v); - - A_ref = A; - b_ref = b; - - A_ref(1,:) = 0; - A_ref(end,:) = 0; - b_ref(1) = 0; - b_ref(end) = 0; - - G = grad(k, m, dx); - - Al = zeros(m + 2); - Al(1,1) = dc1(1); - Bl = zeros(m + 2, m + 1); - Bl(1,1) = -nc1(1); - - Ar = zeros(m + 2); - Ar(end,end) = dc1(2); - Br = zeros(m + 2, m + 1); - Br(end,end) = nc1(2); - - A_ref = A_ref + (Al + Bl*G) + (Ar + Br*G); - - b_ref(1) = v(1); - b_ref(end) = v(2); - - testCase.verifyEqual(A1, A_ref, ... - "Matrix A with boundary conditions doesn't match expected result", ... - 'AbsTol', tol); - testCase.verifyEqual(b1, b_ref, ... - "Vector b with boundary conditions doesn't match expected result", ... - 'AbsTol', tol); - end - - function testBC2D(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')); - - k = 2; - dx = 0.1; - dy = 0.1; - m = 8; - n = m + 2; - tol = 1e-10; - - A = rand((m+2)*(n+2)); - b = rand((m+2)*(n+2),1); - - bclr = 4*ones(n,1); - bcbt = 5*ones(m+2,1); - - v = {bclr; bclr; bcbt; bcbt}; - - % First test case: no boundary conditions - - dc0 = [0; 0; 0; 0]; - nc0 = [0; 0; 0; 0]; - - [A0, b0] = addScalarBC2D(A, b, k, m, dx, n, dy, dc0, nc0, v); - - testCase.verifyEqual(A0, A, ... - "Matrix A with no boundary conditions doesn't match expected result") - testCase.verifyEqual(b0, b, ... - "Vector b with no boundary conditions doesn't match expected result") - - % Second test case: with boundary conditions - - A = rand((m+2)*(n+2)); - b = rand((m+2)*(n+2), 1); - - dc1 = [2; 3; 4; 5]; - nc1 = [6; 7; 8; 9]; - - [A1, b1] = addScalarBC2D(A, b, k, m, dx, n, dy, dc1, nc1, v); - - A_ref = A; - b_ref = b; - - Gx = grad(k, m, dx); - In = eye(n + 2); - In(1,1) = 0; - In(end,end) = 0; - - Al = zeros(m + 2, m + 2); - Al(1,1) = dc1(1); - Bl = zeros(m + 2, m + 1); - Bl(1,1) = -nc1(1); - - Ar = zeros(m + 2, m + 2); - Ar(end,end) = dc1(2); - Br = zeros(m + 2, m + 1); - Br(end,end) = nc1(2); - - Al = Al + Bl*Gx; - Ar = Ar + Br*Gx; - Abclr = kron(In, Al + Ar); - [bcrowslr,~,~] = find(Abclr); - - Gy = grad(k, n, dy); - Im = eye(m + 2); - - Ab = zeros(n + 2, n + 2); - Ab(1,1) = dc1(3); - Bb = zeros(n + 2, n + 1); - Bb(1,1) = -nc1(3); - - At = zeros(n + 2, n + 2); - At(end,end) = dc1(4); - Bt = zeros(n + 2, n + 1); - Bt(end,end) = nc1(4); - - Ab = Ab + Bb*Gy; - At = At + Bt*Gy; - Abcbt = kron(Ab + At, Im); - [bcrowsbt,~,~] = find(Abcbt); - - bcrows = unique([bcrowslr; bcrowsbt]); - A_ref(bcrows, :) = zeros(numel(bcrows), size(A_ref, 2)); - - A_ref = A_ref + Abclr + Abcbt; - - b_ref = reshape(b_ref, m+2, n+2); - b_ref(1,2:n+1) = bclr; - b_ref(end,2:n+1) = bclr; - b_ref(1:m+2,1) = bcbt; - b_ref(1:m+2,end) = bcbt; - b_ref = reshape(b_ref, [], 1); - - testCase.verifyEqual(A1, A_ref, ... - "Matrix A with boundary conditions doesn't match expected result", ... - 'AbsTol', tol); - testCase.verifyEqual(b1, b_ref, ... - "Vector b with boundary conditions doesn't match expected result", ... - 'AbsTol', tol); - end - - function testBC3D(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')); - - k = 2; - dx = 0.1; - dy = 0.1; - dz = 0.1; - m = 8; - n = m + 2; - o = n + 2; - tol = 1e-10; - - A = rand((m+2)*(n+2)*(o+2)); - b = rand((m+2)*(n+2)*(o+2),1); - - bclr = 4*ones(n*o,1); - bcbt = 5*ones((m+2)*o,1); - bcfz = 6*ones((m+2)*(n+2),1); - - v = {bclr; bclr; bcbt; bcbt; bcfz; bcfz}; - - % First test case: no boundary conditions - - dc0 = [0; 0; 0; 0; 0; 0]; - nc0 = [0; 0; 0; 0; 0; 0]; - - [A0, b0] = addScalarBC3D(A, b, k, m, dx, n, dy, o, dz, dc0, nc0, v); - - testCase.verifyEqual(A0, A, ... - "Matrix A with no boundary conditions doesn't match expected result") - testCase.verifyEqual(b0, b, ... - "Vector b with no boundary conditions doesn't match expected result") - - % Second test case: with boundary conditions - - A = rand((m+2)*(n+2)*(o+2)); - b = rand((m+2)*(n+2)*(o+2), 1); - - dc1 = [2; 3; 4; 5; 6; 7]; - nc1 = [8; 9; 10; 11; 12; 13]; - - [A1, b1] = addScalarBC3D(A, b, k, m, dx, n, dy, o, dz, dc1, nc1, v); - - A_ref = A; - b_ref = b; - - Gx = grad(k, m, dx); - Gy = grad(k, n, dy); - Gz = grad(k, o, dz); - - Im = eye(m + 2); - - In = eye(n + 2); - In(1,1) = 0; - In(end,end) = 0; - - Io = eye(o + 2); - Io(1,1) = 0; - Io(end,end) = 0; - - Al = zeros(m+2, m+2); - Al(1,1) = dc1(1); - Bl = zeros(m + 2, m + 1); - Bl(1,1) = -nc1(1); - - Ar = zeros(m+2, m+2); - Ar(end,end) = dc1(2); - Br = zeros(m + 2, m + 1); - Br(end,end) = nc1(2); - - Al = Al + Bl*Gx; - Ar = Ar + Br*Gx; - - Abclr = kron(kron(Io, In), Al+Ar); - [bcrowslr, ~, ~] = find(Abclr); - - Ab = zeros(n + 2, n + 2); - Ab(1,1) = dc1(3); - Bb = zeros(n + 2, n + 1); - Bb(1,1) = -nc1(3); - - At = zeros(n + 2, n + 2); - At(end,end) = dc1(4); - Bt = zeros(n + 2, n + 1); - Bt(end,end) = nc1(4); - - Ab = Ab + Bb*Gy; - At = At + Bt*Gy; - - Abcbt = kron(kron(Io, Ab+At), Im); - [bcrowsbt, ~, ~] = find(Abcbt); - - Af = zeros(o + 2, o + 2); - Af(1,1) = dc1(5); - Bf = zeros(o + 2, o + 1); - Bf(1,1) = -nc1(5); - - Az = zeros(o + 2, o + 2); - Az(end,end) = dc1(6); - Bz = zeros(o + 2, o + 1); - Bz(end,end) = nc1(6); - - Af = Af + Bf*Gz; - Az = Az + Bz*Gz; - In(1,1) = 1; - In(end,end) = 1; - Abcfz = kron(kron(Af+Az, In), Im); - [bcrowsfz, ~, ~] = find(Abcfz); - - bcrows = unique([bcrowslr; bcrowsbt; bcrowsfz]); - A_ref(bcrows, :) = zeros(numel(bcrows), size(A_ref, 2)); - - A_ref = A_ref + Abclr + Abcbt + Abcfz; - - b_ref(1:(m+2)*(n+2)) = bcfz; - b_ref((end-(m+2)*(n+2)+1):end) = bcfz; - - b_ref = reshape(b_ref, m+2, n+2, o+2); - for i=1:o - b_ref(1,2:n+1,i+1) = bclr(((i-1)*n+1):i*n); - b_ref(end,2:n+1,i+1) = bclr(((i-1)*n+1):i*n); - - b_ref(1:m+2,1,i+1) = bcbt(((i-1)*(m+2)+1):i*(m+2)); - b_ref(1:m+2,end,i+1) = bcbt(((i-1)*(m+2)+1):i*(m+2)); - end - b_ref = reshape(b_ref, [], 1); - - testCase.verifyEqual(A1, A_ref, ... - "Matrix A with boundary conditions doesn't match expected result", ... - 'AbsTol', tol); - testCase.verifyEqual(b1, b_ref, ... - "Vector b with boundary conditions doesn't match expected result", ... - 'AbsTol', tol); - end - end -end \ No newline at end of file +function testAddScalarBC () +endfunction + + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave")); +%! +%! k = 2; +%! dx = 0.1; +%! m = 8; +%! tol = 1e-10; +%! +%! A = rand (m + 2); +%! b = rand (m + 2, 1); +%! v = [7; 8]; +%! +%! % No boundary conditions +%! dc0 = [0; 0]; +%! nc0 = [0; 0]; +%! [A0, b0] = addScalarBC1D (A, b, k, m, dx, dc0, nc0, v); +%! +%! assert (isequal (A0, A), ... +%! "Matrix A with no boundary conditions does not match"); +%! assert (isequal (b0, b), ... +%! "Vector b with no boundary conditions does not match"); +%! +%! % Boundary conditions +%! dc1 = [2; 3]; +%! nc1 = [4; 5]; +%! [A1, b1] = addScalarBC1D (A, b, k, m, dx, dc1, nc1, v); +%! +%! A_ref = A; +%! b_ref = b; +%! A_ref (1, :) = 0; +%! A_ref (end, :) = 0; +%! b_ref (1) = 0; +%! b_ref (end) = 0; +%! +%! G = grad (k, m, dx); +%! +%! Al = zeros (m + 2); +%! Al (1, 1) = dc1 (1); +%! Bl = zeros (m + 2, m + 1); +%! Bl (1, 1) = -nc1 (1); +%! +%! Ar = zeros (m + 2); +%! Ar (end, end) = dc1 (2); +%! Br = zeros (m + 2, m + 1); +%! Br (end, end) = nc1 (2); +%! +%! A_ref = A_ref + (Al + Bl * G) + (Ar + Br * G); +%! b_ref (1) = v (1); +%! b_ref (end) = v (2); +%! +%! assert (all (abs (A1(:) - A_ref(:)) <= tol), ... +%! "Matrix A with boundary conditions does not match"); +%! assert (all (abs (b1(:) - b_ref(:)) <= tol), ... +%! "Vector b with boundary conditions does not match"); + + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave")); +%! +%! k = 2; +%! dx = 0.1; +%! dy = 0.1; +%! m = 8; +%! n = m + 2; +%! tol = 1e-10; +%! +%! A = rand ((m + 2) * (n + 2)); +%! b = rand ((m + 2) * (n + 2), 1); +%! +%! bclr = 4 * ones (n, 1); +%! bcbt = 5 * ones (m + 2, 1); +%! v = {bclr; bclr; bcbt; bcbt}; +%! +%! % No boundary conditions +%! dc0 = [0; 0; 0; 0]; +%! nc0 = [0; 0; 0; 0]; +%! +%! [A0, b0] = addScalarBC2D ( ... +%! A, b, k, m, dx, n, dy, dc0, nc0, v); +%! +%! assert (isequal (A0, A), ... +%! "Matrix A with no boundary conditions does not match"); +%! assert (isequal (b0, b), ... +%! "Vector b with no boundary conditions does not match"); +%! +%! % Boundary conditions +%! A = rand ((m + 2) * (n + 2)); +%! b = rand ((m + 2) * (n + 2), 1); +%! +%! dc1 = [2; 3; 4; 5]; +%! nc1 = [6; 7; 8; 9]; +%! +%! [A1, b1] = addScalarBC2D ( ... +%! A, b, k, m, dx, n, dy, dc1, nc1, v); +%! +%! A_ref = A; +%! b_ref = b; +%! +%! Gx = grad (k, m, dx); +%! In = eye (n + 2); +%! In (1, 1) = 0; +%! In (end, end) = 0; +%! +%! Al = zeros (m + 2, m + 2); +%! Al (1, 1) = dc1 (1); +%! Bl = zeros (m + 2, m + 1); +%! Bl (1, 1) = -nc1 (1); +%! +%! Ar = zeros (m + 2, m + 2); +%! Ar (end, end) = dc1 (2); +%! Br = zeros (m + 2, m + 1); +%! Br (end, end) = nc1 (2); +%! +%! Al = Al + Bl * Gx; +%! Ar = Ar + Br * Gx; +%! Abclr = kron (In, Al + Ar); +%! [bcrowslr, ~, ~] = find (Abclr); +%! +%! Gy = grad (k, n, dy); +%! Im = eye (m + 2); +%! +%! Ab = zeros (n + 2, n + 2); +%! Ab (1, 1) = dc1 (3); +%! Bb = zeros (n + 2, n + 1); +%! Bb (1, 1) = -nc1 (3); +%! +%! At = zeros (n + 2, n + 2); +%! At (end, end) = dc1 (4); +%! Bt = zeros (n + 2, n + 1); +%! Bt (end, end) = nc1 (4); +%! +%! Ab = Ab + Bb * Gy; +%! At = At + Bt * Gy; +%! Abcbt = kron (Ab + At, Im); +%! [bcrowsbt, ~, ~] = find (Abcbt); +%! +%! bcrows = unique ([bcrowslr; bcrowsbt]); +%! A_ref (bcrows, :) = 0; +%! A_ref = A_ref + Abclr + Abcbt; +%! +%! b_ref = reshape (b_ref, m + 2, n + 2); +%! b_ref (1, 2:n+1) = bclr; +%! b_ref (end, 2:n+1) = bclr; +%! b_ref (1:m+2, 1) = bcbt; +%! b_ref (1:m+2, end) = bcbt; +%! b_ref = reshape (b_ref, [], 1); +%! +%! assert (all (abs (A1(:) - A_ref(:)) <= tol), ... +%! "Matrix A with boundary conditions does not match"); +%! assert (all (abs (b1(:) - b_ref(:)) <= tol), ... +%! "Vector b with boundary conditions does not match"); + + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave")); +%! +%! k = 2; +%! dx = 0.1; +%! dy = 0.1; +%! dz = 0.1; +%! m = 8; +%! n = m + 2; +%! o = n + 2; +%! tol = 1e-10; +%! +%! A = rand ((m + 2) * (n + 2) * (o + 2)); +%! b = rand ((m + 2) * (n + 2) * (o + 2), 1); +%! +%! bclr = 4 * ones (n * o, 1); +%! bcbt = 5 * ones ((m + 2) * o, 1); +%! bcfz = 6 * ones ((m + 2) * (n + 2), 1); +%! +%! v = {bclr; bclr; bcbt; bcbt; bcfz; bcfz}; +%! +%! % No boundary conditions +%! dc0 = zeros (6, 1); +%! nc0 = zeros (6, 1); +%! +%! [A0, b0] = addScalarBC3D ( ... +%! A, b, k, m, dx, n, dy, o, dz, dc0, nc0, v); +%! +%! assert (isequal (A0, A), ... +%! "Matrix A with no boundary conditions does not match"); +%! assert (isequal (b0, b), ... +%! "Vector b with no boundary conditions does not match"); +%! +%! % Boundary conditions +%! A = rand ((m + 2) * (n + 2) * (o + 2)); +%! b = rand ((m + 2) * (n + 2) * (o + 2), 1); +%! +%! dc1 = [2; 3; 4; 5; 6; 7]; +%! nc1 = [8; 9; 10; 11; 12; 13]; +%! +%! [A1, b1] = addScalarBC3D ( ... +%! A, b, k, m, dx, n, dy, o, dz, dc1, nc1, v); +%! +%! A_ref = A; +%! b_ref = b; +%! +%! Gx = grad (k, m, dx); +%! Gy = grad (k, n, dy); +%! Gz = grad (k, o, dz); +%! +%! Im = eye (m + 2); +%! +%! In = eye (n + 2); +%! In (1, 1) = 0; +%! In (end, end) = 0; +%! +%! Io = eye (o + 2); +%! Io (1, 1) = 0; +%! Io (end, end) = 0; +%! +%! Al = zeros (m + 2, m + 2); +%! Al (1, 1) = dc1 (1); +%! Bl = zeros (m + 2, m + 1); +%! Bl (1, 1) = -nc1 (1); +%! +%! Ar = zeros (m + 2, m + 2); +%! Ar (end, end) = dc1 (2); +%! Br = zeros (m + 2, m + 1); +%! Br (end, end) = nc1 (2); +%! +%! Al = Al + Bl * Gx; +%! Ar = Ar + Br * Gx; +%! Abclr = kron (kron (Io, In), Al + Ar); +%! [bcrowslr, ~, ~] = find (Abclr); +%! +%! Ab = zeros (n + 2, n + 2); +%! Ab (1, 1) = dc1 (3); +%! Bb = zeros (n + 2, n + 1); +%! Bb (1, 1) = -nc1 (3); +%! +%! At = zeros (n + 2, n + 2); +%! At (end, end) = dc1 (4); +%! Bt = zeros (n + 2, n + 1); +%! Bt (end, end) = nc1 (4); +%! +%! Ab = Ab + Bb * Gy; +%! At = At + Bt * Gy; +%! Abcbt = kron (kron (Io, Ab + At), Im); +%! [bcrowsbt, ~, ~] = find (Abcbt); +%! +%! Af = zeros (o + 2, o + 2); +%! Af (1, 1) = dc1 (5); +%! Bf = zeros (o + 2, o + 1); +%! Bf (1, 1) = -nc1 (5); +%! +%! Az = zeros (o + 2, o + 2); +%! Az (end, end) = dc1 (6); +%! Bz = zeros (o + 2, o + 1); +%! Bz (end, end) = nc1 (6); +%! +%! Af = Af + Bf * Gz; +%! Az = Az + Bz * Gz; +%! +%! In (1, 1) = 1; +%! In (end, end) = 1; +%! +%! Abcfz = kron (kron (Af + Az, In), Im); +%! [bcrowsfz, ~, ~] = find (Abcfz); +%! +%! bcrows = unique ([bcrowslr; bcrowsbt; bcrowsfz]); +%! A_ref (bcrows, :) = 0; +%! A_ref = A_ref + Abclr + Abcbt + Abcfz; +%! +%! b_ref (1:(m + 2) * (n + 2)) = bcfz; +%! b_ref (end - (m + 2) * (n + 2) + 1:end) = bcfz; +%! +%! b_ref = reshape (b_ref, m + 2, n + 2, o + 2); +%! +%! for i = 1:o +%! b_ref (1, 2:n+1, i+1) = ... +%! bclr ((i-1) * n + 1:i * n); +%! b_ref (end, 2:n+1, i+1) = ... +%! bclr ((i-1) * n + 1:i * n); +%! +%! b_ref (1:m+2, 1, i+1) = ... +%! bcbt ((i-1) * (m+2) + 1:i * (m+2)); +%! b_ref (1:m+2, end, i+1) = ... +%! bcbt ((i-1) * (m+2) + 1:i * (m+2)); +%! endfor +%! +%! b_ref = reshape (b_ref, [], 1); +%! +%! assert (all (abs (A1(:) - A_ref(:)) <= tol), ... +%! "Matrix A with boundary conditions does not match"); +%! assert (all (abs (b1(:) - b_ref(:)) <= tol), ... +%! "Vector b with boundary conditions does not match"); From 6757e0ce67490421a89f80a8f6a01695165a72be Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 18 Oct 2026 21:22:02 -0700 Subject: [PATCH 05/17] Add an individual Octave testBCConsistency --- tests/octave/Makefile | 1 + tests/octave/testBCConsistency.m | 171 ++++++++++++++++++------------- 2 files changed, 100 insertions(+), 72 deletions(-) diff --git a/tests/octave/Makefile b/tests/octave/Makefile index 639b45ba6..119e52418 100644 --- a/tests/octave/Makefile +++ b/tests/octave/Makefile @@ -3,4 +3,5 @@ run_octave_tests: octave --quiet --eval "addpath(pwd()); test('testDivergence')" octave --quiet --eval "addpath(pwd()); test('testLaplacian')" octave --quiet --eval "addpath(pwd()); test('testAddScalarBC')" + octave --quiet --eval "addpath(pwd()); test('testBCConsistency')" diff --git a/tests/octave/testBCConsistency.m b/tests/octave/testBCConsistency.m index 75fbb8600..8a1314ae7 100644 --- a/tests/octave/testBCConsistency.m +++ b/tests/octave/testBCConsistency.m @@ -1,78 +1,105 @@ -classdef testBCConsistency < matlab.unittest.TestCase - methods(Test) - function test1DCase(testCase) - % 1D Boundary Conditions Consistency Test - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')); - - k = 4; % Order of accuracy - m = 50; % Number of cells - dx = 0.314; % Grid spacing +function testBCConsistency () +endfunction - % Dirichlet-Dirichlet Boundary Conditions - RBC = robinBC(k, m, dx, 1, 0); - MBC = mixedBC(k, m, dx, 'Dirichlet', 1, 'Dirichlet', 1); - testCase.verifyEqual(RBC(1), MBC(1), "1D Dirichlet BCs do not match at west boundary"); - testCase.verifyEqual(RBC(end), MBC(end), "1D Dirichlet BCs do not match at east boundary"); - % Neumann-Neumann Boundary Conditions - RBC = robinBC(k, m, dx, 0, 2); - MBC = mixedBC(k, m, dx, 'Neumann', 2, 'Neumann', 2); - testCase.verifyEqual(RBC(1, :), MBC(1, :), "1D Neumann BCs do not match at west boundary"); - testCase.verifyEqual(RBC(end, :), MBC(end, :), "1D Neumann BCs do not match at east boundary"); +%!test +%! % 1D Boundary Conditions Consistency Test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave")); +%! +%! k = 4; +%! m = 50; +%! dx = 0.314; +%! +%! % Dirichlet-Dirichlet Boundary Conditions +%! RBC = robinBC (k, m, dx, 1, 0); +%! MBC = mixedBC (k, m, dx, "Dirichlet", 1, "Dirichlet", 1); +%! +%! assert (isequal (RBC (1), MBC (1)), ... +%! "1D Dirichlet BCs do not match at west boundary"); +%! assert (isequal (RBC (end), MBC (end)), ... +%! "1D Dirichlet BCs do not match at east boundary"); +%! +%! % Neumann-Neumann Boundary Conditions +%! RBC = robinBC (k, m, dx, 0, 2); +%! MBC = mixedBC (k, m, dx, "Neumann", 2, "Neumann", 2); +%! +%! assert (isequal (RBC (1, :), MBC (1, :)), ... +%! "1D Neumann BCs do not match at west boundary"); +%! assert (isequal (RBC (end, :), MBC (end, :)), ... +%! "1D Neumann BCs do not match at east boundary"); +%! +%! % Robin-Robin Boundary Conditions +%! RBC = robinBC (k, m, dx, 3, 4); +%! MBC = mixedBC (k, m, dx, "Robin", [3, 4], ... +%! "Robin", [3, 4]); +%! +%! assert (isequal (RBC (1, :), MBC (1, :)), ... +%! "1D Robin BCs do not match at west boundary"); +%! assert (isequal (RBC (end, :), MBC (end, :)), ... +%! "1D Robin BCs do not match at east boundary"); +%! +%! % Mixed Robin Boundary Conditions +%! RBC2 = robinBC (k, m, dx, 5, 6); +%! MBC = mixedBC (k, m, dx, "Robin", [3, 4], ... +%! "Robin", [5, 6]); +%! +%! assert (isequal (RBC (1, :), MBC (1, :)), ... +%! "1D mixed Robin BCs do not match at west boundary"); +%! assert (isequal (RBC2 (end, :), MBC (end, :)), ... +%! "1D mixed Robin BCs do not match at east boundary"); - % Robin-Robin Boundary Conditions - RBC = robinBC(k, m, dx, 3, 4); - MBC = mixedBC(k, m, dx, 'Robin', [3, 4], 'Robin', [3, 4]); - testCase.verifyEqual(RBC(1, :), MBC(1, :), "1D Robin BCs do not match at west boundary"); - testCase.verifyEqual(RBC(end, :), MBC(end, :), "1D Robin BCs do not match at east boundary"); - % Mixed Robin BCs - RBC2 = robinBC(k, m, dx, 5, 6); - MBC = mixedBC(k, m, dx, 'Robin', [3, 4], 'Robin', [5, 6]); - testCase.verifyEqual(RBC(1, :), MBC(1, :), "1D Mixed Robin BCs do not match at west boundary"); - testCase.verifyEqual(RBC2(end, :), MBC(end, :), "1D Mixed Robin BCs do not match at east boundary"); - end +%!test +%! % 2D Boundary Conditions Consistency Test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave")); +%! +%! k = 4; +%! m = 50; +%! n = 73; +%! dx = 0.314; +%! dy = 0.123; +%! +%! RBC = robinBC2D (k, m, dx, n, dy, 3, 4); +%! +%! MBC = mixedBC2D ( ... +%! k, m, dx, n, dy, ... +%! "Robin", [3, 4], ... +%! "Robin", [3, 4], ... +%! "Robin", [3, 4], ... +%! "Robin", [3, 4]); +%! +%! assert (isequal (RBC, MBC), ... +%! "2D Robin BCs are not consistent"); - function test2DCase(testCase) - % 2D Boundary Conditions Consistency Test - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')); - - k = 4; % Order of accuracy - m = 50; % Number of cells (x-direction) - n = 73; % Number of cells (y-direction) - dx = 0.314; % x-direction grid spacing - dy = 0.123; % y-direction grid spacing - - % Robin-Robin Boundary Conditions - RBC = robinBC2D(k, m, dx, n, dy, 3, 4); - MBC = mixedBC2D(k, m, dx, n, dy, 'Robin', [3, 4], 'Robin', [3, 4], 'Robin', [3, 4], 'Robin', [3, 4]); - testCase.verifyEqual(RBC, MBC, "2D Robin BCs are not consistent"); - end - - function test3DCase(testCase) - % 3D Boundary Conditions Consistency Test - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')); - - k = 4; % Order of accuracy - m = 50; % Number of cells (x-direction) - n = 73; % Number of cells (y-direction) - o = 18; % Number of cells (z-direction) - dx = 0.314; % x-direction grid spacing - dy = 0.123; % y-direction grid spacing - dz = 0.198; % z-direction grid spacing - - % Robin-Robin Boundary Conditions - RBC = robinBC3D(k, m, dx, n, dy, o, dz, 6, 7); - MBC = mixedBC3D(k, m, dx, n, dy, o, dz, 'Robin', [6, 7], 'Robin', [6, 7], ... - 'Robin', [6, 7], 'Robin', [6, 7], 'Robin', [6, 7], 'Robin', [6, 7]); - testCase.verifyEqual(RBC, MBC, "3D Robin BCs are not consistent"); - end - end -end +%!test +%! % 3D Boundary Conditions Consistency Test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave")); +%! +%! k = 4; +%! m = 50; +%! n = 73; +%! o = 18; +%! dx = 0.314; +%! dy = 0.123; +%! dz = 0.198; +%! +%! RBC = robinBC3D (k, m, dx, n, dy, o, dz, 6, 7); +%! +%! MBC = mixedBC3D ( ... +%! k, m, dx, n, dy, o, dz, ... +%! "Robin", [6, 7], ... +%! "Robin", [6, 7], ... +%! "Robin", [6, 7], ... +%! "Robin", [6, 7], ... +%! "Robin", [6, 7], ... +%! "Robin", [6, 7]); +%! +%! assert (isequal (RBC, MBC), ... +%! "3D Robin BCs are not consistent"); From 4ea794cb32cf45d85a8e484b4ba5f65ebda985aa Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 18 Oct 2026 21:27:48 -0700 Subject: [PATCH 06/17] Add an individual Octave testDivPolynomial --- tests/octave/Makefile | 2 + tests/octave/testDivPolynomial.m | 307 ++++++++++++++++--------------- 2 files changed, 164 insertions(+), 145 deletions(-) diff --git a/tests/octave/Makefile b/tests/octave/Makefile index 119e52418..1a69618ba 100644 --- a/tests/octave/Makefile +++ b/tests/octave/Makefile @@ -4,4 +4,6 @@ run_octave_tests: octave --quiet --eval "addpath(pwd()); test('testLaplacian')" octave --quiet --eval "addpath(pwd()); test('testAddScalarBC')" octave --quiet --eval "addpath(pwd()); test('testBCConsistency')" + octave --quiet --eval "addpath(pwd()); test('testDivPolynomial')" + diff --git a/tests/octave/testDivPolynomial.m b/tests/octave/testDivPolynomial.m index 0264ba886..5169d87d3 100644 --- a/tests/octave/testDivPolynomial.m +++ b/tests/octave/testDivPolynomial.m @@ -1,145 +1,162 @@ -classdef testDivPolynomial < matlab.unittest.TestCase - methods (Test) - function testDivOfPolynomial1D(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - ks = [2, 4, 6, 8]; % Different orders of accuracy - tol = 1e-11; - x1 = 0; - x2 = 1; - - for k = ks - m = 2 * k + 1; - dx = (x2 - x1) / m; - - D = div(k, m, dx); - - sfield = [x1, x1+dx/2:dx:x2-dx/2, x2]'; - vfield = (x1:dx:x2)'; - - approx = D * vfield.^k; - analytic = k * sfield.^(k-1); - - % Note: ignore boundary points for divergence - err = approx(2:end-1) - analytic(2:end-1); - - testCase.verifyLessThan(norm(err), tol, ... - sprintf("Nullity test failed for k = %d", k)); - end - end - - function testDivOfPolynomial2D(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - ks = [2, 4, 6, 8]; % Different orders of accuracy - tol = 1e-11; - x1 = 0; - x2 = 1; - y1 = 0; - y2 = 1; - - for k = ks - m = 2 * k + 1; - dx = (x2 - x1) / m; - n = m + 2; - dy = (y2 - y1) / n; - - D = div2D(k, m, dx, n, dy); - - % staggered grid components - nodes_x = (x1:dx:x2)'; - nodes_y = (y1:dy:y2)'; - centers_x = (x1+dx/2:dx:x2-dx/2)'; - centers_y = (y1+dy/2:dy:y2-dy/2)'; - - % set up vector field - [~, vfield_x] = meshgrid(centers_y, nodes_x); - vfield_x = reshape(vfield_x, [], 1); - [vfield_y, ~] = meshgrid(nodes_y, centers_x); - vfield_y = reshape(vfield_y, [], 1); - - % set up scalar field - [sfield_y, sfield_x] = meshgrid([y1,centers_y',y2]', [x1,centers_x',x2]'); - - % compute approximate and analytical solutions - approx = D * [vfield_x; vfield_y].^k; - analytic = k * sfield_x.^(k-1) + k * sfield_y.^(k-1); - - % remove boundary for divergence; int(erior) only - analytic_int = zeros(m+2, n+2); - analytic_int(2:end-1, 2:end-1) = analytic(2:end-1, 2:end-1); - analytic_int = reshape(analytic_int, [], 1); - - err = approx - analytic_int; - - testCase.verifyLessThan(norm(err), tol, ... - sprintf("Nullity test failed for k = %d", k)); - end - end - - function testDivOfPolynomial3D(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - ks = [2, 4, 6, 8]; % Different orders of accuracy - tol = 1e-11; - x1 = 0; - x2 = 1; - y1 = 0; - y2 = 1; - z1 = 0; - z2 = 1; - - for k = ks - m = 2 * k + 1; - dx = (x2 - x1) / m; - n = m + 2; - dy = (y2 - y1) / n; - o = n + 2; - dz = (z2 - z1) / o; - - D = div3D(k, m, dx, n, dy, o, dz); - - % staggered grid components - nodes_x = (x1:dx:x2)'; - nodes_y = (y1:dy:y2)'; - nodes_z = (z1:dz:z2)'; - centers_x = (x1+dx/2:dx:x2-dx/2)'; - centers_y = (y1+dy/2:dy:y2-dy/2)'; - centers_z = (z1+dz/2:dz:z2-dz/2)'; - - % set up vector field - [~, vfield_x, ~] = meshgrid(centers_y, nodes_x, centers_z); - vfield_x = reshape(vfield_x, [], 1); - [vfield_y, ~, ~] = meshgrid(nodes_y, centers_x, centers_z); - vfield_y = reshape(vfield_y, [], 1); - [~, ~, vfield_z] = meshgrid(centers_y, centers_x, nodes_z); - vfield_z = reshape(vfield_z, [], 1); - - % set up scalar field - [sfield_y, sfield_x, sfield_z] = meshgrid([y1,centers_y',y2]', ... - [x1,centers_x',x2]', [z1,centers_z',z2]'); - - % compute approximate and analytical solutions - approx = D * [vfield_x; vfield_y; vfield_z].^k; - analytic = k * sfield_x.^(k-1) + k * sfield_y.^(k-1) + ... - k * sfield_z.^(k-1); - - % remove boundary for divergence; int(erior) only - analytic_int = zeros(m+2, n+2, o+2); - analytic_int(2:end-1, 2:end-1, 2:end-1) = analytic(2:end-1, 2:end-1, 2:end-1); - analytic_int = reshape(analytic_int, [], 1); - - err = approx - analytic_int; - - testCase.verifyLessThan(norm(err), tol, ... - sprintf("Nullity test failed for k = %d", k)); - end - end - end -end \ No newline at end of file +function testDivPolynomial () +endfunction + + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave")); +%! +%! ks = [2, 4, 6, 8]; +%! tol = 1e-11; +%! x1 = 0; +%! x2 = 1; +%! +%! for k = ks +%! m = 2 * k + 1; +%! dx = (x2 - x1) / m; +%! +%! D = div (k, m, dx); +%! +%! sfield = [x1, x1 + dx/2:dx:x2 - dx/2, x2]'; +%! vfield = (x1:dx:x2)'; +%! +%! approx = D * vfield .^ k; +%! analytic = k * sfield .^ (k - 1); +%! +%! % Ignore boundary points. +%! err = approx (2:end-1) - analytic (2:end-1); +%! +%! assert (norm (err) < tol, ... +%! sprintf ("1D divergence test failed for k = %d", k)); +%! endfor + + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave")); +%! +%! ks = [2, 4, 6, 8]; +%! tol = 1e-11; +%! x1 = 0; +%! x2 = 1; +%! y1 = 0; +%! y2 = 1; +%! +%! for k = ks +%! m = 2 * k + 1; +%! dx = (x2 - x1) / m; +%! n = m + 2; +%! dy = (y2 - y1) / n; +%! +%! D = div2D (k, m, dx, n, dy); +%! +%! % Staggered grid components. +%! nodes_x = (x1:dx:x2)'; +%! nodes_y = (y1:dy:y2)'; +%! centers_x = (x1 + dx/2:dx:x2 - dx/2)'; +%! centers_y = (y1 + dy/2:dy:y2 - dy/2)'; +%! +%! % x-component of vector field. +%! [~, vfield_x] = meshgrid (centers_y, nodes_x); +%! vfield_x = reshape (vfield_x, [], 1); +%! +%! % y-component of vector field. +%! [vfield_y, ~] = meshgrid (nodes_y, centers_x); +%! vfield_y = reshape (vfield_y, [], 1); +%! +%! % Scalar field. +%! [sfield_y, sfield_x] = meshgrid ( ... +%! [y1, centers_y', y2]', ... +%! [x1, centers_x', x2]'); +%! +%! % Approximate and analytical divergence. +%! approx = D * [vfield_x; vfield_y] .^ k; +%! analytic = k * sfield_x .^ (k - 1) + ... +%! k * sfield_y .^ (k - 1); +%! +%! % Retain only interior values. +%! analytic_int = zeros (m + 2, n + 2); +%! analytic_int (2:end-1, 2:end-1) = ... +%! analytic (2:end-1, 2:end-1); +%! analytic_int = reshape (analytic_int, [], 1); +%! +%! err = approx - analytic_int; +%! +%! assert (norm (err) < tol, ... +%! sprintf ("2D divergence test failed for k = %d", k)); +%! endfor + + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave")); +%! +%! ks = [2, 4, 6, 8]; +%! tol = 1e-11; +%! x1 = 0; +%! x2 = 1; +%! y1 = 0; +%! y2 = 1; +%! z1 = 0; +%! z2 = 1; +%! +%! for k = ks +%! m = 2 * k + 1; +%! dx = (x2 - x1) / m; +%! n = m + 2; +%! dy = (y2 - y1) / n; +%! o = n + 2; +%! dz = (z2 - z1) / o; +%! +%! D = div3D (k, m, dx, n, dy, o, dz); +%! +%! % Staggered grid components. +%! nodes_x = (x1:dx:x2)'; +%! nodes_y = (y1:dy:y2)'; +%! nodes_z = (z1:dz:z2)'; +%! +%! centers_x = (x1 + dx/2:dx:x2 - dx/2)'; +%! centers_y = (y1 + dy/2:dy:y2 - dy/2)'; +%! centers_z = (z1 + dz/2:dz:z2 - dz/2)'; +%! +%! % x-component of vector field. +%! [~, vfield_x, ~] = ... +%! meshgrid (centers_y, nodes_x, centers_z); +%! vfield_x = reshape (vfield_x, [], 1); +%! +%! % y-component of vector field. +%! [vfield_y, ~, ~] = ... +%! meshgrid (nodes_y, centers_x, centers_z); +%! vfield_y = reshape (vfield_y, [], 1); +%! +%! % z-component of vector field. +%! [~, ~, vfield_z] = ... +%! meshgrid (centers_y, centers_x, nodes_z); +%! vfield_z = reshape (vfield_z, [], 1); +%! +%! % Scalar field. +%! [sfield_y, sfield_x, sfield_z] = meshgrid ( ... +%! [y1, centers_y', y2]', ... +%! [x1, centers_x', x2]', ... +%! [z1, centers_z', z2]'); +%! +%! % Approximate and analytical divergence. +%! approx = D * [vfield_x; vfield_y; vfield_z] .^ k; +%! analytic = k * sfield_x .^ (k - 1) + ... +%! k * sfield_y .^ (k - 1) + ... +%! k * sfield_z .^ (k - 1); +%! +%! % Retain only interior values. +%! analytic_int = zeros (m + 2, n + 2, o + 2); +%! analytic_int (2:end-1, 2:end-1, 2:end-1) = ... +%! analytic (2:end-1, 2:end-1, 2:end-1); +%! analytic_int = reshape (analytic_int, [], 1); +%! +%! err = approx - analytic_int; +%! +%! assert (norm (err) < tol, ... +%! sprintf ("3D divergence test failed for k = %d", k)); +%! endfor From 99055b156667dc44720a7a5b97ec25851b62444b Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 18 Oct 2026 21:30:58 -0700 Subject: [PATCH 07/17] Add an individual Octave testEnergy --- tests/octave/Makefile | 1 + tests/octave/testEnergy.m | 83 ++++++++++++++++++++------------------- 2 files changed, 43 insertions(+), 41 deletions(-) diff --git a/tests/octave/Makefile b/tests/octave/Makefile index 1a69618ba..a8430278a 100644 --- a/tests/octave/Makefile +++ b/tests/octave/Makefile @@ -5,5 +5,6 @@ run_octave_tests: octave --quiet --eval "addpath(pwd()); test('testAddScalarBC')" octave --quiet --eval "addpath(pwd()); test('testBCConsistency')" octave --quiet --eval "addpath(pwd()); test('testDivPolynomial')" + octave --quiet --eval "addpath(pwd()); test('testEnergy')" diff --git a/tests/octave/testEnergy.m b/tests/octave/testEnergy.m index 25b04d3f8..3b78b453a 100644 --- a/tests/octave/testEnergy.m +++ b/tests/octave/testEnergy.m @@ -1,43 +1,44 @@ -classdef testEnergy < matlab.unittest.TestCase - methods(Test) - function testforEnergy(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - %parameteres - k = 4; - a = -5; - b = 5; - m = 500; - grid = linspace(a, b, m); - dx = grid(2) - grid(1); - tol = 1e-6; - - % Laplacian - L = lap(k, m - 2, dx); +function testEnergy () +endfunction - % Transform grid using an anonymous function - grid = arrayfun(@(x) x^2, grid); - - % Potential matrix V - V = spdiags(grid', 0, m, m); % Transpose grid to match dimension - - % Hamiltonian H - H = -0.5 * L + V; - - % Eigenvalue computation - eigval = eig(full(H)); % Convert sparse matrix to full for eig function - eigval = sort(eigval); - - % Expected eigenvalues - expected = [1, 3, 5, 7, 9]; - - for i= 1:length(expected) - testCase.verifyLessThan(abs(real(eigval(i) / eigval(1)) - expected(i)),tol, ... - sprintf("Nullity test failed for k = %d", k)); - end - end - end -end +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave")); +%! +%! % Parameters +%! k = 4; +%! a = -5; +%! b = 5; +%! m = 500; +%! grid = linspace (a, b, m); +%! dx = grid (2) - grid (1); +%! tol = 1e-6; +%! +%! % Laplacian +%! L = lap (k, m - 2, dx); +%! +%! % Transform the grid +%! grid = arrayfun (@(x) x^2, grid); +%! +%! % Potential matrix +%! V = spdiags (grid', 0, m, m); +%! +%! % Hamiltonian +%! H = -0.5 * L + V; +%! +%! % Eigenvalue computation +%! eigval = eig (full (H)); +%! eigval = sort (eigval); +%! +%! % Expected eigenvalue ratios +%! expected = [1, 3, 5, 7, 9]; +%! +%! for i = 1:length (expected) +%! error_value = abs (real (eigval (i) / eigval (1)) ... +%! - expected (i)); +%! +%! assert (error_value < tol, ... +%! sprintf ("Energy test failed for k = %d", k)); +%! endfor From 59de6dd55d583a1df23242c34a7039356d77d455 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Sun, 18 Oct 2026 22:19:15 -0700 Subject: [PATCH 08/17] Add an individual Octave testGI13Curv --- tests/octave/Makefile | 1 + tests/octave/testGI13Curv.m | 343 ++++++++++++++++++++---------------- 2 files changed, 194 insertions(+), 150 deletions(-) diff --git a/tests/octave/Makefile b/tests/octave/Makefile index a8430278a..3baa366bf 100644 --- a/tests/octave/Makefile +++ b/tests/octave/Makefile @@ -6,5 +6,6 @@ run_octave_tests: octave --quiet --eval "addpath(pwd()); test('testBCConsistency')" octave --quiet --eval "addpath(pwd()); test('testDivPolynomial')" octave --quiet --eval "addpath(pwd()); test('testEnergy')" + octave --quiet --eval "addpath(pwd()); test('testGI13Curv')" diff --git a/tests/octave/testGI13Curv.m b/tests/octave/testGI13Curv.m index de3c59c32..5b71d63d6 100644 --- a/tests/octave/testGI13Curv.m +++ b/tests/octave/testGI13Curv.m @@ -1,159 +1,202 @@ -classdef testGI13Curv < matlab.unittest.TestCase -% Regression tests for GI13, the face-to-face interpolation used by -% grad3DCurvLegacy (which the four-argument grad3DCurv dispatches to). -% -% The defect these guard against was an index-map error: GI13 built -% kron(speye(n*o), I1), walking the source in n*o blocks when the source has -% (n+1)*o blocks, so the map drifted by one eta-plane per zeta-level. -% -% NOTE ON TEST CHOICE. The broken operator had the RIGHT SHAPE and mapped a -% constant to a constant, so shape checks and row-sum checks both passed on it. -% Only an explicit index map (testIndexMap) or a convergence study -% (testGradConvergesOnCurvilinearGrid) detects it. Please keep both. +function testGI13Curv () +endfunction - methods(Static) - function [X, Y, Z, xc, yc, zc] = sineGrid(N, amp) - % Nodes, and the cell centres with a ghost layer, for a smoothly - % distorted grid. Matches examples/matlab_octave conventions. - [a, b, c] = meshgrid(linspace(0,1,N), linspace(0,1,N), linspace(0,1,N)); - X = a + amp*sin(2*pi*b); - Y = b + amp*sin(2*pi*a); - Z = c + 0.5*amp*sin(2*pi*a).*sin(2*pi*b); - s = linspace(0,1,N); - s = [0, s(1:end-1)+0.5/(N-1), 1]; - [p, q, t] = meshgrid(s, s, s); - xc = p + amp*sin(2*pi*q); - yc = q + amp*sin(2*pi*p); - zc = t + 0.5*amp*sin(2*pi*p).*sin(2*pi*q); - end - end +%!function [X, Y, Z, xc, yc, zc] = sineGrid (N, amp) +%! [a, b, c] = meshgrid ( ... +%! linspace (0, 1, N), ... +%! linspace (0, 1, N), ... +%! linspace (0, 1, N)); +%! +%! X = a + amp * sin (2 * pi * b); +%! Y = b + amp * sin (2 * pi * a); +%! Z = c + 0.5 * amp * sin (2 * pi * a) .* sin (2 * pi * b); +%! +%! s = linspace (0, 1, N); +%! s = [0, s(1:end-1) + 0.5 / (N - 1), 1]; +%! +%! [p, q, t] = meshgrid (s, s, s); +%! +%! xc = p + amp * sin (2 * pi * q); +%! yc = q + amp * sin (2 * pi * p); +%! zc = t + 0.5 * amp * sin (2 * pi * p) .* sin (2 * pi * q); +%!endfunction - methods(Test) +%!test +%! % Check the explicit GI13 index map. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! m = 4; +%! n = 3; +%! o = 3; +%! +%! [I, J, K] = ndgrid (1:m, 1:n+1, 1:o); +%! v = I(:) + 100 * J(:) + 10000 * K(:); +%! +%! A = GI13 (speye (m * (n + 1) * o), m, n, o, "Gn"); +%! W = reshape (A * v, m + 1, n, o); +%! +%! bad = 0; +%! +%! for k = 1:o +%! for j = 1:n +%! for i = 1:m+1 +%! val = W (i, j, k); +%! +%! if (val == 0) +%! continue; +%! endif +%! +%! kk = floor (val / 10000); +%! +%! if (kk != k) +%! bad = bad + 1; +%! endif +%! endfor +%! endfor +%! endfor +%! +%! assert (bad == 0, ... +%! sprintf ("%d of %d GI13 outputs read from the wrong zeta-plane", ... +%! bad, (m + 1) * n * o)); - function testIndexMap(testCase) - % Feed GI13 a field whose values encode their own (i,j,k) and - % check every output reads from the correct zeta-plane. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - m = 4; - n = 3; - o = 3; - [I, J, K] = ndgrid(1:m, 1:n+1, 1:o); - v = I(:) + 100*J(:) + 10000*K(:); - W = reshape(GI13(speye(m*(n+1)*o), m, n, o, 'Gn')*v, m+1, n, o); +%!test +%! % Check shapes and preservation of constant fields. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! m = 5; +%! n = 4; +%! o = 3; +%! +%! %ty = {"Gn", "Gc", "Ge", "Gcy", "Gee", "Gnn"}; +%! ty = {'Gn', 'Gc', 'Ge', 'Gcy', 'Gee', 'Gnn'}; +%! +%! src = [ ... +%! m * (n + 1) * o, ... +%! m * n * (o + 1), ... +%! (m + 1) * n * o, ... +%! m * n * (o + 1), ... +%! (m + 1) * n * o, ... +%! m * (n + 1) * o]; +%! +%! tgt = [ ... +%! (m + 1) * n * o, ... +%! (m + 1) * n * o, ... +%! m * (n + 1) * o, ... +%! m * (n + 1) * o, ... +%! m * n * (o + 1), ... +%! m * n * (o + 1)]; +%! +%! for i = 1:numel (ty) +%! clear GI13; +%! A = GI13 (speye (src (i)), m, n, o, ty{i}); +%! +%! expectedSize = [tgt(i), src(i)]; +%! actualSize = size (A); +%! +% fprintf ("GI13 type: %s\n", ty{i}); +% fprintf ("GI13 file: %s\n", which ("GI13")); +% fprintf ("Actual size: [%d %d]\n", actualSize); +% fprintf ("Expected size: [%d %d]\n", expectedSize); +%! +%! assert (all (actualSize == expectedSize), ... +%! sprintf ("GI13 %s has the wrong shape", ty{i})); +%! +%! assert (max (abs (A * ones (src(i), 1) - 1)) < 1e-12, ... +%! sprintf ("GI13 %s does not preserve a constant", ty{i})); +%! endfor - bad = 0; - for k = 1:o - for j = 1:n - for i = 1:m+1 - val = W(i,j,k); - if val == 0 - continue; - end - kk = floor(val/10000); - if kk ~= k - bad = bad + 1; - end - end - end - end - testCase.verifyEqual(bad, 0, ... - sprintf('%d of %d GI13 outputs read from the wrong zeta-plane', ... - bad, (m+1)*n*o)); - end +%!test +%! % Check exactness on a Cartesian grid. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! for N = [11, 17, 25] +%! [X, Y, Z, xc, yc, zc] = sineGrid (N, 0.0); +%! +%! f = reshape ( ... +%! permute (xc.^2 + yc.^2 + zc.^2, [2, 1, 3]), ... +%! [], 1); +%! +%! L = div3DCurv (2, X, Y, Z) * grad3DCurv (2, X, Y, Z); +%! +%! v = reshape (L * f, N + 1, N + 1, N + 1); +%! d = v (3:end-2, 3:end-2, 3:end-2); +%! d = d(:) - 6; +%! +%! assert (sqrt (mean (d.^2)) < 1e-9, ... +%! sprintf ("Cartesian L = D*G is not exact at N = %d", N)); +%! endfor - function testShapesAndConstants(testCase) - % All six shifts must have the right shape and map a constant - % field to a constant field. Weak, but cheap. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - m = 5; - n = 4; - o = 3; - ty = {'Gn','Gc','Ge','Gcy','Gee','Gnn'}; - src = [m*(n+1)*o, m*n*(o+1), (m+1)*n*o, m*n*(o+1), (m+1)*n*o, m*(n+1)*o]; - tgt = [(m+1)*n*o, (m+1)*n*o, m*(n+1)*o, m*(n+1)*o, m*n*(o+1), m*n*(o+1)]; +%!test +%! % Check the gradient of a linear field on a curvilinear grid. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! N = 21; +%! +%! for amp = [0.05, 0.10] +%! [X, Y, Z, xc, yc, zc] = sineGrid (N, amp); +%! +%! f = reshape ( ... +%! permute (xc + 2 * yc + 3 * zc, [2, 1, 3]), ... +%! [], 1); +%! +%! T = grad3DCurv (2, X, Y, Z) * f; +%! +%! gx = T (1:N * (N - 1) * (N - 1)) - 1; +%! +%! assert (sqrt (mean (gx.^2)) < 0.05, ... +%! sprintf ( ... +%! "Gradient of a linear field is badly wrong at amp = %g", ... +%! amp)); +%! endfor - for i = 1:numel(ty) - A = GI13(speye(src(i)), m, n, o, ty{i}); - testCase.verifyEqual(size(A), [tgt(i) src(i)], ... - sprintf('GI13 %s has the wrong shape', ty{i})); - testCase.verifyLessThan(max(abs(A*ones(src(i),1) - 1)), 1e-12, ... - sprintf('GI13 %s does not preserve a constant', ty{i})); - end - end - function testGradExactOnCartesianGrid(testCase) - % On an undistorted grid the composition must be exact to roundoff. - % Guards against a "fix" that trades curvilinear accuracy for - % Cartesian accuracy. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - for N = [11 17 25] - [X, Y, Z, xc, yc, zc] = testGI13Curv.sineGrid(N, 0.0); - f = reshape(permute(xc.^2 + yc.^2 + zc.^2, [2 1 3]), [], 1); - L = div3DCurv(2, X, Y, Z)*grad3DCurv(2, X, Y, Z); - v = reshape(L*f, N+1, N+1, N+1); - d = v(3:end-2, 3:end-2, 3:end-2); - d = d(:) - 6; - testCase.verifyLessThan(sqrt(mean(d.^2)), 1e-9, ... - sprintf('Cartesian L=D*G not exact at n=%d', N)); - end - end - - function testGradExactOnLinearField(testCase) - % The gradient of f = x + 2y + 3z is a constant. The stock - % operator gave rms 0.39 at distortion 0.10; anything near that - % means the face shift is wrong again. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - N = 21; - for amp = [0.05 0.10] - [X, Y, Z, xc, yc, zc] = testGI13Curv.sineGrid(N, amp); - f = reshape(permute(xc + 2*yc + 3*zc, [2 1 3]), [], 1); - T = grad3DCurv(2, X, Y, Z)*f; - gx = T(1:N*(N-1)*(N-1)) - 1; - testCase.verifyLessThan(sqrt(mean(gx.^2)), 0.05, ... - sprintf('grad of a linear field is badly wrong at amp=%g', amp)); - end - end - - function testGradConvergesOnCurvilinearGrid(testCase) - % The headline check. The stock operator gave order 0.03-0.09 -- - % no convergence at all. Require better than first order, which - % the stock operator misses by a wide margin and the corrected one - % clears at ~2.1. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - amp = 0.10; - Ns = [13 21 33]; - errs = zeros(size(Ns)); - for i = 1:numel(Ns) - N = Ns(i); - [X, Y, Z, xc, yc, zc] = testGI13Curv.sineGrid(N, amp); - f = reshape(permute(xc.^2 + yc.^2 + zc.^2, [2 1 3]), [], 1); - T = grad3DCurv(2, X, Y, Z)*f; - cx = linspace(0,1,N); - cx = cx(1:end-1) + 0.5/(N-1); - [aa, bb, ~] = meshgrid(linspace(0,1,N), cx, cx); - xf = aa + amp*sin(2*pi*bb); - gx = T(1:N*(N-1)*(N-1)); - errs(i) = sqrt(mean((gx - reshape(permute(2*xf,[2 1 3]),[],1)).^2)); - end - p = log(errs(1)/errs(end)) / log((Ns(end)-1)/(Ns(1)-1)); - testCase.verifyGreaterThan(p, 1.5, ... - sprintf('grad3DCurv observed order %.2f on a curvilinear grid', p)); - end - - end -end +%!test +%! % Check convergence on a curvilinear grid. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! amp = 0.10; +%! Ns = [13, 21, 33]; +%! errs = zeros (size (Ns)); +%! +%! for i = 1:numel (Ns) +%! N = Ns (i); +%! [X, Y, Z, xc, yc, zc] = sineGrid (N, amp); +%! +%! f = reshape ( ... +%! permute (xc.^2 + yc.^2 + zc.^2, [2, 1, 3]), ... +%! [], 1); +%! +%! T = grad3DCurv (2, X, Y, Z) * f; +%! +%! cx = linspace (0, 1, N); +%! cx = cx (1:end-1) + 0.5 / (N - 1); +%! +%! [aa, bb, ~] = meshgrid (linspace (0, 1, N), cx, cx); +%! xf = aa + amp * sin (2 * pi * bb); +%! +%! gx = T (1:N * (N - 1) * (N - 1)); +%! exact_gx = reshape ( ... +%! permute (2 * xf, [2, 1, 3]), ... +%! [], 1); +%! +%! errs (i) = sqrt (mean ((gx - exact_gx).^2)); +%! endfor +%! +%! p = log (errs (1) / errs (end)) / ... +%! log ((Ns (end) - 1) / (Ns (1) - 1)); +%! +%! assert (p > 1.5, ... +%! sprintf ( ... +%! "Observed gradient convergence order %.2f is too low", p)); From fb60ad3808d2a150381b3f49306e85d901befbe3 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Mon, 19 Oct 2026 11:30:35 -0700 Subject: [PATCH 09/17] Add an individual Octave testGradPolynomial --- .github/workflows/matlab-tests.yml | 2 +- tests/octave/Makefile | 21 +- tests/octave/testGradPolynomial.m | 301 +++++++++++++++-------------- 3 files changed, 171 insertions(+), 153 deletions(-) diff --git a/.github/workflows/matlab-tests.yml b/.github/workflows/matlab-tests.yml index 82578468f..2d88f6d60 100644 --- a/.github/workflows/matlab-tests.yml +++ b/.github/workflows/matlab-tests.yml @@ -34,7 +34,7 @@ jobs: test-results-junit: test-results/results.xml code-coverage-cobertura: code-coverage/coverage.xml source-folder: src/octave - select-by-folder: tests/octave + select-by-folder: tests/matlab strict: false use-parallel: false generate-summary: true diff --git a/tests/octave/Makefile b/tests/octave/Makefile index 3baa366bf..985e30ee8 100644 --- a/tests/octave/Makefile +++ b/tests/octave/Makefile @@ -1,11 +1,14 @@ run_octave_tests: - octave --quiet --eval "addpath(pwd()); test('testGradient')" - octave --quiet --eval "addpath(pwd()); test('testDivergence')" - octave --quiet --eval "addpath(pwd()); test('testLaplacian')" - octave --quiet --eval "addpath(pwd()); test('testAddScalarBC')" - octave --quiet --eval "addpath(pwd()); test('testBCConsistency')" - octave --quiet --eval "addpath(pwd()); test('testDivPolynomial')" - octave --quiet --eval "addpath(pwd()); test('testEnergy')" - octave --quiet --eval "addpath(pwd()); test('testGI13Curv')" - + for test_name in \ + testGradient \ + testDivergence \ + testLaplacian \ + testAddScalarBC \ + testBCConsistency \ + testDivPolynomial \ + testEnergy \ + testGI13Curv \ + testGradPolynomial; do \ + octave --quiet --eval "addpath(pwd()); test('$$test_name')" || exit 1; \ + done diff --git a/tests/octave/testGradPolynomial.m b/tests/octave/testGradPolynomial.m index a7068058b..0c7918d98 100644 --- a/tests/octave/testGradPolynomial.m +++ b/tests/octave/testGradPolynomial.m @@ -1,143 +1,158 @@ -classdef testGradPolynomial < matlab.unittest.TestCase - methods (Test) - function testGradOfPolynomial1D(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - ks = [2, 4, 6, 8]; % Different orders of accuracy - tol = 1e-11; - x1 = 0; - x2 = 1; - - for k = ks - m = 2 * k + 1; - dx = (x2-x1) / m; - - G = grad(k, m, dx); - - sfield = [x1, x1+dx/2:dx:x2-dx/2, x2]'; - vfield = (x1:dx:x2)'; - - approx = G * sfield.^k; - analytic = k * vfield.^(k-1); - - err = approx - analytic; - - testCase.verifyLessThan(norm(err), tol, ... - sprintf("Nullity test failed for k = %d", k)); - end - end - - function testGradOfPolynomial2D(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - ks = [2, 4, 6, 8]; % Different orders of accuracy - tol = 1e-11; - x1 = 0; - x2 = 1; - y1 = 0; - y2 = 1; - - for k = ks - m = 2 * k + 1; - dx = (x2 - x1) / m; - n = m + 2; - dy = (y2 - y1) / n; - - G = grad2D(k, m, dx, n, dy); - - % staggered grid components - nodes_x = (x1:dx:x2)'; - nodes_y = (y1:dy:y2)'; - centers_x = (x1+dx/2:dx:x2-dx/2)'; - centers_y = (y1+dy/2:dy:y2-dy/2)'; - - % set up vector field - [~, vfield_x] = meshgrid(centers_y, nodes_x); - vfield_x = reshape(vfield_x, [], 1); - [vfield_y, ~] = meshgrid(nodes_y, centers_x); - vfield_y = reshape(vfield_y, [], 1); - - % set up scalar field - [sfield_y, sfield_x] = meshgrid([y1,centers_y',y2]', [x1,centers_x',x2]'); - sfield_x = reshape(sfield_x, [], 1); - sfield_y = reshape(sfield_y, [], 1); - - % compute approximate and analytical solutions - approx = G * (sfield_x.^k + sfield_y.^k); - analytic_x = k * vfield_x.^(k-1); - analytic_y = k * vfield_y.^(k-1); - analytic_combined = [analytic_x; analytic_y]; - - err = approx - analytic_combined; - - testCase.verifyLessThan(norm(err), tol, ... - sprintf("Nullity test failed for k = %d", k)); - end - end - - function testGradOfPolynomial3D(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - ks = [2, 4, 6, 8]; % Different orders of accuracy - tol = 1e-11; - x1 = 0; - x2 = 1; - y1 = 0; - y2 = 1; - z1 = 0; - z2 = 1; - - for k = ks - m = 2 * k + 1; - dx = (x2 - x1) / m; - n = m + 2; - dy = (y2 - y1) / n; - o = n + 2; - dz = (z2 - z1) / o; - - G = grad3D(k, m, dx, n, dy, o, dz); - - % staggered grid components - nodes_x = (x1:dx:x2)'; - nodes_y = (y1:dy:y2)'; - nodes_z = (z1:dz:z2)'; - centers_x = (x1+dx/2:dx:x2-dx/2)'; - centers_y = (y1+dy/2:dy:y2-dy/2)'; - centers_z = (z1+dz/2:dz:z2-dz/2)'; - - % set up vector field - [~, vfield_x, ~] = meshgrid(centers_y, nodes_x, centers_z); - vfield_x = reshape(vfield_x, [], 1); - [vfield_y, ~, ~] = meshgrid(nodes_y, centers_x, centers_z); - vfield_y = reshape(vfield_y, [], 1); - [~, ~, vfield_z] = meshgrid(centers_y, centers_x, nodes_z); - vfield_z = reshape(vfield_z, [], 1); - - % set up scalar field - [sfield_y, sfield_x, sfield_z] = meshgrid([y1,centers_y',y2]', ... - [x1,centers_x',x2]', [z1,centers_z',z2]'); - sfield_x = reshape(sfield_x, [], 1); - sfield_y = reshape(sfield_y, [], 1); - sfield_z = reshape(sfield_z, [], 1); - - % compute approximate and analytical solutions - approx = G * (sfield_x.^k + sfield_y.^k + sfield_z.^k); - analytic_x = k * vfield_x.^(k-1); - analytic_y = k * vfield_y.^(k-1); - analytic_z = k * vfield_z.^(k-1); - analytic_combined = [analytic_x; analytic_y; analytic_z]; - - err = approx - analytic_combined; - - testCase.verifyLessThan(norm(err), tol, ... - sprintf("Nullity test failed for k = %d", k)); - end - end - end -end \ No newline at end of file +function testGradPolynomial () +endfunction + + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! ks = [2, 4, 6, 8]; +%! tol = 1e-11; +%! x1 = 0; +%! x2 = 1; +%! +%! for k = ks +%! m = 2 * k + 1; +%! dx = (x2 - x1) / m; +%! +%! G = grad (k, m, dx); +%! +%! sfield = [x1, x1 + dx/2:dx:x2 - dx/2, x2]'; +%! vfield = (x1:dx:x2)'; +%! +%! approx = G * sfield .^ k; +%! analytic = k * vfield .^ (k - 1); +%! +%! err = approx - analytic; +%! +%! assert (norm (err) < tol, ... +%! sprintf ("1D gradient test failed for k = %d", k)); +%! endfor + + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! ks = [2, 4, 6, 8]; +%! tol = 1e-11; +%! x1 = 0; +%! x2 = 1; +%! y1 = 0; +%! y2 = 1; +%! +%! for k = ks +%! m = 2 * k + 1; +%! dx = (x2 - x1) / m; +%! n = m + 2; +%! dy = (y2 - y1) / n; +%! +%! G = grad2D (k, m, dx, n, dy); +%! +%! % Staggered grid components. +%! nodes_x = (x1:dx:x2)'; +%! nodes_y = (y1:dy:y2)'; +%! centers_x = (x1 + dx/2:dx:x2 - dx/2)'; +%! centers_y = (y1 + dy/2:dy:y2 - dy/2)'; +%! +%! % Vector field locations. +%! [~, vfield_x] = meshgrid (centers_y, nodes_x); +%! vfield_x = reshape (vfield_x, [], 1); +%! +%! [vfield_y, ~] = meshgrid (nodes_y, centers_x); +%! vfield_y = reshape (vfield_y, [], 1); +%! +%! % Scalar field locations. +%! [sfield_y, sfield_x] = meshgrid ( ... +%! [y1, centers_y', y2]', ... +%! [x1, centers_x', x2]'); +%! +%! sfield_x = reshape (sfield_x, [], 1); +%! sfield_y = reshape (sfield_y, [], 1); +%! +%! approx = G * (sfield_x .^ k + sfield_y .^ k); +%! +%! analytic_x = k * vfield_x .^ (k - 1); +%! analytic_y = k * vfield_y .^ (k - 1); +%! analytic_combined = [analytic_x; analytic_y]; +%! +%! err = approx - analytic_combined; +%! +%! assert (norm (err) < tol, ... +%! sprintf ("2D gradient test failed for k = %d", k)); +%! endfor + + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! ks = [2, 4, 6, 8]; +%! tol = 1e-11; +%! x1 = 0; +%! x2 = 1; +%! y1 = 0; +%! y2 = 1; +%! z1 = 0; +%! z2 = 1; +%! +%! for k = ks +%! m = 2 * k + 1; +%! dx = (x2 - x1) / m; +%! n = m + 2; +%! dy = (y2 - y1) / n; +%! o = n + 2; +%! dz = (z2 - z1) / o; +%! +%! G = grad3D (k, m, dx, n, dy, o, dz); +%! +%! % Staggered grid components. +%! nodes_x = (x1:dx:x2)'; +%! nodes_y = (y1:dy:y2)'; +%! nodes_z = (z1:dz:z2)'; +%! +%! centers_x = (x1 + dx/2:dx:x2 - dx/2)'; +%! centers_y = (y1 + dy/2:dy:y2 - dy/2)'; +%! centers_z = (z1 + dz/2:dz:z2 - dz/2)'; +%! +%! % Vector field locations. +%! [~, vfield_x, ~] = ... +%! meshgrid (centers_y, nodes_x, centers_z); +%! vfield_x = reshape (vfield_x, [], 1); +%! +%! [vfield_y, ~, ~] = ... +%! meshgrid (nodes_y, centers_x, centers_z); +%! vfield_y = reshape (vfield_y, [], 1); +%! +%! [~, ~, vfield_z] = ... +%! meshgrid (centers_y, centers_x, nodes_z); +%! vfield_z = reshape (vfield_z, [], 1); +%! +%! % Scalar field locations. +%! [sfield_y, sfield_x, sfield_z] = meshgrid ( ... +%! [y1, centers_y', y2]', ... +%! [x1, centers_x', x2]', ... +%! [z1, centers_z', z2]'); +%! +%! sfield_x = reshape (sfield_x, [], 1); +%! sfield_y = reshape (sfield_y, [], 1); +%! sfield_z = reshape (sfield_z, [], 1); +%! +%! approx = G * ( ... +%! sfield_x .^ k + ... +%! sfield_y .^ k + ... +%! sfield_z .^ k); +%! +%! analytic_x = k * vfield_x .^ (k - 1); +%! analytic_y = k * vfield_y .^ (k - 1); +%! analytic_z = k * vfield_z .^ (k - 1); +%! analytic_combined = [analytic_x; analytic_y; analytic_z]; +%! +%! err = approx - analytic_combined; +%! +%! assert (norm (err) < tol, ... +%! sprintf ("3D gradient test failed for k = %d", k)); +%! endfor From f72f326fe267c925f32c599760fd81c0ac41c782 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Mon, 19 Oct 2026 11:45:30 -0700 Subject: [PATCH 10/17] Add an individual Octave testGridOrientation3D --- tests/octave/Makefile | 3 +- tests/octave/testGridOrientation3D.m | 215 ++++++++++++++++----------- 2 files changed, 131 insertions(+), 87 deletions(-) diff --git a/tests/octave/Makefile b/tests/octave/Makefile index 985e30ee8..19059b2da 100644 --- a/tests/octave/Makefile +++ b/tests/octave/Makefile @@ -8,7 +8,8 @@ run_octave_tests: testDivPolynomial \ testEnergy \ testGI13Curv \ - testGradPolynomial; do \ + testGradPolynomial \ + testGridOrientation3D; do \ octave --quiet --eval "addpath(pwd()); test('$$test_name')" || exit 1; \ done diff --git a/tests/octave/testGridOrientation3D.m b/tests/octave/testGridOrientation3D.m index 5cd4c1de9..ff16c2a95 100644 --- a/tests/octave/testGridOrientation3D.m +++ b/tests/octave/testGridOrientation3D.m @@ -1,95 +1,138 @@ -classdef testGridOrientation3D < matlab.unittest.TestCase -% Extends the 2D orientation guard to jacobian3D, and adds detection of -% vanishing Jacobians in both dimensions. -% -% In 3D, swapping two axes is an odd permutation, so handedness inverts and J -% goes uniformly negative -- the same failure as the 2D transpose, with one -% more way to trigger it (mirroring a single axis does it too). -% -% A vanishing Jacobian is a different problem: the mapping is degenerate at -% that cell and is not invertible there, and grad/div divide by J. It is not -% caught by the sign checks, because zero is neither positive nor negative. +function testGridOrientation3D () +endfunction - methods(Static) - function [X, Y, Z] = sineGrid3(N, amp) - [a, b, c] = meshgrid(linspace(0,1,N), linspace(0,1,N), ... - linspace(0,1,N)); - X = a + amp*sin(2*pi*b); - Y = b + amp*sin(2*pi*a); - Z = c + 0.5*amp*sin(2*pi*a).*sin(2*pi*b); - end - end - methods(Test) +%!function [X, Y, Z] = sineGrid3 (N, amp) +%! [a, b, c] = meshgrid ( ... +%! linspace (0, 1, N), ... +%! linspace (0, 1, N), ... +%! linspace (0, 1, N)); +%! +%! X = a + amp * sin (2 * pi * b); +%! Y = b + amp * sin (2 * pi * a); +%! Z = c + 0.5 * amp * sin (2 * pi * a) .* sin (2 * pi * b); +%!endfunction - function testRightHanded3DIsSilent(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - [X, Y, Z] = testGridOrientation3D.sineGrid3(13, 0.08); - testCase.verifyWarningFree(@() jacobian3D(2, X, Y, Z)); - end - function testSwappedAxesWarn(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - [X, Y, Z] = testGridOrientation3D.sineGrid3(13, 0.08); - p = @(A) permute(A, [2 1 3]); - testCase.verifyWarning(@() jacobian3D(2, p(X), p(Y), p(Z)), ... - 'jacobian3D:leftHandedGrid'); - end +%!test +%! % A properly oriented grid should not produce a warning. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! [X, Y, Z] = sineGrid3 (13, 0.08); +%! +%! lastwarn (""); +%! jacobian3D (2, X, Y, Z); +%! [warningMessage, warningId] = lastwarn (); +%! +%! assert (isempty (warningMessage) && isempty (warningId), ... +%! "A right-handed grid generated an unexpected warning"); - function testMirroredAxisWarns(testCase) - % Reversing one axis is also an odd permutation. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - [X, Y, Z] = testGridOrientation3D.sineGrid3(13, 0.08); - f = @(A) A(end:-1:1, :, :); - testCase.verifyWarning(@() jacobian3D(2, f(X), f(Y), f(Z)), ... - 'jacobian3D:leftHandedGrid'); - end - function testWarningReaches3DOperatorUsers(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - [X, Y, Z] = testGridOrientation3D.sineGrid3(13, 0.08); - p = @(A) permute(A, [2 1 3]); - testCase.verifyWarning(@() grad3DCurv(2, p(X), p(Y), p(Z)), ... - 'jacobian3D:leftHandedGrid'); - testCase.verifyWarning(@() div3DCurv(2, p(X), p(Y), p(Z)), ... - 'jacobian3D:leftHandedGrid'); - end +%!test +%! % Swapping two axes should produce a left-handed-grid warning. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! [X, Y, Z] = sineGrid3 (13, 0.08); +%! Xp = permute (X, [2, 1, 3]); +%! Yp = permute (Y, [2, 1, 3]); +%! Zp = permute (Z, [2, 1, 3]); +%! +%! lastwarn (""); +%! jacobian3D (2, Xp, Yp, Zp); +%! [warningMessage, warningId] = lastwarn (); +%! +%! assert (strcmp (warningId, "jacobian3D:leftHandedGrid"), ... +%! sprintf ("Unexpected warning identifier: %s", warningId)); - function testVanishingJacobianWarns(testCase) - % Three identical columns make the central difference across the - % middle one exactly zero, so those cells are degenerate. The sign - % checks cannot see this: zero is neither positive nor negative. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - N = 21; - [X, Y] = meshgrid(linspace(0,1,N), linspace(0,1,N)); - X = X + 0.10*sin(2*pi*Y); - Y = Y + 0.10*sin(2*pi*X); - X(:, 9:11) = repmat(X(:, 10), 1, 3); - Y(:, 9:11) = repmat(Y(:, 10), 1, 3); - testCase.verifyWarning(@() jacobian2D(2, X, Y), ... - 'jacobian2D:vanishingJacobian'); - end - function testHealthyGridHasNoVanishingWarning(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - N = 21; - [X, Y] = meshgrid(linspace(0,1,N), linspace(0,1,N)); - X = X + 0.10*sin(2*pi*Y); - Y = Y + 0.10*sin(2*pi*X); - testCase.verifyWarningFree(@() jacobian2D(2, X, Y)); - end +%!test +%! % Reversing one axis should also produce a left-handed-grid warning. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! [X, Y, Z] = sineGrid3 (13, 0.08); +%! +%! Xf = X (end:-1:1, :, :); +%! Yf = Y (end:-1:1, :, :); +%! Zf = Z (end:-1:1, :, :); +%! +%! lastwarn (""); +%! jacobian3D (2, Xf, Yf, Zf); +%! [warningMessage, warningId] = lastwarn (); +%! +%! assert (strcmp (warningId, "jacobian3D:leftHandedGrid"), ... +%! sprintf ("Unexpected warning identifier: %s", warningId)); - end -end + +%!test +%! % The warning should propagate through grad3DCurv and div3DCurv. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! [X, Y, Z] = sineGrid3 (13, 0.08); +%! Xp = permute (X, [2, 1, 3]); +%! Yp = permute (Y, [2, 1, 3]); +%! Zp = permute (Z, [2, 1, 3]); +%! +%! lastwarn (""); +%! grad3DCurv (2, Xp, Yp, Zp); +%! [warningMessage, warningId] = lastwarn (); +%! +%! assert (strcmp (warningId, "jacobian3D:leftHandedGrid"), ... +%! sprintf ("Unexpected gradient warning identifier: %s", warningId)); +%! +%! lastwarn (""); +%! div3DCurv (2, Xp, Yp, Zp); +%! [warningMessage, warningId] = lastwarn (); +%! +%! assert (strcmp (warningId, "jacobian3D:leftHandedGrid"), ... +%! sprintf ("Unexpected divergence warning identifier: %s", warningId)); + + +%!test +%! % Three identical columns create vanishing Jacobians. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! N = 21; +%! [X, Y] = meshgrid (linspace (0, 1, N), linspace (0, 1, N)); +%! +%! X = X + 0.10 * sin (2 * pi * Y); +%! Y = Y + 0.10 * sin (2 * pi * X); +%! +%! X (:, 9:11) = repmat (X (:, 10), 1, 3); +%! Y (:, 9:11) = repmat (Y (:, 10), 1, 3); +%! +%! lastwarn (""); +%! jacobian2D (2, X, Y); +%! [warningMessage, warningId] = lastwarn (); +%! +%! assert (strcmp (warningId, "jacobian2D:vanishingJacobian"), ... +%! sprintf ("Unexpected warning identifier: %s", warningId)); + + +%!test +%! % A healthy grid should not produce a vanishing-Jacobian warning. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! N = 21; +%! [X, Y] = meshgrid (linspace (0, 1, N), linspace (0, 1, N)); +%! +%! X = X + 0.10 * sin (2 * pi * Y); +%! Y = Y + 0.10 * sin (2 * pi * X); +%! +%! lastwarn (""); +%! jacobian2D (2, X, Y); +%! [warningMessage, warningId] = lastwarn (); +%! +%! assert (isempty (warningMessage) && isempty (warningId), ... +%! "A healthy grid generated an unexpected warning"); From 959aca57a335f4b75bf53033d320956fb75169ef Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Mon, 19 Oct 2026 11:50:02 -0700 Subject: [PATCH 11/17] Add an individual Octave testGridOrientation --- tests/octave/Makefile | 3 +- tests/octave/testGridOrientation.m | 270 +++++++++++++++-------------- 2 files changed, 145 insertions(+), 128 deletions(-) diff --git a/tests/octave/Makefile b/tests/octave/Makefile index 19059b2da..d9c913836 100644 --- a/tests/octave/Makefile +++ b/tests/octave/Makefile @@ -9,7 +9,8 @@ run_octave_tests: testEnergy \ testGI13Curv \ testGradPolynomial \ - testGridOrientation3D; do \ + testGridOrientation3D \ + testGridOrientation; do \ octave --quiet --eval "addpath(pwd()); test('$$test_name')" || exit 1; \ done diff --git a/tests/octave/testGridOrientation.m b/tests/octave/testGridOrientation.m index 2ccf787a7..c66e88226 100644 --- a/tests/octave/testGridOrientation.m +++ b/tests/octave/testGridOrientation.m @@ -1,127 +1,143 @@ -classdef testGridOrientation < matlab.unittest.TestCase -% Tests for the left-handed / tangled grid detection added to jacobian2D. -% -% The curvilinear operators read [n, m] = size(X): rows are eta, columns are -% xi. A grid supplied the other way round is still self-consistent for grad -% and div -- J flips sign, the cofactors flip with it, gradients still come -% out exact -- so nothing fails and nothing used to complain. These tests -% check that it complains now, and, just as importantly, that it stays quiet -% on correctly oriented grids. - - methods(Static) - function [X, Y] = sineGrid(N, amp) - [X, Y] = meshgrid(linspace(0,1,N), linspace(0,1,N)); - X = X + amp*sin(2*pi*Y); - Y = Y + amp*sin(2*pi*X); - end - end - - methods(Test) - - function testRightHandedIsSilent(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - [X, Y] = testGridOrientation.sineGrid(21, 0.10); - testCase.verifyWarningFree(@() jacobian2D(2, X, Y)); - end - - function testTransposedWarns(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - [X, Y] = testGridOrientation.sineGrid(21, 0.10); - testCase.verifyWarning(@() jacobian2D(2, X.', Y.'), ... - 'jacobian2D:leftHandedGrid'); - end - - function testGridGenOutputWarns(testCase) - % The reported case: gridGen returns xi-nodes as ROWS, which is the - % opposite of what the operators read, so passing its output - % straight in yields a globally left-handed grid. - % - % tfi does addpath(['grids/' grid_name]) RELATIVE TO THE CURRENT - % DIRECTORY, so the working directory has to be src/octave, not - % src/octave/grids. And the path is added non-recursively on - % purpose: genpath would put every grid folder on the path at once, - % and the curve functions (right.m, top.m, ...) share names across - % grids, so another grid would silently answer for this one. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - origDir = pwd; - cleanupDir = onCleanup(@() cd(origDir)); - srcDir = fullfile(pwd, '..', '..', 'src', 'octave'); - addpath(srcDir) - cd(srcDir) - - % chevron, not swan: swan is itself tangled (mixed-sign Jacobian - % either way round), which is the other warning. See - % testTangledGridWarns. - [X, Y] = gridGen('TFI', 'chevron', 39, 99, false); - testCase.verifyNotEmpty(strfind(which('right'), 'chevron'), ... - 'the chevron curves are not the ones on the path'); - testCase.verifyWarning(@() jacobian2D(2, X, Y), ... - 'jacobian2D:leftHandedGrid'); - testCase.verifyWarningFree(@() jacobian2D(2, X.', Y.')); - end - - function testMagnitudesAreIdentical(testCase) - % Documents the symptom that makes this hard to notice: same grid, - % same |J| to the last bit, opposite sign. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - [X, Y] = testGridOrientation.sineGrid(21, 0.10); - warning('off', 'jacobian2D:leftHandedGrid'); - restoreWarn = onCleanup(@() warning('on', 'jacobian2D:leftHandedGrid')); - a = sort(abs(nonzeros(jacobian2D(2, X, Y)))); - b = sort(abs(nonzeros(jacobian2D(2, X.', Y.')))); - testCase.verifyEqual(a, b, 'AbsTol', 1e-12); - end - - function testTangledGridWarns(testCase) - % A mesh that folds over itself has mixed-sign J. Different - % problem, different warning. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - [X, Y] = meshgrid(linspace(0,1,15), linspace(0,1,15)); - X(6:10, 6:10) = X(6:10, 6:10) - 0.45; - testCase.verifyWarning(@() jacobian2D(2, X, Y), ... - 'jacobian2D:tangledGrid'); - end - - function testWarningReachesOperatorUsers(testCase) - % grad2DCurv and div2DCurv both route through jacobian2D, so users - % who never call it directly still get told. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - [X, Y] = testGridOrientation.sineGrid(21, 0.10); - testCase.verifyWarning(@() grad2DCurv(2, X.', Y.'), ... - 'jacobian2D:leftHandedGrid'); - testCase.verifyWarning(@() div2DCurv(2, X.', Y.'), ... - 'jacobian2D:leftHandedGrid'); - end - - function testFiveArgPath(testCase) - % jacobian2D has two entry paths: 3 arguments dispatches to - % jacobian2DLegacy, 5 arguments computes inline. Both are checked. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - N = 21; - dc = [1;1;1;1]; - nc = [0;0;0;0]; - s = linspace(0,1,N); - s = [0, s(1:end-1)+0.5/(N-1), 1]; - [X, Y] = meshgrid(s, s); - X = X + 0.08*sin(2*pi*Y); - Y = Y + 0.08*sin(2*pi*X); - testCase.verifyWarningFree(@() jacobian2D(2, X, Y, dc, nc)); - testCase.verifyWarning(@() jacobian2D(2, X.', Y.', dc, nc), ... - 'jacobian2D:leftHandedGrid'); - end - - end -end +function testGridOrientation () +endfunction + + +%!function [X, Y] = sineGrid (N, amp) +%! [X, Y] = meshgrid (linspace (0, 1, N), linspace (0, 1, N)); +%! X = X + amp * sin (2 * pi * Y); +%! Y = Y + amp * sin (2 * pi * X); +%!endfunction + + +%!function assertWarningId (f, expectedId) +%! lastwarn (""); +%! f (); +%! [warningMessage, warningId] = lastwarn (); +%! +%! assert (strcmp (warningId, expectedId), ... +%! sprintf ("Expected warning %s, received %s", ... +%! expectedId, warningId)); +%!endfunction + + +%!function assertNoWarning (f) +%! lastwarn (""); +%! f (); +%! [warningMessage, warningId] = lastwarn (); +%! +%! assert (isempty (warningMessage) && isempty (warningId), ... +%! "An unexpected warning was generated"); +%!endfunction + + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! [X, Y] = sineGrid (21, 0.10); +%! assertNoWarning (@() jacobian2D (2, X, Y)); + + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! [X, Y] = sineGrid (21, 0.10); +%! assertWarningId (@() jacobian2D (2, X.', Y.'), ... +%! "jacobian2D:leftHandedGrid"); + + +%!test +%! % gridGen adds grid curve directories relative to src/octave. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! origDir = pwd (); +%! cleanupDir = onCleanup (@() cd (origDir)); +%! +%! srcDir = fullfile (pwd (), "..", "..", "src", "octave"); +%! addpath (srcDir, "-begin"); +%! cd (srcDir); +%! +%! [X, Y] = gridGen ("TFI", "chevron", 39, 99, false); +%! +%! assert (! isempty (strfind (which ("right"), "chevron")), ... +%! "The chevron curve functions are not on the path"); +%! +%! assertWarningId (@() jacobian2D (2, X, Y), ... +%! "jacobian2D:leftHandedGrid"); +%! +%! assertNoWarning (@() jacobian2D (2, X.', Y.')); + + +%!test +%! % The two orientations have identical absolute Jacobian magnitudes. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! [X, Y] = sineGrid (21, 0.10); +%! +%! warning ("off", "jacobian2D:leftHandedGrid"); +%! restoreWarn = onCleanup ( ... +%! @() warning ("on", "jacobian2D:leftHandedGrid")); +%! +%! a = sort (abs (nonzeros (jacobian2D (2, X, Y)))); +%! b = sort (abs (nonzeros (jacobian2D (2, X.', Y.')))); +%! +%! assert (all (abs (a - b) <= 1e-12), ... +%! "Jacobian magnitudes differ between orientations"); + + +%!test +%! % A folded mesh has mixed-sign Jacobians. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! [X, Y] = meshgrid (linspace (0, 1, 15), ... +%! linspace (0, 1, 15)); +%! +%! X (6:10, 6:10) = X (6:10, 6:10) - 0.45; +%! +%! assertWarningId (@() jacobian2D (2, X, Y), ... +%! "jacobian2D:tangledGrid"); + + +%!test +%! % Warnings propagate through grad2DCurv and div2DCurv. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! [X, Y] = sineGrid (21, 0.10); +%! +%! assertWarningId (@() grad2DCurv (2, X.', Y.'), ... +%! "jacobian2D:leftHandedGrid"); +%! +%! assertWarningId (@() div2DCurv (2, X.', Y.'), ... +%! "jacobian2D:leftHandedGrid"); + + +%!test +%! % Check the five-argument jacobian2D path. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! N = 21; +%! dc = [1; 1; 1; 1]; +%! nc = [0; 0; 0; 0]; +%! +%! s = linspace (0, 1, N); +%! s = [0, s(1:end-1) + 0.5 / (N - 1), 1]; +%! +%! [X, Y] = meshgrid (s, s); +%! X = X + 0.08 * sin (2 * pi * Y); +%! Y = Y + 0.08 * sin (2 * pi * X); +%! +%! assertNoWarning (@() jacobian2D (2, X, Y, dc, nc)); +%! +%! assertWarningId (@() jacobian2D (2, X.', Y.', dc, nc), ... +%! "jacobian2D:leftHandedGrid"); From 82d9f5d11c36d430c98668036640f23aa3e5d62c Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Mon, 19 Oct 2026 11:53:26 -0700 Subject: [PATCH 12/17] Add an individual Octave testLapPolynomial --- tests/octave/Makefile | 3 +- tests/octave/testLapPolynomial.m | 276 ++++++++++++++++--------------- 2 files changed, 144 insertions(+), 135 deletions(-) diff --git a/tests/octave/Makefile b/tests/octave/Makefile index d9c913836..16b203e39 100644 --- a/tests/octave/Makefile +++ b/tests/octave/Makefile @@ -10,7 +10,8 @@ run_octave_tests: testGI13Curv \ testGradPolynomial \ testGridOrientation3D \ - testGridOrientation; do \ + testGridOrientation \ + testLapPolynomial; do \ octave --quiet --eval "addpath(pwd()); test('$$test_name')" || exit 1; \ done diff --git a/tests/octave/testLapPolynomial.m b/tests/octave/testLapPolynomial.m index 68151bfbf..c7f53bcb5 100644 --- a/tests/octave/testLapPolynomial.m +++ b/tests/octave/testLapPolynomial.m @@ -1,134 +1,142 @@ -classdef testLapPolynomial < matlab.unittest.TestCase - methods (Test) - function testLapOfPolynomial1D(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - ks = [2, 4, 6, 8]; % Different orders of accuracy - tol = 1e-10; - x1 = 0; - x2 = 1; - - for k = ks - m = 2 * k + 1; - dx = (x2-x1) / m; - - L = lap(k, m, dx); - - sfield = [x1, x1+dx/2:dx:x2-dx/2, x2]'; - - approx = L * sfield.^k; - analytic = k * (k-1) * sfield.^(k-2); - - % Note: ignore boundary points (because of divergence) - err = approx(2:end-1) - analytic(2:end-1); - - testCase.verifyLessThan(norm(err), tol, ... - sprintf("Nullity test failed for k = %d", k)); - end - end - - function testLapOfPolynomial2D(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - ks = [2, 4, 6, 8]; % Different orders of accuracy - tol = 1e-10; - x1 = 0; - x2 = 1; - y1 = 0; - y2 = 1; - - for k = ks - m = 2 * k + 1; - dx = (x2 - x1) / m; - n = m + 2; - dy = (y2 - y1) / n; - - L = lap2D(k, m, dx, n, dy); - - % staggered grid components (only centers needed) - centers_x = (x1+dx/2:dx:x2-dx/2)'; - centers_y = (y1+dy/2:dy:y2-dy/2)'; - - % set up scalar field only (vector field not needed) - [sfield_y, sfield_x] = meshgrid([y1,centers_y',y2]', [x1,centers_x',x2]'); - sfield_x = reshape(sfield_x, [], 1); - sfield_y = reshape(sfield_y, [], 1); - - % compute approximate and analytical solutions - approx = L * (sfield_x.^k + sfield_y.^k); - analytic = k * (k-1) * sfield_x.^(k-2) + ... - k * (k-1) * sfield_y.^(k-2); - - % remove boundary (because of divergence); int(erior) only - analytic = reshape(analytic, m+2, n+2); - analytic_int = zeros(m+2, n+2); - analytic_int(2:end-1, 2:end-1) = analytic(2:end-1, 2:end-1); - analytic_int = reshape(analytic_int, [], 1); - - err = approx - analytic_int; - - testCase.verifyLessThan(norm(err), tol, ... - sprintf("Nullity test failed for k = %d", k)); - end - end - - function testLapOfPolynomial3D(testCase) - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - addpath(genpath('../../src/octave')) - - ks = [2, 4, 6, 8]; % Different orders of accuracy - tol = 1e-10; - x1 = 0; - x2 = 1; - y1 = 0; - y2 = 1; - z1 = 0; - z2 = 1; - - for k = ks - m = 2 * k + 1; - dx = (x2 - x1) / m; - n = m + 2; - dy = (y2 - y1) / n; - o = n + 2; - dz = (z2 - z1) / o; - - L = lap3D(k, m, dx, n, dy, o, dz); - - % staggered grid components (only centers needed) - centers_x = (x1+dx/2:dx:x2-dx/2)'; - centers_y = (y1+dy/2:dy:y2-dy/2)'; - centers_z = (z1+dz/2:dz:z2-dz/2)'; - - % set up scalar field only (vector field not needed) - [sfield_y, sfield_x, sfield_z] = meshgrid([y1,centers_y',y2]', ... - [x1,centers_x',x2]', [z1,centers_z',z2]'); - sfield_x = reshape(sfield_x, [], 1); - sfield_y = reshape(sfield_y, [], 1); - sfield_z = reshape(sfield_z, [], 1); - - % compute approximate and analytical solutions - approx = L * (sfield_x.^k + sfield_y.^k + sfield_z.^k); - analytic = k * (k-1) * sfield_x.^(k-2) + ... - k * (k-1) * sfield_y.^(k-2) + ... - k * (k-1) * sfield_z.^(k-2); - - % remove boundary (because of divergence); int(erior) only - analytic = reshape(analytic, m+2, n+2, o+2); - analytic_int = zeros(m+2, n+2, o+2); - analytic_int(2:end-1, 2:end-1, 2:end-1) = analytic(2:end-1, 2:end-1, 2:end-1); - analytic_int = reshape(analytic_int, [], 1); - - err = approx - analytic_int; - - testCase.verifyLessThan(norm(err), tol, ... - sprintf("Nullity test failed for k = %d", k)); - end - end - end -end \ No newline at end of file +function testLapPolynomial () +endfunction + + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! ks = [2, 4, 6, 8]; +%! tol = 1e-10; +%! x1 = 0; +%! x2 = 1; +%! +%! for k = ks +%! m = 2 * k + 1; +%! dx = (x2 - x1) / m; +%! +%! L = lap (k, m, dx); +%! +%! sfield = [x1, x1 + dx/2:dx:x2 - dx/2, x2]'; +%! +%! approx = L * sfield .^ k; +%! analytic = k * (k - 1) * sfield .^ (k - 2); +%! +%! % Ignore boundary points. +%! err = approx (2:end-1) - analytic (2:end-1); +%! +%! assert (norm (err) < tol, ... +%! sprintf ("1D Laplacian test failed for k = %d", k)); +%! endfor + + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! ks = [2, 4, 6, 8]; +%! tol = 1e-10; +%! x1 = 0; +%! x2 = 1; +%! y1 = 0; +%! y2 = 1; +%! +%! for k = ks +%! m = 2 * k + 1; +%! dx = (x2 - x1) / m; +%! n = m + 2; +%! dy = (y2 - y1) / n; +%! +%! L = lap2D (k, m, dx, n, dy); +%! +%! centers_x = (x1 + dx/2:dx:x2 - dx/2)'; +%! centers_y = (y1 + dy/2:dy:y2 - dy/2)'; +%! +%! [sfield_y, sfield_x] = meshgrid ( ... +%! [y1, centers_y', y2]', ... +%! [x1, centers_x', x2]'); +%! +%! sfield_x = reshape (sfield_x, [], 1); +%! sfield_y = reshape (sfield_y, [], 1); +%! +%! approx = L * (sfield_x .^ k + sfield_y .^ k); +%! +%! analytic = ... +%! k * (k - 1) * sfield_x .^ (k - 2) + ... +%! k * (k - 1) * sfield_y .^ (k - 2); +%! +%! analytic = reshape (analytic, m + 2, n + 2); +%! +%! analytic_int = zeros (m + 2, n + 2); +%! analytic_int (2:end-1, 2:end-1) = ... +%! analytic (2:end-1, 2:end-1); +%! +%! analytic_int = reshape (analytic_int, [], 1); +%! err = approx - analytic_int; +%! +%! assert (norm (err) < tol, ... +%! sprintf ("2D Laplacian test failed for k = %d", k)); +%! endfor + + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! ks = [2, 4, 6, 8]; +%! tol = 1e-10; +%! x1 = 0; +%! x2 = 1; +%! y1 = 0; +%! y2 = 1; +%! z1 = 0; +%! z2 = 1; +%! +%! for k = ks +%! m = 2 * k + 1; +%! dx = (x2 - x1) / m; +%! n = m + 2; +%! dy = (y2 - y1) / n; +%! o = n + 2; +%! dz = (z2 - z1) / o; +%! +%! L = lap3D (k, m, dx, n, dy, o, dz); +%! +%! centers_x = (x1 + dx/2:dx:x2 - dx/2)'; +%! centers_y = (y1 + dy/2:dy:y2 - dy/2)'; +%! centers_z = (z1 + dz/2:dz:z2 - dz/2)'; +%! +%! [sfield_y, sfield_x, sfield_z] = meshgrid ( ... +%! [y1, centers_y', y2]', ... +%! [x1, centers_x', x2]', ... +%! [z1, centers_z', z2]'); +%! +%! sfield_x = reshape (sfield_x, [], 1); +%! sfield_y = reshape (sfield_y, [], 1); +%! sfield_z = reshape (sfield_z, [], 1); +%! +%! approx = L * ( ... +%! sfield_x .^ k + ... +%! sfield_y .^ k + ... +%! sfield_z .^ k); +%! +%! analytic = ... +%! k * (k - 1) * sfield_x .^ (k - 2) + ... +%! k * (k - 1) * sfield_y .^ (k - 2) + ... +%! k * (k - 1) * sfield_z .^ (k - 2); +%! +%! analytic = reshape (analytic, m + 2, n + 2, o + 2); +%! +%! analytic_int = zeros (m + 2, n + 2, o + 2); +%! analytic_int (2:end-1, 2:end-1, 2:end-1) = ... +%! analytic (2:end-1, 2:end-1, 2:end-1); +%! +%! analytic_int = reshape (analytic_int, [], 1); +%! err = approx - analytic_int; +%! +%! assert (norm (err) < tol, ... +%! sprintf ("3D Laplacian test failed for k = %d", k)); +%! endfor From 6dc36182d9ae195790ee9f5b88bda5060b75a359 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Mon, 19 Oct 2026 11:59:00 -0700 Subject: [PATCH 13/17] Add an individual Octave testPoissonAccuracy --- tests/octave/Makefile | 3 +- tests/octave/testPoissonAccuracy.m | 116 +++++++++++++++-------------- 2 files changed, 62 insertions(+), 57 deletions(-) diff --git a/tests/octave/Makefile b/tests/octave/Makefile index 16b203e39..237d2cc89 100644 --- a/tests/octave/Makefile +++ b/tests/octave/Makefile @@ -11,7 +11,8 @@ run_octave_tests: testGradPolynomial \ testGridOrientation3D \ testGridOrientation \ - testLapPolynomial; do \ + testLapPolynomial \ + testPoissonAccuracy; do \ octave --quiet --eval "addpath(pwd()); test('$$test_name')" || exit 1; \ done diff --git a/tests/octave/testPoissonAccuracy.m b/tests/octave/testPoissonAccuracy.m index 5026236b2..ed8b22d58 100644 --- a/tests/octave/testPoissonAccuracy.m +++ b/tests/octave/testPoissonAccuracy.m @@ -1,56 +1,60 @@ -classdef testPoissonAccuracy < matlab.unittest.TestCase - methods(Test) - function testforEnergy(testCase) - - west = 0; % Domain's limits - east = 1; - - ks = [2, 4, 6]; % Different orders of accuracy - % tests for ks = 8 removed for issues - grid_sizes = [20, 40]; % Grid sizes to test - - for k = ks - - errors = zeros(size(grid_sizes)); - - for i = 1:numel(grid_sizes) - m = grid_sizes(i); % Number of cells - dx = (east - west) / m; % Step length - - L = lap(k, m, dx); % 1D Mimetic Laplacian operator - - % Impose Robin BC on Laplacian operator - a = 1; - b = 1; - L = L + robinBC(k, m, dx, a, b); - - % 1D Staggered grid - grid = [west west+dx/2 : dx : east-dx/2 east]; - - % RHS - U = exp(grid)'; - U(1) = 0; % West BC - U(end) = 2*exp(1); % East BC - - % Solve a linear system of equations - computed_solution = L\U; - - % Compute error using L2 norm - analytical_solution = exp(grid); - errors(i) = max(abs(computed_solution' - analytical_solution)); - end - - % Compute order of accuracy - order = zeros(numel(errors) - 1, 1); - for i = 1:numel(errors) - 1 - order(i) = log2(errors(i) / errors(i + 1)); - - testCase.verifyGreaterThan(order(i),k-0.5, ... - sprintf("Test failed for k = %d", k)); - end - end - end - end -end - - +function testPoissonAccuracy () +endfunction + + +%!test +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! addpath (genpath ("../../src/octave"), "-begin"); +%! +%! west = 0; +%! east = 1; +%! +%! ks = [2, 4, 6]; +%! grid_sizes = [20, 40]; +%! +%! for k = ks +%! errors = zeros (size (grid_sizes)); +%! +%! for i = 1:numel (grid_sizes) +%! m = grid_sizes (i); +%! dx = (east - west) / m; +%! +%! % Mimetic Laplacian operator. +%! L = lap (k, m, dx); +%! +%! % Impose Robin boundary conditions. +%! a = 1; +%! b = 1; +%! L = L + robinBC (k, m, dx, a, b); +%! +%! % One-dimensional staggered grid. +%! grid = [ ... +%! west, ... +%! west + dx/2:dx:east - dx/2, ... +%! east]; +%! +%! % Right-hand side. +%! U = exp (grid)'; +%! U (1) = 0; +%! U (end) = 2 * exp (1); +%! +%! % Solve the linear system. +%! computed_solution = L \ U; +%! +%! % Compute the maximum error. +%! analytical_solution = exp (grid); +%! errors (i) = max ( ... +%! abs (computed_solution' - analytical_solution)); +%! endfor +%! +%! % Compute the observed order of accuracy. +%! order = zeros (numel (errors) - 1, 1); +%! +%! for i = 1:numel (errors) - 1 +%! order (i) = log2 (errors (i) / errors (i + 1)); +%! +%! assert (order (i) > k - 0.5, ... +%! sprintf ("Poisson accuracy test failed for k = %d", k)); +%! endfor +%! endfor From 773842b293211cdf3af6fa5c61b5ebe92e1ae940 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Mon, 19 Oct 2026 12:03:31 -0700 Subject: [PATCH 14/17] Add an individual Octave testTTMGrids --- tests/octave/Makefile | 3 +- tests/octave/testTTMGrids.m | 211 ++++++++++++++++++------------------ 2 files changed, 107 insertions(+), 107 deletions(-) diff --git a/tests/octave/Makefile b/tests/octave/Makefile index 237d2cc89..cd8ce92f0 100644 --- a/tests/octave/Makefile +++ b/tests/octave/Makefile @@ -12,7 +12,8 @@ run_octave_tests: testGridOrientation3D \ testGridOrientation \ testLapPolynomial \ - testPoissonAccuracy; do \ + testPoissonAccuracy \ + testTTMGrids; do \ octave --quiet --eval "addpath(pwd()); test('$$test_name')" || exit 1; \ done diff --git a/tests/octave/testTTMGrids.m b/tests/octave/testTTMGrids.m index f4a3bdc69..e10827ebf 100644 --- a/tests/octave/testTTMGrids.m +++ b/tests/octave/testTTMGrids.m @@ -1,116 +1,115 @@ -classdef testTTMGrids < matlab.unittest.TestCase -% Tests that ttm produces unfolded grids. -% -% Before the TFI initial guess was added, ttm started with the interior at -% zeros(m,n) -- every interior node at the origin -- so the first sweeps -% computed alpha, beta and gamma from a maximally degenerate grid. The -% iteration then had to untangle a fold it had created itself, and on every -% grid shipped with MOLE it failed to: swan, chevron and horseshoe all came -% out with mixed-sign Jacobians at 500 iterations. -% -% A folded grid is locally left-handed. jacobian2D returns the signed -% determinant and grad2DCurv divides by it, so a run on one produces finite, -% plausible output with the metric tensor inverted in a patch. Nothing -% raises. -% -% NOTE ON ORIENTATION. ttm returns X with xi-nodes as rows, which is the -% opposite of what jacobian2D reads, so these tests transpose. Without the -% transpose a correct grid reports ALL NEGATIVE rather than ALL POSITIVE. +function testTTMGrids () +endfunction - methods(Static) - function [npos, nneg] = jacobianSigns(X, Y) - origState = warning('off', 'jacobian2D:leftHandedGrid'); - restore1 = onCleanup(@() warning(origState)); - origState2 = warning('off', 'jacobian2D:tangledGrid'); - restore2 = onCleanup(@() warning(origState2)); - d = nonzeros(jacobian2D(2, X.', Y.')); - npos = sum(d > 0); - nneg = sum(d < 0); - end - end - methods(Test) +%!function [npos, nneg] = jacobianSigns (X, Y) +%! state1 = warning ("off", "jacobian2D:leftHandedGrid"); +%! restore1 = onCleanup (@() warning (state1)); +%! +%! state2 = warning ("off", "jacobian2D:tangledGrid"); +%! restore2 = onCleanup (@() warning (state2)); +%! +%! d = nonzeros (jacobian2D (2, X.', Y.')); +%! npos = sum (d > 0); +%! nneg = sum (d < 0); +%!endfunction - function testShippedGridsAreNotFolded(testCase) - % chevron and horseshoe are clean under TFI, so there is no excuse - % for the elliptic generator to fold them. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - origDir = pwd; - cleanupDir = onCleanup(@() cd(origDir)); - addpath(genpath('../../src/octave')) - cd('../../src/octave') - for g = {'chevron', 'horseshoe'} - [X, Y] = ttm(g{1}, 39, 99, 500, false); - [~, nneg] = testTTMGrids.jacobianSigns(X, Y); - testCase.verifyEqual(nneg, 0, ... - sprintf('ttm folded the %s grid: %d negative Jacobian entries', ... - g{1}, nneg)); - end - end +%!test +%! % Chevron and horseshoe should not be folded by TTM. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! origDir = pwd (); +%! cleanupDir = onCleanup (@() cd (origDir)); +%! +%! srcDir = fullfile (pwd (), "..", "..", "src", "octave"); +%! addpath (srcDir, "-begin"); +%! cd (srcDir); +%! +%! grids = {"chevron", "horseshoe"}; +%! +%! for i = 1:numel (grids) +%! [X, Y] = ttm (grids{i}, 39, 99, 500, false); +%! [~, nneg] = jacobianSigns (X, Y); +%! +%! assert (nneg == 0, ... +%! sprintf ( ... +%! "ttm folded the %s grid: %d negative Jacobian entries", ... +%! grids{i}, nneg)); +%! endfor - function testEllipticSmoothingRepairsAFoldedTFI(testCase) - % swan is folded under TFI itself (113 negative entries). Repairing - % that is exactly what an elliptic generator is for, and with a - % sensible initial guess ttm does it. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - origDir = pwd; - cleanupDir = onCleanup(@() cd(origDir)); - addpath(genpath('../../src/octave')) - cd('../../src/octave') - [Xt, Yt] = tfi('swan', 39, 99, false); - [~, nnegTFI] = testTTMGrids.jacobianSigns(Xt, Yt); - testCase.verifyGreaterThan(nnegTFI, 0, ... - 'swan is expected to be folded under TFI; test premise changed'); +%!test +%! % Elliptic smoothing should repair the folded swan TFI grid. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! origDir = pwd (); +%! cleanupDir = onCleanup (@() cd (origDir)); +%! +%! srcDir = fullfile (pwd (), "..", "..", "src", "octave"); +%! addpath (srcDir, "-begin"); +%! cd (srcDir); +%! +%! [Xt, Yt] = tfi ("swan", 39, 99, false); +%! [~, nnegTFI] = jacobianSigns (Xt, Yt); +%! +%! assert (nnegTFI > 0, ... +%! "swan is expected to be folded under TFI; test premise changed"); +%! +%! [X, Y] = ttm ("swan", 39, 99, 500, false); +%! [~, nneg] = jacobianSigns (X, Y); +%! +%! assert (nneg == 0, ... +%! sprintf ("ttm left swan folded: %d negative entries", nneg)); - [X, Y] = ttm('swan', 39, 99, 500, false); - [~, nneg] = testTTMGrids.jacobianSigns(X, Y); - testCase.verifyEqual(nneg, 0, ... - sprintf('ttm left swan folded: %d negative entries', nneg)); - end - function testFewIterationsDoNotFold(testCase) - % The failure was not slow degradation from over-iteration -- the - % stock version was folded within 5 sweeps. Check the low end. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - origDir = pwd; - cleanupDir = onCleanup(@() cd(origDir)); - addpath(genpath('../../src/octave')) - cd('../../src/octave') +%!test +%! % The grid should remain unfolded after only a few iterations. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! origDir = pwd (); +%! cleanupDir = onCleanup (@() cd (origDir)); +%! +%! srcDir = fullfile (pwd (), "..", "..", "src", "octave"); +%! addpath (srcDir, "-begin"); +%! cd (srcDir); +%! +%! iterations = [5, 25, 100]; +%! +%! for i = 1:numel (iterations) +%! it = iterations (i); +%! [X, Y] = ttm ("chevron", 39, 99, it, false); +%! [~, nneg] = jacobianSigns (X, Y); +%! +%! assert (nneg == 0, ... +%! sprintf ( ... +%! "ttm folded chevron after %d iterations (%d negative)", ... +%! it, nneg)); +%! endfor - for it = [5 25 100] - [X, Y] = ttm('chevron', 39, 99, it, false); - [~, nneg] = testTTMGrids.jacobianSigns(X, Y); - testCase.verifyEqual(nneg, 0, ... - sprintf('ttm folded chevron after %d iterations (%d negative)', ... - it, nneg)); - end - end - function testBoundariesAreUnchanged(testCase) - % The initial guess must only touch the interior. Boundary nodes - % come from the curve definitions and must survive untouched. - origPath = path; - cleanupObj = onCleanup(@() path(origPath)); - origDir = pwd; - cleanupDir = onCleanup(@() cd(origDir)); - addpath(genpath('../../src/octave')) - cd('../../src/octave') - - m = 39; - n = 99; - [X, Y] = ttm('chevron', m, n, 50, false); - [Xt, Yt] = tfi('chevron', m, n, false); - - testCase.verifyEqual(X(1, :), Xt(1, :), 'AbsTol', 1e-12); - testCase.verifyEqual(X(end, :), Xt(end, :), 'AbsTol', 1e-12); - testCase.verifyEqual(Y(:, 1), Yt(:, 1), 'AbsTol', 1e-12); - testCase.verifyEqual(Y(:, end), Yt(:, end), 'AbsTol', 1e-12); - end - - end -end +%!test +%! % TTM should not modify boundary nodes. +%! origPath = path (); +%! cleanupObj = onCleanup (@() path (origPath)); +%! origDir = pwd (); +%! cleanupDir = onCleanup (@() cd (origDir)); +%! +%! srcDir = fullfile (pwd (), "..", "..", "src", "octave"); +%! addpath (srcDir, "-begin"); +%! cd (srcDir); +%! +%! m = 39; +%! n = 99; +%! +%! [X, Y] = ttm ("chevron", m, n, 50, false); +%! [Xt, Yt] = tfi ("chevron", m, n, false); +%! +%! assert (all (abs (X (1, :) - Xt (1, :)) <= 1e-12), ... +%! "Lower X boundary changed"); +%! assert (all (abs (X (end, :) - Xt (end, :)) <= 1e-12), ... +%! "Upper X boundary changed"); +%! assert (all (abs (Y (:, 1) - Yt (:, 1)) <= 1e-12), ... +%! "Left Y boundary changed"); +%! assert (all (abs (Y (:, end) - Yt (:, end)) <= 1e-12), ... +%! "Right Y boundary changed"); From 316671d7a234bcc541111140f50c0f2d1b0e7685 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Mon, 19 Oct 2026 12:27:23 -0700 Subject: [PATCH 15/17] Add Octave test runner invocation --- .github/workflows/matlab-tests.yml | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/.github/workflows/matlab-tests.yml b/.github/workflows/matlab-tests.yml index 2d88f6d60..b6cddc938 100644 --- a/.github/workflows/matlab-tests.yml +++ b/.github/workflows/matlab-tests.yml @@ -54,3 +54,22 @@ jobs: with: name: matlab-code-coverage path: code-coverage/coverage.xml + + - name: Checkout code + uses: actions/checkout@v4 + with: + # Full git history is needed to get a proper + # list of changed files within `super-linter` + fetch-depth: 0 + + - name: Install system deps + run: | + wget -O - https://apt.llvm.org/llvm-snapshot.gpg.key | sudo apt-key add - + sudo add-apt-repository "deb http://apt.llvm.org/noble/ llvm-toolchain-noble-20 main" + sudo apt-get update + sudo apt-get install -y octave + + - name: Run Octave tests + run: | + cd tests + octave --quiet --eval "oruntests('octave')" From dd13226c4e63d15c5df69c9ebf60c4816b5cd401 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Mon, 19 Oct 2026 12:43:05 -0700 Subject: [PATCH 16/17] Add Octave test runner invocation, bug fix --- .github/workflows/matlab-tests.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/matlab-tests.yml b/.github/workflows/matlab-tests.yml index b6cddc938..d01ced943 100644 --- a/.github/workflows/matlab-tests.yml +++ b/.github/workflows/matlab-tests.yml @@ -71,5 +71,5 @@ jobs: - name: Run Octave tests run: | - cd tests - octave --quiet --eval "oruntests('octave')" + cd tests/octave + octave --quiet --eval "oruntests('.')" From d4e917fe836ad7d423e55439aa33e609fed9d598 Mon Sep 17 00:00:00 2001 From: Christopher Paolini Date: Mon, 26 Oct 2026 19:23:13 +0000 Subject: [PATCH 17/17] Enhance README with installation and test instructions Updated installation instructions for Ubuntu/Debian and added verification steps for Octave. Signed-off-by: Christopher Paolini --- README.md | 33 ++++++++++++++++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index e7ea109b7..43ec5567d 100644 --- a/README.md +++ b/README.md @@ -48,7 +48,8 @@ For documentation build requirements, please refer to the [Documentation Guide]( ```bash # Install all required packages -sudo apt install cmake libopenblas-dev libeigen3-dev +sudo apt update +sudo apt install octave make cmake libopenblas-dev libeigen3-dev ``` #### macOS Systems @@ -125,10 +126,40 @@ make run_tests Octave equivalent of the C++ test suite. We recommend running these tests before using MOLE to ensure proper setup. +First, verify the installation: + +```bash +octave --version +make --version +``` + +GNU Octave is the supported open-source environment for running the tests. MATLAB is not required. The general contribution guidelines list GNU Octave as an accepted environment for MATLAB/Octave development. + ```matlab +cd tests/octave make run_octave_tests ``` +The test suite includes tests for: +- 1D, 2D, and 3D gradient operators; +- 1D, 2D, and 3D divergence operators; +- 1D, 2D, and 3D Laplacian operators; +- polynomial exactness and convergence; +- boundary-condition operators; +- Poisson-solution accuracy; +- curvilinear-grid operators; +- Jacobian orientation and degeneracy detection; and +- transfinite and elliptic grid generation. + +If GNU Make is unavailable, run the complete test directory directly: + +```bash +cd tests/octave +octave --quiet --eval "addpath(pwd()); oruntests(pwd())" +``` + +A successful test produces no test failed message. A failed test reports the test block and the associated assertion message. + ## Examples Many of the examples require 'gnuplot' to visualize the results. You can get gnuplot on macOSX with