Shiftedmetric updates

include/bout/mesh.hxx

@@ -609,13 +609,19 @@ class Mesh {
   /// Derivative functions of a velocity field, and field stencil v, f
   typedef BoutReal (*flux_func)(stencil&, stencil &);
 
+  /// Calculate yup/ydown fields in case they can be computed for an
+  /// intermediate variable without communicating
+  void calcYUpDown(Field3D &f) {
+    getParallelTransform().calcYUpDown(f);
+  }
+
   /// Transform a field into field-aligned coordinates
-  const Field3D toFieldAligned(const Field3D &f) {
-    return getParallelTransform().toFieldAligned(f);
+  const Field3D toFieldAligned(const Field3D &f, const REGION region = RGN_NOX) {
+    return getParallelTransform().toFieldAligned(f, region);
   }
   /// Convert back into standard form
-  const Field3D fromFieldAligned(const Field3D &f) {
-    return getParallelTransform().fromFieldAligned(f);
+  const Field3D fromFieldAligned(const Field3D &f, const REGION region = RGN_NOX) {
+    return getParallelTransform().fromFieldAligned(f, region);
   }
 
   bool canToFromFieldAligned() {

include/bout/paralleltransform.hxx

@@ -38,11 +38,11 @@ public:
 
   /// Convert a 3D field into field-aligned coordinates
   /// so that the y index is along the magnetic field
-  virtual const Field3D toFieldAligned(const Field3D &f) = 0;
+  virtual const Field3D toFieldAligned(const Field3D &f, const REGION region = RGN_NOX) = 0;
 
   /// Convert back from field-aligned coordinates
   /// into standard form
-  virtual const Field3D fromFieldAligned(const Field3D &f) = 0;
+  virtual const Field3D fromFieldAligned(const Field3D &f, const REGION region = RGN_NOX) = 0;
 
   virtual bool canToFromFieldAligned() = 0;
 };
@@ -59,21 +59,24 @@ public:
    * Merges the yup and ydown() fields of f, so that
    * f.yup() = f.ydown() = f
    */ 
-  void calcYUpDown(Field3D &f) override {f.mergeYupYdown();}
+  void calcYUpDown(Field3D &f) override {
+    f.mergeYupYdown();
+    f.setHasValidYUpDown(true);
+  }
 
   /*!
    * The field is already aligned in Y, so this
    * does nothing
    */ 
-  const Field3D toFieldAligned(const Field3D &f) override {
+  const Field3D toFieldAligned(const Field3D &f, const REGION UNUSED(region)) override {
     return f;
   }
 
   /*!
    * The field is already aligned in Y, so this
    * does nothing
    */
-  const Field3D fromFieldAligned(const Field3D &f) override {
+  const Field3D fromFieldAligned(const Field3D &f, const REGION UNUSED(region)) override {
     return f;
   }
 
@@ -109,13 +112,13 @@ public:
    * in X-Z, and the metric tensor will need to be changed 
    * if X derivatives are used.
    */
-  const Field3D toFieldAligned(const Field3D &f) override;
+  const Field3D toFieldAligned(const Field3D &f, const REGION region=RGN_NOX) override;
 
   /*!
    * Converts a field back to X-Z orthogonal coordinates
    * from field aligned coordinates.
    */
-  const Field3D fromFieldAligned(const Field3D &f) override;
+  const Field3D fromFieldAligned(const Field3D &f, const REGION region=RGN_NOX) override;
 
   bool canToFromFieldAligned() override{
     return true;
@@ -142,7 +145,7 @@ private:
    * Shift a 2D field in Z. 
    * Since 2D fields are constant in Z, this has no effect
    */
-  const Field2D shiftZ(const Field2D &f, const Field2D &UNUSED(zangle)){return f;};
+  const Field2D shiftZ(const Field2D &f, const Field2D &UNUSED(zangle), const REGION UNUSED(region)=RGN_NOX){return f;};
 
   /*!
    * Shift a 3D field \p f in Z by the given \p zangle
@@ -151,7 +154,7 @@ private:
    * @param[in] zangle   Toroidal angle (z)
    *
    */ 
-  const Field3D shiftZ(const Field3D &f, const Field2D &zangle);
+  const Field3D shiftZ(const Field3D &f, const Field2D &zangle, const REGION region=RGN_NOX);
 
   /*!
    * Shift a 3D field \p f by the given phase \p phs in Z
@@ -162,7 +165,7 @@ private:
    * @param[in] f  The field to shift
    * @param[in] phs  The phase to shift by
    */
-  const Field3D shiftZ(const Field3D &f, const arr3Dvec &phs);
+  const Field3D shiftZ(const Field3D &f, const arr3Dvec &phs, const REGION region=RGN_NOX);
 
   /*!
    * Shift a given 1D array, assumed to be in Z, by the given \p zangle

include/field3d.hxx

@@ -223,10 +223,23 @@ class Field3D : public Field, public FieldData {
    * Ensure that yup and ydown refer to this field
    */
   void mergeYupYdown();
+
+  /*!
+   * Delete all yup/ydown fields and field_fa
+   */
+  void deleteYupYdown();
 
   /// Check if this field has yup and ydown fields
   bool hasYupYdown() const {
-    return (yup_field != nullptr) && (ydown_field != nullptr);
+    return hasValidYUpDown && (yup_field != nullptr) && (ydown_field != nullptr);
+  }
+
+  /// Set validity state of yup/ydown fields
+  void setHasValidYUpDown(bool set) {
+    if (set) {
+      ASSERT1(yup_field != nullptr && ydown_field != nullptr);
+    }
+    hasValidYUpDown = set;
   }
 
   /// Return reference to yup field
@@ -457,6 +470,9 @@ private:
 
   Field3D *deriv; ///< Time derivative (may be NULL)
 
+  /// flag to record whether yup_field and ydown_field have been set correctly
+  bool hasValidYUpDown;
+
   /// Pointers to fields containing values along Y
   Field3D *yup_field, *ydown_field;
 };

src/field/field3d.cxx

@@ -146,12 +146,8 @@ Field3D::~Field3D() {
     // Now delete them as part of the deriv vector
     delete deriv;
   }
-
-  if((yup_field != this) && (yup_field != nullptr))
-    delete yup_field;
-
-  if((ydown_field != this) && (ydown_field != nullptr))
-    delete ydown_field;
+
+  deleteYupYdown();
 }
 
 void Field3D::allocate() {
@@ -195,18 +191,25 @@ void Field3D::mergeYupYdown() {
   if(yup_field == this && ydown_field == this)
     return;
 
-  if(yup_field != nullptr){
-    delete yup_field;
-  }
-
-  if(ydown_field != nullptr) {
-    delete ydown_field;
-  }
+  deleteYupYdown();
 
   yup_field = this;
   ydown_field = this;
 }
 
+void Field3D::deleteYupYdown() {
+  // Delete auxiliary fields if they have been set
+  if (yup_field == this && ydown_field == this) {
+    return;
+  }
+
+  delete yup_field;
+  yup_field = nullptr;
+
+  delete ydown_field;
+  ydown_field = nullptr;
+}
+
 Field3D& Field3D::ynext(int dir) {
   switch(dir) {
   case +1:
@@ -307,6 +310,7 @@ Field3D & Field3D::operator=(const Field3D &rhs) {
 
   setLocation(rhs.location);
 
+  setHasValidYUpDown(false);
   return *this;
 }
 
@@ -328,7 +332,9 @@ Field3D & Field3D::operator=(const Field2D &rhs) {
 
   /// Only 3D fields have locations for now
   //location = CELL_CENTRE;
-
+
+  setHasValidYUpDown(false);
+
   return *this;
 }
 
@@ -347,6 +353,8 @@ void Field3D::operator=(const FieldPerp &rhs) {
   BOUT_FOR(i, region_all) {
     (*this)(i, rhs.getIndex()) = rhs[i];
   }
+
+  setHasValidYUpDown(false);
 }
 
 Field3D & Field3D::operator=(const BoutReal val) {
@@ -366,6 +374,8 @@ Field3D & Field3D::operator=(const BoutReal val) {
   //location = CELL_CENTRE;
   // DON'T RE-SET LOCATION
 
+  setHasValidYUpDown(false);
+
   return *this;
 }
 

src/field/gen_fieldops.jinja

@@ -95,6 +95,11 @@
   } else {
     (*this) = (*this) {{operator}} {{rhs.name}};
   }
+
+  {% if (out == "Field3D") %}
+    setHasValidYUpDown(false);
+  {% endif %}
+
   return *this;
 }
 {% endif %}