diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index 2d48107412..cce4136b3b 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -39,6 +39,8 @@ #include #include +#include + class Mesh; /*! @@ -249,6 +251,8 @@ private: void checkCovariant(); // check that contravariant tensors are positive (if expected) and finite (always) void checkContravariant(); + /// Read quantities with given suffix from `Mesh` + void readFromMesh(Options* options, const std::string& suffix); }; #endif // BOUT_COORDINATES_H diff --git a/include/bout/griddata.hxx b/include/bout/griddata.hxx index 29a32e5779..fd913a27f1 100644 --- a/include/bout/griddata.hxx +++ b/include/bout/griddata.hxx @@ -28,12 +28,15 @@ class GridDataSource; #ifndef BOUT_GRIDDATA_H #define BOUT_GRIDDATA_H -#include "mesh.hxx" -#include "bout/bout_types.hxx" -#include "bout/options.hxx" - +#include #include #include +#include +#include +#include + +#include +#include /// Interface class to serve grid data /*! @@ -46,8 +49,8 @@ public: GridDataSource(const bool source_is_file = false) : is_file(source_is_file) {} virtual ~GridDataSource() = default; - virtual bool - hasVar(const std::string& name) = 0; ///< Test if source can supply a variable + /// Test if source can supply a variable + virtual bool hasVar(const std::string& name) const = 0; /// Get a string virtual bool get(Mesh* m, std::string& sval, const std::string& name, @@ -95,9 +98,9 @@ class GridFile : public GridDataSource { public: GridFile() = delete; GridFile(std::string gridfilename); - ~GridFile() = default; + ~GridFile() override = default; - bool hasVar(const std::string& name) override; + bool hasVar(const std::string& name) const override; /// Get a string bool get(Mesh* m, std::string& sval, const std::string& name, @@ -178,7 +181,7 @@ public: /*! * Checks if the options has a given variable */ - bool hasVar(const std::string& name) override; + bool hasVar(const std::string& name) const override; /*! * Reads strings from options. Uses Options::get to handle diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index 25e66fe5d2..17d3edbe31 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -479,7 +479,7 @@ public: // Boundary regions /// Return a vector containing all the boundary regions on this processor - virtual std::vector getBoundaries() = 0; + virtual std::vector getBoundaries() const = 0; /// Get the set of all possible boundaries in this configuration virtual std::set getPossibleBoundaries() const { return {}; } @@ -496,7 +496,7 @@ public: /// get only xout: /// mesh->getBoundariesPar(Mesh::BoundaryParType::xout) virtual std::vector> - getBoundariesPar(BoundaryParType type = BoundaryParType::all) = 0; + getBoundariesPar(BoundaryParType type = BoundaryParType::all) const = 0; /// Add a parallel(Y) boundary to this processor virtual void addBoundaryPar(std::shared_ptr UNUSED(bndry), diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index d2634d3f88..fb239519bc 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -5,66 +5,75 @@ **************************************************************************/ #include +#include +#include +#include #include #include #include -#include -#include -#include - +#include #include #include +#include +#include +#include #include +#include +#include +#include #include +#include +#include +#include +#include +#include +#include +#include -#include +#include +#include +#include +#include +#include +#include +#include + +#include #include "invert3x3.hxx" #include "parallel/fci.hxx" #include "parallel/shiftedmetricinterp.hxx" namespace { -/// Interpolate a Field2D to a new CELL_LOC with interp_to. -/// Communicates to set internal guard cells. -/// Boundary guard cells are set by extrapolating from the grid, like -/// 'free_o3' boundary conditions -/// Corner guard cells are set to BoutNaN -Field2D interpolateAndExtrapolate(const Field2D& f, CELL_LOC location, bool extrapolate_x, - bool extrapolate_y, bool no_extra_interpolate, - ParallelTransform* UNUSED(pt) = nullptr) { +// Extrapolate into boundaries (if requested) so that differential geometry +// terms can be interpolated if necessary +// Note: cannot use applyBoundary("free_o3") here because applyBoundary() +// would try to create a new Coordinates object since we have not finished +// initializing yet, leading to an infinite recursion. +// Also, here we interpolate for the boundary points at xstart/ystart and +// (xend+1)/(yend+1) instead of extrapolating. +template > +void fillGuards_impl(T& result, CELL_LOC location, const T& f, bool extrapolate_x, + bool extrapolate_y, bool no_extra_interpolate = false) { + const auto* localmesh = result.getMesh(); + + const auto zstart = bout::build::use_metric_3d ? localmesh->zstart : 0; + const auto zend = bout::build::use_metric_3d ? localmesh->zend : 0; - Mesh* localmesh = f.getMesh(); - Field2D result = interp_to(f, location, "RGN_NOBNDRY"); - // Ensure result's data is unique. Otherwise result might be a duplicate of - // f (if no interpolation is needed, e.g. if interpolation is in the - // z-direction); then f would be communicated. Since this function is used - // on geometrical quantities that might not be periodic in y even on closed - // field lines (due to dependence on integrated shear), we don't want to - // communicate f. We will sort out result's boundary guard cells below, but - // not f's so we don't want to change f. - result.allocate(); - localmesh->communicate_no_slices(result); - - // Extrapolate into boundaries (if requested) so that differential geometry - // terms can be interpolated if necessary - // Note: cannot use applyBoundary("free_o3") here because applyBoundary() - // would try to create a new Coordinates object since we have not finished - // initializing yet, leading to an infinite recursion. - // Also, here we interpolate for the boundary points at xstart/ystart and - // (xend+1)/(yend+1) instead of extrapolating. for (auto& bndry : localmesh->getBoundaries()) { if ((extrapolate_x and bndry->bx != 0) or (extrapolate_y and bndry->by != 0)) { - int extrap_start = 0; - if (not no_extra_interpolate) { - // Can use no_extra_interpolate argument to skip the extra interpolation when we - // want to extrapolate the Christoffel symbol terms which come from derivatives so - // don't have the extra point set already - if ((location == CELL_XLOW) && (bndry->bx > 0)) { - extrap_start = 1; - } else if ((location == CELL_YLOW) && (bndry->by > 0)) { - extrap_start = 1; - } - } + // Can use no_extra_interpolate argument to skip the extra interpolation when we + // want to extrapolate the Christoffel symbol terms which come from derivatives so + // don't have the extra point set already + const bool interpolating_into_x_boundary = + ((location == CELL_XLOW) and (bndry->bx > 0)); + const bool interpolating_into_y_boundary = + ((location == CELL_YLOW) and (bndry->by > 0)); + const int extrap_start = + (not no_extra_interpolate) + and (interpolating_into_x_boundary or interpolating_into_y_boundary) + ? 1 + : 0; for (bndry->first(); !bndry->isDone(); bndry->next1d()) { // interpolate extra boundary point that is missed by interp_to, if // necessary. @@ -76,19 +85,21 @@ Field2D interpolateAndExtrapolate(const Field2D& f, CELL_LOC location, bool extr ASSERT1(bndry->bx == 0 or localmesh->xstart > 1); ASSERT1(bndry->by == 0 or localmesh->ystart > 1); // note that either bx or by is >0 here - result(bndry->x, bndry->y) = - (9. - * (f(bndry->x - bndry->bx, bndry->y - bndry->by) - + f(bndry->x, bndry->y)) - - f(bndry->x - 2 * bndry->bx, bndry->y - 2 * bndry->by) - - f(bndry->x + bndry->bx, bndry->y + bndry->by)) - / 16.; + for (int zi = zstart; zi <= zend; ++zi) { + result(bndry->x, bndry->y, zi) = + (9. + * (f(bndry->x - bndry->bx, bndry->y - bndry->by, zi) + + f(bndry->x, bndry->y, zi)) + - f(bndry->x - (2 * bndry->bx), bndry->y - (2 * bndry->by), zi) + - f(bndry->x + bndry->bx, bndry->y + bndry->by, zi)) + / 16.; + } } // set boundary guard cells - if ((bndry->bx != 0 && localmesh->GlobalNx - 2 * bndry->width >= 3) + if ((bndry->bx != 0 && localmesh->GlobalNx - (2 * bndry->width) >= 3) || (bndry->by != 0 - && localmesh->GlobalNy - localmesh->numberOfYBoundaries() * bndry->width + && localmesh->GlobalNy - (localmesh->numberOfYBoundaries() * bndry->width) >= 3)) { if (bndry->bx != 0 && localmesh->LocalNx == 1 && bndry->width == 1) { throw BoutException( @@ -106,17 +117,22 @@ Field2D interpolateAndExtrapolate(const Field2D& f, CELL_LOC location, bool extr } // extrapolate into boundary guard cells if there are enough grid points for (int i = extrap_start; i < bndry->width; i++) { - int xi = bndry->x + i * bndry->bx; - int yi = bndry->y + i * bndry->by; - result(xi, yi) = 3.0 * result(xi - bndry->bx, yi - bndry->by) - - 3.0 * result(xi - 2 * bndry->bx, yi - 2 * bndry->by) - + result(xi - 3 * bndry->bx, yi - 3 * bndry->by); + const int xi = bndry->x + (i * bndry->bx); + const int yi = bndry->y + (i * bndry->by); + for (int zi = zstart; zi <= zend; ++zi) { + result(xi, yi, zi) = + (3.0 * result(xi - bndry->bx, yi - bndry->by, zi)) + - (3.0 * result(xi - (2 * bndry->bx), yi - (2 * bndry->by), zi)) + + result(xi - (3 * bndry->bx), yi - (3 * bndry->by), zi); + } } } else { // not enough grid points to extrapolate, set equal to last grid point for (int i = extrap_start; i < bndry->width; i++) { - result(bndry->x + i * bndry->bx, bndry->y + i * bndry->by) = - result(bndry->x - bndry->bx, bndry->y - bndry->by); + for (int zi = zstart; zi <= zend; ++zi) { + result(bndry->x + (i * bndry->bx), bndry->y + (i * bndry->by), zi) = + result(bndry->x - bndry->bx, bndry->y - bndry->by, zi); + } } } } @@ -134,36 +150,40 @@ Field2D interpolateAndExtrapolate(const Field2D& f, CELL_LOC location, bool extr // Invalidate corner guard cells for (int i = 0; i < localmesh->xstart; i++) { for (int j = 0; j < localmesh->ystart; j++) { - result(i, j) = BoutNaN; - result(i, localmesh->LocalNy - 1 - j) = BoutNaN; - result(localmesh->LocalNx - 1 - i, j) = BoutNaN; - result(localmesh->LocalNx - 1 - i, localmesh->LocalNy - 1 - j) = BoutNaN; + for (int k = 0; k < f.getNz(); ++k) { + result(i, j, k) = BoutNaN; + result(i, localmesh->LocalNy - 1 - j, k) = BoutNaN; + result(localmesh->LocalNx - 1 - i, j, k) = BoutNaN; + result(localmesh->LocalNx - 1 - i, localmesh->LocalNy - 1 - j, k) = BoutNaN; + } } } } #endif - - return result; } -#if BOUT_USE_METRIC_3D -Field3D interpolateAndExtrapolate(const Field3D& f_, CELL_LOC location, - bool extrapolate_x, bool extrapolate_y, - bool no_extra_interpolate, ParallelTransform* pt_) { +/// Interpolate a Field2D to a new CELL_LOC with interp_to. +/// Communicates to set internal guard cells. +/// Boundary guard cells are set by extrapolating from the grid, like +/// 'free_o3' boundary conditions +/// Corner guard cells are set to BoutNaN +template > +auto interpolateAndExtrapolate(const T& f_, CELL_LOC location, bool extrapolate_x, + bool extrapolate_y, bool no_extra_interpolate, + ParallelTransform* pt_ = nullptr) -> T { Mesh* localmesh = f_.getMesh(); - Field3D result; - Field3D f = f_; - ParallelTransform* pt_f; - if (f.getCoordinates() == nullptr) { - // if input f is member of the Coordinates we are currently constructing, it will not - // have Coordinates and needs to use the passed-in ParallelTransform - pt_f = pt_; - } else { - // if input f is from Coordinates at a different location, it will have its own - // Coordinates, and we should use its ParallelTransform - pt_f = &f.getCoordinates()->getParallelTransform(); - } + T f = f_; + + // If input f is member of the Coordinates we are currently constructing, it will not + // have Coordinates and needs to use the passed-in ParallelTransform. + // Otherwise, if input f is from Coordinates at a different location, it will have its own + // Coordinates, and we should use its ParallelTransform + ParallelTransform* pt_f = + (not bout::build::use_metric_3d) or f.getCoordinates() == nullptr + ? pt_ + : &f.getCoordinates()->getParallelTransform(); + if (f.getDirectionY() != YDirectionType::Standard) { if (pt_f->canToFromFieldAligned()) { f = pt_f->fromFieldAligned(f); @@ -171,15 +191,15 @@ Field3D interpolateAndExtrapolate(const Field3D& f_, CELL_LOC location, f.setDirectionY(YDirectionType::Standard); } } + + T result; if (location == CELL_YLOW and f.getLocation() != CELL_YLOW) { auto f_aligned = pt_f->toFieldAligned(f, "RGN_NOX"); result = interp_to(f_aligned, location, "RGN_NOBNDRY"); - ParallelTransform* pt_result; - if (result.getCoordinates() == nullptr) { - pt_result = pt_; - } else { - pt_result = &result.getCoordinates()->getParallelTransform(); - } + const auto* result_coords = result.getCoordinates(); + ParallelTransform* pt_result = + result_coords == nullptr ? pt_ : &result_coords->getParallelTransform(); + result = pt_result->fromFieldAligned(result, "RGN_NOBNDRY"); } else { result = interp_to(f, location, "RGN_NOBNDRY"); @@ -194,144 +214,22 @@ Field3D interpolateAndExtrapolate(const Field3D& f_, CELL_LOC location, result.allocate(); localmesh->communicate_no_slices(result); - // Extrapolate into boundaries (if requested) so that differential geometry - // terms can be interpolated if necessary - // Note: cannot use applyBoundary("free_o3") here because applyBoundary() - // would try to create a new Coordinates object since we have not finished - // initializing yet, leading to an infinite recursion. - // Also, here we interpolate for the boundary points at xstart/ystart and - // (xend+1)/(yend+1) instead of extrapolating. - for (auto& bndry : localmesh->getBoundaries()) { - if ((extrapolate_x and bndry->bx != 0) or (extrapolate_y and bndry->by != 0)) { - int extrap_start = 0; - if (not no_extra_interpolate) { - // Can use no_extra_interpolate argument to skip the extra interpolation when we - // want to extrapolate the Christoffel symbol terms which come from derivatives so - // don't have the extra point set already - if ((location == CELL_XLOW) && (bndry->bx > 0)) { - extrap_start = 1; - } else if ((location == CELL_YLOW) && (bndry->by > 0)) { - extrap_start = 1; - } - } - for (bndry->first(); !bndry->isDone(); bndry->next1d()) { - // interpolate extra boundary point that is missed by interp_to, if - // necessary. - // Only interpolate this point if we are actually changing location. E.g. - // when we use this function to extrapolate J and Bxy on staggered grids, - // this point should already be set correctly because the metric - // components have been interpolated to here. - if (extrap_start > 0 and f.getLocation() != location) { - ASSERT1(bndry->bx == 0 or localmesh->xstart > 1); - ASSERT1(bndry->by == 0 or localmesh->ystart > 1); - // note that either bx or by is >0 here - for (int zi = localmesh->zstart; zi <= localmesh->zend; ++zi) { - result(bndry->x, bndry->y, zi) = - (9. - * (f(bndry->x - bndry->bx, bndry->y - bndry->by, zi) - + f(bndry->x, bndry->y, zi)) - - f(bndry->x - 2 * bndry->bx, bndry->y - 2 * bndry->by, zi) - - f(bndry->x + bndry->bx, bndry->y + bndry->by, zi)) - / 16.; - } - } - // set boundary guard cells - if ((bndry->bx != 0 && localmesh->GlobalNx - 2 * bndry->width >= 3) - || (bndry->by != 0 && localmesh->GlobalNy - 2 * bndry->width >= 3)) { - if (bndry->bx != 0 && localmesh->LocalNx == 1 && bndry->width == 1) { - throw BoutException( - "Not enough points in the x-direction on this " - "processor for extrapolation needed to use staggered grids. " - "Increase number of x-guard cells MXG or decrease number of " - "processors in the x-direction NXPE."); - } - if (bndry->by != 0 && localmesh->LocalNy == 1 && bndry->width == 1) { - throw BoutException( - "Not enough points in the y-direction on this " - "processor for extrapolation needed to use staggered grids. " - "Increase number of y-guard cells MYG or decrease number of " - "processors in the y-direction NYPE."); - } - // extrapolate into boundary guard cells if there are enough grid points - for (int i = extrap_start; i < bndry->width; i++) { - int xi = bndry->x + i * bndry->bx; - int yi = bndry->y + i * bndry->by; - for (int zi = localmesh->zstart; zi <= localmesh->zend; ++zi) { - result(xi, yi, zi) = - 3.0 * result(xi - bndry->bx, yi - bndry->by, zi) - - 3.0 * result(xi - 2 * bndry->bx, yi - 2 * bndry->by, zi) - + result(xi - 3 * bndry->bx, yi - 3 * bndry->by, zi); - } - } - } else { - // not enough grid points to extrapolate, set equal to last grid point - for (int i = extrap_start; i < bndry->width; i++) { - for (int zi = localmesh->zstart; zi <= localmesh->zend; ++zi) { - result(bndry->x + i * bndry->bx, bndry->y + i * bndry->by, zi) = - result(bndry->x - bndry->bx, bndry->y - bndry->by, zi); - } - } - } - } - } - } -#if CHECK > 0 - if (not( - // if include_corner_cells=true, then we extrapolate valid data into the - // corner cells if they are not already filled - localmesh->include_corner_cells - - // if we are not extrapolating at all, the corner cells should contain valid - // data - or (not extrapolate_x and not extrapolate_y))) { - // Invalidate corner guard cells - for (int i = 0; i < localmesh->xstart; i++) { - for (int j = 0; j < localmesh->ystart; j++) { - for (int k = 0; k < localmesh->LocalNz; ++k) { - result(i, j, k) = BoutNaN; - result(i, localmesh->LocalNy - 1 - j, k) = BoutNaN; - result(localmesh->LocalNx - 1 - i, j, k) = BoutNaN; - result(localmesh->LocalNx - 1 - i, localmesh->LocalNy - 1 - j, k) = BoutNaN; - } - } - } - } -#endif // CHECK > 0 + fillGuards_impl(result, location, f, extrapolate_x, extrapolate_y, + no_extra_interpolate); return result; } -#endif // BOUT_USE_METRIC_3D // If the CELL_CENTRE variable was read, the staggered version is required to // also exist for consistency void checkStaggeredGet(Mesh* mesh, const std::string& name, const std::string& suffix) { if (mesh->sourceHasVar(name) != mesh->sourceHasVar(name + suffix)) { - throw BoutException("Attempting to read staggered fields from grid, but " + name - + " is not present in both CELL_CENTRE and staggered versions."); + throw BoutException("Attempting to read staggered fields from grid, but '{}'" + " is not present in both CELL_CENTRE and staggered versions.", + name); } } -// convenience function for repeated code -int getAtLoc(Mesh* mesh, Coordinates::FieldMetric& var, const std::string& name, - const std::string& suffix, CELL_LOC location, BoutReal default_value = 0.) { - - checkStaggeredGet(mesh, name, suffix); - int result = mesh->get(var, name + suffix, default_value, false, location); - - return result; -} - -int getAtLocAndFillGuards(Mesh* mesh, Coordinates::FieldMetric& var, - const std::string& name, const std::string& suffix, - CELL_LOC location, BoutReal default_value, bool extrapolate_x, - bool extrapolate_y, bool no_extra_interpolate, - ParallelTransform* pt) { - auto ret = getAtLoc(mesh, var, name, suffix, location, default_value); - var = interpolateAndExtrapolate(var, location, extrapolate_x, extrapolate_y, - no_extra_interpolate, pt); - return ret; -} - std::string getLocationSuffix(CELL_LOC location) { switch (location) { case CELL_CENTRE: { @@ -381,8 +279,92 @@ Coordinates::Coordinates(Mesh* mesh, Options* options) G1_13(mesh), G1_23(mesh), G2_11(mesh), G2_22(mesh), G2_33(mesh), G2_12(mesh), G2_13(mesh), G2_23(mesh), G3_11(mesh), G3_22(mesh), G3_33(mesh), G3_12(mesh), G3_13(mesh), G3_23(mesh), G1(mesh), G2(mesh), G3(mesh), ShiftTorsion(mesh), - IntShiftTorsion(mesh), localmesh(mesh), location(CELL_CENTRE) { + IntShiftTorsion(mesh), nz(mesh->LocalNz), localmesh(mesh), location(CELL_CENTRE) { + + readFromMesh(options, ""); +} + +Coordinates::Coordinates(Mesh* mesh, Options* options, const CELL_LOC loc, + const Coordinates* coords_in, bool force_interpolate_from_centre) + : dx(1., mesh), dy(1., mesh), dz(1., mesh), d1_dx(mesh), d1_dy(mesh), d1_dz(mesh), + J(1., mesh), Bxy(1., mesh), + // Identity metric tensor + g11(1., mesh), g22(1., mesh), g33(1., mesh), g12(0, mesh), g13(0, mesh), + g23(0, mesh), g_11(1., mesh), g_22(1., mesh), g_33(1., mesh), g_12(0, mesh), + g_13(0, mesh), g_23(0, mesh), G1_11(mesh), G1_22(mesh), G1_33(mesh), G1_12(mesh), + G1_13(mesh), G1_23(mesh), G2_11(mesh), G2_22(mesh), G2_33(mesh), G2_12(mesh), + G2_13(mesh), G2_23(mesh), G3_11(mesh), G3_22(mesh), G3_33(mesh), G3_12(mesh), + G3_13(mesh), G3_23(mesh), G1(mesh), G2(mesh), G3(mesh), ShiftTorsion(mesh), + IntShiftTorsion(mesh), nz(mesh->LocalNz), localmesh(mesh), location(loc) { + + const std::string suffix = getLocationSuffix(location); + + if (not force_interpolate_from_centre and mesh->sourceHasVar("dx" + suffix)) { + readFromMesh(options, suffix); + } else { + // Interpolate fields from coords_in + + if (isUniform(coords_in->dz)) { + dz = coords_in->dz; + dz.setLocation(location); + } else { + throw BoutException( + "We are asked to transform dz to get dz before we have a transform, which " + "might require dz!\nPlease provide a dz for the staggered quantity!"); + } + setParallelTransform(options); + + auto interpField = [&, this](const FieldMetric& f) -> FieldMetric { + return interpolateAndExtrapolate(f, location, true, true, false, transform.get()); + }; + + dx = interpField(coords_in->dx); + dy = interpField(coords_in->dy); + // not really needed - we have used dz already ... + dz = interpField(coords_in->dz); + // Diagonal components of metric tensor g^{ij} + g11 = interpField(coords_in->g11); + g22 = interpField(coords_in->g22); + g33 = interpField(coords_in->g33); + + // Off-diagonal elements. + g12 = interpField(coords_in->g12); + g13 = interpField(coords_in->g13); + g23 = interpField(coords_in->g23); + + // 3x3 matrix inversion can exaggerate small interpolation errors, so it is + // more robust to interpolate and extrapolate derived quantities directly, + // rather than deriving from interpolated/extrapolated covariant metric + // components + g_11 = interpField(coords_in->g_11); + g_22 = interpField(coords_in->g_22); + g_33 = interpField(coords_in->g_33); + g_12 = interpField(coords_in->g_12); + g_13 = interpField(coords_in->g_13); + g_23 = interpField(coords_in->g_23); + + // Check input metrics + checkContravariant(); + checkCovariant(); + + J = interpField(coords_in->J); + Bxy = interpField(coords_in->Bxy); + + bout::checkFinite(J, "The Jacobian", "RGN_NOCORNERS"); + bout::checkPositive(J, "The Jacobian", "RGN_NOCORNERS"); + bout::checkFinite(Bxy, "Bxy", "RGN_NOCORNERS"); + bout::checkPositive(Bxy, "Bxy", "RGN_NOCORNERS"); + + ShiftTorsion = interpField(coords_in->ShiftTorsion); + + if (mesh->IncIntShear) { + IntShiftTorsion = interpField(coords_in->IntShiftTorsion); + } + } +} + +void Coordinates::readFromMesh(Options* options, const std::string& suffix) { if (options == nullptr) { options = Options::getRoot()->getSection("mesh"); } @@ -392,9 +374,9 @@ Coordinates::Coordinates(Mesh* mesh, Options* options) // smooth at all the boundaries. const bool extrapolate_x = - (*options)["extrapolate_x"].withDefault(not mesh->sourceHasXBoundaryGuards()); + (*options)["extrapolate_x"].withDefault(not localmesh->sourceHasXBoundaryGuards()); const bool extrapolate_y = - (*options)["extrapolate_y"].withDefault(not mesh->sourceHasYBoundaryGuards()); + (*options)["extrapolate_y"].withDefault(not localmesh->sourceHasYBoundaryGuards()); if (extrapolate_x) { output_warn.write(_("WARNING: extrapolating input mesh quantities into x-boundary " @@ -406,12 +388,39 @@ Coordinates::Coordinates(Mesh* mesh, Options* options) "cells. Set option extrapolate_y=false to disable this.\n")); } - mesh->get(dx, "dx", 1.0, false); - mesh->get(dy, "dy", 1.0, false); + // Some helper functions that let us avoid passing the same arguments every time - nz = mesh->LocalNz; + // Read from the mesh and transform from field aligned if required + auto readField = [this, &suffix](const std::string& name, + BoutReal def_value) -> FieldMetric { + checkStaggeredGet(localmesh, name, suffix); + + FieldMetric field{localmesh}; + localmesh->get(field, name + suffix, def_value, false, location); + if (field.getDirectionY() == YDirectionType::Aligned + and transform->canToFromFieldAligned()) { + return transform->fromFieldAligned(field); + } + return field.setDirectionY(YDirectionType::Standard); + }; + + // Just fill in the guard cells (via extrapolation) for an existing field + auto fillGuards = [&, this](FieldMetric& field) { + // Passing `field` in here twice is gross but ok because the second argument + // is only used when interpolating, and we're not interpolating here + fillGuards_impl(field, location, field, extrapolate_x, extrapolate_y); + }; + + // Read the field, transform if required, and fill in the guards + auto readAndFillGuards = [&](const std::string& name, + BoutReal def_value) -> FieldMetric { + auto field = readField(name, def_value); + fillGuards(field); + return field; + }; { + // dz auto& options = Options::root(); const bool has_zperiod = options.isSet("zperiod"); const auto zmin = has_zperiod ? 0.0 : options["ZMIN"].withDefault(0.0); @@ -420,60 +429,45 @@ Coordinates::Coordinates(Mesh* mesh, Options* options) const auto default_dz = (zmax - zmin) * TWOPI / nz; - mesh->get(dz, "dz", default_dz, false); + // We can't use the helper functions here because in 3D we might (always?) + // need to transform from field aligned -- which requires dz! So we have to + // read it "plain" first... + localmesh->get(dz, "dz" + suffix, default_dz, false, location); } - // required early for differentiation. + // ...then we can set the transform (required early for differentiation)... setParallelTransform(options); + // ...and finally we can transform/interpolate/fill in guards dz = interpolateAndExtrapolate(dz, location, extrapolate_x, extrapolate_y, false, transform.get()); - dx = interpolateAndExtrapolate(dx, location, extrapolate_x, extrapolate_y, false, - transform.get()); - - if (mesh->periodicX) { - mesh->communicate_no_slices(dx); - } - dy = interpolateAndExtrapolate(dy, location, extrapolate_x, extrapolate_y, false, - transform.get()); + // everything else from this point on can use our helper functions + dx = readAndFillGuards("dx", 1.0); - auto getUnaligned = [this](auto& field, const std::string& name, - BoutReal default_value) { - localmesh->get(field, name, default_value, false); - if (field.getDirectionY() == YDirectionType::Aligned - and transform->canToFromFieldAligned()) { - return transform->fromFieldAligned(field); - } else { - return field.setDirectionY(YDirectionType::Standard); - } - }; + if (localmesh->periodicX) { + localmesh->communicate_no_slices(dx); + } - auto getUnalignedAtLocation = [this, extrapolate_x, extrapolate_y, - getUnaligned](auto& field, const std::string& name, - BoutReal default_value) { - field = getUnaligned(field, name, default_value); - return interpolateAndExtrapolate(field, location, extrapolate_x, extrapolate_y, false, - transform.get()); - }; + dy = readAndFillGuards("dy", 1.0); // Diagonal components of metric tensor g^{ij} (default to 1) - g11 = getUnalignedAtLocation(g11, "g11", 1.0); - g22 = getUnalignedAtLocation(g22, "g22", 1.0); - g33 = getUnalignedAtLocation(g33, "g33", 1.0); + g11 = readAndFillGuards("g11", 1.0); + g22 = readAndFillGuards("g22", 1.0); + g33 = readAndFillGuards("g33", 1.0); // Off-diagonal elements. Default to 0 - g12 = getUnalignedAtLocation(g12, "g12", 0.0); - g13 = getUnalignedAtLocation(g13, "g13", 0.0); - g23 = getUnalignedAtLocation(g23, "g23", 0.0); + g12 = readAndFillGuards("g12", 0.0); + g13 = readAndFillGuards("g13", 0.0); + g23 = readAndFillGuards("g23", 0.0); // Check input metrics checkContravariant(); - /// Find covariant metric components + // Find covariant metric components auto covariant_component_names = {"g_11", "g_22", "g_33", "g_12", "g_13", "g_23"}; - auto source_has_component = [&mesh](const std::string& name) { - return mesh->sourceHasVar(name); + auto source_has_component = [&suffix, this](const std::string& name) { + return localmesh->sourceHasVar(name + suffix); }; // Check if any of the components are present if (std::any_of(begin(covariant_component_names), end(covariant_component_names), @@ -481,12 +475,12 @@ Coordinates::Coordinates(Mesh* mesh, Options* options) // Check that all components are present if (std::all_of(begin(covariant_component_names), end(covariant_component_names), source_has_component)) { - g_11 = getUnaligned(g_11, "g_11", 1.0); - g_22 = getUnaligned(g_22, "g_22", 1.0); - g_33 = getUnaligned(g_33, "g_33", 1.0); - g_12 = getUnaligned(g_12, "g_12", 0.0); - g_13 = getUnaligned(g_13, "g_13", 0.0); - g_23 = getUnaligned(g_23, "g_23", 0.0); + g_11 = readField("g_11", 1.0); + g_22 = readField("g_22", 1.0); + g_33 = readField("g_33", 1.0); + g_12 = readField("g_12", 0.0); + g_13 = readField("g_13", 0.0); + g_23 = readField("g_23", 0.0); output_warn.write("\tWARNING! Covariant components of metric tensor set manually. " "Contravariant components NOT recalculated\n"); @@ -495,56 +489,47 @@ Coordinates::Coordinates(Mesh* mesh, Options* options) output_warn.write("Not all covariant components of metric tensor found. " "Calculating all from the contravariant tensor\n"); /// Calculate contravariant metric components if not found - if (calcCovariant("RGN_NOCORNERS")) { + if (calcCovariant("RGN_NOCORNERS") != 0) { throw BoutException("Error in calcCovariant call"); } } } else { /// Calculate contravariant metric components if not found - if (calcCovariant("RGN_NOCORNERS")) { + if (calcCovariant("RGN_NOCORNERS") != 0) { throw BoutException("Error in calcCovariant call"); } } // More robust to extrapolate derived quantities directly, rather than // deriving from extrapolated covariant metric components - g_11 = interpolateAndExtrapolate(g_11, location, extrapolate_x, extrapolate_y, false, - transform.get()); - g_22 = interpolateAndExtrapolate(g_22, location, extrapolate_x, extrapolate_y, false, - transform.get()); - g_33 = interpolateAndExtrapolate(g_33, location, extrapolate_x, extrapolate_y, false, - transform.get()); - g_12 = interpolateAndExtrapolate(g_12, location, extrapolate_x, extrapolate_y, false, - transform.get()); - g_13 = interpolateAndExtrapolate(g_13, location, extrapolate_x, extrapolate_y, false, - transform.get()); - g_23 = interpolateAndExtrapolate(g_23, location, extrapolate_x, extrapolate_y, false, - transform.get()); + fillGuards(g_11); + fillGuards(g_22); + fillGuards(g_33); + fillGuards(g_12); + fillGuards(g_13); + fillGuards(g_23); // Check covariant metrics checkCovariant(); - /// Calculate Jacobian and Bxy - if (jacobian()) { - throw BoutException("Error in jacobian call"); - } + // Calculate Jacobian and Bxy + jacobian(); // Attempt to read J from the grid file - auto Jcalc = J; - if (mesh->get(J, "J", 0.0, false)) { - output_warn.write( - "\tWARNING: Jacobian 'J' not found. Calculating from metric tensor\n"); - J = Jcalc; - } else { - J = interpolateAndExtrapolate(J, location, extrapolate_x, extrapolate_y, false, - transform.get()); - - // Compare calculated and loaded values + if (localmesh->sourceHasVar("J" + suffix)) { + // Copy value previously calculated from metric components, use to check + // read in value + const auto Jcalc = J; + J = readAndFillGuards("J", 0.0); output_warn.write("\tMaximum difference in J is {:e}\n", max(abs(J - Jcalc))); - mesh->communicate_no_slices(J); + localmesh->communicate_no_slices(J); // Re-evaluate Bxy using new J Bxy = sqrt(g_22) / J; + fillGuards(Bxy); + } else { + output_warn.write( + "\tWARNING: Jacobian 'J' not found. Calculating from metric tensor\n"); } // Check jacobian @@ -555,334 +540,43 @@ Coordinates::Coordinates(Mesh* mesh, Options* options) } // Attempt to read Bxy from the grid file - auto Bcalc = Bxy; - if (mesh->get(Bxy, "Bxy", 0.0, false)) { + if (localmesh->sourceHasVar("Bxy" + suffix)) { + // Copy value previously calculated from J, use to check read in value + const auto Bcalc = Bxy; + Bxy = readAndFillGuards("Bxy", 0.0); + output_warn.write("\tMaximum difference in Bxy is {:e}\n", max(abs(Bxy - Bcalc))); + + } else { output_warn.write("\tWARNING: Magnitude of B field 'Bxy' not found. Calculating from " "metric tensor\n"); - Bxy = Bcalc; - } else { - - Bxy = interpolateAndExtrapolate(Bxy, location, extrapolate_x, extrapolate_y, false, - transform.get()); - output_warn.write("\tMaximum difference in Bxy is {:e}\n", max(abs(Bxy - Bcalc))); } // Check Bxy bout::checkFinite(Bxy, "Bxy", "RGN_NOCORNERS"); bout::checkPositive(Bxy, "Bxy", "RGN_NOCORNERS"); - if (mesh->get(ShiftTorsion, "ShiftTorsion", 0.0, false)) { + if (localmesh->get(ShiftTorsion, "ShiftTorsion" + suffix, 0.0, false) != 0) { output_warn.write( "\tWARNING: No Torsion specified for zShift. Derivatives may not be correct\n"); ShiftTorsion = 0.0; } - ShiftTorsion = interpolateAndExtrapolate(ShiftTorsion, location, extrapolate_x, - extrapolate_y, false, transform.get()); + fillGuards(ShiftTorsion); ////////////////////////////////////////////////////// - if (mesh->IncIntShear) { - if (mesh->get(IntShiftTorsion, "IntShiftTorsion", 0.0, false)) { + if (localmesh->IncIntShear) { + if (localmesh->get(IntShiftTorsion, "IntShiftTorsion", 0.0, false) != 0) { output_warn.write("\tWARNING: No Integrated torsion specified\n"); } - IntShiftTorsion = interpolateAndExtrapolate(IntShiftTorsion, location, extrapolate_x, - extrapolate_y, false, transform.get()); + fillGuards(IntShiftTorsion); } else { // IntShiftTorsion will not be used, but set to zero to avoid uninitialized field IntShiftTorsion = 0.; } } -Coordinates::Coordinates(Mesh* mesh, Options* options, const CELL_LOC loc, - const Coordinates* coords_in, bool force_interpolate_from_centre) - : dx(1., mesh), dy(1., mesh), dz(1., mesh), d1_dx(mesh), d1_dy(mesh), d1_dz(mesh), - J(1., mesh), Bxy(1., mesh), - // Identity metric tensor - g11(1., mesh), g22(1., mesh), g33(1., mesh), g12(0, mesh), g13(0, mesh), - g23(0, mesh), g_11(1., mesh), g_22(1., mesh), g_33(1., mesh), g_12(0, mesh), - g_13(0, mesh), g_23(0, mesh), G1_11(mesh), G1_22(mesh), G1_33(mesh), G1_12(mesh), - G1_13(mesh), G1_23(mesh), G2_11(mesh), G2_22(mesh), G2_33(mesh), G2_12(mesh), - G2_13(mesh), G2_23(mesh), G3_11(mesh), G3_22(mesh), G3_33(mesh), G3_12(mesh), - G3_13(mesh), G3_23(mesh), G1(mesh), G2(mesh), G3(mesh), ShiftTorsion(mesh), - IntShiftTorsion(mesh), localmesh(mesh), location(loc) { - - std::string suffix = getLocationSuffix(location); - - nz = mesh->LocalNz; - - // Default to true in case staggered quantities are not read from file - bool extrapolate_x = true; - bool extrapolate_y = true; - - if (!force_interpolate_from_centre && mesh->sourceHasVar("dx" + suffix)) { - - extrapolate_x = not mesh->sourceHasXBoundaryGuards(); - extrapolate_y = not mesh->sourceHasYBoundaryGuards(); - - if (extrapolate_x) { - output_warn.write(_("WARNING: extrapolating input mesh quantities into x-boundary " - "cells\n")); - } - - if (extrapolate_y) { - output_warn.write(_("WARNING: extrapolating input mesh quantities into y-boundary " - "cells\n")); - } - - { - auto& options = Options::root(); - const bool has_zperiod = options.isSet("zperiod"); - const auto zmin = has_zperiod ? 0.0 : options["ZMIN"].withDefault(0.0); - const auto zmax = has_zperiod ? 1.0 / options["zperiod"].withDefault(1.0) - : options["ZMAX"].withDefault(1.0); - - const auto default_dz = (zmax - zmin) * TWOPI / nz; - getAtLoc(mesh, dz, "dz", suffix, location, default_dz); - } - setParallelTransform(options); - - dz = interpolateAndExtrapolate(dz, location, extrapolate_x, extrapolate_y, false, - transform.get()); - - getAtLocAndFillGuards(mesh, dx, "dx", suffix, location, 1.0, extrapolate_x, - extrapolate_y, false, transform.get()); - - if (mesh->periodicX) { - mesh->communicate_no_slices(dx); - } - - getAtLocAndFillGuards(mesh, dy, "dy", suffix, location, 1.0, extrapolate_x, - extrapolate_y, false, transform.get()); - - // grid data source has staggered fields, so read instead of interpolating - // Diagonal components of metric tensor g^{ij} (default to 1) - getAtLocAndFillGuards(mesh, g11, "g11", suffix, location, 1.0, extrapolate_x, - extrapolate_y, false, transform.get()); - getAtLocAndFillGuards(mesh, g22, "g22", suffix, location, 1.0, extrapolate_x, - extrapolate_y, false, transform.get()); - getAtLocAndFillGuards(mesh, g33, "g33", suffix, location, 1.0, extrapolate_x, - extrapolate_y, false, transform.get()); - getAtLocAndFillGuards(mesh, g12, "g12", suffix, location, 0.0, extrapolate_x, - extrapolate_y, false, transform.get()); - getAtLocAndFillGuards(mesh, g13, "g13", suffix, location, 0.0, extrapolate_x, - extrapolate_y, false, transform.get()); - getAtLocAndFillGuards(mesh, g23, "g23", suffix, location, 0.0, extrapolate_x, - extrapolate_y, false, transform.get()); - - // Check input metrics - checkContravariant(); - - /// Find covariant metric components - auto covariant_component_names = {"g_11", "g_22", "g_33", "g_12", "g_13", "g_23"}; - auto source_has_component = [&suffix, &mesh](const std::string& name) { - return mesh->sourceHasVar(name + suffix); - }; - // Check if any of the components are present - if (std::any_of(begin(covariant_component_names), end(covariant_component_names), - source_has_component)) { - // Check that all components are present - if (std::all_of(begin(covariant_component_names), end(covariant_component_names), - source_has_component)) { - - getAtLoc(mesh, g_11, "g_11", suffix, location); - getAtLoc(mesh, g_22, "g_22", suffix, location); - getAtLoc(mesh, g_33, "g_33", suffix, location); - getAtLoc(mesh, g_12, "g_12", suffix, location); - getAtLoc(mesh, g_13, "g_13", suffix, location); - getAtLoc(mesh, g_23, "g_23", suffix, location); - - output_warn.write( - "\tWARNING! Staggered covariant components of metric tensor set manually. " - "Contravariant components NOT recalculated\n"); - - } else { - output_warn.write( - "Not all staggered covariant components of metric tensor found. " - "Calculating all from the contravariant tensor\n"); - /// Calculate contravariant metric components if not found - if (calcCovariant("RGN_NOCORNERS")) { - throw BoutException("Error in staggered calcCovariant call"); - } - } - } else { - /// Calculate contravariant metric components if not found - if (calcCovariant("RGN_NOCORNERS")) { - throw BoutException("Error in staggered calcCovariant call"); - } - } - // More robust to extrapolate derived quantities directly, rather than - // deriving from extrapolated covariant metric components - g_11 = interpolateAndExtrapolate(g_11, location, extrapolate_x, extrapolate_y, false, - transform.get()); - g_22 = interpolateAndExtrapolate(g_22, location, extrapolate_x, extrapolate_y, false, - transform.get()); - g_33 = interpolateAndExtrapolate(g_33, location, extrapolate_x, extrapolate_y, false, - transform.get()); - g_12 = interpolateAndExtrapolate(g_12, location, extrapolate_x, extrapolate_y, false, - transform.get()); - g_13 = interpolateAndExtrapolate(g_13, location, extrapolate_x, extrapolate_y, false, - transform.get()); - g_23 = interpolateAndExtrapolate(g_23, location, extrapolate_x, extrapolate_y, false, - transform.get()); - - // Check covariant metrics - checkCovariant(); - - /// Calculate Jacobian and Bxy - if (jacobian()) { - throw BoutException("Error in jacobian call while constructing staggered " - "Coordinates"); - } - - // Attempt to read J from the grid file - auto Jcalc = J; - if (getAtLoc(mesh, J, "J", suffix, location)) { - output_warn.write( - "\tWARNING: Jacobian 'J_{:s}' not found. Calculating from metric tensor\n", - suffix); - J = Jcalc; - } else { - J = interpolateAndExtrapolate(J, location, extrapolate_x, extrapolate_y, false, - transform.get()); - - // Compare calculated and loaded values - output_warn.write("\tMaximum difference in J is %e\n", max(abs(J - Jcalc))); - - // Re-evaluate Bxy using new J - Bxy = sqrt(g_22) / J; - } - - // Check jacobian - bout::checkFinite(J, "J" + suffix, "RGN_NOCORNERS"); - bout::checkPositive(J, "J" + suffix, "RGN_NOCORNERS"); - if (min(abs(J)) < 1.0e-10) { - throw BoutException("\tERROR: Jacobian{:s} becomes very small\n", suffix); - } - - // Attempt to read Bxy from the grid file - auto Bcalc = Bxy; - if (getAtLoc(mesh, Bxy, "Bxy", suffix, location)) { - output_warn.write( - "\tWARNING: Magnitude of B field 'Bxy_{:s}' not found. Calculating " - " from metric tensor\n", - suffix); - Bxy = Bcalc; - } else { - Bxy = interpolateAndExtrapolate(Bxy, location, extrapolate_x, extrapolate_y, false, - transform.get()); - - output_warn.write("\tMaximum difference in Bxy is %e\n", max(abs(Bxy - Bcalc))); - } - - // Check Bxy - bout::checkFinite(Bxy, "Bxy" + suffix, "RGN_NOCORNERS"); - bout::checkPositive(Bxy, "Bxy" + suffix, "RGN_NOCORNERS"); - - checkStaggeredGet(mesh, "ShiftTorsion", suffix); - if (mesh->get(ShiftTorsion, "ShiftTorsion" + suffix, 0.0, false)) { - output_warn.write( - "\tWARNING: No Torsion specified for zShift. Derivatives may not be correct\n"); - ShiftTorsion = 0.0; - } - ShiftTorsion.setLocation(location); - ShiftTorsion = interpolateAndExtrapolate(ShiftTorsion, location, extrapolate_x, - extrapolate_y, false, transform.get()); - - ////////////////////////////////////////////////////// - - if (mesh->IncIntShear) { - checkStaggeredGet(mesh, "IntShiftTorsion", suffix); - if (mesh->get(IntShiftTorsion, "IntShiftTorsion" + suffix, 0.0, false)) { - output_warn.write("\tWARNING: No Integrated torsion specified\n"); - IntShiftTorsion = 0.0; - } - IntShiftTorsion.setLocation(location); - IntShiftTorsion = - interpolateAndExtrapolate(IntShiftTorsion, location, extrapolate_x, - extrapolate_y, false, transform.get()); - } else { - // IntShiftTorsion will not be used, but set to zero to avoid uninitialized field - IntShiftTorsion = 0.; - } - } else { - // Interpolate fields from coords_in - - if (isUniform(coords_in->dz)) { - dz = coords_in->dz; - dz.setLocation(location); - } else { - throw BoutException( - "We are asked to transform dz to get dz before we have a transform, which " - "might require dz!\nPlease provide a dz for the staggered quantity!"); - } - setParallelTransform(options); - dx = interpolateAndExtrapolate(coords_in->dx, location, true, true, false, - transform.get()); - dy = interpolateAndExtrapolate(coords_in->dy, location, true, true, false, - transform.get()); - // not really needed - we have used dz already ... - dz = interpolateAndExtrapolate(coords_in->dz, location, true, true, false, - transform.get()); - - // Diagonal components of metric tensor g^{ij} - g11 = interpolateAndExtrapolate(coords_in->g11, location, true, true, false, - transform.get()); - g22 = interpolateAndExtrapolate(coords_in->g22, location, true, true, false, - transform.get()); - g33 = interpolateAndExtrapolate(coords_in->g33, location, true, true, false, - transform.get()); - - // Off-diagonal elements. - g12 = interpolateAndExtrapolate(coords_in->g12, location, true, true, false, - transform.get()); - g13 = interpolateAndExtrapolate(coords_in->g13, location, true, true, false, - transform.get()); - g23 = interpolateAndExtrapolate(coords_in->g23, location, true, true, false, - transform.get()); - - // 3x3 matrix inversion can exaggerate small interpolation errors, so it is - // more robust to interpolate and extrapolate derived quantities directly, - // rather than deriving from interpolated/extrapolated covariant metric - // components - g_11 = interpolateAndExtrapolate(coords_in->g_11, location, true, true, false, - transform.get()); - g_22 = interpolateAndExtrapolate(coords_in->g_22, location, true, true, false, - transform.get()); - g_33 = interpolateAndExtrapolate(coords_in->g_33, location, true, true, false, - transform.get()); - g_12 = interpolateAndExtrapolate(coords_in->g_12, location, true, true, false, - transform.get()); - g_13 = interpolateAndExtrapolate(coords_in->g_13, location, true, true, false, - transform.get()); - g_23 = interpolateAndExtrapolate(coords_in->g_23, location, true, true, false, - transform.get()); - - // Check input metrics - checkContravariant(); - checkCovariant(); - - J = interpolateAndExtrapolate(coords_in->J, location, true, true, false, - transform.get()); - Bxy = interpolateAndExtrapolate(coords_in->Bxy, location, true, true, false, - transform.get()); - - bout::checkFinite(J, "The Jacobian", "RGN_NOCORNERS"); - bout::checkPositive(J, "The Jacobian", "RGN_NOCORNERS"); - bout::checkFinite(Bxy, "Bxy", "RGN_NOCORNERS"); - bout::checkPositive(Bxy, "Bxy", "RGN_NOCORNERS"); - - ShiftTorsion = interpolateAndExtrapolate(coords_in->ShiftTorsion, location, true, - true, false, transform.get()); - - if (mesh->IncIntShear) { - IntShiftTorsion = interpolateAndExtrapolate(coords_in->IntShiftTorsion, location, - true, true, false, transform.get()); - } - } -} - void Coordinates::outputVars(Options& output_options) { - Timer time("io"); + const Timer time("io"); const std::string loc_string = (location == CELL_CENTRE) ? "" : "_" + toString(location); diff --git a/src/mesh/data/gridfromfile.cxx b/src/mesh/data/gridfromfile.cxx index 10f6f6926e..414f29c7ac 100644 --- a/src/mesh/data/gridfromfile.cxx +++ b/src/mesh/data/gridfromfile.cxx @@ -32,7 +32,7 @@ GridFile::GridFile(std::string gridfilename) * Tests whether a variable exists in the file * */ -bool GridFile::hasVar(const std::string& name) { return data.isSet(name); } +bool GridFile::hasVar(const std::string& name) const { return data.isSet(name); } /*! * Read a string from file. If the string is not @@ -120,7 +120,7 @@ bool GridFile::get(Mesh* UNUSED(m), BoutReal& rval, const std::string& name, /*! * Reads a 2D, 3D or FieldPerp field variable from a file - * + * * Successfully reads Field2D or FieldPerp if the variable in the file is 0-D or 2-D. * Successfully reads Field3D if the variable in the file is 0-D, 2-D or 3-D. */ @@ -378,7 +378,7 @@ void GridFile::readField(Mesh* m, const std::string& name, int ys, int yd, int n for (int x = xs; x < xs + nx_to_read; ++x) { for (int y = ys; y < ys + ny_to_read; ++y) { - BoutReal const value = full_var(x, y); + const BoutReal value = full_var(x, y); for (int z = 0; z < var.getNz(); z++) { var(x - xs + xd, y - ys + yd, z) = value; } diff --git a/src/mesh/data/gridfromoptions.cxx b/src/mesh/data/gridfromoptions.cxx index 662329c526..7d21548340 100644 --- a/src/mesh/data/gridfromoptions.cxx +++ b/src/mesh/data/gridfromoptions.cxx @@ -7,7 +7,9 @@ using bout::generator::Context; -bool GridFromOptions::hasVar(const std::string& name) { return options->isSet(name); } +bool GridFromOptions::hasVar(const std::string& name) const { + return options->isSet(name); +} namespace { /// Return value of \p name in \p options, using \p def as a default diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index e1e34d2f4f..16e90a7a98 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -3201,10 +3201,10 @@ RangeIterator BoutMesh::iterateBndryUpperY() const { return RangeIterator(xs, xe); } -std::vector BoutMesh::getBoundaries() { return boundary; } +std::vector BoutMesh::getBoundaries() const { return boundary; } std::vector> -BoutMesh::getBoundariesPar(BoundaryParType type) { +BoutMesh::getBoundariesPar(BoundaryParType type) const { return par_boundary[static_cast(type)]; } diff --git a/src/mesh/impls/bout/boutmesh.hxx b/src/mesh/impls/bout/boutmesh.hxx index 07453c2676..38d6dc4d91 100644 --- a/src/mesh/impls/bout/boutmesh.hxx +++ b/src/mesh/impls/bout/boutmesh.hxx @@ -166,9 +166,9 @@ public: bool hasBndryUpperY() const override { return has_boundary_upper_y; } // Boundary regions - std::vector getBoundaries() override; + std::vector getBoundaries() const override; std::vector> - getBoundariesPar(BoundaryParType type) override; + getBoundariesPar(BoundaryParType type) const override; void addBoundaryPar(std::shared_ptr bndry, BoundaryParType type) override; std::set getPossibleBoundaries() const override; diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx index 6dbbd6200b..3258496b3f 100644 --- a/tests/unit/fake_mesh.hxx +++ b/tests/unit/fake_mesh.hxx @@ -172,9 +172,9 @@ public: bool hasBndryLowerY() const override { return false; } bool hasBndryUpperY() const override { return false; } void addBoundary(BoundaryRegion* region) override { boundaries.push_back(region); } - std::vector getBoundaries() override { return boundaries; } + std::vector getBoundaries() const override { return boundaries; } std::vector> - getBoundariesPar(BoundaryParType UNUSED(type)) override { + getBoundariesPar(BoundaryParType UNUSED(type)) const override { return std::vector>(); } BoutReal GlobalX(int jx) const override { return jx; } @@ -280,7 +280,7 @@ public: /// Take an rvalue (e.g. initializer list), convert to lvalue and delegate constructor FakeGridDataSource(Options&& values) : FakeGridDataSource(values) {} - bool hasVar(const std::string& UNUSED(name)) override { return false; } + bool hasVar(const std::string& name) const override { return values.isSet(name); } bool get([[maybe_unused]] Mesh* m, std::string& sval, const std::string& name, const std::string& def = "") override {