Remove mergeYupYdown()

include/bout/paralleltransform.hxx

@@ -56,10 +56,9 @@ public:
 class ParallelTransformIdentity : public ParallelTransform {
 public:
   /*!
-   * Merges the yup and ydown() fields of f, so that
-   * f.yup() = f.ydown() = f
+   * Does nothing, so field-aligned derivatives, etc. are used
    */ 
-  void calcYUpDown(Field3D &f) override {f.mergeYupYdown();}
+  void calcYUpDown(Field3D &UNUSED(f)) override {}
 
   /*!
    * The field is already aligned in Y, so this

include/field3d.hxx

@@ -152,12 +152,9 @@ class Mesh;  // #include "bout/mesh.hxx"
 
       f.yup() // error; f.yup not allocated
 
-      f.mergeYupYdown(); // f.yup() and f.ydown() now point to f
-      f.yup()(0,1,0)  // ok, gives value of f at (0,1,0)
-
       To have separate fields for yup and ydown, first call
 
-      f.splitYupYdown(); // f.yup() and f.ydown() separate
+      f.createYupYdown(); // create f.yup() and f.ydown() fields
 
       f.yup(); // ok
       f.yup()(0,1,0) // error; f.yup not allocated
@@ -227,13 +224,8 @@ class Field3D : public Field, public FieldData {
    * Ensure that this field has separate fields
    * for yup and ydown.
    */
-  void splitYupYdown();
+  void createYupYdown();
 
-  /*!
-   * Ensure that yup and ydown refer to this field
-   */
-  void mergeYupYdown();
-
   /// Check if this field has yup and ydown fields
   bool hasYupYdown() const {
     return (yup_field != nullptr) && (ydown_field != nullptr);

src/field/field3d.cxx

@@ -178,8 +178,8 @@ Field3D* Field3D::timeDeriv() {
   return deriv;
 }
 
-void Field3D::splitYupYdown() {
-  TRACE("Field3D::splitYupYdown");
+void Field3D::createYupYdown() {
+  TRACE("Field3D::createYupYdown");
 
   if((yup_field != this) && (yup_field != nullptr))
     return;
@@ -189,24 +189,6 @@ void Field3D::splitYupYdown() {
   ydown_field = new Field3D(fieldmesh);
 }
 
-void Field3D::mergeYupYdown() {
-  TRACE("Field3D::mergeYupYdown");
-
-  if(yup_field == this && ydown_field == this)
-    return;
-
-  if(yup_field != nullptr){
-    delete yup_field;
-  }
-
-  if(ydown_field != nullptr) {
-    delete ydown_field;
-  }
-
-  yup_field = this;
-  ydown_field = this;
-}
-
 Field3D& Field3D::ynext(int dir) {
   switch(dir) {
   case +1:

src/mesh/coordinates.cxx

@@ -680,15 +680,9 @@ const Field3D Coordinates::Div_par(const Field3D &f, CELL_LOC outloc,
   if (f.hasYupYdown()) {
     // Need to modify yup and ydown fields
     Field3D f_B = f / Bxy;
-    if (&f.yup() == &f) {
-      // Identity, yup and ydown point to same field
-      f_B.mergeYupYdown();
-    } else {
-      // Distinct fields
-      f_B.splitYupYdown();
-      f_B.yup() = f.yup() / Bxy;
-      f_B.ydown() = f.ydown() / Bxy;
-    }
+    f_B.createYupYdown();
+    f_B.yup() = f.yup() / Bxy;
+    f_B.ydown() = f.ydown() / Bxy;
     return Bxy * Grad_par(f_B, outloc, method);
   }
 

src/mesh/difops.cxx

@@ -315,8 +315,8 @@ const Field3D Vpar_Grad_par_LCtoC(const Field3D &v, const Field3D &f, REGION reg
 
   result.allocate();
 
-  bool vUseUpDown = (v.hasYupYdown() && ((&v.yup() != &v) || (&v.ydown() != &v)));
-  bool fUseUpDown = (f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f)));
+  bool vUseUpDown = v.hasYupYdown();
+  bool fUseUpDown = f.hasYupYdown();
 
   if (vUseUpDown && fUseUpDown) {
     // Both v and f have up/down fields

src/mesh/index_derivs.cxx

@@ -928,7 +928,7 @@ const Field3D Mesh::applyYdiff(const Field3D &var, Mesh::deriv_func func, CELL_L
   Field3D result(this);
   result.allocate(); // Make sure data allocated
 
-  if (var.hasYupYdown() && ((&var.yup() != &var) || (&var.ydown() != &var))) {
+  if (var.hasYupYdown()) {
     // Field "var" has distinct yup and ydown fields which
     // will be used to calculate a derivative along
     // the magnetic field
@@ -2230,8 +2230,8 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
     // If vUseUpDown is true, field "v" has distinct yup and ydown fields which
     // will be used to calculate a derivative along
     // the magnetic field
-    bool vUseUpDown = (v.hasYupYdown() && ((&v.yup() != &v) || (&v.ydown() != &v)));
-    bool fUseUpDown = (f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f)));
+    bool vUseUpDown = v.hasYupYdown();
+    bool fUseUpDown = f.hasYupYdown();
 
     if (vUseUpDown && fUseUpDown) {
       // Both v and f have up/down fields
@@ -2379,7 +2379,7 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
       func = lookupFunc(table, method);
     }
 
-    if (f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f))) {
+    if (f.hasYupYdown()) {
       // f has yup and ydown fields which are distinct
 
       stencil fs;
@@ -2993,8 +2993,8 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo
   // If vUseUpDown is true, field "v" has distinct yup and ydown fields which
   // will be used to calculate a derivative along
   // the magnetic field
-  bool vUseUpDown = (v.hasYupYdown() && ((&v.yup() != &v) || (&v.ydown() != &v)));
-  bool fUseUpDown = (f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f)));
+  bool vUseUpDown = v.hasYupYdown();
+  bool fUseUpDown = f.hasYupYdown();
 
   if (vUseUpDown && fUseUpDown) {
     // Both v and f have up/down fields

src/mesh/interpolation.cxx

@@ -116,7 +116,7 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) {
       case CELL_YLOW: {
         ASSERT0(fieldmesh->ystart >= 2); // At least 2 boundary cells needed for interpolation in y-direction
 
-        if (var.hasYupYdown() && ((&var.yup() != &var) || (&var.ydown() != &var))) {
+        if (var.hasYupYdown()) {
           // Field "var" has distinct yup and ydown fields which
           // will be used to calculate a derivative along
           // the magnetic field

src/mesh/parallel/fci.cxx

@@ -296,8 +296,8 @@ const Field3D FCIMap::integrate(Field3D &f) const {
 void FCITransform::calcYUpDown(Field3D &f) {
   TRACE("FCITransform::calcYUpDown");
 
-  // Ensure that yup and ydown are different fields
-  f.splitYupYdown();
+  // Ensure that yup and ydown are not null
+  f.createYupYdown();
 
   // Interpolate f onto yup and ydown fields
   f.ynext(forward_map.dir) = forward_map.interpolate(f);
@@ -307,8 +307,8 @@ void FCITransform::calcYUpDown(Field3D &f) {
 void FCITransform::integrateYUpDown(Field3D &f) {
   TRACE("FCITransform::integrateYUpDown");
 
-  // Ensure that yup and ydown are different fields
-  f.splitYupYdown();
+  // Ensure that yup and ydown are not null
+  f.createYupYdown();
 
   // Integrate f onto yup and ydown fields
   f.ynext(forward_map.dir) = forward_map.integrate(f);

src/mesh/parallel/shiftedmetric.cxx

@@ -95,7 +95,8 @@ ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) {
  * Calculate the Y up and down fields
  */
 void ShiftedMetric::calcYUpDown(Field3D &f) {
-  f.splitYupYdown();
+  // Ensure that yup and ydown are not null
+  f.createYupYdown();
 
   Field3D& yup = f.yup();
   yup.allocate();

src/physics/smoothing.cxx

@@ -74,6 +74,7 @@ const Field3D smooth_y(const Field3D &f) {
   Mesh *mesh = f.getMesh();
   Field3D result(mesh);
   result.allocate();
+  result.setLocation(f.getLocation());
 
   // Copy boundary region
   for(int jx=0;jx<mesh->LocalNx;jx++)
@@ -84,11 +85,21 @@ const Field3D smooth_y(const Field3D &f) {
 
   // Smooth using simple 1-2-1 filter
 
-  for(int jx=0;jx<mesh->LocalNx;jx++)
-    for(int jy=1;jy<mesh->LocalNy-1;jy++)
-      for(int jz=0;jz<mesh->LocalNz;jz++) {
-	result(jx,jy,jz) = 0.5*f(jx,jy,jz) + 0.25*( f.ydown()(jx,jy-1,jz) + f.yup()(jx,jy+1,jz) );
-      }
+  if (f.hasYupYdown()) {
+    for(int jx=0;jx<mesh->LocalNx;jx++)
+      for(int jy=1;jy<mesh->LocalNy-1;jy++)
+        for(int jz=0;jz<mesh->LocalNz;jz++) {
+          result(jx,jy,jz) = 0.5*f(jx,jy,jz) + 0.25*( f.ydown()(jx,jy-1,jz) + f.yup()(jx,jy+1,jz) );
+        }
+  } else {
+    Field3D f_fa = mesh->toFieldAligned(f);
+    for(int jx=0;jx<mesh->LocalNx;jx++)
+      for(int jy=1;jy<mesh->LocalNy-1;jy++)
+        for(int jz=0;jz<mesh->LocalNz;jz++) {
+          result(jx,jy,jz) = 0.5*f_fa(jx,jy,jz) + 0.25*( f_fa(jx,jy-1,jz) + f_fa(jx,jy+1,jz) );
+        }
+    result = mesh->fromFieldAligned(result);
+  }
 
   // Need to communicate boundaries
   mesh->communicate(result);

tests/integrated/test-smooth/test_smooth.cxx

@@ -17,8 +17,6 @@ int main(int argc, char **argv) {
   Field2D input2d = f.create2D("1 + sin(2*y)");
   Field3D input3d = f.create3D("gauss(x-0.5,0.2)*gauss(y-pi)*sin(3*y - z)");
 
-  input3d.mergeYupYdown();
-
   SAVE_ONCE2(input2d, input3d);
 
   // Average in 3D

tests/unit/field/test_field3d.cxx

@@ -213,7 +213,7 @@ TEST_F(Field3DTest, SplitYupYDown) {
 
   EXPECT_FALSE(field.hasYupYdown());
 
-  field.splitYupYdown();
+  field.createYupYdown();
 
   EXPECT_TRUE(field.hasYupYdown());
 
@@ -223,7 +223,7 @@ TEST_F(Field3DTest, SplitYupYDown) {
   EXPECT_NE(&field, &ydown);
 
   // Should be able to split again without any problems
-  field.splitYupYdown();
+  field.createYupYdown();
 
   // Would be nice to check yup2 != yup, but not sure this is possible
   // to do in general
@@ -233,55 +233,11 @@ TEST_F(Field3DTest, SplitYupYDown) {
   EXPECT_NE(&field, &ydown2);
 }
 
-TEST_F(Field3DTest, MergeYupYDown) {
-  Field3D field;
-
-  field = 0.;
-
-  EXPECT_FALSE(field.hasYupYdown());
-
-  field.mergeYupYdown();
-
-  EXPECT_TRUE(field.hasYupYdown());
-
-  auto& yup = field.yup();
-  EXPECT_EQ(&field, &yup);
-  auto& ydown = field.ydown();
-  EXPECT_EQ(&field, &ydown);
-
-  // Should be able to merge again without any problems
-  field.mergeYupYdown();
-
-  auto& yup2 = field.yup();
-  EXPECT_EQ(&field, &yup2);
-  auto& ydown2 = field.ydown();
-  EXPECT_EQ(&field, &ydown2);
-}
-
-TEST_F(Field3DTest, SplitThenMergeYupYDown) {
-  Field3D field;
-
-  field = 0.;
-  field.splitYupYdown();
-
-  auto& yup = field.yup();
-  EXPECT_NE(&field, &yup);
-  auto& ydown = field.ydown();
-  EXPECT_NE(&field, &ydown);
-
-  field.mergeYupYdown();
-
-  auto& yup2 = field.yup();
-  EXPECT_EQ(&field, &yup2);
-  auto& ydown2 = field.ydown();
-  EXPECT_EQ(&field, &ydown2);
-}
-
 TEST_F(Field3DTest, Ynext) {
   Field3D field;
 
   field = 0.;
-  field.splitYupYdown();
+  field.createYupYdown();
 
   auto& yup = field.ynext(1);
   EXPECT_NE(&field, &yup);
@@ -295,7 +251,7 @@ TEST_F(Field3DTest, Ynext) {
 TEST_F(Field3DTest, ConstYnext) {
   Field3D field(0.);
 
-  field.splitYupYdown();
+  field.createYupYdown();
 
   const Field3D& field2 = field;