Use Mesh::addCoordinates(location) to initialize Coordinates

include/bout/coordinates.hxx

@@ -143,6 +143,9 @@ public:
   const Field3D Laplace(const Field3D &f, CELL_LOC outloc=CELL_DEFAULT);
 
 private:
+  // temporary work-around method to allow 'geometry' to be called without
+  // trying to re-calculate the non-CELL_CENTRE fields
+  int geometryNoRecalculate();
   int nz; // Size of mesh in Z. This is mesh->ngz-1
   Mesh * localmesh;
   CELL_LOC location;

include/bout/mesh.hxx

@@ -434,21 +434,26 @@ 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();
   }
 
   Coordinates *DEPRECATED(coordinates(const CELL_LOC location = CELL_CENTRE)) {
     return getCoordinates(location);
   }
 
+  /// 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);
+
   // First derivatives in index space
   // Implemented in src/mesh/index_derivs.hxx
 
@@ -750,9 +755,6 @@ class Mesh {
                            REGION region = RGN_NOBNDRY);
 
 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;

src/field/field3d.cxx

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

src/mesh/coordinates.cxx

@@ -151,7 +151,7 @@ Coordinates::Coordinates(Mesh *mesh)
 
   //////////////////////////////////////////////////////
   /// Calculate Christoffel symbols. Needs communication
-  if (geometry()) {
+  if (geometryNoRecalculate()) {
     throw BoutException("Differential geometry failed\n");
   }
 
@@ -269,7 +269,7 @@ Coordinates::Coordinates(Mesh *mesh, const CELL_LOC loc, const Coordinates* coor
 
   //////////////////////////////////////////////////////
   /// Calculate Christoffel symbols. Needs communication
-  if (geometry()) {
+  if (geometryNoRecalculate()) {
     throw BoutException("Differential geometry failed\n");
   }
 
@@ -304,7 +304,7 @@ void Coordinates::outputVars(Datafile &file) {
   file.add(J, "J", false);
 }
 
-int Coordinates::geometry() {
+int Coordinates::geometryNoRecalculate() {
   TRACE("Coordinates::geometry");
 
   output_progress.write("Calculating differential geometry terms\n");
@@ -442,13 +442,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 = localmesh->indexDDX(1. / dx); // d/di(1/dx)
   } else {
+    // Shift d2x to our location
+    d2x = interp_to(d2x, location);
+
     d1_dx = -d2x / (dx * dx);
   }
 
@@ -457,12 +460,50 @@ int Coordinates::geometry() {
         "\tWARNING: differencing quantity 'd2y' not found. Calculating from dy\n");
     d1_dy = localmesh->indexDDY(1. / dy); // d/di(1/dy)
   } else {
+    // Shift d2y to our location
+    d2y = interp_to(d2y, location);
+
     d1_dy = -d2y / (dy * dy);
   }
 
   return 0;
 }
 
+int Coordinates::geometry() {
+
+  geometryNoRecalculate();
+
+  if (location != CELL_CENTRE) {
+    throw BoutException("geometry() called from a location other than "
+      "CELL_CENTRE. This is an error as the other Coordinates are calculated "
+      "from the CELL_CENTRE version, so the changes you have made to this "
+      "object would be overwritten.");
+  }
+  // Coordinates objects at staggered location were calculated from
+  // CELL_CENTRE ones. geometry() has been called on the CELL_CENTRE
+  // Coordinates, so they must have changed; we need to re-calculate the
+  // staggered location Coordinates objects.
+
+  if (localmesh->StaggerGrids) {
+    // Replace Coordinates objects at staggered locations, if there are
+    // enough grid points.
+    // This is a temporary workaround. In v4.3 we will require users to call
+    // Mesh::addCoordinates(location) for each location that is needed and
+    // change this so that we don't waste memory on unneeded Coordinates.
+    if (localmesh->LocalNx >= 4) {
+      localmesh->addCoordinates(CELL_XLOW, true);
+    }
+    if (localmesh->LocalNy >= 4) {
+      localmesh->addCoordinates(CELL_YLOW, true);
+    }
+    // Can always add ZLOW Coordinates, since z-interpolation on Field2D is a
+    // null operation
+    localmesh->addCoordinates(CELL_ZLOW, true);
+  }
+
+  return 0;
+}
+
 int Coordinates::calcCovariant() {
   TRACE("Coordinates::calcCovariant");
 

src/mesh/impls/bout/boutmesh.cxx

@@ -839,6 +839,26 @@ int BoutMesh::load() {
   // Add boundary regions
   addBoundaryRegions();
 
+  // Create CELL_CENTRE Coordinates object
+  addCoordinates(CELL_CENTRE);
+
+  if (StaggerGrids) {
+    // Add Coordinates objects at staggered locations, if there are enough grid
+    // points.
+    // This is a temporary workaround. In v4.3 we will require users to call
+    // Mesh::addCoordinates(location) for each location that is needed and
+    // change this so that we don't waste memory on unneeded Coordinates.
+    if (LocalNx >= 4) {
+      addCoordinates(CELL_XLOW);
+    }
+    if (LocalNy >= 4) {
+      addCoordinates(CELL_YLOW);
+    }
+    // Can always add ZLOW Coordinates, since z-interpolation on Field2D is a
+    // null operation
+    addCoordinates(CELL_ZLOW);
+  }
+
   output_info.write("\tdone\n");
 
   return 0;
@@ -2208,7 +2228,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 +2245,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,38 @@ 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, std::make_shared<Coordinates>(this));
+      } else if (replace_coords) {
+        // location does already exists 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 (!coords_map.count(location)) {
+        // location does not exist in coords_map, so create new entry
+        coords_map.emplace(location, std::make_shared<Coordinates>(this, location, getCoordinates(CELL_CENTRE)));
+      } else if (replace_coords) {
+        // location does already exists in coords_map, so reset it
+        coords_map.at(location).reset(new Coordinates(this, location, getCoordinates(CELL_CENTRE)));
+      }
+    }
+  }
 }
 
 

tests/unit/test_extras.hxx

@@ -73,6 +73,11 @@ public:
     StaggerGrids = false;
     IncIntShear = false;
     maxregionblocksize = MAXREGIONBLOCKSIZE;
+
+    coords_map.emplace(CELL_CENTRE, std::shared_ptr<Coordinates>(nullptr));
+    coords_map.emplace(CELL_XLOW, std::shared_ptr<Coordinates>(nullptr));
+    coords_map.emplace(CELL_YLOW, std::shared_ptr<Coordinates>(nullptr));
+    coords_map.emplace(CELL_ZLOW, std::shared_ptr<Coordinates>(nullptr));
   }
 
   comm_handle send(FieldGroup &UNUSED(g)) { return nullptr; };