From f9f16c3dc1a4301f5d82b314965bded8f5c9d664 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 24 Jul 2018 12:04:18 +0100 Subject: [PATCH 1/7] Record which coordinate system Field3D is on Define enum class COORDINATE_SYSTEM. Add a COORDINATE_SYSTEM member to Field3D and write as an attribute in output files. The default is set by the ParallelTransform but changed to FieldAligned by 'toFieldAligned()' and changed back by 'fromFieldAligned()'. --- include/bout/mesh.hxx | 6 +- include/bout/paralleltransform.hxx | 13 ++ include/field3d.hxx | 34 +++++ src/field/field3d.cxx | 31 +++- src/field/field_factory.cxx | 1 + src/field/gen_fieldops.jinja | 14 +- src/field/generated_fieldops.cxx | 228 +++++++++++++++++----------- src/fileio/datafile.cxx | 6 + src/mesh/index_derivs.cxx | 1 - src/mesh/interpolation.cxx | 1 + src/mesh/parallel/fci.hxx | 4 + src/mesh/parallel/shiftedmetric.cxx | 10 +- 12 files changed, 250 insertions(+), 99 deletions(-) diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index 432cf0fa22..cec5a6eacb 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -695,6 +695,10 @@ class Mesh { return f; } + COORDINATE_SYSTEM getCoordinateSystem() const { + return getParallelTransform().getCoordinateSystem(); + } + bool canToFromFieldAligned() { return getParallelTransform().canToFromFieldAligned(); } @@ -717,7 +721,7 @@ class Mesh { void setParallelTransform(); /*! - * Return the parallel transform, setting it if need be + * Return the parallel transform */ ParallelTransform& getParallelTransform(); diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx index f2f40ed62f..aace91f21f 100644 --- a/include/bout/paralleltransform.hxx +++ b/include/bout/paralleltransform.hxx @@ -43,6 +43,11 @@ public: /// into standard form virtual const Field3D fromFieldAligned(const Field3D &f, const REGION region = RGN_NOX) = 0; + /*! + * Return the default coordinate system for Field3Ds when using this ParallelTransform + */ + virtual COORDINATE_SYSTEM getCoordinateSystem() const = 0; + virtual bool canToFromFieldAligned() = 0; }; @@ -76,6 +81,10 @@ public: return f; } + COORDINATE_SYSTEM getCoordinateSystem() const override { + return COORDINATE_SYSTEM::FieldAligned; + } + bool canToFromFieldAligned() override{ return true; } @@ -116,6 +125,10 @@ public: */ const Field3D fromFieldAligned(const Field3D &f, const REGION region=RGN_NOX) override; + COORDINATE_SYSTEM getCoordinateSystem() const override { + return COORDINATE_SYSTEM::Orthogonal; + } + bool canToFromFieldAligned() override{ return true; } diff --git a/include/field3d.hxx b/include/field3d.hxx index f72d21f502..3b70a4e95e 100644 --- a/include/field3d.hxx +++ b/include/field3d.hxx @@ -40,6 +40,28 @@ class Mesh; // #include "bout/mesh.hxx" #include "bout/field_visitor.hxx" +/*! + * Coordinate systems that Field3Ds can be stored in + * + * None - special value to be used for Field3Ds initialized before the Mesh and + * ParallelTransform have finished being initialized, because then we cannot + * get the COORDINATE_SYSTEM to use from the ParallelTransform, but Field3Ds + * created then must not need to know their COORDINATE_SYSTEM, or must set it + * explicitly later. + */ +enum class COORDINATE_SYSTEM { None, FieldAligned, Orthogonal, FCI }; + +#define COORDENUMSTR(val) {COORDINATE_SYSTEM::val, #val} + +const std::map COORDINATE_SYSTEMtoString = { + COORDENUMSTR(None), + COORDENUMSTR(FieldAligned), + COORDENUMSTR(Orthogonal), + COORDENUMSTR(FCI) +}; + +#undef COORDENUMSTR + /// Class for 3D X-Y-Z scalar fields /*! This class represents a scalar field defined over the mesh. @@ -215,6 +237,16 @@ class Field3D : public Field, public FieldData { */ int getNz() const override {return nz;}; + /*! + * Return the coordinate system of this field + */ + COORDINATE_SYSTEM getCoordinateSystem() const { return coordinate_system; } + + /*! + * Reset the coordinate system of this field + */ + void setCoordinateSystem(COORDINATE_SYSTEM new_coords) { coordinate_system = new_coords; } + /*! * Ensure that this field has separate fields * for yup and ydown. @@ -478,6 +510,8 @@ private: /// Array sizes (from fieldmesh). These are valid only if fieldmesh is not null int nx{-1}, ny{-1}, nz{-1}; + COORDINATE_SYSTEM coordinate_system; ///< Coordinate system used by this field, e.g. fieldaligned, orthogonal + /// Internal data array. Handles allocation/freeing of memory Array data; diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 6ab7a4ef79..6f0650dac7 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -54,11 +54,19 @@ Field3D::Field3D(Mesh* localmesh) : Field(localmesh) { ny = fieldmesh->LocalNy; nz = fieldmesh->LocalNz; } + if (fieldmesh) { + coordinate_system = fieldmesh->getCoordinateSystem(); + } else { + coordinate_system = COORDINATE_SYSTEM::None; + } + } /// Doesn't copy any data, just create a new reference to the same data (copy on change /// later) -Field3D::Field3D(const Field3D& f) : Field(f.fieldmesh), data(f.data) { +Field3D::Field3D(const Field3D& f) + : Field(f.fieldmesh), coordinate_system(f.coordinate_system), data(f.data), + location(f.location) { TRACE("Field3D(Field3D&)"); @@ -72,7 +80,6 @@ Field3D::Field3D(const Field3D& f) : Field(f.fieldmesh), data(f.data) { nz = fieldmesh->LocalNz; } - location = f.location; fieldCoordinates = f.fieldCoordinates; } @@ -80,6 +87,8 @@ Field3D::Field3D(const Field2D& f) : Field(f.getMesh()) { TRACE("Field3D: Copy constructor from Field2D"); + coordinate_system = fieldmesh->getCoordinateSystem(); + nx = fieldmesh->LocalNx; ny = fieldmesh->LocalNy; nz = fieldmesh->LocalNz; @@ -94,6 +103,8 @@ Field3D::Field3D(const BoutReal val, Mesh* localmesh) : Field(localmesh) { TRACE("Field3D: Copy constructor from value"); + coordinate_system = fieldmesh->getCoordinateSystem(); + nx = fieldmesh->LocalNx; ny = fieldmesh->LocalNy; nz = fieldmesh->LocalNz; @@ -275,6 +286,8 @@ Field3D & Field3D::operator=(const Field3D &rhs) { fieldmesh = rhs.fieldmesh; nx = rhs.nx; ny = rhs.ny; nz = rhs.nz; + coordinate_system = rhs.coordinate_system; + data = rhs.data; setLocation(rhs.location); @@ -614,6 +627,7 @@ Field3D pow(const Field3D &lhs, const Field3D &rhs, REGION rgn) { TRACE("pow(Field3D, Field3D)"); ASSERT1(lhs.getLocation() == rhs.getLocation()); + ASSERT1(lhs.getCoordinateSystem() == rhs.getCoordinateSystem()); ASSERT1(lhs.getMesh() == rhs.getMesh()); Field3D result(lhs.getMesh()); @@ -622,6 +636,7 @@ Field3D pow(const Field3D &lhs, const Field3D &rhs, REGION rgn) { BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs[i]); } result.setLocation( lhs.getLocation() ); + result.setCoordinateSystem( lhs.getCoordinateSystem() ); checkData(result); return result; @@ -641,6 +656,7 @@ Field3D pow(const Field3D &lhs, const Field2D &rhs, REGION rgn) { BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs[i]); } result.setLocation( lhs.getLocation() ); + result.setCoordinateSystem( lhs.getCoordinateSystem() ); checkData(result); return result; @@ -653,6 +669,11 @@ FieldPerp pow(const Field3D &lhs, const FieldPerp &rhs, REGION rgn) { checkData(rhs); ASSERT1(lhs.getMesh() == rhs.getMesh()); + // Assume FieldPerp is always on the default coordinate system of its Mesh, + // since at the moment FieldPerp cannot switch between orthogonal and field + // aligned coordinates + ASSERT2(lhs.getCoordinateSystem() == rhs.getMesh()->getCoordinateSystem()); + FieldPerp result{rhs.getMesh()}; result.allocate(); result.setIndex(rhs.getIndex()); @@ -679,6 +700,7 @@ Field3D pow(const Field3D &lhs, BoutReal rhs, REGION rgn) { BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs); } result.setLocation( lhs.getLocation() ); + result.setCoordinateSystem( lhs.getCoordinateSystem() ); checkData(result); return result; @@ -697,6 +719,7 @@ Field3D pow(BoutReal lhs, const Field3D &rhs, REGION rgn) { BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs, rhs[i]); } result.setLocation( rhs.getLocation() ); + result.setCoordinateSystem( rhs.getCoordinateSystem() ); checkData(result); return result; @@ -802,6 +825,7 @@ BoutReal mean(const Field3D &f, bool allpe, REGION rgn) { result.allocate(); \ BOUT_FOR(d, result.getRegion(rgn)) { result[d] = func(f[d]); } \ result.setLocation(f.getLocation()); \ + result.setCoordinateSystem(f.getCoordinateSystem()); \ checkData(result); \ return result; \ } @@ -865,6 +889,7 @@ const Field3D filter(const Field3D &var, int N0, REGION rgn) { #endif result.setLocation(var.getLocation()); + result.setCoordinateSystem(var.getCoordinateSystem()); checkData(result); return result; @@ -912,6 +937,7 @@ const Field3D lowPass(const Field3D &var, int zmax, REGION rgn) { } } result.setLocation(var.getLocation()); + result.setCoordinateSystem(var.getCoordinateSystem()); checkData(result); return result; @@ -962,6 +988,7 @@ const Field3D lowPass(const Field3D &var, int zmax, int zmin, REGION rgn) { } result.setLocation(var.getLocation()); + result.setCoordinateSystem(var.getCoordinateSystem()); checkData(result); return result; diff --git a/src/field/field_factory.cxx b/src/field/field_factory.cxx index bc679e8a30..3daf96054d 100644 --- a/src/field/field_factory.cxx +++ b/src/field/field_factory.cxx @@ -229,6 +229,7 @@ const Field3D FieldFactory::create3D(const std::string &value, const Options *op if (localmesh->canToFromFieldAligned()){ // Ask wheter it is possible // Transform from field aligned coordinates, to be compatible with // older BOUT++ inputs. This is not a particularly "nice" solution. + result.setCoordinateSystem(COORDINATE_SYSTEM::FieldAligned); result = localmesh->fromFieldAligned(result, RGN_ALL); } diff --git a/src/field/gen_fieldops.jinja b/src/field/gen_fieldops.jinja index 58f07d0796..82dc3cc95e 100644 --- a/src/field/gen_fieldops.jinja +++ b/src/field/gen_fieldops.jinja @@ -16,6 +16,10 @@ ASSERT1(localmesh == {{rhs.name}}.getMesh()); {% endif %} + {% if lhs == "Field3D" and rhs == "Field3D" %} + ASSERT1({{lhs.name}}.getCoordinateSystem() == {{rhs.name}}.getCoordinateSystem()); + {% endif %} + {{out.field_type}} {{out.name}}(localmesh); {{out.name}}.allocate(); checkData({{lhs.name}}); @@ -46,8 +50,12 @@ {% if out in ['Field3D','Field2D'] %} {{out.name}}.setLocation({{rhs.name if rhs in ["Field3D","Field2D"] else lhs.name}}.getLocation()); + {% endif %} + {% if out in ['Field3D'] %} + {{out.name}}.setCoordinateSystem({{rhs.name if rhs == "Field3D" else lhs.name}}.getCoordinateSystem()); + {% endif %} checkData({{out.name}}); return {{out.name}}; } @@ -58,7 +66,11 @@ // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { - {% if ((lhs in ["Field3D", "Field2D"]) and (rhs in ["Field3D", "Field2D"])) %} + {% if ((lhs in ["Field3D"]) and (rhs in ["Field3D"])) %} + ASSERT1(this->getCoordinateSystem() == {{rhs.name}}.getCoordinateSystem()) + + {% endif %} + {% if ((lhs in ["Field3D", "Field2D"]) and (rhs in ["Field3D", "Field2D"])) %} #if CHECK > 0 if (this->getLocation() != rhs.getLocation()) { throw BoutException( diff --git a/src/field/generated_fieldops.cxx b/src/field/generated_fieldops.cxx index 0aa552060f..f11e207659 100644 --- a/src/field/generated_fieldops.cxx +++ b/src/field/generated_fieldops.cxx @@ -7,7 +7,7 @@ #include // Provide the C++ wrapper for multiplication of Field3D and Field3D -Field3D operator*(const Field3D &lhs, const Field3D &rhs) { +Field3D operator*(const Field3D& lhs, const Field3D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator*(Field3D, Field3D): fields at different " @@ -16,10 +16,12 @@ Field3D operator*(const Field3D &lhs, const Field3D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); + ASSERT1(lhs.getCoordinateSystem() == rhs.getCoordinateSystem()); + Field3D result(localmesh); result.allocate(); checkData(lhs); @@ -31,15 +33,18 @@ Field3D operator*(const Field3D &lhs, const Field3D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(rhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ operator to update Field3D by multiplication with Field3D -Field3D &Field3D::operator*=(const Field3D &rhs) { +Field3D& Field3D::operator*=(const Field3D& rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { + ASSERT1(this->getCoordinateSystem() == rhs.getCoordinateSystem()) #if CHECK > 0 if (this->getLocation() != rhs.getLocation()) { @@ -66,7 +71,7 @@ Field3D &Field3D::operator*=(const Field3D &rhs) { } // Provide the C++ wrapper for division of Field3D and Field3D -Field3D operator/(const Field3D &lhs, const Field3D &rhs) { +Field3D operator/(const Field3D& lhs, const Field3D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator/(Field3D, Field3D): fields at different " @@ -75,10 +80,12 @@ Field3D operator/(const Field3D &lhs, const Field3D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); + ASSERT1(lhs.getCoordinateSystem() == rhs.getCoordinateSystem()); + Field3D result(localmesh); result.allocate(); checkData(lhs); @@ -90,15 +97,18 @@ Field3D operator/(const Field3D &lhs, const Field3D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(rhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ operator to update Field3D by division with Field3D -Field3D &Field3D::operator/=(const Field3D &rhs) { +Field3D& Field3D::operator/=(const Field3D& rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { + ASSERT1(this->getCoordinateSystem() == rhs.getCoordinateSystem()) #if CHECK > 0 if (this->getLocation() != rhs.getLocation()) { @@ -125,7 +135,7 @@ Field3D &Field3D::operator/=(const Field3D &rhs) { } // Provide the C++ wrapper for addition of Field3D and Field3D -Field3D operator+(const Field3D &lhs, const Field3D &rhs) { +Field3D operator+(const Field3D& lhs, const Field3D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator+(Field3D, Field3D): fields at different " @@ -134,10 +144,12 @@ Field3D operator+(const Field3D &lhs, const Field3D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); + ASSERT1(lhs.getCoordinateSystem() == rhs.getCoordinateSystem()); + Field3D result(localmesh); result.allocate(); checkData(lhs); @@ -149,15 +161,18 @@ Field3D operator+(const Field3D &lhs, const Field3D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(rhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ operator to update Field3D by addition with Field3D -Field3D &Field3D::operator+=(const Field3D &rhs) { +Field3D& Field3D::operator+=(const Field3D& rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { + ASSERT1(this->getCoordinateSystem() == rhs.getCoordinateSystem()) #if CHECK > 0 if (this->getLocation() != rhs.getLocation()) { @@ -184,7 +199,7 @@ Field3D &Field3D::operator+=(const Field3D &rhs) { } // Provide the C++ wrapper for subtraction of Field3D and Field3D -Field3D operator-(const Field3D &lhs, const Field3D &rhs) { +Field3D operator-(const Field3D& lhs, const Field3D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator-(Field3D, Field3D): fields at different " @@ -193,10 +208,12 @@ Field3D operator-(const Field3D &lhs, const Field3D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); + ASSERT1(lhs.getCoordinateSystem() == rhs.getCoordinateSystem()); + Field3D result(localmesh); result.allocate(); checkData(lhs); @@ -208,15 +225,18 @@ Field3D operator-(const Field3D &lhs, const Field3D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(rhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ operator to update Field3D by subtraction with Field3D -Field3D &Field3D::operator-=(const Field3D &rhs) { +Field3D& Field3D::operator-=(const Field3D& rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { + ASSERT1(this->getCoordinateSystem() == rhs.getCoordinateSystem()) #if CHECK > 0 if (this->getLocation() != rhs.getLocation()) { @@ -243,7 +263,7 @@ Field3D &Field3D::operator-=(const Field3D &rhs) { } // Provide the C++ wrapper for multiplication of Field3D and Field2D -Field3D operator*(const Field3D &lhs, const Field2D &rhs) { +Field3D operator*(const Field3D& lhs, const Field2D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator*(Field3D, Field2D): fields at different " @@ -252,7 +272,7 @@ Field3D operator*(const Field3D &lhs, const Field2D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); @@ -270,16 +290,17 @@ Field3D operator*(const Field3D &lhs, const Field2D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(lhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ operator to update Field3D by multiplication with Field2D -Field3D &Field3D::operator*=(const Field2D &rhs) { +Field3D& Field3D::operator*=(const Field2D& rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { - #if CHECK > 0 if (this->getLocation() != rhs.getLocation()) { throw BoutException("Error in Field3D::operator*=(Field2D): fields at different " @@ -310,7 +331,7 @@ Field3D &Field3D::operator*=(const Field2D &rhs) { } // Provide the C++ wrapper for division of Field3D and Field2D -Field3D operator/(const Field3D &lhs, const Field2D &rhs) { +Field3D operator/(const Field3D& lhs, const Field2D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator/(Field3D, Field2D): fields at different " @@ -319,7 +340,7 @@ Field3D operator/(const Field3D &lhs, const Field2D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); @@ -338,16 +359,17 @@ Field3D operator/(const Field3D &lhs, const Field2D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(lhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ operator to update Field3D by division with Field2D -Field3D &Field3D::operator/=(const Field2D &rhs) { +Field3D& Field3D::operator/=(const Field2D& rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { - #if CHECK > 0 if (this->getLocation() != rhs.getLocation()) { throw BoutException("Error in Field3D::operator/=(Field2D): fields at different " @@ -379,7 +401,7 @@ Field3D &Field3D::operator/=(const Field2D &rhs) { } // Provide the C++ wrapper for addition of Field3D and Field2D -Field3D operator+(const Field3D &lhs, const Field2D &rhs) { +Field3D operator+(const Field3D& lhs, const Field2D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator+(Field3D, Field2D): fields at different " @@ -388,7 +410,7 @@ Field3D operator+(const Field3D &lhs, const Field2D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); @@ -406,16 +428,17 @@ Field3D operator+(const Field3D &lhs, const Field2D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(lhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ operator to update Field3D by addition with Field2D -Field3D &Field3D::operator+=(const Field2D &rhs) { +Field3D& Field3D::operator+=(const Field2D& rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { - #if CHECK > 0 if (this->getLocation() != rhs.getLocation()) { throw BoutException("Error in Field3D::operator+=(Field2D): fields at different " @@ -446,7 +469,7 @@ Field3D &Field3D::operator+=(const Field2D &rhs) { } // Provide the C++ wrapper for subtraction of Field3D and Field2D -Field3D operator-(const Field3D &lhs, const Field2D &rhs) { +Field3D operator-(const Field3D& lhs, const Field2D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator-(Field3D, Field2D): fields at different " @@ -455,7 +478,7 @@ Field3D operator-(const Field3D &lhs, const Field2D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); @@ -473,16 +496,17 @@ Field3D operator-(const Field3D &lhs, const Field2D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(lhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ operator to update Field3D by subtraction with Field2D -Field3D &Field3D::operator-=(const Field2D &rhs) { +Field3D& Field3D::operator-=(const Field2D& rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { - #if CHECK > 0 if (this->getLocation() != rhs.getLocation()) { throw BoutException("Error in Field3D::operator-=(Field2D): fields at different " @@ -513,9 +537,9 @@ Field3D &Field3D::operator-=(const Field2D &rhs) { } // Provide the C++ wrapper for multiplication of Field3D and BoutReal -Field3D operator*(const Field3D &lhs, const BoutReal rhs) { +Field3D operator*(const Field3D& lhs, const BoutReal rhs) { - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); Field3D result(localmesh); result.allocate(); @@ -526,12 +550,14 @@ Field3D operator*(const Field3D &lhs, const BoutReal rhs) { result.setLocation(lhs.getLocation()); + result.setCoordinateSystem(lhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ operator to update Field3D by multiplication with BoutReal -Field3D &Field3D::operator*=(const BoutReal rhs) { +Field3D& Field3D::operator*=(const BoutReal rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { @@ -550,9 +576,9 @@ Field3D &Field3D::operator*=(const BoutReal rhs) { } // Provide the C++ wrapper for division of Field3D and BoutReal -Field3D operator/(const Field3D &lhs, const BoutReal rhs) { +Field3D operator/(const Field3D& lhs, const BoutReal rhs) { - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); Field3D result(localmesh); result.allocate(); @@ -563,12 +589,14 @@ Field3D operator/(const Field3D &lhs, const BoutReal rhs) { result.setLocation(lhs.getLocation()); + result.setCoordinateSystem(lhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ operator to update Field3D by division with BoutReal -Field3D &Field3D::operator/=(const BoutReal rhs) { +Field3D& Field3D::operator/=(const BoutReal rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { @@ -587,9 +615,9 @@ Field3D &Field3D::operator/=(const BoutReal rhs) { } // Provide the C++ wrapper for addition of Field3D and BoutReal -Field3D operator+(const Field3D &lhs, const BoutReal rhs) { +Field3D operator+(const Field3D& lhs, const BoutReal rhs) { - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); Field3D result(localmesh); result.allocate(); @@ -600,12 +628,14 @@ Field3D operator+(const Field3D &lhs, const BoutReal rhs) { result.setLocation(lhs.getLocation()); + result.setCoordinateSystem(lhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ operator to update Field3D by addition with BoutReal -Field3D &Field3D::operator+=(const BoutReal rhs) { +Field3D& Field3D::operator+=(const BoutReal rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { @@ -624,9 +654,9 @@ Field3D &Field3D::operator+=(const BoutReal rhs) { } // Provide the C++ wrapper for subtraction of Field3D and BoutReal -Field3D operator-(const Field3D &lhs, const BoutReal rhs) { +Field3D operator-(const Field3D& lhs, const BoutReal rhs) { - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); Field3D result(localmesh); result.allocate(); @@ -637,12 +667,14 @@ Field3D operator-(const Field3D &lhs, const BoutReal rhs) { result.setLocation(lhs.getLocation()); + result.setCoordinateSystem(lhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ operator to update Field3D by subtraction with BoutReal -Field3D &Field3D::operator-=(const BoutReal rhs) { +Field3D& Field3D::operator-=(const BoutReal rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { @@ -661,7 +693,7 @@ Field3D &Field3D::operator-=(const BoutReal rhs) { } // Provide the C++ wrapper for multiplication of Field2D and Field3D -Field3D operator*(const Field2D &lhs, const Field3D &rhs) { +Field3D operator*(const Field2D& lhs, const Field3D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator*(Field2D, Field3D): fields at different " @@ -670,7 +702,7 @@ Field3D operator*(const Field2D &lhs, const Field3D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); @@ -688,12 +720,14 @@ Field3D operator*(const Field2D &lhs, const Field3D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(rhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ wrapper for division of Field2D and Field3D -Field3D operator/(const Field2D &lhs, const Field3D &rhs) { +Field3D operator/(const Field2D& lhs, const Field3D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator/(Field2D, Field3D): fields at different " @@ -702,7 +736,7 @@ Field3D operator/(const Field2D &lhs, const Field3D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); @@ -720,12 +754,14 @@ Field3D operator/(const Field2D &lhs, const Field3D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(rhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ wrapper for addition of Field2D and Field3D -Field3D operator+(const Field2D &lhs, const Field3D &rhs) { +Field3D operator+(const Field2D& lhs, const Field3D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator+(Field2D, Field3D): fields at different " @@ -734,7 +770,7 @@ Field3D operator+(const Field2D &lhs, const Field3D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); @@ -752,12 +788,14 @@ Field3D operator+(const Field2D &lhs, const Field3D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(rhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ wrapper for subtraction of Field2D and Field3D -Field3D operator-(const Field2D &lhs, const Field3D &rhs) { +Field3D operator-(const Field2D& lhs, const Field3D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator-(Field2D, Field3D): fields at different " @@ -766,7 +804,7 @@ Field3D operator-(const Field2D &lhs, const Field3D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); @@ -784,12 +822,14 @@ Field3D operator-(const Field2D &lhs, const Field3D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(rhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ wrapper for multiplication of Field2D and Field2D -Field2D operator*(const Field2D &lhs, const Field2D &rhs) { +Field2D operator*(const Field2D& lhs, const Field2D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator*(Field2D, Field2D): fields at different " @@ -798,7 +838,7 @@ Field2D operator*(const Field2D &lhs, const Field2D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); @@ -818,11 +858,10 @@ Field2D operator*(const Field2D &lhs, const Field2D &rhs) { } // Provide the C++ operator to update Field2D by multiplication with Field2D -Field2D &Field2D::operator*=(const Field2D &rhs) { +Field2D& Field2D::operator*=(const Field2D& rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { - #if CHECK > 0 if (this->getLocation() != rhs.getLocation()) { throw BoutException("Error in Field2D::operator*=(Field2D): fields at different " @@ -848,7 +887,7 @@ Field2D &Field2D::operator*=(const Field2D &rhs) { } // Provide the C++ wrapper for division of Field2D and Field2D -Field2D operator/(const Field2D &lhs, const Field2D &rhs) { +Field2D operator/(const Field2D& lhs, const Field2D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator/(Field2D, Field2D): fields at different " @@ -857,7 +896,7 @@ Field2D operator/(const Field2D &lhs, const Field2D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); @@ -877,11 +916,10 @@ Field2D operator/(const Field2D &lhs, const Field2D &rhs) { } // Provide the C++ operator to update Field2D by division with Field2D -Field2D &Field2D::operator/=(const Field2D &rhs) { +Field2D& Field2D::operator/=(const Field2D& rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { - #if CHECK > 0 if (this->getLocation() != rhs.getLocation()) { throw BoutException("Error in Field2D::operator/=(Field2D): fields at different " @@ -907,7 +945,7 @@ Field2D &Field2D::operator/=(const Field2D &rhs) { } // Provide the C++ wrapper for addition of Field2D and Field2D -Field2D operator+(const Field2D &lhs, const Field2D &rhs) { +Field2D operator+(const Field2D& lhs, const Field2D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator+(Field2D, Field2D): fields at different " @@ -916,7 +954,7 @@ Field2D operator+(const Field2D &lhs, const Field2D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); @@ -936,11 +974,10 @@ Field2D operator+(const Field2D &lhs, const Field2D &rhs) { } // Provide the C++ operator to update Field2D by addition with Field2D -Field2D &Field2D::operator+=(const Field2D &rhs) { +Field2D& Field2D::operator+=(const Field2D& rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { - #if CHECK > 0 if (this->getLocation() != rhs.getLocation()) { throw BoutException("Error in Field2D::operator+=(Field2D): fields at different " @@ -966,7 +1003,7 @@ Field2D &Field2D::operator+=(const Field2D &rhs) { } // Provide the C++ wrapper for subtraction of Field2D and Field2D -Field2D operator-(const Field2D &lhs, const Field2D &rhs) { +Field2D operator-(const Field2D& lhs, const Field2D& rhs) { #if CHECK > 0 if (lhs.getLocation() != rhs.getLocation()) { throw BoutException("Error in operator-(Field2D, Field2D): fields at different " @@ -975,7 +1012,7 @@ Field2D operator-(const Field2D &lhs, const Field2D &rhs) { } #endif - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); ASSERT1(localmesh == rhs.getMesh()); @@ -995,11 +1032,10 @@ Field2D operator-(const Field2D &lhs, const Field2D &rhs) { } // Provide the C++ operator to update Field2D by subtraction with Field2D -Field2D &Field2D::operator-=(const Field2D &rhs) { +Field2D& Field2D::operator-=(const Field2D& rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { - #if CHECK > 0 if (this->getLocation() != rhs.getLocation()) { throw BoutException("Error in Field2D::operator-=(Field2D): fields at different " @@ -1025,9 +1061,9 @@ Field2D &Field2D::operator-=(const Field2D &rhs) { } // Provide the C++ wrapper for multiplication of Field2D and BoutReal -Field2D operator*(const Field2D &lhs, const BoutReal rhs) { +Field2D operator*(const Field2D& lhs, const BoutReal rhs) { - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); Field2D result(localmesh); result.allocate(); @@ -1043,7 +1079,7 @@ Field2D operator*(const Field2D &lhs, const BoutReal rhs) { } // Provide the C++ operator to update Field2D by multiplication with BoutReal -Field2D &Field2D::operator*=(const BoutReal rhs) { +Field2D& Field2D::operator*=(const BoutReal rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { @@ -1062,9 +1098,9 @@ Field2D &Field2D::operator*=(const BoutReal rhs) { } // Provide the C++ wrapper for division of Field2D and BoutReal -Field2D operator/(const Field2D &lhs, const BoutReal rhs) { +Field2D operator/(const Field2D& lhs, const BoutReal rhs) { - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); Field2D result(localmesh); result.allocate(); @@ -1080,7 +1116,7 @@ Field2D operator/(const Field2D &lhs, const BoutReal rhs) { } // Provide the C++ operator to update Field2D by division with BoutReal -Field2D &Field2D::operator/=(const BoutReal rhs) { +Field2D& Field2D::operator/=(const BoutReal rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { @@ -1099,9 +1135,9 @@ Field2D &Field2D::operator/=(const BoutReal rhs) { } // Provide the C++ wrapper for addition of Field2D and BoutReal -Field2D operator+(const Field2D &lhs, const BoutReal rhs) { +Field2D operator+(const Field2D& lhs, const BoutReal rhs) { - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); Field2D result(localmesh); result.allocate(); @@ -1117,7 +1153,7 @@ Field2D operator+(const Field2D &lhs, const BoutReal rhs) { } // Provide the C++ operator to update Field2D by addition with BoutReal -Field2D &Field2D::operator+=(const BoutReal rhs) { +Field2D& Field2D::operator+=(const BoutReal rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { @@ -1136,9 +1172,9 @@ Field2D &Field2D::operator+=(const BoutReal rhs) { } // Provide the C++ wrapper for subtraction of Field2D and BoutReal -Field2D operator-(const Field2D &lhs, const BoutReal rhs) { +Field2D operator-(const Field2D& lhs, const BoutReal rhs) { - Mesh *localmesh = lhs.getMesh(); + Mesh* localmesh = lhs.getMesh(); Field2D result(localmesh); result.allocate(); @@ -1154,7 +1190,7 @@ Field2D operator-(const Field2D &lhs, const BoutReal rhs) { } // Provide the C++ operator to update Field2D by subtraction with BoutReal -Field2D &Field2D::operator-=(const BoutReal rhs) { +Field2D& Field2D::operator-=(const BoutReal rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { @@ -1173,9 +1209,9 @@ Field2D &Field2D::operator-=(const BoutReal rhs) { } // Provide the C++ wrapper for multiplication of BoutReal and Field3D -Field3D operator*(const BoutReal lhs, const Field3D &rhs) { +Field3D operator*(const BoutReal lhs, const Field3D& rhs) { - Mesh *localmesh = rhs.getMesh(); + Mesh* localmesh = rhs.getMesh(); Field3D result(localmesh); result.allocate(); @@ -1186,14 +1222,16 @@ Field3D operator*(const BoutReal lhs, const Field3D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(rhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ wrapper for division of BoutReal and Field3D -Field3D operator/(const BoutReal lhs, const Field3D &rhs) { +Field3D operator/(const BoutReal lhs, const Field3D& rhs) { - Mesh *localmesh = rhs.getMesh(); + Mesh* localmesh = rhs.getMesh(); Field3D result(localmesh); result.allocate(); @@ -1204,14 +1242,16 @@ Field3D operator/(const BoutReal lhs, const Field3D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(rhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ wrapper for addition of BoutReal and Field3D -Field3D operator+(const BoutReal lhs, const Field3D &rhs) { +Field3D operator+(const BoutReal lhs, const Field3D& rhs) { - Mesh *localmesh = rhs.getMesh(); + Mesh* localmesh = rhs.getMesh(); Field3D result(localmesh); result.allocate(); @@ -1222,14 +1262,16 @@ Field3D operator+(const BoutReal lhs, const Field3D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(rhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ wrapper for subtraction of BoutReal and Field3D -Field3D operator-(const BoutReal lhs, const Field3D &rhs) { +Field3D operator-(const BoutReal lhs, const Field3D& rhs) { - Mesh *localmesh = rhs.getMesh(); + Mesh* localmesh = rhs.getMesh(); Field3D result(localmesh); result.allocate(); @@ -1240,14 +1282,16 @@ Field3D operator-(const BoutReal lhs, const Field3D &rhs) { result.setLocation(rhs.getLocation()); + result.setCoordinateSystem(rhs.getCoordinateSystem()); + checkData(result); return result; } // Provide the C++ wrapper for multiplication of BoutReal and Field2D -Field2D operator*(const BoutReal lhs, const Field2D &rhs) { +Field2D operator*(const BoutReal lhs, const Field2D& rhs) { - Mesh *localmesh = rhs.getMesh(); + Mesh* localmesh = rhs.getMesh(); Field2D result(localmesh); result.allocate(); @@ -1263,9 +1307,9 @@ Field2D operator*(const BoutReal lhs, const Field2D &rhs) { } // Provide the C++ wrapper for division of BoutReal and Field2D -Field2D operator/(const BoutReal lhs, const Field2D &rhs) { +Field2D operator/(const BoutReal lhs, const Field2D& rhs) { - Mesh *localmesh = rhs.getMesh(); + Mesh* localmesh = rhs.getMesh(); Field2D result(localmesh); result.allocate(); @@ -1281,9 +1325,9 @@ Field2D operator/(const BoutReal lhs, const Field2D &rhs) { } // Provide the C++ wrapper for addition of BoutReal and Field2D -Field2D operator+(const BoutReal lhs, const Field2D &rhs) { +Field2D operator+(const BoutReal lhs, const Field2D& rhs) { - Mesh *localmesh = rhs.getMesh(); + Mesh* localmesh = rhs.getMesh(); Field2D result(localmesh); result.allocate(); @@ -1299,9 +1343,9 @@ Field2D operator+(const BoutReal lhs, const Field2D &rhs) { } // Provide the C++ wrapper for subtraction of BoutReal and Field2D -Field2D operator-(const BoutReal lhs, const Field2D &rhs) { +Field2D operator-(const BoutReal lhs, const Field2D& rhs) { - Mesh *localmesh = rhs.getMesh(); + Mesh* localmesh = rhs.getMesh(); Field2D result(localmesh); result.allocate(); diff --git a/src/fileio/datafile.cxx b/src/fileio/datafile.cxx index b6ef8a1298..41dcb1332e 100644 --- a/src/fileio/datafile.cxx +++ b/src/fileio/datafile.cxx @@ -930,6 +930,12 @@ bool Datafile::write() { // Write 3D fields for (const auto& var : f3d_arr) { write_f3d(var.name, var.ptr, var.save_repeat); + // Add coordinate system + if (shiftOutput) { + file->setAttribute(var.name, "coordinate_system", COORDINATE_SYSTEMtoString.at(COORDINATE_SYSTEM::FieldAligned)); + } else { + file->setAttribute(var.name, "coordinate_system", COORDINATE_SYSTEMtoString.at(var.ptr->getCoordinateSystem())); + } } // 2D vectors diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx index ea06fa66e7..71090eca61 100644 --- a/src/mesh/index_derivs.cxx +++ b/src/mesh/index_derivs.cxx @@ -45,7 +45,6 @@ STAGGER Mesh::getStagger(const CELL_LOC inloc, const CELL_LOC outloc, TRACE("Mesh::getStagger -- three arguments"); ASSERT1(outloc == inloc || (outloc == CELL_CENTRE && inloc == allowedStaggerLoc) || (outloc == allowedStaggerLoc && inloc == CELL_CENTRE)); - if ((!StaggerGrids) || outloc == inloc) return STAGGER::None; if (outloc == allowedStaggerLoc) { diff --git a/src/mesh/interpolation.cxx b/src/mesh/interpolation.cxx index 45b7d47d9b..b687030617 100644 --- a/src/mesh/interpolation.cxx +++ b/src/mesh/interpolation.cxx @@ -161,6 +161,7 @@ const Field3D interp_to(const Field3D &var, CELL_LOC loc, REGION region) { result = var_fa; // NOTE: This is just for boundaries. FIX! result.allocate(); } + result.setCoordinateSystem(COORDINATE_SYSTEM::FieldAligned); if (fieldmesh->ystart > 1) { // More than one guard cell, so set pp and mm values diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx index 747f311c03..6e3e7c90a3 100644 --- a/src/mesh/parallel/fci.hxx +++ b/src/mesh/parallel/fci.hxx @@ -81,6 +81,10 @@ public: throw BoutException("FCI method cannot transform into field aligned grid"); } + COORDINATE_SYSTEM getCoordinateSystem() const override { + return COORDINATE_SYSTEM::FCI; + } + bool canToFromFieldAligned() override{ return false; } diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx index 7a7c439c14..b3ddcb8676 100644 --- a/src/mesh/parallel/shiftedmetric.cxx +++ b/src/mesh/parallel/shiftedmetric.cxx @@ -106,7 +106,10 @@ void ShiftedMetric::calcYUpDown(Field3D &f) { * and Y is then field aligned. */ const Field3D ShiftedMetric::toFieldAligned(const Field3D &f, const REGION region) { - return shiftZ(f, toAlignedPhs, region); + ASSERT2(f.getCoordinateSystem() == COORDINATE_SYSTEM::Orthogonal); + Field3D result = shiftZ(f, toAlignedPhs, region); + result.setCoordinateSystem(COORDINATE_SYSTEM::FieldAligned); + return result; } /*! @@ -114,7 +117,10 @@ const Field3D ShiftedMetric::toFieldAligned(const Field3D &f, const REGION regio * but Y is not field aligned. */ const Field3D ShiftedMetric::fromFieldAligned(const Field3D &f, const REGION region) { - return shiftZ(f, fromAlignedPhs, region); + ASSERT2(f.getCoordinateSystem() == COORDINATE_SYSTEM::FieldAligned); + Field3D result = shiftZ(f, fromAlignedPhs, region); + result.setCoordinateSystem(f.getMesh()->getCoordinateSystem()); + return result; } const Field3D ShiftedMetric::shiftZ(const Field3D& f, const Tensor& phs, From 5fb0d7da2707974f7b6fe7cdf3c52edb7f37cd46 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Fri, 14 Dec 2018 09:10:52 +0000 Subject: [PATCH 2/7] Add coordinate_system to swap(Field3D,Field3D) --- include/field3d.hxx | 1 + 1 file changed, 1 insertion(+) diff --git a/include/field3d.hxx b/include/field3d.hxx index 3b70a4e95e..bd5d8eee53 100644 --- a/include/field3d.hxx +++ b/include/field3d.hxx @@ -492,6 +492,7 @@ class Field3D : public Field, public FieldData { swap(first.nx, second.nx); swap(first.ny, second.ny); swap(first.nz, second.nz); + swap(first.coordinate_system, second.coordinate_system); swap(first.location, second.location); swap(first.deriv, second.deriv); swap(first.yup_field, second.yup_field); From a1699f6d75b4b696cb94551b27d92b6975ffb43f Mon Sep 17 00:00:00 2001 From: John Omotani Date: Sun, 16 Dec 2018 14:05:58 +0000 Subject: [PATCH 3/7] Set CoordinateSystem of result in derivatives Derivative functions are templates, so Field2D needs get/setCoordinates methods (even if these are never called). --- include/bout/index_derivs_interface.hxx | 6 ++++++ include/field.hxx | 22 ++++++++++++++++++++++ include/field2d.hxx | 18 ++++++++++++++++++ include/field3d.hxx | 22 ---------------------- 4 files changed, 46 insertions(+), 22 deletions(-) diff --git a/include/bout/index_derivs_interface.hxx b/include/bout/index_derivs_interface.hxx index 0bbfad8cfc..971d866c16 100644 --- a/include/bout/index_derivs_interface.hxx +++ b/include/bout/index_derivs_interface.hxx @@ -96,6 +96,9 @@ T flowDerivative(const T& vel, const T& f, CELL_LOC outloc, const std::string& m T result(localmesh); result.allocate(); // Make sure data allocated result.setLocation(outloc); + if (std::is_base_of::value) { + result.setCoordinateSystem(f.getCoordinateSystem()); + } // Apply method derivativeMethod(vel, f, result, region); @@ -162,6 +165,9 @@ T standardDerivative(const T& f, CELL_LOC outloc, const std::string& method, T result(localmesh); result.allocate(); // Make sure data allocated result.setLocation(outloc); + if (std::is_base_of::value) { + result.setCoordinateSystem(f.getCoordinateSystem()); + } // Apply method derivativeMethod(f, result, region); diff --git a/include/field.hxx b/include/field.hxx index 13c288f635..9c17b055ac 100644 --- a/include/field.hxx +++ b/include/field.hxx @@ -46,6 +46,28 @@ extern Mesh * mesh; ///< Global mesh #include #endif +/*! + * Coordinate systems that Field3Ds can be stored in + * + * None - special value to be used for Field3Ds initialized before the Mesh and + * ParallelTransform have finished being initialized, because then we cannot + * get the COORDINATE_SYSTEM to use from the ParallelTransform, but Field3Ds + * created then must not need to know their COORDINATE_SYSTEM, or must set it + * explicitly later. + */ +enum class COORDINATE_SYSTEM { None, FieldAligned, Orthogonal, FCI }; + +#define COORDENUMSTR(val) {COORDINATE_SYSTEM::val, #val} + +const std::map COORDINATE_SYSTEMtoString = { + COORDENUMSTR(None), + COORDENUMSTR(FieldAligned), + COORDENUMSTR(Orthogonal), + COORDENUMSTR(FCI) +}; + +#undef COORDENUMSTR + /*! * \brief Base class for fields * diff --git a/include/field2d.hxx b/include/field2d.hxx index adbf8b0800..d5174750c4 100644 --- a/include/field2d.hxx +++ b/include/field2d.hxx @@ -111,6 +111,24 @@ class Field2D : public Field, public FieldData { */ int getNz() const override {return 1;}; + /*! + * Return the coordinate system of this field + */ + COORDINATE_SYSTEM getCoordinateSystem() const { +#if CHECK > 0 + throw BoutException("getCoordinateSystem() should not be called for Field2D"); +#endif + } + + /*! + * Reset the coordinate system of this field + */ + void setCoordinateSystem(COORDINATE_SYSTEM UNUSED(new_coords)) { +#if CHECK > 0 + throw BoutException("setCoordinateSystem() should not be called for Field2D"); +#endif + } + /// Check if this field has yup and ydown fields bool hasYupYdown() const { return true; diff --git a/include/field3d.hxx b/include/field3d.hxx index bd5d8eee53..56269fece9 100644 --- a/include/field3d.hxx +++ b/include/field3d.hxx @@ -40,28 +40,6 @@ class Mesh; // #include "bout/mesh.hxx" #include "bout/field_visitor.hxx" -/*! - * Coordinate systems that Field3Ds can be stored in - * - * None - special value to be used for Field3Ds initialized before the Mesh and - * ParallelTransform have finished being initialized, because then we cannot - * get the COORDINATE_SYSTEM to use from the ParallelTransform, but Field3Ds - * created then must not need to know their COORDINATE_SYSTEM, or must set it - * explicitly later. - */ -enum class COORDINATE_SYSTEM { None, FieldAligned, Orthogonal, FCI }; - -#define COORDENUMSTR(val) {COORDINATE_SYSTEM::val, #val} - -const std::map COORDINATE_SYSTEMtoString = { - COORDENUMSTR(None), - COORDENUMSTR(FieldAligned), - COORDENUMSTR(Orthogonal), - COORDENUMSTR(FCI) -}; - -#undef COORDENUMSTR - /// Class for 3D X-Y-Z scalar fields /*! This class represents a scalar field defined over the mesh. From a8e1e99e837e40a0bc0575393233d449d6776801 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 24 Jul 2018 19:27:01 +0100 Subject: [PATCH 4/7] Avoid circular dependencies Initialize ParallelTransform in BoutMesh::load(), so it is almost always available. Creating the FCITransform requires nx/ny/nz to be set in the mesh, so cannot call setParallelTransform() in the Mesh constructor, since it is called before these are initialized. --- include/bout/mesh.hxx | 14 ++++++++++---- src/field/field3d.cxx | 7 +++---- src/mesh/impls/bout/boutmesh.cxx | 5 ++++- src/mesh/mesh.cxx | 5 +---- 4 files changed, 18 insertions(+), 13 deletions(-) diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index cec5a6eacb..61bb21e1ea 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -681,7 +681,7 @@ class Mesh { /// Transform a field into field-aligned coordinates const Field3D toFieldAligned(const Field3D &f, const REGION region = RGN_NOX) { - return getParallelTransform().toFieldAligned(f, region); + return transform->toFieldAligned(f, region); } const Field2D toFieldAligned(const Field2D &f, const REGION UNUSED(region) = RGN_NOX) { return f; @@ -689,18 +689,24 @@ class Mesh { /// Convert back into standard form const Field3D fromFieldAligned(const Field3D &f, const REGION region = RGN_NOX) { - return getParallelTransform().fromFieldAligned(f, region); + return transform->fromFieldAligned(f, region); } const Field2D fromFieldAligned(const Field2D &f, const REGION UNUSED(region) = RGN_NOX) { return f; } COORDINATE_SYSTEM getCoordinateSystem() const { - return getParallelTransform().getCoordinateSystem(); + if (transform) { + return transform->getCoordinateSystem(); + } else { + // ParallelTransform not initialized yet. Any Field3D created better not + // need to know its coordinate system. + return COORDINATE_SYSTEM::None; + } } bool canToFromFieldAligned() { - return getParallelTransform().canToFromFieldAligned(); + return transform->canToFromFieldAligned(); } /*! diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 6f0650dac7..428db8e274 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -87,14 +87,13 @@ Field3D::Field3D(const Field2D& f) : Field(f.getMesh()) { TRACE("Field3D: Copy constructor from Field2D"); - coordinate_system = fieldmesh->getCoordinateSystem(); - nx = fieldmesh->LocalNx; ny = fieldmesh->LocalNy; nz = fieldmesh->LocalNz; location = f.getLocation(); fieldCoordinates = nullptr; + coordinate_system = fieldmesh->getCoordinateSystem(); *this = f; } @@ -103,12 +102,12 @@ Field3D::Field3D(const BoutReal val, Mesh* localmesh) : Field(localmesh) { TRACE("Field3D: Copy constructor from value"); - coordinate_system = fieldmesh->getCoordinateSystem(); - nx = fieldmesh->LocalNx; ny = fieldmesh->LocalNy; nz = fieldmesh->LocalNz; + coordinate_system = fieldmesh->getCoordinateSystem(); + *this = val; } diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index 0bf59de09a..4548a9368f 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -839,7 +839,10 @@ int BoutMesh::load() { // Add boundary regions addBoundaryRegions(); - output_info.write(_("\tdone\n")); + // Set the ParallelTransform + setParallelTransform(); + + output_info.write("\tdone\n"); return 0; } diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx index f645b89f0c..2152052d30 100644 --- a/src/mesh/mesh.cxx +++ b/src/mesh/mesh.cxx @@ -316,10 +316,7 @@ void Mesh::setParallelTransform() { } ParallelTransform& Mesh::getParallelTransform() { - if(!transform) { - // No ParallelTransform object yet. Set from options - setParallelTransform(); - } + ASSERT1(transform); // Return a reference to the ParallelTransform object return *transform; From 15e40bba648620738ed627818e0915c60be0b4d6 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 25 Jul 2018 17:44:05 +0100 Subject: [PATCH 5/7] Fix coordinate system tracking in test-yupdown test-yupdown has a custom field-aligned derivative that did not set the coordinate system, causing an exception when checking was turned on. --- tests/integrated/test-yupdown/test_yupdown.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/integrated/test-yupdown/test_yupdown.cxx b/tests/integrated/test-yupdown/test_yupdown.cxx index dbadc6372f..cd5ea6ed3b 100644 --- a/tests/integrated/test-yupdown/test_yupdown.cxx +++ b/tests/integrated/test-yupdown/test_yupdown.cxx @@ -24,6 +24,7 @@ const Field3D DDY_aligned(const Field3D &f) { Field3D result; result.allocate(); result = 0.0; + result.setCoordinateSystem(f.getCoordinateSystem()); for(int i=0;iLocalNx;i++) for(int j=mesh->ystart;j<=mesh->yend;j++) From bfd1f7ab17275c386c36c17a76658461e644bfd5 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 17 Dec 2018 22:59:35 +0000 Subject: [PATCH 6/7] Add COORDINATE_SYSTEM setting and checks to ShiftToFieldAligned Also reduce CHECK level of COORDINATE_SYSTEM asserts in ShiftedMetric from 2 to 1. --- src/mesh/parallel/shiftedmetric.cxx | 4 ++-- src/mesh/parallel/shifttofieldaligned.cxx | 14 ++++++++++++-- src/mesh/parallel/shifttofieldaligned.hxx | 2 ++ 3 files changed, 16 insertions(+), 4 deletions(-) diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx index b3ddcb8676..231d63188e 100644 --- a/src/mesh/parallel/shiftedmetric.cxx +++ b/src/mesh/parallel/shiftedmetric.cxx @@ -106,7 +106,7 @@ void ShiftedMetric::calcYUpDown(Field3D &f) { * and Y is then field aligned. */ const Field3D ShiftedMetric::toFieldAligned(const Field3D &f, const REGION region) { - ASSERT2(f.getCoordinateSystem() == COORDINATE_SYSTEM::Orthogonal); + ASSERT1(f.getCoordinateSystem() == COORDINATE_SYSTEM::Orthogonal); Field3D result = shiftZ(f, toAlignedPhs, region); result.setCoordinateSystem(COORDINATE_SYSTEM::FieldAligned); return result; @@ -117,7 +117,7 @@ const Field3D ShiftedMetric::toFieldAligned(const Field3D &f, const REGION regio * but Y is not field aligned. */ const Field3D ShiftedMetric::fromFieldAligned(const Field3D &f, const REGION region) { - ASSERT2(f.getCoordinateSystem() == COORDINATE_SYSTEM::FieldAligned); + ASSERT1(f.getCoordinateSystem() == COORDINATE_SYSTEM::FieldAligned); Field3D result = shiftZ(f, fromAlignedPhs, region); result.setCoordinateSystem(f.getMesh()->getCoordinateSystem()); return result; diff --git a/src/mesh/parallel/shifttofieldaligned.cxx b/src/mesh/parallel/shifttofieldaligned.cxx index 9b760fd929..8bb06f02c8 100644 --- a/src/mesh/parallel/shifttofieldaligned.cxx +++ b/src/mesh/parallel/shifttofieldaligned.cxx @@ -120,7 +120,14 @@ void ShiftToFieldAligned::Implementation::cachePhases() { */ const Field3D ShiftToFieldAligned::Implementation::toFieldAligned(const Field3D& f, const REGION region) { - return shiftZ(f, toAlignedPhs, region); + if (f.getCoordinateSystem() == COORDINATE_SYSTEM::FieldAligned) { + return f; + } else { + ASSERT1(f.getCoordinateSystem() == COORDINATE_SYSTEM::Orthogonal); + Field3D result = shiftZ(f, toAlignedPhs, region); + result.setCoordinateSystem(COORDINATE_SYSTEM::FieldAligned); + return result; + } } /*! @@ -129,7 +136,10 @@ const Field3D ShiftToFieldAligned::Implementation::toFieldAligned(const Field3D& */ const Field3D ShiftToFieldAligned::Implementation::fromFieldAligned(const Field3D& f, const REGION region) { - return shiftZ(f, fromAlignedPhs, region); + ASSERT1(f.getCoordinateSystem() == COORDINATE_SYSTEM::FieldAligned); + Field3D result = shiftZ(f, fromAlignedPhs, region); + result.setCoordinateSystem(COORDINATE_SYSTEM::Orthogonal); + return result; } const Field3D ShiftToFieldAligned::Implementation::shiftZ( diff --git a/src/mesh/parallel/shifttofieldaligned.hxx b/src/mesh/parallel/shifttofieldaligned.hxx index 57b8fcfc6c..31b404819e 100644 --- a/src/mesh/parallel/shifttofieldaligned.hxx +++ b/src/mesh/parallel/shifttofieldaligned.hxx @@ -61,6 +61,8 @@ public: return implementations.at(f.getLocation()).fromFieldAligned(f, region); } + COORDINATE_SYSTEM getCoordinateSystem() const { return COORDINATE_SYSTEM::Orthogonal; } + bool canToFromFieldAligned() override { return true; } private: From 7eb5c0f8c8257b9a0fe8cf5a6d4c7dcb4c4b6c7a Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 17 Dec 2018 23:02:03 +0000 Subject: [PATCH 7/7] Fix test-yupdown-shifttofieldaligned by setting coordinate system --- .../integrated/test-yupdown-shifttofieldaligned/test_yupdown.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/integrated/test-yupdown-shifttofieldaligned/test_yupdown.cxx b/tests/integrated/test-yupdown-shifttofieldaligned/test_yupdown.cxx index f7cef52403..9fa03bb49a 100644 --- a/tests/integrated/test-yupdown-shifttofieldaligned/test_yupdown.cxx +++ b/tests/integrated/test-yupdown-shifttofieldaligned/test_yupdown.cxx @@ -9,6 +9,7 @@ const Field3D DDY_aligned(const Field3D &f) { Field3D result; result.allocate(); result = 0.0; + result.setCoordinateSystem(f.getCoordinateSystem()); for(int i=0;iLocalNx;i++) for(int j=mesh->ystart;j<=mesh->yend;j++)