From bdb56d88f8fe52fa1715db31b6db1949820b4dd0 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 20 Nov 2018 14:29:13 +0000 Subject: [PATCH 1/7] Use localmesh for d2x and d2y, and interpolate them to location Previously the localmesh was not passed through to the constructor of d2x and d2y in Coordinates::geometry(), so they used the global 'mesh'. Also, when d2x/d2y are read from the mesh, they are at CELL_CENTRE, so we need to interpolate them to the location of the Coordinates. --- src/mesh/coordinates.cxx | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 81c0e921e4..8b287638df 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -442,13 +442,16 @@ int Coordinates::geometry() { OPTION(Options::getRoot(), non_uniform, true); - Field2D d2x, d2y; // d^2 x / d i^2 + Field2D d2x(localmesh), d2y(localmesh); // d^2 x / d i^2 // Read correction for non-uniform meshes if (localmesh->get(d2x, "d2x")) { output_warn.write( "\tWARNING: differencing quantity 'd2x' not found. Calculating from dx\n"); d1_dx = localmesh->indexDDX(1. / dx); // d/di(1/dx) } else { + // Shift d2x to our location + d2x = interp_to(d2x, location); + d1_dx = -d2x / (dx * dx); } @@ -457,6 +460,9 @@ int Coordinates::geometry() { "\tWARNING: differencing quantity 'd2y' not found. Calculating from dy\n"); d1_dy = localmesh->indexDDY(1. / dy); // d/di(1/dy) } else { + // Shift d2y to our location + d2y = interp_to(d2y, location); + d1_dy = -d2y / (dy * dy); } From 8795acad2b5ec34f53bf815604113348849e502a Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 20 Nov 2018 14:32:15 +0000 Subject: [PATCH 2/7] Remove uses of global 'mesh' in BoutMesh::addBoundaryRegions() Could cause bugs, for example in unit tests of BoutMesh where mesh==nullptr. --- src/mesh/impls/bout/boutmesh.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index e4f3c487a0..e2516fd771 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -2208,7 +2208,7 @@ void BoutMesh::addBoundaryRegions() { all_boundaries.emplace_back("RGN_UPPER_Y"); // Inner X - if(mesh->firstX() && !mesh->periodicX) { + if(firstX() && !periodicX) { addRegion3D("RGN_INNER_X", Region(0, xstart-1, ystart, yend, 0, LocalNz-1, LocalNy, LocalNz, maxregionblocksize)); addRegion2D("RGN_INNER_X", Region(0, xstart-1, ystart, yend, 0, 0, @@ -2225,7 +2225,7 @@ void BoutMesh::addBoundaryRegions() { } // Outer X - if(mesh->firstX() && !mesh->periodicX) { + if(firstX() && !periodicX) { addRegion3D("RGN_OUTER_X", Region(xend+1, LocalNx-1, ystart, yend, 0, LocalNz-1, LocalNy, LocalNz, maxregionblocksize)); addRegion2D("RGN_OUTER_X", Region(xend+1, LocalNx-1, ystart, yend, 0, 0, From 21a0fc548b3221e30869c3a3b17e170f82b5ca5e Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 14 Nov 2018 15:07:23 +0000 Subject: [PATCH 3/7] Allow initialization of Coordinates with Mesh::addCoordinates(location) Instead of creating Coordinates objects when Mesh::getCoordinates is called, create them during initialization. As a temporary workaround to maintain backward compatibility, initialize all locations in BoutMesh::load() so that they are always available. Allow the user to explicitly request a location to be available, by calling Mesh::addCoordinates(location) for each location required. This is optional in v4.2 but will be required from v4.3. When Coordinates::geometry() is called, assume that staggered location Coordinates objects need to be recalculated because the user has changed the CELL_CENTRE version. As a temporary workaround, replace geometry() with geometryNoRecalculate() and geometry() now calls geometryNoRecalculate and then recalculates the staggered versions. This means geometry() can be reset back to not recalculating in next (which is OK because we will add the requirement to call addCoordinates() after calling geometry()) without changing the function signature. --- include/bout/coordinates.hxx | 3 +++ include/bout/mesh.hxx | 25 +++++++++++-------- src/field/field3d.cxx | 6 +++++ src/mesh/coordinates.cxx | 41 +++++++++++++++++++++++++++++--- src/mesh/impls/bout/boutmesh.cxx | 20 ++++++++++++++++ src/mesh/mesh.cxx | 39 ++++++++++++++++++++++++------ 6 files changed, 114 insertions(+), 20 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index e649cdfcc7..624b15b745 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -143,6 +143,9 @@ public: const Field3D Laplace(const Field3D &f, CELL_LOC outloc=CELL_DEFAULT); private: + // temporary work-around method to allow 'geometry' to be called without + // trying to re-calculate the non-CELL_CENTRE fields + int geometryNoRecalculate(); int nz; // Size of mesh in Z. This is mesh->ngz-1 Mesh * localmesh; CELL_LOC location; diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index fcee0e58e3..2c76177640 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -434,21 +434,26 @@ class Mesh { ASSERT1(location != CELL_DEFAULT); ASSERT1(location != CELL_VSHIFT); - if (coords_map.count(location)) { // True branch most common, returns immediately - return coords_map[location].get(); - } else { - // No coordinate system set. Create default - // Note that this can't be allocated here due to incomplete type - // (circular dependency between Mesh and Coordinates) - coords_map.emplace(location, createDefaultCoordinates(location)); - return coords_map[location].get(); +#if CHECK > 0 + if (!coords_map.count(location)) { + throw BoutException("Error: Coordinates for %s have not been added to " + "this Mesh. You should use REQUEST_STAGGER(location) before " + "initializing fields staggered to 'CELL_LOC location'.", + CELL_LOC_STRING(location).c_str()); } +#endif + return coords_map.at(location).get(); } Coordinates *DEPRECATED(coordinates(const CELL_LOC location = CELL_CENTRE)) { return getCoordinates(location); } + /// Add Coordinates object at a certain location. + /// If replace_coords is set to true, reset the object in coords_map if it + /// already exists, otherwise add a new one + void addCoordinates(const CELL_LOC location, bool replace_coords = false); + // First derivatives in index space // Implemented in src/mesh/index_derivs.hxx @@ -712,7 +717,7 @@ class Mesh { GridDataSource *source; ///< Source for grid data - std::map > coords_map; ///< Coordinate systems at different CELL_LOCs + std::map > coords_map; ///< Coordinate systems at different CELL_LOCs Options *options; ///< Mesh options section @@ -751,7 +756,7 @@ class Mesh { private: /// Allocates default Coordinates objects - std::shared_ptr createDefaultCoordinates(const CELL_LOC location); + std::unique_ptr createDefaultCoordinates(const CELL_LOC location); //Internal region related information std::map> regionMap3D; diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 1143eb9fdc..022fd06aca 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -240,6 +240,12 @@ void Field3D::setLocation(CELL_LOC new_location) { } location = new_location; +#if CHECK > 1 + // Check Coordinates for location have been added + // For CHECK > 0, getCoordinates will throw if location has not been added. + fieldmesh->getCoordinates(location); +#endif + // Invalidate the coordinates pointer if (new_location != location) fieldCoordinates = nullptr; diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 8b287638df..11075b37df 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -151,7 +151,7 @@ Coordinates::Coordinates(Mesh *mesh) ////////////////////////////////////////////////////// /// Calculate Christoffel symbols. Needs communication - if (geometry()) { + if (geometryNoRecalculate()) { throw BoutException("Differential geometry failed\n"); } @@ -269,7 +269,7 @@ Coordinates::Coordinates(Mesh *mesh, const CELL_LOC loc, const Coordinates* coor ////////////////////////////////////////////////////// /// Calculate Christoffel symbols. Needs communication - if (geometry()) { + if (geometryNoRecalculate()) { throw BoutException("Differential geometry failed\n"); } @@ -304,7 +304,7 @@ void Coordinates::outputVars(Datafile &file) { file.add(J, "J", false); } -int Coordinates::geometry() { +int Coordinates::geometryNoRecalculate() { TRACE("Coordinates::geometry"); output_progress.write("Calculating differential geometry terms\n"); @@ -469,6 +469,41 @@ int Coordinates::geometry() { return 0; } +int Coordinates::geometry() { + + geometryNoRecalculate(); + + if (location != CELL_CENTRE) { + throw BoutException("geometry() called from a location other than " + "CELL_CENTRE. This is an error as the other Coordinates are calculated " + "from the CELL_CENTRE version, so the changes you have made to this " + "object would be overwritten."); + } + // Coordinates objects at staggered location were calculated from + // CELL_CENTRE ones. geometry() has been called on the CELL_CENTRE + // Coordinates, so they must have changed; we need to re-calculate the + // staggered location Coordinates objects. + + if (localmesh->StaggerGrids) { + // Replace Coordinates objects at staggered locations, if there are + // enough grid points. + // This is a temporary workaround. In v4.3 we will users to call + // REQUEST_LOCATION(location) for each location that is needed and change + // this so that we don't waste memory on unneeded Coordinates. + if (localmesh->LocalNx >= 4) { + localmesh->addCoordinates(CELL_XLOW, true); + } + if (localmesh->LocalNy >= 4) { + localmesh->addCoordinates(CELL_YLOW, true); + } + // Can always add ZLOW Coordinates, since z-interpolation on Field2D is a + // null operation + localmesh->addCoordinates(CELL_ZLOW, true); + } + + return 0; +} + int Coordinates::calcCovariant() { TRACE("Coordinates::calcCovariant"); diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index e2516fd771..62b6a11979 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -839,6 +839,26 @@ int BoutMesh::load() { // Add boundary regions addBoundaryRegions(); + // Create CELL_CENTRE Coordinates object + addCoordinates(CELL_CENTRE); + + if (StaggerGrids) { + // Add Coordinates objects at staggered locations, if there are enough grid + // points. + // This is a temporary workaround. In v4.3 we will users to call + // REQUEST_LOCATION(location) for each location that is needed and remove + // this so that we don't waste memory on unneeded Coordinates. + if (LocalNx >= 4) { + addCoordinates(CELL_XLOW); + } + if (LocalNy >= 4) { + addCoordinates(CELL_YLOW); + } + // Can always add ZLOW Coordinates, since z-interpolation on Field2D is a + // null operation + addCoordinates(CELL_ZLOW); + } + output_info.write("\tdone\n"); return 0; diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx index 6a857f57a7..023beabccc 100644 --- a/src/mesh/mesh.cxx +++ b/src/mesh/mesh.cxx @@ -325,13 +325,38 @@ ParallelTransform& Mesh::getParallelTransform() { return *transform; } -std::shared_ptr Mesh::createDefaultCoordinates(const CELL_LOC location) { - if (location == CELL_CENTRE || location == CELL_DEFAULT) - // Initialize coordinates from input - return std::make_shared(this); - else - // Interpolate coordinates from CELL_CENTRE version - return std::make_shared(this, location, getCoordinates(CELL_CENTRE)); +void Mesh::addCoordinates(const CELL_LOC location, bool replace_coords) { + ASSERT1(location != CELL_DEFAULT); + + if (location == CELL_VSHIFT) { + // CELL_VSHIFT puts vector components at CELL_XLOW, CELL_YLOW and + // CELL_ZLOW, so require Coordinates at all three. + addCoordinates(CELL_XLOW, replace_coords); + addCoordinates(CELL_YLOW, replace_coords); + addCoordinates(CELL_ZLOW, replace_coords); + } else { + // No coordinate system set. Create default + if (location == CELL_CENTRE) { + // Initialize coordinates from input + if (!coords_map.count(location)) { + // location does not exist in coords_map, so create new entry + coords_map.emplace(location, std::unique_ptr(new Coordinates(this))); + } else if (replace_coords) { + // location does already exists in coords_map, so reset it + coords_map.at(location).reset(new Coordinates(this)); + } + } else { + // Interpolate coordinates from CELL_CENTRE version + ASSERT1(StaggerGrids); // If StaggerGrids==false, it doesn't make sense to have non-CELL_CENTRE Coordinates + if (!coords_map.count(location)) { + // location does not exist in coords_map, so create new entry + coords_map.emplace(location, std::unique_ptr(new Coordinates(this, location, getCoordinates(CELL_CENTRE)))); + } else if (replace_coords) { + // location does already exists in coords_map, so reset it + coords_map.at(location).reset(new Coordinates(this, location, getCoordinates(CELL_CENTRE))); + } + } + } } From f528e6aba79be79f13f09cc213c937fa2b4e2e05 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Sun, 18 Nov 2018 20:50:22 +0000 Subject: [PATCH 4/7] Add nullptr entries to coords_map in FakeMesh of unit tests Allows Coordinates pointers to be retrieved where necessary in the unit tests. --- tests/unit/test_extras.hxx | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tests/unit/test_extras.hxx b/tests/unit/test_extras.hxx index 6e6c9dc2ad..b143edef1f 100644 --- a/tests/unit/test_extras.hxx +++ b/tests/unit/test_extras.hxx @@ -73,6 +73,11 @@ public: StaggerGrids = false; IncIntShear = false; maxregionblocksize = MAXREGIONBLOCKSIZE; + + coords_map.emplace(CELL_CENTRE, std::unique_ptr(nullptr)); + coords_map.emplace(CELL_XLOW, std::unique_ptr(nullptr)); + coords_map.emplace(CELL_YLOW, std::unique_ptr(nullptr)); + coords_map.emplace(CELL_ZLOW, std::unique_ptr(nullptr)); } comm_handle send(FieldGroup &UNUSED(g)) { return nullptr; }; From b5b80f24caae6b6c1cc6136c2159444015c2e2c2 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 9 Jan 2019 19:40:43 +0000 Subject: [PATCH 5/7] Fix error message, replacing reference to 'REQUEST_STAGGER' Also similarly fix comments in Coordinates::geometry() and BoutMesh::load(). --- include/bout/mesh.hxx | 4 ++-- src/mesh/coordinates.cxx | 6 +++--- src/mesh/impls/bout/boutmesh.cxx | 6 +++--- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index 2c76177640..d8780deffc 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -437,8 +437,8 @@ class Mesh { #if CHECK > 0 if (!coords_map.count(location)) { throw BoutException("Error: Coordinates for %s have not been added to " - "this Mesh. You should use REQUEST_STAGGER(location) before " - "initializing fields staggered to 'CELL_LOC location'.", + "this Mesh. You should call the method Mesh::addCoordinates(location) " + "before initializing fields staggered to 'CELL_LOC location'.", CELL_LOC_STRING(location).c_str()); } #endif diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 11075b37df..0dd32d2d33 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -487,9 +487,9 @@ int Coordinates::geometry() { if (localmesh->StaggerGrids) { // Replace Coordinates objects at staggered locations, if there are // enough grid points. - // This is a temporary workaround. In v4.3 we will users to call - // REQUEST_LOCATION(location) for each location that is needed and change - // this so that we don't waste memory on unneeded Coordinates. + // This is a temporary workaround. In v4.3 we will require users to call + // Mesh::addCoordinates(location) for each location that is needed and + // change this so that we don't waste memory on unneeded Coordinates. if (localmesh->LocalNx >= 4) { localmesh->addCoordinates(CELL_XLOW, true); } diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index 62b6a11979..54359833aa 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -845,9 +845,9 @@ int BoutMesh::load() { if (StaggerGrids) { // Add Coordinates objects at staggered locations, if there are enough grid // points. - // This is a temporary workaround. In v4.3 we will users to call - // REQUEST_LOCATION(location) for each location that is needed and remove - // this so that we don't waste memory on unneeded Coordinates. + // This is a temporary workaround. In v4.3 we will require users to call + // Mesh::addCoordinates(location) for each location that is needed and + // change this so that we don't waste memory on unneeded Coordinates. if (LocalNx >= 4) { addCoordinates(CELL_XLOW); } From 0bf16afe372c6506ad714602c37b5929c515f545 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 9 Jan 2019 19:42:27 +0000 Subject: [PATCH 6/7] Remove createDefaultCoordinates() method from header Implementation was already removed from Mesh. --- include/bout/mesh.hxx | 3 --- 1 file changed, 3 deletions(-) diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index d8780deffc..705281727f 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -755,9 +755,6 @@ class Mesh { REGION region = RGN_NOBNDRY); private: - /// Allocates default Coordinates objects - std::unique_ptr createDefaultCoordinates(const CELL_LOC location); - //Internal region related information std::map> regionMap3D; std::map> regionMap2D; From 75ce52256d2815fbd674afddfd4434270a045466 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 9 Jan 2019 20:37:05 +0000 Subject: [PATCH 7/7] Revert change of Mesh::coords_map to use std::unique_ptr --- include/bout/mesh.hxx | 2 +- src/mesh/mesh.cxx | 4 ++-- tests/unit/test_extras.hxx | 8 ++++---- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index 705281727f..f5b459d675 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -717,7 +717,7 @@ class Mesh { GridDataSource *source; ///< Source for grid data - std::map > coords_map; ///< Coordinate systems at different CELL_LOCs + std::map > coords_map; ///< Coordinate systems at different CELL_LOCs Options *options; ///< Mesh options section diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx index 023beabccc..6941d10e33 100644 --- a/src/mesh/mesh.cxx +++ b/src/mesh/mesh.cxx @@ -340,7 +340,7 @@ void Mesh::addCoordinates(const CELL_LOC location, bool replace_coords) { // Initialize coordinates from input if (!coords_map.count(location)) { // location does not exist in coords_map, so create new entry - coords_map.emplace(location, std::unique_ptr(new Coordinates(this))); + coords_map.emplace(location, std::make_shared(this)); } else if (replace_coords) { // location does already exists in coords_map, so reset it coords_map.at(location).reset(new Coordinates(this)); @@ -350,7 +350,7 @@ void Mesh::addCoordinates(const CELL_LOC location, bool replace_coords) { ASSERT1(StaggerGrids); // If StaggerGrids==false, it doesn't make sense to have non-CELL_CENTRE Coordinates if (!coords_map.count(location)) { // location does not exist in coords_map, so create new entry - coords_map.emplace(location, std::unique_ptr(new Coordinates(this, location, getCoordinates(CELL_CENTRE)))); + coords_map.emplace(location, std::make_shared(this, location, getCoordinates(CELL_CENTRE))); } else if (replace_coords) { // location does already exists in coords_map, so reset it coords_map.at(location).reset(new Coordinates(this, location, getCoordinates(CELL_CENTRE))); diff --git a/tests/unit/test_extras.hxx b/tests/unit/test_extras.hxx index b143edef1f..b644a3bac8 100644 --- a/tests/unit/test_extras.hxx +++ b/tests/unit/test_extras.hxx @@ -74,10 +74,10 @@ public: IncIntShear = false; maxregionblocksize = MAXREGIONBLOCKSIZE; - coords_map.emplace(CELL_CENTRE, std::unique_ptr(nullptr)); - coords_map.emplace(CELL_XLOW, std::unique_ptr(nullptr)); - coords_map.emplace(CELL_YLOW, std::unique_ptr(nullptr)); - coords_map.emplace(CELL_ZLOW, std::unique_ptr(nullptr)); + coords_map.emplace(CELL_CENTRE, std::shared_ptr(nullptr)); + coords_map.emplace(CELL_XLOW, std::shared_ptr(nullptr)); + coords_map.emplace(CELL_YLOW, std::shared_ptr(nullptr)); + coords_map.emplace(CELL_ZLOW, std::shared_ptr(nullptr)); } comm_handle send(FieldGroup &UNUSED(g)) { return nullptr; };