diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx index a25540354a..0c0a8c7efa 100644 --- a/include/bout/paralleltransform.hxx +++ b/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 diff --git a/include/field3d.hxx b/include/field3d.hxx index 2ae253086d..9cc526a058 100644 --- a/include/field3d.hxx +++ b/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); diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 876dc925bf..121037ae74 100644 --- a/src/field/field3d.cxx +++ b/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: diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 29f27c13bf..d259dd7d0a 100644 --- a/src/mesh/coordinates.cxx +++ b/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); } diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 859b6cb5d7..cdd658443a 100644 --- a/src/mesh/difops.cxx +++ b/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 diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx index e3cf4af6ce..7d6a997610 100644 --- a/src/mesh/index_derivs.cxx +++ b/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 diff --git a/src/mesh/interpolation.cxx b/src/mesh/interpolation.cxx index e5d5e6da97..00231d3188 100644 --- a/src/mesh/interpolation.cxx +++ b/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 diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 63f53856a7..e7cf67ee7a 100644 --- a/src/mesh/parallel/fci.cxx +++ b/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); diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx index 13067fe07c..cfedf9150a 100644 --- a/src/mesh/parallel/shiftedmetric.cxx +++ b/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(); diff --git a/src/physics/smoothing.cxx b/src/physics/smoothing.cxx index 7386a51037..bf03217773 100644 --- a/src/physics/smoothing.cxx +++ b/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;jxLocalNx;jx++) @@ -84,11 +85,21 @@ const Field3D smooth_y(const Field3D &f) { // Smooth using simple 1-2-1 filter - for(int jx=0;jxLocalNx;jx++) - for(int jy=1;jyLocalNy-1;jy++) - for(int jz=0;jzLocalNz;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;jxLocalNx;jx++) + for(int jy=1;jyLocalNy-1;jy++) + for(int jz=0;jzLocalNz;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;jxLocalNx;jx++) + for(int jy=1;jyLocalNy-1;jy++) + for(int jz=0;jzLocalNz;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); diff --git a/tests/integrated/test-smooth/test_smooth.cxx b/tests/integrated/test-smooth/test_smooth.cxx index f1246d0fed..d0884d2734 100644 --- a/tests/integrated/test-smooth/test_smooth.cxx +++ b/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 diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index 71f39f34a6..135d7fa3e8 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/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;