Require Mesh::addCoordinates(location) to initialize Coordinates at location

examples/staggered_grid/data/BOUT.inp

@@ -6,6 +6,8 @@ MZ = 1
 
 grid = "test-staggered.nc"
 
+[mesh]
+
 StaggerGrids = true
 
 [mesh:ddy]

examples/staggered_grid/test/BOUT.inp

@@ -4,9 +4,9 @@ timestep = 0.02
 
 MZ = 1
 
+[mesh]
 StaggerGrids = true
 
-[mesh]
 nx = 5
 ny = 16
 

examples/staggered_grid/test_staggered.cxx

@@ -8,10 +8,17 @@
 #include <derivs.hxx>
 
 Field3D n, v;
+CELL_LOC maybe_ylow = CELL_CENTRE;
 
 int physics_init(bool restart) {
 
-  v.setLocation(CELL_YLOW); // Staggered relative to n
+  if (mesh->StaggerGrids) {
+    maybe_ylow = CELL_YLOW;
+
+    mesh->addCoordinates(CELL_YLOW);
+
+    v.setLocation(CELL_YLOW); // Staggered relative to n
+  }
 
   SOLVE_FOR(n, v);
 
@@ -24,7 +31,7 @@ int physics_run(BoutReal time) {
   //ddt(n) = -Div_par_flux(v, n, CELL_CENTRE);
   ddt(n) = -n*Grad_par(v, CELL_CENTRE) - Vpar_Grad_par(v, n, CELL_CENTRE);
 
-  ddt(v) = -Grad_par(n, CELL_YLOW);
+  ddt(v) = -Grad_par(n, maybe_ylow);
 
   // Have to manually apply the lower Y boundary region, using a width of 3
   for( RangeIterator rlow = mesh->iterateBndryLowerY(); !rlow.isDone(); rlow++)

include/bout/coordinates.hxx

@@ -102,7 +102,7 @@ public:
   Field2D IntShiftTorsion; ///< Integrated shear (I in BOUT notation)
 
   /// Calculate differential geometry quantities from the metric tensor
-  int geometry();
+  int geometry(bool allow_geometry_without_recalculate_staggered = false);
   int calcCovariant(); ///< Inverts contravatiant metric to get covariant
   int calcContravariant(); ///< Invert covariant metric to get contravariant
   int jacobian(); ///< Calculate J and Bxy

include/bout/mesh.hxx

@@ -64,6 +64,7 @@ class Mesh;
 #include "sys/range.hxx" // RangeIterator
 
 #include <bout/griddata.hxx>
+#include <bout/macro_for_each.hxx>
 
 #include "coordinates.hxx"    // Coordinates class
 
@@ -438,15 +439,30 @@ class Mesh {
     ASSERT1(location != CELL_DEFAULT);
     ASSERT1(location != CELL_VSHIFT);
 
-    if (coords_map.count(location)) { // True branch most common, returns immediately
-      return coords_map[location].get();
-    } else {
-      // No coordinate system set. Create default
-      // Note that this can't be allocated here due to incomplete type
-      // (circular dependency between Mesh and Coordinates)
-      coords_map.emplace(location, createDefaultCoordinates(location));
-      return coords_map[location].get();
+#if CHECK > 0
+    if (!coords_map.count(location)) {
+      throw BoutException("Error: Coordinates for %s have not been added to "
+          "this Mesh. You should call the method Mesh::addCoordinates(location) "
+          "before initializing fields staggered to 'CELL_LOC location'.",
+          CELL_LOC_STRING(location).c_str());
     }
+#endif
+    return coords_map.at(location).get();
+  }
+
+  /// Add Coordinates object at a certain location.
+  /// If replace_coords is set to true, reset the object in coords_map if it
+  /// already exists, otherwise add a new one
+  void addCoordinates(const CELL_LOC location, bool replace_coords = false);
+
+  /// Check if Coordinates object at location has been added
+  bool hasCoordinates(const CELL_LOC location) {
+    return coords_map.count(location);
+  }
+
+  /// Count how many Coordinates objects have been added
+  int countCoordinates() {
+    return coords_map.size();
   }
 
   /// Returns the non-CELL_CENTRE location
@@ -790,7 +806,7 @@ class Mesh {
 
   GridDataSource *source; ///< Source for grid data
 
-  std::map<CELL_LOC, std::shared_ptr<Coordinates> > coords_map; ///< Coordinate systems at different CELL_LOCs
+  std::map<CELL_LOC, std::unique_ptr<Coordinates> > coords_map; ///< Coordinate systems at different CELL_LOCs
 
   Options *options; ///< Mesh options section
 
@@ -808,9 +824,6 @@ class Mesh {
 
 private:
 
-  /// Allocates default Coordinates objects
-  std::shared_ptr<Coordinates> createDefaultCoordinates(const CELL_LOC location);
-
   //Internal region related information
   std::map<std::string, Region<Ind3D>> regionMap3D;
   std::map<std::string, Region<Ind2D>> regionMap2D;

manual/sphinx/user_docs/staggered_grids.rst

@@ -16,8 +16,9 @@ staggered grids, set::
 
     StaggerGrids = true
 
-in the top section of the ``BOUT.inp`` file. The **test-staggered**
-example illustrates how to use staggered grids in BOUT++.
+in the top section of the ``BOUT.inp`` file and enable the locations you will
+use with a call to mesh->addCoordinates(location) (see below).  The
+**test-staggered** example illustrates how to use staggered grids in BOUT++.
 
 There are four possible locations in a grid cell where a quantity can be
 defined in BOUT++: centre, lower X, lower Y, and lower Z. These are
@@ -29,17 +30,22 @@ illustrated in :numref:`staggergrids-location`.
 
    The four possible cell locations for defining quantities
 
-To specify the location of a variable, use the method
-`Field3D::setLocation` with one of the `CELL_LOC` locations
-`CELL_CENTRE`, `CELL_XLOW`, `CELL_YLOW`, or `CELL_ZLOW`.
+The possible locations are specified with the `CELL_LOC` type, which has the
+possible values `CELL_CENTRE`, `CELL_XLOW`, `CELL_YLOW`, or `CELL_ZLOW`.
+`CELL_CENTRE` is enabled by default, but a call to
+mesh->addCoordinates(location) is required to enable the others.
 
 The key lines in the **staggered_grid** example which specify the
 locations of the evolving variables are::
 
     Field3D n, v;
 
     int init(bool restart) {
+
+      mesh->addCoordinates(CELL_YLOW);
+
       v.setLocation(CELL_YLOW); // Staggered relative to n
+
       SOLVE_FOR(n, v);
       ...
 
@@ -62,6 +68,33 @@ in Y, whilst the density :math:`n` remains cell centred.
           `CELL_CENTRE` if staggered grids are off, regardless of what
           you pass it.
 
+.. note:: For advanced users:
+          If you change members of the Coordinates object manually, you should
+          change the CELL_CENTRE Coordinates and only call addCoordinates() for
+          other locations after you call Coordinates::geometry() on the
+          CELL_CENTRE Coordinates. Then the Coordinates at staggered locations
+          will be interpolated from the correct, final CELL_CENTRE version.
+
+          The example uses the global Mesh object 'mesh'. If you are using any
+          other Mesh objects, you need to initialize the Coordinates objects in
+          their coords_map members by calling the addCoordinates(CELL_LOC
+          location) method for each location you will use.
+
+          An exception will be thrown if you call Coordinates::geometry() from
+          any Coordinates object at a staggered location, since these are
+          expected to be consistent with (and calculated from) the CELL_CENTRE
+          Coordinates. If you need to change them, you can set the option
+          mesh:allow_geometry_without_recalculate_staggered=true to disable
+          this check; you must then ensure that all the Coordinates objects in
+          Mesh::coords_map are consistent with each other. Setting this option
+          also allows changes to be made and geometry() to be called on the
+          CELL_CENTRE Coordinates after other locations have been added to
+          coords_map; in this case you will need to update the other locations
+          explicitly by calling Mesh::addCoordinates(location, true) - the
+          optional second argument causes addCoordinates to overwrite any
+          existing Coordinates object at location and replace it with one
+          calculated from the current CELL_CENTRE Coordinates.
+
 
 Arithmetic operations can only be performed between variables with the same
 location. When performing a calculation at one location, to include a variable

src/field/field3d.cxx

@@ -209,6 +209,12 @@ void Field3D::setLocation(CELL_LOC new_location) {
       new_location = CELL_CENTRE;
     }
 
+#if CHECK > 1
+    // Check Coordinates for location have been added
+    // For CHECK > 0, getCoordinates will throw if location has not been added.
+    getMesh()->getCoordinates(location);
+#endif
+
     // Invalidate the coordinates pointer
     if (new_location != location) {
       fieldCoordinates = nullptr;

src/mesh/coordinates.cxx

@@ -324,7 +324,7 @@ void Coordinates::outputVars(Datafile &file) {
   file.add(J, "J", false);
 }
 
-int Coordinates::geometry() {
+int Coordinates::geometry(bool allow_geometry_without_recalculate_staggered) {
   TRACE("Coordinates::geometry");
 
   output_progress.write("Calculating differential geometry terms\n");
@@ -462,13 +462,16 @@ int Coordinates::geometry() {
 
   OPTION(Options::getRoot(), non_uniform, true);
 
-  Field2D d2x, d2y; // d^2 x / d i^2
+  Field2D d2x(localmesh), d2y(localmesh); // d^2 x / d i^2
   // Read correction for non-uniform meshes
   if (localmesh->get(d2x, "d2x")) {
     output_warn.write(
         "\tWARNING: differencing quantity 'd2x' not found. Calculating from dx\n");
     d1_dx = bout::derivatives::index::DDX(1. / dx); // d/di(1/dx)
   } else {
+    // Shift d2x to our location
+    d2x = interp_to(d2x, location);
+
     d1_dx = -d2x / (dx * dx);
   }
 
@@ -477,9 +480,43 @@ int Coordinates::geometry() {
         "\tWARNING: differencing quantity 'd2y' not found. Calculating from dy\n");
     d1_dy = bout::derivatives::index::DDY(1. / dy); // d/di(1/dy)
   } else {
+    // Shift d2y to our location
+    d2y = interp_to(d2y, location);
+
     d1_dy = -d2y / (dy * dy);
   }
 
+  if (!allow_geometry_without_recalculate_staggered) {
+    // Only CELL_CENTRE Coordinates should ever be changed (and therefore need
+    // geometry() calling). If CELL_CENTRE Coordinates are changed they should
+    // be changed before other Coordinates are added to localmesh->coords_map,
+    // since Coordinates at other locations are calculated from the CELL_CENTRE
+    // ones. The allow_geometry_without_recalculate_staggered option overrides
+    // this check, in case for some reason the user does need to call
+    // geometry() in other situations.
+    if ( (localmesh->hasCoordinates(location))                                            // this location already added
+         && !(localmesh->countCoordinates()==1 && localmesh->hasCoordinates(CELL_CENTRE)) // OK to be added already only if CELL_CENTRE and no other locations added
+       ) {
+      if (location == CELL_CENTRE) {
+        throw BoutException("Coordinates::geometry() called at CELL_CENTRE, but "
+          "other locations have already been added to the Mesh. These would "
+          "need recalculating. If possible call Mesh::addCoordinates() after "
+          "this call to geometry(). If you need to recalculate multiple "
+          "Coordinates objects explicitly, call geometry(true) [i.e. setting the "
+          "argument allow_geometry_without_recalculate_staggered=true] file to "
+          "disable this check.");
+      } else {
+        throw BoutException("Coordinates::geometry() called at location %s, but "
+          "this location has already been initialized and added to the Mesh. "
+          "If you need to recalculate multiple Coordinates objects explicitly, "
+          "call geometry(true) [i.e. setting the argument "
+          "allow_geometry_without_recalculate_staggered=true] input file to "
+          "disable this check.", CELL_LOC_STRING(location).c_str());
+      }
+    }
+
+  }
+
   return 0;
 }
 

src/mesh/impls/bout/boutmesh.cxx

@@ -839,7 +839,10 @@ int BoutMesh::load() {
   // Add boundary regions
   addBoundaryRegions();
 
-  output_info.write(_("\tdone\n"));
+  // Create CELL_CENTRE Coordinates object
+  addCoordinates(CELL_CENTRE);
+
+  output_info.write("\tdone\n");
 
   return 0;
 }
@@ -2208,7 +2211,7 @@ void BoutMesh::addBoundaryRegions() {
   all_boundaries.emplace_back("RGN_UPPER_Y");
 
   // Inner X
-  if(mesh->firstX() && !mesh->periodicX) {
+  if(firstX() && !periodicX) {
     addRegion3D("RGN_INNER_X", Region<Ind3D>(0, xstart-1, ystart, yend, 0, LocalNz-1,
                                              LocalNy, LocalNz, maxregionblocksize));
     addRegion2D("RGN_INNER_X", Region<Ind2D>(0, xstart-1, ystart, yend, 0, 0,
@@ -2225,7 +2228,7 @@ void BoutMesh::addBoundaryRegions() {
   }
 
   // Outer X
-  if(mesh->firstX() && !mesh->periodicX) {
+  if(firstX() && !periodicX) {
     addRegion3D("RGN_OUTER_X", Region<Ind3D>(xend+1, LocalNx-1, ystart, yend, 0, LocalNz-1,
                                              LocalNy, LocalNz, maxregionblocksize));
     addRegion2D("RGN_OUTER_X", Region<Ind2D>(xend+1, LocalNx-1, ystart, yend, 0, 0,

src/mesh/mesh.cxx

@@ -325,13 +325,44 @@ ParallelTransform& Mesh::getParallelTransform() {
   return *transform;
 }
 
-std::shared_ptr<Coordinates> Mesh::createDefaultCoordinates(const CELL_LOC location) {
-  if (location == CELL_CENTRE || location == CELL_DEFAULT)
-    // Initialize coordinates from input
-    return std::make_shared<Coordinates>(this);
-  else
-    // Interpolate coordinates from CELL_CENTRE version
-    return std::make_shared<Coordinates>(this, location, getCoordinates(CELL_CENTRE));
+void Mesh::addCoordinates(const CELL_LOC location, bool replace_coords) {
+  ASSERT1(location != CELL_DEFAULT);
+
+  if (location == CELL_VSHIFT) {
+    // CELL_VSHIFT puts vector components at CELL_XLOW, CELL_YLOW and
+    // CELL_ZLOW, so require Coordinates at all three.
+    addCoordinates(CELL_XLOW, replace_coords);
+    addCoordinates(CELL_YLOW, replace_coords);
+    addCoordinates(CELL_ZLOW, replace_coords);
+  } else {
+    // No coordinate system set. Create default
+    if (location == CELL_CENTRE) {
+      // Initialize coordinates from input
+      if (!coords_map.count(location)) {
+        // location does not exist in coords_map, so create new entry
+        coords_map.emplace(location, bout::utils::make_unique<Coordinates>(this));
+      } else if (replace_coords) {
+        // location does already exist in coords_map, so reset it
+        coords_map.at(location).reset(new Coordinates(this));
+      }
+    } else {
+      // Interpolate coordinates from CELL_CENTRE version
+      ASSERT1(StaggerGrids); // If StaggerGrids==false, it doesn't make sense to have non-CELL_CENTRE Coordinates
+
+      if (!replace_coords and (coords_map.count(location) > 0)) {
+        throw BoutException("Coordinates at %s already added to Mesh",
+            CELL_LOC_STRING(location).c_str());
+      }
+      if (replace_coords) {
+        // first erase the existing entry to avoid throwing an exception from
+        // Coordinates::geometry(): replacement of the Coordinates object has
+        // been explicitly requested.
+        coords_map.erase(location);
+      }
+
+      coords_map.emplace(location, bout::utils::make_unique<Coordinates>(this, location, getCoordinates(CELL_CENTRE)));
+    }
+  }
 }
 
 

tests/MMS/derivatives3/runtest

@@ -33,6 +33,8 @@ def runtests(functions,derivatives,directions,stag,msg):
                                ,"2*pi/(%d)"%(nz),force=True)
             dirnfac=direction+"*"+fac
             mesh=boutcore.Mesh(section="mesh"+direction)
+            for loc in locations[1:]:
+                mesh.addCoordinates(loc)
             f=boutcore.create3D(infunc.replace("%s",dirnfac),mesh
                                 ,outloc=inloc)
             sim=diff_func(f,method=diff,outloc=outloc)

tests/MMS/upwinding3/runtest

@@ -34,6 +34,8 @@ def runtests(functions,derivatives,directions,stag,msg):
                                ,"2*pi/(%d)"%(nz),force=True)
             dirnfac=direction+"*"+fac
             mesh=boutcore.Mesh(section="mesh"+direction)
+            for loc in locations[1:]:
+                mesh.addCoordinates(loc)
             f=boutcore.create3D(ffunc.replace("%s",dirnfac),mesh
                                 ,outloc=floc)
             v=boutcore.create3D(vfunc.replace("%s",dirnfac),mesh

tests/MMS/wave-1d-y/wave.cxx

@@ -9,6 +9,7 @@ class Wave1D : public PhysicsModel {
 protected:
   int init(bool restarting) {
 
+    mesh->addCoordinates(CELL_YLOW);
     g.setLocation(CELL_YLOW); // g staggered 
 
     // Tell BOUT++ to solve f and g

tests/MMS/wave-1d/wave.cxx

@@ -83,6 +83,7 @@ class Wave1D : public PhysicsModel {
     coord->g_23 = 0.0;
     coord->geometry();
 
+    mesh->addCoordinates(CELL_XLOW);
     g.setLocation(CELL_XLOW); // g staggered to the left of f
 
     //Dirichlet everywhere except inner x-boundary Neumann