Merge pull request #2713 from boutproject/petsc-interface-array

Refactor `PetscVector` to use buffered writes

.clang-tidy

@@ -359,5 +359,11 @@ CheckOptions:
     value:           'false'
   - key:             readability-simplify-subscript-expr.Types
     value:           '::std::basic_string;::std::basic_string_view;::std::vector;::std::array'
+  - key:             readability-identifier-length.IgnoredVariableNames
+    value:           '^n?[xyz]$'
+  - key:             readability-identifier-length.IgnoredParameterNames
+    value:           '^[ijkxyz][01]?$'
+  - key:             readability-identifier-length.IgnoredLoopCounterNames
+    value:           '^[ijkxyz_]$'
 ...
 

include/bout/field.hxx

@@ -122,6 +122,9 @@ public:
    */
   virtual int getNz() const;
 
+  /// Get the total number of points
+  virtual int size() const = 0;
+
   friend void swap(Field& first, Field& second) noexcept {
     using std::swap;
     swap(static_cast<FieldData&>(first), static_cast<FieldData&>(second));

include/bout/field2d.hxx

@@ -277,10 +277,12 @@ public:
 
   friend void swap(Field2D& first, Field2D& second) noexcept;
 
+  int size() const override { return nx * ny; };
+
+private:
   /// Internal data array. Handles allocation/freeing of memory
   Array<BoutReal> data;
 
-private:
   /// Array sizes (from fieldmesh). These are valid only if fieldmesh is not null
   int nx{-1}, ny{-1};
 

include/bout/field3d.hxx

@@ -489,6 +489,8 @@ public:
 
   friend void swap(Field3D& first, Field3D& second) noexcept;
 
+  int size() const override { return nx * ny * nz; };
+
 private:
   /// Array sizes (from fieldmesh). These are valid only if fieldmesh is not null
   int nx{-1}, ny{-1}, nz{-1};

include/bout/fieldperp.hxx

@@ -287,6 +287,8 @@ public:
 
   friend void swap(FieldPerp& first, FieldPerp& second) noexcept;
 
+  int size() const override { return nx * nz; };
+
 private:
   /// The Y index at which this FieldPerp is defined
   int yindex{-1};

include/bout/globalindexer.hxx

@@ -45,7 +45,8 @@ public:
                          bool autoInitialise = true)
       : fieldmesh(localmesh), indices(-1., localmesh), stencils(std::move(stencil)) {
 
-    Region<ind_type> allCandidate, bndryCandidate;
+    Region<ind_type> allCandidate;
+    Region<ind_type> bndryCandidate;
     if (stencils.getNumParts() > 0) {
       std::set<ind_type> allIndices(getRegionNobndry().getIndices().begin(),
                                     getRegionNobndry().getIndices().end());
@@ -105,7 +106,18 @@ public:
   /// Finish setting up the indexer, communicating indices across
   /// processes and, if possible, calculating the sparsity pattern of
   /// any matrices.
-  void initialise() { fieldmesh->communicate(indices); }
+  void initialise() {
+    // We need to ensure any _guard_ cells are -1 so we don't include them
+    int_indices.reallocate(indices.size());
+    BOUT_FOR(index, indices.getRegion("RGN_GUARDS")) { int_indices[index.ind] = -1; }
+    // Now we can communicate to get global indices from neighbouring processes
+    fieldmesh->communicate(indices);
+    // Finally, we fill in the global indices including in the
+    // _boundaries_ (*not* guards)
+    BOUT_FOR(index, regionAll) {
+      int_indices[index.ind] = static_cast<int>(indices[index]);
+    }
+  }
 
   Mesh* getMesh() const { return fieldmesh; }
 
@@ -157,6 +169,8 @@ public:
 
   int size() const { return regionAll.size(); }
 
+  const Array<int>& getIntIndices() const { return int_indices; }
+
 protected:
   // Must not be const as the index field needs to be mutable in order
   // to fake parallel communication in the unit tests.
@@ -206,6 +220,8 @@ private:
 
   /// Fields containing the indices for each element (as reals)
   T indices;
+  /// Indices as integers
+  Array<int> int_indices;
   /// The first and last global index on this processor (inclusive in
   /// both cases)
   int globalStart, globalEnd;