From 132c14e897d63926953241f536a7ed7f8572cc7e Mon Sep 17 00:00:00 2001 From: John Omotani Date: Sun, 16 Dec 2018 23:07:37 +0000 Subject: [PATCH 1/5] Add hasBranchCutUp/Down methods to Mesh Useful for knowing where to correct y-periodic-boundary cells of zShift using ShiftAngle. --- include/bout/mesh.hxx | 8 ++++++++ src/mesh/impls/bout/boutmesh.hxx | 14 ++++++++++++++ 2 files changed, 22 insertions(+) diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index 432cf0fa22..0508de4cc3 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -322,7 +322,15 @@ class Mesh { /// \param[in] jx The local (on this processor) index in X /// \param[out] ts The Twist-Shift angle if periodic virtual bool periodicY(int jx, BoutReal &ts) const = 0; + + /// Is there a branch cut at the lower y edge of this processor's grid points + /// at x-index jx? + virtual bool hasBranchCutDown(int jx) const = 0; + /// Is there a branch cut at the upper y edge of this processor's grid points + /// at x-index jx? + virtual bool hasBranchCutUp(int jx) const = 0; + virtual int ySize(int jx) const; ///< The number of points in Y at fixed X index \p jx // Y communications diff --git a/src/mesh/impls/bout/boutmesh.hxx b/src/mesh/impls/bout/boutmesh.hxx index 6964ebf2dc..c8babd1226 100644 --- a/src/mesh/impls/bout/boutmesh.hxx +++ b/src/mesh/impls/bout/boutmesh.hxx @@ -105,6 +105,20 @@ class BoutMesh : public Mesh { /// \param[in] jx The local (on this processor) index in X bool periodicY(int jx) const; + /// Is there a branch cut at the lower y edge of this processor's grid points + /// at x-index jx? + bool hasBranchCutDown(int jx) const { + return (TS_down_in && (DDATA_INDEST != -1) && (jx < DDATA_XSPLIT)) + || (TS_down_out && (DDATA_OUTDEST != -1) && (jx >= DDATA_XSPLIT)); + } + + /// Is there a branch cut at the upper y edge of this processor's grid points + /// at x-index jx? + bool hasBranchCutUp(int jx) const { + return (TS_up_in && (UDATA_INDEST != -1) && (jx < UDATA_XSPLIT)) + || (TS_up_out && (UDATA_OUTDEST != -1) && (jx >= UDATA_XSPLIT)); + } + int ySize(int jx) const; ///< The number of points in Y at fixed X index \p jx ///////////////////////////////////////////// From 68a27d89f2361ec0b9c65ec56e73677636975b44 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 17 Dec 2018 15:09:45 +0000 Subject: [PATCH 2/5] Add hasBranchCutDown() and hasBranchCutUp() methods to FakeMesh --- tests/unit/test_extras.hxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/unit/test_extras.hxx b/tests/unit/test_extras.hxx index b042e2a933..a36b842280 100644 --- a/tests/unit/test_extras.hxx +++ b/tests/unit/test_extras.hxx @@ -104,6 +104,8 @@ public: MPI_Comm getYcomm(int UNUSED(jx)) const { return MPI_COMM_NULL; } bool periodicY(int UNUSED(jx)) const { return true; } bool periodicY(int UNUSED(jx), BoutReal &UNUSED(ts)) const { return true; } + bool hasBranchCutDown(int UNUSED(jx)) const { return false; }; + bool hasBranchCutUp(int UNUSED(jx)) const { return false; }; bool firstY() const { return true; } bool lastY() const { return true; } bool firstY(int UNUSED(xpos)) const { return true; } From 46a36a90be9cd672f4263e49c6354e9bd56265ed Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 11 Dec 2018 15:59:46 +0000 Subject: [PATCH 3/5] Add ShiftToFieldAligned implementation of ParallelTransform Similar to ShiftedMetric, but transforms to globally field-aligned coordinates instead of calculating yup/ydown fields. --- include/bout/paralleltransform.hxx | 6 +- src/mesh/mesh.cxx | 7 +- src/mesh/parallel/makefile | 2 +- src/mesh/parallel/shifttofieldaligned.cxx | 183 ++++++++++++++++++++++ src/mesh/parallel/shifttofieldaligned.hxx | 137 ++++++++++++++++ 5 files changed, 330 insertions(+), 5 deletions(-) create mode 100644 src/mesh/parallel/shifttofieldaligned.cxx create mode 100644 src/mesh/parallel/shifttofieldaligned.hxx diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx index f2f40ed62f..f313c978c3 100644 --- a/include/bout/paralleltransform.hxx +++ b/include/bout/paralleltransform.hxx @@ -173,17 +173,17 @@ private: * @param[in] zangle The angle (z coordinate) to shift by * @param[out] out A 1D array of length \p len, already allocated */ - void shiftZ(const BoutReal *in, int len, BoutReal zangle, BoutReal *out); + void shiftZ(const BoutReal* in, int len, BoutReal zangle, BoutReal* out); /*! * Shift a given 1D array, assumed to be in Z, by the given \p zangle * * @param[in] in A 1D array of length mesh.LocalNz - * @param[in] phs Phase shift, assumed to have length (mesh.LocalNz/2 + 1) i.e. the number of modes + * @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 dcomplex* phs, BoutReal* out); }; - #endif // __PARALLELTRANSFORM_H__ diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx index f645b89f0c..ea4147fa3a 100644 --- a/src/mesh/mesh.cxx +++ b/src/mesh/mesh.cxx @@ -13,6 +13,7 @@ #include #include "parallel/fci.hxx" +#include "parallel/shifttofieldaligned.hxx" Mesh* Mesh::create(GridDataSource *s, Options *opt) { return MeshFactory::getInstance()->createMesh(s, opt); @@ -299,8 +300,12 @@ void Mesh::setParallelTransform() { }else if(ptstr == "shifted") { // Shifted metric method - transform = bout::utils::make_unique(*this); + transform = bout::utils::make_unique(*this); + }else if(ptstr == "shifttofieldaligned") { + // Shifted metric method, shifting to globally field-aligned coordinates + transform = bout::utils::make_unique(*this); + }else if(ptstr == "fci") { Options *fci_options = Options::getRoot()->getSection("fci"); diff --git a/src/mesh/parallel/makefile b/src/mesh/parallel/makefile index 9c4817b008..97728dd3c7 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 shifttofieldaligned.cxx TARGET = lib include $(BOUT_TOP)/make.config diff --git a/src/mesh/parallel/shifttofieldaligned.cxx b/src/mesh/parallel/shifttofieldaligned.cxx new file mode 100644 index 0000000000..a3b84f774f --- /dev/null +++ b/src/mesh/parallel/shifttofieldaligned.cxx @@ -0,0 +1,183 @@ +/* + * Implements the shifted metric method for parallel derivatives by + * transforming to globally field-aligned coordinates (in contrast to + * ShiftedMetric where the coordinates are only locally field-aligned). + * + * By default fields are stored so that X-Z are orthogonal, + * and so not aligned in Y. + * + */ + +#include "shifttofieldaligned.hxx" + +#include +#include +#include +#include + +namespace { +void fixzShiftBoundaries(Field2D& zShift) { + Mesh* localmesh = zShift.getMesh(); + + localmesh->communicate(zShift); + + for (int x = 0; x < localmesh->LocalNx; x++) { + if (localmesh->hasBranchCutDown(x)) { + // ShiftAngle is the total toroidal shift of a field line after one poloidal turn + BoutReal ShiftAngle = 0.; + localmesh->periodicY(x, ShiftAngle); + for (int y = 0; y < localmesh->ystart; y++) { + zShift(x, y) -= ShiftAngle; + } + } + if (localmesh->hasBranchCutUp(x)) { + // ShiftAngle is the total toroidal shift of a field line after one poloidal turn + BoutReal ShiftAngle = 0.; + localmesh->periodicY(x, ShiftAngle); + for (int y = localmesh->yend + 1; y < localmesh->LocalNy; y++) { + zShift(x, y) += ShiftAngle; + } + } + } + + zShift.applyBoundary("free_o3"); +} +} // namespace + +ShiftToFieldAligned::ShiftToFieldAligned(Mesh& mesh_in) { + // Must *not* twist-shift when communicating + bool twistshift = Options::root()["TwistShift"].withDefault(false); + if (twistshift) { + throw BoutException("ShiftToFieldAligned requires TwistShift=false"); + } + + Field2D zShift(&mesh_in); + // Read the zShift angle from the mesh + if (mesh_in.get(zShift, "zShift")) { + // No zShift variable. Try qinty in BOUT grid files + mesh_in.get(zShift, "qinty"); + } + + init(mesh_in, zShift); +} + +void ShiftToFieldAligned::init(Mesh& mesh_in, const Field2D& zShift) { + // always include CELL_CENTRE implementation + implementations.emplace(CELL_CENTRE, Implementation(mesh_in, CELL_CENTRE, zShift)); + + if (mesh_in.StaggerGrids) { + // Don't know which locations we will need, so populate all those that can + // be interpolated (for now at least) + implementations.emplace(CELL_XLOW, Implementation(mesh_in, CELL_XLOW, zShift)); + implementations.emplace(CELL_XLOW, Implementation(mesh_in, CELL_YLOW, zShift)); + implementations.emplace(CELL_XLOW, Implementation(mesh_in, CELL_ZLOW, zShift)); + } +} + +ShiftToFieldAligned::Implementation::Implementation(Mesh& mesh_in, + const CELL_LOC location_in, + const Field2D& zShift_in) + : mesh(mesh_in), location(location_in), zShift(&mesh_in) { + + zShift = interp_to(zShift_in, location); + fixzShiftBoundaries(zShift); + cachePhases(); +} + +void ShiftToFieldAligned::Implementation::cachePhases() { + // If we wanted to be efficient we could move the following cached phase setup + // into the relevant shifting routines (with a bool to record if the phases + // have already been cached) so that we only calculate the phase if we + // actually call a relevant shift routine -- however that would make calls to + // toFieldAligned and fromFieldAligned non-const, so it's better to keep this + // in the constructor + + // As we're attached to a mesh we can expect the z direction to + // not change once we've been created so precalculate the complex + // phases used in transformations + int nmodes = mesh.LocalNz / 2 + 1; + BoutReal zlength = mesh.getCoordinates()->zlength(); + + // Allocate storage for complex intermediate + cmplx = Array(nmodes); + std::fill(cmplx.begin(), cmplx.end(), 0.0); + + // Allocate storage for our 3d vector structures. + fromAlignedPhs = Matrix>(mesh.LocalNx, mesh.LocalNy); + for (auto& element : fromAlignedPhs) { + element = Array(mesh.LocalNz); + } + toAlignedPhs = Matrix>(mesh.LocalNx, mesh.LocalNy); + for (auto& element : toAlignedPhs) { + element = Array(mesh.LocalNz); + } + + // To/From field aligned phases + for (int jx = 0; jx < mesh.LocalNx; jx++) { + for (int jy = 0; jy < mesh.LocalNy; jy++) { + for (int jz = 0; jz < nmodes; jz++) { + BoutReal kwave = jz * 2.0 * PI / zlength; // wave number is 1/[rad] + fromAlignedPhs(jx, jy)[jz] = + dcomplex(cos(kwave * zShift(jx, jy)), -sin(kwave * zShift(jx, jy))); + toAlignedPhs(jx, jy)[jz] = + dcomplex(cos(kwave * zShift(jx, jy)), sin(kwave * zShift(jx, jy))); + } + } + } +} + +/*! + * Get the shifted field so that X-Z is not orthogonal, + * and Y is then field aligned. + */ +const Field3D ShiftToFieldAligned::Implementation::toFieldAligned(const Field3D& f, + const REGION region) { + return shiftZ(f, toAlignedPhs, region); +} + +/*! + * Shift back, so that X-Z is orthogonal, + * but Y is not field aligned. + */ +const Field3D ShiftToFieldAligned::Implementation::fromFieldAligned(const Field3D& f, + const REGION region) { + return shiftZ(f, fromAlignedPhs, region); +} + +const Field3D ShiftToFieldAligned::Implementation::shiftZ( + const Field3D& f, const Matrix>& phs, const REGION region) { + ASSERT1(&mesh == f.getMesh()); + ASSERT1(f.getLocation() == location); + + if (mesh.LocalNz == 1) + return f; // Shifting makes no difference + + Field3D result(&mesh); + result.allocate(); + result.setLocation(location); + + for (const auto& i : mesh.getRegion2D(REGION_STRING(region))) { + shiftZ(f(i.x(), i.y()), phs(i.x(), i.y()), result(i.x(), i.y())); + } + + return result; +} + +void ShiftToFieldAligned::Implementation::shiftZ(const BoutReal* in, + const Array& phs, + BoutReal* out) { + // Take forward FFT + rfft(in, mesh.LocalNz, cmplx.begin()); + + // Following is an algorithm approach to write a = a*b where a and b are + // vectors of dcomplex. + // std::transform(cmplxOneOff.begin(),cmplxOneOff.end(), ptr.begin(), + // cmplxOneOff.begin(), std::multiplies()); + + const int nmodes = cmplx.size(); + for (int jz = 1; jz < nmodes; jz++) { + cmplx[jz] *= phs[jz]; + } + + irfft(cmplx.begin(), mesh.LocalNz, out); // Reverse FFT +} diff --git a/src/mesh/parallel/shifttofieldaligned.hxx b/src/mesh/parallel/shifttofieldaligned.hxx new file mode 100644 index 0000000000..ea963f7868 --- /dev/null +++ b/src/mesh/parallel/shifttofieldaligned.hxx @@ -0,0 +1,137 @@ +/************************************************************************** + * Parallel derivatives use shifted-metric scheme implemented by shifting + * to globally field-aligned coordinates + * + ************************************************************************** + * Copyright 2018 B.D.Dudson, P. Hill, J. Omotani + * + * Contact: Ben Dudson, bd512@york.ac.uk + * + * This file is part of BOUT++. + * + * BOUT++ is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * BOUT++ is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with BOUT++. If not, see . + * + **************************************************************************/ + +#ifndef __SHIFTTOFIELDALIGNED_H__ +#define __SHIFTTOFIELDALIGNED_H__ + +#include + +class ShiftToFieldAligned : public ParallelTransform { +public: + ShiftToFieldAligned(Mesh& mesh_in); + /// Use an existing zShift + ShiftToFieldAligned(Mesh& mesh_in, const Field2D& zShift) { init(mesh_in, zShift); } + + /*! + * Do nothing here since we don't cache field-aligned fields + */ + void calcYUpDown(Field3D& UNUSED(f)) override {} + + /*! + * Uses FFTs and a phase shift to align the grid points + * with the y coordinate (along magnetic field usually). + * + * Note that the returned field will no longer be orthogonal + * in X-Z, and the metric tensor will need to be changed + * if X derivatives are used. + */ + const Field3D toFieldAligned(const Field3D& f, const REGION region = RGN_NOX) override { + return implementations.at(f.getLocation()).toFieldAligned(f, region); + } + + /*! + * Converts a field back to X-Z orthogonal coordinates + * from field aligned coordinates. + */ + const Field3D fromFieldAligned(const Field3D& f, + const REGION region = RGN_NOX) override { + return implementations.at(f.getLocation()).fromFieldAligned(f, region); + } + + bool canToFromFieldAligned() override { return true; } + +private: + class Implementation { + public: + Implementation(Mesh& mesh_in, const CELL_LOC location_in, const Field2D& zShift_in); + + /*! + * Uses FFTs and a phase shift to align the grid points + * with the y coordinate (along magnetic field usually). + * + * Note that the returned field will no longer be orthogonal + * in X-Z, and the metric tensor will need to be changed + * if X derivatives are used. + */ + const Field3D toFieldAligned(const Field3D& f, const REGION region = RGN_NOX); + + /*! + * Converts a field back to X-Z orthogonal coordinates + * from field aligned coordinates. + */ + const Field3D fromFieldAligned(const Field3D& f, const REGION region = RGN_NOX); + + private: + Mesh& mesh; + const CELL_LOC location; + Field2D zShift; + + /// Calculate and store the phases for to/from field aligned using zShift + void cachePhases(); + + Array cmplx; ///< A temporary array, used for input/output to fft routines + Array + cmplxLoc; ///< A temporary array, used for input/output to fft routines + + Matrix> toAlignedPhs; ///< Cache of phase shifts for + /// transforming from X-Z orthogonal + /// coordinates to field-aligned + /// coordinates. + Matrix> fromAlignedPhs; ///< Cache of phase shifts for + /// transforming from field-aligned + /// coordinates to X-Z orthogonal + /// coordinates. + + /*! + * Shift a 3D field \p f by the given phase \p phs in Z + * + * Calculates FFT in Z, multiplies by the complex phase + * and inverse FFTS. + * + * @param[in] f The field to shift + * @param[in] phs The phase to shift by + */ + const Field3D shiftZ(const Field3D& f, const Matrix>& phs, + const REGION region = RGN_NOX); + + /*! + * Shift a given 1D array, assumed to be in Z, by the given \p zangle + * + * @param[in] in A 1D array of length mesh.LocalNz + * @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 Array& phs, BoutReal* out); + }; + + std::map implementations; + + // populate the implementations map + void init(Mesh& mesh_in, const Field2D& zShift); +}; + +#endif // __SHIFTTOFIELDALIGNED_H__ From f3038b9512756deaaf53b1616163129a54c9685d Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 11 Dec 2018 17:47:12 +0000 Subject: [PATCH 4/5] Basic test comparing ShiftToFieldAligned with methods from ShiftedMetric --- .../data/BOUT.inp | 12 +++++ .../test-yupdown-shifttofieldaligned/makefile | 6 +++ .../test-yupdown-shifttofieldaligned/runtest | 33 ++++++++++++ .../test_yupdown.cxx | 51 +++++++++++++++++++ .../integrated/test-yupdown/test_yupdown.cxx | 34 ++++++------- 5 files changed, 119 insertions(+), 17 deletions(-) create mode 100644 tests/integrated/test-yupdown-shifttofieldaligned/data/BOUT.inp create mode 100644 tests/integrated/test-yupdown-shifttofieldaligned/makefile create mode 100755 tests/integrated/test-yupdown-shifttofieldaligned/runtest create mode 100644 tests/integrated/test-yupdown-shifttofieldaligned/test_yupdown.cxx diff --git a/tests/integrated/test-yupdown-shifttofieldaligned/data/BOUT.inp b/tests/integrated/test-yupdown-shifttofieldaligned/data/BOUT.inp new file mode 100644 index 0000000000..fdc5cae5d0 --- /dev/null +++ b/tests/integrated/test-yupdown-shifttofieldaligned/data/BOUT.inp @@ -0,0 +1,12 @@ +ShiftWithoutTwist = true + +[mesh] +nx = 12 +ny = 12 + +zShift = y + +var = sin(y - z) + +paralleltransform = shifttofieldaligned + diff --git a/tests/integrated/test-yupdown-shifttofieldaligned/makefile b/tests/integrated/test-yupdown-shifttofieldaligned/makefile new file mode 100644 index 0000000000..26728c595f --- /dev/null +++ b/tests/integrated/test-yupdown-shifttofieldaligned/makefile @@ -0,0 +1,6 @@ + +BOUT_TOP = ../../.. + +SOURCEC = test_yupdown.cxx + +include $(BOUT_TOP)/make.config diff --git a/tests/integrated/test-yupdown-shifttofieldaligned/runtest b/tests/integrated/test-yupdown-shifttofieldaligned/runtest new file mode 100755 index 0000000000..881beb6eb2 --- /dev/null +++ b/tests/integrated/test-yupdown-shifttofieldaligned/runtest @@ -0,0 +1,33 @@ +#!/usr/bin/env python3 + +from boututils.run_wrapper import shell, shell_safe, launch_safe, getmpirun +from boutdata.collect import collect +from sys import stdout, exit + +from numpy import max, abs + +MPIRUN=getmpirun() + +shell_safe("make > make.log") + + +s, out = launch_safe("./test_yupdown", runcmd=MPIRUN, nproc=1, pipe=True, verbose=True) + +with open("run.log", "w") as f: + f.write(out) + +vars = [ ("ddy", "ddy2") ] +for v1, v2 in vars: + stdout.write("Testing %s and %s ... " % (v1, v2) ) + ddy = collect(v1, path="data", xguards=False, yguards=False, info=False) + ddy2 = collect(v2, path="data", xguards=False, yguards=False, info=False) + + diff = max(abs(ddy - ddy2)) + + if diff < 1e-8: + print("Passed (Max difference %e)" % (diff)) + else: + print("Failed (Max difference %e)" % (diff)) + exit(1) + +exit(0) diff --git a/tests/integrated/test-yupdown-shifttofieldaligned/test_yupdown.cxx b/tests/integrated/test-yupdown-shifttofieldaligned/test_yupdown.cxx new file mode 100644 index 0000000000..f7cef52403 --- /dev/null +++ b/tests/integrated/test-yupdown-shifttofieldaligned/test_yupdown.cxx @@ -0,0 +1,51 @@ +#include + +#include + +#include + +// Y derivative assuming field is aligned in Y +const Field3D DDY_aligned(const Field3D &f) { + Field3D result; + result.allocate(); + result = 0.0; + + for(int i=0;iLocalNx;i++) + for(int j=mesh->ystart;j<=mesh->yend;j++) + for(int k=0;kLocalNz;k++) + result(i,j,k) = 0.5*(f(i,j+1,k) - f(i,j-1,k)); + + return result; +} + +int main(int argc, char** argv) { + + BoutInitialise(argc, argv); + + ShiftedMetric s(*mesh); + + // Read variable from mesh + Field3D var; + mesh->get(var, "var"); + + // Var starts in orthogonal X-Z coordinates + + // Calculate d/dy using ShiftToFieldAligned + Field3D ddy = DDY(var); + + // Change into field-aligned coordinates + Field3D var_aligned = s.toFieldAligned(var); + + // var now field aligned + Field3D ddy2 = DDY_aligned(var_aligned); + + // Shift back to orthogonal X-Z coordinates + ddy2 = s.fromFieldAligned(ddy2); + + SAVE_ONCE2(ddy, ddy2); + dump.write(); + + BoutFinalise(); + + return 0; +} diff --git a/tests/integrated/test-yupdown/test_yupdown.cxx b/tests/integrated/test-yupdown/test_yupdown.cxx index dbadc6372f..e05f5db587 100644 --- a/tests/integrated/test-yupdown/test_yupdown.cxx +++ b/tests/integrated/test-yupdown/test_yupdown.cxx @@ -5,31 +5,31 @@ #include // Y derivative using yup() and ydown() fields -const Field3D DDY_yud(const Field3D &f) { +const Field3D DDY_yud(const Field3D& f) { Field3D result; result.allocate(); result = 0.0; - - for(int i=0;iLocalNx;i++) - for(int j=mesh->ystart;j<=mesh->yend;j++) - for(int k=0;kLocalNz;k++) - result(i,j,k) = 0.5*(f.yup()(i,j+1,k) - f.ydown()(i,j-1,k)); + + for (int i = 0; i < mesh->LocalNx; i++) + for (int j = mesh->ystart; j <= mesh->yend; j++) + for (int k = 0; k < mesh->LocalNz; k++) + result(i, j, k) = 0.5 * (f.yup()(i, j + 1, k) - f.ydown()(i, j - 1, k)); return result; } // Y derivative assuming field is aligned in Y -const Field3D DDY_aligned(const Field3D &f) { +const Field3D DDY_aligned(const Field3D& f) { Field3D result; result.allocate(); result = 0.0; - - for(int i=0;iLocalNx;i++) - for(int j=mesh->ystart;j<=mesh->yend;j++) - for(int k=0;kLocalNz;k++) - result(i,j,k) = 0.5*(f(i,j+1,k) - f(i,j-1,k)); - + + for (int i = 0; i < mesh->LocalNx; i++) + for (int j = mesh->ystart; j <= mesh->yend; j++) + for (int k = 0; k < mesh->LocalNz; k++) + result(i, j, k) = 0.5 * (f(i, j + 1, k) - f(i, j - 1, k)); + return result; } @@ -49,7 +49,7 @@ int main(int argc, char** argv) { // Calculate yup and ydown s.calcYUpDown(var); - + // Calculate d/dy using yup() and ydown() fields Field3D ddy = DDY(var); @@ -59,13 +59,13 @@ int main(int argc, char** argv) { // Change into field-aligned coordinates Field3D var_aligned = mesh->toFieldAligned(var); - + // var now field aligned Field3D ddy_check = DDY_aligned(var_aligned); - + // Shift back to orthogonal X-Z coordinates ddy_check = mesh->fromFieldAligned(ddy_check); - + SAVE_ONCE3(ddy, ddy2, ddy_check); dump.write(); From 7a73f8ce7dca5dcd56d7a0a81fe23a4edcaf5535 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 17 Dec 2018 15:28:08 +0000 Subject: [PATCH 5/5] Use Tensor and avoid race conditions in ShiftToFieldAligned Copies changes made in #1462 to ShiftedMetric --- src/mesh/parallel/shifttofieldaligned.cxx | 47 +++++++++-------------- src/mesh/parallel/shifttofieldaligned.hxx | 24 +++++------- 2 files changed, 28 insertions(+), 43 deletions(-) diff --git a/src/mesh/parallel/shifttofieldaligned.cxx b/src/mesh/parallel/shifttofieldaligned.cxx index a3b84f774f..9b760fd929 100644 --- a/src/mesh/parallel/shifttofieldaligned.cxx +++ b/src/mesh/parallel/shifttofieldaligned.cxx @@ -95,33 +95,21 @@ void ShiftToFieldAligned::Implementation::cachePhases() { // As we're attached to a mesh we can expect the z direction to // not change once we've been created so precalculate the complex // phases used in transformations - int nmodes = mesh.LocalNz / 2 + 1; + nmodes = mesh.LocalNz / 2 + 1; BoutReal zlength = mesh.getCoordinates()->zlength(); - // Allocate storage for complex intermediate - cmplx = Array(nmodes); - std::fill(cmplx.begin(), cmplx.end(), 0.0); - // Allocate storage for our 3d vector structures. - fromAlignedPhs = Matrix>(mesh.LocalNx, mesh.LocalNy); - for (auto& element : fromAlignedPhs) { - element = Array(mesh.LocalNz); - } - toAlignedPhs = Matrix>(mesh.LocalNx, mesh.LocalNy); - for (auto& element : toAlignedPhs) { - element = Array(mesh.LocalNz); - } + fromAlignedPhs = Tensor(mesh.LocalNx, mesh.LocalNy, mesh.LocalNz); + toAlignedPhs = Tensor(mesh.LocalNx, mesh.LocalNy, mesh.LocalNz); // To/From field aligned phases - for (int jx = 0; jx < mesh.LocalNx; jx++) { - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = 0; jz < nmodes; jz++) { - BoutReal kwave = jz * 2.0 * PI / zlength; // wave number is 1/[rad] - fromAlignedPhs(jx, jy)[jz] = - dcomplex(cos(kwave * zShift(jx, jy)), -sin(kwave * zShift(jx, jy))); - toAlignedPhs(jx, jy)[jz] = - dcomplex(cos(kwave * zShift(jx, jy)), sin(kwave * zShift(jx, jy))); - } + BOUT_FOR(i, mesh.getRegion2D("RGN_ALL")) { + for (int jz = 0; jz < nmodes; jz++) { + BoutReal kwave = jz * 2.0 * PI / zlength; // wave number is 1/[rad] + fromAlignedPhs(i.x(), i.y(), jz) = + dcomplex(cos(kwave * zShift[i]), -sin(kwave * zShift[i])); + toAlignedPhs(i.x(), i.y(), jz) = + dcomplex(cos(kwave * zShift[i]), sin(kwave * zShift[i])); } } } @@ -145,7 +133,7 @@ const Field3D ShiftToFieldAligned::Implementation::fromFieldAligned(const Field3 } const Field3D ShiftToFieldAligned::Implementation::shiftZ( - const Field3D& f, const Matrix>& phs, const REGION region) { + const Field3D& f, const Tensor& phs, const REGION region) { ASSERT1(&mesh == f.getMesh()); ASSERT1(f.getLocation() == location); @@ -156,28 +144,29 @@ const Field3D ShiftToFieldAligned::Implementation::shiftZ( result.allocate(); result.setLocation(location); - for (const auto& i : mesh.getRegion2D(REGION_STRING(region))) { - shiftZ(f(i.x(), i.y()), phs(i.x(), i.y()), result(i.x(), i.y())); + BOUT_FOR(i, mesh.getRegion2D(REGION_STRING(region))) { + shiftZ(&f(i, 0), &phs(i.x(), i.y(), 0), &result(i, 0)); } return result; } void ShiftToFieldAligned::Implementation::shiftZ(const BoutReal* in, - const Array& phs, + const dcomplex* phs, BoutReal* out) { + Array cmplx(nmodes); + // Take forward FFT - rfft(in, mesh.LocalNz, cmplx.begin()); + rfft(in, mesh.LocalNz, &cmplx[0]); // Following is an algorithm approach to write a = a*b where a and b are // vectors of dcomplex. // std::transform(cmplxOneOff.begin(),cmplxOneOff.end(), ptr.begin(), // cmplxOneOff.begin(), std::multiplies()); - const int nmodes = cmplx.size(); for (int jz = 1; jz < nmodes; jz++) { cmplx[jz] *= phs[jz]; } - irfft(cmplx.begin(), mesh.LocalNz, out); // Reverse FFT + irfft(&cmplx[0], mesh.LocalNz, out); // Reverse FFT } diff --git a/src/mesh/parallel/shifttofieldaligned.hxx b/src/mesh/parallel/shifttofieldaligned.hxx index ea963f7868..57b8fcfc6c 100644 --- a/src/mesh/parallel/shifttofieldaligned.hxx +++ b/src/mesh/parallel/shifttofieldaligned.hxx @@ -92,18 +92,14 @@ private: /// Calculate and store the phases for to/from field aligned using zShift void cachePhases(); - Array cmplx; ///< A temporary array, used for input/output to fft routines - Array - cmplxLoc; ///< A temporary array, used for input/output to fft routines - - Matrix> toAlignedPhs; ///< Cache of phase shifts for - /// transforming from X-Z orthogonal - /// coordinates to field-aligned - /// coordinates. - Matrix> fromAlignedPhs; ///< Cache of phase shifts for - /// transforming from field-aligned - /// coordinates to X-Z orthogonal - /// coordinates. + int nmodes; + + Tensor toAlignedPhs; ///< Cache of phase shifts for transforming from X-Z + /// orthogonal coordinates to field-aligned + /// coordinates. + Tensor fromAlignedPhs; ///< Cache of phase shifts for transforming from + /// field-aligned coordinates to X-Z orthogonal + /// coordinates. /*! * Shift a 3D field \p f by the given phase \p phs in Z @@ -114,7 +110,7 @@ private: * @param[in] f The field to shift * @param[in] phs The phase to shift by */ - const Field3D shiftZ(const Field3D& f, const Matrix>& phs, + const Field3D shiftZ(const Field3D& f, const Tensor& phs, const REGION region = RGN_NOX); /*! @@ -125,7 +121,7 @@ private: * i.e. the number of modes * @param[out] out A 1D array of length mesh.LocalNz, already allocated */ - void shiftZ(const BoutReal* in, const Array& phs, BoutReal* out); + void shiftZ(const BoutReal* in, const dcomplex* phs, BoutReal* out); }; std::map implementations;