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..f5b459d675 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 call the method Mesh::addCoordinates(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 @@ -750,9 +755,6 @@ class Mesh { REGION region = RGN_NOBNDRY); private: - /// Allocates default Coordinates objects - std::shared_ptr createDefaultCoordinates(const CELL_LOC location); - //Internal region related information std::map> regionMap3D; std::map> regionMap2D; 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 81c0e921e4..0dd32d2d33 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"); @@ -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,12 +460,50 @@ 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); } 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 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); + } + 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 e4f3c487a0..54359833aa 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 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); + } + 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; @@ -2208,7 +2228,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 +2245,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, diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx index 6a857f57a7..6941d10e33 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::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)); + } + } 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::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 6e6c9dc2ad..b644a3bac8 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::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; };