diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index 8ae0e90419..74a2e21371 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -614,12 +614,12 @@ class Mesh { typedef BoutReal (*flux_func)(stencil&, stencil &); /// 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() { diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx index afd296844d..03d280454b 100644 --- a/include/bout/paralleltransform.hxx +++ b/include/bout/paralleltransform.hxx @@ -6,10 +6,13 @@ #ifndef __PARALLELTRANSFORM_H__ #define __PARALLELTRANSFORM_H__ +#include #include +#include #include #include #include +#include class Mesh; @@ -38,13 +41,16 @@ 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; + + /// Write out ParallelTransform variables to file + virtual void outputVars(Datafile &UNUSED(file)) {}; }; @@ -65,7 +71,7 @@ public: * 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; } @@ -73,13 +79,16 @@ public: * 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; } bool canToFromFieldAligned() override{ return true; } + + /// Write out ParallelTransform variables to file + virtual void outputVars(Datafile &UNUSED(file)) {}; }; /*! @@ -108,20 +117,18 @@ 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; } - /// A 3D array, implemented as nested vectors - typedef std::vector>> arr3Dvec; private: ShiftedMetric(); @@ -129,21 +136,49 @@ private: /// This is the shift in toroidal angle (z) which takes a point from /// X-Z orthogonal to field-aligned along Y. - Field2D zShift; - std::vector cmplx; ///< A temporary array, used for input/output to fft routines - std::vector cmplxLoc; ///< A temporary array, used for input/output to fft routines - - arr3Dvec toAlignedPhs; ///< Cache of phase shifts for transforming from X-Z orthogonal coordinates to field-aligned coordinates - arr3Dvec fromAlignedPhs; ///< Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates - - arr3Dvec yupPhs; ///< Cache of phase shifts for calculating yup fields - arr3Dvec ydownPhs; ///< Cache of phase shifts for calculating ydown fields + Flexible zShift; + Array cmplx; ///< A temporary array, used for input/output to fft routines + Array cmplxLoc; ///< A temporary array, used for input/output to fft routines + + Matrix< Array > getToAlignedPhs(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary; + Matrix< Array > getFromAlignedPhs(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary; + Matrix< Array > getYupPhs1(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary; + Matrix< Array > getYdownPhs1(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary; + Matrix< Array > getYupPhs2(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary; + Matrix< Array > getYdownPhs2(CELL_LOC location = CELL_CENTRE); ///< Get phase shifts, calculating if necessary; + + bool has_toAligned_CENTRE, has_toAligned_XLOW, has_toAligned_YLOW; ///< Have phase shifts for shift to field aligned coordinates been calculated + bool has_fromAligned_CENTRE, has_fromAligned_XLOW, has_fromAligned_YLOW; ///< Have phase shifts for shift from field aligned coordinates been calculated + bool has_yupPhs1_CENTRE, has_yupPhs1_XLOW, has_yupPhs1_YLOW; ///< Have phase shifts for yup1 been calculated + bool has_ydownPhs1_CENTRE, has_ydownPhs1_XLOW, has_ydownPhs1_YLOW; ///< Have phase shifts for ydown1 been calculated + bool has_yupPhs2_CENTRE, has_yupPhs2_XLOW, has_yupPhs2_YLOW; ///< Have phase shifts for yup2 been calculated + bool has_ydownPhs2_CENTRE, has_ydownPhs2_XLOW, has_ydownPhs2_YLOW; ///< Have phase shifts for ydown2 been calculated + + Matrix< Array > toAlignedPhs_CENTRE; ///< Cache of phase shifts for transforming from X-Z orthogonal coordinates to field-aligned coordinates. Cell centre version. + Matrix< Array > fromAlignedPhs_CENTRE; ///< Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates. Cell centre version. + Matrix< Array > toAlignedPhs_XLOW; ///< Cache of phase shifts for transforming from X-Z orthogonal coordinates to field-aligned coordinates. Interpolated to CELL_XLOW. + Matrix< Array > fromAlignedPhs_XLOW; ///< Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates. Interpolated to CELL_XLOW. + Matrix< Array > toAlignedPhs_YLOW; ///< Cache of phase shifts for transforming from X-Z orthogonal coordinates to field-aligned coordinates. Interpolated to CELL_YLOW. + Matrix< Array > fromAlignedPhs_YLOW; ///< Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates. Interpolated to CELL_YLOW. + + Matrix< Array > yupPhs1_CENTRE; ///< Cache of phase shifts for calculating yup1 fields. Cell centre version. + Matrix< Array > ydownPhs1_CENTRE; ///< Cache of phase shifts for calculating ydown1 fields. Cell centre version. + Matrix< Array > yupPhs2_CENTRE; ///< Cache of phase shifts for calculating yup2 fields. Cell centre version. + Matrix< Array > ydownPhs2_CENTRE; ///< Cache of phase shifts for calculating ydown2 fields. Cell centre version. + Matrix< Array > yupPhs1_XLOW; ///< Cache of phase shifts for calculating yup1 fields. Interpolated to CELL_XLOW. + Matrix< Array > ydownPhs1_XLOW; ///< Cache of phase shifts for calculating ydown1 fields. Interpolated to CELL_XLOW. + Matrix< Array > yupPhs2_XLOW; ///< Cache of phase shifts for calculating yup2 fields. Interpolated to CELL_XLOW. + Matrix< Array > ydownPhs2_XLOW; ///< Cache of phase shifts for calculating ydown2 fields. Interpolated to CELL_XLOW. + Matrix< Array > yupPhs1_YLOW; ///< Cache of phase shifts for calculating yup1 fields. Interpolated to CELL_YLOW. + Matrix< Array > ydownPhs1_YLOW; ///< Cache of phase shifts for calculating ydown1 fields. Interpolated to CELL_YLOW. + Matrix< Array > yupPhs2_YLOW; ///< Cache of phase shifts for calculating yup2 fields. Interpolated to CELL_YLOW. + Matrix< Array > ydownPhs2_YLOW; ///< Cache of phase shifts for calculating ydown2 fields. Interpolated to CELL_YLOW. /*! * 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 @@ -152,7 +187,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 @@ -163,7 +198,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 Matrix< Array > &phs, const REGION region=RGN_NOX); /*! * Shift a given 1D array, assumed to be in Z, by the given \p zangle @@ -182,7 +217,10 @@ private: * @param[in] phs Phase shift, assumed to have length (mesh.LocalNz/2 + 1) i.e. the number of modes * @param[out] out A 1D array of length mesh.LocalNz, already allocated */ - void shiftZ(const BoutReal *in, const std::vector &phs, BoutReal *out); + void shiftZ(const BoutReal *in, const Array &phs, BoutReal *out); + + /// Write out ParallelTransform variables to file + void outputVars(Datafile &file); }; diff --git a/include/field3d.hxx b/include/field3d.hxx index 6e09d6f583..7ea5d5e6b5 100644 --- a/include/field3d.hxx +++ b/include/field3d.hxx @@ -151,21 +151,21 @@ template class Flexible; Field3D f(0.0); // f allocated, set to zero - f.yup() // error; f.yup not allocated + f.yup(1) // 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) + f.mergeYupYdown(); // f.yup(i) and f.ydown(i) now point to f + f.yup(1)(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.splitYupYdown(); // f.yup(i) and f.ydown(i) separate - f.yup(); // ok - f.yup()(0,1,0) // error; f.yup not allocated + f.yup(1); // ok + f.yup(1)(0,1,0) // error; f.yup not allocated - f.yup() = 1.0; // Set f.yup() field to 1.0 + f.yup(1) = 1.0; // Set f.yup() field to 1.0 - f.yup()(0,1,0) // ok + f.yup(1)(0,1,0) // ok */ class Field3D : public Field, public FieldData { @@ -234,36 +234,30 @@ class Field3D : public Field, public FieldData { * Ensure that yup and ydown refer to this field */ void mergeYupYdown(); + + /*! + * Clear all yup/ydown fields and field_fa + */ + void clearYupYdown(); /// Check if this field has yup and ydown fields bool hasYupYdown() const { - return (yup_field != nullptr) && (ydown_field != nullptr); + return (yup1_field != nullptr) && (ydown1_field != nullptr); } /// Return reference to yup field - Field3D& yup() { - ASSERT2(yup_field != nullptr); // Check for communicate - return *yup_field; - } - /// Return const reference to yup field - const Field3D& yup() const { - ASSERT2(yup_field != nullptr); - return *yup_field; - } + Field3D& yup(const int i = 1); + /// Return const reference to yup field + const Field3D& yup(const int i = 1) const; + /// Return reference to ydown field - Field3D& ydown() { - ASSERT2(ydown_field != nullptr); - return *ydown_field; - } + Field3D& ydown(const int i = 1); /// Return const reference to ydown field - const Field3D& ydown() const { - ASSERT2(ydown_field != nullptr); - return *ydown_field; - } + const Field3D& ydown(const int i = 1) const; - /// Return yup if dir=+1, and ydown if dir=-1 + /// Return yup(dir) if dir>0, and ydown(-dir) if dir<0 Field3D& ynext(int dir); const Field3D& ynext(int dir) const; @@ -528,8 +522,8 @@ private: Field3D *deriv; ///< Time derivative (may be NULL) - /// Pointers to fields containing values along Y - Field3D *yup_field, *ydown_field; + /// Arrays of pointers to fields containing values along Y + Field3D *yup1_field, *ydown1_field, *yup2_field, *ydown2_field; }; // Non-member overloaded operators diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index a0ea01e69b..e76d28a353 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -45,8 +45,8 @@ /// Constructor Field3D::Field3D(Mesh *localmesh) - : Field(localmesh), background(nullptr), deriv(nullptr), yup_field(nullptr), - ydown_field(nullptr) { + : Field(localmesh), background(nullptr), deriv(nullptr), yup1_field(nullptr), + ydown1_field(nullptr), yup2_field(nullptr), ydown2_field(nullptr) { #ifdef TRACK name = ""; #endif @@ -74,7 +74,8 @@ Field3D::Field3D(Mesh *localmesh) Field3D::Field3D(const Field3D &f) : Field(f.fieldmesh), // The mesh containing array sizes background(nullptr), data(f.data), // This handles references to the data array - deriv(nullptr), yup_field(nullptr), ydown_field(nullptr) { + deriv(nullptr), yup1_field(nullptr), + ydown1_field(nullptr), yup2_field(nullptr), ydown2_field(nullptr) { TRACE("Field3D(Field3D&)"); @@ -101,8 +102,8 @@ Field3D::Field3D(const Field3D &f) } Field3D::Field3D(const Field2D &f) - : Field(f.getMesh()), background(nullptr), deriv(nullptr), yup_field(nullptr), - ydown_field(nullptr) { + : Field(f.getMesh()), background(nullptr), deriv(nullptr), yup1_field(nullptr), + ydown1_field(nullptr), yup2_field(nullptr), ydown2_field(nullptr) { TRACE("Field3D: Copy constructor from Field2D"); @@ -119,8 +120,8 @@ Field3D::Field3D(const Field2D &f) } Field3D::Field3D(const BoutReal val, Mesh *localmesh) - : Field(localmesh), background(nullptr), deriv(nullptr), yup_field(nullptr), - ydown_field(nullptr) { + : Field(localmesh), background(nullptr), deriv(nullptr), yup1_field(nullptr), + ydown1_field(nullptr), yup2_field(nullptr), ydown2_field(nullptr) { TRACE("Field3D: Copy constructor from value"); @@ -141,20 +142,20 @@ Field3D::~Field3D() { // The ddt of the yup/ydown_fields point to the same place as ddt.yup_field // only delete once // Also need to check that separate yup_field exists - if ((yup_field != this) && (yup_field != nullptr)) - yup_field->deriv = nullptr; - if ((ydown_field != this) && (ydown_field != nullptr)) - ydown_field->deriv = nullptr; + if ((yup1_field != this) && (yup1_field != nullptr)) + yup1_field->deriv = nullptr; + if ((ydown1_field != this) && (ydown1_field != nullptr)) + ydown1_field->deriv = nullptr; + if ((yup2_field != this) && (yup2_field != nullptr)) + yup2_field->deriv = nullptr; + if ((ydown2_field != this) && (ydown2_field != nullptr)) + ydown2_field->deriv = nullptr; // 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; + + clearYupYdown(); } void Field3D::allocate() { @@ -184,51 +185,154 @@ Field3D* Field3D::timeDeriv() { void Field3D::splitYupYdown() { TRACE("Field3D::splitYupYdown"); - if((yup_field != this) && (yup_field != nullptr)) + if((yup1_field != this) && (yup1_field != nullptr)) return; - // yup_field and ydown_field null - yup_field = new Field3D(fieldmesh); - ydown_field = new Field3D(fieldmesh); + // yup_array and ydown_array null + yup1_field = new Field3D(fieldmesh); + yup1_field->setLocation(location); + ydown1_field = new Field3D(fieldmesh); + ydown1_field->setLocation(location); + if (fieldmesh->ystart>1) { + yup2_field = new Field3D(fieldmesh); + yup2_field->setLocation(location); + ydown2_field = new Field3D(fieldmesh); + ydown2_field->setLocation(location); + } } void Field3D::mergeYupYdown() { TRACE("Field3D::mergeYupYdown"); - if(yup_field == this && ydown_field == this) + if (yup1_field == this && ydown1_field == this) return; - if(yup_field != nullptr){ - delete yup_field; + clearYupYdown(); + + yup1_field = this; + ydown1_field = this; + if (fieldmesh->ystart>1) { + yup2_field = this; + ydown2_field = this; + } +} + +void Field3D::clearYupYdown() { + // Delete auxiliary fields if they have been set + if (yup1_field != nullptr && yup1_field != this) { + delete yup1_field; + } + yup1_field = nullptr; + + if (ydown1_field != nullptr && ydown1_field != this) { + delete ydown1_field; + } + ydown1_field = nullptr; + + if (yup2_field != nullptr && yup2_field != this) { + delete yup2_field; + } + yup2_field = nullptr; + + if (ydown2_field != nullptr && ydown2_field != this) { + delete ydown2_field; } + ydown2_field = nullptr; +} - if(ydown_field != nullptr) { - delete ydown_field; +Field3D& Field3D::yup(const int i) { + switch (i) { + case 1: { + ASSERT2(yup1_field != nullptr); // Check for communicate + return *yup1_field; } + case 2: { + ASSERT2(yup2_field != nullptr); // Check for communicate + return *yup2_field; + } +#if CHECK > 1 + default: { + throw BoutException("There is no yup field for i=%i", i); + } +#endif + } +} - yup_field = this; - ydown_field = this; +const Field3D& Field3D::yup(const int i) const { + switch (i) { + case 1: { + ASSERT2(yup1_field != nullptr); // Check for communicate + return *yup1_field; + } + case 2: { + ASSERT2(yup2_field != nullptr); // Check for communicate + return *yup2_field; + } +#if CHECK > 1 + default: { + throw BoutException("There is no yup field for i=%i", i); + } +#endif + } +} + +Field3D& Field3D::ydown(const int i) { + switch (i) { + case 1: { + ASSERT2(ydown1_field != nullptr); // Check for communicate + return *ydown1_field; + } + case 2: { + ASSERT2(ydown2_field != nullptr); // Check for communicate + return *ydown2_field; + } +#if CHECK > 1 + default: { + throw BoutException("There is no ydown field for i=%i", i); + } +#endif + } +} + +const Field3D& Field3D::ydown(const int i) const { + switch (i) { + case 1: { + ASSERT2(ydown1_field != nullptr); // Check for communicate + return *ydown1_field; + } + case 2: { + ASSERT2(ydown2_field != nullptr); // Check for communicate + return *ydown2_field; + } +#if CHECK > 1 + default: { + throw BoutException("There is no ydown field for i=%i", i); + } +#endif + } } Field3D& Field3D::ynext(int dir) { - switch(dir) { - case +1: - return yup(); - case -1: - return ydown(); - default: - throw BoutException("Field3D: Call to ynext with strange direction %d. Only +/-1 currently supported", dir); + if (dir>0) { + return yup(dir); + } else if (dir<0) { + return ydown(-dir); + } else { +#if CHECK > 1 + throw BoutException("Field3D: Call to ynext with strange direction %d. Should not be zero.", dir); +#endif } } const Field3D& Field3D::ynext(int dir) const { - switch(dir) { - case +1: - return yup(); - case -1: - return ydown(); - default: - throw BoutException("Field3D: Call to ynext with strange direction %d. Only +/-1 currently supported", dir); + if (dir>0) { + return yup(dir); + } else if (dir<0) { + return ydown(-dir); + } else { +#if CHECK > 1 + throw BoutException("Field3D: Call to ynext with strange direction %d. Should not be zero.", dir); +#endif } } @@ -357,7 +461,9 @@ Field3D & Field3D::operator=(const Field3D &rhs) { location = rhs.location; - return *this; + clearYupYdown(); + + return *this; } Field3D & Field3D::operator=(const Field2D &rhs) { @@ -377,7 +483,9 @@ Field3D & Field3D::operator=(const Field2D &rhs) { /// Only 3D fields have locations for now //location = CELL_CENTRE; - + + clearYupYdown(); + return *this; } @@ -408,6 +516,8 @@ Field3D & Field3D::operator=(const BoutReal val) { //location = CELL_CENTRE; // DON'T RE-SET LOCATION + clearYupYdown(); + return *this; } diff --git a/src/field/gen_fieldops.jinja b/src/field/gen_fieldops.jinja index 0f529132c8..3aa4975570 100644 --- a/src/field/gen_fieldops.jinja +++ b/src/field/gen_fieldops.jinja @@ -95,6 +95,11 @@ } else { (*this) = (*this) {{operator}} {{rhs.name}}; } + + {% if (out == "Field3D") %} + clearYupYdown(); + {% endif %} + return *this; } {% endif %} diff --git a/src/field/generated_fieldops.cxx b/src/field/generated_fieldops.cxx index 5a3a938065..99c227dbba 100644 --- a/src/field/generated_fieldops.cxx +++ b/src/field/generated_fieldops.cxx @@ -61,6 +61,9 @@ Field3D &Field3D::operator*=(const Field3D &rhs) { } else { (*this) = (*this) * rhs; } + + clearYupYdown(); + return *this; } @@ -119,6 +122,9 @@ Field3D &Field3D::operator/=(const Field3D &rhs) { } else { (*this) = (*this) / rhs; } + + clearYupYdown(); + return *this; } @@ -177,6 +183,9 @@ Field3D &Field3D::operator+=(const Field3D &rhs) { } else { (*this) = (*this) + rhs; } + + clearYupYdown(); + return *this; } @@ -235,6 +244,9 @@ Field3D &Field3D::operator-=(const Field3D &rhs) { } else { (*this) = (*this) - rhs; } + + clearYupYdown(); + return *this; } @@ -298,6 +310,9 @@ Field3D &Field3D::operator*=(const Field2D &rhs) { } else { (*this) = (*this) * rhs; } + + clearYupYdown(); + return *this; } @@ -362,6 +377,9 @@ Field3D &Field3D::operator/=(const Field2D &rhs) { } else { (*this) = (*this) / rhs; } + + clearYupYdown(); + return *this; } @@ -425,6 +443,9 @@ Field3D &Field3D::operator+=(const Field2D &rhs) { } else { (*this) = (*this) + rhs; } + + clearYupYdown(); + return *this; } @@ -488,6 +509,9 @@ Field3D &Field3D::operator-=(const Field2D &rhs) { } else { (*this) = (*this) - rhs; } + + clearYupYdown(); + return *this; } @@ -526,6 +550,9 @@ Field3D &Field3D::operator*=(const BoutReal rhs) { } else { (*this) = (*this) * rhs; } + + clearYupYdown(); + return *this; } @@ -564,6 +591,9 @@ Field3D &Field3D::operator/=(const BoutReal rhs) { } else { (*this) = (*this) / rhs; } + + clearYupYdown(); + return *this; } @@ -602,6 +632,9 @@ Field3D &Field3D::operator+=(const BoutReal rhs) { } else { (*this) = (*this) + rhs; } + + clearYupYdown(); + return *this; } @@ -640,6 +673,9 @@ Field3D &Field3D::operator-=(const BoutReal rhs) { } else { (*this) = (*this) - rhs; } + + clearYupYdown(); + return *this; } @@ -822,6 +858,7 @@ Field2D &Field2D::operator*=(const Field2D &rhs) { } else { (*this) = (*this) * rhs; } + return *this; } @@ -880,6 +917,7 @@ Field2D &Field2D::operator/=(const Field2D &rhs) { } else { (*this) = (*this) / rhs; } + return *this; } @@ -938,6 +976,7 @@ Field2D &Field2D::operator+=(const Field2D &rhs) { } else { (*this) = (*this) + rhs; } + return *this; } @@ -996,6 +1035,7 @@ Field2D &Field2D::operator-=(const Field2D &rhs) { } else { (*this) = (*this) - rhs; } + return *this; } @@ -1034,6 +1074,7 @@ Field2D &Field2D::operator*=(const BoutReal rhs) { } else { (*this) = (*this) * rhs; } + return *this; } @@ -1072,6 +1113,7 @@ Field2D &Field2D::operator/=(const BoutReal rhs) { } else { (*this) = (*this) / rhs; } + return *this; } @@ -1110,6 +1152,7 @@ Field2D &Field2D::operator+=(const BoutReal rhs) { } else { (*this) = (*this) + rhs; } + return *this; } @@ -1148,6 +1191,7 @@ Field2D &Field2D::operator-=(const BoutReal rhs) { } else { (*this) = (*this) - rhs; } + return *this; } diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 7a8ba661a5..954230331e 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -796,6 +796,11 @@ const Field3D Coordinates::Div_par(const Field3D &f, CELL_LOC outloc, f_B.splitYupYdown(); f_B.yup() = f.yup() / Bxy; f_B.ydown() = f.ydown() / Bxy; + if (mesh->ystart > 1) { + // Have a second yup/ydown field + f_B.yup(2) = f.yup(2) / Bxy; + f_B.ydown(2) = f.ydown(2) / Bxy; + } } return Bxy * Grad_par(f_B, outloc, method); } diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index f688d19cea..cf4aef9acc 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -341,63 +341,11 @@ const Field3D Vpar_Grad_par_LCtoC(const Field3D &v, const Field3D &f, REGION reg result[i] -= (vval.p >= 0.0) ? vval.p * fval.c : vval.p * fval.p; } } - else if (vUseUpDown) { - // Only v has up/down fields - // f must shift to field aligned coordinates - Field3D f_fa = mesh->toFieldAligned(f); - - vval.mm = nan(""); - vval.pp = nan(""); - - for (const auto &i : result.region(region)) { - - fval.mm = f_fa[i.offset(0, -2, 0)]; - fval.m = f_fa[i.ym()]; - fval.c = f_fa[i]; - fval.p = f_fa[i.yp()]; - fval.pp = f_fa[i.offset(0, 2, 0)]; - - vval.m = v.ydown()[i.ym()]; - vval.c = v[i]; - vval.p = v.yup()[i.yp()]; - - // Left side - result[i] = (vval.c >= 0.0) ? vval.c * fval.m : vval.c * fval.c; - // Right side - result[i] -= (vval.p >= 0.0) ? vval.p * fval.c : vval.p * fval.p; - } - } - else if (fUseUpDown) { - // Only f has up/down fields - // v must shift to field aligned coordinates - Field3D v_fa = mesh->toFieldAligned(v); - - stencil vval; - - stencil fval; - fval.mm = nan(""); - fval.pp = nan(""); - - for (const auto &i : result.region(region)) { - - fval.m = f.ydown()[i.ym()]; - fval.c = f[i]; - fval.p = f.yup()[i.yp()]; - - vval.mm = v_fa[i.offset(0,-2,0)]; - vval.m = v_fa[i.ym()]; - vval.c = v_fa[i]; - vval.p = v_fa[i.yp()]; - vval.pp = v_fa[i.offset(0,2,0)]; - - // Left side - result[i] = (vval.c >= 0.0) ? vval.c * fval.m : vval.c * fval.c; - // Right side - result[i] -= (vval.p >= 0.0) ? vval.p * fval.c : vval.p * fval.p; - } - } else { // Both must shift to field aligned + // (even if one of v and f has yup/ydown fields, it doesn't make sense to + // multiply them with one in field-aligned and one in non-field-aligned + // coordinates) Field3D v_fa = mesh->toFieldAligned(v); Field3D f_fa = mesh->toFieldAligned(f); @@ -420,6 +368,8 @@ const Field3D Vpar_Grad_par_LCtoC(const Field3D &v, const Field3D &f, REGION reg // Right side result[i] -= (vval.p >= 0.0) ? vval.p * fval.c : vval.p * fval.p; } + + result = mesh->fromFieldAligned(result); } result.setLocation(CELL_CENTRE); @@ -431,10 +381,11 @@ const Field3D Grad_par_LtoC(const Field3D &var) { ASSERT1(var.getLocation() == CELL_YLOW); } - Field3D result(var.getMesh()); + Mesh* fieldmesh = var.getMesh(); + Field3D result(fieldmesh); result.allocate(); - Coordinates *metric = var.getMesh()->coordinates(); + Coordinates *metric = fieldmesh->coordinates(); if (var.hasYupYdown()) { for (auto &i : result.region(RGN_NOBNDRY)) { @@ -443,12 +394,13 @@ const Field3D Grad_par_LtoC(const Field3D &var) { } else { // No yup/ydown field, so transform to field aligned - Field3D var_fa = var.getMesh()->toFieldAligned(var); + Field3D var_fa = fieldmesh->toFieldAligned(var); for(auto &i : result.region(RGN_NOBNDRY)) { result[i] = (var_fa[i.yp()] - var_fa[i]) / (metric->dy[i]*sqrt(metric->g_22.get(CELL_CENTRE).operator[](i))); } - result = var.getMesh()->fromFieldAligned(result); + + result = fieldmesh->fromFieldAligned(result); } result.setLocation(CELL_CENTRE); diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index 91eca272bf..5e758e85e1 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -2552,4 +2552,5 @@ void BoutMesh::outputVars(Datafile &file) { file.add(jyseps2_2, "jyseps2_2", 0); coordinates()->outputVars(file); + getParallelTransform().outputVars(file); } diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx index ace0fe0af8..c50f43b180 100644 --- a/src/mesh/index_derivs.cxx +++ b/src/mesh/index_derivs.cxx @@ -677,10 +677,9 @@ const Field2D Mesh::applyXdiff(const Field2D &var, Mesh::deriv_func func, return Field2D(0., this); } - CELL_LOC diffloc = var.getLocation(); - Field2D result(this); result.allocate(); // Make sure data allocated + result.setLocation(loc); if (this->StaggerGrids && (loc != CELL_DEFAULT) && (loc != var.getLocation())) { // Staggered differencing @@ -765,8 +764,6 @@ const Field2D Mesh::applyXdiff(const Field2D &var, Mesh::deriv_func func, } } - result.setLocation(diffloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -786,10 +783,9 @@ const Field3D Mesh::applyXdiff(const Field3D &var, Mesh::deriv_func func, return Field3D(0., this); } - CELL_LOC diffloc = var.getLocation(); - Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(loc); if (this->StaggerGrids && (loc != CELL_DEFAULT) && (loc != var.getLocation())) { // Staggered differencing @@ -874,8 +870,6 @@ const Field3D Mesh::applyXdiff(const Field3D &var, Mesh::deriv_func func, } } - result.setLocation(diffloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -900,6 +894,7 @@ const Field2D Mesh::applyYdiff(const Field2D &var, Mesh::deriv_func func, CELL_L Field2D result(this); result.allocate(); // Make sure data allocated + result.setLocation(diffloc); if (this->ystart > 1) { // More than one guard cell, so set pp and mm values @@ -931,8 +926,6 @@ const Field2D Mesh::applyYdiff(const Field2D &var, Mesh::deriv_func func, CELL_L } } - result.setLocation(diffloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_yup = result.bndry_ydown = false; @@ -951,10 +944,9 @@ const Field3D Mesh::applyYdiff(const Field3D &var, Mesh::deriv_func func, CELL_L return Field3D(0., this); } - CELL_LOC diffloc = var.getLocation(); - Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(loc); if (var.hasYupYdown() && ((&var.yup() != &var) || (&var.ydown() != &var))) { // Field "var" has distinct yup and ydown fields which @@ -968,38 +960,78 @@ const Field3D Mesh::applyYdiff(const Field3D &var, Mesh::deriv_func func, CELL_L CELL_LOC location = var.getLocation(); stencil s; - s.pp = nan(""); - s.mm = nan(""); - for (const auto &i : result.region(region)) { - // Set stencils - s.c = var[i]; - s.p = var.yup()[i.yp()]; - s.m = var.ydown()[i.ym()]; + if (mesh->ystart==1) { + // Only one guard cell, so can only use 3-point stencils + s.pp = nan(""); + s.mm = nan(""); + for (const auto &i : result.region(region)) { + // Set stencils + s.c = var[i]; + s.p = var.yup()[i.yp()]; + s.m = var.ydown()[i.ym()]; - if ((location == CELL_CENTRE) && (loc == CELL_YLOW)) { - // Producing a stencil centred around a lower Y value - s.pp = s.p; - s.p = s.c; - } else if (location == CELL_YLOW) { - // Stencil centred around a cell centre - s.mm = s.m; - s.m = s.c; + if ((location == CELL_CENTRE) && (loc == CELL_YLOW)) { + // Producing a stencil centred around a lower Y value + s.pp = s.p; + s.p = s.c; + } else if (location == CELL_YLOW) { + // Stencil centred around a cell centre + s.mm = s.m; + s.m = s.c; + } + + result[i] = func(s); } + } else { + // Can use 5-point stencils + for (const auto &i : result.region(region)) { + // Set stencils + s.c = var[i]; + s.p = var.yup()[i.yp()]; + s.m = var.ydown()[i.ym()]; + s.pp = var.yup(2)[i.offset(0, 2, 0)]; + s.mm = var.ydown(2)[i.offset(0, -2, 0)]; - result[i] = func(s); + if ((location == CELL_CENTRE) && (loc == CELL_YLOW)) { + // Producing a stencil centred around a lower Y value + s.pp = s.p; + s.p = s.c; + } else if (location == CELL_YLOW) { + // Stencil centred around a cell centre + s.mm = s.m; + s.m = s.c; + } + + result[i] = func(s); + } } } else { // Non-staggered stencil s; - s.pp = nan(""); - s.mm = nan(""); - for (const auto &i : result.region(region)) { - // Set stencils - s.c = var[i]; - s.p = var.yup()[i.yp()]; - s.m = var.ydown()[i.ym()]; + if (mesh->ystart == 1) { + // Only one guard cell, so can only use 3-point stencils + s.pp = nan(""); + s.mm = nan(""); + for (const auto &i : result.region(region)) { + // Set stencils + s.c = var[i]; + s.p = var.yup()[i.yp()]; + s.m = var.ydown()[i.ym()]; - result[i] = func(s); + result[i] = func(s); + } + } else { + // Can use 5-point stencils + for (const auto &i : result.region(region)) { + // Set stencils + s.c = var[i]; + s.p = var.yup()[i.yp()]; + s.m = var.ydown()[i.ym()]; + s.pp = var.yup(2)[i.offset(0, 2, 0)]; + s.mm = var.ydown(2)[i.offset(0, -2, 0)]; + + result[i] = func(s); + } } } } else { @@ -1100,11 +1132,10 @@ const Field3D Mesh::applyYdiff(const Field3D &var, Mesh::deriv_func func, CELL_L result = this->fromFieldAligned(result); } - result.setLocation(diffloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; + invalidateGuards(result); // extra check: set guard cells to NaN if CHECK>2 #endif return result; @@ -1122,9 +1153,6 @@ const Field3D Mesh::applyZdiff(const Field3D &var, Mesh::deriv_func func, CELL_L return Field3D(0., this); } - - CELL_LOC diffloc = var.getLocation(); - if (this->StaggerGrids && (loc != CELL_DEFAULT) && (loc != var.getLocation())) { // Staggered differencing throw BoutException("No one used this before. And no one implemented it."); @@ -1132,6 +1160,7 @@ const Field3D Mesh::applyZdiff(const Field3D &var, Mesh::deriv_func func, CELL_L Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(loc); // Check that the input variable has data ASSERT1(var.isAllocated()); @@ -1147,8 +1176,6 @@ const Field3D Mesh::applyZdiff(const Field3D &var, Mesh::deriv_func func, CELL_L result[i] = func(s); } - result.setLocation(diffloc); - return result; } @@ -1166,8 +1193,6 @@ const Field3D Mesh::indexDDX(const Field3D &f, CELL_LOC outloc, DIFF_METHOD meth CELL_LOC inloc = f.getLocation(); // Input location CELL_LOC diffloc = inloc; // Location of differential result - Field3D result(this); - if (this->StaggerGrids && (outloc == CELL_DEFAULT)) { // Take care of CELL_DEFAULT case outloc = diffloc; // No shift (i.e. same as no stagger case) @@ -1199,7 +1224,7 @@ const Field3D Mesh::indexDDX(const Field3D &f, CELL_LOC outloc, DIFF_METHOD meth throw BoutException("Cannot use FFT for X derivatives"); } - result = applyXdiff(f, func, diffloc, region); + Field3D result = applyXdiff(f, func, diffloc, region); result.setLocation(diffloc); // Set the result location @@ -1223,8 +1248,6 @@ const Field3D Mesh::indexDDY(const Field3D &f, CELL_LOC outloc, CELL_LOC inloc = f.getLocation(); // Input location CELL_LOC diffloc = inloc; // Location of differential result - Field3D result(this); - if (this->StaggerGrids && (outloc == CELL_DEFAULT)) { // Take care of CELL_DEFAULT case outloc = diffloc; // No shift (i.e. same as no stagger case) @@ -1255,7 +1278,7 @@ const Field3D Mesh::indexDDY(const Field3D &f, CELL_LOC outloc, throw BoutException("Cannot use FFT for Y derivatives"); } - result = applyYdiff(f, func, diffloc, region); + Field3D result = applyYdiff(f, func, diffloc, region); result.setLocation(diffloc); // Set the result location @@ -1426,8 +1449,6 @@ const Field3D Mesh::indexD2DX2(const Field3D &f, CELL_LOC outloc, ASSERT1(this == f.getMesh()); - Field3D result(this); - if (StaggerGrids && (outloc == CELL_DEFAULT)) { // Take care of CELL_DEFAULT case outloc = diffloc; // No shift (i.e. same as no stagger case) @@ -1459,7 +1480,7 @@ const Field3D Mesh::indexD2DX2(const Field3D &f, CELL_LOC outloc, throw BoutException("Cannot use FFT for X derivatives"); } - result = applyXdiff(f, func, diffloc, region); + Field3D result = applyXdiff(f, func, diffloc, region); result.setLocation(diffloc); @@ -1505,8 +1526,6 @@ const Field3D Mesh::indexD2DY2(const Field3D &f, CELL_LOC outloc, ASSERT1(this == f.getMesh()); - Field3D result(this); - if (StaggerGrids && (outloc == CELL_DEFAULT)) { // Take care of CELL_DEFAULT case outloc = diffloc; // No shift (i.e. same as no stagger case) @@ -1538,7 +1557,7 @@ const Field3D Mesh::indexD2DY2(const Field3D &f, CELL_LOC outloc, throw BoutException("Cannot use FFT for Y derivatives"); } - result = applyYdiff(f, func, diffloc, region); + Field3D result = applyYdiff(f, func, diffloc, region); result.setLocation(diffloc); @@ -1761,6 +1780,7 @@ const Field2D Mesh::indexVDDX(const Field2D &v, const Field2D &f, CELL_LOC outlo Field2D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (this->xstart > 1) { // Two or more guard cells @@ -1800,8 +1820,6 @@ const Field2D Mesh::indexVDDX(const Field2D &v, const Field2D &f, CELL_LOC outlo result.bndry_xin = result.bndry_xout = false; #endif - result.setLocation(diffloc); - return result; } @@ -1816,6 +1834,7 @@ const Field3D Mesh::indexVDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); CELL_LOC vloc = v.getLocation(); CELL_LOC inloc = f.getLocation(); // Input location @@ -1989,8 +2008,6 @@ const Field3D Mesh::indexVDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo } } - result.setLocation(diffloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -2012,6 +2029,7 @@ const Field2D Mesh::indexVDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo Field2D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); CELL_LOC vloc = v.getLocation(); CELL_LOC inloc = f.getLocation(); // Input location @@ -2024,7 +2042,6 @@ const Field2D Mesh::indexVDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo if (this->LocalNy == 1){ result=0; - result.setLocation(outloc); return result; } @@ -2187,8 +2204,6 @@ const Field2D Mesh::indexVDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo } } - result.setLocation(diffloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -2207,6 +2222,7 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); CELL_LOC vloc = v.getLocation(); CELL_LOC inloc = f.getLocation(); // Input location @@ -2219,7 +2235,6 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo if (this->LocalNy == 1){ result=0; - result.setLocation(outloc); return result; } @@ -2252,12 +2267,8 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo func = lookupFluxFunc(table, method); } - // There are four cases, corresponding to whether or not f and v - // have yup, ydown fields. - - // If vUseUpDown is true, field "v" has distinct yup and ydown fields which - // will be used to calculate a derivative along - // the magnetic field + // If *UseUpDown is true, field "*" 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))); @@ -2265,104 +2276,72 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo // Both v and f have up/down fields stencil vval, fval; - vval.pp = nan(""); - vval.mm = nan(""); - fval.pp = nan(""); - fval.mm = nan(""); - for (const auto &i : result.region(region)) { - vval.c = v[i]; - vval.p = v.yup()[i.yp()]; - vval.m = v.ydown()[i.ym()]; - fval.c = f[i]; - fval.p = f.yup()[i.yp()]; - fval.m = f.ydown()[i.ym()]; - - if (diffloc != CELL_DEFAULT) { - // Non-centred stencil - if ((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) { - // Producing a stencil centred around a lower Y value - vval.pp = vval.p; - vval.p = vval.c; - } else if (vloc == CELL_YLOW) { - // Stencil centred around a cell centre - vval.mm = vval.m; - vval.m = vval.c; - } - // Shifted in one direction -> shift in another - // Could produce warning - } - result[i] = func(vval, fval); - } - } else if (vUseUpDown) { - // Only v has up/down fields - // f must shift to field aligned coordinates - Field3D f_fa = this->toFieldAligned(f); - - stencil vval, fval; - vval.pp = nan(""); - vval.mm = nan(""); - for (const auto &i : result.region(region)) { - vval.c = v[i]; - vval.p = v.yup()[i.yp()]; - vval.m = v.ydown()[i.ym()]; - fval.c = f_fa[i]; - fval.p = f_fa[i.yp()]; - fval.m = f_fa[i.ym()]; - fval.pp = f_fa[i.offset(0, 2, 0)]; - fval.mm = f_fa[i.offset(0, -2, 0)]; - - if (diffloc != CELL_DEFAULT) { - // Non-centred stencil - if ((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) { - // Producing a stencil centred around a lower Y value - vval.pp = vval.p; - vval.p = vval.c; - } else if (vloc == CELL_YLOW) { - // Stencil centred around a cell centre - vval.mm = vval.m; - vval.m = vval.c; + if (mesh->ystart==1) { + // Only one guard cell, so can only use 3-point stencils + vval.pp = nan(""); + vval.mm = nan(""); + fval.pp = nan(""); + fval.mm = nan(""); + for (const auto &i : result.region(region)) { + vval.c = v[i]; + vval.p = v.yup()[i.yp()]; + vval.m = v.ydown()[i.ym()]; + fval.c = f[i]; + fval.p = f.yup()[i.yp()]; + fval.m = f.ydown()[i.ym()]; + + if (diffloc != CELL_DEFAULT) { + // Non-centred stencil + if ((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) { + // Producing a stencil centred around a lower Y value + vval.pp = vval.p; + vval.p = vval.c; + } else if (vloc == CELL_YLOW) { + // Stencil centred around a cell centre + vval.mm = vval.m; + vval.m = vval.c; + } + // Shifted in one direction -> shift in another + // Could produce warning } - // Shifted in one direction -> shift in another - // Could produce warning + result[i] = func(vval, fval); } - result[i] = func(vval, fval); - } - } else if (fUseUpDown) { - // Only f has up/down fields - // v must shift to field aligned coordinates - Field3D v_fa = this->toFieldAligned(v); - - stencil vval, fval; - fval.pp = nan(""); - fval.mm = nan(""); - for (const auto &i : result.region(region)) { - vval.c = v_fa[i]; - vval.p = v_fa[i.yp()]; - vval.m = v_fa[i.ym()]; - vval.pp = v_fa[i.offset(0, 2, 0)]; - vval.mm = v_fa[i.offset(0, -2, 0)]; - fval.c = f[i]; - fval.p = f.yup()[i.yp()]; - fval.m = f.ydown()[i.ym()]; - - if (diffloc != CELL_DEFAULT) { - // Non-centred stencil - if ((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) { - // Producing a stencil centred around a lower Y value - vval.pp = vval.p; - vval.p = vval.c; - } else if (vloc == CELL_YLOW) { - // Stencil centred around a cell centre - vval.mm = vval.m; - vval.m = vval.c; + } else { + // Can use 5-point stencils + for (const auto &i : result.region(region)) { + vval.c = v[i]; + vval.p = v.yup()[i.yp()]; + vval.m = v.ydown()[i.ym()]; + vval.pp = v.yup(2)[i.offset(0, 2, 0)]; + vval.mm = v.ydown(2)[i.offset(0, -2, 0)]; + fval.c = f[i]; + fval.p = f.yup()[i.yp()]; + fval.m = f.ydown()[i.ym()]; + fval.pp = f.yup(2)[i.offset(0, 2, 0)]; + fval.mm = f.ydown(2)[i.offset(0, -2, 0)]; + + if (diffloc != CELL_DEFAULT) { + // Non-centred stencil + if ((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) { + // Producing a stencil centred around a lower Y value + vval.pp = vval.p; + vval.p = vval.c; + } else if (vloc == CELL_YLOW) { + // Stencil centred around a cell centre + vval.mm = vval.m; + vval.m = vval.c; + } + // Shifted in one direction -> shift in another + // Could produce warning } - // Shifted in one direction -> shift in another - // Could produce warning + result[i] = func(vval, fval); } - result[i] = func(vval, fval); } } else { // Both must shift to field aligned + // (even if one of v and f has yup/ydown fields, it doesn't make sense to + // multiply them with one in field-aligned and one in non-field-aligned + // coordinates) Field3D v_fa = this->toFieldAligned(v); Field3D f_fa = this->toFieldAligned(f); @@ -2373,7 +2352,7 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo vval.m = v_fa[i.ym()]; vval.pp = v_fa[i.offset(0, 2, 0)]; vval.mm = v_fa[i.offset(0, -2, 0)]; - fval.c = f[i]; + fval.c = f_fa[i]; fval.p = f_fa[i.yp()]; fval.m = f_fa[i.ym()]; fval.pp = f_fa[i.offset(0, 2, 0)]; @@ -2395,6 +2374,8 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo } result[i] = func(vval, fval); } + + result = this->fromFieldAligned(result, RGN_NOBNDRY); } } else { // Non-staggered case @@ -2411,23 +2392,37 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo // f has yup and ydown fields which are distinct stencil fs; - fs.pp = nan(""); - fs.mm = nan(""); + if (mesh->ystart==1) { + // Only one guard cell, so can only use 3-point stencils + fs.pp = nan(""); + fs.mm = nan(""); - Field3D f_yup = f.yup(); - Field3D f_ydown = f.ydown(); + for (const auto &i : result.region(region)) { - for (const auto &i : result.region(region)) { + fs.c = f[i]; + fs.p = f.yup()[i.yp()]; + fs.m = f.ydown()[i.ym()]; - fs.c = f[i]; - fs.p = f_yup[i.yp()]; - fs.m = f_ydown[i.ym()]; + result[i] = func(v[i], fs); + } + } else { + // Can use 5-point stencils + for (const auto &i : result.region(region)) { - result[i] = func(v[i], fs); - } + fs.c = f[i]; + fs.p = f.yup()[i.yp()]; + fs.m = f.ydown()[i.ym()]; + fs.pp = f.yup(2)[i.offset(0, 2, 0)]; + fs.mm = f.ydown(2)[i.offset(0, -2, 0)]; + result[i] = func(v[i], fs); + } + } } else { // Not using yup/ydown fields, so first transform to field-aligned coordinates + // (even if one of v and f has yup/ydown fields, it doesn't make sense to + // multiply them with one in field-aligned and one in non-field-aligned + // coordinates) Field3D f_fa = this->toFieldAligned(f); Field3D v_fa = this->toFieldAligned(v); @@ -2458,15 +2453,14 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo } } // Shift result back - result = this->fromFieldAligned(result); + result = this->fromFieldAligned(result, RGN_NOBNDRY); } } - result.setLocation(diffloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; + invalidateGuards(result); // extra check: set guard cells to NaN if CHECK>2 #endif return result; @@ -2484,6 +2478,7 @@ const Field3D Mesh::indexVDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); CELL_LOC vloc = v.getLocation(); CELL_LOC inloc = f.getLocation(); // Input location @@ -2575,8 +2570,6 @@ const Field3D Mesh::indexVDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo } } - result.setLocation(diffloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -2593,6 +2586,8 @@ const Field2D Mesh::indexFDDX(const Field2D &v, const Field2D &f, CELL_LOC outlo DIFF_METHOD method, REGION region) { TRACE("Mesh::::indexFDDX(Field2D, Field2D)"); + CELL_LOC diffloc = f.getLocation(); + if ((method == DIFF_SPLIT) || ((method == DIFF_DEFAULT) && (fFDDX == nullptr))) { // Split into an upwind and a central differencing part // d/dx(v*f) = v*d/dx(f) + f*d/dx(v) @@ -2724,6 +2719,7 @@ const Field3D Mesh::indexFDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(diffloc); if (this->xstart > 1) { // Two or more guard cells @@ -2860,8 +2856,6 @@ const Field3D Mesh::indexFDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo throw BoutException("Error: Derivatives in X requires at least one guard cell"); } - result.setLocation(diffloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -2895,6 +2889,7 @@ const Field2D Mesh::indexFDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo Field2D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (StaggerGrids && ((v.getLocation() != CELL_CENTRE) || (f.getLocation() != CELL_CENTRE))) { @@ -2950,8 +2945,6 @@ const Field2D Mesh::indexFDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo throw BoutException("Error: Derivatives in Y requires at least one guard cell"); } - result.setLocation(diffloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = false; @@ -3014,131 +3007,88 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(diffloc); - // There are four cases, corresponding to whether or not f and v - // have yup, ydown fields. - - // If vUseUpDown is true, field "v" has distinct yup and ydown fields which - // will be used to calculate a derivative along - // the magnetic field + // If *UseUpDown is true, field "*" 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))); if (vUseUpDown && fUseUpDown) { // Both v and f have up/down fields stencil vval, fval; - vval.mm = nan(""); - vval.pp = nan(""); - fval.mm = nan(""); - fval.pp = nan(""); - for (const auto &i : result.region(region)) { + if (mesh->ystart==1) { + // Only one guard cell, so can only use 3-point stencils + vval.mm = nan(""); + vval.pp = nan(""); + fval.mm = nan(""); + fval.pp = nan(""); + for (const auto &i : result.region(region)) { - fval.m = f.ydown()[i.ym()]; - fval.c = f[i]; - fval.p = f.yup()[i.yp()]; + fval.m = f.ydown()[i.ym()]; + fval.c = f[i]; + fval.p = f.yup()[i.yp()]; - vval.m = v.ydown()[i.ym()]; - vval.c = v[i]; - vval.p = v.yup()[i.yp()]; + vval.m = v.ydown()[i.ym()]; + vval.c = v[i]; + vval.p = v.yup()[i.yp()]; - if(StaggerGrids && (diffloc != CELL_DEFAULT) && (diffloc != vloc)) { - // Non-centred stencil - if((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) { - // Producing a stencil centred around a lower Y value - vval.pp = vval.p; - vval.p = vval.c; - }else if(vloc == CELL_YLOW) { - // Stencil centred around a cell centre - vval.mm = vval.m; - vval.m = vval.c; + if(StaggerGrids && (diffloc != CELL_DEFAULT) && (diffloc != vloc)) { + // Non-centred stencil + if((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) { + // Producing a stencil centred around a lower Y value + vval.pp = vval.p; + vval.p = vval.c; + }else if(vloc == CELL_YLOW) { + // Stencil centred around a cell centre + vval.mm = vval.m; + vval.m = vval.c; + } + // Shifted in one direction -> shift in another + // Could produce warning } - // Shifted in one direction -> shift in another - // Could produce warning + result[i] = func(vval, fval); } - result[i] = func(vval, fval); - } - } - else if (vUseUpDown) { - // Only v has up/down fields - // f must shift to field aligned coordinates - Field3D f_fa = this->toFieldAligned(f); - - stencil vval; - vval.mm = nan(""); - vval.pp = nan(""); - - stencil fval; - for (const auto &i : result.region(region)) { + } else { + // Can use 5-point stencils + for (const auto &i : result.region(region)) { - fval.mm = f_fa[i.offset(0, -2, 0)]; - fval.m = f_fa[i.ym()]; - fval.c = f_fa[i]; - fval.p = f_fa[i.yp()]; - fval.pp = f_fa[i.offset(0, 2, 0)]; + fval.mm = f.ydown(2)[i.offset(0, -2, 0)]; + fval.m = f.ydown()[i.ym()]; + fval.c = f[i]; + fval.p = f.yup()[i.yp()]; + fval.pp = f.yup(2)[i.offset(0, 2, 0)]; - vval.m = v.ydown()[i.ym()]; - vval.c = v[i]; - vval.p = v.yup()[i.yp()]; + vval.mm = v.ydown(2)[i.offset(0, -2, 0)]; + vval.m = v.ydown()[i.ym()]; + vval.c = v[i]; + vval.p = v.yup()[i.yp()]; + vval.pp = v.yup(2)[i.offset(0, 2, 0)]; - if(StaggerGrids && (diffloc != CELL_DEFAULT) && (diffloc != vloc)) { - // Non-centred stencil - if((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) { - // Producing a stencil centred around a lower Y value - vval.pp = vval.p; - vval.p = vval.c; - }else if(vloc == CELL_YLOW) { - // Stencil centred around a cell centre - vval.mm = vval.m; - vval.m = vval.c; + if(StaggerGrids && (diffloc != CELL_DEFAULT) && (diffloc != vloc)) { + // Non-centred stencil + if((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) { + // Producing a stencil centred around a lower Y value + vval.pp = vval.p; + vval.p = vval.c; + }else if(vloc == CELL_YLOW) { + // Stencil centred around a cell centre + vval.mm = vval.m; + vval.m = vval.c; + } + // Shifted in one direction -> shift in another + // Could produce warning } - // Shifted in one direction -> shift in another - // Could produce warning + result[i] = func(vval, fval); } - result[i] = func(vval, fval); - } - } - else if (fUseUpDown) { - // Only f has up/down fields - // v must shift to field aligned coordinates - Field3D v_fa = this->toFieldAligned(v); - - stencil vval; - - stencil fval; - fval.pp = nan(""); - fval.mm = nan(""); - - for (const auto &i : result.region(region)) { - fval.m = f.ydown()[i.ym()]; - fval.c = f[i]; - fval.p = f.yup()[i.yp()]; - - vval.mm = v_fa[i.offset(0,-2,0)]; - vval.m = v_fa[i.ym()]; - vval.c = v_fa[i]; - vval.p = v_fa[i.yp()]; - vval.pp = v_fa[i.offset(0,2,0)]; - - if(StaggerGrids && (diffloc != CELL_DEFAULT) && (diffloc != vloc)) { - // Non-centred stencil - if((vloc == CELL_CENTRE) && (diffloc == CELL_YLOW)) { - // Producing a stencil centred around a lower Y value - vval.pp = vval.p; - vval.p = vval.c; - }else if(vloc == CELL_YLOW) { - // Stencil centred around a cell centre - vval.mm = vval.m; - vval.m = vval.c; - } - // Shifted in one direction -> shift in another - // Could produce warning - } - result[i] = func(vval, fval); } } else { // Both must shift to field aligned + // (even if one of v and f has yup/ydown fields, it doesn't make sense to + // multiply them with one in field-aligned and one in non-field-aligned + // coordinates) Field3D v_fa = this->toFieldAligned(v); Field3D f_fa = this->toFieldAligned(f); @@ -3174,9 +3124,9 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo } result[i] = func(vval, fval); } - } - result.setLocation(diffloc); + result = this->fromFieldAligned(result); + } #if CHECK > 0 // Mark boundaries as invalid @@ -3238,6 +3188,7 @@ const Field3D Mesh::indexFDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(diffloc); stencil vval, fval; for (const auto &i : result.region(region)) { @@ -3274,8 +3225,6 @@ const Field3D Mesh::indexFDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo result[i] = func(vval, fval); } - result.setLocation(diffloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; diff --git a/src/mesh/interpolation.cxx b/src/mesh/interpolation.cxx index 9ade01b1ab..684316c6e8 100644 --- a/src/mesh/interpolation.cxx +++ b/src/mesh/interpolation.cxx @@ -56,6 +56,7 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) { Mesh *fieldmesh = var.getMesh(); Field3D result(fieldmesh); + result.setLocation(loc); if ((loc != CELL_CENTRE && loc != CELL_DEFAULT) && (fieldmesh->StaggerGrids == false)) { throw BoutException("Asked to interpolate, but StaggerGrids is disabled!"); @@ -120,37 +121,59 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) { // Field "var" has distinct yup and ydown fields which // will be used to calculate a derivative along // the magnetic field - throw BoutException("At the moment, fields with yup/ydown cannot use interp_to.\n" - "If we implement a 3-point stencil for interpolate or double-up\n" - "/double-down fields, then we can use this case."); - s.pp = nan(""); - s.mm = nan(""); + // More than one guard cell, so set pp and mm values + // This allows higher-order methods to be used + if (fieldmesh->ystart > 1) { + for(const auto &i : result.region(RGN_NOY)) { + // Set stencils + s.c = var[i]; + s.p = var.yup()[i.yp()]; + s.m = var.ydown()[i.ym()]; + s.pp = var.yup(2)[i.offset(0,2,0)]; + s.mm = var.ydown(2)[i.offset(0,-2,0)]; - for (const auto &i : result.region(RGN_NOBNDRY)) { - // Set stencils - s.c = var[i]; - s.p = var.yup()[i.yp()]; - s.m = var.ydown()[i.ym()]; + if (location == CELL_CENTRE) { + // Producing a stencil centred around a lower Y value + s.pp = s.p; + s.p = s.c; + } else { + // Stencil centred around a cell centre + s.mm = s.m; + s.m = s.c; + } - if ((location == CELL_CENTRE) && (loc == CELL_YLOW)) { - // Producing a stencil centred around a lower Y value - s.pp = s.p; - s.p = s.c; - } else if (location == CELL_YLOW) { - // Stencil centred around a cell centre - s.mm = s.m; - s.m = s.c; + result[i] = interp(s); } + } else { + // Note: at the moment we cannot reach this case because of the + // 'ASSERT0(mesh->ystart >=2)' above, but if we implement a 3-point + // stencil for interp, then this will be useful + s.pp = nan(""); + s.mm = nan(""); + for(const auto &i : result.region(RGN_NOY)) { + // Set stencils + s.c = var[i]; + s.p = var.yup()[i.yp()]; + s.m = var.ydown()[i.ym()]; - result[i] = interp(s); + if (location == CELL_CENTRE) { + // Producing a stencil centred around a lower Y value + s.pp = s.p; + s.p = s.c; + } else { + // Stencil centred around a cell centre + s.mm = s.m; + s.m = s.c; + } + + result[i] = interp(s); + } } } else { // var has no yup/ydown fields, so we need to shift into field-aligned // coordinates Field3D var_fa = fieldmesh->toFieldAligned(var); - Field3D result_fa; - result_fa.allocate(); if (fieldmesh->ystart > 1) { // More than one guard cell, so set pp and mm values @@ -173,7 +196,7 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) { s.m = s.c; } - result_fa[i] = interp(s); + result[i] = interp(s); } } else { // Only one guard cell, so no pp or mm values @@ -198,11 +221,11 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) { s.m = s.c; } - result_fa[i] = interp(s); + result[i] = interp(s); } } - result = fieldmesh->fromFieldAligned(result_fa); + result = fieldmesh->fromFieldAligned(result, RGN_NOBNDRY); } break; } @@ -235,6 +258,8 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) { } }; + invalidateGuards(result); // Fill guard cells with NaN so we can check they are not used when unset. + if ((dir != CELL_ZLOW) && (region != RGN_NOBNDRY)) { fieldmesh->communicate(result); } diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx index 05bd0093d6..915356b50e 100644 --- a/src/mesh/parallel/fci.hxx +++ b/src/mesh/parallel/fci.hxx @@ -69,20 +69,34 @@ public: class FCITransform : public ParallelTransform { public: DEPRECATED(FCITransform(Mesh &mesh, bool UNUSED(yperiodic), bool zperiodic)) - : FCITransform(mesh, zperiodic) {} + : FCITransform(mesh, zperiodic) { + if (mesh.ystart > 1) + // FCITransform can only use myg=1 because it only loads grid + // information for one point forward or back along the magnetic + // field, so it cannot set Field3D::yup2_field or + // Field3D::ydown2_field + throw BoutException("FCI method must use only one y-guard cell: set option myg=1"); + } FCITransform(Mesh &mesh, bool zperiodic = true) : mesh(mesh), forward_map(mesh, +1, zperiodic), backward_map(mesh, -1, zperiodic), - zperiodic(zperiodic) {} + zperiodic(zperiodic) { + if (mesh.ystart > 1) + // FCITransform can only use myg=1 because it only loads grid + // information for one point forward or back along the magnetic + // field, so it cannot set Field3D::yup2_field or + // Field3D::ydown2_field + throw BoutException("FCI method must use only one y-guard cell: set option myg=1"); + } void calcYUpDown(Field3D &f) override; void integrateYUpDown(Field3D &f) override; - const Field3D toFieldAligned(const Field3D &UNUSED(f)) override { + const Field3D toFieldAligned(const Field3D &UNUSED(f), const REGION UNUSED(region)) override { throw BoutException("FCI method cannot transform into field aligned grid"); } - const Field3D fromFieldAligned(const Field3D &UNUSED(f)) override { + const Field3D fromFieldAligned(const Field3D &UNUSED(f), const REGION UNUSED(region)) override { throw BoutException("FCI method cannot transform into field aligned grid"); } diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx index 13067fe07c..cb998d38c7 100644 --- a/src/mesh/parallel/shiftedmetric.cxx +++ b/src/mesh/parallel/shiftedmetric.cxx @@ -8,6 +8,7 @@ #include #include +#include #include #include @@ -15,7 +16,14 @@ #include -ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) { +ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m), + has_toAligned_CENTRE(false), has_toAligned_XLOW(false), has_toAligned_YLOW(false), + has_fromAligned_CENTRE(false), has_fromAligned_XLOW(false), has_fromAligned_YLOW(false), + has_yupPhs1_CENTRE(false), has_yupPhs1_XLOW(false), has_yupPhs1_YLOW(false), + has_ydownPhs1_CENTRE(false), has_ydownPhs1_XLOW(false), has_ydownPhs1_YLOW(false), + has_yupPhs2_CENTRE(false), has_yupPhs2_XLOW(false), has_yupPhs2_YLOW(false), + has_ydownPhs2_CENTRE(false), has_ydownPhs2_XLOW(false), has_ydownPhs2_YLOW(false) +{ // Read the zShift angle from the mesh if(mesh.get(zShift, "zShift")) { @@ -23,95 +31,670 @@ ShiftedMetric::ShiftedMetric(Mesh &m) : mesh(m), zShift(&m) { mesh.get(zShift, "qinty"); } - //If we wanted to be efficient we could move the following cached phase setup - //into the relevant shifting routines (with static bool first protection) - //so that we only calculate the phase if we actually call a relevant shift - //routine -- however as we're only going to do this initialisation once I - //think it's cleaner to put it in the constructor here. + if (mesh.StaggerGrids) { + if (mesh.xstart >=2) { + // Can interpolate in x-direction + // Calculate staggered field for zShift and apply boundary conditions + Field2D zShift_XLOW = interp_to(zShift, CELL_XLOW, RGN_ALL); + zShift_XLOW.applyBoundary("neumann"); // Set boundary guard cells to closest grid cell value + zShift.set(zShift_XLOW); + } + if (mesh.ystart >=2) { + // Can interpolate in y-direction + // Calculate staggered field for zShift and apply boundary conditions + Field2D zShift_YLOW = interp_to(zShift, CELL_YLOW, RGN_ALL); + zShift_YLOW.applyBoundary("neumann"); // Set boundary guard cells to closest grid cell value + zShift.set(zShift_YLOW); + } + } - //As we're attached to a mesh we can expect the z direction to - //not change once we've been created so precalculate the complex - //phases used in transformations int nmodes = mesh.LocalNz/2 + 1; - BoutReal zlength = mesh.coordinates()->zlength(); - //Allocate storage for complex intermediate - cmplx.resize(nmodes); + cmplx = Array(nmodes); std::fill(cmplx.begin(), cmplx.end(), 0.0); +} - //Allocate storage for our 3d vector structures. - //This could be made more succinct but this approach is fairly - //verbose --> transparent - fromAlignedPhs.resize(mesh.LocalNx); - toAlignedPhs.resize(mesh.LocalNx); - - yupPhs.resize(mesh.LocalNx); - ydownPhs.resize(mesh.LocalNx); +//As we're attached to a mesh we can expect the z direction to not change +//once we've been created so cache the complex phases used in transformations +//the first time they are needed +Matrix< Array > ShiftedMetric::getFromAlignedPhs(CELL_LOC location) { + switch (location) { + case CELL_CENTRE: { + if (!has_fromAligned_CENTRE) { + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + + has_fromAligned_CENTRE = true; + fromAlignedPhs_CENTRE = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : fromAlignedPhs_CENTRE) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=0;jy=2); //otherwise we cannot interpolate in the x-direction + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + + // get zShift at XLOW + Field2D zShift_XLOW = zShift.get(CELL_XLOW); + + has_fromAligned_XLOW = true; + fromAlignedPhs_XLOW = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : fromAlignedPhs_XLOW) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=0;jy=2); //otherwise we cannot interpolate in the y-direction + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + + // get zShift at YLOW + Field2D zShift_YLOW = zShift.get(CELL_YLOW); + + has_fromAligned_YLOW = true; + fromAlignedPhs_YLOW = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : fromAlignedPhs_YLOW) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=0;jy > ShiftedMetric::getToAlignedPhs(CELL_LOC location) { + switch (location) { + case CELL_CENTRE: { + if (!has_toAligned_CENTRE) { + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + + has_toAligned_CENTRE = true; + toAlignedPhs_CENTRE = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : toAlignedPhs_CENTRE) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=0;jy=2); //otherwise we cannot interpolate in the x-direction + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + + // get zShift at XLOW + Field2D zShift_XLOW = zShift.get(CELL_XLOW); + + has_toAligned_XLOW = true; + toAlignedPhs_XLOW = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : toAlignedPhs_XLOW) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=0;jy=2); //otherwise we cannot interpolate in the y-direction + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + + // get zShift at YLOW + Field2D zShift_YLOW = zShift.get(CELL_YLOW); + + has_toAligned_YLOW = true; + toAlignedPhs_YLOW = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : toAlignedPhs_YLOW) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=0;jy > ShiftedMetric::getYupPhs1(CELL_LOC location) { + switch (location) { + case CELL_CENTRE: { + if (!has_yupPhs1_CENTRE) { + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + + has_yupPhs1_CENTRE = true; + yupPhs1_CENTRE = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : yupPhs1_CENTRE) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=mesh.ystart; jy<=mesh.yend; jy++) { + BoutReal yupShift1 = zShift(jx,jy) - zShift(jx,jy+1); + for(int jz=0;jz=2); //otherwise we cannot interpolate in the x-direction + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + + // get zShift at XLOW + Field2D zShift_XLOW = zShift.get(CELL_XLOW); + + has_yupPhs1_XLOW = true; + yupPhs1_XLOW = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : yupPhs1_XLOW) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=mesh.ystart; jy<=mesh.yend; jy++) { + BoutReal yupShift1 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy+1); + for(int jz=0;jz=2); //otherwise we cannot interpolate in the y-direction + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + + // get zShift at YLOW + Field2D zShift_YLOW = zShift.get(CELL_YLOW); + + has_yupPhs1_YLOW = true; + yupPhs1_YLOW = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : yupPhs1_YLOW) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=mesh.ystart; jy<=mesh.yend; jy++) { + BoutReal yupShift1 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy+1); + for(int jz=0;jz > ShiftedMetric::getYupPhs2(CELL_LOC location) { + switch (location) { + case CELL_CENTRE: { + if (!has_yupPhs2_CENTRE) { + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + + has_yupPhs2_CENTRE = true; + yupPhs2_CENTRE = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : yupPhs2_CENTRE) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=mesh.ystart; jy<=mesh.yend; jy++) { + BoutReal yupShift2 = zShift(jx,jy) - zShift(jx,jy+2); + for(int jz=0;jz=2); //otherwise we cannot interpolate in the x-direction + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + + // get zShift at XLOW + Field2D zShift_XLOW = zShift.get(CELL_XLOW); + + has_yupPhs2_XLOW = true; + yupPhs2_XLOW = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : yupPhs2_XLOW) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=mesh.ystart; jy<=mesh.yend; jy++) { + BoutReal yupShift2 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy+2); + for(int jz=0;jz=2); //otherwise we cannot interpolate in the y-direction + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + + // get zShift at YLOW + Field2D zShift_YLOW = zShift.get(CELL_YLOW); + + has_yupPhs2_YLOW = true; + yupPhs2_YLOW = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : yupPhs2_YLOW) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=mesh.ystart; jy<=mesh.yend; jy++) { + BoutReal yupShift2 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy+2); + for(int jz=0;jz > ShiftedMetric::getYdownPhs1(CELL_LOC location) { + switch (location) { + case CELL_CENTRE: { + if (has_ydownPhs1_CENTRE) { + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + + has_ydownPhs1_CENTRE = true; + ydownPhs1_CENTRE = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : ydownPhs1_CENTRE) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=mesh.ystart; jy<=mesh.yend; jy++) { + BoutReal ydownShift1 = zShift(jx,jy) - zShift(jx,jy-1); + for(int jz=0;jz=2); //otherwise we cannot interpolate in the x-direction + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + + // get zShift at XLOW + Field2D zShift_XLOW = zShift.get(CELL_XLOW); + + has_ydownPhs1_XLOW = true; + ydownPhs1_XLOW = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : ydownPhs1_XLOW) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=mesh.ystart; jy<=mesh.yend; jy++) { + BoutReal ydownShift1 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy-1); + for(int jz=0;jz=2); //otherwise we cannot interpolate in the y-direction + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); - for(int jx=0;jx >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : ydownPhs1_YLOW) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=mesh.ystart; jy<=mesh.yend; jy++) { + BoutReal ydownShift1 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy-1); + for(int jz=0;jz > ShiftedMetric::getYdownPhs2(CELL_LOC location) { + switch (location) { + case CELL_CENTRE: { + if (!has_ydownPhs2_CENTRE) { + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + + has_ydownPhs2_CENTRE = true; + ydownPhs2_CENTRE = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : ydownPhs2_CENTRE) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=mesh.ystart; jy<=mesh.yend; jy++) { + BoutReal ydownShift2 = zShift(jx,jy) - zShift(jx,jy-2); + for(int jz=0;jz=2); //otherwise we cannot interpolate in the x-direction + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); - //Yup/Ydown phases -- note we don't shift in the boundaries/guards - for(int jx=0;jx >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : ydownPhs2_XLOW) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=mesh.ystart; jy<=mesh.yend; jy++) { + BoutReal ydownShift2 = zShift_XLOW(jx,jy) - zShift_XLOW(jx,jy-2); + for(int jz=0;jz=2); //otherwise we cannot interpolate in the y-direction + int nmodes = mesh.LocalNz/2 + 1; + BoutReal zlength = mesh.coordinates()->zlength(); + // get zShift at YLOW + Field2D zShift_YLOW = zShift.get(CELL_YLOW); + + has_ydownPhs2_YLOW = true; + ydownPhs2_YLOW = Matrix< Array >(mesh.LocalNx, mesh.LocalNy); + for (auto &element : ydownPhs2_YLOW) { + element = Array(mesh.LocalNz); + } + + //To/From field aligned phases + for(int jx=mesh.xstart; jx<=mesh.xend; jx++) { + for(int jy=mesh.ystart; jy<=mesh.yend; jy++) { + BoutReal ydownShift2 = zShift_YLOW(jx,jy) - zShift_YLOW(jx,jy-2); + for(int jz=0;jz > phases = getYupPhs1(location); + for(auto i : f.region2D(RGN_NOBNDRY)) { + shiftZ(f(i.x, i.y+1), phases(i.x, i.y), yup1(i.x, i.y+1)); + } + if (mesh.ystart>1) { + Field3D& yup2 = f.yup(2); + yup2.allocate(); + invalidateGuards(yup2); // Won't set x-guard cells, so allow checking to throw exception if they are used. + phases = getYupPhs2(location); + for(auto i : f.region2D(RGN_NOBNDRY)) { + shiftZ(f(i.x, i.y+2), phases(i.x, i.y), yup2(i.x, i.y+2)); } } - Field3D& ydown = f.ydown(); - ydown.allocate(); - - for(int jx=0;jx 1) { + Field3D& ydown2 = f.ydown(2); + ydown2.allocate(); + invalidateGuards(ydown2); // Won't set x-guard cells, so allow checking to throw exception if they are used. + phases = getYdownPhs2(location); + for(auto i : f.region2D(RGN_NOBNDRY)) { + shiftZ(f(i.x, i.y-2), phases(i.x, i.y), ydown2(i.x, i.y-2)); } } } @@ -120,39 +703,37 @@ void ShiftedMetric::calcYUpDown(Field3D &f) { * Shift the field so that X-Z is not orthogonal, * and Y is then field aligned. */ -const Field3D ShiftedMetric::toFieldAligned(const Field3D &f) { - return shiftZ(f, toAlignedPhs); +const Field3D ShiftedMetric::toFieldAligned(const Field3D &f, const REGION region) { + return shiftZ(f, getToAlignedPhs(f.getLocation()), region); } /*! * Shift back, so that X-Z is orthogonal, * but Y is not field aligned. */ -const Field3D ShiftedMetric::fromFieldAligned(const Field3D &f) { - return shiftZ(f, fromAlignedPhs); +const Field3D ShiftedMetric::fromFieldAligned(const Field3D &f, const REGION region) { + return shiftZ(f, getFromAlignedPhs(f.getLocation()), region); } -const Field3D ShiftedMetric::shiftZ(const Field3D &f, const arr3Dvec &phs) { +const Field3D ShiftedMetric::shiftZ(const Field3D &f, const Matrix< Array > &phs, const REGION region) { ASSERT1(&mesh == f.getMesh()); + ASSERT1(region == RGN_NOX || region == RGN_NOBNDRY); // Never calculate x-guard cells here if(mesh.LocalNz == 1) return f; // Shifting makes no difference - Field3D result(&mesh); - result.allocate(); - - for(int jx=0;jx &phs, BoutReal *out) { +void ShiftedMetric::shiftZ(const BoutReal *in, const Array &phs, BoutReal *out) { // Take forward FFT - rfft(in, mesh.LocalNz, &cmplx[0]); + rfft(in, mesh.LocalNz, cmplx.begin()); //Following is an algorithm approach to write a = a*b where a and b are //vectors of dcomplex. @@ -164,22 +745,29 @@ void ShiftedMetric::shiftZ(const BoutReal *in, const std::vector &phs, cmplx[jz] *= phs[jz]; } - irfft(&cmplx[0], mesh.LocalNz, out); // Reverse FFT + irfft(cmplx.begin(), mesh.LocalNz, out); // Reverse FFT } //Old approach retained so we can still specify a general zShift -const Field3D ShiftedMetric::shiftZ(const Field3D &f, const Field2D &zangle) { +const Field3D ShiftedMetric::shiftZ(const Field3D &f, const Field2D &zangle, const REGION region) { ASSERT1(&mesh == f.getMesh()); + ASSERT1(region == RGN_NOX || region == RGN_NOBNDRY); // Never calculate x-guard cells here + ASSERT1(f.getLocation() == zangle.getLocation()); if(mesh.LocalNz == 1) return f; // Shifting makes no difference Field3D result(&mesh); result.allocate(); + invalidateGuards(result); // Won't set x-guard cells, so allow checking to throw exception if they are used. - for(int jx=0;jx(nmodes); // Take forward FFT - rfft(in, len, &cmplxLoc[0]); + rfft(in, len, cmplxLoc.begin()); // Apply phase shift BoutReal zlength = mesh.coordinates()->zlength(); @@ -201,5 +789,9 @@ void ShiftedMetric::shiftZ(const BoutReal *in, int len, BoutReal zangle, BoutRe cmplxLoc[jz] *= dcomplex(cos(kwave*zangle) , -sin(kwave*zangle)); } - irfft(&cmplxLoc[0], len, out); // Reverse FFT + irfft(cmplxLoc.begin(), len, out); // Reverse FFT +} + +void ShiftedMetric::outputVars(Datafile &file) { + file.add(zShift, "zShift", 0); } diff --git a/tests/integrated/test-fci-slab/data/BOUT.inp b/tests/integrated/test-fci-slab/data/BOUT.inp index 218946e676..dc0ff71b81 100644 --- a/tests/integrated/test-fci-slab/data/BOUT.inp +++ b/tests/integrated/test-fci-slab/data/BOUT.inp @@ -6,6 +6,8 @@ timestep = 0.2 MZ = 64 +myg = 1 + [mesh] paralleltransform = fci diff --git a/tests/integrated/test-fci-slab/mms/BOUT.inp b/tests/integrated/test-fci-slab/mms/BOUT.inp index 3c1e5f9559..5280f748be 100644 --- a/tests/integrated/test-fci-slab/mms/BOUT.inp +++ b/tests/integrated/test-fci-slab/mms/BOUT.inp @@ -7,6 +7,8 @@ MZ = 64 NXPE = 1 +myg = 1 + [mesh] paralleltransform = fci diff --git a/tests/integrated/test-yupdown/test_yupdown.cxx b/tests/integrated/test-yupdown/test_yupdown.cxx index e414c0ba88..f7c9e02b79 100644 --- a/tests/integrated/test-yupdown/test_yupdown.cxx +++ b/tests/integrated/test-yupdown/test_yupdown.cxx @@ -37,8 +37,6 @@ int main(int argc, char** argv) { BoutInitialise(argc, argv); - ShiftedMetric s(*mesh); - // Read variable from mesh Field3D var; mesh->get(var, "var"); @@ -46,7 +44,7 @@ int main(int argc, char** argv) { // Var starts in orthogonal X-Z coordinates // Calculate yup and ydown - s.calcYUpDown(var); + mesh->communicate(var); // Calculate d/dy ysing yup() and ydown() fields Field3D ddy = DDY_yud(var); diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index 2fb7e612c5..bd9a1b4872 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -12,8 +12,9 @@ #include #include -/// Global mesh +/// Global meshes extern Mesh *mesh; +Mesh *mesh2; // mesh2 has 2 guard cells /// Test fixture to make sure the global mesh is our fake one class Field3DTest : public ::testing::Test { @@ -26,11 +27,24 @@ class Field3DTest : public ::testing::Test { } mesh = new FakeMesh(nx, ny, nz); mesh->createDefaultRegions(); + + if (mesh2 != nullptr) { + delete mesh2; + mesh2 = nullptr; + } + mesh2 = new FakeMesh(nx, ny, nz); + mesh2->createDefaultRegions(); + mesh2->xstart += 1; + mesh2->xend -= 1; + mesh2->ystart += 1; + mesh2->yend -=1; } static void TearDownTestCase() { delete mesh; mesh = nullptr; + delete mesh2; + mesh2 = nullptr; } public: @@ -275,6 +289,34 @@ TEST_F(Field3DTest, SplitThenMergeYupYDown) { EXPECT_EQ(&field, &ydown2); } +TEST_F(Field3DTest, SplitThenMergeYupYDown2) { + Field3D field(mesh2); + + field = 0.; + field.splitYupYdown(); + + auto& yup1 = field.yup(); + EXPECT_NE(&field, &yup1); + auto& ydown1 = field.ydown(); + EXPECT_NE(&field, &ydown1); + auto& yup2 = field.yup(2); + EXPECT_NE(&field, &yup2); + auto& ydown2 = field.ydown(2); + EXPECT_NE(&field, &ydown2); + + field.mergeYupYdown(); + + auto& yup1_2 = field.yup(); + EXPECT_EQ(&field, &yup1_2); + auto& ydown1_2 = field.ydown(); + EXPECT_EQ(&field, &ydown1_2); + + auto& yup2_2 = field.yup(2); + EXPECT_EQ(&field, &yup2_2); + auto& ydown2_2 = field.ydown(2); + EXPECT_EQ(&field, &ydown2_2); +} + TEST_F(Field3DTest, Ynext) { Field3D field; @@ -287,7 +329,9 @@ TEST_F(Field3DTest, Ynext) { EXPECT_NE(&field, &ydown); EXPECT_NE(&yup, &ydown); +#if CHECK > 1 EXPECT_THROW(field.ynext(99), BoutException); +#endif } TEST_F(Field3DTest, ConstYnext) { @@ -303,7 +347,56 @@ TEST_F(Field3DTest, ConstYnext) { EXPECT_NE(&field2, &ydown); EXPECT_NE(&yup, &ydown); +#if CHECK > 1 EXPECT_THROW(field2.ynext(99), BoutException); +#endif +} + +TEST_F(Field3DTest, Ynext2) { + Field3D field(mesh2); + + field = 0.; + field.splitYupYdown(); + + auto& yup = field.ynext(1); + EXPECT_NE(&field, &yup); + auto& ydown = field.ynext(-1); + EXPECT_NE(&field, &ydown); + EXPECT_NE(&yup, &ydown); + + auto& yup2 = field.ynext(2); + EXPECT_NE(&field, &yup2); + auto& ydown2 = field.ynext(-2); + EXPECT_NE(&field, &ydown2); + EXPECT_NE(&yup2, &ydown2); + +#if CHECK > 1 + EXPECT_THROW(field.ynext(99), BoutException); +#endif +} + +TEST_F(Field3DTest, ConstYnext2) { + Field3D field(0., mesh2); + + field.splitYupYdown(); + + const Field3D& field2 = field; + + auto& yup = field2.ynext(1); + EXPECT_NE(&field2, &yup); + auto& ydown = field2.ynext(-1); + EXPECT_NE(&field2, &ydown); + EXPECT_NE(&yup, &ydown); + + auto& yup2 = field2.ynext(2); + EXPECT_NE(&field2, &yup2); + auto& ydown2 = field2.ynext(-2); + EXPECT_NE(&field2, &ydown2); + EXPECT_NE(&yup2, &ydown2); + +#if CHECK > 1 + EXPECT_THROW(field2.ynext(99), BoutException); +#endif } TEST_F(Field3DTest, GetGlobalMesh) {