New shifted metric 2: coordinate system labelling

include/bout/index_derivs_interface.hxx

@@ -96,6 +96,9 @@ T flowDerivative(const T& vel, const T& f, CELL_LOC outloc, const std::string& m
   T result(localmesh);
   result.allocate(); // Make sure data allocated
   result.setLocation(outloc);
+  if (std::is_base_of<Field3D, T>::value) {
+    result.setCoordinateSystem(f.getCoordinateSystem());
+  }
 
   // Apply method
   derivativeMethod(vel, f, result, region);
@@ -162,6 +165,9 @@ T standardDerivative(const T& f, CELL_LOC outloc, const std::string& method,
   T result(localmesh);
   result.allocate(); // Make sure data allocated
   result.setLocation(outloc);
+  if (std::is_base_of<Field3D, T>::value) {
+    result.setCoordinateSystem(f.getCoordinateSystem());
+  }
 
   // Apply method
   derivativeMethod(f, result, region);

include/bout/mesh.hxx

@@ -689,22 +689,32 @@ class Mesh {
 
   /// Transform a field into field-aligned coordinates
   const Field3D toFieldAligned(const Field3D &f, const REGION region = RGN_NOX) {
-    return getParallelTransform().toFieldAligned(f, region);
+    return transform->toFieldAligned(f, region);
   }
   const Field2D toFieldAligned(const Field2D &f, const REGION UNUSED(region) = RGN_NOX) {
     return f;
   }
 
   /// Convert back into standard form
   const Field3D fromFieldAligned(const Field3D &f, const REGION region = RGN_NOX) {
-    return getParallelTransform().fromFieldAligned(f, region);
+    return transform->fromFieldAligned(f, region);
   }
   const Field2D fromFieldAligned(const Field2D &f, const REGION UNUSED(region) = RGN_NOX) {
     return f;
   }
 
+  COORDINATE_SYSTEM getCoordinateSystem() const {
+    if (transform) {
+      return transform->getCoordinateSystem();
+    } else {
+      // ParallelTransform not initialized yet. Any Field3D created better not
+      // need to know its coordinate system.
+      return COORDINATE_SYSTEM::None;
+    }
+  }
+
   bool canToFromFieldAligned() {
-    return getParallelTransform().canToFromFieldAligned();
+    return transform->canToFromFieldAligned();
   }
 
   /*!
@@ -725,7 +735,7 @@ class Mesh {
   void setParallelTransform();
 
   /*!
-   * Return the parallel transform, setting it if need be
+   * Return the parallel transform
    */
   ParallelTransform& getParallelTransform();
 

include/bout/paralleltransform.hxx

@@ -43,6 +43,11 @@ public:
   /// into standard form
   virtual const Field3D fromFieldAligned(const Field3D &f, const REGION region = RGN_NOX) = 0;
 
+  /*!
+   * Return the default coordinate system for Field3Ds when using this ParallelTransform
+   */
+  virtual COORDINATE_SYSTEM getCoordinateSystem() const = 0;
+
   virtual bool canToFromFieldAligned() = 0;
 };
 
@@ -76,6 +81,10 @@ public:
     return f;
   }
 
+  COORDINATE_SYSTEM getCoordinateSystem() const override {
+    return COORDINATE_SYSTEM::FieldAligned;
+  }
+
   bool canToFromFieldAligned() override{
     return true;
   }
@@ -116,6 +125,10 @@ public:
    */
   const Field3D fromFieldAligned(const Field3D &f, const REGION region=RGN_NOX) override;
 
+  COORDINATE_SYSTEM getCoordinateSystem() const override {
+    return COORDINATE_SYSTEM::Orthogonal;
+  }
+
   bool canToFromFieldAligned() override{
     return true;
   }

include/field.hxx

@@ -46,6 +46,28 @@ extern Mesh * mesh; ///< Global mesh
 #include <string>
 #endif
 
+/*!
+ * Coordinate systems that Field3Ds can be stored in
+ *
+ * None - special value to be used for Field3Ds initialized before the Mesh and
+ * ParallelTransform have finished being initialized, because then we cannot
+ * get the COORDINATE_SYSTEM to use from the ParallelTransform, but Field3Ds
+ * created then must not need to know their COORDINATE_SYSTEM, or must set it
+ * explicitly later.
+ */
+enum class COORDINATE_SYSTEM { None, FieldAligned, Orthogonal, FCI };
+
+#define COORDENUMSTR(val) {COORDINATE_SYSTEM::val, #val}
+
+const std::map<COORDINATE_SYSTEM, std::string> COORDINATE_SYSTEMtoString = {
+  COORDENUMSTR(None),
+  COORDENUMSTR(FieldAligned),
+  COORDENUMSTR(Orthogonal),
+  COORDENUMSTR(FCI)
+};
+
+#undef COORDENUMSTR
+
 /*!
  * \brief Base class for fields
  *

include/field2d.hxx

@@ -111,6 +111,24 @@ class Field2D : public Field, public FieldData {
    */
   int getNz() const override {return 1;};
 
+  /*!
+   * Return the coordinate system of this field
+   */
+  COORDINATE_SYSTEM getCoordinateSystem() const {
+#if CHECK > 0
+    throw BoutException("getCoordinateSystem() should not be called for Field2D");
+#endif
+  }
+
+  /*!
+   * Reset the coordinate system of this field
+   */
+  void setCoordinateSystem(COORDINATE_SYSTEM UNUSED(new_coords)) {
+#if CHECK > 0
+    throw BoutException("setCoordinateSystem() should not be called for Field2D");
+#endif
+  }
+
   /// Check if this field has yup and ydown fields
   bool hasYupYdown() const {
     return true;

include/field3d.hxx

@@ -215,6 +215,16 @@ class Field3D : public Field, public FieldData {
    */
   int getNz() const override {return nz;};
 
+  /*!
+   * Return the coordinate system of this field
+   */
+  COORDINATE_SYSTEM getCoordinateSystem() const { return coordinate_system; }
+
+  /*!
+   * Reset the coordinate system of this field
+   */
+  void setCoordinateSystem(COORDINATE_SYSTEM new_coords) { coordinate_system = new_coords; }
+
   /*!
    * Ensure that this field has separate fields
    * for yup and ydown.
@@ -460,6 +470,7 @@ class Field3D : public Field, public FieldData {
     swap(first.nx, second.nx);
     swap(first.ny, second.ny);
     swap(first.nz, second.nz);
+    swap(first.coordinate_system, second.coordinate_system);
     swap(first.location, second.location);
     swap(first.deriv, second.deriv);
     swap(first.yup_field, second.yup_field);
@@ -478,6 +489,8 @@ private:
   /// Array sizes (from fieldmesh). These are valid only if fieldmesh is not null
   int nx{-1}, ny{-1}, nz{-1};
 
+  COORDINATE_SYSTEM coordinate_system; ///< Coordinate system used by this field, e.g. fieldaligned, orthogonal
+
   /// Internal data array. Handles allocation/freeing of memory
   Array<BoutReal> data;
 

src/field/field3d.cxx

@@ -54,11 +54,19 @@ Field3D::Field3D(Mesh* localmesh) : Field(localmesh) {
     ny = fieldmesh->LocalNy;
     nz = fieldmesh->LocalNz;
   }
+  if (fieldmesh) {
+    coordinate_system = fieldmesh->getCoordinateSystem();
+  } else {
+    coordinate_system = COORDINATE_SYSTEM::None;
+  }
+
 }
 
 /// Doesn't copy any data, just create a new reference to the same data (copy on change
 /// later)
-Field3D::Field3D(const Field3D& f) : Field(f.fieldmesh), data(f.data) {
+Field3D::Field3D(const Field3D& f)
+  : Field(f.fieldmesh), coordinate_system(f.coordinate_system), data(f.data),
+    location(f.location) {
 
   TRACE("Field3D(Field3D&)");
 
@@ -72,7 +80,6 @@ Field3D::Field3D(const Field3D& f) : Field(f.fieldmesh), data(f.data) {
     nz = fieldmesh->LocalNz;
   }
 
-  location = f.location;
   fieldCoordinates = f.fieldCoordinates;
 }
 
@@ -86,6 +93,7 @@ Field3D::Field3D(const Field2D& f) : Field(f.getMesh()) {
 
   location = f.getLocation();
   fieldCoordinates = nullptr;
+  coordinate_system = fieldmesh->getCoordinateSystem();
 
   *this = f;
 }
@@ -98,6 +106,8 @@ Field3D::Field3D(const BoutReal val, Mesh* localmesh) : Field(localmesh) {
   ny = fieldmesh->LocalNy;
   nz = fieldmesh->LocalNz;
 
+  coordinate_system = fieldmesh->getCoordinateSystem();
+
   *this = val;
 }
 
@@ -275,6 +285,8 @@ Field3D & Field3D::operator=(const Field3D &rhs) {
   fieldmesh = rhs.fieldmesh;
   nx = rhs.nx; ny = rhs.ny; nz = rhs.nz; 
 
+  coordinate_system = rhs.coordinate_system;
+
   data = rhs.data;
 
   setLocation(rhs.location);
@@ -614,6 +626,7 @@ Field3D pow(const Field3D &lhs, const Field3D &rhs, REGION rgn) {
   TRACE("pow(Field3D, Field3D)");
 
   ASSERT1(lhs.getLocation() == rhs.getLocation());
+  ASSERT1(lhs.getCoordinateSystem() == rhs.getCoordinateSystem());
 
   ASSERT1(lhs.getMesh() == rhs.getMesh());
   Field3D result(lhs.getMesh());
@@ -622,6 +635,7 @@ Field3D pow(const Field3D &lhs, const Field3D &rhs, REGION rgn) {
   BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs[i]); }
 
   result.setLocation( lhs.getLocation() );
+  result.setCoordinateSystem( lhs.getCoordinateSystem() );
 
   checkData(result);
   return result;
@@ -641,6 +655,7 @@ Field3D pow(const Field3D &lhs, const Field2D &rhs, REGION rgn) {
   BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs[i]); }
 
   result.setLocation( lhs.getLocation() );
+  result.setCoordinateSystem( lhs.getCoordinateSystem() );
 
   checkData(result);
   return result;
@@ -653,6 +668,11 @@ FieldPerp pow(const Field3D &lhs, const FieldPerp &rhs, REGION rgn) {
   checkData(rhs);
   ASSERT1(lhs.getMesh() == rhs.getMesh());
 
+  // Assume FieldPerp is always on the default coordinate system of its Mesh,
+  // since at the moment FieldPerp cannot switch between orthogonal and field
+  // aligned coordinates
+  ASSERT2(lhs.getCoordinateSystem() == rhs.getMesh()->getCoordinateSystem());
+
   FieldPerp result{rhs.getMesh()};
   result.allocate();
   result.setIndex(rhs.getIndex());
@@ -679,6 +699,7 @@ Field3D pow(const Field3D &lhs, BoutReal rhs, REGION rgn) {
   BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs); }
 
   result.setLocation( lhs.getLocation() );
+  result.setCoordinateSystem( lhs.getCoordinateSystem() );
 
   checkData(result);
   return result;
@@ -697,6 +718,7 @@ Field3D pow(BoutReal lhs, const Field3D &rhs, REGION rgn) {
   BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs, rhs[i]); }
 
   result.setLocation( rhs.getLocation() );
+  result.setCoordinateSystem( rhs.getCoordinateSystem() );
 
   checkData(result);
   return result;
@@ -802,6 +824,7 @@ BoutReal mean(const Field3D &f, bool allpe, REGION rgn) {
     result.allocate();                                                                   \
     BOUT_FOR(d, result.getRegion(rgn)) { result[d] = func(f[d]); }                       \
     result.setLocation(f.getLocation());                                                 \
+    result.setCoordinateSystem(f.getCoordinateSystem());                                 \
     checkData(result);                                                                   \
     return result;                                                                       \
   }
@@ -865,6 +888,7 @@ const Field3D filter(const Field3D &var, int N0, REGION rgn) {
 #endif
 
   result.setLocation(var.getLocation());
+  result.setCoordinateSystem(var.getCoordinateSystem());
 
   checkData(result);
   return result;
@@ -912,6 +936,7 @@ const Field3D lowPass(const Field3D &var, int zmax, REGION rgn) {
     }
   }
   result.setLocation(var.getLocation());
+  result.setCoordinateSystem(var.getCoordinateSystem());
 
   checkData(result);
   return result;
@@ -962,6 +987,7 @@ const Field3D lowPass(const Field3D &var, int zmax, int zmin, REGION rgn) {
   }
 
   result.setLocation(var.getLocation());
+  result.setCoordinateSystem(var.getCoordinateSystem());
 
   checkData(result);
   return result;

src/field/field_factory.cxx

@@ -229,6 +229,7 @@ const Field3D FieldFactory::create3D(const std::string &value, const Options *op
   if (localmesh->canToFromFieldAligned()){ // Ask wheter it is possible
     // Transform from field aligned coordinates, to be compatible with
     // older BOUT++ inputs. This is not a particularly "nice" solution.
+    result.setCoordinateSystem(COORDINATE_SYSTEM::FieldAligned);
     result = localmesh->fromFieldAligned(result, RGN_ALL);
   }
 

src/field/gen_fieldops.jinja

@@ -16,6 +16,10 @@
     ASSERT1(localmesh == {{rhs.name}}.getMesh());
   {% endif %}
 
+  {% if lhs == "Field3D" and rhs == "Field3D" %}
+    ASSERT1({{lhs.name}}.getCoordinateSystem() == {{rhs.name}}.getCoordinateSystem());
+  {% endif %}
+
   {{out.field_type}} {{out.name}}(localmesh);
   {{out.name}}.allocate();
   checkData({{lhs.name}});
@@ -46,8 +50,12 @@
 
   {% if out in ['Field3D','Field2D'] %}
     {{out.name}}.setLocation({{rhs.name if rhs in ["Field3D","Field2D"] else lhs.name}}.getLocation());
+
   {% endif %}
+  {% if out in ['Field3D'] %}
+    {{out.name}}.setCoordinateSystem({{rhs.name if rhs == "Field3D" else lhs.name}}.getCoordinateSystem());
 
+  {% endif %}
   checkData({{out.name}});
   return {{out.name}};
 }
@@ -58,7 +66,11 @@
   // only if data is unique we update the field
   // otherwise just call the non-inplace version
   if (data.unique()) {
-    {% if ((lhs in ["Field3D", "Field2D"]) and (rhs in ["Field3D", "Field2D"])) %}    
+    {% if ((lhs in ["Field3D"]) and (rhs in ["Field3D"])) %}
+    ASSERT1(this->getCoordinateSystem() == {{rhs.name}}.getCoordinateSystem())
+
+    {% endif %}
+    {% if ((lhs in ["Field3D", "Field2D"]) and (rhs in ["Field3D", "Field2D"])) %}
 #if CHECK > 0
     if (this->getLocation() != rhs.getLocation()) {
       throw BoutException(