From 0386f676712594aea5c4df308804b39139a5a884 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Wed, 7 Jan 2026 17:09:12 +0000 Subject: [PATCH] Use `zstart/zend` in more loops --- examples/elm-pb/elm_pb.cxx | 18 +-- .../iterator-offsets/iterator-offsets.cxx | 4 +- examples/performance/iterator/iterator.cxx | 4 +- include/bout/fv_ops.hxx | 2 +- manual/sphinx/developer_docs/data_types.rst | 6 +- manual/sphinx/user_docs/boundary_options.rst | 2 +- .../laplace/impls/naulin/naulin_laplace.cxx | 4 +- .../laplace/impls/petsc/petsc_laplace.cxx | 40 +++--- src/invert/laplace/impls/spt/spt.cxx | 12 +- .../laplacexz/impls/petsc/laplacexz-petsc.cxx | 50 +++---- src/mesh/boundary_standard.cxx | 136 +++++++++--------- src/mesh/coordinates.cxx | 6 +- src/mesh/difops.cxx | 4 +- src/mesh/fv_ops.cxx | 16 +-- src/solver/impls/imex-bdf2/imex-bdf2.cxx | 20 +-- src/solver/impls/petsc/petsc.cxx | 6 +- src/solver/impls/snes/snes.cxx | 10 +- tests/integrated/test-invpar/test_invpar.cxx | 2 +- .../test_multigrid_laplace.cxx | 10 +- .../test_naulin_laplace.cxx | 10 +- .../test-petsc_laplace/test_petsc_laplace.cxx | 54 +++---- .../test_petsc_laplace_MAST_grid.cxx | 66 ++++----- 22 files changed, 241 insertions(+), 241 deletions(-) diff --git a/examples/elm-pb/elm_pb.cxx b/examples/elm-pb/elm_pb.cxx index fe18a756ae..57f0ec0778 100644 --- a/examples/elm-pb/elm_pb.cxx +++ b/examples/elm-pb/elm_pb.cxx @@ -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); } } @@ -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 @@ -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); } @@ -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 @@ -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); } @@ -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 @@ -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)); } @@ -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); } } @@ -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); } } diff --git a/examples/performance/iterator-offsets/iterator-offsets.cxx b/examples/performance/iterator-offsets/iterator-offsets.cxx index 455df38ea5..48012f1793 100644 --- a/examples/performance/iterator-offsets/iterator-offsets.cxx +++ b/examples/performance/iterator-offsets/iterator-offsets.cxx @@ -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)); } } @@ -76,7 +76,7 @@ int main(int argc, char** argv) { BOUT_OMP_PERF(parallel for) for(int i=0;iLocalNx;++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)); } } diff --git a/examples/performance/iterator/iterator.cxx b/examples/performance/iterator/iterator.cxx index 7a29b00298..e1fc59d067 100644 --- a/examples/performance/iterator/iterator.cxx +++ b/examples/performance/iterator/iterator.cxx @@ -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); } } @@ -88,7 +88,7 @@ int main(int argc, char** argv) { BOUT_OMP_PERF(parallel for) for(int i=0;iLocalNx;++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); } } diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 94007a57a2..df61c6aa24 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -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 = diff --git a/manual/sphinx/developer_docs/data_types.rst b/manual/sphinx/developer_docs/data_types.rst index fa8e9e6ea6..f9411dcb39 100644 --- a/manual/sphinx/developer_docs/data_types.rst +++ b/manual/sphinx/developer_docs/data_types.rst @@ -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) } } diff --git a/manual/sphinx/user_docs/boundary_options.rst b/manual/sphinx/user_docs/boundary_options.rst index a3cdf0078b..8493049516 100644 --- a/manual/sphinx/user_docs/boundary_options.rst +++ b/manual/sphinx/user_docs/boundary_options.rst @@ -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;zLocalNz;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]; } diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.cxx b/src/invert/laplace/impls/naulin/naulin_laplace.cxx index e6f68d850d..f61dc93908 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.cxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.cxx @@ -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); } } @@ -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); } } diff --git a/src/invert/laplace/impls/petsc/petsc_laplace.cxx b/src/invert/laplace/impls/petsc/petsc_laplace.cxx index 5f417d4faa..4170c88e80 100644 --- a/src/invert/laplace/impls/petsc/petsc_laplace.cxx +++ b/src/invert/laplace/impls/petsc/petsc_laplace.cxx @@ -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; @@ -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; @@ -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; @@ -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; @@ -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; @@ -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)) { @@ -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); @@ -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); @@ -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; @@ -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; @@ -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; @@ -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; @@ -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; @@ -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; @@ -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 @@ -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 @@ -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 diff --git a/src/invert/laplace/impls/spt/spt.cxx b/src/invert/laplace/impls/spt/spt.cxx index 2e4c844c94..e39ca7e89f 100644 --- a/src/invert/laplace/impls/spt/spt.cxx +++ b/src/invert/laplace/impls/spt/spt.cxx @@ -95,7 +95,7 @@ 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]; } } @@ -103,7 +103,7 @@ FieldPerp LaplaceSPT::solve(const FieldPerp& b, const FieldPerp& x0) { 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]; } } @@ -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); } } @@ -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); } } @@ -519,7 +519,7 @@ 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; } } @@ -527,7 +527,7 @@ void LaplaceSPT::finish(SPT_data& data, FieldPerp& x) { 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; } } diff --git a/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx b/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx index 769e797352..c61db141b8 100644 --- a/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx +++ b/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx @@ -202,25 +202,25 @@ LaplaceXZpetsc::LaplaceXZpetsc(Mesh* m, Options* opt, const CELL_LOC loc) // X boundaries if (localmesh->firstX()) { - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { d_nnz[z] = 2; } } else { // One point on another processor - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { d_nnz[z] -= 1; o_nnz[z] += 1; } } if (localmesh->lastX()) { - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { int ind = localN - (localmesh->LocalNz) + z; d_nnz[ind] = 2; } } else { // One point on another processor - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { int ind = localN - (localmesh->LocalNz) + z; d_nnz[ind] -= 1; o_nnz[ind] += 1; @@ -335,7 +335,7 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) { /* NOTE: Sign of the elements are opposite of what one might expect, * see note about BC in LaplaceXZ constructor for more details */ - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { PetscScalar val = 1.0; MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES); @@ -347,7 +347,7 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) { } } else if (inner_boundary_flags & INVERT_SET) { // Setting BC from x0 - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { PetscScalar val = 1.0; MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES); @@ -359,7 +359,7 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) { } } else if (inner_boundary_flags & INVERT_RHS) { // Setting BC from b - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { PetscScalar val = 1.0; MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES); @@ -374,7 +374,7 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) { /* NOTE: Sign of the elements are opposite of what one might expect, * see note about BC in LaplaceXZ constructor for more details */ - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { PetscScalar val = 0.5; MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES); @@ -393,7 +393,7 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) { Coordinates* coords = localmesh->getCoordinates(location); 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++) { // stencil entries PetscScalar c, xm, xp, zm, zp; // Diagonal entries @@ -655,7 +655,7 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) { if (localmesh->lastX()) { if (outer_boundary_flags & INVERT_AC_GRAD) { // Neumann 0 - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { PetscScalar val = 1.0; MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES); @@ -667,7 +667,7 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) { } } else if (outer_boundary_flags & INVERT_SET) { // Setting BC from x0 - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { PetscScalar val = 1.0; MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES); @@ -679,7 +679,7 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) { } } else if (outer_boundary_flags & INVERT_RHS) { // Setting BC from b - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { PetscScalar val = 1.0; MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES); @@ -693,7 +693,7 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) { //Default: Dirichlet on outer X boundary PetscScalar val = 0.5; - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES); int col = row - (localmesh->LocalNz); // -1 in X @@ -804,7 +804,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) { if (localmesh->firstX()) { if (inner_boundary_flags & INVERT_AC_GRAD) { // Neumann 0 - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { // Setting the initial guess x0 PetscScalar val = x0(localmesh->xstart - 1, y, z); VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); @@ -816,7 +816,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) { } } else if (inner_boundary_flags & INVERT_SET) { // Setting BC from x0 - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { // Setting the initial guess x0 PetscScalar val = x0(localmesh->xstart - 1, y, z); VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); @@ -828,7 +828,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) { } } else if (inner_boundary_flags & INVERT_RHS) { // Setting BC from b - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { // Setting the initial guess x0 PetscScalar val = x0(localmesh->xstart - 1, y, z); VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); @@ -840,7 +840,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) { } } else { // Default: Neumann on inner x boundary - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { // Setting the initial guess x0 PetscScalar val = x0(localmesh->xstart - 1, y, z); VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); @@ -855,7 +855,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) { // Set the inner points 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 = x0(x, y, z); VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); @@ -869,7 +869,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) { if (localmesh->lastX()) { if (outer_boundary_flags & INVERT_AC_GRAD) { // Neumann 0 - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { // Setting the initial guess x0 PetscScalar val = x0(localmesh->xend + 1, y, z); VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); @@ -882,7 +882,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) { } } else if (outer_boundary_flags & INVERT_SET) { // Setting BC from x0 - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { // Setting the initial guess x0 PetscScalar val = x0(localmesh->xend + 1, y, z); VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); @@ -895,7 +895,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) { } } else if (outer_boundary_flags & INVERT_RHS) { // Setting BC from b - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { // Setting the initial guess x0 PetscScalar val = x0(localmesh->xend + 1, y, z); VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); @@ -908,7 +908,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) { } } else { //Default: Dirichlet on outer X boundary - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { // Setting the initial guess x0 PetscScalar val = x0(localmesh->xend + 1, y, z); VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); @@ -952,7 +952,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) { ind = Istart; // Inner X boundary if (localmesh->firstX()) { - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { PetscScalar val; VecGetValues(xs, 1, &ind, &val); for (int x = localmesh->xstart - 1; x >= 0; --x) { @@ -963,7 +963,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) { } 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, &ind, &val); result(x, y, z) = val; @@ -973,7 +973,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) { // Outer X boundary if (localmesh->lastX()) { - for (int z = 0; z < localmesh->LocalNz; z++) { + for (int z = localmesh->zstart; z <= localmesh->zend; z++) { PetscScalar val; VecGetValues(xs, 1, &ind, &val); for (int x = localmesh->xend + 1; x < localmesh->LocalNx; ++x) { diff --git a/src/mesh/boundary_standard.cxx b/src/mesh/boundary_standard.cxx index 8f2c7df5ec..c3fd5f8f2e 100644 --- a/src/mesh/boundary_standard.cxx +++ b/src/mesh/boundary_standard.cxx @@ -325,7 +325,7 @@ void BoundaryDirichlet::apply(Field3D& f, BoutReal t) { // Outer x boundary for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -346,7 +346,7 @@ void BoundaryDirichlet::apply(Field3D& f, BoutReal t) { if (bndry->bx < 0) { // Inner x boundary. Set one point inwards for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -368,7 +368,7 @@ void BoundaryDirichlet::apply(Field3D& f, BoutReal t) { if (bndry->by != 0) { // y boundaries for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -394,7 +394,7 @@ void BoundaryDirichlet::apply(Field3D& f, BoutReal t) { // Upper y boundary boundary for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -415,7 +415,7 @@ void BoundaryDirichlet::apply(Field3D& f, BoutReal t) { if (bndry->by < 0) { // Lower y boundary. Set one point inwards for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -436,7 +436,7 @@ void BoundaryDirichlet::apply(Field3D& f, BoutReal t) { if (bndry->bx != 0) { // x boundaries for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -469,7 +469,7 @@ void BoundaryDirichlet::apply(Field3D& f, BoutReal t) { * (mesh->GlobalY(bndry->y) // In the guard cell + mesh->GlobalY(bndry->y - bndry->by)); // the grid cell - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate( Context(bndry, zk, loc, t, mesh).set("x", xnorm, "y", ynorm)); @@ -513,7 +513,7 @@ void BoundaryDirichlet::apply(Field3D& f, BoutReal t) { const int yi = bndry->y + (i * bndry->by); const auto xnorm = mesh->GlobalX(xi); const auto ynorm = mesh->GlobalY(yi) * TWOPI; - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate( Context(bndry, zk, loc, t, mesh).set("x", xnorm, "y", ynorm)); @@ -528,7 +528,7 @@ void BoundaryDirichlet::apply(Field3D& f, BoutReal t) { } else { // Standard (non-staggered) case for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -572,7 +572,7 @@ void BoundaryDirichlet::apply(Field3D& f, BoutReal t) { // Set any other guard cells using the values on the cells int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -596,7 +596,7 @@ void BoundaryDirichlet::apply_ddt(Field3D& f) { Field3D* dt = f.timeDeriv(); for (bndry->first(); !bndry->isDone(); bndry->next()) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { (*dt)(bndry->x, bndry->y, z) = 0.; // Set time derivative to zero } } @@ -827,7 +827,7 @@ void BoundaryDirichlet_O3::apply(Field3D& f, BoutReal t) { // Outer x boundary for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -848,7 +848,7 @@ void BoundaryDirichlet_O3::apply(Field3D& f, BoutReal t) { if (bndry->bx < 0) { // Inner x boundary. Set one point inwards for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -870,7 +870,7 @@ void BoundaryDirichlet_O3::apply(Field3D& f, BoutReal t) { // y boundaries for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -898,7 +898,7 @@ void BoundaryDirichlet_O3::apply(Field3D& f, BoutReal t) { // Upper y boundary for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -919,7 +919,7 @@ void BoundaryDirichlet_O3::apply(Field3D& f, BoutReal t) { if (bndry->by < 0) { // Lower y boundary. Set one point inwards for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -940,7 +940,7 @@ void BoundaryDirichlet_O3::apply(Field3D& f, BoutReal t) { if (bndry->bx != 0) { // x boundaries for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -974,7 +974,7 @@ void BoundaryDirichlet_O3::apply(Field3D& f, BoutReal t) { * (mesh->GlobalY(bndry->y) // In the guard cell + mesh->GlobalY(bndry->y - bndry->by)); // the grid cell - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate( Context(bndry, zk, loc, t, mesh).set("x", xnorm, "y", ynorm)); @@ -1000,7 +1000,7 @@ void BoundaryDirichlet_O3::apply(Field3D& f, BoutReal t) { } else { // Standard (non-staggered) case for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -1038,7 +1038,7 @@ void BoundaryDirichlet_O3::apply_ddt(Field3D& f) { bndry->first(); for (bndry->first(); !bndry->isDone(); bndry->next()) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { (*dt)(bndry->x, bndry->y, z) = 0.; // Set time derivative to zero } } @@ -1284,7 +1284,7 @@ void BoundaryDirichlet_O4::apply(Field3D& f, BoutReal t) { // Outer x boundary for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -1306,7 +1306,7 @@ void BoundaryDirichlet_O4::apply(Field3D& f, BoutReal t) { if (bndry->bx < 0) { // Inner x boundary. Set one point inwards for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -1329,7 +1329,7 @@ void BoundaryDirichlet_O4::apply(Field3D& f, BoutReal t) { // y boundaries for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -1358,7 +1358,7 @@ void BoundaryDirichlet_O4::apply(Field3D& f, BoutReal t) { // Outer y boundary for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -1381,7 +1381,7 @@ void BoundaryDirichlet_O4::apply(Field3D& f, BoutReal t) { if (bndry->by < 0) { // Inner y boundary. Set one point inwards for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -1405,7 +1405,7 @@ void BoundaryDirichlet_O4::apply(Field3D& f, BoutReal t) { // x boundaries for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -1440,7 +1440,7 @@ void BoundaryDirichlet_O4::apply(Field3D& f, BoutReal t) { * (mesh->GlobalY(bndry->y) // In the guard cell + mesh->GlobalY(bndry->y - bndry->by)); // the grid cell - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate( Context(bndry, zk, loc, t, mesh).set("x", xnorm, "y", ynorm)); @@ -1468,7 +1468,7 @@ void BoundaryDirichlet_O4::apply(Field3D& f, BoutReal t) { } else { // Standard (non-staggered) case for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); } @@ -1505,7 +1505,7 @@ void BoundaryDirichlet_O4::apply_ddt(Field3D& f) { ASSERT1(mesh == f.getMesh()); Field3D* dt = f.timeDeriv(); for (bndry->first(); !bndry->isDone(); bndry->next()) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { (*dt)(bndry->x, bndry->y, z) = 0.; // Set time derivative to zero } } @@ -1547,7 +1547,7 @@ void BoundaryDirichlet_4thOrder::apply(Field3D& f) { // Set (at 4th order) the value at the mid-point between the guard cell and the grid // cell to be val for (bndry->first(); !bndry->isDone(); bndry->next1d()) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { f(bndry->x, bndry->y, z) = 128. / 35. * val - 4. * f(bndry->x - bndry->bx, bndry->y - bndry->by, z) + 2. * f(bndry->x - 2 * bndry->bx, bndry->y - 2 * bndry->by, z) @@ -1574,7 +1574,7 @@ void BoundaryDirichlet_4thOrder::apply_ddt(Field3D& f) { ASSERT1(mesh == f.getMesh()); Field3D* dt = f.timeDeriv(); for (bndry->first(); !bndry->isDone(); bndry->next()) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { (*dt)(bndry->x, bndry->y, z) = 0.; // Set time derivative to zero } } @@ -1661,7 +1661,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { // Loop over all elements and set equal to the next point in for (bndry->first(); !bndry->isDone(); bndry->next1d()) { #if BOUT_USE_METRIC_3D - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { #else int z = 0; #endif @@ -1680,7 +1680,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { // because derivative values don't exist in boundary region // NOTE: should be fixed to interpolate to boundary line #if not(BOUT_USE_METRIC_3D) - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { #endif BoutReal xshift = g12shift * dfdy(bndry->x - bndry->bx, bndry->y, z) + g13shift * dfdz(bndry->x - bndry->bx, bndry->y, z); @@ -1968,7 +1968,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { if (bndry->bx > 0) { // Outer x boundary for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)) * metric->dx(bndry->x, bndry->y, zk); @@ -1997,7 +1997,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { if (bndry->bx < 0) { // Inner x boundary for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)) * metric->dx(bndry->x - bndry->bx, bndry->y, zk); @@ -2029,7 +2029,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { BoutReal delta = bndry->bx * metric->dx(bndry->x, bndry->y) + bndry->by * metric->dy(bndry->x, bndry->y); #endif - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { #if BOUT_USE_METRIC_3D BoutReal delta = bndry->bx * metric->dx(bndry->x, bndry->y, zk) + bndry->by * metric->dy(bndry->x, bndry->y, zk); @@ -2053,7 +2053,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { if (bndry->by > 0) { // Outer y boundary for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)) * metric->dy(bndry->x, bndry->y, zk); @@ -2081,7 +2081,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { if (bndry->by < 0) { // Inner y boundary. Set one point inwards for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)) * metric->dy(bndry->x, bndry->y - bndry->by, zk); @@ -2115,7 +2115,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { BoutReal delta = bndry->bx * metric->dx(bndry->x, bndry->y, zk) + bndry->by * metric->dy(bndry->x, bndry->y, zk); #endif - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { #if BOUT_USE_METRIC_3D BoutReal delta = bndry->bx * metric->dx(bndry->x, bndry->y, zk) + bndry->by * metric->dy(bndry->x, bndry->y, zk); @@ -2148,7 +2148,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { * (mesh->GlobalY(bndry->y) // In the guard cell + mesh->GlobalY(bndry->y - bndry->by)); // the grid cell - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { BoutReal delta = bndry->bx * metric->dx(bndry->x, bndry->y, zk) + bndry->by * metric->dy(bndry->x, bndry->y, zk); if (fg) { @@ -2170,13 +2170,13 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { } else { for (; !bndry->isDone(); bndry->next1d()) { #if BOUT_USE_METRIC_3D - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { BoutReal delta = bndry->bx * metric->dx(bndry->x, bndry->y, zk) + bndry->by * metric->dy(bndry->x, bndry->y, zk); #else BoutReal delta = bndry->bx * metric->dx(bndry->x, bndry->y) + bndry->by * metric->dy(bndry->x, bndry->y); - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { #endif if (fg) { val = fg->generate(Context(bndry, zk, loc, t, mesh)); @@ -2205,7 +2205,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { ASSERT1(mesh == f.getMesh()); Field3D* dt = f.timeDeriv(); for (bndry->first(); !bndry->isDone(); bndry->next()) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { (*dt)(bndry->x, bndry->y, z) = 0.; // Set time derivative to zero } } @@ -2303,7 +2303,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { } else { Coordinates* coords = f.getCoordinates(); for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { BoutReal delta = bndry->bx * coords->dx(bndry->x, bndry->y, zk) + bndry->by * coords->dy(bndry->x, bndry->y, zk); if (fg) { @@ -2339,7 +2339,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { ASSERT1(mesh == f.getMesh()); Field3D* dt = f.timeDeriv(); for (bndry->first(); !bndry->isDone(); bndry->next()) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { (*dt)(bndry->x, bndry->y, z) = 0.; // Set time derivative to zero } } @@ -2398,7 +2398,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { // grid cell to be val This sets the value of the co-ordinate derivative, i.e. DDX/DDY // not Grad_par/Grad_perp.x for (bndry->first(); !bndry->isDone(); bndry->next1d()) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { BoutReal delta = -(bndry->bx * metric->dx(bndry->x, bndry->y, z) + bndry->by * metric->dy(bndry->x, bndry->y, z)); f(bndry->x, bndry->y, z) = @@ -2431,7 +2431,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { ASSERT1(mesh == f.getMesh()); Field3D* dt = f.timeDeriv(); for (bndry->first(); !bndry->isDone(); bndry->next()) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { (*dt)(bndry->x, bndry->y, z) = 0.; // Set time derivative to zero } } @@ -2469,7 +2469,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { ASSERT1(mesh == f.getMesh()); Coordinates* metric = f.getCoordinates(); 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) * sqrt(metric->g_22(bndry->x, bndry->y, z) @@ -2536,7 +2536,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { ASSERT1(mesh == f.getMesh()); if (fabs(bval) < 1.e-12) { 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) = gval / aval; } } @@ -2546,7 +2546,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { sign = -1.; } 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) + sign * (gval - aval * f(bndry->x - bndry->bx, bndry->y - bndry->by, z)) @@ -2570,7 +2570,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { Mesh* mesh = bndry->localmesh; ASSERT1(mesh == f.getMesh()); 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) = 2. * f(bndry->x - bndry->bx, bndry->y - bndry->by, z) - f(bndry->x - 2 * bndry->bx, bndry->y - 2 * bndry->by, z); @@ -3191,7 +3191,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { for (int i = 0; i < bndry->width; i++) { int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; @@ -3205,7 +3205,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { // Inner x boundary. Set one point inwards for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { for (int i = -1; i < bndry->width; i++) { int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; @@ -3220,7 +3220,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { for (int i = 0; i < bndry->width; i++) { int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; @@ -3236,7 +3236,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { if (bndry->by > 0) { // Upper y boundary for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { for (int i = 0; i < bndry->width; i++) { int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; @@ -3250,7 +3250,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { // Lower y boundary. Set one point inwards for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { for (int i = -1; i < bndry->width; i++) { int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; @@ -3264,7 +3264,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { // x boundaries for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { for (int i = 0; i < bndry->width; i++) { int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; @@ -3279,7 +3279,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { // Standard (non-staggered) case for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { for (int i = 0; i < bndry->width; i++) { int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; @@ -3303,7 +3303,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { ASSERT1(mesh == f.getMesh()); Field3D* dt = f.timeDeriv(); for (bndry->first(); !bndry->isDone(); bndry->next()) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { (*dt)(bndry->x, bndry->y, z) = 0.; // Set time derivative to zero } } @@ -3453,7 +3453,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { for (int i = 0; i < bndry->width; i++) { int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; @@ -3468,7 +3468,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { // Inner x boundary. Set one point inwards for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { for (int i = -1; i < bndry->width; i++) { int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; @@ -3484,7 +3484,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { for (int i = 0; i < bndry->width; i++) { int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; @@ -3501,7 +3501,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { if (bndry->by > 0) { // Upper y boundary for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { for (int i = 0; i < bndry->width; i++) { int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; @@ -3516,7 +3516,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { // Lower y boundary. Set one point inwards for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { for (int i = -1; i < bndry->width; i++) { int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; @@ -3531,7 +3531,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { // x boundaries for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { for (int i = 0; i < bndry->width; i++) { int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; @@ -3547,7 +3547,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { // Standard (non-staggered) case for (; !bndry->isDone(); bndry->next1d()) { - for (int zk = 0; zk < mesh->LocalNz; zk++) { + for (int zk = mesh->zstart; zk <= mesh->zend; zk++) { for (int i = 0; i < bndry->width; i++) { int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; @@ -3572,7 +3572,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { ASSERT1(mesh == f.getMesh()); Field3D* dt = f.timeDeriv(); for (bndry->first(); !bndry->isDone(); bndry->next()) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { (*dt)(bndry->x, bndry->y, z) = 0.; // Set time derivative to zero } } @@ -3632,7 +3632,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { op->apply(g); // Set time-derivatives for (bndry->first(); !bndry->isDone(); bndry->next()) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { ddt(f)(bndry->x, bndry->y, z) = r * (g(bndry->x, bndry->y, z) - f(bndry->x, bndry->y, z)); } diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 3dfee6a553..986707a79e 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -225,7 +225,7 @@ Field3D interpolateAndExtrapolate(const Field3D& f_, CELL_LOC location, ASSERT1(bndry->bx == 0 or localmesh->xstart > 1); ASSERT1(bndry->by == 0 or localmesh->ystart > 1); // note that either bx or by is >0 here - for (int zi = 0; zi < localmesh->LocalNz; ++zi) { + for (int zi = localmesh->zstart; zi <= localmesh->zend; ++zi) { result(bndry->x, bndry->y, zi) = (9. * (f(bndry->x - bndry->bx, bndry->y - bndry->by, zi) @@ -256,7 +256,7 @@ Field3D interpolateAndExtrapolate(const Field3D& f_, CELL_LOC location, for (int i = extrap_start; i < bndry->width; i++) { int xi = bndry->x + i * bndry->bx; int yi = bndry->y + i * bndry->by; - for (int zi = 0; zi < localmesh->LocalNz; ++zi) { + for (int zi = localmesh->zstart; zi <= localmesh->zend; ++zi) { result(xi, yi, zi) = 3.0 * result(xi - bndry->bx, yi - bndry->by, zi) - 3.0 * result(xi - 2 * bndry->bx, yi - 2 * bndry->by, zi) @@ -266,7 +266,7 @@ Field3D interpolateAndExtrapolate(const Field3D& f_, CELL_LOC location, } else { // not enough grid points to extrapolate, set equal to last grid point for (int i = extrap_start; i < bndry->width; i++) { - for (int zi = 0; zi < localmesh->LocalNz; ++zi) { + for (int zi = localmesh->zstart; zi <= localmesh->zend; ++zi) { result(bndry->x + i * bndry->bx, bndry->y + i * bndry->by, zi) = result(bndry->x - bndry->bx, bndry->y - bndry->by, zi); } diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 42fa4d6ca5..64b072e5ec 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -779,7 +779,7 @@ Field3D bracket(const Field3D& f, const Field2D& g, BRACKET_METHOD method, for (int jy = mesh->ystart; jy <= mesh->yend; jy++) { const BoutReal partialFactor = 1.0 / (12 * metric->dz(jx, jy)); const BoutReal spacingFactor = partialFactor / metric->dx(jx, jy); - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { const int jzp = jz + 1 < ncz ? jz + 1 : 0; // Above is alternative to const int jzp = (jz + 1) % ncz; const int jzm = jz - 1 >= 0 ? jz - 1 : ncz - 1; @@ -904,7 +904,7 @@ Field3D bracket(const Field3D& f, const Field3D& g, BRACKET_METHOD method, int ncz = mesh->LocalNz; for (int y = mesh->ystart; y <= mesh->yend; y++) { for (int x = 1; x <= mesh->LocalNx - 2; x++) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { int zm = (z - 1 + ncz) % ncz; int zp = (z + 1) % ncz; diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx index fe5422b4d1..723e1d45da 100644 --- a/src/mesh/fv_ops.cxx +++ b/src/mesh/fv_ops.cxx @@ -48,7 +48,7 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) { for (int i = xs; i <= xe; 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++) { // Calculate flux from i to i+1 BoutReal fout = 0.5 * (a(i, j, k) + a(i + 1, j, k)) @@ -282,7 +282,7 @@ const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { mesh->yend; for (int j = ystart; j <= yend; j++) { - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = mesh->zstart; k <= mesh->zend; k++) { BoutReal dy3 = SQ(coord->dy(i, j, k)) * coord->dy(i, j, k); // 3rd derivative at upper boundary @@ -329,7 +329,7 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { if (j != mesh->yend || !has_upper_boundary) { - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = mesh->zstart; k <= mesh->zend; k++) { // Right boundary common factors const BoutReal common_factor = 0.25 * (coord->dy(i, j, k) + coord->dy(i, j + 1, k)) @@ -353,7 +353,7 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { // At a domain boundary // Use a one-sided difference formula - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = mesh->zstart; k <= mesh->zend; k++) { // Right boundary common factors const BoutReal common_factor = 0.25 * (coord->dy(i, j, k) + coord->dy(i, j + 1, k)) @@ -383,7 +383,7 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { // Calculate the fluxes if (j != mesh->ystart || !has_lower_boundary) { - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = mesh->zstart; k <= mesh->zend; k++) { const BoutReal common_factor = 0.25 * (coord->dy(i, j, k) + coord->dy(i, j + 1, k)) * (coord->J(i, j, k) + coord->J(i, j - 1, k)); @@ -402,7 +402,7 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { } } else { // On a domain (Y) boundary - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = mesh->zstart; k <= mesh->zend; k++) { const BoutReal common_factor = 0.25 * (coord->dy(i, j, k) + coord->dy(i, j + 1, k)) * (coord->J(i, j, k) + coord->J(i, j - 1, k)); @@ -463,7 +463,7 @@ void communicateFluxes(Field3D& f) { mesh->wait(xin); // Add to cells for (int y = mesh->ystart; y <= mesh->yend; y++) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { f(2, y, z) += f(0, y, z); } } @@ -472,7 +472,7 @@ void communicateFluxes(Field3D& f) { mesh->wait(xout); // Add to cells for (int y = mesh->ystart; y <= mesh->yend; y++) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { f(mesh->LocalNx - 3, y, z) += f(mesh->LocalNx - 1, y, z); } } diff --git a/src/solver/impls/imex-bdf2/imex-bdf2.cxx b/src/solver/impls/imex-bdf2/imex-bdf2.cxx index 07188110ab..5ad70a9b8b 100644 --- a/src/solver/impls/imex-bdf2/imex-bdf2.cxx +++ b/src/solver/impls/imex-bdf2/imex-bdf2.cxx @@ -331,7 +331,7 @@ void IMEXBDF2::constructSNES(SNES* snesIn) { if (mesh->firstX()) { // Lower X boundary for (int y = mesh->ystart; y <= mesh->yend; y++) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { int localIndex = ROUND(index(mesh->xstart, y, z)); ASSERT2((localIndex >= 0) && (localIndex < localN)); if (z == 0) { @@ -350,7 +350,7 @@ void IMEXBDF2::constructSNES(SNES* snesIn) { } else { // On another processor for (int y = mesh->ystart; y <= mesh->yend; y++) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { int localIndex = ROUND(index(mesh->xstart, y, z)); ASSERT2((localIndex >= 0) && (localIndex < localN)); if (z == 0) { @@ -373,7 +373,7 @@ void IMEXBDF2::constructSNES(SNES* snesIn) { if (mesh->lastX()) { // Upper X boundary for (int y = mesh->ystart; y <= mesh->yend; y++) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { int localIndex = ROUND(index(mesh->xend, y, z)); ASSERT2((localIndex >= 0) && (localIndex < localN)); if (z == 0) { @@ -392,7 +392,7 @@ void IMEXBDF2::constructSNES(SNES* snesIn) { } else { // On another processor for (int y = mesh->ystart; y <= mesh->yend; y++) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { int localIndex = ROUND(index(mesh->xend, y, z)); ASSERT2((localIndex >= 0) && (localIndex < localN)); if (z == 0) { @@ -559,7 +559,7 @@ void IMEXBDF2::constructSNES(SNES* snesIn) { } // 3D fields - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { int ind = ROUND(index(x, y, z)); @@ -1420,7 +1420,7 @@ void IMEXBDF2::loopVars(BoutReal* u) { if (mesh->firstX() && !mesh->periodicX) { for (int jx = 0; jx < mesh->xstart; ++jx) { for (int jy = mesh->ystart; jy <= mesh->yend; ++jy) { - for (int jz = 0; jz < mesh->LocalNz; ++jz) { + for (int jz = mesh->zstart; jz <= mesh->zend; ++jz) { op.run(jx, jy, jz, u); ++u; } @@ -1432,7 +1432,7 @@ void IMEXBDF2::loopVars(BoutReal* u) { if (mesh->lastX() && !mesh->periodicX) { for (int jx = mesh->xend + 1; jx < mesh->LocalNx; ++jx) { for (int jy = mesh->ystart; jy <= mesh->yend; ++jy) { - for (int jz = 0; jz < mesh->LocalNz; ++jz) { + for (int jz = mesh->zstart; jz <= mesh->zend; ++jz) { op.run(jx, jy, jz, u); ++u; } @@ -1442,7 +1442,7 @@ void IMEXBDF2::loopVars(BoutReal* u) { // Lower Y for (RangeIterator xi = mesh->iterateBndryLowerY(); !xi.isDone(); ++xi) { for (int jy = 0; jy < mesh->ystart; ++jy) { - for (int jz = 0; jz < mesh->LocalNz; ++jz) { + for (int jz = mesh->zstart; jz <= mesh->zend; ++jz) { op.run(*xi, jy, jz, u); ++u; } @@ -1452,7 +1452,7 @@ void IMEXBDF2::loopVars(BoutReal* u) { // Upper Y for (RangeIterator xi = mesh->iterateBndryUpperY(); !xi.isDone(); ++xi) { for (int jy = mesh->yend + 1; jy < mesh->LocalNy; ++jy) { - for (int jz = 0; jz < mesh->LocalNz; ++jz) { + for (int jz = mesh->zstart; jz <= mesh->zend; ++jz) { op.run(*xi, jy, jz, u); ++u; } @@ -1463,7 +1463,7 @@ void IMEXBDF2::loopVars(BoutReal* u) { // Bulk of points for (int jx = mesh->xstart; jx <= mesh->xend; ++jx) { for (int jy = mesh->ystart; jy <= mesh->yend; ++jy) { - for (int jz = 0; jz < mesh->LocalNz; ++jz) { + for (int jz = mesh->zstart; jz <= mesh->zend; ++jz) { op.run(jx, jy, jz, u); ++u; } diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index d0535f069d..81f1cc7951 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -660,7 +660,7 @@ int PetscSolver::init() { } } // 3D fields - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { const int ind = ROUND(index(x, y, z)) - Istart; for (int i = 0; i < n3d; i++) { @@ -768,7 +768,7 @@ int PetscSolver::init() { } } // 3D fields - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { int const ind = ROUND(index(x, y, z)); for (int i = 0; i < n3d; i++) { @@ -792,7 +792,7 @@ int PetscSolver::init() { || (yi >= mesh->LocalNy)) { continue; } - for (int zi = 0; zi < mesh->LocalNz; ++zi) { + for (int zi = mesh->zstart; zi <= mesh->zend; ++zi) { int ind2 = ROUND(index(xi, yi, zi)); if (ind2 < 0) { continue; // Boundary point diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 311bc371fe..446903dd4f 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -200,7 +200,7 @@ PetscErrorCode SNESSolver::FDJinitialise() { } } // 3D fields - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { const int ind = ROUND(index(x, y, z)) - Istart; for (int i = 0; i < n3d; i++) { @@ -305,7 +305,7 @@ PetscErrorCode SNESSolver::FDJinitialise() { } } // 3D fields - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { int ind = ROUND(index(x, y, z)); for (int i = 0; i < n3d; i++) { @@ -329,7 +329,7 @@ PetscErrorCode SNESSolver::FDJinitialise() { || (yi >= mesh->LocalNy)) { continue; } - for (int zi = 0; zi < mesh->LocalNz; ++zi) { + for (int zi = mesh->zstart; zi <= mesh->zend; ++zi) { int ind2 = ROUND(index(xi, yi, zi)); if (ind2 < 0) { continue; // Boundary point @@ -1331,7 +1331,7 @@ PetscErrorCode SNESSolver::updatePseudoTimestepping(Vec x) { } // Field3D quantities evolved together within a cell - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { count = 0; residual = 0.0; for (const auto& f : f3d) { @@ -1380,7 +1380,7 @@ PetscErrorCode SNESSolver::updatePseudoTimestepping(Vec x) { // Field3D quantities evolved together within a cell if (!f3d.empty()) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { auto i3d = mesh->ind2Dto3D(i2d, jz); BoutReal residual = 0.0; diff --git a/tests/integrated/test-invpar/test_invpar.cxx b/tests/integrated/test-invpar/test_invpar.cxx index 350437d66b..99b4e99c5c 100644 --- a/tests/integrated/test-invpar/test_invpar.cxx +++ b/tests/integrated/test-invpar/test_invpar.cxx @@ -47,7 +47,7 @@ int test(const std::string& acoef, const std::string& bcoef, const std::string& local_ystart = mesh->ystart + 1; } for (int y = local_ystart; y < mesh->yend; y++) { - for (int z = 0; z < mesh->LocalNz; z++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { output.write("result: [{:d},{:d}] : {:e}, {:e}, {:e}\n", y, z, input(mesh->xstart, y, z), result(mesh->xstart, y, z), deriv(mesh->xstart, y, z)); diff --git a/tests/integrated/test-multigrid_laplace/test_multigrid_laplace.cxx b/tests/integrated/test-multigrid_laplace/test_multigrid_laplace.cxx index 2b03f78049..a3e128bfdc 100644 --- a/tests/integrated/test-multigrid_laplace/test_multigrid_laplace.cxx +++ b/tests/integrated/test-multigrid_laplace/test_multigrid_laplace.cxx @@ -173,14 +173,14 @@ int main(int argc, char** argv) { // make field to pass in boundary conditions Field3D x0 = 0.; if (mesh->firstX()) { - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = mesh->zstart; k <= mesh->zend; k++) { x0(mesh->xstart - 1, mesh->ystart, k) = 0.5 * (f3(mesh->xstart - 1, mesh->ystart, k) + f3(mesh->xstart, mesh->ystart, k)); } } if (mesh->lastX()) { - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = mesh->zstart; k <= mesh->zend; k++) { x0(mesh->xend + 1, mesh->ystart, k) = 0.5 * (f3(mesh->xend + 1, mesh->ystart, k) + f3(mesh->xend, mesh->ystart, k)); } @@ -238,7 +238,7 @@ int main(int argc, char** argv) { // make field to pass in boundary conditions x0 = 0.; if (mesh->firstX()) { - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = mesh->zstart; k <= mesh->zend; k++) { x0(mesh->xstart - 1, mesh->ystart, k) = (f4(mesh->xstart, mesh->ystart, k) - f4(mesh->xstart - 1, mesh->ystart, k)) / mesh->getCoordinates()->dx(mesh->xstart, mesh->ystart, k) @@ -246,7 +246,7 @@ int main(int argc, char** argv) { } } if (mesh->lastX()) { - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = mesh->zstart; k <= mesh->zend; k++) { x0(mesh->xend + 1, mesh->ystart, k) = (f4(mesh->xend + 1, mesh->ystart, k) - f4(mesh->xend, mesh->ystart, k)) / mesh->getCoordinates()->dx(mesh->xend, mesh->ystart, k) @@ -311,7 +311,7 @@ BoutReal max_error_at_ystart(const Field3D& error) { BoutReal local_max_error = error(mesh->xstart, mesh->ystart, 0); for (int jx = mesh->xstart; jx <= mesh->xend; jx++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { if (local_max_error < error(jx, mesh->ystart, jz)) { local_max_error = error(jx, mesh->ystart, jz); } diff --git a/tests/integrated/test-naulin-laplace/test_naulin_laplace.cxx b/tests/integrated/test-naulin-laplace/test_naulin_laplace.cxx index 07a403e2e2..e9af4e5b36 100644 --- a/tests/integrated/test-naulin-laplace/test_naulin_laplace.cxx +++ b/tests/integrated/test-naulin-laplace/test_naulin_laplace.cxx @@ -175,14 +175,14 @@ int main(int argc, char** argv) { // make field to pass in boundary conditions Field3D x0 = 0.; if (mesh->firstX()) { - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = mesh->zstart; k <= mesh->zend; k++) { x0(mesh->xstart - 1, mesh->ystart, k) = 0.5 * (f3(mesh->xstart - 1, mesh->ystart, k) + f3(mesh->xstart, mesh->ystart, k)); } } if (mesh->lastX()) { - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = mesh->zstart; k <= mesh->zend; k++) { x0(mesh->xend + 1, mesh->ystart, k) = 0.5 * (f3(mesh->xend + 1, mesh->ystart, k) + f3(mesh->xend, mesh->ystart, k)); } @@ -242,7 +242,7 @@ int main(int argc, char** argv) { // make field to pass in boundary conditions x0 = 0.; if (mesh->firstX()) { - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = mesh->zstart; k <= mesh->zend; k++) { x0(mesh->xstart - 1, mesh->ystart, k) = (f4(mesh->xstart, mesh->ystart, k) - f4(mesh->xstart - 1, mesh->ystart, k)) / mesh->getCoordinates()->dx(mesh->xstart, mesh->ystart, k) @@ -250,7 +250,7 @@ int main(int argc, char** argv) { } } if (mesh->lastX()) { - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = mesh->zstart; k <= mesh->zend; k++) { x0(mesh->xend + 1, mesh->ystart, k) = (f4(mesh->xend + 1, mesh->ystart, k) - f4(mesh->xend, mesh->ystart, k)) / mesh->getCoordinates()->dx(mesh->xend, mesh->ystart, k) @@ -315,7 +315,7 @@ BoutReal max_error_at_ystart(const Field3D& error) { BoutReal local_max_error = error(mesh->xstart, mesh->ystart, 0); for (int jx = mesh->xstart; jx <= mesh->xend; jx++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { if (local_max_error < error(jx, mesh->ystart, jz)) { local_max_error = error(jx, mesh->ystart, jz); } diff --git a/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx b/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx index 1e3cdde310..c42c55d8d6 100644 --- a/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx +++ b/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx @@ -222,7 +222,7 @@ BoutReal max_error_at_ystart(const Field3D& error) { BoutReal local_max_error = error(mesh->xstart, mesh->ystart, 0); for (int jx = mesh->xstart; jx <= mesh->xend; jx++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { if (local_max_error < error(jx, mesh->ystart, jz)) { local_max_error = error(jx, mesh->ystart, jz); } @@ -241,7 +241,7 @@ void apply_flat_boundary(Field3D& bcoef) { if (mesh.firstX()) { for (int jx = mesh.xstart - 1; jx >= 0; jx--) { for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { bcoef(jx, jy, jz) = bcoef(jx + 1, jy, jz); } } @@ -250,7 +250,7 @@ void apply_flat_boundary(Field3D& bcoef) { if (mesh.lastX()) { for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { bcoef(jx, jy, jz) = bcoef(jx - 1, jy, jz); } } @@ -270,7 +270,7 @@ Field3D generate_f1(const Mesh& mesh) { for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; //make the gradients zero at both x-boundaries result(jx, jy, jz) = 0. @@ -289,7 +289,7 @@ Field3D generate_f1(const Mesh& mesh) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; //make the gradients zero at both x-boundaries result(jx, jy, jz) = 0. @@ -308,7 +308,7 @@ Field3D generate_f1(const Mesh& mesh) { for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; //make the gradients zero at both x-boundaries result(jx, jy, jz) = 0. @@ -340,7 +340,7 @@ Field3D generate_d1(const Mesh& mesh) { for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = 1. + 0.2 * exp(-50. * pow(x - p, 2) / 4.) * sin(2. * PI * (z - q) * 3.); @@ -351,7 +351,7 @@ Field3D generate_d1(const Mesh& mesh) { for (int jx = mesh.xstart - 1; jx >= 0; jx--) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = 1. + 0.2 * exp(-50. * pow(x - p, 2) / 4.) * sin(2. * PI * (z - q) * 3.); @@ -363,7 +363,7 @@ Field3D generate_d1(const Mesh& mesh) { for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = 1. + 0.2 * exp(-50. * pow(x - p, 2) / 4.) * sin(2. * PI * (z - q) * 3.); @@ -386,7 +386,7 @@ Field3D generate_c1(const Mesh& mesh) { for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = 1. + 0.15 * exp(-50. * pow(x - p, 2) * 2.) * sin(2. * PI * (z - q) * 2.); @@ -397,7 +397,7 @@ Field3D generate_c1(const Mesh& mesh) { for (int jx = mesh.xstart - 1; jx >= 0; jx--) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = 1. + 0.15 * exp(-50. * pow(x - p, 2) * 2.) * sin(2. * PI * (z - q) * 2.); @@ -409,7 +409,7 @@ Field3D generate_c1(const Mesh& mesh) { for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = 1. + 0.15 * exp(-50. * pow(x - p, 2) * 2.) * sin(2. * PI * (z - q) * 2.); @@ -433,7 +433,7 @@ Field3D generate_a1(const Mesh& mesh) { for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = -1. + 0.1 * exp(-50. * pow(x - p, 2) * 2.5) * sin(2. * PI * (z - q) * 7.); @@ -444,7 +444,7 @@ Field3D generate_a1(const Mesh& mesh) { for (int jx = mesh.xstart - 1; jx >= 0; jx--) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = -1. + 0.1 * exp(-50. * pow(x - p, 2) * 2.5) * sin(2. * PI * (z - q) * 7.); @@ -456,7 +456,7 @@ Field3D generate_a1(const Mesh& mesh) { for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = -1. + 0.1 * exp(-50. * pow(x - p, 2) * 2.5) * sin(2. * PI * (z - q) * 7.); @@ -479,7 +479,7 @@ Field3D generate_f5(const Mesh& mesh) { for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; //make the gradients zero at both x-boundaries result(jx, jy, jz) = @@ -496,7 +496,7 @@ Field3D generate_f5(const Mesh& mesh) { for (int jx = mesh.xstart - 1; jx >= 0; jx--) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; //make the gradients zero at both x-boundaries result(jx, jy, jz) = 0. @@ -515,7 +515,7 @@ Field3D generate_f5(const Mesh& mesh) { for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; //make the gradients zero at both x-boundaries result(jx, jy, jz) = 0. @@ -545,7 +545,7 @@ Field3D generate_d5(const Mesh& mesh) { for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = 1. + p * cos(2. * PI * x) * sin(2. * PI * (z - q) * 3.); } @@ -555,7 +555,7 @@ Field3D generate_d5(const Mesh& mesh) { for (int jx = mesh.xstart - 1; jx >= 0; jx--) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = 1. + p * cos(2. * PI * x) * sin(2. * PI * (z - q) * 3.); } @@ -566,7 +566,7 @@ Field3D generate_d5(const Mesh& mesh) { for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = 1. + p * cos(2. * PI * x) * sin(2. * PI * (z - q) * 3.); } @@ -589,7 +589,7 @@ Field3D generate_c5(const Mesh& mesh) { for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = 1. + p * cos(2. * PI * x * 5) * sin(2. * PI * (z - q) * 2.); } @@ -599,7 +599,7 @@ Field3D generate_c5(const Mesh& mesh) { for (int jx = mesh.xstart - 1; jx >= 0; jx--) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = 1. + p * cos(2. * PI * x * 5) * sin(2. * PI * (z - q) * 2.); @@ -611,7 +611,7 @@ Field3D generate_c5(const Mesh& mesh) { for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = 1. + p * cos(2. * PI * x * 5) * sin(2. * PI * (z - q) * 2.); @@ -633,7 +633,7 @@ Field3D generate_a5(const Mesh& mesh) { for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = -1. + p * cos(2. * PI * x * 2.) * sin(2. * PI * (z - q) * 7.); @@ -644,7 +644,7 @@ Field3D generate_a5(const Mesh& mesh) { for (int jx = mesh.xstart - 1; jx >= 0; jx--) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = -1. + p * cos(2. * PI * x * 2.) * sin(2. * PI * (z - q) * 7.); @@ -656,7 +656,7 @@ Field3D generate_a5(const Mesh& mesh) { for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < mesh.LocalNz; jz++) { + for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = -1. + p * cos(2. * PI * x * 2.) * sin(2. * PI * (z - q) * 7.); diff --git a/tests/integrated/test-petsc_laplace_MAST-grid/test_petsc_laplace_MAST_grid.cxx b/tests/integrated/test-petsc_laplace_MAST-grid/test_petsc_laplace_MAST_grid.cxx index 7646b915a7..945d5c6443 100644 --- a/tests/integrated/test-petsc_laplace_MAST-grid/test_petsc_laplace_MAST_grid.cxx +++ b/tests/integrated/test-petsc_laplace_MAST-grid/test_petsc_laplace_MAST_grid.cxx @@ -63,7 +63,7 @@ int main(int argc, char** argv) { f1.allocate(); for (int jx = mesh->xstart; jx <= mesh->xend; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; f1(jx, jy, jz) = @@ -84,7 +84,7 @@ int main(int argc, char** argv) { if (mesh->firstX()) { for (int jx = mesh->xstart - 1; jx >= 0; jx--) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; f1(jx, jy, jz) = @@ -106,7 +106,7 @@ int main(int argc, char** argv) { if (mesh->lastX()) { for (int jx = mesh->xend + 1; jx < mesh->LocalNx; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; f1(jx, jy, jz) = @@ -131,7 +131,7 @@ int main(int argc, char** argv) { d1.allocate(); for (int jx = mesh->xstart; jx <= mesh->xend; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; d1(jx, jy, jz) = @@ -143,7 +143,7 @@ int main(int argc, char** argv) { if (mesh->firstX()) { for (int jx = mesh->xstart - 1; jx >= 0; jx--) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; d1(jx, jy, jz) = @@ -157,7 +157,7 @@ int main(int argc, char** argv) { if (mesh->lastX()) { for (int jx = mesh->xend + 1; jx < mesh->LocalNx; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; d1(jx, jy, jz) = @@ -174,7 +174,7 @@ int main(int argc, char** argv) { c1.allocate(); for (int jx = mesh->xstart; jx <= mesh->xend; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; c1(jx, jy, jz) = 1. @@ -186,7 +186,7 @@ int main(int argc, char** argv) { if (mesh->firstX()) { for (int jx = mesh->xstart - 1; jx >= 0; jx--) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; c1(jx, jy, jz) = 1. @@ -199,7 +199,7 @@ int main(int argc, char** argv) { if (mesh->lastX()) { for (int jx = mesh->xend + 1; jx < mesh->LocalNx; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; c1(jx, jy, jz) = 1. @@ -215,7 +215,7 @@ int main(int argc, char** argv) { a1.allocate(); for (int jx = mesh->xstart; jx <= mesh->xend; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; a1(jx, jy, jz) = @@ -226,7 +226,7 @@ int main(int argc, char** argv) { if (mesh->firstX()) { for (int jx = mesh->xstart - 1; jx >= 0; jx--) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; a1(jx, jy, jz) = @@ -238,7 +238,7 @@ int main(int argc, char** argv) { if (mesh->lastX()) { for (int jx = mesh->xend + 1; jx < mesh->LocalNx; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; a1(jx, jy, jz) = @@ -256,7 +256,7 @@ int main(int argc, char** argv) { if (mesh->firstX()) { for (int jx = mesh->xstart - 1; jx >= 0; jx--) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { b1(jx, jy, jz) = b1(jx + 1, jy, jz); } } @@ -265,7 +265,7 @@ int main(int argc, char** argv) { if (mesh->lastX()) { for (int jx = mesh->xend + 1; jx < mesh->LocalNx; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { b1(jx, jy, jz) = b1(jx - 1, jy, jz); } } @@ -371,7 +371,7 @@ int main(int argc, char** argv) { if (mesh->firstX()) { for (int jx = mesh->xstart - 1; jx >= 0; jx--) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { b3(jx, jy, jz) = b3(jx + 1, jy, jz); } } @@ -380,7 +380,7 @@ int main(int argc, char** argv) { if (mesh->lastX()) { for (int jx = mesh->xend + 1; jx < mesh->LocalNx; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { b3(jx, jy, jz) = b3(jx - 1, jy, jz); } } @@ -468,7 +468,7 @@ int main(int argc, char** argv) { f5.allocate(); for (int jx = mesh->xstart; jx <= mesh->xend; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; f5(jx, jy, jz) = @@ -489,7 +489,7 @@ int main(int argc, char** argv) { if (mesh->firstX()) { for (int jx = mesh->xstart - 1; jx >= 0; jx--) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; f5(jx, jy, jz) = @@ -511,7 +511,7 @@ int main(int argc, char** argv) { if (mesh->lastX()) { for (int jx = mesh->xend + 1; jx < mesh->LocalNx; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; f5(jx, jy, jz) = @@ -536,7 +536,7 @@ int main(int argc, char** argv) { d5.allocate(); for (int jx = mesh->xstart; jx <= mesh->xend; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; d5(jx, jy, jz) = @@ -547,7 +547,7 @@ int main(int argc, char** argv) { if (mesh->firstX()) { for (int jx = mesh->xstart - 1; jx >= 0; jx--) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; d5(jx, jy, jz) = @@ -559,7 +559,7 @@ int main(int argc, char** argv) { if (mesh->lastX()) { for (int jx = mesh->xend + 1; jx < mesh->LocalNx; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; d5(jx, jy, jz) = @@ -574,7 +574,7 @@ int main(int argc, char** argv) { c5.allocate(); for (int jx = mesh->xstart; jx <= mesh->xend; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; c5(jx, jy, jz) = @@ -585,7 +585,7 @@ int main(int argc, char** argv) { if (mesh->firstX()) { for (int jx = mesh->xstart - 1; jx >= 0; jx--) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; c5(jx, jy, jz) = @@ -597,7 +597,7 @@ int main(int argc, char** argv) { if (mesh->lastX()) { for (int jx = mesh->xend + 1; jx < mesh->LocalNx; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; c5(jx, jy, jz) = @@ -612,7 +612,7 @@ int main(int argc, char** argv) { a5.allocate(); for (int jx = mesh->xstart; jx <= mesh->xend; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; a5(jx, jy, jz) = -1. + p * cos(2. * PI * x * 2.) * sin(2. * PI * (z - q) * 7.); @@ -622,7 +622,7 @@ int main(int argc, char** argv) { if (mesh->firstX()) { for (int jx = mesh->xstart - 1; jx >= 0; jx--) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; a5(jx, jy, jz) = @@ -634,7 +634,7 @@ int main(int argc, char** argv) { if (mesh->lastX()) { for (int jx = mesh->xend + 1; jx < mesh->LocalNx; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { BoutReal x = BoutReal(mesh->getGlobalXIndex(jx) - mesh->xstart) / nx; BoutReal z = BoutReal(jz) / nz; a5(jx, jy, jz) = @@ -650,7 +650,7 @@ int main(int argc, char** argv) { if (mesh->firstX()) { for (int jx = mesh->xstart - 1; jx >= 0; jx--) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { b5(jx, jy, jz) = b5(jx + 1, jy, jz); } } @@ -659,7 +659,7 @@ int main(int argc, char** argv) { if (mesh->lastX()) { for (int jx = mesh->xend + 1; jx < mesh->LocalNx; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { b5(jx, jy, jz) = b5(jx - 1, jy, jz); } } @@ -762,7 +762,7 @@ int main(int argc, char** argv) { if (mesh->firstX()) { for (int jx = mesh->xstart - 1; jx >= 0; jx--) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { b7(jx, jy, jz) = b7(jx + 1, jy, jz); } } @@ -771,7 +771,7 @@ int main(int argc, char** argv) { if (mesh->lastX()) { for (int jx = mesh->xend + 1; jx < mesh->LocalNx; jx++) { for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { b7(jx, jy, jz) = b7(jx - 1, jy, jz); } } @@ -864,7 +864,7 @@ BoutReal max_error_at_ystart(const Field3D& error) { BoutReal local_max_error = error(mesh->xstart, mesh->ystart, 0); for (int jx = mesh->xstart; jx <= mesh->xend; jx++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { if (local_max_error < error(jx, mesh->ystart, jz)) { local_max_error = error(jx, mesh->ystart, jz); }