From f667ea6106860b058f13bd534db0386eef70f120 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 25 Oct 2018 13:30:05 +0100 Subject: [PATCH 01/34] Add multiple yup/ydown fields --- include/field3d.hxx | 48 +++++++++------ src/field/field3d.cxx | 97 +++++++++++++------------------ tests/unit/field/test_field3d.cxx | 28 ++++++++- 3 files changed, 100 insertions(+), 73 deletions(-) diff --git a/include/field3d.hxx b/include/field3d.hxx index e9f090bea6..3a774e2d48 100644 --- a/include/field3d.hxx +++ b/include/field3d.hxx @@ -40,6 +40,8 @@ class Mesh; // #include "bout/mesh.hxx" #include "bout/field_visitor.hxx" +#include + /// Class for 3D X-Y-Z scalar fields /*! This class represents a scalar field defined over the mesh. @@ -226,30 +228,42 @@ class Field3D : public Field, public FieldData { /// Check if this field has yup and ydown fields bool hasYupYdown() const { - return (yup_field != nullptr) && (ydown_field != nullptr); + return !yup_fields.empty() and !ydown_fields.empty(); } /// Return reference to yup field - Field3D& yup() { - ASSERT2(yup_field != nullptr); // Check for communicate - return *yup_field; + Field3D &yup(std::vector::size_type index = 0) { + if (yup_fields.empty()) { + return *this; + } + ASSERT2(index < yup_fields.size()); // Check for communicate + return yup_fields[index]; } /// Return const reference to yup field - const Field3D& yup() const { - ASSERT2(yup_field != nullptr); - return *yup_field; + const Field3D &yup(std::vector::size_type index = 0) const { + if (yup_fields.empty()) { + return *this; + } + ASSERT2(index < yup_fields.size()); + return yup_fields[index]; } - + /// Return reference to ydown field - Field3D& ydown() { - ASSERT2(ydown_field != nullptr); - return *ydown_field; + Field3D &ydown(std::vector::size_type index = 0) { + if (ydown_fields.empty()) { + return *this; + } + ASSERT2(index < ydown_fields.size()); + return ydown_fields[index]; } - + /// Return const reference to ydown field - const Field3D& ydown() const { - ASSERT2(ydown_field != nullptr); - return *ydown_field; + const Field3D &ydown(std::vector::size_type index = 0) const { + if (ydown_fields.empty()) { + return *this; + } + ASSERT2(index < ydown_fields.size()); + return ydown_fields[index]; } /// Return yup if dir=+1, and ydown if dir=-1 @@ -457,8 +471,8 @@ private: Field3D *deriv; ///< Time derivative (may be NULL) - /// Pointers to fields containing values along Y - Field3D *yup_field, *ydown_field; + /// Fields containing values along Y + std::vector yup_fields{}, ydown_fields{}; }; // Non-member overloaded operators diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 7d9dda7cea..519db35b08 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -45,8 +45,7 @@ /// Constructor Field3D::Field3D(Mesh *localmesh) - : Field(localmesh), background(nullptr), deriv(nullptr), yup_field(nullptr), - ydown_field(nullptr) { + : Field(localmesh), background(nullptr), deriv(nullptr) { #ifdef TRACK name = ""; #endif @@ -55,6 +54,9 @@ Field3D::Field3D(Mesh *localmesh) nx = fieldmesh->LocalNx; ny = fieldmesh->LocalNy; nz = fieldmesh->LocalNz; + + yup_fields.reserve(fieldmesh->ystart); + ydown_fields.reserve(fieldmesh->ystart); } #if CHECK > 0 else { @@ -72,7 +74,7 @@ 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) { TRACE("Field3D(Field3D&)"); @@ -84,6 +86,9 @@ Field3D::Field3D(const Field3D &f) nx = fieldmesh->LocalNx; ny = fieldmesh->LocalNy; nz = fieldmesh->LocalNz; + + yup_fields.reserve(fieldmesh->ystart); + ydown_fields.reserve(fieldmesh->ystart); } #if CHECK > 0 else { @@ -95,13 +100,12 @@ Field3D::Field3D(const Field3D &f) location = f.location; fieldCoordinates = f.fieldCoordinates; - + boundaryIsSet = false; } 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) { TRACE("Field3D: Copy constructor from Field2D"); @@ -118,8 +122,7 @@ 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) { TRACE("Field3D: Copy constructor from value"); @@ -135,23 +138,8 @@ Field3D::Field3D(const BoutReal val, Mesh *localmesh) Field3D::~Field3D() { /// Delete the time derivative variable if allocated if (deriv != nullptr) { - // 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; - - // 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; } void Field3D::allocate() { @@ -181,54 +169,53 @@ Field3D* Field3D::timeDeriv() { void Field3D::splitYupYdown() { TRACE("Field3D::splitYupYdown"); - if((yup_field != this) && (yup_field != nullptr)) + if (!yup_fields.empty()) { return; + } - // yup_field and ydown_field null - yup_field = new Field3D(fieldmesh); - ydown_field = new Field3D(fieldmesh); + for (int i = 0; i < fieldmesh->ystart; ++i) { + yup_fields.emplace_back(fieldmesh); + ydown_fields.emplace_back(fieldmesh); + } } void Field3D::mergeYupYdown() { TRACE("Field3D::mergeYupYdown"); - - if(yup_field == this && ydown_field == this) - return; - if(yup_field != nullptr){ - delete yup_field; - } - - if(ydown_field != nullptr) { - delete ydown_field; + if (yup_fields.empty() && ydown_fields.empty()) { + return; } - yup_field = this; - ydown_field = this; + yup_fields.clear(); + ydown_fields.clear(); } -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); +const Field3D& Field3D::ynext(int dir) const { + // Asked for more than yguards + if (std::abs(dir) > fieldmesh->ystart) { + throw BoutException( + "Field3D: Call to ynext with %d which is more than number of yguards (%d)", dir, + fieldmesh->ystart); } -} -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); + // ynext uses 1-indexing, but yup wants 0-indexing + if (dir > 0) { + return yup(dir - 1); + } else if (dir < 0) { + return ydown(std::abs(dir) - 1); + } else { + return *this; } } +Field3D &Field3D::ynext(int dir) { + // Call the `const` version: need to add `const` to `this` to call + // it, then throw it away after. This is ok because `this` wasn't + // `const` to begin with. + // See Effective C++, Scott Meyers, p23, for a better explanation + return const_cast(static_cast(*this).ynext(dir)); +} + void Field3D::setLocation(CELL_LOC new_location) { if (getMesh()->StaggerGrids) { if (new_location == CELL_VSHIFT) { diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index 649a3abb1b..775a383621 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -243,7 +243,7 @@ TEST_F(Field3DTest, MergeYupYDown) { field.mergeYupYdown(); - EXPECT_TRUE(field.hasYupYdown()); + EXPECT_FALSE(field.hasYupYdown()); auto& yup = field.yup(); EXPECT_EQ(&field, &yup); @@ -278,6 +278,32 @@ TEST_F(Field3DTest, SplitThenMergeYupYDown) { EXPECT_EQ(&field, &ydown2); } +TEST_F(Field3DTest, MultipleYupYdown) { + FakeMesh newmesh{3, 5, 7}; + newmesh.ystart = 2; + + Field3D field{&newmesh}; + + field.splitYupYdown(); + + EXPECT_TRUE(field.hasYupYdown()); + + auto &yup = field.yup(); + EXPECT_NE(&field, &yup); + auto &ydown = field.ydown(); + EXPECT_NE(&field, &ydown); + auto &yup1 = field.yup(1); + EXPECT_NE(&field, &yup1); + EXPECT_NE(&yup, &yup1); + auto &ydown1 = field.ydown(1); + EXPECT_NE(&field, &ydown1); + EXPECT_NE(&ydown, &ydown1); + +#if CHECK > 1 + EXPECT_THROW(field.yup(2), BoutException); +#endif +} + TEST_F(Field3DTest, Ynext) { Field3D field; From a972f59b980a5ad91b2fac9b46bdb99f2c2adfe7 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Wed, 7 Nov 2018 17:27:15 +0000 Subject: [PATCH 02/34] Don't return *this from yup/ydown if there are no parallel slices This changes the behaviour of mergeYupYdown: afterwards, the parallel slices are not valid --- include/field3d.hxx | 12 ------------ tests/unit/field/test_field3d.cxx | 19 +++++-------------- 2 files changed, 5 insertions(+), 26 deletions(-) diff --git a/include/field3d.hxx b/include/field3d.hxx index 3a774e2d48..b662225860 100644 --- a/include/field3d.hxx +++ b/include/field3d.hxx @@ -233,35 +233,23 @@ class Field3D : public Field, public FieldData { /// Return reference to yup field Field3D &yup(std::vector::size_type index = 0) { - if (yup_fields.empty()) { - return *this; - } ASSERT2(index < yup_fields.size()); // Check for communicate return yup_fields[index]; } /// Return const reference to yup field const Field3D &yup(std::vector::size_type index = 0) const { - if (yup_fields.empty()) { - return *this; - } ASSERT2(index < yup_fields.size()); return yup_fields[index]; } /// Return reference to ydown field Field3D &ydown(std::vector::size_type index = 0) { - if (ydown_fields.empty()) { - return *this; - } ASSERT2(index < ydown_fields.size()); return ydown_fields[index]; } /// Return const reference to ydown field const Field3D &ydown(std::vector::size_type index = 0) const { - if (ydown_fields.empty()) { - return *this; - } ASSERT2(index < ydown_fields.size()); return ydown_fields[index]; } diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index 775a383621..842c279a61 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -245,18 +245,11 @@ TEST_F(Field3DTest, MergeYupYDown) { EXPECT_FALSE(field.hasYupYdown()); - auto& yup = field.yup(); - EXPECT_EQ(&field, &yup); - auto& ydown = field.ydown(); - EXPECT_EQ(&field, &ydown); + EXPECT_THROW(field.yup(), BoutException); + EXPECT_THROW(field.ydown(), BoutException); // Should be able to merge again without any problems - field.mergeYupYdown(); - - auto& yup2 = field.yup(); - EXPECT_EQ(&field, &yup2); - auto& ydown2 = field.ydown(); - EXPECT_EQ(&field, &ydown2); + EXPECT_NO_THROW(field.mergeYupYdown()); } TEST_F(Field3DTest, SplitThenMergeYupYDown) { @@ -272,10 +265,8 @@ TEST_F(Field3DTest, SplitThenMergeYupYDown) { field.mergeYupYdown(); - auto& yup2 = field.yup(); - EXPECT_EQ(&field, &yup2); - auto& ydown2 = field.ydown(); - EXPECT_EQ(&field, &ydown2); + EXPECT_THROW(field.yup(), BoutException); + EXPECT_THROW(field.ydown(), BoutException); } TEST_F(Field3DTest, MultipleYupYdown) { From a66fb8496dda744211f308b60c795873c72de092 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Wed, 21 Nov 2018 10:33:23 +0000 Subject: [PATCH 03/34] Don't reserve yup/ydown_fields vectors now that we have move ctors --- src/field/field3d.cxx | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 02740a12c5..c59aa6f7db 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -53,9 +53,6 @@ Field3D::Field3D(Mesh* localmesh) : Field(localmesh) { nx = fieldmesh->LocalNx; ny = fieldmesh->LocalNy; nz = fieldmesh->LocalNz; - - yup_fields.reserve(fieldmesh->ystart); - ydown_fields.reserve(fieldmesh->ystart); } } @@ -73,9 +70,6 @@ Field3D::Field3D(const Field3D& f) : Field(f.fieldmesh), data(f.data) { nx = fieldmesh->LocalNx; ny = fieldmesh->LocalNy; nz = fieldmesh->LocalNz; - - yup_fields.reserve(fieldmesh->ystart); - ydown_fields.reserve(fieldmesh->ystart); } location = f.location; From 21d400cc5997406bed37d688d99111f12c99217d Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Wed, 21 Nov 2018 15:14:38 +0000 Subject: [PATCH 04/34] Add implementation for Identity parallel transform + tests --- include/bout/paralleltransform.hxx | 2 +- src/mesh/parallel/identity.cxx | 11 +++++++ src/mesh/parallel/makefile | 2 +- tests/unit/mesh/test_paralleltransform.cxx | 37 ++++++++++++++++++++++ 4 files changed, 50 insertions(+), 2 deletions(-) create mode 100644 src/mesh/parallel/identity.cxx create mode 100644 tests/unit/mesh/test_paralleltransform.cxx diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx index a25540354a..9308b8dad9 100644 --- a/include/bout/paralleltransform.hxx +++ b/include/bout/paralleltransform.hxx @@ -59,7 +59,7 @@ public: * Merges the yup and ydown() fields of f, so that * f.yup() = f.ydown() = f */ - void calcYUpDown(Field3D &f) override {f.mergeYupYdown();} + void calcYUpDown(Field3D &f) override; /*! * The field is already aligned in Y, so this diff --git a/src/mesh/parallel/identity.cxx b/src/mesh/parallel/identity.cxx new file mode 100644 index 0000000000..437d04178b --- /dev/null +++ b/src/mesh/parallel/identity.cxx @@ -0,0 +1,11 @@ +#include "bout/paralleltransform.hxx" +#include "bout/mesh.hxx" + +void ParallelTransformIdentity::calcYUpDown(Field3D& f) { + f.splitYupYdown(); + + for (int i = 0; i < f.getMesh()->ystart; ++i) { + f.yup(i) = f; + f.ydown(i) = f; + } +} diff --git a/src/mesh/parallel/makefile b/src/mesh/parallel/makefile index 9c4817b008..2e522318de 100644 --- a/src/mesh/parallel/makefile +++ b/src/mesh/parallel/makefile @@ -2,7 +2,7 @@ BOUT_TOP = ../../.. DIRS = -SOURCEC = shiftedmetric.cxx fci.cxx +SOURCEC = shiftedmetric.cxx fci.cxx identity.cxx TARGET = lib include $(BOUT_TOP)/make.config diff --git a/tests/unit/mesh/test_paralleltransform.cxx b/tests/unit/mesh/test_paralleltransform.cxx new file mode 100644 index 0000000000..6550243b2a --- /dev/null +++ b/tests/unit/mesh/test_paralleltransform.cxx @@ -0,0 +1,37 @@ +#include "gtest/gtest.h" + +#include "test_extras.hxx" +#include "bout/paralleltransform.hxx" + +extern Mesh* mesh; + +using ParallelTransformTest = FakeMeshFixture; + +TEST_F(ParallelTransformTest, IdentityCalcYUpDown) { + + ParallelTransformIdentity transform{}; + + Field3D field{1.0}; + + transform.calcYUpDown(field); + + EXPECT_TRUE(IsField3DEqualBoutReal(field.yup(), 1.0)); + EXPECT_TRUE(IsField3DEqualBoutReal(field.ydown(), 1.0)); +} + +TEST_F(ParallelTransformTest, IdentityCalcYUpDownTwoSlices) { + + ParallelTransformIdentity transform{}; + + mesh->ystart = 2; + + Field3D field{1.0}; + + transform.calcYUpDown(field); + + EXPECT_TRUE(IsField3DEqualBoutReal(field.yup(0), 1.0)); + EXPECT_TRUE(IsField3DEqualBoutReal(field.yup(1), 1.0)); + + EXPECT_TRUE(IsField3DEqualBoutReal(field.ydown(0), 1.0)); + EXPECT_TRUE(IsField3DEqualBoutReal(field.ydown(1), 1.0)); +} From d1449bbcb0f8f920d850852a2be420d8a793608e Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 14 Jan 2019 11:16:37 +0000 Subject: [PATCH 05/34] Implement semi-optimisation for ShiftedMetric::calcYUpDown Forward-FFT whole field, then inverse-FFT for each parallel slice. Will make it easier to generalise to multiple parallel slices --- include/bout/paralleltransform.hxx | 7 ++- src/mesh/parallel/shiftedmetric.cxx | 69 +++++++++++++++++++++-------- 2 files changed, 55 insertions(+), 21 deletions(-) diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx index 34fabd061f..38b1604703 100644 --- a/include/bout/paralleltransform.hxx +++ b/include/bout/paralleltransform.hxx @@ -125,7 +125,7 @@ public: } /// A 3D array, implemented as nested vectors - using arr3Dvec = std::vector>>; + using arr3Dvec = std::vector>>; private: Mesh &mesh; ///< The mesh this paralleltransform is part of @@ -185,11 +185,14 @@ 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) const; + void shiftZ(const BoutReal *in, const Array &phs, BoutReal *out) const; /// Calculate and store the phases for to/from field aligned and for /// the parallel slices using zShift void cachePhases(); + + std::vector shiftZ(const Field3D& f, const std::vector& phases, + const std::vector& y_offsets) const; }; diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx index 55eee670bd..6da967135b 100644 --- a/src/mesh/parallel/shiftedmetric.cxx +++ b/src/mesh/parallel/shiftedmetric.cxx @@ -109,28 +109,15 @@ void ShiftedMetric::cachePhases() { /*! * Calculate the Y up and down fields */ -void ShiftedMetric::calcYUpDown(Field3D &f) { +void ShiftedMetric::calcYUpDown(Field3D& f) { f.splitYupYdown(); - - Field3D& yup = f.yup(); - yup.allocate(); - - for(int jx=0;jx& phs, +void ShiftedMetric::shiftZ(const BoutReal* in, const Array& phs, BoutReal* out) const { int nmodes = mesh.LocalNz / 2 + 1; @@ -187,6 +174,50 @@ void ShiftedMetric::shiftZ(const BoutReal* in, const std::vector& phs, irfft(&cmplx[0], mesh.LocalNz, out); // Reverse FFT } +std::vector ShiftedMetric::shiftZ(const Field3D& f, + const std::vector& phases, + const std::vector& y_offsets) const { + + ASSERT1(phases.size() == y_offsets.size()); + + const int nmodes = mesh.LocalNz / 2 + 1; + + // FFT in Z of input field at each (x, y) point + arr3Dvec f_fft(mesh.LocalNx, + std::vector>(mesh.LocalNy, Array(nmodes))); + + std::vector results{}; + + for (int jx = 0; jx < mesh.LocalNx; jx++) { + for (int jy = 0; jy < mesh.LocalNy; jy++) { + rfft(f(jx, jy), mesh.LocalNz, &f_fft[jx][jy][0]); + } + } + + for (std::size_t i = 0; i < phases.size(); ++i) { + + results.emplace_back(&mesh); + results[i].allocate(); + results[i].setLocation(f.getLocation()); + + for (int jx = 0; jx < mesh.LocalNx; jx++) { + for (int jy = mesh.ystart; jy <= mesh.yend; jy++) { + + Array shifted_temp(f_fft[jx][jy + y_offsets[i]]); + shifted_temp.ensureUnique(); + + for (int jz = 1; jz < nmodes; ++jz) { + shifted_temp[jz] *= phases[i][jx][jy][jz]; + } + + irfft(&shifted_temp[0], mesh.LocalNz, results[i](jx, jy + y_offsets[i])); + } + } + } + + return results; +} + //Old approach retained so we can still specify a general zShift const Field3D ShiftedMetric::shiftZ(const Field3D &f, const Field2D &zangle) const { ASSERT1(&mesh == f.getMesh()); From f7fc1d861b7ebeb530f2dd2ddf4d3a1fc1952bc5 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 14 Jan 2019 15:51:57 +0000 Subject: [PATCH 06/34] Implement multiple parallel slices for ShiftedMetric --- include/bout/paralleltransform.hxx | 23 ++- src/mesh/parallel/shiftedmetric.cxx | 112 ++++++++------ .../unit/mesh/parallel/test_shiftedmetric.cxx | 146 ++++++++++++------ 3 files changed, 191 insertions(+), 90 deletions(-) diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx index 38b1604703..5ed14cd4f2 100644 --- a/include/bout/paralleltransform.hxx +++ b/include/bout/paralleltransform.hxx @@ -137,6 +137,20 @@ private: /// Cache of phase shifts for transforming from field-aligned coordinates to X-Z orthogonal coordinates arr3Dvec fromAlignedPhs; + /// Helper POD for parallel slice phase shifts + struct ParallelSlicePhase { + arr3Dvec phase_shift; + int y_offset; + }; + + /// Cache of phase shifts for the parallel slices. Slices are stored + /// in the following order: + /// {+1, ..., +n, -1, ..., -n} + /// slice[i] stores offset i+1 + /// slice[2*i + 1] stores offset -(i+1) + /// where i goes from 0 to (n-1), with n the number of y guard cells + std::vector parallel_slice_phases; + arr3Dvec yupPhs; ///< Cache of phase shifts for calculating yup fields arr3Dvec ydownPhs; ///< Cache of phase shifts for calculating ydown fields @@ -191,8 +205,13 @@ private: /// the parallel slices using zShift void cachePhases(); - std::vector shiftZ(const Field3D& f, const std::vector& phases, - const std::vector& y_offsets) const; + /// Shift a 3D field \p f in Z to all the parallel slices in \p phases + /// + /// @param[in] f The field to shift + /// @param[in] phases The phase and offset information for each parallel slice + /// @return The shifted parallel slices + std::vector shiftZ(const Field3D& f, + const std::vector& phases) const; }; diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx index 6da967135b..7abf316c7e 100644 --- a/src/mesh/parallel/shiftedmetric.cxx +++ b/src/mesh/parallel/shiftedmetric.cxx @@ -58,21 +58,13 @@ void ShiftedMetric::cachePhases() { fromAlignedPhs.resize(mesh.LocalNx); toAlignedPhs.resize(mesh.LocalNx); - yupPhs.resize(mesh.LocalNx); - ydownPhs.resize(mesh.LocalNx); - for (int jx = 0; jx < mesh.LocalNx; jx++) { fromAlignedPhs[jx].resize(mesh.LocalNy); toAlignedPhs[jx].resize(mesh.LocalNy); - yupPhs[jx].resize(mesh.LocalNy); - ydownPhs[jx].resize(mesh.LocalNy); for (int jy = 0; jy < mesh.LocalNy; jy++) { fromAlignedPhs[jx][jy].resize(nmodes); toAlignedPhs[jx][jy].resize(nmodes); - - yupPhs[jx][jy].resize(nmodes); - ydownPhs[jx][jy].resize(nmodes); } } @@ -89,33 +81,46 @@ void ShiftedMetric::cachePhases() { } } - // Yup/Ydown phases -- note we don't shift in the boundaries/guards - for (int jx = 0; jx < mesh.LocalNx; jx++) { - for (int jy = mesh.ystart; jy <= mesh.yend; jy++) { - BoutReal yupShift = zShift(jx, jy) - zShift(jx, jy + 1); - BoutReal ydownShift = zShift(jx, jy) - zShift(jx, jy - 1); - - for (int jz = 0; jz < nmodes; jz++) { - BoutReal kwave = jz * 2.0 * PI / zlength; // wave number is 1/[rad] + // Allocate space for parallel slice caches: y-guard cells in each + // direction + parallel_slice_phases.resize(mesh.ystart * 2); + + // Careful with the indices/offsets! Offsets are 1-indexed (as 0 + // would be the original slice), and Mesh::ystart is the number of + // guard cells. The parallel slice vector stores the offsets as + // {+1, ..., +n, -1, ..., -n} + // Once parallel_slice_phases is initialised though, each element + // stores its phase and offset, so we don't need to faff about after + // this + for (int i = 0; i < mesh.ystart; ++i) { + parallel_slice_phases[i].phase_shift = + arr3Dvec(mesh.LocalNx, + std::vector>(mesh.LocalNy, Array(nmodes))); + parallel_slice_phases[i].y_offset = i + 1; - yupPhs[jx][jy][jz] = dcomplex(cos(kwave * yupShift), -sin(kwave * yupShift)); - ydownPhs[jx][jy][jz] = - dcomplex(cos(kwave * ydownShift), -sin(kwave * ydownShift)); - } - } + parallel_slice_phases[mesh.ystart + i].phase_shift = + arr3Dvec(mesh.LocalNx, + std::vector>(mesh.LocalNy, Array(nmodes))); + parallel_slice_phases[mesh.ystart + i].y_offset = -(i + 1); } -} -/*! - * Calculate the Y up and down fields - */ -void ShiftedMetric::calcYUpDown(Field3D& f) { - f.splitYupYdown(); + // Parallel slice phases -- note we don't shift in the boundaries/guards + for (auto& slice : parallel_slice_phases) { + for (int jx = 0; jx < mesh.LocalNx; jx++) { + for (int jy = mesh.ystart; jy <= mesh.yend; jy++) { + + BoutReal slice_shift = zShift(jx, jy) - zShift(jx, jy + slice.y_offset); - auto results = shiftZ(f, {yupPhs, ydownPhs}, {+1, -1}); + for (int jz = 0; jz < nmodes; jz++) { + // wave number is 1/[rad] + BoutReal kwave = jz * 2.0 * PI / zlength; - f.yup() = results[0]; - f.ydown() = results[1]; + slice.phase_shift[jx][jy][jz] = + dcomplex(cos(kwave * slice_shift), -sin(kwave * slice_shift)); + } + } + } + } } /*! @@ -174,11 +179,23 @@ void ShiftedMetric::shiftZ(const BoutReal* in, const Array& phs, irfft(&cmplx[0], mesh.LocalNz, out); // Reverse FFT } -std::vector ShiftedMetric::shiftZ(const Field3D& f, - const std::vector& phases, - const std::vector& y_offsets) const { - ASSERT1(phases.size() == y_offsets.size()); +void ShiftedMetric::calcYUpDown(Field3D& f) { + + auto results = shiftZ(f, parallel_slice_phases); + + ASSERT3(results.size() == parallel_slice_phases.size()); + + f.splitYupYdown(); + + for (std::size_t i = 0; i < results.size(); ++i) { + f.ynext(parallel_slice_phases[i].y_offset) = std::move(results[i]); + } +} + +std::vector +ShiftedMetric::shiftZ(const Field3D& f, + const std::vector& phases) const { const int nmodes = mesh.LocalNz / 2 + 1; @@ -186,31 +203,40 @@ std::vector ShiftedMetric::shiftZ(const Field3D& f, arr3Dvec f_fft(mesh.LocalNx, std::vector>(mesh.LocalNy, Array(nmodes))); - std::vector results{}; - for (int jx = 0; jx < mesh.LocalNx; jx++) { for (int jy = 0; jy < mesh.LocalNy; jy++) { rfft(f(jx, jy), mesh.LocalNz, &f_fft[jx][jy][0]); } } - for (std::size_t i = 0; i < phases.size(); ++i) { + std::vector results{}; + + for (auto& phase : phases) { + // In C++17 std::vector::emplace_back returns a reference, which + // would be very useful here! - results.emplace_back(&mesh); - results[i].allocate(); - results[i].setLocation(f.getLocation()); + // FIXME: initialisation to -1 to avoid checkData choking on the + // uninitialised regions in assignment into the Field parallel + // slices in calcYUpDown + results.emplace_back(-1.0, &mesh); + auto& current_result = results.back(); + // FIXME: uncomment the following after fixing Field3D::operator= + // current_result.allocate(); + current_result.setLocation(f.getLocation()); for (int jx = 0; jx < mesh.LocalNx; jx++) { for (int jy = mesh.ystart; jy <= mesh.yend; jy++) { - Array shifted_temp(f_fft[jx][jy + y_offsets[i]]); + // Deep copy the FFT'd field + Array shifted_temp(f_fft[jx][jy + phase.y_offset]); shifted_temp.ensureUnique(); for (int jz = 1; jz < nmodes; ++jz) { - shifted_temp[jz] *= phases[i][jx][jy][jz]; + shifted_temp[jz] *= phase.phase_shift[jx][jy][jz]; } - irfft(&shifted_temp[0], mesh.LocalNz, results[i](jx, jy + y_offsets[i])); + irfft(shifted_temp.begin(), mesh.LocalNz, + current_result(jx, jy + phase.y_offset)); } } } diff --git a/tests/unit/mesh/parallel/test_shiftedmetric.cxx b/tests/unit/mesh/parallel/test_shiftedmetric.cxx index 4d8c2885a7..a03cd9a593 100644 --- a/tests/unit/mesh/parallel/test_shiftedmetric.cxx +++ b/tests/unit/mesh/parallel/test_shiftedmetric.cxx @@ -136,20 +136,29 @@ TEST_F(ShiftedMetricTest, FromFieldAligned) { } TEST_F(ShiftedMetricTest, CalcYUpDown) { + // Use two y-guards to test multiple parallel slices + mesh->ystart = 2; + mesh->yend = mesh->LocalNy - 3; + + // We don't shift in the guard cells, and the parallel slices are + // stored offset in y, therefore we need to make new regions that we + // can compare the expected and actual outputs over output_info.disable(); - auto region_yup = mesh->getRegion("RGN_NOY"); - region_yup.periodicShift(ShiftedMetricTest::nz, - ShiftedMetricTest::ny * ShiftedMetricTest::nz); - mesh->addRegion("RGN_YUP", region_yup); - - auto region_ydown = mesh->getRegion("RGN_NOY"); - region_ydown.periodicShift(-ShiftedMetricTest::nz, - ShiftedMetricTest::ny * ShiftedMetricTest::nz); - mesh->addRegion("RGN_YDOWN", region_ydown); + mesh->addRegion3D("RGN_YUP", + Region(0, mesh->LocalNx - 1, mesh->ystart + 1, mesh->yend + 1, + 0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz)); + mesh->addRegion3D("RGN_YUP2", + Region(0, mesh->LocalNx - 1, mesh->ystart + 2, mesh->yend + 2, + 0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz)); + + mesh->addRegion3D("RGN_YDOWN", + Region(0, mesh->LocalNx - 1, mesh->ystart - 1, mesh->yend - 1, + 0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz)); + mesh->addRegion3D("RGN_YDOWN2", + Region(0, mesh->LocalNx - 1, mesh->ystart - 2, mesh->yend - 2, + 0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz)); output_info.enable(); - ShiftedMetric shifted{*mesh, zShift}; - Field3D input{mesh}; fillField(input, {{{1., 2., 3., 4., 5., 6., 7.}, @@ -170,48 +179,95 @@ TEST_F(ShiftedMetricTest, CalcYUpDown) { {1., 2., 3., 4., 5., 6., 7.}, {1., 2., 3., 4., 5., 6., 7.}}}); + // Actual interesting bit here! + ShiftedMetric shifted{*mesh, zShift}; shifted.calcYUpDown(input); - Field3D expected_up{mesh}; - - fillField(expected_up, {{{0., 0., 0., 0., 0., 0., 0.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}}, - - {{0., 0., 0., 0., 0., 0., 0.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}}, - - {{0., 0., 0., 0., 0., 0., 0.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}, - {2., 3., 4., 5., 6., 7., 1.}}}); + // Expected output values - Field3D expected_down{mesh}; + Field3D expected_up_1{mesh}; - fillField(expected_down, {{{7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, + // Note: here zeroes are for values we don't expect to read + fillField(expected_up_1, {{{0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {2., 3., 4., 5., 6., 7., 1.}, {0., 0., 0., 0., 0., 0., 0.}}, - {{7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, + {{0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {2., 3., 4., 5., 6., 7., 1.}, {0., 0., 0., 0., 0., 0., 0.}}, - {{7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, - {7., 1., 2., 3., 4., 5., 6.}, + {{0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {2., 3., 4., 5., 6., 7., 1.}, {0., 0., 0., 0., 0., 0., 0.}}}); - EXPECT_TRUE(IsField3DEqualField3D(input.yup(), expected_up, "RGN_YUP")); - EXPECT_TRUE(IsField3DEqualField3D(input.ydown(), expected_down, "RGN_YDOWN")); + Field3D expected_up_2{mesh}; + + fillField(expected_up_2, {{{0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {3., 4., 5., 6., 7., 1., 2.}}, + + {{0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {3., 4., 5., 6., 7., 1., 2.}}, + + {{0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {3., 4., 5., 6., 7., 1., 2.}}}); + + Field3D expected_down_1{mesh}; + + fillField(expected_down_1, {{{0., 0., 0., 0., 0., 0., 0.}, + {7., 1., 2., 3., 4., 5., 6.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}}, + + {{0., 0., 0., 0., 0., 0., 0.}, + {7., 1., 2., 3., 4., 5., 6.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}}, + + {{0., 0., 0., 0., 0., 0., 0.}, + {7., 1., 2., 3., 4., 5., 6.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}}}); + + Field3D expected_down2{mesh}; + + fillField(expected_down2, {{{6., 7., 1., 2., 3., 4., 5.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}}, + + {{6., 7., 1., 2., 3., 4., 5.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}}, + + {{6., 7., 1., 2., 3., 4., 5.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0., 0.}}}); + + EXPECT_TRUE(IsField3DEqualField3D(input.ynext(1), expected_up_1, "RGN_YUP")); + EXPECT_TRUE(IsField3DEqualField3D(input.ynext(2), expected_up_2, "RGN_YUP2")); + EXPECT_TRUE(IsField3DEqualField3D(input.ynext(-1), expected_down_1, "RGN_YDOWN")); + EXPECT_TRUE(IsField3DEqualField3D(input.ynext(-2), expected_down2, "RGN_YDOWN2")); } From 4064d23e6fd1daa223fbb87361ba95df32571723 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Wed, 21 Nov 2018 15:37:29 +0000 Subject: [PATCH 07/34] Guard tests that might not throw at low CHECK --- tests/unit/field/test_field3d.cxx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index 2b21b284a5..c6fe34f8c3 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -218,8 +218,10 @@ TEST_F(Field3DTest, MergeYupYDown) { EXPECT_FALSE(field.hasYupYdown()); +#if CHECK > 2 EXPECT_THROW(field.yup(), BoutException); EXPECT_THROW(field.ydown(), BoutException); +#endif // Should be able to merge again without any problems EXPECT_NO_THROW(field.mergeYupYdown()); @@ -238,8 +240,10 @@ TEST_F(Field3DTest, SplitThenMergeYupYDown) { field.mergeYupYdown(); +#if CHECK > 2 EXPECT_THROW(field.yup(), BoutException); EXPECT_THROW(field.ydown(), BoutException); +#endif } TEST_F(Field3DTest, MultipleYupYdown) { From 901b8c119e21f92ab0fafc6485ecef4b73d0b6eb Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 14 Jan 2019 16:38:34 +0000 Subject: [PATCH 08/34] Fix test-smooth for multiple parallel slices --- tests/integrated/test-smooth/test_smooth.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/integrated/test-smooth/test_smooth.cxx b/tests/integrated/test-smooth/test_smooth.cxx index f1246d0fed..15f7c656ec 100644 --- a/tests/integrated/test-smooth/test_smooth.cxx +++ b/tests/integrated/test-smooth/test_smooth.cxx @@ -17,7 +17,8 @@ int main(int argc, char **argv) { Field2D input2d = f.create2D("1 + sin(2*y)"); Field3D input3d = f.create3D("gauss(x-0.5,0.2)*gauss(y-pi)*sin(3*y - z)"); - input3d.mergeYupYdown(); + input3d.splitYupYdown(); + mesh->getParallelTransform().calcYUpDown(input3d); SAVE_ONCE2(input2d, input3d); From a215003ade26096122280c162fb1b18306f55ebb Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 14 Jan 2019 16:42:16 +0000 Subject: [PATCH 09/34] Let interp_to in y work with multiple parallel slices --- include/interpolation.hxx | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/include/interpolation.hxx b/include/interpolation.hxx index ad72870878..038285f824 100644 --- a/include/interpolation.hxx +++ b/include/interpolation.hxx @@ -118,15 +118,7 @@ const T interp_to(const T& var, CELL_LOC loc, REGION region = RGN_ALL) { // At least 2 boundary cells needed for interpolation in y-direction ASSERT0(fieldmesh->ystart >= 2); - if (var.hasYupYdown() && ((&var.yup() != &var) || (&var.ydown() != &var))) { - // 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."); - + if (var.hasYupYdown()) { if ((location == CELL_CENTRE) && (loc == CELL_YLOW)) { // C2L BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { // Producing a stencil centred around a lower X value From 1b0d149ca5afb394bf15cda06d45a3db612b48fe Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 14 Jan 2019 17:29:52 +0000 Subject: [PATCH 10/34] Enable higher-order derivatives with multiple parallel slices --- include/bout/index_derivs.hxx | 31 ++--------------------- include/bout/index_derivs_interface.hxx | 33 ++++++++++--------------- src/mesh/difops.cxx | 4 +-- 3 files changed, 17 insertions(+), 51 deletions(-) diff --git a/include/bout/index_derivs.hxx b/include/bout/index_derivs.hxx index 6a4c2f5293..af453a5f83 100644 --- a/include/bout/index_derivs.hxx +++ b/include/bout/index_derivs.hxx @@ -619,7 +619,7 @@ struct registerMethod { ///////////////////////////////////////////////////////////////////////////////// produceCombinations, + WRAP_ENUM(DIRECTION, YOrthogonal), WRAP_ENUM(DIRECTION, Z)>, Set, Set, TypeContainer>, Set< @@ -640,22 +640,8 @@ produceCombinations>> registerDerivatives(registerMethod{}); -produceCombinations, Set, - Set, TypeContainer>, - Set< - // Standard - DerivativeType, DerivativeType, - // Standard 2nd order - DerivativeType, - // Standard 4th order - // Upwind - DerivativeType, DerivativeType, - // Flux - DerivativeType>> - registerDerivativesYOrtho(registerMethod{}); - produceCombinations, + WRAP_ENUM(DIRECTION, YOrthogonal), WRAP_ENUM(DIRECTION, Z)>, Set, Set, TypeContainer>, Set< @@ -670,19 +656,6 @@ produceCombinations>> registerStaggeredDerivatives(registerMethod{}); -produceCombinations, - Set, - Set, TypeContainer>, - Set< - // Standard - DerivativeType, - // Standard 2nd order - // Upwind - DerivativeType, DerivativeType, - // Flux - DerivativeType>> - registerStaggeredDerivativesYOrtho(registerMethod{}); - class FFTDerivativeType { public: template diff --git a/include/bout/index_derivs_interface.hxx b/include/bout/index_derivs_interface.hxx index 154a6b9fd6..4ebdc2ac59 100644 --- a/include/bout/index_derivs_interface.hxx +++ b/include/bout/index_derivs_interface.hxx @@ -207,14 +207,13 @@ template T DDY(const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", REGION region = RGN_NOBNDRY) { AUTO_TRACE(); - if (std::is_base_of::value && f.hasYupYdown() - && ((&f.yup() != &f) || (&f.ydown() != &f))) { + if (std::is_base_of::value && f.hasYupYdown()) { return standardDerivative(f, outloc, method, region); } else { const T f_aligned = f.getMesh()->toFieldAligned(f); - T result = - standardDerivative(f_aligned, outloc, method, region); + T result = standardDerivative(f_aligned, outloc, + method, region); return f.getMesh()->fromFieldAligned(result); } } @@ -223,14 +222,13 @@ template T D2DY2(const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", REGION region = RGN_NOBNDRY) { AUTO_TRACE(); - if (std::is_base_of::value && f.hasYupYdown() - && ((&f.yup() != &f) || (&f.ydown() != &f))) { + if (std::is_base_of::value && f.hasYupYdown()) { return standardDerivative( f, outloc, method, region); } else { const T f_aligned = f.getMesh()->toFieldAligned(f); - T result = standardDerivative(f_aligned, outloc, - method, region); + T result = standardDerivative( + f_aligned, outloc, method, region); return f.getMesh()->fromFieldAligned(result); } } @@ -239,14 +237,13 @@ template T D4DY4(const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", REGION region = RGN_NOBNDRY) { AUTO_TRACE(); - if (std::is_base_of::value && f.hasYupYdown() - && ((&f.yup() != &f) || (&f.ydown() != &f))) { + if (std::is_base_of::value && f.hasYupYdown()) { return standardDerivative( f, outloc, method, region); } else { const T f_aligned = f.getMesh()->toFieldAligned(f); - T result = standardDerivative(f_aligned, outloc, - method, region); + T result = standardDerivative( + f_aligned, outloc, method, region); return f.getMesh()->fromFieldAligned(result); } } @@ -312,10 +309,8 @@ template T VDDY(const T& vel, const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", REGION region = RGN_NOBNDRY) { AUTO_TRACE(); - bool fHasParallelSlices = (std::is_base_of::value && f.hasYupYdown() - && ((&f.yup() != &f) || (&f.ydown() != &f))); - bool velHasParallelSlices = (std::is_base_of::value && vel.hasYupYdown() - && ((&vel.yup() != &vel) || (&vel.ydown() != &vel))); + bool fHasParallelSlices = (std::is_base_of::value && f.hasYupYdown()); + bool velHasParallelSlices = (std::is_base_of::value && vel.hasYupYdown()); if (fHasParallelSlices && velHasParallelSlices) { return flowDerivative(vel, f, outloc, method, region); @@ -332,10 +327,8 @@ template T FDDY(const T& vel, const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", REGION region = RGN_NOBNDRY) { AUTO_TRACE(); - bool fHasParallelSlices = (std::is_base_of::value && f.hasYupYdown() - && ((&f.yup() != &f) || (&f.ydown() != &f))); - bool velHasParallelSlices = (std::is_base_of::value && vel.hasYupYdown() - && ((&vel.yup() != &vel) || (&vel.ydown() != &vel))); + bool fHasParallelSlices = (std::is_base_of::value && f.hasYupYdown()); + bool velHasParallelSlices = (std::is_base_of::value && vel.hasYupYdown()); if (fHasParallelSlices && velHasParallelSlices) { return flowDerivative(vel, f, outloc, method, region); diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index c924ae11e0..fa88ec6868 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -340,8 +340,8 @@ const Field3D Vpar_Grad_par_LCtoC(const Field3D &v, const Field3D &f, REGION reg result.allocate(); - bool vUseUpDown = (v.hasYupYdown() && ((&v.yup() != &v) || (&v.ydown() != &v))); - bool fUseUpDown = (f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f))); + bool vUseUpDown = v.hasYupYdown(); + bool fUseUpDown = f.hasYupYdown(); if (vUseUpDown && fUseUpDown) { // Both v and f have up/down fields From 4d065b3a022579282b21ed23a414d3ed9cfa3f06 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 15 Jan 2019 11:57:44 +0000 Subject: [PATCH 11/34] Fix bug with new ShiftedMetric implementation: Arrays weren't unique Tighten up the unit tests to catch this --- src/mesh/parallel/shiftedmetric.cxx | 8 +- .../unit/mesh/parallel/test_shiftedmetric.cxx | 342 +++++++++--------- 2 files changed, 182 insertions(+), 168 deletions(-) diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx index 7abf316c7e..cf6124bd29 100644 --- a/src/mesh/parallel/shiftedmetric.cxx +++ b/src/mesh/parallel/shiftedmetric.cxx @@ -93,11 +93,15 @@ void ShiftedMetric::cachePhases() { // stores its phase and offset, so we don't need to faff about after // this for (int i = 0; i < mesh.ystart; ++i) { + // NOTE: std::vector constructor here takes a **copy** of the + // Array! We *must* call `Array::ensureUnique` on each element + // before using it! parallel_slice_phases[i].phase_shift = arr3Dvec(mesh.LocalNx, std::vector>(mesh.LocalNy, Array(nmodes))); parallel_slice_phases[i].y_offset = i + 1; + // Backwards parallel slices parallel_slice_phases[mesh.ystart + i].phase_shift = arr3Dvec(mesh.LocalNx, std::vector>(mesh.LocalNy, Array(nmodes))); @@ -109,6 +113,7 @@ void ShiftedMetric::cachePhases() { for (int jx = 0; jx < mesh.LocalNx; jx++) { for (int jy = mesh.ystart; jy <= mesh.yend; jy++) { + slice.phase_shift[jx][jy].ensureUnique(); BoutReal slice_shift = zShift(jx, jy) - zShift(jx, jy + slice.y_offset); for (int jz = 0; jz < nmodes; jz++) { @@ -205,7 +210,8 @@ ShiftedMetric::shiftZ(const Field3D& f, for (int jx = 0; jx < mesh.LocalNx; jx++) { for (int jy = 0; jy < mesh.LocalNy; jy++) { - rfft(f(jx, jy), mesh.LocalNz, &f_fft[jx][jy][0]); + f_fft[jx][jy].ensureUnique(); + rfft(f(jx, jy), mesh.LocalNz, f_fft[jx][jy].begin()); } } diff --git a/tests/unit/mesh/parallel/test_shiftedmetric.cxx b/tests/unit/mesh/parallel/test_shiftedmetric.cxx index a03cd9a593..8a806506c3 100644 --- a/tests/unit/mesh/parallel/test_shiftedmetric.cxx +++ b/tests/unit/mesh/parallel/test_shiftedmetric.cxx @@ -20,7 +20,37 @@ class ShiftedMetricTest : public ::testing::Test { zShift = Field2D{mesh}; - fillField(zShift, {{1., 2., 3., 4., 5.}, {1., 2., 3., 4., 5.}, {1., 2., 3., 4., 5.}}); + fillField(zShift, {{1., 2., 3., 4., 5., 6., 7.}, + {2., 4., 6., 8., 10., 12., 14.}, + {3., 6., 9., 12., 15., 18., 21.}}); + + Field3D input_temp{mesh}; + + fillField(input_temp, {{{1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}}, + + {{1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}}, + + {{1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}, + {1., 2., 3., 4., 5.}}}); + + input = std::move(input_temp); dynamic_cast(mesh)->setCoordinates(std::make_shared( mesh, Field2D{1.0}, Field2D{1.0}, BoutReal{1.0}, Field2D{1.0}, Field2D{0.0}, @@ -35,100 +65,74 @@ class ShiftedMetricTest : public ::testing::Test { } static constexpr int nx = 3; - static constexpr int ny = 5; - static constexpr int nz = 7; + static constexpr int ny = 7; + static constexpr int nz = 5; Field2D zShift; + Field3D input; }; TEST_F(ShiftedMetricTest, ToFieldAligned) { ShiftedMetric shifted{*mesh, zShift}; - Field3D input{mesh}; - - fillField(input, {{{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}, - - {{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}, - - {{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}}); - Field3D expected{mesh}; - fillField(expected, {{{2., 3., 4., 5., 6., 7., 1.}, - {3., 4., 5., 6., 7., 1., 2.}, - {4., 5., 6., 7., 1., 2., 3.}, - {5., 6., 7., 1., 2., 3., 4.}, - {6., 7., 1., 2., 3., 4., 5.}}, - - {{2., 3., 4., 5., 6., 7., 1.}, - {3., 4., 5., 6., 7., 1., 2.}, - {4., 5., 6., 7., 1., 2., 3.}, - {5., 6., 7., 1., 2., 3., 4.}, - {6., 7., 1., 2., 3., 4., 5.}}, - - {{2., 3., 4., 5., 6., 7., 1.}, - {3., 4., 5., 6., 7., 1., 2.}, - {4., 5., 6., 7., 1., 2., 3.}, - {5., 6., 7., 1., 2., 3., 4.}, - {6., 7., 1., 2., 3., 4., 5.}}}); - - EXPECT_TRUE(IsField3DEqualField3D(shifted.toFieldAligned(input), expected)); + fillField(expected, {{{2., 3., 4., 5., 1.}, + {3., 4., 5., 1., 2.}, + {4., 5., 1., 2., 3.}, + {5., 1., 2., 3., 4.}, + {1., 2., 3., 4., 5.}, + {2., 3., 4., 5., 1.}, + {3., 4., 5., 1., 2.}}, + + {{3., 4., 5., 1., 2.}, + {5., 1., 2., 3., 4.}, + {2., 3., 4., 5., 1.}, + {4., 5., 1., 2., 3.}, + {1., 2., 3., 4., 5.}, + {3., 4., 5., 1., 2.}, + {5., 1., 2., 3., 4.}}, + + {{4., 5., 1., 2., 3.}, + {2., 3., 4., 5., 1.}, + {5., 1., 2., 3., 4.}, + {3., 4., 5., 1., 2.}, + {1., 2., 3., 4., 5.}, + {4., 5., 1., 2., 3.}, + {2., 3., 4., 5., 1.}}}); + + EXPECT_TRUE( + IsField3DEqualField3D(shifted.toFieldAligned(input), expected, "RGN_ALL", 1.e-12)); } TEST_F(ShiftedMetricTest, FromFieldAligned) { ShiftedMetric shifted{*mesh, zShift}; - Field3D input{mesh}; - - fillField(input, {{{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}, - - {{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}, - - {{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}}); - Field3D expected{mesh}; - fillField(expected, {{{7., 1., 2., 3., 4., 5., 6.}, - {6., 7., 1., 2., 3., 4., 5.}, - {5., 6., 7., 1., 2., 3., 4.}, - {4., 5., 6., 7., 1., 2., 3.}, - {3., 4., 5., 6., 7., 1., 2.}}, - - {{7., 1., 2., 3., 4., 5., 6.}, - {6., 7., 1., 2., 3., 4., 5.}, - {5., 6., 7., 1., 2., 3., 4.}, - {4., 5., 6., 7., 1., 2., 3.}, - {3., 4., 5., 6., 7., 1., 2.}}, - - {{7., 1., 2., 3., 4., 5., 6.}, - {6., 7., 1., 2., 3., 4., 5.}, - {5., 6., 7., 1., 2., 3., 4.}, - {4., 5., 6., 7., 1., 2., 3.}, - {3., 4., 5., 6., 7., 1., 2.}}}); + fillField(expected, {{{5., 1., 2., 3., 4.}, + {4., 5., 1., 2., 3.}, + {3., 4., 5., 1., 2.}, + {2., 3., 4., 5., 1.}, + {1., 2., 3., 4., 5.}, + {5., 1., 2., 3., 4.}, + {4., 5., 1., 2., 3.}}, + + {{4., 5., 1., 2., 3.}, + {2., 3., 4., 5., 1.}, + {5., 1., 2., 3., 4.}, + {3., 4., 5., 1., 2.}, + {1., 2., 3., 4., 5.}, + {4., 5., 1., 2., 3.}, + {2., 3., 4., 5., 1.}}, + + {{3., 4., 5., 1., 2.}, + {5., 1., 2., 3., 4.}, + {2., 3., 4., 5., 1.}, + {4., 5., 1., 2., 3.}, + {1., 2., 3., 4., 5.}, + {3., 4., 5., 1., 2.}, + {5., 1., 2., 3., 4.}}}); // Loosen tolerance a bit due to FFTs EXPECT_TRUE(IsField3DEqualField3D(shifted.fromFieldAligned(input), expected, "RGN_ALL", @@ -159,26 +163,6 @@ TEST_F(ShiftedMetricTest, CalcYUpDown) { 0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz)); output_info.enable(); - Field3D input{mesh}; - - fillField(input, {{{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}, - - {{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}, - - {{1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}, - {1., 2., 3., 4., 5., 6., 7.}}}); - // Actual interesting bit here! ShiftedMetric shifted{*mesh, zShift}; shifted.calcYUpDown(input); @@ -188,83 +172,107 @@ TEST_F(ShiftedMetricTest, CalcYUpDown) { Field3D expected_up_1{mesh}; // Note: here zeroes are for values we don't expect to read - fillField(expected_up_1, {{{0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {2., 3., 4., 5., 6., 7., 1.}, - {0., 0., 0., 0., 0., 0., 0.}}, - - {{0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {2., 3., 4., 5., 6., 7., 1.}, - {0., 0., 0., 0., 0., 0., 0.}}, - - {{0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {2., 3., 4., 5., 6., 7., 1.}, - {0., 0., 0., 0., 0., 0., 0.}}}); + fillField(expected_up_1, {{{0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {2., 3., 4., 5., 1.}, + {2., 3., 4., 5., 1.}, + {2., 3., 4., 5., 1.}, + {0., 0., 0., 0., 0.}}, + + {{0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {3., 4., 5., 1., 2.}, + {3., 4., 5., 1., 2.}, + {3., 4., 5., 1., 2.}, + {0., 0., 0., 0., 0.}}, + + {{0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {4., 5., 1., 2., 3.}, + {4., 5., 1., 2., 3.}, + {4., 5., 1., 2., 3.}, + {0., 0., 0., 0., 0.}}}); Field3D expected_up_2{mesh}; - fillField(expected_up_2, {{{0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {3., 4., 5., 6., 7., 1., 2.}}, - - {{0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {3., 4., 5., 6., 7., 1., 2.}}, - - {{0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {3., 4., 5., 6., 7., 1., 2.}}}); + fillField(expected_up_2, {{{0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {3., 4., 5., 1., 2.}, + {3., 4., 5., 1., 2.}, + {3., 4., 5., 1., 2.}}, + + {{0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {5., 1., 2., 3., 4.}, + {5., 1., 2., 3., 4.}, + {5., 1., 2., 3., 4.}}, + + {{0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {2., 3., 4., 5., 1.}, + {2., 3., 4., 5., 1.}, + {2., 3., 4., 5., 1.}}}); Field3D expected_down_1{mesh}; - fillField(expected_down_1, {{{0., 0., 0., 0., 0., 0., 0.}, - {7., 1., 2., 3., 4., 5., 6.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}}, - - {{0., 0., 0., 0., 0., 0., 0.}, - {7., 1., 2., 3., 4., 5., 6.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}}, - - {{0., 0., 0., 0., 0., 0., 0.}, - {7., 1., 2., 3., 4., 5., 6.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}}}); + fillField(expected_down_1, {{{0., 0., 0., 0., 0.}, + {5., 1., 2., 3., 4.}, + {5., 1., 2., 3., 4.}, + {5., 1., 2., 3., 4.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}}, + + {{0., 0., 0., 0., 0.}, + {4., 5., 1., 2., 3.}, + {4., 5., 1., 2., 3.}, + {4., 5., 1., 2., 3.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}}, + + {{0., 0., 0., 0., 0.}, + {3., 4., 5., 1., 2.}, + {3., 4., 5., 1., 2.}, + {3., 4., 5., 1., 2.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}}}); Field3D expected_down2{mesh}; - fillField(expected_down2, {{{6., 7., 1., 2., 3., 4., 5.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}}, - - {{6., 7., 1., 2., 3., 4., 5.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}}, - - {{6., 7., 1., 2., 3., 4., 5.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0., 0.}}}); + fillField(expected_down2, {{{4., 5., 1., 2., 3.}, + {4., 5., 1., 2., 3.}, + {4., 5., 1., 2., 3.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}}, + + {{2., 3., 4., 5., 1.}, + {2., 3., 4., 5., 1.}, + {2., 3., 4., 5., 1.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}}, + + {{5., 1., 2., 3., 4.}, + {5., 1., 2., 3., 4.}, + {5., 1., 2., 3., 4.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0.}}}); EXPECT_TRUE(IsField3DEqualField3D(input.ynext(1), expected_up_1, "RGN_YUP")); EXPECT_TRUE(IsField3DEqualField3D(input.ynext(2), expected_up_2, "RGN_YUP2")); From bc133eca6de3083f33f42ad4c6d838b39ad89008 Mon Sep 17 00:00:00 2001 From: David Dickinson Date: Tue, 15 Jan 2019 14:28:28 +0000 Subject: [PATCH 12/34] Ensure populateStencil fills in extra guard cells from parallel slices --- include/stencils.hxx | 44 ++++++++++++++++++++++++++++++++++---------- 1 file changed, 34 insertions(+), 10 deletions(-) diff --git a/include/stencils.hxx b/include/stencils.hxx index 517929b1dc..fe3b8b1160 100644 --- a/include/stencils.hxx +++ b/include/stencils.hxx @@ -46,49 +46,73 @@ void inline populateStencil(stencil &s, const FieldType& f, const typename Field switch(stagger) { case(STAGGER::None): - if (nGuard == 2) s.mm = f[i.template minus<2, direction>()]; + if (nGuard == 2) { + if (direction == DIRECTION::YOrthogonal) { + s.mm = f.ynext(-2)[i.template minus<2, direction>()]; + } else { + s.mm = f[i.template minus<2, direction>()]; + } + } if (direction == DIRECTION::YOrthogonal) { - s.m = f.ydown()[i.template minus<1, direction>()]; + s.m = f.ynext(-1)[i.template minus<1, direction>()]; } else { s.m = f[i.template minus<1, direction>()]; } s.c = f[i]; if (direction == DIRECTION::YOrthogonal) { - s.p = f.yup()[i.template plus<1, direction>()]; + s.p = f.ynext(1)[i.template plus<1, direction>()]; } else { s.p = f[i.template plus<1, direction>()]; } - if (nGuard == 2) s.pp = f[i.template plus<2, direction>()]; + if (nGuard == 2) { + if (direction == DIRECTION::YOrthogonal) { + s.pp = f.ynext(2)[i.template plus<2, direction>()]; + } else { + s.pp = f[i.template plus<2, direction>()]; + } + } break; case(STAGGER::C2L): - if (nGuard == 2) s.mm = f[i.template minus<2, direction>()]; + if (nGuard == 2) { + if (direction == DIRECTION::YOrthogonal) { + s.mm = f.ynext(-2)[i.template minus<2, direction>()]; + } else { + s.mm = f[i.template minus<2, direction>()]; + } + } if (direction == DIRECTION::YOrthogonal) { - s.m = f.ydown()[i.template minus<1, direction>()]; + s.m = f.ynext(-1)[i.template minus<1, direction>()]; } else { s.m = f[i.template minus<1, direction>()]; } s.c = f[i]; s.p = s.c; if (direction == DIRECTION::YOrthogonal) { - s.pp = f.yup()[i.template plus<1, direction>()]; + s.pp = f.ynext(1)[i.template plus<1, direction>()]; } else { s.pp = f[i.template plus<1, direction>()]; } break; case(STAGGER::L2C): if (direction == DIRECTION::YOrthogonal) { - s.mm = f.ydown()[i.template minus<1, direction>()]; + s.mm = f.ynext(-1)[i.template minus<1, direction>()]; } else { s.mm = f[i.template minus<1, direction>()]; } s.m = f[i]; s.c = s.m; if (direction == DIRECTION::YOrthogonal) { - s.p = f.yup()[i.template plus<1, direction>()]; + s.p = f.ynext(1)[i.template plus<1, direction>()]; } else { s.p = f[i.template plus<1, direction>()]; } - if (nGuard == 2) s.pp = f[i.template plus<2, direction>()]; + if (nGuard == 2) { + if (direction == DIRECTION::YOrthogonal) { + s.pp = f.ynext(2)[i.template plus<2, direction>()]; + } else { + s.pp = f[i.template plus<2, direction>()]; + } + } break; } return; From ebcf7dd607b2795f41332a0ab8f4b631a1fde9cc Mon Sep 17 00:00:00 2001 From: David Dickinson Date: Tue, 15 Jan 2019 14:36:51 +0000 Subject: [PATCH 13/34] Add `ynext` to `Field2D` --- include/field2d.hxx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/include/field2d.hxx b/include/field2d.hxx index adbf8b0800..eae5befe62 100644 --- a/include/field2d.hxx +++ b/include/field2d.hxx @@ -129,7 +129,10 @@ class Field2D : public Field, public FieldData { const Field2D& ydown() const { return *this; } - + + Field2D& ynext(int UNUSED(dir)) { return *this; } + const Field2D& ynext(int UNUSED(dir)) const { return *this; } + // Operators /*! From 6e6dcd6fc0357d8526900d2a0215e2afd86c43da Mon Sep 17 00:00:00 2001 From: David Dickinson Date: Tue, 15 Jan 2019 14:37:14 +0000 Subject: [PATCH 14/34] Move macros and add YOrthogonal to list of automatic directions to register --- include/bout/index_derivs.hxx | 108 +++++++++++++++++----------------- 1 file changed, 54 insertions(+), 54 deletions(-) diff --git a/include/bout/index_derivs.hxx b/include/bout/index_derivs.hxx index af453a5f83..88522b780b 100644 --- a/include/bout/index_derivs.hxx +++ b/include/bout/index_derivs.hxx @@ -171,6 +171,60 @@ public: #define DEFINE_FLUX_DERIV_STAGGERED(name, key, nGuards, type) \ DEFINE_FLUX_DERIV(name, key, nGuards, type) +/// Some helper defines for now that allow us to wrap up enums +/// and the specific methods. +#define WRAP_ENUM(family, value) enumWrapper + +#define REGISTER_DERIVATIVE(name) \ + namespace { \ + produceCombinations, \ + Set, \ + Set, TypeContainer>, \ + Set>> \ + reg(registerMethod{}); \ + } +#define REGISTER_STAGGERED_DERIVATIVE(name) \ + namespace { \ + produceCombinations, \ + Set, \ + Set, TypeContainer>, \ + Set>> \ + reg(registerMethod{}); \ + } + +#define REGISTER_STANDARD_DERIVATIVE(name, key, nGuards, type) \ + DEFINE_STANDARD_DERIV_CORE(name, key, nGuards, type) \ + REGISTER_DERIVATIVE(name) \ + BoutReal name::operator()(const stencil& f) const + +#define REGISTER_UPWIND_DERIVATIVE(name, key, nGuards, type) \ + DEFINE_UPWIND_DERIV_CORE(name, key, nGuards, type) \ + REGISTER_DERIVATIVE(name) \ + BoutReal name::operator()(BoutReal vc, const stencil& f) const + +#define REGISTER_FLUX_DERIVATIVE(name, key, nGuards, type) \ + DEFINE_FLUX_DERIV_CORE(name, key, nGuards, type) \ + REGISTER_DERIVATIVE(name) \ + BoutReal name::operator()(const stencil& v, const stencil& f) const + +#define REGISTER_STANDARD_STAGGERED_DERIVATIVE(name, key, nGuards, type) \ + DEFINE_STANDARD_DERIV_CORE(name, key, nGuards, type) \ + REGISTER_STAGGERED_DERIVATIVE(name) \ + BoutReal name::operator()(const stencil& f) const + +#define REGISTER_UPWIND_STAGGERED_DERIVATIVE(name, key, nGuards, type) \ + /*Note staggered upwind looks like flux*/ \ + DEFINE_FLUX_DERIV_CORE(name, key, nGuards, type) \ + REGISTER_STAGGERED_DERIVATIVE(name) \ + BoutReal name::operator()(const stencil& v, const stencil& f) const + +#define REGISTER_FLUX_STAGGERED_DERIVATIVE(name, key, nGuards, type) \ + DEFINE_FLUX_DERIV_CORE(name, key, nGuards, type) \ + REGISTER_STAGGERED_DERIVATIVE(name) \ + BoutReal name::operator()(const stencil& v, const stencil& f) const + ////////////////////// FIRST DERIVATIVES ///////////////////// /// central, 2nd order @@ -559,60 +613,6 @@ struct registerMethod { } }; -/// Some helper defines for now that allow us to wrap up enums -/// and the specific methods. -#define WRAP_ENUM(family, value) enumWrapper - -#define REGISTER_DERIVATIVE(name) \ - namespace { \ - produceCombinations, \ - Set, \ - Set, TypeContainer>, \ - Set>> \ - reg(registerMethod{}); \ - } -#define REGISTER_STAGGERED_DERIVATIVE(name) \ - namespace { \ - produceCombinations, \ - Set, \ - Set, TypeContainer>, \ - Set>> \ - reg(registerMethod{}); \ - } - -#define REGISTER_STANDARD_DERIVATIVE(name, key, nGuards, type) \ - DEFINE_STANDARD_DERIV_CORE(name, key, nGuards, type) \ - REGISTER_DERIVATIVE(name) \ - BoutReal name::operator()(const stencil& f) const - -#define REGISTER_UPWIND_DERIVATIVE(name, key, nGuards, type) \ - DEFINE_UPWIND_DERIV_CORE(name, key, nGuards, type) \ - REGISTER_DERIVATIVE(name) \ - BoutReal name::operator()(BoutReal vc, const stencil& f) const - -#define REGISTER_FLUX_DERIVATIVE(name, key, nGuards, type) \ - DEFINE_FLUX_DERIV_CORE(name, key, nGuards, type) \ - REGISTER_DERIVATIVE(name) \ - BoutReal name::operator()(const stencil& v, const stencil& f) const - -#define REGISTER_STANDARD_STAGGERED_DERIVATIVE(name, key, nGuards, type) \ - DEFINE_STANDARD_DERIV_CORE(name, key, nGuards, type) \ - REGISTER_STAGGERED_DERIVATIVE(name) \ - BoutReal name::operator()(const stencil& f) const - -#define REGISTER_UPWIND_STAGGERED_DERIVATIVE(name, key, nGuards, type) \ - /*Note staggered upwind looks like flux*/ \ - DEFINE_FLUX_DERIV_CORE(name, key, nGuards, type) \ - REGISTER_STAGGERED_DERIVATIVE(name) \ - BoutReal name::operator()(const stencil& v, const stencil& f) const - -#define REGISTER_FLUX_STAGGERED_DERIVATIVE(name, key, nGuards, type) \ - DEFINE_FLUX_DERIV_CORE(name, key, nGuards, type) \ - REGISTER_STAGGERED_DERIVATIVE(name) \ - BoutReal name::operator()(const stencil& v, const stencil& f) const - ///////////////////////////////////////////////////////////////////////////////// /// Here's an example of registering a couple of DerivativeType methods /// at once for no staggering From c24dc664f8d2c2f917b2e8b82c28cdd659d0512a Mon Sep 17 00:00:00 2001 From: David Dickinson Date: Tue, 15 Jan 2019 14:51:07 +0000 Subject: [PATCH 15/34] Switch to define and register in one step for in-build derivative methods --- include/bout/index_derivs.hxx | 296 ++++++++++++++++------------------ 1 file changed, 136 insertions(+), 160 deletions(-) diff --git a/include/bout/index_derivs.hxx b/include/bout/index_derivs.hxx index 88522b780b..80e63dac79 100644 --- a/include/bout/index_derivs.hxx +++ b/include/bout/index_derivs.hxx @@ -123,6 +123,85 @@ public: const metaData meta = func.meta; }; +///////////////////////////////////////////////////////////////////////////////// +/// Following code is for dealing with registering a method/methods for all +/// template combinations, in conjunction with the template_combinations code. +///////////////////////////////////////////////////////////////////////////////// + +struct registerMethod { + template + void operator()(Direction, Stagger, FieldTypeContainer, Method) { + AUTO_TRACE(); + using namespace std::placeholders; + + // Now we want to get the actual field type out of the TypeContainer + // used to pass this around + using FieldType = typename FieldTypeContainer::type; + + Method method{}; + + // Note whilst this should be known at compile time using this directly in the + // template parameters below causes problems for old versions of gcc/libstdc++ + // (tested with 4.8.3) so we currently use a hacky workaround. Once we drop + // support for these versions the branching in the case statement below can be + // removed and we can use nGuard directly in the template statement. + const int nGuards = method.meta.nGuards; + + auto& derivativeRegister = DerivativeStore::getInstance(); + + switch (method.meta.derivType) { + case (DERIV::Standard): + case (DERIV::StandardSecond): + case (DERIV::StandardFourth): { + if (nGuards == 1) { + const auto theFunc = std::bind( + // Method to store in function + &Method::template standard, + // Arguments -- first is hidden this of type-bound, others are placeholders + // for input field, output field, region + method, _1, _2, _3); + derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); + } else { + const auto theFunc = std::bind( + // Method to store in function + &Method::template standard, + // Arguments -- first is hidden this of type-bound, others are placeholders + // for input field, output field, region + method, _1, _2, _3); + derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); + } + break; + } + case (DERIV::Upwind): + case (DERIV::Flux): { + if (nGuards == 1) { + const auto theFunc = std::bind( + // Method to store in function + &Method::template upwindOrFlux, + // Arguments -- first is hidden this of type-bound, others are placeholders + // for input field, output field, region + method, _1, _2, _3, _4); + derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); + } else { + const auto theFunc = std::bind( + // Method to store in function + &Method::template upwindOrFlux, + // Arguments -- first is hidden this of type-bound, others are placeholders + // for input field, output field, region + method, _1, _2, _3, _4); + derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); + } + break; + } + default: + throw BoutException("Unhandled derivative method in registerMethod."); + }; + } +}; + #define DEFINE_STANDARD_DERIV_CORE(name, key, nGuards, type) \ struct name { \ BoutReal operator()(const stencil& f) const; \ @@ -175,23 +254,23 @@ public: /// and the specific methods. #define WRAP_ENUM(family, value) enumWrapper -#define REGISTER_DERIVATIVE(name) \ - namespace { \ - produceCombinations, \ - Set, \ - Set, TypeContainer>, \ - Set>> \ - reg(registerMethod{}); \ + Set, \ + Set, TypeContainer>, \ + Set>> \ + reg##name(registerMethod{}); \ } -#define REGISTER_STAGGERED_DERIVATIVE(name) \ - namespace { \ - produceCombinations, \ - Set, \ - Set, TypeContainer>, \ - Set>> \ - reg(registerMethod{}); \ + Set, \ + Set, TypeContainer>, \ + Set>> \ + reg##name(registerMethod{}); \ } #define REGISTER_STANDARD_DERIVATIVE(name, key, nGuards, type) \ @@ -214,29 +293,39 @@ public: REGISTER_STAGGERED_DERIVATIVE(name) \ BoutReal name::operator()(const stencil& f) const +#define REGISTER_STANDARD_DERIVATIVE_STAGGERED(name, key, nGuards, type) \ + REGISTER_STANDARD_STAGGERED_DERIVATIVE(name, key, nGuards, type) + #define REGISTER_UPWIND_STAGGERED_DERIVATIVE(name, key, nGuards, type) \ /*Note staggered upwind looks like flux*/ \ DEFINE_FLUX_DERIV_CORE(name, key, nGuards, type) \ REGISTER_STAGGERED_DERIVATIVE(name) \ BoutReal name::operator()(const stencil& v, const stencil& f) const +#define REGISTER_UPWIND_DERIVATIVE_STAGGERED(name, key, nGuards, type) \ + REGISTER_UPWIND_STAGGERED_DERIVATIVE(name, key, nGuards, type) + #define REGISTER_FLUX_STAGGERED_DERIVATIVE(name, key, nGuards, type) \ DEFINE_FLUX_DERIV_CORE(name, key, nGuards, type) \ REGISTER_STAGGERED_DERIVATIVE(name) \ BoutReal name::operator()(const stencil& v, const stencil& f) const +#define REGISTER_FLUX_DERIVATIVE_STAGGERED(name, key, nGuards, type) \ + REGISTER_FLUX_STAGGERED_DERIVATIVE(name, key, nGuards, type) ////////////////////// FIRST DERIVATIVES ///////////////////// /// central, 2nd order -DEFINE_STANDARD_DERIV(DDX_C2, "C2", 1, DERIV::Standard) { return 0.5 * (f.p - f.m); }; +REGISTER_STANDARD_DERIVATIVE(DDX_C2, "C2", 1, DERIV::Standard) { + return 0.5 * (f.p - f.m); +}; /// central, 4th order -DEFINE_STANDARD_DERIV(DDX_C4, "C4", 2, DERIV::Standard) { +REGISTER_STANDARD_DERIVATIVE(DDX_C4, "C4", 2, DERIV::Standard) { return (8. * f.p - 8. * f.m + f.mm - f.pp) / 12.; } /// Central WENO method, 2nd order (reverts to 1st order near shocks) -DEFINE_STANDARD_DERIV(DDX_CWENO2, "W2", 1, DERIV::Standard) { +REGISTER_STANDARD_DERIVATIVE(DDX_CWENO2, "W2", 1, DERIV::Standard) { BoutReal isl, isr, isc; // Smoothness indicators BoutReal al, ar, ac, sa; // Un-normalised weights BoutReal dl, dr, dc; // Derivatives using different stencils @@ -258,7 +347,7 @@ DEFINE_STANDARD_DERIV(DDX_CWENO2, "W2", 1, DERIV::Standard) { } // Smoothing 2nd order derivative -DEFINE_STANDARD_DERIV(DDX_S2, "S2", 2, DERIV::Standard) { +REGISTER_STANDARD_DERIVATIVE(DDX_S2, "S2", 2, DERIV::Standard) { // 4th-order differencing BoutReal result = (8. * f.p - 8. * f.m + f.mm - f.pp) / 12.; @@ -275,19 +364,19 @@ DEFINE_STANDARD_DERIV(DDX_S2, "S2", 2, DERIV::Standard) { ////////////////////////////// /// Second derivative: Central, 2nd order -DEFINE_STANDARD_DERIV(D2DX2_C2, "C2", 1, DERIV::StandardSecond) { +REGISTER_STANDARD_DERIVATIVE(D2DX2_C2, "C2", 1, DERIV::StandardSecond) { return f.p + f.m - 2. * f.c; } /// Second derivative: Central, 4th order -DEFINE_STANDARD_DERIV(D2DX2_C4, "C4", 2, DERIV::StandardSecond) { +REGISTER_STANDARD_DERIVATIVE(D2DX2_C4, "C4", 2, DERIV::StandardSecond) { return (-f.pp + 16. * f.p - 30. * f.c + 16. * f.m - f.mm) / 12.; } ////////////////////////////// //--- Fourth order derivatives ////////////////////////////// -DEFINE_STANDARD_DERIV(D4DX4_C2, "C2", 2, DERIV::StandardFourth) { +REGISTER_STANDARD_DERIVATIVE(D4DX4_C2, "C2", 2, DERIV::StandardFourth) { return (f.pp - 4. * f.p + 6. * f.c - 4. * f.m + f.mm); } @@ -304,15 +393,17 @@ std::tuple vUpDown(BoutReal v) { } /// Upwinding: Central, 2nd order -DEFINE_UPWIND_DERIV(VDDX_C2, "C2", 1, DERIV::Upwind) { return vc * 0.5 * (f.p - f.m); } +REGISTER_UPWIND_DERIVATIVE(VDDX_C2, "C2", 1, DERIV::Upwind) { + return vc * 0.5 * (f.p - f.m); +} /// Upwinding: Central, 4th order -DEFINE_UPWIND_DERIV(VDDX_C4, "C4", 2, DERIV::Upwind) { +REGISTER_UPWIND_DERIVATIVE(VDDX_C4, "C4", 2, DERIV::Upwind) { return vc * (8. * f.p - 8. * f.m + f.mm - f.pp) / 12.; } /// upwind, 1st order -DEFINE_UPWIND_DERIV(VDDX_U1, "U1", 1, DERIV::Upwind) { // No vec +REGISTER_UPWIND_DERIVATIVE(VDDX_U1, "U1", 1, DERIV::Upwind) { // No vec // Existing form doesn't vectorise due to branching return vc >= 0.0 ? vc * (f.c - f.m) : vc * (f.p - f.c); // Alternative form would but may involve more operations @@ -321,7 +412,7 @@ DEFINE_UPWIND_DERIV(VDDX_U1, "U1", 1, DERIV::Upwind) { // No vec } /// upwind, 2nd order -DEFINE_UPWIND_DERIV(VDDX_U2, "U2", 2, DERIV::Upwind) { // No vec +REGISTER_UPWIND_DERIVATIVE(VDDX_U2, "U2", 2, DERIV::Upwind) { // No vec // Existing form doesn't vectorise due to branching return vc >= 0.0 ? vc * (1.5 * f.c - 2.0 * f.m + 0.5 * f.mm) : vc * (-0.5 * f.pp + 2.0 * f.p - 1.5 * f.c); @@ -332,7 +423,7 @@ DEFINE_UPWIND_DERIV(VDDX_U2, "U2", 2, DERIV::Upwind) { // No vec } /// upwind, 3rd order -DEFINE_UPWIND_DERIV(VDDX_U3, "U3", 2, DERIV::Upwind) { // No vec +REGISTER_UPWIND_DERIVATIVE(VDDX_U3, "U3", 2, DERIV::Upwind) { // No vec // Existing form doesn't vectorise due to branching return vc >= 0.0 ? vc * (4. * f.p - 12. * f.m + 2. * f.mm + 6. * f.c) / 12. : vc * (-4. * f.m + 12. * f.p - 2. * f.pp - 6. * f.c) / 12.; @@ -344,7 +435,7 @@ DEFINE_UPWIND_DERIV(VDDX_U3, "U3", 2, DERIV::Upwind) { // No vec } /// 3rd-order WENO scheme -DEFINE_UPWIND_DERIV(VDDX_WENO3, "W3", 2, DERIV::Upwind) { // No vec +REGISTER_UPWIND_DERIVATIVE(VDDX_WENO3, "W3", 2, DERIV::Upwind) { // No vec BoutReal deriv, w, r; // Existing form doesn't vectorise due to branching @@ -373,7 +464,7 @@ DEFINE_UPWIND_DERIV(VDDX_WENO3, "W3", 2, DERIV::Upwind) { // No vec ///----------------------------------------------------------------- /// 3rd-order CWENO. Uses the upwinding code and split flux -DEFINE_STANDARD_DERIV(DDX_CWENO3, "W3", 2, DERIV::Standard) { +REGISTER_STANDARD_DERIVATIVE(DDX_CWENO3, "W3", 2, DERIV::Standard) { BoutReal a, ma = fabs(f.c); // Split flux a = fabs(f.m); @@ -416,7 +507,7 @@ DEFINE_STANDARD_DERIV(DDX_CWENO3, "W3", 2, DERIV::Standard) { /// //////////////////////////////////////////////////////////////////////////////// -DEFINE_FLUX_DERIV(FDDX_U1, "U1", 1, DERIV::Flux) { // No vec +REGISTER_FLUX_DERIVATIVE(FDDX_U1, "U1", 1, DERIV::Flux) { // No vec // Velocity at lower end BoutReal vs = 0.5 * (v.m + v.c); @@ -432,9 +523,11 @@ DEFINE_FLUX_DERIV(FDDX_U1, "U1", 1, DERIV::Flux) { // No vec return result - std::get<0>(vSplit) * f.c + std::get<1>(vSplit) * f.p; } -DEFINE_FLUX_DERIV(FDDX_C2, "C2", 2, DERIV::Flux) { return 0.5 * (v.p * f.p - v.m * f.m); } +REGISTER_FLUX_DERIVATIVE(FDDX_C2, "C2", 2, DERIV::Flux) { + return 0.5 * (v.p * f.p - v.m * f.m); +} -DEFINE_FLUX_DERIV(FDDX_C4, "C4", 2, DERIV::Flux) { +REGISTER_FLUX_DERIVATIVE(FDDX_C4, "C4", 2, DERIV::Flux) { return (8. * v.p * f.p - 8. * v.m * f.m + v.mm * f.mm - v.pp * f.pp) / 12.; } @@ -461,25 +554,25 @@ DEFINE_FLUX_DERIV(FDDX_C4, "C4", 2, DERIV::Flux) { //////////////////////////////////////////////////////////////////////////////// /// Standard methods -- first order //////////////////////////////////////////////////////////////////////////////// -DEFINE_STANDARD_DERIV_STAGGERED(DDX_C2_stag, "C2", 1, DERIV::Standard) { +REGISTER_STANDARD_DERIVATIVE_STAGGERED(DDX_C2_stag, "C2", 1, DERIV::Standard) { return f.p - f.m; } -DEFINE_STANDARD_DERIV_STAGGERED(DDX_C4_stag, "C4", 2, DERIV::Standard) { +REGISTER_STANDARD_DERIVATIVE_STAGGERED(DDX_C4_stag, "C4", 2, DERIV::Standard) { return (27. * (f.p - f.m) - (f.pp - f.mm)) / 24.; } //////////////////////////////////////////////////////////////////////////////// /// Standard methods -- second order //////////////////////////////////////////////////////////////////////////////// -DEFINE_STANDARD_DERIV_STAGGERED(D2DX2_C2_stag, "C2", 2, DERIV::StandardSecond) { +REGISTER_STANDARD_DERIVATIVE_STAGGERED(D2DX2_C2_stag, "C2", 2, DERIV::StandardSecond) { return (f.pp + f.mm - f.p - f.m) / 2.; } //////////////////////////////////////////////////////////////////////////////// /// Upwind methods //////////////////////////////////////////////////////////////////////////////// -DEFINE_UPWIND_DERIV_STAGGERED(VDDX_U1_stag, "U1", 1, DERIV::Upwind) { +REGISTER_UPWIND_DERIVATIVE_STAGGERED(VDDX_U1_stag, "U1", 1, DERIV::Upwind) { // Lower cell boundary BoutReal result = (v.m >= 0) ? v.m * f.m : v.m * f.c; @@ -492,7 +585,7 @@ DEFINE_UPWIND_DERIV_STAGGERED(VDDX_U1_stag, "U1", 1, DERIV::Upwind) { return result; } -DEFINE_UPWIND_DERIV_STAGGERED(VDDX_U2_stag, "U2", 2, DERIV::Upwind) { +REGISTER_UPWIND_DERIVATIVE_STAGGERED(VDDX_U2_stag, "U2", 2, DERIV::Upwind) { // Calculate d(v*f)/dx = (v*f)[i+1/2] - (v*f)[i-1/2] // Upper cell boundary @@ -508,13 +601,13 @@ DEFINE_UPWIND_DERIV_STAGGERED(VDDX_U2_stag, "U2", 2, DERIV::Upwind) { return result; } -DEFINE_UPWIND_DERIV_STAGGERED(VDDX_C2_stag, "C2", 1, DERIV::Upwind) { +REGISTER_UPWIND_DERIVATIVE_STAGGERED(VDDX_C2_stag, "C2", 1, DERIV::Upwind) { // Result is needed at location of f: interpolate v to f's location and take an // unstaggered derivative of f return 0.5 * (v.p + v.m) * 0.5 * (f.p - f.m); } -DEFINE_UPWIND_DERIV_STAGGERED(VDDX_C4_stag, "C4", 2, DERIV::Upwind) { +REGISTER_UPWIND_DERIVATIVE_STAGGERED(VDDX_C4_stag, "C4", 2, DERIV::Upwind) { // Result is needed at location of f: interpolate v to f's location and take an // unstaggered derivative of f return (9. * (v.m + v.p) - v.mm - v.pp) / 16. * (8. * f.p - 8. * f.m + f.mm - f.pp) @@ -524,7 +617,7 @@ DEFINE_UPWIND_DERIV_STAGGERED(VDDX_C4_stag, "C4", 2, DERIV::Upwind) { //////////////////////////////////////////////////////////////////////////////// /// Flux methods //////////////////////////////////////////////////////////////////////////////// -DEFINE_FLUX_DERIV_STAGGERED(FDDX_U1_stag, "U1", 1, DERIV::Flux) { +REGISTER_FLUX_DERIVATIVE_STAGGERED(FDDX_U1_stag, "U1", 1, DERIV::Flux) { // Lower cell boundary BoutReal result = (v.m >= 0) ? v.m * f.m : v.m * f.c; @@ -534,127 +627,10 @@ DEFINE_FLUX_DERIV_STAGGERED(FDDX_U1_stag, "U1", 1, DERIV::Flux) { return -result; } -///////////////////////////////////////////////////////////////////////////////// -/// Following code is for dealing with registering a method/methods for all -/// template combinations, in conjunction with the template_combinations code. -///////////////////////////////////////////////////////////////////////////////// - -struct registerMethod { - template - void operator()(Direction, Stagger, FieldTypeContainer, Method) { - AUTO_TRACE(); - using namespace std::placeholders; - - // Now we want to get the actual field type out of the TypeContainer - // used to pass this around - using FieldType = typename FieldTypeContainer::type; - - Method method{}; - - // Note whilst this should be known at compile time using this directly in the - // template parameters below causes problems for old versions of gcc/libstdc++ - // (tested with 4.8.3) so we currently use a hacky workaround. Once we drop - // support for these versions the branching in the case statement below can be - // removed and we can use nGuard directly in the template statement. - const int nGuards = method.meta.nGuards; - - auto& derivativeRegister = DerivativeStore::getInstance(); - - switch (method.meta.derivType) { - case (DERIV::Standard): - case (DERIV::StandardSecond): - case (DERIV::StandardFourth): { - if (nGuards == 1) { - const auto theFunc = std::bind( - // Method to store in function - &Method::template standard, - // Arguments -- first is hidden this of type-bound, others are placeholders - // for input field, output field, region - method, _1, _2, _3); - derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); - } else { - const auto theFunc = std::bind( - // Method to store in function - &Method::template standard, - // Arguments -- first is hidden this of type-bound, others are placeholders - // for input field, output field, region - method, _1, _2, _3); - derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); - } - break; - } - case (DERIV::Upwind): - case (DERIV::Flux): { - if (nGuards == 1) { - const auto theFunc = std::bind( - // Method to store in function - &Method::template upwindOrFlux, - // Arguments -- first is hidden this of type-bound, others are placeholders - // for input field, output field, region - method, _1, _2, _3, _4); - derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); - } else { - const auto theFunc = std::bind( - // Method to store in function - &Method::template upwindOrFlux, - // Arguments -- first is hidden this of type-bound, others are placeholders - // for input field, output field, region - method, _1, _2, _3, _4); - derivativeRegister.registerDerivative(theFunc, Direction{}, Stagger{}, method); - } - break; - } - default: - throw BoutException("Unhandled derivative method in registerMethod."); - }; - } -}; - -///////////////////////////////////////////////////////////////////////////////// -/// Here's an example of registering a couple of DerivativeType methods -/// at once for no staggering -///////////////////////////////////////////////////////////////////////////////// - -produceCombinations, - Set, - Set, TypeContainer>, - Set< - // Standard - DerivativeType, DerivativeType, - DerivativeType, DerivativeType, - DerivativeType, - // Standard 2nd order - DerivativeType, DerivativeType, - // Standard 4th order - DerivativeType, - // Upwind - DerivativeType, DerivativeType, - DerivativeType, DerivativeType, - DerivativeType, DerivativeType, - // Flux - DerivativeType, DerivativeType, - DerivativeType>> - registerDerivatives(registerMethod{}); - -produceCombinations, - Set, - Set, TypeContainer>, - Set< - // Standard - DerivativeType, DerivativeType, - // Standard 2nd order - DerivativeType, - // Upwind - DerivativeType, DerivativeType, - DerivativeType, DerivativeType, - // Flux - DerivativeType>> - registerStaggeredDerivatives(registerMethod{}); +///////////////////////////////////////////////////////////////////////////////////// +/// Here's an example of defining and registering a custom method that doesn't fit +/// into the standard stencil based approach. +// ///////////////////////////////////////////////////////////////////////////////// class FFTDerivativeType { public: From 7c34454f097ab6573523846ada89d03e26c94198 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 15 Jan 2019 14:16:11 +0000 Subject: [PATCH 16/34] Silence output from fft_init in ShiftedMetric tests --- tests/unit/mesh/parallel/test_shiftedmetric.cxx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/unit/mesh/parallel/test_shiftedmetric.cxx b/tests/unit/mesh/parallel/test_shiftedmetric.cxx index 8a806506c3..38f368ce6a 100644 --- a/tests/unit/mesh/parallel/test_shiftedmetric.cxx +++ b/tests/unit/mesh/parallel/test_shiftedmetric.cxx @@ -1,6 +1,7 @@ #include "gtest/gtest.h" #include "bout/paralleltransform.hxx" +#include "fft.hxx" #include "test_extras.hxx" extern Mesh* mesh; @@ -18,6 +19,9 @@ class ShiftedMetricTest : public ::testing::Test { mesh->createDefaultRegions(); output_info.enable(); + // Make sure fft functions are quiet by setting fft_measure to false + bout::fft::fft_init(false); + zShift = Field2D{mesh}; fillField(zShift, {{1., 2., 3., 4., 5., 6., 7.}, From 5e6672fd87cd9b8012c25eb37b46b1a37a5e7408 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 15 Jan 2019 14:39:38 +0000 Subject: [PATCH 17/34] Use consistent tolerance for FFT-based functions across unit tests --- tests/unit/invert/test_fft.cxx | 17 +++++++--------- .../unit/mesh/parallel/test_shiftedmetric.cxx | 20 +++++++++++-------- tests/unit/test_extras.hxx | 4 +++- 3 files changed, 22 insertions(+), 19 deletions(-) diff --git a/tests/unit/invert/test_fft.cxx b/tests/unit/invert/test_fft.cxx index 37ed3cb9b8..735fcedb93 100644 --- a/tests/unit/invert/test_fft.cxx +++ b/tests/unit/invert/test_fft.cxx @@ -39,9 +39,6 @@ class FFTTest : public ::testing::TestWithParam { Array real_signal; Array fft_signal; - - // FFTs have a slightly looser tolerance than other functions - static constexpr BoutReal fft_tolerance{1.e-12}; }; // Test the FFT functions with both even- and odd-length real signals @@ -57,8 +54,8 @@ TEST_P(FFTTest, rfft) { EXPECT_EQ(output.size(), nmodes); for (int i = 0; i < nmodes; ++i) { - EXPECT_NEAR(real(output[i]), real(fft_signal[i]), fft_tolerance); - EXPECT_NEAR(imag(output[i]), imag(fft_signal[i]), fft_tolerance); + EXPECT_NEAR(real(output[i]), real(fft_signal[i]), FFTTolerance); + EXPECT_NEAR(imag(output[i]), imag(fft_signal[i]), FFTTolerance); } } @@ -72,7 +69,7 @@ TEST_P(FFTTest, irfft) { EXPECT_EQ(output.size(), size); for (int i = 0; i < size; ++i) { - EXPECT_NEAR(output[i], real_signal[i], fft_tolerance); + EXPECT_NEAR(output[i], real_signal[i], FFTTolerance); } } @@ -84,8 +81,8 @@ TEST_P(FFTTest, rfftWithArray) { EXPECT_EQ(output.size(), nmodes); for (int i = 0; i < nmodes; ++i) { - EXPECT_NEAR(real(output[i]), real(fft_signal[i]), fft_tolerance); - EXPECT_NEAR(imag(output[i]), imag(fft_signal[i]), fft_tolerance); + EXPECT_NEAR(real(output[i]), real(fft_signal[i]), FFTTolerance); + EXPECT_NEAR(imag(output[i]), imag(fft_signal[i]), FFTTolerance); } } @@ -97,7 +94,7 @@ TEST_P(FFTTest, irfftWithArray) { EXPECT_EQ(output.size(), size); for (int i = 0; i < size; ++i) { - EXPECT_NEAR(output[i], real_signal[i], fft_tolerance); + EXPECT_NEAR(output[i], real_signal[i], FFTTolerance); } } @@ -109,6 +106,6 @@ TEST_P(FFTTest, RoundTrip) { EXPECT_EQ(output.size(), real_signal.size()); for (int i = 0; i < size; ++i) { - EXPECT_NEAR(output[i], real_signal[i], fft_tolerance); + EXPECT_NEAR(output[i], real_signal[i], FFTTolerance); } } diff --git a/tests/unit/mesh/parallel/test_shiftedmetric.cxx b/tests/unit/mesh/parallel/test_shiftedmetric.cxx index 38f368ce6a..cd224592e2 100644 --- a/tests/unit/mesh/parallel/test_shiftedmetric.cxx +++ b/tests/unit/mesh/parallel/test_shiftedmetric.cxx @@ -1,8 +1,8 @@ #include "gtest/gtest.h" -#include "bout/paralleltransform.hxx" #include "fft.hxx" #include "test_extras.hxx" +#include "bout/paralleltransform.hxx" extern Mesh* mesh; @@ -105,8 +105,8 @@ TEST_F(ShiftedMetricTest, ToFieldAligned) { {4., 5., 1., 2., 3.}, {2., 3., 4., 5., 1.}}}); - EXPECT_TRUE( - IsField3DEqualField3D(shifted.toFieldAligned(input), expected, "RGN_ALL", 1.e-12)); + EXPECT_TRUE(IsField3DEqualField3D(shifted.toFieldAligned(input), expected, "RGN_ALL", + FFTTolerance)); } TEST_F(ShiftedMetricTest, FromFieldAligned) { @@ -140,7 +140,7 @@ TEST_F(ShiftedMetricTest, FromFieldAligned) { // Loosen tolerance a bit due to FFTs EXPECT_TRUE(IsField3DEqualField3D(shifted.fromFieldAligned(input), expected, "RGN_ALL", - 1.e-12)); + FFTTolerance)); } TEST_F(ShiftedMetricTest, CalcYUpDown) { @@ -278,8 +278,12 @@ TEST_F(ShiftedMetricTest, CalcYUpDown) { {0., 0., 0., 0., 0.}, {0., 0., 0., 0., 0.}}}); - EXPECT_TRUE(IsField3DEqualField3D(input.ynext(1), expected_up_1, "RGN_YUP")); - EXPECT_TRUE(IsField3DEqualField3D(input.ynext(2), expected_up_2, "RGN_YUP2")); - EXPECT_TRUE(IsField3DEqualField3D(input.ynext(-1), expected_down_1, "RGN_YDOWN")); - EXPECT_TRUE(IsField3DEqualField3D(input.ynext(-2), expected_down2, "RGN_YDOWN2")); + EXPECT_TRUE( + IsField3DEqualField3D(input.ynext(1), expected_up_1, "RGN_YUP", FFTTolerance)); + EXPECT_TRUE( + IsField3DEqualField3D(input.ynext(2), expected_up_2, "RGN_YUP2", FFTTolerance)); + EXPECT_TRUE( + IsField3DEqualField3D(input.ynext(-1), expected_down_1, "RGN_YDOWN", FFTTolerance)); + EXPECT_TRUE( + IsField3DEqualField3D(input.ynext(-2), expected_down2, "RGN_YDOWN2", FFTTolerance)); } diff --git a/tests/unit/test_extras.hxx b/tests/unit/test_extras.hxx index a4593c8657..6b0c711a8d 100644 --- a/tests/unit/test_extras.hxx +++ b/tests/unit/test_extras.hxx @@ -12,7 +12,9 @@ #include "field3d.hxx" #include "unused.hxx" -const BoutReal BoutRealTolerance = 1e-15; +static constexpr BoutReal BoutRealTolerance{1e-15}; +// FFTs have a slightly looser tolerance than other functions +static constexpr BoutReal FFTTolerance{1.e-12}; /// Does \p str contain \p substring? ::testing::AssertionResult IsSubString(const std::string &str, From 1cbf1338879ad03e898305fa35a8ebc0f008b52e Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 24 Jan 2019 09:55:36 +0000 Subject: [PATCH 18/34] Generalise Zoidberg to multiple parallel slices --- tools/pylib/zoidberg/zoidberg.py | 209 +++++++++++++++---------------- 1 file changed, 101 insertions(+), 108 deletions(-) diff --git a/tools/pylib/zoidberg/zoidberg.py b/tools/pylib/zoidberg/zoidberg.py index a7333d6762..e287132b5b 100644 --- a/tools/pylib/zoidberg/zoidberg.py +++ b/tools/pylib/zoidberg/zoidberg.py @@ -2,6 +2,11 @@ import numpy as np from boututils import datafile as bdata +from collections import namedtuple +from itertools import chain + +from . import fieldtracer +from .progress import update_progress # PyEVTK might be called pyevtk or evtk, depending on where it was # installed from @@ -14,13 +19,24 @@ except ImportError: have_evtk = False -# from . import grid -# from . import field -from . import fieldtracer -from .progress import update_progress + +def parallel_slice_field_name(field, offset): + """Form a unique, backwards-compatible name for field at a given offset + + Parameters + ---------- + field : str + Name of the field to convert + offset : int + Parallel slice offset + + """ + prefix = 'forward' if offset > 0 else 'backward' + suffix = "_{}".format(abs(offset)) if abs(offset) > 1 else "" + return "{}_{}{}".format(prefix, field, suffix) -def make_maps(grid, magnetic_field, quiet=False, **kwargs): +def make_maps(grid, magnetic_field, nslice=1, quiet=False, **kwargs): """Make the forward and backward FCI maps Parameters @@ -29,6 +45,8 @@ def make_maps(grid, magnetic_field, quiet=False, **kwargs): Grid generated by Zoidberg magnetic_field : :py:obj:`zoidberg.field.MagneticField` Zoidberg magnetic field object + nslice : int + Number of parallel slices in each direction quiet : bool Don't display progress bar kwargs @@ -51,110 +69,85 @@ def make_maps(grid, magnetic_field, quiet=False, **kwargs): shape = (nx, ny, nz) # Coordinates of each grid point - R = np.zeros( shape ) - Z = np.zeros( shape ) - - # Arrays to store X index at end of field-line - # starting from (x,y,z) and going forward in toroidal angle (y) - forward_xt_prime = np.zeros( shape ) - forward_zt_prime = np.zeros( shape ) - - forward_R = np.zeros( shape ) - forward_Z = np.zeros( shape ) - - # Same but going backwards in toroidal angle - backward_xt_prime = np.zeros( shape ) - backward_zt_prime = np.zeros( shape ) - - backward_R = np.zeros( shape ) - backward_Z = np.zeros( shape ) - - field_tracer = fieldtracer.FieldTracer(magnetic_field) + R = np.zeros(shape) + Z = np.zeros(shape) - try: - rtol = kwargs["rtol"] - except KeyError: - rtol = None - - # TODO: if axisymmetric, don't loop, do one slice and copy for j in range(ny): - if (not quiet) and (ny > 1): - update_progress(float(j)/float(ny-1), **kwargs) - - # Get this poloidal grid - pol, ycoord = grid.getPoloidalGrid(j) - - # Store coordinates - R[:,j,:] = pol.R - Z[:,j,:] = pol.Z - - # Get the next (forward) poloidal grid - pol_forward, y_forward = grid.getPoloidalGrid(j+1) - - # We only want the end point, as [0,...] is the initial position - coord = field_tracer.follow_field_lines(pol.R, pol.Z, [ycoord, y_forward], rtol=rtol)[1,...] - - # Store the coordinates in real space - forward_R[:,j,:] = coord[:,:,0] - forward_Z[:,j,:] = coord[:,:,1] - - # Get the indices into the forward poloidal grid - if pol_forward is None: - # No forward grid, so hit a boundary - xind = -1 - zind = -1 - else: - # Find the indices for these new locations on the forward poloidal grid - xcoord = coord[:,:,0] - zcoord = coord[:,:,1] - xind, zind = pol_forward.findIndex(xcoord, zcoord) - - # Check boundary defined by the field - outside = magnetic_field.boundary.outside(xcoord, y_forward, zcoord) - xind[outside] = -1 - zind[outside] = -1 - - forward_xt_prime[:,j,:] = xind - forward_zt_prime[:,j,:] = zind - - # Go backwards one poloidal grid - pol_back, y_back = grid.getPoloidalGrid(j-1) - - # We only want the end point, as [0,...] is the initial position - coord = field_tracer.follow_field_lines(pol.R, pol.Z, [ycoord, y_back], rtol=rtol)[1,...] - - # Store the coordinates in real space - backward_R[:,j,:] = coord[:,:,0] - backward_Z[:,j,:] = coord[:,:,1] - - if pol_back is None: - # Hit boundary - xind = -1 - zind = -1 - else: - # Find the indices for these new locations on the backward poloidal grid - xcoord = coord[:,:,0] - zcoord = coord[:,:,1] - xind, zind = pol_back.findIndex(xcoord, zcoord) - - # Check boundary defined by the field - outside = magnetic_field.boundary.outside(xcoord, y_back, zcoord) - xind[outside] = -1 - zind[outside] = -1 - - backward_xt_prime[:,j,:] = xind - backward_zt_prime[:,j,:] = zind + pol, _ = grid.getPoloidalGrid(j) + R[:, j, :] = pol.R + Z[:, j, :] = pol.Z + + field_tracer = fieldtracer.FieldTracer(magnetic_field) + + rtol = kwargs.get("rtol", None) + # The field line maps and coordinates, etc. maps = { - 'R' : R, 'Z':Z, - 'forward_R':forward_R, 'forward_Z':forward_Z, - 'forward_xt_prime' : forward_xt_prime, - 'forward_zt_prime' : forward_zt_prime, - 'backward_R':backward_R, 'backward_Z':backward_Z, - 'backward_xt_prime' : backward_xt_prime, - 'backward_zt_prime' : backward_zt_prime + 'R': R, + 'Z': Z, } + # A helper data structure that groups the various field line maps along with the offset + ParallelSlice = namedtuple('ParallelSlice', ['R', 'Z', 'xt_prime', 'zt_prime', 'offset']) + # A list of the above data structures for each offset we want + parallel_slices = [] + + # Loop over offsets {1, ... nslice, -1, ... -nslice} + for offset in chain(range(1, nslice + 1), range(-1, -(nslice + 1), -1)): + # Unique names of the field line maps for this offset + field_names = [parallel_slice_field_name(field, offset) + for field in ['R', 'Z', 'xt_prime', 'zt_prime']] + + # Initialise the field arrays -- puts them straight into the result dict + for field in field_names: + maps[field] = np.zeros(shape) + + # Get the field arrays we just made and wrap them up in our helper tuple + fields = map(lambda x: maps[x], field_names) + parallel_slices.append(ParallelSlice(*fields, offset)) + + # Total size of the progress bar + total_work = float((len(parallel_slices) - 1) * (ny-1)) + + # TODO: if axisymmetric, don't loop, do one slice and copy + # TODO: restart tracing for adjacent offsets + for slice_index, parallel_slice in enumerate(parallel_slices): + for j in range(ny): + if (not quiet) and (ny > 1): + update_progress(float(slice_index * j) / total_work, **kwargs) + + # Get this poloidal grid + pol, ycoord = grid.getPoloidalGrid(j) + + # Get the next poloidal grid + pol_slice, y_slice = grid.getPoloidalGrid(j + parallel_slice.offset) + + # We only want the end point, as [0,...] is the initial position + coord = field_tracer.follow_field_lines(pol.R, pol.Z, [ycoord, y_slice], rtol=rtol)[1, ...] + + # Store the coordinates in real space + parallel_slice.R[:, j, :] = coord[:, :, 0] + parallel_slice.Z[:, j, :] = coord[:, :, 1] + + # Get the indices into the slice poloidal grid + if pol_slice is None: + # No slice grid, so hit a boundary + xind = -1 + zind = -1 + else: + # Find the indices for these new locations on the slice poloidal grid + xcoord = coord[:, :, 0] + zcoord = coord[:, :, 1] + xind, zind = pol_slice.findIndex(xcoord, zcoord) + + # Check boundary defined by the field + outside = magnetic_field.boundary.outside(xcoord, y_slice, zcoord) + xind[outside] = -1 + zind[outside] = -1 + + parallel_slice.xt_prime[:, j, :] = xind + parallel_slice.zt_prime[:, j, :] = zind + return maps @@ -203,7 +196,7 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', Rmaj[:,yindex,:] = magnetic_field.Rfunc(pol_grid.R, pol_grid.Z, ypos) metric["gyy"] = 1./Rmaj**2 metric["g_yy"] = Rmaj**2 - + # Get magnetic field and pressure Bmag = np.zeros(grid.shape) pressure = np.zeros(grid.shape) @@ -220,7 +213,7 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', pol_grid, ypos = grid.getPoloidalGrid(yindex) attribute[:,yindex,:] = magnetic_field.attributes[name](pol_grid.R, pol_grid.Z, ypos) attributes[name] = attribute - + # Metric is now 3D if metric2d: # Remove the Z dimension from metric components @@ -250,7 +243,7 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', f.write("ixseps2",ixseps) # Metric tensor - + if new_names: for key, val in metric.items(): f.write(key, val) @@ -271,7 +264,7 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', if name in name_changes: name = name_changes[name] f.write(name, metric[key]) - + # Magnetic field f.write("B", Bmag) @@ -281,7 +274,7 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', # Attributes for name in attributes: f.write(name, attributes[name]) - + # Maps - write everything to file for key in maps: f.write(key, maps[key]) From 1e0a6148a2356f19f281aba3b64e55287caf5c95 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 24 Jan 2019 14:35:33 +0000 Subject: [PATCH 19/34] Remove unneeded yup/ydown phase caches from ShiftedMetric --- include/bout/paralleltransform.hxx | 3 --- 1 file changed, 3 deletions(-) diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx index 5ed14cd4f2..c7b6ebf7d8 100644 --- a/include/bout/paralleltransform.hxx +++ b/include/bout/paralleltransform.hxx @@ -151,9 +151,6 @@ private: /// where i goes from 0 to (n-1), with n the number of y guard cells std::vector parallel_slice_phases; - arr3Dvec yupPhs; ///< Cache of phase shifts for calculating yup fields - arr3Dvec ydownPhs; ///< Cache of phase shifts for calculating ydown fields - /*! * Shift a 2D field in Z. * Since 2D fields are constant in Z, this has no effect From 497838772abfe4373f9232e1a89ca7de44f18f14 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 24 Jan 2019 16:46:10 +0000 Subject: [PATCH 20/34] Generalise Zoidberg test to multiple parallel slices --- tools/pylib/zoidberg/test_zoidberg.py | 67 +++++++++++++-------------- 1 file changed, 33 insertions(+), 34 deletions(-) diff --git a/tools/pylib/zoidberg/test_zoidberg.py b/tools/pylib/zoidberg/test_zoidberg.py index d82c30f0e2..23117c8dad 100644 --- a/tools/pylib/zoidberg/test_zoidberg.py +++ b/tools/pylib/zoidberg/test_zoidberg.py @@ -1,69 +1,68 @@ - +from itertools import chain, product import numpy as np from . import zoidberg, grid, field + def test_make_maps_slab(): nx = 5 - ny = 6 + ny = 8 nz = 7 # Create a straight magnetic field in a slab straight_field = field.Slab(By=1.0, Bz=0.0, Bzprime=0.0) - + # Create a rectangular grid in (x,y,z) - rectangle = grid.rectangular_grid(nx,ny,nz) + rectangle = grid.rectangular_grid(nx, ny, nz) + + # Two parallel slices in each direction + nslice = 2 # Calculate forwards and backwards maps - maps = zoidberg.make_maps(rectangle, straight_field) - - # Check that maps has the required forward and backward index variables - for var in ['forward_xt_prime', 'forward_zt_prime', 'backward_xt_prime', 'backward_zt_prime']: - assert var in maps + maps = zoidberg.make_maps(rectangle, straight_field, nslice=nslice) - # Each map should have the same shape as the grid - assert maps['forward_xt_prime'].shape == (nx,ny,nz) - assert maps['backward_xt_prime'].shape == (nx,ny,nz) - assert maps['forward_zt_prime'].shape == (nx,ny,nz) - assert maps['backward_zt_prime'].shape == (nx,ny,nz) - - # Since this is a straight magnetic field in a simple rectangle, + # Since this is a straight magnetic field in a simple rectangle, # all the maps should be the same, and should be the identity - identity_map_x, identity_map_z = np.meshgrid(np.arange(nx), np.arange(nz), indexing='ij') - for y in range(ny-1): - assert np.allclose(maps['forward_xt_prime'][:,y,:], identity_map_x) - assert np.allclose(maps['forward_zt_prime'][:,y,:], identity_map_z) + # Check that maps has the required forward and backward index variables + offsets = chain(range(1, nslice + 1), range(-1, -(nslice + 1), -1)) + field_line_maps = ["xt_prime", "zt_prime"] + + for field_line_map, offset in product(field_line_maps, offsets): + var = zoidberg.parallel_slice_field_name(field_line_map, offset) + print("Current field: ", var) + assert var in maps + + # Each map should have the same shape as the grid + maps[var].shape == (nx, ny, nz) - for y in range(1,ny): - assert np.allclose(maps['backward_xt_prime'][:,y,:], identity_map_x) - assert np.allclose(maps['backward_zt_prime'][:,y,:], identity_map_z) + # The first/last abs(offset) points are not valid, so ignore those + interior_range = range(ny-abs(offset)) if offset > 0 else range(abs(offset), ny) + # Those invalid points should be set to -1 + end_slice = slice(-1, -(offset + 1), -1) if offset > 0 else slice(0, -offset) + identity_map = identity_map_x if "x" in var else identity_map_z - # The last forward map should hit a boundary - assert np.allclose(maps['forward_xt_prime'][:,-1,:], -1.0) - assert np.allclose(maps['forward_zt_prime'][:,-1,:], -1.0) + for y in interior_range: + assert np.allclose(maps[var][:, y, :], identity_map) - # First backward map hits boundary - assert np.allclose(maps['backward_xt_prime'][:,0,:], -1.0) - assert np.allclose(maps['backward_zt_prime'][:,0,:], -1.0) + # The end slice should hit a boundary + assert np.allclose(maps[var][:, end_slice, :], -1.0) def test_make_maps_straight_stellarator(): nx = 5 ny = 6 nz = 7 - + # Create magnetic field magnetic_field = field.StraightStellarator(radius = np.sqrt(2.0)) - + # Create a rectangular grid in (x,y,z) rectangle = grid.rectangular_grid(nx,ny,nz, Lx = 1.0, Lz = 1.0, Ly = 10.0, yperiodic = True) - + # Here both the field and and grid are centred at (x,z) = (0,0) # and the rectangular grid here fits entirely within the coils maps = zoidberg.make_maps(rectangle, magnetic_field) - - From af2ca9b1e78c83f35fd5bc79fa9db16e56c4ad4c Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 25 Jan 2019 14:04:53 +0000 Subject: [PATCH 21/34] Add proper scaling test for FCI with C2 and C4 --- tests/MMS/.gitignore | 3 + tests/MMS/spatial/fci/data/BOUT.inp | 23 ++++ tests/MMS/spatial/fci/fci_mms.cxx | 28 +++++ tests/MMS/spatial/fci/makefile | 6 + tests/MMS/spatial/fci/mms.py | 35 ++++++ tests/MMS/spatial/fci/runtest | 176 ++++++++++++++++++++++++++++ tools/pylib/zoidberg/zoidberg.py | 54 +++++---- 7 files changed, 300 insertions(+), 25 deletions(-) create mode 100644 tests/MMS/spatial/fci/data/BOUT.inp create mode 100644 tests/MMS/spatial/fci/fci_mms.cxx create mode 100644 tests/MMS/spatial/fci/makefile create mode 100755 tests/MMS/spatial/fci/mms.py create mode 100755 tests/MMS/spatial/fci/runtest diff --git a/tests/MMS/.gitignore b/tests/MMS/.gitignore index 2946406b9b..c8f94147b3 100644 --- a/tests/MMS/.gitignore +++ b/tests/MMS/.gitignore @@ -30,6 +30,9 @@ BOUT.settings /spatial/d2dx2/test_d2dx2 /spatial/d2dz2/test_d2dz2 /spatial/diffusion/diffusion +/spatial/fci/fci_mms +/spatial/fci/fci.grid.nc +/spatial/fci/fci_mms.pkl /time/time /tokamak/tokamak /tokamak/tokamak.pkl diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp new file mode 100644 index 0000000000..6ad5d66ca2 --- /dev/null +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -0,0 +1,23 @@ +grid = fci.grid.nc + +input = sin(y - 2*z) + sin(y - z) + +solution = 6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z) + +MXG = 1 +NXPE = 1 + +[mesh] +paralleltransform = fci +symmetricglobalx = true + +[mesh:ddy] +first = C2 +second = C2 + +[fci] +y_periodic = true +z_periodic = true + +[interpolation] +type = lagrange4pt diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx new file mode 100644 index 0000000000..32b668f575 --- /dev/null +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -0,0 +1,28 @@ +#include "bout.hxx" +#include "derivs.hxx" +#include "field_factory.hxx" + +int main(int argc, char** argv) { + BoutInitialise(argc, argv); + + Field3D input{FieldFactory::get()->create3D("input", Options::getRoot(), mesh)}; + Field3D solution{FieldFactory::get()->create3D("solution", Options::getRoot(), mesh)}; + + // Communicate to calculate parallel transform + mesh->communicate(input); + + Field3D result{DDY(input)}; + Field3D error{result - solution}; + BoutReal l_2{sqrt(mean(SQ(error)))}; + BoutReal l_inf{max(abs(error), true)}; + + SAVE_ONCE6(input, solution, result, error, l_2, l_inf); + + for (int slice = 1; slice < mesh->ystart; ++slice) { + SAVE_ONCE2(input.ynext(-slice), input.ynext(slice)); + } + + dump.write(); + + BoutFinalise(); +} diff --git a/tests/MMS/spatial/fci/makefile b/tests/MMS/spatial/fci/makefile new file mode 100644 index 0000000000..88ba6c77e7 --- /dev/null +++ b/tests/MMS/spatial/fci/makefile @@ -0,0 +1,6 @@ + +BOUT_TOP = ../../../.. + +SOURCEC = fci_mms.cxx + +include $(BOUT_TOP)/make.config diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py new file mode 100755 index 0000000000..1154cfc9cc --- /dev/null +++ b/tests/MMS/spatial/fci/mms.py @@ -0,0 +1,35 @@ +#!/usr/bin/env python3 +# +# Generate manufactured solution and sources for FCI test +# + +from __future__ import division +from __future__ import print_function + +from boutdata.mms import * + +from sympy import sin, cos, sqrt + +from math import pi + +f = sin(y - z) + sin(y - 2*z) + +Lx = 0.1 +Ly = 10. +Lz = 1. + +Bt = 1.0 +Bp = 0.05 +Bpprime = 0.1 + +Bpx = Bp + (x-0.5)*Lx * Bpprime # Note: x in range [0,1] +B = sqrt(Bpx**2 + Bt**2) + +def FCI_ddy(f): + return ( Bt * diff(f, y)*2.*pi/Ly + Bpx * diff(f, z)*2.*pi/Lz ) + +############################################ +# Equations solved + +print("input = " + exprToStr(f)) +print("solution = " + exprToStr(FCI_ddy(f))) diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest new file mode 100755 index 0000000000..36eb96c7d7 --- /dev/null +++ b/tests/MMS/spatial/fci/runtest @@ -0,0 +1,176 @@ +#!/usr/bin/env python3 +# +# Python script to run and analyse MMS test +# +from __future__ import division +from __future__ import print_function + +from boututils.run_wrapper import shell_safe, launch_safe, getmpirun +from boutdata.collect import collect + +from numpy import array, log, polyfit, linspace, arange + +import pickle + +from sys import stdout + +import zoidberg as zb + +nx = 3 # Not changed for these tests + +# Resolution in y and z +nlist = [8, 16, 32, 64, 128] + +# Number of parallel slices (in each direction) +nslices = [1, 2] + +directory = "data" + +nproc = 2 +mthread = 2 + +MPIRUN = getmpirun() + +success = True + +error_2 = {} +error_inf = {} +method_orders = {} + +# Run with periodic Y? +yperiodic = True + +failures = [] + +print("Making fci MMS test") +shell_safe("make > make.log") + +for nslice in nslices: + error_2[nslice] = [] + error_inf[nslice] = [] + + # Which central difference scheme to use and its expected order + order = nslice * 2 + method_orders[nslice] = { + "name": "C{}".format(order), + "order": order + } + + for n in nlist: + # Define the magnetic field using new poloidal gridding method + # Note that the Bz and Bzprime parameters here must be the same as in mms.py + field = zb.field.Slab(Bz=0.05, Bzprime=0.1) + # Create rectangular poloidal grids + poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid(nx, n, 1., 1.) + # Set the ylength and y locations + ylength = 10. + + if yperiodic: + ycoords = linspace(0.0, ylength, n, endpoint=False) + else: + # Doesn't include the end points + ycoords = (arange(n) + 0.5)*ylength/float(n) + + # Create the grid + grid = zb.grid.Grid(poloidal_grid, ycoords, ylength, yperiodic=yperiodic) + # Make and write maps + maps = zb.make_maps(grid, field, nslice=nslice, quiet=True) + zb.write_maps(grid, field, maps, new_names=False, metric2d=True, quiet=True) + + args = (" MZ={} MYG={} fci:y_periodic={} mesh:ddy:first={}" + .format(n, nslice, yperiodic, method_orders[nslice]["name"])) + + # Command to run + cmd = "./fci_mms "+args + + print("Running command: "+cmd) + + # Launch using MPI + s, out = launch_safe(cmd, runcmd=MPIRUN, nproc=nproc, mthread=mthread, pipe=True) + + # Save output to log file + with open("run.log."+str(n), "w") as f: + f.write(out) + + if s: + print("Run failed!\nOutput was:\n") + print(out) + exit(s) + + # Collect data + l_2 = collect("l_2", tind=[1, 1], info=False, + path=directory, xguards=False, yguards=False) + l_inf = collect("l_inf", tind=[1, 1], info=False, + path=directory, xguards=False, yguards=False) + + error_2[nslice].append(l_2) + error_inf[nslice].append(l_inf) + + print("Errors : l-2 {:f} l-inf {:f}".format(l_2, l_inf)) + + dx = 1. / array(nlist) + + # Calculate convergence order + fit = polyfit(log(dx), log(error_2[nslice]), 1) + order = fit[0] + stdout.write("Convergence order = {:f} (fit)".format(order)) + + order = log(error_2[nslice][-2]/error_2[nslice][-1])/log(dx[-2]/dx[-1]) + stdout.write(", {:f} (small spacing)".format(order)) + + # Should be close to the expected order + if order > method_orders[nslice]["order"] * 0.95: + print("............ PASS\n") + else: + print("............ FAIL\n") + success = False + failures.append(method_orders[nslice]["name"]) + + +with open("fci_mms.pkl", "wb") as output: + pickle.dump(nlist, output) + for nslice in nslices: + pickle.dump(error_2[nslice], output) + pickle.dump(error_inf[nslice], output) + +# Do we want to show the plot as well as save it to file. +showPlot = True + +if False: + try: + # Plot using matplotlib if available + import matplotlib.pyplot as plt + + fig, ax = plt.subplots(1, 1) + + for nslice in nslices: + ax.plot(dx, error_2[nslice], '-', + label="{} $l_2$".format(method_orders[nslice]["name"])) + ax.plot(dx, error_inf[nslice], '--', + label="{} $l_\inf$".format(method_orders[nslice]["name"])) + ax.legend(loc="upper left") + ax.grid() + ax.set_yscale('log') + ax.set_xscale('log') + ax.set_title('error scaling') + ax.set_xlabel(r'Mesh spacing $\delta x$') + ax.set_ylabel("Error norm") + + plt.savefig("fci_mms.pdf") + + print("Plot saved to fci_mms.pdf") + + if showPlot: + plt.show() + plt.close() + except ImportError: + print("No matplotlib") + +if success: + print("All tests passed") + exit(0) +else: + print("Some tests failed:") + for failure in failures: + print("\t" + failure) + exit(1) diff --git a/tools/pylib/zoidberg/zoidberg.py b/tools/pylib/zoidberg/zoidberg.py index e287132b5b..3217b70c76 100644 --- a/tools/pylib/zoidberg/zoidberg.py +++ b/tools/pylib/zoidberg/zoidberg.py @@ -152,7 +152,8 @@ def make_maps(grid, magnetic_field, nslice=1, quiet=False, **kwargs): def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', - new_names=False, metric2d=True, format="NETCDF3_64BIT"): + new_names=False, metric2d=True, format="NETCDF3_64BIT", + quiet=False): """Write FCI maps to BOUT++ grid file Parameters @@ -169,8 +170,10 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', Write "g_yy" rather than "g_22" metric2d : bool, optional Output only 2D metrics - format : str + format : str, optional Specifies file format to use, passed to boutdata.DataFile + quiet : bool, optional + Don't warn about 2D metrics Returns ------- @@ -186,14 +189,14 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', # Check if the magnetic field is in cylindrical coordinates # If so, we need to change the gyy and g_yy metrics - pol_grid,ypos = grid.getPoloidalGrid(0) + pol_grid, ypos = grid.getPoloidalGrid(0) Rmaj = magnetic_field.Rfunc(pol_grid.R, pol_grid.Z, ypos) if Rmaj is not None: # In cylindrical coordinates Rmaj = np.zeros(grid.shape) for yindex in range(grid.numberOfPoloidalGrids()): - pol_grid,ypos = grid.getPoloidalGrid(yindex) - Rmaj[:,yindex,:] = magnetic_field.Rfunc(pol_grid.R, pol_grid.Z, ypos) + pol_grid, ypos = grid.getPoloidalGrid(yindex) + Rmaj[:, yindex, :] = magnetic_field.Rfunc(pol_grid.R, pol_grid.Z, ypos) metric["gyy"] = 1./Rmaj**2 metric["g_yy"] = Rmaj**2 @@ -201,9 +204,9 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', Bmag = np.zeros(grid.shape) pressure = np.zeros(grid.shape) for yindex in range(grid.numberOfPoloidalGrids()): - pol_grid,ypos = grid.getPoloidalGrid(yindex) - Bmag[:,yindex,:] = magnetic_field.Bmag(pol_grid.R, pol_grid.Z, ypos) - pressure[:,yindex,:] = magnetic_field.pressure(pol_grid.R, pol_grid.Z, ypos) + pol_grid, ypos = grid.getPoloidalGrid(yindex) + Bmag[:, yindex, :] = magnetic_field.Bmag(pol_grid.R, pol_grid.Z, ypos) + pressure[:, yindex, :] = magnetic_field.pressure(pol_grid.R, pol_grid.Z, ypos) # Get attributes from magnetic field (e.g. psi) attributes = {} @@ -211,23 +214,24 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', attribute = np.zeros(grid.shape) for yindex in range(grid.numberOfPoloidalGrids()): pol_grid, ypos = grid.getPoloidalGrid(yindex) - attribute[:,yindex,:] = magnetic_field.attributes[name](pol_grid.R, pol_grid.Z, ypos) + attribute[:, yindex, :] = magnetic_field.attributes[name](pol_grid.R, pol_grid.Z, ypos) attributes[name] = attribute # Metric is now 3D if metric2d: # Remove the Z dimension from metric components - print("WARNING: Outputting 2D metrics, discarding metric information.") + if not quiet: + print("WARNING: Outputting 2D metrics, discarding metric information.") for key in metric: try: - metric[key] = metric[key][:,:,0] - except: + metric[key] = metric[key][:, :, 0] + except TypeError: pass # Make dz a constant - metric["dz"] = metric["dz"][0,0] + metric["dz"] = metric["dz"][0, 0] # Add Rxy, Bxy - metric["Rxy"] = maps["R"][:,:,0] - metric["Bxy"] = Bmag[:,:,0] + metric["Rxy"] = maps["R"][:, :, 0] + metric["Bxy"] = Bmag[:, :, 0] with bdata.DataFile(gridfile, write=True, create=True, format=format) as f: ixseps = nx+1 @@ -239,8 +243,8 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', f.write("dy", metric["dy"]) f.write("dz", metric["dz"]) - f.write("ixseps1",ixseps) - f.write("ixseps2",ixseps) + f.write("ixseps1", ixseps) + f.write("ixseps2", ixseps) # Metric tensor @@ -251,14 +255,14 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', # Translate between output variable names and metric names # Map from new to old names. Anything not in this dict # is output unchanged - name_changes = {"g_yy":"g_22", - "gyy":"g22", - "gxx":"g11", - "gxz":"g13", - "gzz":"g33", - "g_xx":"g_11", - "g_xz":"g_13", - "g_zz":"g_33"} + name_changes = {"g_yy": "g_22", + "gyy": "g22", + "gxx": "g11", + "gxz": "g13", + "gzz": "g33", + "g_xx": "g_11", + "g_xz": "g_13", + "g_zz": "g_33"} for key in metric: name = key if name in name_changes: From 7c41f598d12ec87a23558a99cdcc4932829c3426 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 25 Jan 2019 14:14:23 +0000 Subject: [PATCH 22/34] Generalise FCI transform to multiple parallel slices --- src/mesh/parallel/fci.cxx | 72 +++++++++++++++++++-------------------- src/mesh/parallel/fci.hxx | 44 +++++++++++++++--------- 2 files changed, 63 insertions(+), 53 deletions(-) diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index eb62510406..89558fc628 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -42,10 +42,12 @@ #include "parallel_boundary_region.hxx" #include #include -#include // See this for codes +#include #include #include +#include + /** * Return the sign of val */ @@ -53,14 +55,18 @@ inline BoutReal sgn(BoutReal val) { return (BoutReal(0) < val) - (val < BoutReal // Calculate all the coefficients needed for the spline interpolation // dir MUST be either +1 or -1 -FCIMap::FCIMap(Mesh &mesh, int dir, bool zperiodic) - : dir(dir), boundary_mask(mesh), corner_boundary_mask(mesh), y_prime(&mesh) { +FCIMap::FCIMap(Mesh &mesh, int offset_, BoundaryRegionPar* boundary, bool zperiodic) + : offset(offset_), boundary_mask(mesh), corner_boundary_mask(mesh), y_prime(&mesh) { + + if (offset == 0) { + throw BoutException("FCIMap called with offset = 0; You probably didn't mean to do that"); + } interp = InterpolationFactory::getInstance()->create(&mesh); - interp->setYOffset(dir); + interp->setYOffset(offset); interp_corner = InterpolationFactory::getInstance()->create(&mesh); - interp_corner->setYOffset(dir); + interp_corner->setYOffset(offset); // Index arrays contain guard cells in order to get subscripts right // x-index of bottom-left grid point @@ -76,28 +82,20 @@ FCIMap::FCIMap(Mesh &mesh, int dir, bool zperiodic) mesh.get(R, "R", 0.0, false); mesh.get(Z, "Z", 0.0, false); - // Load the floating point indices from the grid file - // Future, higher order parallel derivatives could require maps to +/-2 slices - if (dir == +1) { - mesh.get(xt_prime, "forward_xt_prime", 0.0, false); - mesh.get(zt_prime, "forward_zt_prime", 0.0, false); - mesh.get(R_prime, "forward_R", 0.0, false); - mesh.get(Z_prime, "forward_Z", 0.0, false); - boundary = new BoundaryRegionPar("FCI_forward", BNDRY_PAR_FWD, dir, &mesh); - } else if (dir == -1) { - mesh.get(xt_prime, "backward_xt_prime", 0.0, false); - mesh.get(zt_prime, "backward_zt_prime", 0.0, false); - mesh.get(R_prime, "backward_R", 0.0, false); - mesh.get(Z_prime, "backward_Z", 0.0, false); - boundary = new BoundaryRegionPar("FCI_backward", BNDRY_PAR_BKWD, dir, &mesh); - } else { - // Definitely shouldn't be called - throw BoutException("FCIMap called with strange direction: %d. Only +/-1 currently supported.", dir); - } + const auto parallel_slice_field_name = [&](std::string field) -> std::string { + const std::string direction = (offset > 0) ? "forward" : "backward"; + // We only have a suffix for parallel slices beyond the first + // This is for backwards compatibility + const std::string slice_suffix = + (std::abs(offset) > 1) ? "_" + std::to_string(std::abs(offset)) : ""; + return direction + "_" + field + slice_suffix; + }; + + mesh.get(xt_prime, parallel_slice_field_name("xt_prime"), 0.0, false); + mesh.get(zt_prime, parallel_slice_field_name("zt_prime"), 0.0, false); + mesh.get(R_prime, parallel_slice_field_name("R"), 0.0, false); + mesh.get(Z_prime, parallel_slice_field_name("Z"), 0.0, false); - // Add the boundary region to the mesh's vector of parallel boundaries - mesh.addBoundaryPar(boundary); - // Cell corners Field3D xt_prime_corner(&mesh), zt_prime_corner(&mesh); xt_prime_corner.allocate(); @@ -215,7 +213,7 @@ FCIMap::FCIMap(Mesh &mesh, int dir, bool zperiodic) BoutReal dx = (dZ_dz * dR - dR_dz * dZ) / det; BoutReal dz = (dR_dx * dZ - dZ_dx * dR) / det; boundary->add_point(x, y, z, - x + dx, y + 0.5*dir, z + dz, // Intersection point in local index space + x + dx, y + 0.5*offset, z + dz, // Intersection point in local index space 0.5*coord.dy(x,y), //sqrt( SQ(dR) + SQ(dZ) ), // Distance to intersection PI // Right-angle intersection ); @@ -236,7 +234,7 @@ FCIMap::FCIMap(Mesh &mesh, int dir, bool zperiodic) interp->setMask(boundary_mask); } -const Field3D FCIMap::integrate(Field3D &f) const { +Field3D FCIMap::integrate(Field3D &f) const { TRACE("FCIMap::integrate"); // Cell centre values @@ -254,7 +252,7 @@ const Field3D FCIMap::integrate(Field3D &f) const { for(int x = mesh->xstart; x <= mesh->xend; x++) { for(int y = mesh->ystart; y <= mesh->yend; y++) { - int ynext = y+dir; + int ynext = y+offset; for(int z = 0; z < nz; z++) { if (boundary_mask(x,y,z)) @@ -294,24 +292,26 @@ const Field3D FCIMap::integrate(Field3D &f) const { return result; } -void FCITransform::calcYUpDown(Field3D &f) { +void FCITransform::calcYUpDown(Field3D& f) { TRACE("FCITransform::calcYUpDown"); // Ensure that yup and ydown are different fields f.splitYupYdown(); // Interpolate f onto yup and ydown fields - f.ynext(forward_map.dir) = forward_map.interpolate(f); - f.ynext(backward_map.dir) = backward_map.interpolate(f); + for (const auto& map : field_line_maps) { + f.ynext(map.offset) = map.interpolate(f); + } } -void FCITransform::integrateYUpDown(Field3D &f) { +void FCITransform::integrateYUpDown(Field3D& f) { TRACE("FCITransform::integrateYUpDown"); - + // Ensure that yup and ydown are different fields f.splitYupYdown(); // Integrate f onto yup and ydown fields - f.ynext(forward_map.dir) = forward_map.integrate(f); - f.ynext(backward_map.dir) = backward_map.integrate(f); + for (const auto& map : field_line_maps) { + f.ynext(map.offset) = map.integrate(f); + } } diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx index 293113b91d..765fc9fc83 100644 --- a/src/mesh/parallel/fci.hxx +++ b/src/mesh/parallel/fci.hxx @@ -32,6 +32,8 @@ #include #include +#include + /*! * Field line map - contains the coefficients for interpolation */ @@ -40,24 +42,21 @@ class FCIMap { Interpolation *interp; // Cell centre Interpolation *interp_corner; // Cell corner at (x+1, z+1) - /// Private constructor - must be initialised with mesh - FCIMap(); public: - /// dir MUST be either +1 or -1 - FCIMap(Mesh& mesh, int dir, bool zperiodic); + FCIMap() = delete; + FCIMap(Mesh& mesh, int offset, BoundaryRegionPar* boundary, bool zperiodic); - int dir; /**< Direction of map */ + /// Direction of map + const int offset; BoutMask boundary_mask; /**< boundary mask - has the field line left the domain */ BoutMask corner_boundary_mask; ///< If any of the integration area has left the domain Field3D y_prime; /**< distance to intersection with boundary */ - BoundaryRegionPar* boundary; /**< boundary region */ - - const Field3D interpolate(Field3D &f) const { return interp->interpolate(f); } + Field3D interpolate(Field3D &f) const { return interp->interpolate(f); } - const Field3D integrate(Field3D &f) const; + Field3D integrate(Field3D &f) const; }; /*! @@ -65,9 +64,21 @@ public: */ class FCITransform : public ParallelTransform { public: - FCITransform(Mesh &mesh, bool zperiodic = true) - : mesh(mesh), forward_map(mesh, +1, zperiodic), backward_map(mesh, -1, zperiodic), - zperiodic(zperiodic) {} + FCITransform() = delete; + FCITransform(Mesh& mesh, bool zperiodic = true) : mesh(mesh), zperiodic(zperiodic) { + auto forward_boundary = new BoundaryRegionPar("FCI_forward", BNDRY_PAR_FWD, +1, &mesh); + auto backward_boundary = new BoundaryRegionPar("FCI_backward", BNDRY_PAR_BKWD, -1, &mesh); + + // Add the boundary region to the mesh's vector of parallel boundaries + mesh.addBoundaryPar(forward_boundary); + mesh.addBoundaryPar(backward_boundary); + + field_line_maps.reserve(mesh.ystart * 2); + for (int offset = 1; offset < mesh.ystart + 1; ++offset) { + field_line_maps.emplace_back(mesh, offset, forward_boundary, zperiodic); + field_line_maps.emplace_back(mesh, -offset, backward_boundary, zperiodic); + } + } void calcYUpDown(Field3D &f) override; @@ -85,14 +96,13 @@ public: return false; } private: - FCITransform(); - Mesh& mesh; - FCIMap forward_map; /**< FCI map for field lines in +ve y */ - FCIMap backward_map; /**< FCI map for field lines in -ve y */ + /// FCI maps for field lines in +ve y + std::vector field_line_maps; - bool zperiodic; /**< Is the z-direction periodic? */ + /// Is the z-direction periodic? + bool zperiodic; }; #endif // __FCITRANSFORM_H__ From 13b68c12fd21d98a6378c260d30eaeef981587fe Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 25 Jan 2019 14:28:53 +0000 Subject: [PATCH 23/34] Throw if not enough parallel slices in FCI grid file --- src/mesh/parallel/fci.cxx | 56 ++++++++++++++++++++++++++------------- src/mesh/parallel/fci.hxx | 8 ++---- 2 files changed, 39 insertions(+), 25 deletions(-) diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 89558fc628..54b9152602 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -58,6 +58,8 @@ inline BoutReal sgn(BoutReal val) { return (BoutReal(0) < val) - (val < BoutReal FCIMap::FCIMap(Mesh &mesh, int offset_, BoundaryRegionPar* boundary, bool zperiodic) : offset(offset_), boundary_mask(mesh), corner_boundary_mask(mesh), y_prime(&mesh) { + TRACE("Creating FCIMAP for direction %d", offset); + if (offset == 0) { throw BoutException("FCIMap called with offset = 0; You probably didn't mean to do that"); } @@ -67,7 +69,7 @@ FCIMap::FCIMap(Mesh &mesh, int offset_, BoundaryRegionPar* boundary, bool zperio interp_corner = InterpolationFactory::getInstance()->create(&mesh); interp_corner->setYOffset(offset); - + // Index arrays contain guard cells in order to get subscripts right // x-index of bottom-left grid point auto i_corner = Tensor(mesh.LocalNx, mesh.LocalNy, mesh.LocalNz); @@ -91,10 +93,26 @@ FCIMap::FCIMap(Mesh &mesh, int offset_, BoundaryRegionPar* boundary, bool zperio return direction + "_" + field + slice_suffix; }; - mesh.get(xt_prime, parallel_slice_field_name("xt_prime"), 0.0, false); - mesh.get(zt_prime, parallel_slice_field_name("zt_prime"), 0.0, false); - mesh.get(R_prime, parallel_slice_field_name("R"), 0.0, false); - mesh.get(Z_prime, parallel_slice_field_name("Z"), 0.0, false); + if (mesh.get(xt_prime, parallel_slice_field_name("xt_prime"), 0.0, false) != 0) { + throw BoutException("Could not read %s from grid file!\n" + " Either add it to the grid file, or reduce MYG", + parallel_slice_field_name("xt_prime").c_str()); + } + if (mesh.get(zt_prime, parallel_slice_field_name("zt_prime"), 0.0, false) != 0) { + throw BoutException("Could not read %s from grid file!\n" + " Either add it to the grid file, or reduce MYG", + parallel_slice_field_name("zt_prime").c_str()); + } + if (mesh.get(R_prime, parallel_slice_field_name("R"), 0.0, false) != 0) { + throw BoutException("Could not read %s from grid file!\n" + " Either add it to the grid file, or reduce MYG", + parallel_slice_field_name("R").c_str()); + } + if (mesh.get(Z_prime, parallel_slice_field_name("Z"), 0.0, false) != 0) { + throw BoutException("Could not read %s from grid file!\n" + " Either add it to the grid file, or reduce MYG", + parallel_slice_field_name("Z").c_str()); + } // Cell corners Field3D xt_prime_corner(&mesh), zt_prime_corner(&mesh); @@ -110,12 +128,12 @@ FCIMap::FCIMap(Mesh &mesh, int offset_, BoundaryRegionPar* boundary, bool zperio (xt_prime(x + 1, y, z + 1) < 0.0) || (xt_prime(x, y, z + 1) < 0.0)) { // Hit a boundary corner_boundary_mask(x, y, z) = true; - + xt_prime_corner(x, y, z) = -1.0; zt_prime_corner(x, y, z) = -1.0; continue; } - + xt_prime_corner(x, y, z) = 0.25 * (xt_prime(x, y, z) + xt_prime(x + 1, y, z) + xt_prime(x, y, z + 1) + xt_prime(x + 1, y, z + 1)); @@ -126,12 +144,12 @@ FCIMap::FCIMap(Mesh &mesh, int offset_, BoundaryRegionPar* boundary, bool zperio } } } - + interp_corner->setMask(corner_boundary_mask); interp_corner->calcWeights(xt_prime_corner, zt_prime_corner); - + interp->calcWeights(xt_prime, zt_prime); - + int ncz = mesh.LocalNz; BoutReal t_x, t_z; @@ -212,7 +230,7 @@ FCIMap::FCIMap(Mesh &mesh, int offset_, BoundaryRegionPar* boundary, bool zperio // Invert 2x2 matrix to get change in index BoutReal dx = (dZ_dz * dR - dR_dz * dZ) / det; BoutReal dz = (dR_dx * dZ - dZ_dx * dR) / det; - boundary->add_point(x, y, z, + boundary->add_point(x, y, z, x + dx, y + 0.5*offset, z + dz, // Intersection point in local index space 0.5*coord.dy(x,y), //sqrt( SQ(dR) + SQ(dZ) ), // Distance to intersection PI // Right-angle intersection @@ -236,10 +254,10 @@ FCIMap::FCIMap(Mesh &mesh, int offset_, BoundaryRegionPar* boundary, bool zperio Field3D FCIMap::integrate(Field3D &f) const { TRACE("FCIMap::integrate"); - + // Cell centre values Field3D centre = interp->interpolate(f); - + // Cell corner values (x+1/2, z+1/2) Field3D corner = interp_corner->interpolate(f); @@ -248,23 +266,23 @@ Field3D FCIMap::integrate(Field3D &f) const { result.setLocation(f.getLocation()); int nz = mesh->LocalNz; - + for(int x = mesh->xstart; x <= mesh->xend; x++) { for(int y = mesh->ystart; y <= mesh->yend; y++) { - + int ynext = y+offset; - + for(int z = 0; z < nz; z++) { if (boundary_mask(x,y,z)) continue; - + int zm = z - 1; if (z == 0) { zm = nz-1; } - + BoutReal f_c = centre(x,ynext,z); - + if (corner_boundary_mask(x, y, z) || corner_boundary_mask(x - 1, y, z) || corner_boundary_mask(x, y, zm) || corner_boundary_mask(x - 1, y, zm) || (x == mesh->xstart)) { diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx index 765fc9fc83..95fa241996 100644 --- a/src/mesh/parallel/fci.hxx +++ b/src/mesh/parallel/fci.hxx @@ -65,7 +65,8 @@ public: class FCITransform : public ParallelTransform { public: FCITransform() = delete; - FCITransform(Mesh& mesh, bool zperiodic = true) : mesh(mesh), zperiodic(zperiodic) { + FCITransform(Mesh& mesh, bool zperiodic = true) { + auto forward_boundary = new BoundaryRegionPar("FCI_forward", BNDRY_PAR_FWD, +1, &mesh); auto backward_boundary = new BoundaryRegionPar("FCI_backward", BNDRY_PAR_BKWD, -1, &mesh); @@ -96,13 +97,8 @@ public: return false; } private: - Mesh& mesh; - /// FCI maps for field lines in +ve y std::vector field_line_maps; - - /// Is the z-direction periodic? - bool zperiodic; }; #endif // __FCITRANSFORM_H__ From 33155db7bdccc014edeba3883ed9b5581782a3e2 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 25 Jan 2019 17:36:36 +0000 Subject: [PATCH 24/34] Clean up FCI a bit - Remove unused variables - Store mesh and check it's the same as the field to be interpolated/integrated - Store interpolations as unique_ptrs to avoid potential memory leaks --- src/mesh/parallel/fci.cxx | 75 ++++++++++++++++++++------------------- src/mesh/parallel/fci.hxx | 40 +++++++++++---------- 2 files changed, 61 insertions(+), 54 deletions(-) diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 54b9152602..d85f7df1e4 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -48,15 +48,9 @@ #include -/** - * Return the sign of val - */ -inline BoutReal sgn(BoutReal val) { return (BoutReal(0) < val) - (val < BoutReal(0)); } - -// Calculate all the coefficients needed for the spline interpolation -// dir MUST be either +1 or -1 -FCIMap::FCIMap(Mesh &mesh, int offset_, BoundaryRegionPar* boundary, bool zperiodic) - : offset(offset_), boundary_mask(mesh), corner_boundary_mask(mesh), y_prime(&mesh) { +FCIMap::FCIMap(Mesh& mesh, int offset_, BoundaryRegionPar* boundary, bool zperiodic) + : map_mesh(mesh), offset(offset_), boundary_mask(map_mesh), + corner_boundary_mask(map_mesh) { TRACE("Creating FCIMAP for direction %d", offset); @@ -64,26 +58,31 @@ FCIMap::FCIMap(Mesh &mesh, int offset_, BoundaryRegionPar* boundary, bool zperio throw BoutException("FCIMap called with offset = 0; You probably didn't mean to do that"); } - interp = InterpolationFactory::getInstance()->create(&mesh); + interp = + std::unique_ptr(InterpolationFactory::getInstance()->create(&map_mesh)); interp->setYOffset(offset); - interp_corner = InterpolationFactory::getInstance()->create(&mesh); + interp_corner = + std::unique_ptr(InterpolationFactory::getInstance()->create(&map_mesh)); interp_corner->setYOffset(offset); // Index arrays contain guard cells in order to get subscripts right // x-index of bottom-left grid point - auto i_corner = Tensor(mesh.LocalNx, mesh.LocalNy, mesh.LocalNz); + auto i_corner = Tensor(map_mesh.LocalNx, map_mesh.LocalNy, map_mesh.LocalNz); // z-index of bottom-left grid point - auto k_corner = Tensor(mesh.LocalNx, mesh.LocalNy, mesh.LocalNz); + auto k_corner = Tensor(map_mesh.LocalNx, map_mesh.LocalNy, map_mesh.LocalNz); - Field3D xt_prime(&mesh), zt_prime(&mesh); - Field3D R(&mesh), Z(&mesh); // Real-space coordinates of grid points - Field3D R_prime(&mesh), - Z_prime(&mesh); // Real-space coordinates of forward/backward points + // Index-space coordinates of forward/backward points + Field3D xt_prime(&map_mesh), zt_prime(&map_mesh); + // Real-space coordinates of grid points + Field3D R(&map_mesh), Z(&map_mesh); + // Real-space coordinates of forward/backward points + Field3D R_prime(&map_mesh), Z_prime(&map_mesh); - mesh.get(R, "R", 0.0, false); - mesh.get(Z, "Z", 0.0, false); + map_mesh.get(R, "R", 0.0, false); + map_mesh.get(Z, "Z", 0.0, false); + // Get a unique name for a field based on the sign/magnitude of the offset const auto parallel_slice_field_name = [&](std::string field) -> std::string { const std::string direction = (offset > 0) ? "forward" : "backward"; // We only have a suffix for parallel slices beyond the first @@ -93,35 +92,37 @@ FCIMap::FCIMap(Mesh &mesh, int offset_, BoundaryRegionPar* boundary, bool zperio return direction + "_" + field + slice_suffix; }; - if (mesh.get(xt_prime, parallel_slice_field_name("xt_prime"), 0.0, false) != 0) { + // If we can't read in any of these fields, things will silently not + // work, so best throw + if (map_mesh.get(xt_prime, parallel_slice_field_name("xt_prime"), 0.0, false) != 0) { throw BoutException("Could not read %s from grid file!\n" " Either add it to the grid file, or reduce MYG", parallel_slice_field_name("xt_prime").c_str()); } - if (mesh.get(zt_prime, parallel_slice_field_name("zt_prime"), 0.0, false) != 0) { + if (map_mesh.get(zt_prime, parallel_slice_field_name("zt_prime"), 0.0, false) != 0) { throw BoutException("Could not read %s from grid file!\n" " Either add it to the grid file, or reduce MYG", parallel_slice_field_name("zt_prime").c_str()); } - if (mesh.get(R_prime, parallel_slice_field_name("R"), 0.0, false) != 0) { + if (map_mesh.get(R_prime, parallel_slice_field_name("R"), 0.0, false) != 0) { throw BoutException("Could not read %s from grid file!\n" " Either add it to the grid file, or reduce MYG", parallel_slice_field_name("R").c_str()); } - if (mesh.get(Z_prime, parallel_slice_field_name("Z"), 0.0, false) != 0) { + if (map_mesh.get(Z_prime, parallel_slice_field_name("Z"), 0.0, false) != 0) { throw BoutException("Could not read %s from grid file!\n" " Either add it to the grid file, or reduce MYG", parallel_slice_field_name("Z").c_str()); } // Cell corners - Field3D xt_prime_corner(&mesh), zt_prime_corner(&mesh); + Field3D xt_prime_corner(&map_mesh), zt_prime_corner(&map_mesh); xt_prime_corner.allocate(); zt_prime_corner.allocate(); - for (int x = mesh.xstart; x <= mesh.xend; x++) { - for (int y = mesh.ystart; y <= mesh.yend; y++) { - for (int z = 0; z < mesh.LocalNz - 1; z++) { + for (int x = map_mesh.xstart; x <= map_mesh.xend; x++) { + for (int y = map_mesh.ystart; y <= map_mesh.yend; y++) { + for (int z = 0; z < map_mesh.LocalNz - 1; z++) { // Point interpolated from (x+1/2, z+1/2) if ((xt_prime(x, y, z) < 0.0) || (xt_prime(x + 1, y, z) < 0.0) || @@ -150,13 +151,13 @@ FCIMap::FCIMap(Mesh &mesh, int offset_, BoundaryRegionPar* boundary, bool zperio interp->calcWeights(xt_prime, zt_prime); - int ncz = mesh.LocalNz; + int ncz = map_mesh.LocalNz; BoutReal t_x, t_z; - Coordinates &coord = *(mesh.getCoordinates()); + Coordinates &coord = *(map_mesh.getCoordinates()); - for (int x = mesh.xstart; x <= mesh.xend; x++) { - for (int y = mesh.ystart; y <= mesh.yend; y++) { + for (int x = map_mesh.xstart; x <= map_mesh.xend; x++) { + for (int y = map_mesh.ystart; y <= map_mesh.yend; y++) { for (int z = 0; z < ncz; z++) { // The integer part of xt_prime, zt_prime are the indices of the cell @@ -213,7 +214,7 @@ FCIMap::FCIMap(Mesh &mesh, int offset_, BoundaryRegionPar* boundary, bool zperio dR_dz = R(x, y, z + 1) - R(x, y, z); dZ_dz = Z(x, y, z + 1) - Z(x, y, z); - } else if (z == mesh.LocalNz - 1) { + } else if (z == map_mesh.LocalNz - 1) { dR_dz = R(x, y, z) - R(x, y, z - 1); dZ_dz = Z(x, y, z) - Z(x, y, z - 1); @@ -255,6 +256,8 @@ FCIMap::FCIMap(Mesh &mesh, int offset_, BoundaryRegionPar* boundary, bool zperio Field3D FCIMap::integrate(Field3D &f) const { TRACE("FCIMap::integrate"); + ASSERT3(&map_mesh == f.getMesh()); + // Cell centre values Field3D centre = interp->interpolate(f); @@ -265,10 +268,10 @@ Field3D FCIMap::integrate(Field3D &f) const { result.allocate(); result.setLocation(f.getLocation()); - int nz = mesh->LocalNz; + int nz = map_mesh.LocalNz; - for(int x = mesh->xstart; x <= mesh->xend; x++) { - for(int y = mesh->ystart; y <= mesh->yend; y++) { + for(int x = map_mesh.xstart; x <= map_mesh.xend; x++) { + for(int y = map_mesh.ystart; y <= map_mesh.yend; y++) { int ynext = y+offset; @@ -285,7 +288,7 @@ Field3D FCIMap::integrate(Field3D &f) const { if (corner_boundary_mask(x, y, z) || corner_boundary_mask(x - 1, y, z) || corner_boundary_mask(x, y, zm) || corner_boundary_mask(x - 1, y, zm) || - (x == mesh->xstart)) { + (x == map_mesh.xstart)) { // One of the corners leaves the domain. // Use the cell centre value, since boundary conditions are not // currently applied to corners. diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx index 95fa241996..37a87a96e6 100644 --- a/src/mesh/parallel/fci.hxx +++ b/src/mesh/parallel/fci.hxx @@ -32,36 +32,41 @@ #include #include +#include #include -/*! - * Field line map - contains the coefficients for interpolation - */ + +/// Field line map - contains the coefficients for interpolation class FCIMap { - /// Interpolation object - Interpolation *interp; // Cell centre - Interpolation *interp_corner; // Cell corner at (x+1, z+1) + /// Interpolation objects + std::unique_ptr interp; // Cell centre + std::unique_ptr interp_corner; // Cell corner at (x+1, z+1) public: FCIMap() = delete; FCIMap(Mesh& mesh, int offset, BoundaryRegionPar* boundary, bool zperiodic); + // The mesh this map was created on + Mesh& map_mesh; + /// Direction of map const int offset; - BoutMask boundary_mask; /**< boundary mask - has the field line left the domain */ - BoutMask corner_boundary_mask; ///< If any of the integration area has left the domain + /// boundary mask - has the field line left the domain + BoutMask boundary_mask; + /// If any of the integration area has left the domain + BoutMask corner_boundary_mask; - Field3D y_prime; /**< distance to intersection with boundary */ - - Field3D interpolate(Field3D &f) const { return interp->interpolate(f); } + Field3D interpolate(Field3D& f) const { + ASSERT3(&map_mesh == f.getMesh()); + return interp->interpolate(f); + } Field3D integrate(Field3D &f) const; }; -/*! - * Flux Coordinate Independent method for parallel derivatives - */ + +/// Flux Coordinate Independent method for parallel derivatives class FCITransform : public ParallelTransform { public: FCITransform() = delete; @@ -93,11 +98,10 @@ public: throw BoutException("FCI method cannot transform into field aligned grid"); } - bool canToFromFieldAligned() override{ - return false; - } + bool canToFromFieldAligned() override { return false; } + private: - /// FCI maps for field lines in +ve y + /// FCI maps for each of the parallel slices std::vector field_line_maps; }; From e743a525b6917eb150f4c57985f41a1319839447 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 25 Jan 2019 17:39:47 +0000 Subject: [PATCH 25/34] Remove old fci-slab test (replaced with better MMS test) --- tests/integrated/test-fci-slab/data/BOUT.inp | 26 --- tests/integrated/test-fci-slab/fci.py | 110 ------------ tests/integrated/test-fci-slab/fci_slab.cxx | 53 ------ tests/integrated/test-fci-slab/generate.py | 126 ------------- tests/integrated/test-fci-slab/makefile | 6 - tests/integrated/test-fci-slab/mms.py | 47 ----- tests/integrated/test-fci-slab/mms/BOUT.inp | 39 ----- tests/integrated/test-fci-slab/plot_funcs.py | 28 --- tests/integrated/test-fci-slab/plot_interp.py | 17 -- tests/integrated/test-fci-slab/runtest | 165 ------------------ tests/integrated/test-fci-slab/simple_test.py | 32 ---- 11 files changed, 649 deletions(-) delete mode 100644 tests/integrated/test-fci-slab/data/BOUT.inp delete mode 100644 tests/integrated/test-fci-slab/fci.py delete mode 100644 tests/integrated/test-fci-slab/fci_slab.cxx delete mode 100644 tests/integrated/test-fci-slab/generate.py delete mode 100644 tests/integrated/test-fci-slab/makefile delete mode 100644 tests/integrated/test-fci-slab/mms.py delete mode 100644 tests/integrated/test-fci-slab/mms/BOUT.inp delete mode 100644 tests/integrated/test-fci-slab/plot_funcs.py delete mode 100644 tests/integrated/test-fci-slab/plot_interp.py delete mode 100755 tests/integrated/test-fci-slab/runtest delete mode 100644 tests/integrated/test-fci-slab/simple_test.py diff --git a/tests/integrated/test-fci-slab/data/BOUT.inp b/tests/integrated/test-fci-slab/data/BOUT.inp deleted file mode 100644 index 218946e676..0000000000 --- a/tests/integrated/test-fci-slab/data/BOUT.inp +++ /dev/null @@ -1,26 +0,0 @@ -grid = fci.grid.nc -#grid = simple_test.nc - -nout = 50 -timestep = 0.2 - -MZ = 64 - -[mesh] -paralleltransform = fci - -symmetricglobalx = true - -[interpolation] -type=lagrange4pt - -[fci] -y_periodic = false -z_periodic = false - -[f] -scale = 1.0 -function = cos(y-z) - -bndry_par_yup = parallel_dirichlet(0.0) -bndry_par_ydown = parallel_dirichlet(0.0) diff --git a/tests/integrated/test-fci-slab/fci.py b/tests/integrated/test-fci-slab/fci.py deleted file mode 100644 index 1c505ffa20..0000000000 --- a/tests/integrated/test-fci-slab/fci.py +++ /dev/null @@ -1,110 +0,0 @@ -from __future__ import division -from builtins import object -from past.utils import old_div -import numpy as np -from math import pi -from scipy.integrate import odeint -import boututils.datafile as bdata -from boutdata.input import transform3D - -# Parameters -nx = 34 -########## y is toroidal! -ny = 64 -########## z is poloidal! -nz = 64 - -Lx = 0.1 # Radial domain size [m] -Ltor = 10. # "Toroidal" length [m] -Lpol = 1. # "Poloidal" length [m] - -delta_x = old_div(Lx,(nx)) -delta_pol = old_div(Lpol,(nz)) -delta_tor = old_div(Ltor,(ny)) - -Bt = 1.0 # Magnetic field [T] -Bp = 0.1 # Poloidal field at the middle of the domain [T] -Bpprime = 1.0 # Bp gradient [T/m] Bp(x) = Bp + Bpprime * x - -# Coord arrays -x = np.linspace(0,Lx,nx) -y = np.linspace(0,Ltor,ny) -z = np.linspace(0,Lpol,nz,endpoint=False) - -############################################################ - -# Effective major radius -R = old_div(Ltor, (2.*pi)) - -# Set poloidal magnetic field - -Bpx = Bp + (x-old_div(Lx,2)) * Bpprime - -Bpxy = np.transpose(np.resize(Bpx, (nz, ny, nx)), (2,1,0)) - -Bxy = np.sqrt(Bpxy**2 + Bt**2) - -############################################################ - -class Mappoint(object): - def __init__(self, xt, zt): - self.xt = xt - self.zt = zt - - self.xt_prime = old_div(xt,delta_x) - self.zt_prime = old_div(zt,delta_pol) - -def unroll_map_coeff(map_list, coeff): - coeff_array = np.transpose(np.resize(np.array([getattr(f, coeff) for f in map_list]).reshape( (nx,nz) ), (ny, nx, nz) ), (1, 0, 2) ) - return coeff_array - -def b_field(vector, y): - x0 = 0.05 # Centre of box, where bz = 0. - x, z = vector; - bx = 0. - bz = Bp + (x-x0) * Bpprime - - return [bx, bz] - -def field_line_tracer(direction, map_list): - - result = np.zeros( (nx, nz, 2) ) - - for i in np.arange(0,nx): - for k in np.arange(0,nz): - result[i,k,:] = odeint(b_field, [x[i], z[k]], [0, delta_tor*direction])[1,:] - result[i,k,1] = np.mod(result[i,k,1], Lpol) - - map_list.append(Mappoint(result[i,k,0],result[i,k,1])) - - return result - -if __name__ == "__main__": - - forward_map = [] - forward_coords = field_line_tracer(+1, forward_map) - backward_map = [] - backward_coords = field_line_tracer(-1, backward_map) - - X,Y = np.meshgrid(x,y,indexing='ij') - x0 = 0.5 - g_22 = np.sqrt(((Bp + (X-x0) * Lx * Bpprime)**2 + 1)) - - with bdata.DataFile('fci.grid.nc', write=True, create=True) as f: - f.write('nx', nx) - f.write('ny', ny) - f.write('nz', nz) - f.write("dx", delta_x) - f.write("dy", delta_tor) - f.write("g_22", g_22) - f.write("Bxy", transform3D(Bxy)) - - xt_prime = unroll_map_coeff(forward_map, 'xt_prime') - f.write('forward_xt_prime', transform3D(xt_prime)) - zt_prime = unroll_map_coeff(forward_map, 'zt_prime') - f.write('forward_zt_prime', transform3D(zt_prime)) - - xt_prime = unroll_map_coeff(backward_map, 'xt_prime') - f.write('backward_xt_prime', transform3D(xt_prime)) - zt_prime = unroll_map_coeff(backward_map, 'zt_prime') - f.write('backward_zt_prime', transform3D(zt_prime)) diff --git a/tests/integrated/test-fci-slab/fci_slab.cxx b/tests/integrated/test-fci-slab/fci_slab.cxx deleted file mode 100644 index 60a6bfd4ca..0000000000 --- a/tests/integrated/test-fci-slab/fci_slab.cxx +++ /dev/null @@ -1,53 +0,0 @@ -#include -#include -#include - -class FCISlab : public PhysicsModel { -public: - - // We need to initialise the FCI object with the mesh - FCISlab() {} - - int init(bool UNUSED(restarting)) { - - D = 10; - - Coordinates *coord = mesh->getCoordinates(); - - mesh->get(coord->g_22, "g_22"); - - coord->geometry(); - - solver->add(f, "f"); - solver->add(g, "g"); - - f.applyBoundary("dirichlet"); - g.applyBoundary("dirichlet"); - - return 0; - } - - int rhs(BoutReal time); - -private: - Field3D f, g; - - BoutReal D; -}; - -BOUTMAIN(FCISlab); - -int FCISlab::rhs(BoutReal time) { - mesh->communicate(f,g); - - Coordinates *coord = mesh->getCoordinates(); - - f.applyParallelBoundary(time); - g.applyParallelBoundary(time); - - ddt(f) = Grad_par(g) + D*SQ(coord->dy)*Grad2_par2(f); - - ddt(g) = Grad_par(f) + D*SQ(coord->dy)*Grad2_par2(g); - - return 0; -} diff --git a/tests/integrated/test-fci-slab/generate.py b/tests/integrated/test-fci-slab/generate.py deleted file mode 100644 index 03767db3a6..0000000000 --- a/tests/integrated/test-fci-slab/generate.py +++ /dev/null @@ -1,126 +0,0 @@ -from __future__ import division -from builtins import object -from past.utils import old_div -# -# Routines to generate slab meshes for FCI -# - -import numpy as np -from math import pi -from scipy.integrate import odeint -import boututils.datafile as bdata - - -def slab(nx, ny, nz, - filename="fci.grid.nc", - Lx=0.1, Ly=10., Lz = 1., - Bt=1.0, Bp = 0.1, Bpprime = 1.0): - """ - nx - Number of radial points - ny - Number of toroidal points (NOTE: Different to BOUT++ standard) - nz - Number of poloidal points - - Lx - Radial domain size [m] - Ly - Toroidal domain size [m] - Lz - Poloidal domain size [m] - - Bt - Toroidal magnetic field [T] - Bp - Poloidal magnetic field [T] - Bpprime - Gradient of Bp [T/m] Bp(x) = Bp + Bpprime * x - """ - - MXG = 2 - - # Make sure input types are sane - nx = int(nx) - ny = int(ny) - nz = int(nz) - - Lx = float(Lx) - Ly = float(Ly) - Lz = float(Lz) - - delta_x = old_div(Lx,(nx-2.*MXG)) - delta_pol = old_div(Lz,(nz)) - delta_tor = old_div(Ly,(ny)) - - # Coord arrays - x = Lx * (np.arange(nx) - MXG + 0.5)/(nx - 2.*MXG) # 0 and 1 half-way between cells - y = np.linspace(0,Ly,ny) - z = np.linspace(0,Lz,nz,endpoint=False) - - ############################################################ - - # Effective major radius - R = old_div(Ly, (2.*pi)) - - # Set poloidal magnetic field - - Bpx = Bp + (x-old_div(Lx,2)) * Bpprime - - Bpxy = np.transpose(np.resize(Bpx, (nz, ny, nx)), (2,1,0)) - - Bxy = np.sqrt(Bpxy**2 + Bt**2)[:,:,0] - - class Mappoint(object): - def __init__(self, xt, zt): - self.xt = xt - self.zt = zt - - self.xt_prime = old_div(xt,delta_x) + MXG - 0.5 - self.zt_prime = old_div(zt,delta_pol) - - def unroll_map_coeff(map_list, coeff): - coeff_array = np.transpose(np.resize(np.array([getattr(f, coeff) for f in map_list]).reshape( (nx,nz) ), (ny, nx, nz) ), (1, 0, 2) ) - return coeff_array - - def b_field(vector, y): - x0 = old_div(Lx,2.) # Centre of box, where bz = 0. - x, z = vector; - bx = 0. - bz = Bp + (x-x0) * Bpprime - - return [bx, bz] - - def field_line_tracer(direction, map_list): - - result = np.zeros( (nx, nz, 2) ) - - for i in np.arange(0,nx): - for k in np.arange(0,nz): - result[i,k,:] = odeint(b_field, [x[i], z[k]], [0, delta_tor*direction])[1,:] - map_list.append(Mappoint(result[i,k,0],result[i,k,1])) - - return result - - forward_map = [] - forward_coords = field_line_tracer(+1, forward_map) - backward_map = [] - backward_coords = field_line_tracer(-1, backward_map) - - X,Y = np.meshgrid(x,y,indexing='ij') - x0 = 0.5 - g_22 = old_div(((Bp + (X-x0) * Lx * Bpprime)**2 + Bt**2), Bt**2) - - with bdata.DataFile(filename, write=True, create=True) as f: - f.write('nx', nx) - f.write('ny', ny) - f.write('nz', nz) - f.write("dx", delta_x) - f.write("dy", delta_tor) - f.write("g_22", g_22) - f.write("Bxy", (Bxy)) - - xt_prime = unroll_map_coeff(forward_map, 'xt_prime') - f.write('forward_xt_prime', (xt_prime)) - zt_prime = unroll_map_coeff(forward_map, 'zt_prime') - f.write('forward_zt_prime', (zt_prime)) - - xt_prime = unroll_map_coeff(backward_map, 'xt_prime') - f.write('backward_xt_prime', (xt_prime)) - zt_prime = unroll_map_coeff(backward_map, 'zt_prime') - f.write('backward_zt_prime', (zt_prime)) - - -if __name__ == "__main__": - slab(34, 64, 64, filename="fci.grid.nc") diff --git a/tests/integrated/test-fci-slab/makefile b/tests/integrated/test-fci-slab/makefile deleted file mode 100644 index b0fbe93385..0000000000 --- a/tests/integrated/test-fci-slab/makefile +++ /dev/null @@ -1,6 +0,0 @@ - -BOUT_TOP = ../../.. - -SOURCEC = fci_slab.cxx - -include $(BOUT_TOP)/make.config diff --git a/tests/integrated/test-fci-slab/mms.py b/tests/integrated/test-fci-slab/mms.py deleted file mode 100644 index df00c5d68d..0000000000 --- a/tests/integrated/test-fci-slab/mms.py +++ /dev/null @@ -1,47 +0,0 @@ -# -# Generate manufactured solution and sources for FCI test -# - -from __future__ import division -from __future__ import print_function - -from boutdata.mms import * - -from sympy import sin, cos, sqrt - -from math import pi - -f = sin(y - z) + cos(t)*sin(y - 2*z) - -g = cos(y - z) - cos(t)*sin(y - 2*z) - -Lx = 0.1 -Ly = 10. -Lz = 1. - -Bt = 1.0 -Bp = 0.05 -Bpprime = 0.1 - -Bpx = Bp + (x-0.5)*Lx * Bpprime # Note: x in range [0,1] -B = sqrt(Bpx**2 + Bt**2) - -def FCI_Grad_par(f): - return ( Bt * diff(f, y)*2.*pi/Ly + Bpx * diff(f, z)*2.*pi/Lz ) / B - -############################################ -# Equations solved - -dfdt = FCI_Grad_par(g) -dgdt = FCI_Grad_par(f) - -# Loop over variables and print solution, source etc. -for v, dvdt, name in [ (f, dfdt, "f"), (g, dgdt, "g") ]: - # Calculate source - S = diff(v, t) - dvdt - - print("\n["+name+"]") - print("solution = "+exprToStr(v)) - print("\nsource = "+exprToStr(S)) - print("\nbndry_par_all = parallel_dirichlet("+name+":solution)") - diff --git a/tests/integrated/test-fci-slab/mms/BOUT.inp b/tests/integrated/test-fci-slab/mms/BOUT.inp deleted file mode 100644 index 3c1e5f9559..0000000000 --- a/tests/integrated/test-fci-slab/mms/BOUT.inp +++ /dev/null @@ -1,39 +0,0 @@ -grid = fci.grid.nc - -nout = 1 -timestep = 0.01 - -MZ = 64 - -NXPE = 1 - -[mesh] -paralleltransform = fci - -symmetricglobalx = true - -[fci] -y_periodic = false -z_periodic = false - -[interpolation] -type=lagrange4pt - -[solver] -ATOL = 1e-12 -RTOL = 1e-8 -mms = true - -[f] -solution = sin(y - 2*z)*cos(t) + sin(y - z) - -source = -sin(t)*sin(y - 2*z) - (6.28318530717959*(0.01*x + 0.045)*(sin(y - z) + 2*cos(t)*cos(y - 2*z)) - 0.628318530717959*sin(y - z) - 0.628318530717959*cos(t)*cos(y - 2*z))/sqrt((0.01*x + 0.045)^2 + 1.0) - -bndry_par_all = parallel_dirichlet(f:solution) - -[g] -solution = -sin(y - 2*z)*cos(t) + cos(y - z) - -source = sin(t)*sin(y - 2*z) - (6.28318530717959*(0.01*x + 0.045)*(-2*cos(t)*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(t)*cos(y - 2*z) + 0.628318530717959*cos(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) - -bndry_par_all = parallel_dirichlet(g:solution) diff --git a/tests/integrated/test-fci-slab/plot_funcs.py b/tests/integrated/test-fci-slab/plot_funcs.py deleted file mode 100644 index d233858f8f..0000000000 --- a/tests/integrated/test-fci-slab/plot_funcs.py +++ /dev/null @@ -1,28 +0,0 @@ -from builtins import str -from builtins import range -# Plot interpolating functions -# Input generated by simple_test.py - -from numpy import linspace -import matplotlib.pyplot as plt -from boutdata.collect import collect - -f = collect("f", path="data") -yup = collect("yup", path="data") - -ny = 20 -nz = 8 - -# Note: yup[y=0] is never set -y = linspace(-1, 1, ny-1) - -plt.plot(f[0,4,4,:], 'o', label="f") - -for z in range(nz): - plt.plot(y+z, yup[4,1:,z], label="z = "+str(z)) - -plt.legend(loc='upper center') - -plt.savefig("plot_funcs.pdf") - -plt.show() diff --git a/tests/integrated/test-fci-slab/plot_interp.py b/tests/integrated/test-fci-slab/plot_interp.py deleted file mode 100644 index 2278684fba..0000000000 --- a/tests/integrated/test-fci-slab/plot_interp.py +++ /dev/null @@ -1,17 +0,0 @@ - -import matplotlib.pyplot as plt - -from boutdata.collect import collect - -f = collect("f", path="data") -yup = collect("yup", path="data") -ydown = collect("ydown", path="data") - -plt.plot(f[0,4,4,:], label="f") -plt.plot(yup[4,4,:], label="f.yup") -plt.plot(ydown[4,4,:], label="f.ydown") - -plt.legend() - -plt.savefig("plot_interp.pdf") -plt.show() diff --git a/tests/integrated/test-fci-slab/runtest b/tests/integrated/test-fci-slab/runtest deleted file mode 100755 index 196b2471f0..0000000000 --- a/tests/integrated/test-fci-slab/runtest +++ /dev/null @@ -1,165 +0,0 @@ -#!/usr/bin/env python3 -# -# Python script to run and analyse MMS test -# -from __future__ import division -from __future__ import print_function -from builtins import zip -from builtins import str - -from boututils.run_wrapper import shell, shell_safe, launch_safe, getmpirun -from boututils.datafile import DataFile -from boutdata.collect import collect - -from numpy import sqrt, max, abs, mean, array, log, pi, polyfit, linspace, arange - -import pickle - -from sys import stdout - -import zoidberg as zb - -showPlot = False #Do we want to show the plot as well as save it to file. - -nx = 5 # Not changed for these tests - -# Resolution in y and z -nlist = [64,128] #[8,16,32,64,128,256] - -nproc = 2 - -directory = "mms" - -varlist = ["f", "g"] -markers = ['bo', 'r^'] -labels = [r'$f$', r'$g$'] - -MPIRUN = getmpirun() - -success=True - -print("Making fci-slab test") -shell_safe("make > make.log") - -error_2 = {} -error_inf = {} -for var in varlist: - error_2[var] = [] # The L2 error (RMS) - error_inf[var] = [] # The maximum error - -yperiodic=False # Run with periodic Y? - -for n in nlist: - - # Define the magnetic field using new poloidal gridding method - # Note that the Bz and Bzprime parameters here must be the same as in mms.py - field = zb.field.Slab(Bz=0.05, Bzprime=0.1) - # Create rectangular poloidal grids - poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid(nx,n,1.,1.) - # Set the ylength and y locations - ylength = 10. - - if yperiodic: - ycoords = linspace(0.0, ylength, n, endpoint=False) - else: - # Doesn't include the end points - ycoords = (arange(n) + 0.5)*ylength/float(n) - - # Create the grid - grid = zb.grid.Grid(poloidal_grid, ycoords, ylength, yperiodic=yperiodic) - # Make and write maps - maps = zb.make_maps(grid, field) - zb.write_maps(grid, field, maps, new_names=False, metric2d=True) - - args = " -d "+directory+" MZ="+str(n)+ " fci:y_periodic="+str(yperiodic) - - # Command to run - cmd = "./fci_slab "+args - - print("Running command: "+cmd) - - # Launch using MPI - s, out = launch_safe(cmd, runcmd=MPIRUN, nproc=nproc, mthread=1, pipe=True) - - # Save output to log file - with open("run.log."+str(n), "w") as f: - f.write(out) - - if s: - print("Run failed!\nOutput was:\n") - print(out) - exit(s) - - for var in varlist: - # Collect data - E = collect("E_"+var, tind=[1,1], info=False, path=directory) - E = E[:,2:-2, :,:] - - # Average error over domain - l2 = sqrt(mean(E**2)) - linf = max(abs( E )) - - error_2[var].append( l2 ) - error_inf[var].append( linf ) - - print("%s : l-2 %f l-inf %f" % (var, l2, linf)) - -dx = 1. / array(nlist) - -# Save data -with open("fci_mms.pkl", "wb") as output: - pickle.dump(nlist, output) - pickle.dump(error_2, output) - pickle.dump(error_inf, output) - -# Calculate convergence order -for var,mark,label in zip(varlist, markers, labels): - fit = polyfit(log(dx), log(error_2[var]), 1) - order = fit[0] - stdout.write("%s Convergence order = %f (fit)" % (var, order)) - - order = log(error_2[var][-2]/error_2[var][-1])/log(dx[-2]/dx[-1]) - stdout.write(", %f (small spacing)" % (order,)) - - if order > 1.5: # Should be second order accurate - print("............ PASS") - else: - print("............ FAIL") - success = False - -if False: - try: - # Plot using matplotlib if available - import matplotlib.pyplot as plt - - plt.figure() - - for var,mark,label in zip(varlist, markers, labels): - plt.plot(dx, error_2[var], '-'+mark, label=label) - plt.plot(dx, error_inf[var], '--'+mark) - - plt.legend(loc="upper left") - plt.grid() - - plt.yscale('log') - plt.xscale('log') - - plt.xlabel(r'Mesh spacing $\delta x$') - plt.ylabel("Error norm") - - plt.savefig("fci-norm.pdf") - - print("Plot saved to fci-norm.pdf") - - if showPlot: - plt.show() - plt.close() - except ImportError: - print("No matplotlib") -else: - print("Plotting disabled") - -if success: - exit(0) -else: - exit(1) diff --git a/tests/integrated/test-fci-slab/simple_test.py b/tests/integrated/test-fci-slab/simple_test.py deleted file mode 100644 index 95f245c8c2..0000000000 --- a/tests/integrated/test-fci-slab/simple_test.py +++ /dev/null @@ -1,32 +0,0 @@ -from builtins import range -from numpy import zeros, linspace, concatenate -import boututils.datafile as bdata -from boutdata.input import transform3D - -# Parameters -nx = 10 -ny = 20 -nz = 8 - -shape = [nx, ny, nz] - -xt_prime = zeros(shape) -zt_prime = zeros(shape) - -for x in range(nx): - # No interpolation in x - xt_prime[x,:,:] = x - - # Each y slice scans between neighbouring z points - for z in range(nz): - zt_prime[x,:,z] = z + concatenate([linspace(-1, 1, ny-1), [0]]) - - -with bdata.DataFile('simple_test.nc', write=True, create=True) as f: - f.write('nx',nx) - f.write('ny',ny) - - for direction_name in ['forward', 'backward']: - f.write(direction_name + '_xt_prime', transform3D(xt_prime)) - f.write(direction_name + '_zt_prime', transform3D(zt_prime)) - From 6d95e434272794b6f32ec3f614ce17a6ad181e53 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 25 Jan 2019 17:42:50 +0000 Subject: [PATCH 26/34] Fix outdated comments --- include/field3d.hxx | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/include/field3d.hxx b/include/field3d.hxx index 2274df81a6..f9b23400ba 100644 --- a/include/field3d.hxx +++ b/include/field3d.hxx @@ -235,7 +235,7 @@ class Field3D : public Field, public FieldData { /// Return reference to yup field Field3D &yup(std::vector::size_type index = 0) { - ASSERT2(index < yup_fields.size()); // Check for communicate + ASSERT2(index < yup_fields.size()); return yup_fields[index]; } /// Return const reference to yup field @@ -256,9 +256,11 @@ class Field3D : public Field, public FieldData { return ydown_fields[index]; } - /// Return yup if dir=+1, and ydown if dir=-1 - Field3D& ynext(int dir); - const Field3D& ynext(int dir) const; + /// Return the parallel slice at \p offset + /// + /// \p offset of 0 returns the main field itself + Field3D& ynext(int offset); + const Field3D& ynext(int offset) const; /// Set variable location for staggered grids to @param new_location /// From 49de89c40459071d45940899fd78ee0c03259cef Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 29 Jan 2019 10:20:14 +0000 Subject: [PATCH 27/34] Reorder namedtuple arguments for compatibility with Python < 3.5 --- tools/pylib/zoidberg/zoidberg.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/pylib/zoidberg/zoidberg.py b/tools/pylib/zoidberg/zoidberg.py index 3217b70c76..5833d01134 100644 --- a/tools/pylib/zoidberg/zoidberg.py +++ b/tools/pylib/zoidberg/zoidberg.py @@ -88,7 +88,7 @@ def make_maps(grid, magnetic_field, nslice=1, quiet=False, **kwargs): } # A helper data structure that groups the various field line maps along with the offset - ParallelSlice = namedtuple('ParallelSlice', ['R', 'Z', 'xt_prime', 'zt_prime', 'offset']) + ParallelSlice = namedtuple('ParallelSlice', ['offset', 'R', 'Z', 'xt_prime', 'zt_prime']) # A list of the above data structures for each offset we want parallel_slices = [] @@ -104,7 +104,7 @@ def make_maps(grid, magnetic_field, nslice=1, quiet=False, **kwargs): # Get the field arrays we just made and wrap them up in our helper tuple fields = map(lambda x: maps[x], field_names) - parallel_slices.append(ParallelSlice(*fields, offset)) + parallel_slices.append(ParallelSlice(offset, *fields)) # Total size of the progress bar total_work = float((len(parallel_slices) - 1) * (ny-1)) From 1e3b5469f03f6788dd42189f489998bd1d72d689 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 29 Jan 2019 10:34:11 +0000 Subject: [PATCH 28/34] Add some runtime checks for parallel slice consistency --- src/field/field3d.cxx | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index c0c85b193a..43b2768bf9 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -135,6 +135,13 @@ Field3D* Field3D::timeDeriv() { void Field3D::splitYupYdown() { TRACE("Field3D::splitYupYdown"); +#if CHECK > 2 + if (yup_fields.size() != ydown_fields.size()) { + throw BoutException("Field3D::splitYupYdown: forward/backward parallel slices not in sync.\n" + " This is an internal library error"); + } +#endif + if (!yup_fields.empty()) { return; } @@ -148,6 +155,13 @@ void Field3D::splitYupYdown() { void Field3D::mergeYupYdown() { TRACE("Field3D::mergeYupYdown"); +#if CHECK > 2 + if (yup_fields.size() != ydown_fields.size()) { + throw BoutException("Field3D::mergeYupYdown: forward/backward parallel slices not in sync.\n" + " This is an internal library error"); + } +#endif + if (yup_fields.empty() && ydown_fields.empty()) { return; } @@ -157,12 +171,14 @@ void Field3D::mergeYupYdown() { } const Field3D& Field3D::ynext(int dir) const { +#if CHECK > 0 // Asked for more than yguards if (std::abs(dir) > fieldmesh->ystart) { throw BoutException( "Field3D: Call to ynext with %d which is more than number of yguards (%d)", dir, fieldmesh->ystart); } +#endif // ynext uses 1-indexing, but yup wants 0-indexing if (dir > 0) { From 682060459577f6e02e089415839739c3751b90b2 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 29 Jan 2019 11:31:55 +0000 Subject: [PATCH 29/34] Guard unit tests for Field3D::ynext --- tests/unit/field/test_field3d.cxx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index c6fe34f8c3..b9b9aac33b 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -284,7 +284,9 @@ TEST_F(Field3DTest, Ynext) { EXPECT_NE(&field, &ydown); EXPECT_NE(&yup, &ydown); +#if CHECK > 0 EXPECT_THROW(field.ynext(99), BoutException); +#endif } TEST_F(Field3DTest, ConstYnext) { @@ -300,7 +302,9 @@ TEST_F(Field3DTest, ConstYnext) { EXPECT_NE(&field2, &ydown); EXPECT_NE(&yup, &ydown); +#if CHECK > 0 EXPECT_THROW(field2.ynext(99), BoutException); +#endif } TEST_F(Field3DTest, GetGlobalMesh) { From a7b1bd6c4d9c7617dc37f69be0ad2472451d6821 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 4 Feb 2019 13:42:17 +0000 Subject: [PATCH 30/34] Add missing assert in test_zoidberg --- tools/pylib/zoidberg/test_zoidberg.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/pylib/zoidberg/test_zoidberg.py b/tools/pylib/zoidberg/test_zoidberg.py index 23117c8dad..45c1626381 100644 --- a/tools/pylib/zoidberg/test_zoidberg.py +++ b/tools/pylib/zoidberg/test_zoidberg.py @@ -34,7 +34,7 @@ def test_make_maps_slab(): assert var in maps # Each map should have the same shape as the grid - maps[var].shape == (nx, ny, nz) + assert maps[var].shape == (nx, ny, nz) # The first/last abs(offset) points are not valid, so ignore those interior_range = range(ny-abs(offset)) if offset > 0 else range(abs(offset), ny) From 6fb831965576cc91cc199cf97453c640b92deec2 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 4 Feb 2019 13:43:20 +0000 Subject: [PATCH 31/34] Don't store fields with placeholder values in shiftedmetric --- src/mesh/parallel/shiftedmetric.cxx | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx index cf6124bd29..d6171fc117 100644 --- a/src/mesh/parallel/shiftedmetric.cxx +++ b/src/mesh/parallel/shiftedmetric.cxx @@ -220,14 +220,9 @@ ShiftedMetric::shiftZ(const Field3D& f, for (auto& phase : phases) { // In C++17 std::vector::emplace_back returns a reference, which // would be very useful here! - - // FIXME: initialisation to -1 to avoid checkData choking on the - // uninitialised regions in assignment into the Field parallel - // slices in calcYUpDown - results.emplace_back(-1.0, &mesh); + results.emplace_back(&mesh); auto& current_result = results.back(); - // FIXME: uncomment the following after fixing Field3D::operator= - // current_result.allocate(); + current_result.allocate(); current_result.setLocation(f.getLocation()); for (int jx = 0; jx < mesh.LocalNx; jx++) { From 3e62a35942631ba2fa312fe4b8229c6ccf9b1a70 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 4 Feb 2019 17:23:03 +0000 Subject: [PATCH 32/34] Fix FCI MMS: g_22 should depend on the length of the field line Should be done better by actually calculating the arc length of the field line --- tests/MMS/spatial/fci/data/BOUT.inp | 2 +- tests/MMS/spatial/fci/fci_mms.cxx | 6 +++--- tests/MMS/spatial/fci/mms.py | 2 +- tests/MMS/spatial/fci/runtest | 2 +- tools/pylib/zoidberg/grid.py | 4 ++-- tools/pylib/zoidberg/zoidberg.py | 7 +++++++ 6 files changed, 15 insertions(+), 8 deletions(-) diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index 6ad5d66ca2..cf2208e3b3 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -2,7 +2,7 @@ grid = fci.grid.nc input = sin(y - 2*z) + sin(y - z) -solution = 6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z) +solution = (6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) MXG = 1 NXPE = 1 diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index 32b668f575..dc250d55b7 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -11,10 +11,10 @@ int main(int argc, char** argv) { // Communicate to calculate parallel transform mesh->communicate(input); - Field3D result{DDY(input)}; + Field3D result{Grad_par(input)}; Field3D error{result - solution}; - BoutReal l_2{sqrt(mean(SQ(error)))}; - BoutReal l_inf{max(abs(error), true)}; + BoutReal l_2{sqrt(mean(SQ(error), true, RGN_NOBNDRY))}; + BoutReal l_inf{max(abs(error), true, RGN_NOBNDRY)}; SAVE_ONCE6(input, solution, result, error, l_2, l_inf); diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py index 1154cfc9cc..477c605e84 100755 --- a/tests/MMS/spatial/fci/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -26,7 +26,7 @@ B = sqrt(Bpx**2 + Bt**2) def FCI_ddy(f): - return ( Bt * diff(f, y)*2.*pi/Ly + Bpx * diff(f, z)*2.*pi/Lz ) + return ( Bt * diff(f, y)*2.*pi/Ly + Bpx * diff(f, z)*2.*pi/Lz ) / B ############################################ # Equations solved diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 36eb96c7d7..99ec92edc2 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -61,7 +61,7 @@ for nslice in nslices: # Note that the Bz and Bzprime parameters here must be the same as in mms.py field = zb.field.Slab(Bz=0.05, Bzprime=0.1) # Create rectangular poloidal grids - poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid(nx, n, 1., 1.) + poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid(nx, n, 0.1, 1.) # Set the ylength and y locations ylength = 10. diff --git a/tools/pylib/zoidberg/grid.py b/tools/pylib/zoidberg/grid.py index 99e247d43f..8ca964517b 100644 --- a/tools/pylib/zoidberg/grid.py +++ b/tools/pylib/zoidberg/grid.py @@ -184,8 +184,8 @@ def metric(self): # Note: These y metrics are for Cartesian coordinates # If in cylindrical coordinates then these should be different - g_yy = 1.0 # Rmaj**2 - gyy = 1.0 # 1/Rmaj**2 + g_yy = np.ones(self.shape) + gyy = np.ones(self.shape) return {"dx":dx, "dy":dy3d, "dz": dz, "gyy": gyy, "g_yy":g_yy, diff --git a/tools/pylib/zoidberg/zoidberg.py b/tools/pylib/zoidberg/zoidberg.py index 5833d01134..7d11882fb2 100644 --- a/tools/pylib/zoidberg/zoidberg.py +++ b/tools/pylib/zoidberg/zoidberg.py @@ -208,6 +208,13 @@ def write_maps(grid, magnetic_field, maps, gridfile='fci.grid.nc', Bmag[:, yindex, :] = magnetic_field.Bmag(pol_grid.R, pol_grid.Z, ypos) pressure[:, yindex, :] = magnetic_field.pressure(pol_grid.R, pol_grid.Z, ypos) + metric["g_yy"][:, yindex, :] = (metric["g_yy"][:, yindex, :] + * (Bmag[:, yindex, :] + / magnetic_field.Byfunc(pol_grid.R, pol_grid.Z, ypos))**2) + metric["gyy"][:, yindex, :] = (metric["gyy"][:, yindex, :] + * (magnetic_field.Byfunc(pol_grid.R, pol_grid.Z, ypos) + / Bmag[:, yindex, :])**2) + # Get attributes from magnetic field (e.g. psi) attributes = {} for name in magnetic_field.attributes: From 35cdc728776c2fb98cc478ebc4cdc413e00e7162 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 4 Feb 2019 17:29:54 +0000 Subject: [PATCH 33/34] Remove some unnecessary code from test-smooth --- tests/integrated/test-smooth/test_smooth.cxx | 6 ------ 1 file changed, 6 deletions(-) diff --git a/tests/integrated/test-smooth/test_smooth.cxx b/tests/integrated/test-smooth/test_smooth.cxx index 15f7c656ec..c4670656c1 100644 --- a/tests/integrated/test-smooth/test_smooth.cxx +++ b/tests/integrated/test-smooth/test_smooth.cxx @@ -17,7 +17,6 @@ int main(int argc, char **argv) { Field2D input2d = f.create2D("1 + sin(2*y)"); Field3D input3d = f.create3D("gauss(x-0.5,0.2)*gauss(y-pi)*sin(3*y - z)"); - input3d.splitYupYdown(); mesh->getParallelTransform().calcYUpDown(input3d); SAVE_ONCE2(input2d, input3d); @@ -32,11 +31,6 @@ int main(int argc, char **argv) { // Output data dump.write(); - dump.close(); - - output << "\nFinished running test. Triggering error to quit\n\n"; - - MPI_Barrier(BoutComm::get()); // Wait for all processors to write data BoutFinalise(); return 0; From 69aea2521fa91024329a68ff86099bc050046194 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 5 Feb 2019 11:13:51 +0000 Subject: [PATCH 34/34] Fix: do y-interpolation in field-aligned space --- include/interpolation.hxx | 71 +++++++++++++++------------------------ 1 file changed, 28 insertions(+), 43 deletions(-) diff --git a/include/interpolation.hxx b/include/interpolation.hxx index 038285f824..2fe51191d9 100644 --- a/include/interpolation.hxx +++ b/include/interpolation.hxx @@ -118,53 +118,38 @@ const T interp_to(const T& var, CELL_LOC loc, REGION region = RGN_ALL) { // At least 2 boundary cells needed for interpolation in y-direction ASSERT0(fieldmesh->ystart >= 2); - if (var.hasYupYdown()) { - if ((location == CELL_CENTRE) && (loc == CELL_YLOW)) { // C2L - BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { - // Producing a stencil centred around a lower X value - result[i] = interp( - populateStencil(var, i)); - } - } else if (location == CELL_YLOW) { // L2C - BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { - // Stencil centred around a cell centre - result[i] = interp( - populateStencil(var, i)); - } + // We can't interpolate in y unless we're field-aligned + // FIXME: Add check once we label fields as orthogonal/aligned + + const T var_fa = fieldmesh->toFieldAligned(var); + if (region != RGN_NOBNDRY) { + // repeat the hack above for boundary points + // this avoids a duplicate toFieldAligned call if we had called + // result = toFieldAligned(result) + // to get the boundary cells + // + // result is requested in some boundary region(s) + result = var_fa; // NOTE: This is just for boundaries. FIX! + result.allocate(); + result.setLocation(loc); + } + + if ((location == CELL_CENTRE) && (loc == CELL_YLOW)) { // C2L + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + // Producing a stencil centred around a lower X value + result[i] = + interp(populateStencil(var_fa, i)); } - } else { - // var has no yup/ydown fields, so we need to shift into field-aligned - // coordinates - - const T var_fa = fieldmesh->toFieldAligned(var); - if (region != RGN_NOBNDRY) { - // repeat the hack above for boundary points - // this avoids a duplicate toFieldAligned call if we had called - // result = toFieldAligned(result) - // to get the boundary cells - // - // result is requested in some boundary region(s) - result = var_fa; // NOTE: This is just for boundaries. FIX! - result.allocate(); - result.setLocation(loc); + } else if (location == CELL_YLOW) { // L2C + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + // Stencil centred around a cell centre + result[i] = + interp(populateStencil(var_fa, i)); } + } - if ((location == CELL_CENTRE) && (loc == CELL_YLOW)) { // C2L - BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { - // Producing a stencil centred around a lower X value - result[i] = interp( - populateStencil(var_fa, i)); - } - } else if (location == CELL_YLOW) { // L2C - BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { - // Stencil centred around a cell centre - result[i] = interp( - populateStencil(var_fa, i)); - } - } + result = fieldmesh->fromFieldAligned(result); - result = fieldmesh->fromFieldAligned(result); - } break; } case CELL_ZLOW: {