Skip to content

Use zstart/zend in more loops #3244

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

Merged
ZedThree merged 2 commits into from
Jan 9, 2026
Merged

Use zstart/zend in more loops #3244

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0386f67
Use `zstart/zend` in more loops
ZedThree Jan 7, 2026
b20d19d
Merge branch 'next' into more-zstart-zend-loops
ZedThree Jan 9, 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
18 changes: 9 additions & 9 deletions examples/elm-pb/elm_pb.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1366,7 +1366,7 @@ class ELMpb : public PhysicsModel {
// Calculate a single phi boundary value for all Y slices
BoutReal philocal = 0.0;
for (int j = mesh->ystart; j <= mesh->yend; j++) {
for (int k = 0; k < mesh->LocalNz; k++) {
for (int k = mesh->zstart; k <= mesh->zend; k++) {
philocal += phi(mesh->xstart, j, k);
}
}
Expand All @@ -1379,7 +1379,7 @@ class ELMpb : public PhysicsModel {
for (int j = mesh->ystart; j <= mesh->yend; j++) {
if (!phi_core_averagey) {
phivalue = 0.0; // Calculate phi boundary for each Y index separately
for (int k = 0; k < mesh->LocalNz; k++) {
for (int k = mesh->zstart; k <= mesh->zend; k++) {
phivalue += phi(mesh->xstart, j, k);
}
phivalue /= mesh->LocalNz; // Average in Z of point next to boundary
Expand All @@ -1393,7 +1393,7 @@ class ELMpb : public PhysicsModel {
BoutReal const newvalue = weight * oldvalue + (1. - weight) * phivalue;

// Set phi at the boundary to this value
for (int k = 0; k < mesh->LocalNz; k++) {
for (int k = mesh->zstart; k <= mesh->zend; k++) {
phi(mesh->xstart - 1, j, k) = 2. * newvalue - phi(mesh->xstart, j, k);
phi(mesh->xstart - 2, j, k) = phi(mesh->xstart - 1, j, k);
}
Expand All @@ -1403,7 +1403,7 @@ class ELMpb : public PhysicsModel {
if (mesh->lastX()) {
for (int j = mesh->ystart; j <= mesh->yend; j++) {
BoutReal phivalue = 0.0;
for (int k = 0; k < mesh->LocalNz; k++) {
for (int k = mesh->zstart; k <= mesh->zend; k++) {
phivalue += phi(mesh->xend, j, k);
}
phivalue /= mesh->LocalNz; // Average in Z of point next to boundary
Expand All @@ -1416,7 +1416,7 @@ class ELMpb : public PhysicsModel {
BoutReal const newvalue = weight * oldvalue + (1. - weight) * phivalue;

// Set phi at the boundary to this value
for (int k = 0; k < mesh->LocalNz; k++) {
for (int k = mesh->zstart; k <= mesh->zend; k++) {
phi(mesh->xend + 1, j, k) = 2. * newvalue - phi(mesh->xend, j, k);
phi(mesh->xend + 2, j, k) = phi(mesh->xend + 1, j, k);
}
Expand All @@ -1441,7 +1441,7 @@ class ELMpb : public PhysicsModel {

if (mesh->firstX()) {
for (int j = mesh->ystart; j <= mesh->yend; j++) {
for (int k = 0; k < mesh->LocalNz; k++) {
for (int k = mesh->zstart; k <= mesh->zend; k++) {
// Average phi + Pi at the boundary, and set the boundary cell
// to this value. The phi solver will then put the value back
// onto the cell mid-point
Expand All @@ -1453,7 +1453,7 @@ class ELMpb : public PhysicsModel {

if (mesh->lastX()) {
for (int j = mesh->ystart; j <= mesh->yend; j++) {
for (int k = 0; k < mesh->LocalNz; k++) {
for (int k = mesh->zstart; k <= mesh->zend; k++) {
phi_shift(mesh->xend + 1, j, k) =
0.5 * (phi_shift(mesh->xend + 1, j, k) + phi_shift(mesh->xend, j, k));
}
Expand Down Expand Up @@ -1513,7 +1513,7 @@ class ELMpb : public PhysicsModel {
if (mesh->firstX()) {
for (int i = mesh->xstart - 2; i >= 0; --i) {
for (int j = mesh->ystart; j <= mesh->yend; ++j) {
for (int k = 0; k < mesh->LocalNz; ++k) {
for (int k = mesh->zstart; k <= mesh->zend; ++k) {
phi(i, j, k) = phi(i + 1, j, k);
}
}
Expand All @@ -1523,7 +1523,7 @@ class ELMpb : public PhysicsModel {
if (mesh->lastX()) {
for (int i = mesh->xend + 2; i < mesh->LocalNx; ++i) {
for (int j = mesh->ystart; j <= mesh->yend; ++j) {
for (int k = 0; k < mesh->LocalNz; ++k) {
for (int k = mesh->zstart; k <= mesh->zend; ++k) {
phi(i, j, k) = phi(i - 1, j, k);
}
}
Expand Down
4 changes: 2 additions & 2 deletions examples/performance/iterator-offsets/iterator-offsets.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ int main(int argc, char** argv) {
ITERATOR_TEST_BLOCK(
"Nested loop", for (int i = 0; i < mesh->LocalNx; ++i) {
for (int j = mesh->ystart; j < mesh->yend; ++j) {
for (int k = 0; k < mesh->LocalNz; ++k) {
for (int k = mesh->zstart; k <= mesh->zend; ++k) {
result(i, j, k) = (a(i, j + 1, k) - a(i, j - 1, k));
}
}
Expand All @@ -76,7 +76,7 @@ int main(int argc, char** argv) {
BOUT_OMP_PERF(parallel for)
for(int i=0;i<mesh->LocalNx;++i) {
for (int j = mesh->ystart; j < mesh->yend; ++j) {
for (int k = 0; k < mesh->LocalNz; ++k) {
for (int k = mesh->zstart; k <= mesh->zend; ++k) {
result(i, j, k) = (a(i, j + 1, k) - a(i, j - 1, k));
}
}
Expand Down
4 changes: 2 additions & 2 deletions examples/performance/iterator/iterator.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ int main(int argc, char** argv) {
ITERATOR_TEST_BLOCK(
"Nested loop", for (int i = 0; i < mesh->LocalNx; ++i) {
for (int j = 0; j < mesh->LocalNy; ++j) {
for (int k = 0; k < mesh->LocalNz; ++k) {
for (int k = mesh->zstart; k <= mesh->zend; ++k) {
result(i, j, k) = a(i, j, k) + b(i, j, k);
}
}
Expand All @@ -88,7 +88,7 @@ int main(int argc, char** argv) {
BOUT_OMP_PERF(parallel for)
for(int i=0;i<mesh->LocalNx;++i) {
for (int j = 0; j < mesh->LocalNy; ++j) {
for (int k = 0; k < mesh->LocalNz; ++k) {
for (int k = mesh->zstart; k <= mesh->zend; ++k) {
result(i, j, k) = a(i, j, k) + b(i, j, k);
}
}
Expand Down
2 changes: 1 addition & 1 deletion include/bout/fv_ops.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -258,7 +258,7 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in,
BoutReal flux_factor_lm =
common_factor / (coord->dy(i, j - 1) * coord->J(i, j - 1));
#endif
for (int k = 0; k < mesh->LocalNz; k++) {
for (int k = mesh->zstart; k <= mesh->zend; k++) {
#if BOUT_USE_METRIC_3D
// For right cell boundaries
BoutReal common_factor =
Expand Down
6 changes: 3 additions & 3 deletions manual/sphinx/developer_docs/data_types.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -255,9 +255,9 @@ parallelise and vectorise. Some tuning of this is possible, see below
for details. It replaces the C-style triple-nested loop::

Field3D f(0.0);
for (int i = mesh->xstart; i < mesh->xend; ++i) {
for (int j = mesh->ystart; j < mesh->yend; ++j) {
for (int k = 0; k < mesh->LocalNz; ++k) {
for (int i = mesh->xstart; i <= mesh->xend; ++i) {
for (int j = mesh->ystart; j <= mesh->yend; ++j) {
for (int k = mesh->zstart; k <= mesh->zend; ++k) {
f(i,j,k) = a(i,j,k) + b(i,j,k)
}
}
Expand Down
2 changes: 1 addition & 1 deletion manual/sphinx/user_docs/boundary_options.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -575,7 +575,7 @@ and implemented in ``boundary_standard.cxx``

void BoundaryNeumann::apply(Field3D &f) {
for(bndry->first(); !bndry->isDone(); bndry->next())
for(int z=0;z<mesh->LocalNz;z++)
for(int z= mesh->zstart; z <= mesh->zend;z++)
f[bndry->x][bndry->y][z] = f[bndry->x - bndry->bx][bndry->y -
bndry->by][z];
}
Expand Down
4 changes: 2 additions & 2 deletions src/invert/laplace/impls/naulin/naulin_laplace.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -402,7 +402,7 @@ void LaplaceNaulin::copy_x_boundaries(Field3D& x, const Field3D& x0, Mesh* local
if (localmesh->firstX()) {
for (int i = localmesh->xstart - 1; i >= 0; --i) {
for (int j = localmesh->ystart; j <= localmesh->yend; ++j) {
for (int k = 0; k < localmesh->LocalNz; ++k) {
for (int k = localmesh->zstart; k <= localmesh->zend; ++k) {
x(i, j, k) = x0(i, j, k);
}
}
Expand All @@ -411,7 +411,7 @@ void LaplaceNaulin::copy_x_boundaries(Field3D& x, const Field3D& x0, Mesh* local
if (localmesh->lastX()) {
for (int i = localmesh->xend + 1; i < localmesh->LocalNx; ++i) {
for (int j = localmesh->ystart; j <= localmesh->yend; ++j) {
for (int k = 0; k < localmesh->LocalNz; ++k) {
for (int k = localmesh->zstart; k <= localmesh->zend; ++k) {
x(i, j, k) = x0(i, j, k);
}
}
Expand Down
40 changes: 20 additions & 20 deletions src/invert/laplace/impls/petsc/petsc_laplace.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -145,7 +145,7 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in,
if (fourth_order) {
// first and last 2*localmesh-LocalNz entries are the edge x-values that (may) have 'off-diagonal' components (i.e. on another processor)
if (localmesh->firstX() && localmesh->lastX()) {
for (int i = 0; i < localmesh->LocalNz; i++) {
for (int i = localmesh->zstart; i <= localmesh->zend; i++) {
d_nnz[i] = 15;
d_nnz[localN - 1 - i] = 15;
o_nnz[i] = 0;
Expand All @@ -158,7 +158,7 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in,
o_nnz[localN - 1 - i] = 0;
}
} else if (localmesh->firstX()) {
for (int i = 0; i < localmesh->LocalNz; i++) {
for (int i = localmesh->zstart; i <= localmesh->zend; i++) {
d_nnz[i] = 15;
d_nnz[localN - 1 - i] = 15;
o_nnz[i] = 0;
Expand All @@ -171,7 +171,7 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in,
o_nnz[localN - 1 - i] = 5;
}
} else if (localmesh->lastX()) {
for (int i = 0; i < localmesh->LocalNz; i++) {
for (int i = localmesh->zstart; i <= localmesh->zend; i++) {
d_nnz[i] = 15;
d_nnz[localN - 1 - i] = 15;
o_nnz[i] = 10;
Expand All @@ -184,7 +184,7 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in,
o_nnz[localN - 1 - i] = 0;
}
} else {
for (int i = 0; i < localmesh->LocalNz; i++) {
for (int i = localmesh->zstart; i <= localmesh->zend; i++) {
d_nnz[i] = 15;
d_nnz[localN - 1 - i] = 15;
o_nnz[i] = 10;
Expand Down Expand Up @@ -215,28 +215,28 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in,
} else {
// first and last localmesh->LocalNz entries are the edge x-values that (may) have 'off-diagonal' components (i.e. on another processor)
if (localmesh->firstX() && localmesh->lastX()) {
for (int i = 0; i < localmesh->LocalNz; i++) {
for (int i = localmesh->zstart; i <= localmesh->zend; i++) {
d_nnz[i] = 6;
d_nnz[localN - 1 - i] = 6;
o_nnz[i] = 0;
o_nnz[localN - 1 - i] = 0;
}
} else if (localmesh->firstX()) {
for (int i = 0; i < localmesh->LocalNz; i++) {
for (int i = localmesh->zstart; i <= localmesh->zend; i++) {
d_nnz[i] = 6;
d_nnz[localN - 1 - i] = 6;
o_nnz[i] = 0;
o_nnz[localN - 1 - i] = 3;
}
} else if (localmesh->lastX()) {
for (int i = 0; i < localmesh->LocalNz; i++) {
for (int i = localmesh->zstart; i <= localmesh->zend; i++) {
d_nnz[i] = 6;
d_nnz[localN - 1 - i] = 6;
o_nnz[i] = 3;
o_nnz[localN - 1 - i] = 0;
}
} else {
for (int i = 0; i < localmesh->LocalNz; i++) {
for (int i = localmesh->zstart; i <= localmesh->zend; i++) {
d_nnz[i] = 6;
d_nnz[localN - 1 - i] = 6;
o_nnz[i] = 3;
Expand Down Expand Up @@ -377,7 +377,7 @@ FieldPerp LaplacePetsc::solve(const FieldPerp& b, const FieldPerp& x0) {
// Set the values for the inner boundary region
if (localmesh->firstX()) {
for (int x = 0; x < localmesh->xstart; x++) {
for (int z = 0; z < localmesh->LocalNz; z++) {
for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
PetscScalar val; // Value of element to be set in the matrix
// If Neumann Boundary Conditions are set.
if (isInnerBoundaryFlagSet(INVERT_AC_GRAD)) {
Expand Down Expand Up @@ -461,7 +461,7 @@ FieldPerp LaplacePetsc::solve(const FieldPerp& b, const FieldPerp& x0) {

// Set the values for the main domain
for (int x = localmesh->xstart; x <= localmesh->xend; x++) {
for (int z = 0; z < localmesh->LocalNz; z++) {
for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
// NOTE: Only A0 is the A from setCoefA ()
BoutReal A0, A1, A2, A3, A4, A5;
A0 = A(x, y, z);
Expand Down Expand Up @@ -639,7 +639,7 @@ FieldPerp LaplacePetsc::solve(const FieldPerp& b, const FieldPerp& x0) {
// Set the values for the outer boundary region
if (localmesh->lastX()) {
for (int x = localmesh->xend + 1; x < localmesh->LocalNx; x++) {
for (int z = 0; z < localmesh->LocalNz; z++) {
for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
// Set Diagonal Values to 1
PetscScalar val = 1;
Element(i, x, z, 0, 0, val, MatA);
Expand Down Expand Up @@ -829,7 +829,7 @@ FieldPerp LaplacePetsc::solve(const FieldPerp& b, const FieldPerp& x0) {
// Set the inner boundary values
if (localmesh->firstX()) {
for (int x = 0; x < localmesh->xstart; x++) {
for (int z = 0; z < localmesh->LocalNz; z++) {
for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
PetscScalar val = 0;
VecGetValues(xs, 1, &i, &val);
sol[x][z] = val;
Expand All @@ -840,7 +840,7 @@ FieldPerp LaplacePetsc::solve(const FieldPerp& b, const FieldPerp& x0) {

// Set the main domain values
for (int x = localmesh->xstart; x <= localmesh->xend; x++) {
for (int z = 0; z < localmesh->LocalNz; z++) {
for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
PetscScalar val = 0;
VecGetValues(xs, 1, &i, &val);
sol[x][z] = val;
Expand All @@ -851,7 +851,7 @@ FieldPerp LaplacePetsc::solve(const FieldPerp& b, const FieldPerp& x0) {
// Set the outer boundary values
if (localmesh->lastX()) {
for (int x = localmesh->xend + 1; x < localmesh->LocalNx; x++) {
for (int z = 0; z < localmesh->LocalNz; z++) {
for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
PetscScalar val = 0;
VecGetValues(xs, 1, &i, &val);
sol[x][z] = val;
Expand Down Expand Up @@ -1076,7 +1076,7 @@ void LaplacePetsc::vecToField(Vec xs, FieldPerp& f) {
int i = Istart;
if (localmesh->firstX()) {
for (int x = 0; x < localmesh->xstart; x++) {
for (int z = 0; z < localmesh->LocalNz; z++) {
for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
PetscScalar val;
VecGetValues(xs, 1, &i, &val);
f[x][z] = val;
Expand All @@ -1086,7 +1086,7 @@ void LaplacePetsc::vecToField(Vec xs, FieldPerp& f) {
}

for (int x = localmesh->xstart; x <= localmesh->xend; x++) {
for (int z = 0; z < localmesh->LocalNz; z++) {
for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
PetscScalar val;
VecGetValues(xs, 1, &i, &val);
f[x][z] = val;
Expand All @@ -1096,7 +1096,7 @@ void LaplacePetsc::vecToField(Vec xs, FieldPerp& f) {

if (localmesh->lastX()) {
for (int x = localmesh->xend + 1; x < localmesh->LocalNx; x++) {
for (int z = 0; z < localmesh->LocalNz; z++) {
for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
PetscScalar val;
VecGetValues(xs, 1, &i, &val);
f[x][z] = val;
Expand All @@ -1113,7 +1113,7 @@ void LaplacePetsc::fieldToVec(const FieldPerp& f, Vec bs) {
int i = Istart;
if (localmesh->firstX()) {
for (int x = 0; x < localmesh->xstart; x++) {
for (int z = 0; z < localmesh->LocalNz; z++) {
for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
PetscScalar val = f[x][z];
VecSetValues(bs, 1, &i, &val, INSERT_VALUES);
i++; // Increment row in Petsc matrix
Expand All @@ -1122,7 +1122,7 @@ void LaplacePetsc::fieldToVec(const FieldPerp& f, Vec bs) {
}

for (int x = localmesh->xstart; x <= localmesh->xend; x++) {
for (int z = 0; z < localmesh->LocalNz; z++) {
for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
PetscScalar val = f[x][z];
VecSetValues(bs, 1, &i, &val, INSERT_VALUES);
i++; // Increment row in Petsc matrix
Expand All @@ -1131,7 +1131,7 @@ void LaplacePetsc::fieldToVec(const FieldPerp& f, Vec bs) {

if (localmesh->lastX()) {
for (int x = localmesh->xend + 1; x < localmesh->LocalNx; x++) {
for (int z = 0; z < localmesh->LocalNz; z++) {
for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
PetscScalar val = f[x][z];
VecSetValues(bs, 1, &i, &val, INSERT_VALUES);
i++; // Increment row in Petsc matrix
Expand Down
12 changes: 6 additions & 6 deletions src/invert/laplace/impls/spt/spt.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,15 +95,15 @@ FieldPerp LaplaceSPT::solve(const FieldPerp& b, const FieldPerp& x0) {
if (isInnerBoundaryFlagSetOnFirstX(INVERT_SET)) {
// Copy x0 inner boundary into bs
for (int ix = 0; ix < xbndry; ix++) {
for (int iz = 0; iz < localmesh->LocalNz; iz++) {
for (int iz = localmesh->zstart; iz <= localmesh->zend; iz++) {
bs[ix][iz] = x0[ix][iz];
}
}
}
if (isOuterBoundaryFlagSetOnLastX(INVERT_SET)) {
// Copy x0 outer boundary into bs
for (int ix = localmesh->LocalNx - 1; ix >= localmesh->LocalNx - xbndry; ix--) {
for (int iz = 0; iz < localmesh->LocalNz; iz++) {
for (int iz = localmesh->zstart; iz <= localmesh->zend; iz++) {
bs[ix][iz] = x0[ix][iz];
}
}
Expand Down Expand Up @@ -181,7 +181,7 @@ Field3D LaplaceSPT::solve(const Field3D& b, const Field3D& x0) {
// Copy x0 inner boundary into bs
for (int ix = 0; ix < xbndry; ix++) {
for (int iy = 0; iy < localmesh->LocalNy; iy++) {
for (int iz = 0; iz < localmesh->LocalNz; iz++) {
for (int iz = localmesh->zstart; iz <= localmesh->zend; iz++) {
bs(ix, iy, iz) = x0(ix, iy, iz);
}
}
Expand All @@ -191,7 +191,7 @@ Field3D LaplaceSPT::solve(const Field3D& b, const Field3D& x0) {
// Copy x0 outer boundary into bs
for (int ix = localmesh->LocalNx - 1; ix >= localmesh->LocalNx - xbndry; ix--) {
for (int iy = 0; iy < localmesh->LocalNy; iy++) {
for (int iz = 0; iz < localmesh->LocalNz; iz++) {
for (int iz = localmesh->zstart; iz <= localmesh->zend; iz++) {
bs(ix, iy, iz) = x0(ix, iy, iz);
}
}
Expand Down Expand Up @@ -519,15 +519,15 @@ void LaplaceSPT::finish(SPT_data& data, FieldPerp& x) {
if (!localmesh->firstX()) {
// Set left boundary to zero (Prevent unassigned values in corners)
for (int ix = 0; ix < localmesh->xstart; ix++) {
for (int kz = 0; kz < localmesh->LocalNz; kz++) {
for (int kz = localmesh->zstart; kz <= localmesh->zend; kz++) {
x(ix, kz) = 0.0;
}
}
}
if (!localmesh->lastX()) {
// Same for right boundary
for (int ix = localmesh->xend + 1; ix < localmesh->LocalNx; ix++) {
for (int kz = 0; kz < localmesh->LocalNz; kz++) {
for (int kz = localmesh->zstart; kz <= localmesh->zend; kz++) {
x(ix, kz) = 0.0;
}
}
Expand Down
Loading
Loading