Skip to content

Cpp addscalarBC bug fix #504

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
bewagner1 wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Cpp addscalarBC bug fix #504

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
46f8f36
addscalarbc zeroRows fix
bewagner1 Oct 20, 2026
1fc0e37
test addscalarbc update
bewagner1 Oct 20, 2026
9104509
density tweaks
bewagner1 Oct 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 4 additions & 7 deletions src/cpp/addscalarbc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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();
}

/**
Expand Down Expand Up @@ -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
} // namespace AddScalarBC
98 changes: 61 additions & 37 deletions tests/cpp/test_addscalarbc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 = 1.0 - 2.0 / static_cast<Real>(ncenters * ncenters) + 0.01;

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;
Expand Down Expand Up @@ -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 = 1.0 - 2.0 / static_cast<Real>(ncenters * ncenters) + 0.01;

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
Expand Down Expand Up @@ -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 = 1.0 - 2.0 / static_cast<Real>(ncenters * ncenters) + 0.01;

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
Expand Down Expand Up @@ -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 = 1.0 - 2.0 / static_cast<Real>(ncenters * ncenters) + 0.01;

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
Expand Down Expand Up @@ -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 = 1.0 - 2.0 / static_cast<Real>(ncenters * ncenters) + 0.01;

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
Expand Down Expand Up @@ -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 = 1.0 - 2.0 / static_cast<Real>(ncenters * ncenters);

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
Expand Down Expand Up @@ -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 = 1.0 - 2.0 / static_cast<Real>(ncenters * ncenters);

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};
Expand Down Expand Up @@ -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 = 1.0 - 2.0 / static_cast<Real>(ncenters * ncenters);

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
Expand Down Expand Up @@ -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 = 1.0 - 2.0 / static_cast<Real>(ncenters * ncenters);

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};
Expand Down Expand Up @@ -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 = 1.0 - 2.0 / static_cast<Real>(ncenters * ncenters);

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};
Expand Down Expand Up @@ -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 = 1.0 - 2.0 / static_cast<Real>(ncenters * ncenters);

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
Expand Down Expand Up @@ -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 = 1.0 - 2.0 / static_cast<Real>(ncenters * ncenters) + 0.01;

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;
Expand Down Expand Up @@ -656,4 +680,4 @@ int main() {
}

return 0;
}
}
Loading