From 46f8f36b6e7354973bcb2e318b0ce8dfe950abcd Mon Sep 17 00:00:00 2001 From: Ben Wagner Date: Tue, 20 Oct 2026 12:58:52 -0700 Subject: [PATCH 1/3] addscalarbc zeroRows fix --- src/cpp/addscalarbc.cpp | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/cpp/addscalarbc.cpp b/src/cpp/addscalarbc.cpp index 9bb4affd5..cf089212a 100644 --- a/src/cpp/addscalarbc.cpp +++ b/src/cpp/addscalarbc.cpp @@ -51,13 +51,10 @@ uvec collectUniqueRows(const sp_mat &M) { * replacing them with boundary condition equations. */ void zeroRows(sp_mat &A, const uvec &rows) { - for (uword i = 0; i < rows.n_elem; i++) { - uword row_idx = rows(i); - for (sp_mat::const_row_iterator it = A.begin_row(row_idx); - it != A.end_row(row_idx); ++it) { - A(it.row(), it.col()) = 0; - } + for (uword idx : rows) { + A.row(idx).zeros(); } + A.sync(); } /** @@ -362,4 +359,4 @@ void addScalarBC(sp_mat &A, vec &b, u16 k, u32 m, Real dx, addScalarBCrhs(b, bc.dc, bc.nc, bc.v, rl, rr, rb, rt, rf, rk); } -} // namespace AddScalarBC \ No newline at end of file +} // namespace AddScalarBC From 1fc0e374e908e6086e976964f75791d6410fe2e2 Mon Sep 17 00:00:00 2001 From: Ben Wagner Date: Tue, 20 Oct 2026 12:59:29 -0700 Subject: [PATCH 2/3] test addscalarbc update --- tests/cpp/test_addscalarbc.cpp | 98 +++++++++++++++++++++------------- 1 file changed, 61 insertions(+), 37 deletions(-) diff --git a/tests/cpp/test_addscalarbc.cpp b/tests/cpp/test_addscalarbc.cpp index 5b0a171a6..f99d4f859 100644 --- a/tests/cpp/test_addscalarbc.cpp +++ b/tests/cpp/test_addscalarbc.cpp @@ -122,11 +122,13 @@ void test_1D_dirichlet() { const u32 m = 10; // number of interior points (ghost points add +2 in this codebase) const Real dx = 0.1; // grid spacing - Laplacian L(k, m, dx); - sp_mat A = sp_mat(L); + const uword ncenters = m + 2; + const Real density = 0.1; + + sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; - vec b(m + 2, fill::ones); + vec b(ncenters, fill::ones); vec b_orig = b; BC1D bc; @@ -165,11 +167,13 @@ void test_1D_neumann() { const u32 m = 10; const Real dx = 0.1; - Laplacian L(k, m, dx); - sp_mat A = sp_mat(L); + const uword ncenters = m + 2; + const Real density = 0.1; + + sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; - vec b(m + 2, fill::ones); + vec b(ncenters, fill::ones); BC1D bc; bc.dc = {0.0, 0.0}; // no Dirichlet component @@ -206,11 +210,13 @@ void test_1D_mixed_dirichlet_neumann() { const u32 m = 10; const Real dx = 0.1; - Laplacian L(k, m, dx); - sp_mat A = sp_mat(L); + const uword ncenters = m + 2; + const Real density = 0.1; + + sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; - vec b(m + 2, fill::ones); + vec b(ncenters, fill::ones); BC1D bc; bc.dc = {1.0, 0.0}; // Dirichlet on left only @@ -246,11 +252,13 @@ void test_1D_mixed_neumann_dirichlet() { const u32 m = 10; const Real dx = 0.1; - Laplacian L(k, m, dx); - sp_mat A = sp_mat(L); + const uword ncenters = m + 2; + const Real density = 0.1; + + sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; - vec b(m + 2, fill::ones); + vec b(ncenters, fill::ones); BC1D bc; bc.dc = {0.0, 1.0}; // Dirichlet on right only @@ -292,11 +300,13 @@ void test_1D_robin() { const u32 m = 10; const Real dx = 0.1; - Laplacian L(k, m, dx); - sp_mat A = sp_mat(L); + const uword ncenters = m + 2; + const Real density = 0.1; + + sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; - vec b(m + 2, fill::ones); + vec b(ncenters, fill::ones); BC1D bc; bc.dc = {2.0, 3.0}; // value term weights @@ -336,11 +346,13 @@ void test_2D_dirichlet_all() { const Real dx = 0.1; const Real dy = 0.2; - Laplacian L(k, m, n, dx, dy); - sp_mat A = sp_mat(L); + const uword ncenters = (m + 2) * (n + 2); + const Real density = 0.1; + + sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; - vec b((m + 2) * (n + 2), fill::ones); + vec b(ncenters, fill::ones); BC2D bc; bc.dc = {1.0, 1.0, 1.0, 1.0}; // Left, Right, Bottom, Top => Dirichlet @@ -378,11 +390,13 @@ void test_2D_neumann_all() { const Real dx = 0.1; const Real dy = 0.2; - Laplacian L(k, m, n, dx, dy); - sp_mat A = sp_mat(L); + const uword ncenters = (m + 2) * (n + 2); + const Real density = 0.1; + + sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; - vec b((m + 2) * (n + 2), fill::ones); + vec b(ncenters, fill::ones); BC2D bc; bc.dc = {0.0, 0.0, 0.0, 0.0}; @@ -419,11 +433,13 @@ void test_2D_mixed() { const Real dx = 0.1; const Real dy = 0.1; - Laplacian L(k, m, n, dx, dy); - sp_mat A = sp_mat(L); + const uword ncenters = (m + 2) * (n + 2); + const Real density = 0.1; + + sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; - vec b((m + 2) * (n + 2), fill::ones); + vec b(ncenters, fill::ones); BC2D bc; bc.dc = {1.0, 1.0, 0.0, 0.0}; // Dirichlet L/R @@ -463,11 +479,13 @@ void test_3D_dirichlet_all() { const Real dy = 0.1; const Real dz = 0.1; - Laplacian L(k, m, n, o, dx, dy, dz); - sp_mat A = sp_mat(L); + const uword ncenters = (m + 2) * (n + 2) * (o + 2); + const Real density = 0.1; + + sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; - vec b((m + 2) * (n + 2) * (o + 2), fill::ones); + vec b(ncenters, fill::ones); BC3D bc; bc.dc = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; @@ -508,11 +526,13 @@ void test_3D_neumann_all() { const Real dy = 0.1; const Real dz = 0.1; - Laplacian L(k, m, n, o, dx, dy, dz); - sp_mat A = sp_mat(L); + const uword ncenters = (m + 2) * (n + 2) * (o + 2); + const Real density = 0.1; + + sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; - vec b((m + 2) * (n + 2) * (o + 2), fill::ones); + vec b(ncenters, fill::ones); BC3D bc; bc.dc = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; @@ -550,11 +570,13 @@ void test_3D_mixed() { const Real dy = 0.1; const Real dz = 0.1; - Laplacian L(k, m, n, o, dx, dy, dz); - sp_mat A = sp_mat(L); + const uword ncenters = (m + 2) * (n + 2) * (o + 2); + const Real density = 0.1; + + sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; - vec b((m + 2) * (n + 2) * (o + 2), fill::ones); + vec b(ncenters, fill::ones); BC3D bc; bc.dc = {1.0, 1.0, 0.0, 0.0, 0.0, 0.0}; // Dirichlet on x faces @@ -596,11 +618,13 @@ void test_periodic_1D_do_nothing() { const u32 m = 10; const Real dx = 0.1; - Laplacian L(k, m, dx); - sp_mat A_orig = sp_mat(L); + const uword ncenters = m; + const Real density = 0.1; + + sp_mat A_orig = sprandn(ncenters, ncenters, density); sp_mat A = A_orig; - vec b_orig(m + 2, fill::ones); + vec b_orig(ncenters, fill::ones); vec b = b_orig; BC1D bc; @@ -656,4 +680,4 @@ int main() { } return 0; -} \ No newline at end of file +} From 91045090bb104319edf91d3b605c46d6b5c5cd72 Mon Sep 17 00:00:00 2001 From: Ben Wagner Date: Tue, 20 Oct 2026 13:35:02 -0700 Subject: [PATCH 3/3] density tweaks --- tests/cpp/test_addscalarbc.cpp | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/tests/cpp/test_addscalarbc.cpp b/tests/cpp/test_addscalarbc.cpp index f99d4f859..b92f69a04 100644 --- a/tests/cpp/test_addscalarbc.cpp +++ b/tests/cpp/test_addscalarbc.cpp @@ -123,7 +123,7 @@ void test_1D_dirichlet() { const Real dx = 0.1; // grid spacing const uword ncenters = m + 2; - const Real density = 0.1; + const Real density = 1.0 - 2.0 / static_cast(ncenters * ncenters) + 0.01; sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; @@ -168,7 +168,7 @@ void test_1D_neumann() { const Real dx = 0.1; const uword ncenters = m + 2; - const Real density = 0.1; + const Real density = 1.0 - 2.0 / static_cast(ncenters * ncenters) + 0.01; sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; @@ -211,7 +211,7 @@ void test_1D_mixed_dirichlet_neumann() { const Real dx = 0.1; const uword ncenters = m + 2; - const Real density = 0.1; + const Real density = 1.0 - 2.0 / static_cast(ncenters * ncenters) + 0.01; sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; @@ -253,7 +253,7 @@ void test_1D_mixed_neumann_dirichlet() { const Real dx = 0.1; const uword ncenters = m + 2; - const Real density = 0.1; + const Real density = 1.0 - 2.0 / static_cast(ncenters * ncenters) + 0.01; sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; @@ -301,7 +301,7 @@ void test_1D_robin() { const Real dx = 0.1; const uword ncenters = m + 2; - const Real density = 0.1; + const Real density = 1.0 - 2.0 / static_cast(ncenters * ncenters) + 0.01; sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; @@ -347,7 +347,7 @@ void test_2D_dirichlet_all() { const Real dy = 0.2; const uword ncenters = (m + 2) * (n + 2); - const Real density = 0.1; + const Real density = 1.0 - 2.0 / static_cast(ncenters * ncenters); sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; @@ -391,7 +391,7 @@ void test_2D_neumann_all() { const Real dy = 0.2; const uword ncenters = (m + 2) * (n + 2); - const Real density = 0.1; + const Real density = 1.0 - 2.0 / static_cast(ncenters * ncenters); sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; @@ -434,7 +434,7 @@ void test_2D_mixed() { const Real dy = 0.1; const uword ncenters = (m + 2) * (n + 2); - const Real density = 0.1; + const Real density = 1.0 - 2.0 / static_cast(ncenters * ncenters); sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; @@ -480,7 +480,7 @@ void test_3D_dirichlet_all() { const Real dz = 0.1; const uword ncenters = (m + 2) * (n + 2) * (o + 2); - const Real density = 0.1; + const Real density = 1.0 - 2.0 / static_cast(ncenters * ncenters); sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; @@ -527,7 +527,7 @@ void test_3D_neumann_all() { const Real dz = 0.1; const uword ncenters = (m + 2) * (n + 2) * (o + 2); - const Real density = 0.1; + const Real density = 1.0 - 2.0 / static_cast(ncenters * ncenters); sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; @@ -571,7 +571,7 @@ void test_3D_mixed() { const Real dz = 0.1; const uword ncenters = (m + 2) * (n + 2) * (o + 2); - const Real density = 0.1; + const Real density = 1.0 - 2.0 / static_cast(ncenters * ncenters); sp_mat A = sprandn(ncenters, ncenters, density); sp_mat A_orig = A; @@ -619,7 +619,7 @@ void test_periodic_1D_do_nothing() { const Real dx = 0.1; const uword ncenters = m; - const Real density = 0.1; + const Real density = 1.0 - 2.0 / static_cast(ncenters * ncenters) + 0.01; sp_mat A_orig = sprandn(ncenters, ncenters, density); sp_mat A = A_orig;