Boundary conditions fixes

CITATION.cff

@@ -139,11 +139,11 @@ authors:
 
   - family-names: Wang
     given-names: Zhanhui
-version: 4.1.2
-date-released: 2017-12-01
+version: 4.2.0
+date-released: TBC
 repository-code: https://github.com/boutproject/BOUT-dev
 url: http://boutproject.github.io/
-doi: 10.5281/zenodo.1423213
+doi: TBC
 license: 'LGPL-3.0-or-later'
 references:
   - type: article

bin/bout-squashoutput

@@ -9,37 +9,42 @@ from sys import exit
 try:
     import argcomplete
 except ImportError:
-    argcomplete=None
+    argcomplete = None
 import boutdata.squashoutput as squash
 
 # Parse command line arguments
-parser = argparse.ArgumentParser(squash.__doc__+"\n\n"+squash.squashoutput.__doc__)
+parser = argparse.ArgumentParser(
+    squash.__doc__ + "\n\n" + squash.squashoutput.__doc__)
+
 
 def str_to_bool(string):
-    return string.lower()=="true" or string.lower()=="t"
+    return string.lower() == "true" or string.lower() == "t"
+
 
 def int_or_none(string):
     try:
         return int(string)
     except ValueError:
-        if string.lower()=='none' or string.lower()=='n':
+        if string.lower() == 'none' or string.lower() == 'n':
             return None
         else:
             raise
 
 parser.add_argument("datadir", nargs='?', default=".")
-parser.add_argument("--outputname",default="BOUT.dmp.nc")
+parser.add_argument("--outputname", default="BOUT.dmp.nc")
 parser.add_argument("--tind", type=int_or_none, nargs='*', default=[None])
 parser.add_argument("--xind", type=int_or_none, nargs='*', default=[None])
 parser.add_argument("--yind", type=int_or_none, nargs='*', default=[None])
 parser.add_argument("--zind", type=int_or_none, nargs='*', default=[None])
-parser.add_argument("-s","--singleprecision", action="store_true", default=False)
-parser.add_argument("-c","--compress", action="store_true", default=False)
-parser.add_argument("-l","--complevel", type=int_or_none, default=None)
-parser.add_argument("-i","--least-significant-digit", type=int_or_none, default=None)
-parser.add_argument("-q","--quiet", action="store_true", default=False)
-parser.add_argument("-a","--append", action="store_true", default=False)
-parser.add_argument("-d","--delete", action="store_true", default=False)
+parser.add_argument("-s", "--singleprecision",
+                    action="store_true", default=False)
+parser.add_argument("-c", "--compress", action="store_true", default=False)
+parser.add_argument("-l", "--complevel", type=int_or_none, default=None)
+parser.add_argument("-i", "--least-significant-digit",
+                    type=int_or_none, default=None)
+parser.add_argument("-q", "--quiet", action="store_true", default=False)
+parser.add_argument("-a", "--append", action="store_true", default=False)
+parser.add_argument("-d", "--delete", action="store_true", default=False)
 
 if argcomplete:
     argcomplete.autocomplete(parser)
@@ -49,7 +54,7 @@ args = parser.parse_args()
 # Late imports to not slow down bash completion
 
 for ind in "txyz":
-    args.__dict__[ind+"ind"]=slice(*args.__dict__[ind+"ind"])
+    args.__dict__[ind + "ind"] = slice(*args.__dict__[ind + "ind"])
 # Call the function, using command line arguments
 squash.squashoutput(**args.__dict__)
 

configure

@@ -1,6 +1,6 @@
 #! /bin/sh
 # Guess values for system-dependent variables and create Makefiles.
-# Generated by GNU Autoconf 2.69 for BOUT++ 4.1.2.
+# Generated by GNU Autoconf 2.69 for BOUT++ 4.2.0.
 #
 # Report bugs to <bd512@york.ac.uk>.
 #
@@ -580,8 +580,8 @@ MAKEFLAGS=
 # Identity of this package.
 PACKAGE_NAME='BOUT++'
 PACKAGE_TARNAME='bout--'
-PACKAGE_VERSION='4.1.2'
-PACKAGE_STRING='BOUT++ 4.1.2'
+PACKAGE_VERSION='4.2.0'
+PACKAGE_STRING='BOUT++ 4.2.0'
 PACKAGE_BUGREPORT='bd512@york.ac.uk'
 PACKAGE_URL=''
 
@@ -1359,7 +1359,7 @@ if test "$ac_init_help" = "long"; then
   # Omit some internal or obsolete options to make the list less imposing.
   # This message is too long to be a string in the A/UX 3.1 sh.
   cat <<_ACEOF
-\`configure' configures BOUT++ 4.1.2 to adapt to many kinds of systems.
+\`configure' configures BOUT++ 4.2.0 to adapt to many kinds of systems.
 
 Usage: $0 [OPTION]... [VAR=VALUE]...
 
@@ -1421,7 +1421,7 @@ fi
 
 if test -n "$ac_init_help"; then
   case $ac_init_help in
-     short | recursive ) echo "Configuration of BOUT++ 4.1.2:";;
+     short | recursive ) echo "Configuration of BOUT++ 4.2.0:";;
    esac
   cat <<\_ACEOF
 
@@ -1550,7 +1550,7 @@ fi
 test -n "$ac_init_help" && exit $ac_status
 if $ac_init_version; then
   cat <<\_ACEOF
-BOUT++ configure 4.1.2
+BOUT++ configure 4.2.0
 generated by GNU Autoconf 2.69
 
 Copyright (C) 2012 Free Software Foundation, Inc.
@@ -2131,7 +2131,7 @@ cat >config.log <<_ACEOF
 This file contains any messages produced by compilers while
 running configure, to aid debugging if configure makes a mistake.
 
-It was created by BOUT++ $as_me 4.1.2, which was
+It was created by BOUT++ $as_me 4.2.0, which was
 generated by GNU Autoconf 2.69.  Invocation command line was
 
   $ $0 $@
@@ -12705,7 +12705,7 @@ cat >>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1
 # report actual input values of CONFIG_FILES etc. instead of their
 # values after options handling.
 ac_log="
-This file was extended by BOUT++ $as_me 4.1.2, which was
+This file was extended by BOUT++ $as_me 4.2.0, which was
 generated by GNU Autoconf 2.69.  Invocation command line was
 
   CONFIG_FILES    = $CONFIG_FILES
@@ -12758,7 +12758,7 @@ _ACEOF
 cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
 ac_cs_config="`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`"
 ac_cs_version="\\
-BOUT++ config.status 4.1.2
+BOUT++ config.status 4.2.0
 configured by $0, generated by GNU Autoconf 2.69,
   with options \\"\$ac_cs_config\\"
 
@@ -13933,7 +13933,7 @@ cat >>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1
 # report actual input values of CONFIG_FILES etc. instead of their
 # values after options handling.
 ac_log="
-This file was extended by BOUT++ $as_me 4.1.2, which was
+This file was extended by BOUT++ $as_me 4.2.0, which was
 generated by GNU Autoconf 2.69.  Invocation command line was
 
   CONFIG_FILES    = $CONFIG_FILES
@@ -13986,7 +13986,7 @@ _ACEOF
 cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
 ac_cs_config="`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`"
 ac_cs_version="\\
-BOUT++ config.status 4.1.2
+BOUT++ config.status 4.2.0
 configured by $0, generated by GNU Autoconf 2.69,
   with options \\"\$ac_cs_config\\"
 

configure.ac

@@ -32,7 +32,7 @@
 #
 
 AC_PREREQ([2.69])
-AC_INIT([BOUT++],[4.1.2],[bd512@york.ac.uk])
+AC_INIT([BOUT++],[4.2.0],[bd512@york.ac.uk])
 AC_CONFIG_AUX_DIR([build-aux])
 AC_CONFIG_MACRO_DIR([m4])
 

examples/bout_runners_example/diffusion_3D.cxx

@@ -60,8 +60,8 @@ int physics_init(bool restarting) {
         // The boundary lies (1/2)*dx away from the last point As there
         // are 2 boundaries there will effectively add one more line
         // segment in the domain. Hence
-        mesh->coordinates()->dx = Lx/(internal_x_points);
-        mesh->coordinates()->dy = Ly/(internal_y_points);
+        mesh->getCoordinates()->dx = Lx/(internal_x_points);
+        mesh->getCoordinates()->dy = Ly/(internal_y_points);
     }
     // ************************************************************************
 

examples/laplacexy/alfven-wave/alfven.cxx

@@ -170,7 +170,7 @@ class Alfven : public PhysicsModel {
     Field2D Rxy, Bpxy, Btxy, hthe, sinty;
     GRID_LOAD5(Rxy, Bpxy, Btxy, hthe, sinty); // Load metrics
 
-    Coordinates *coord = mesh->coordinates(); // Metric tensor
+    Coordinates *coord = mesh->getCoordinates(); // Metric tensor
 
     // Checking for dpsi and qinty used in BOUT grids
     Field2D dx;

examples/laplacexy/laplace_perp/test.cxx

@@ -20,7 +20,7 @@ int main(int argc, char** argv) {
     mesh->get(hthe, "hthe"); // m
     mesh->get(I,    "sinty");// m^-2 T^-1
 
-    Coordinates *coord = mesh->coordinates();
+    Coordinates *coord = mesh->getCoordinates();
 
     // Calculate metrics
     coord->g11 = SQ(Rxy * Bpxy);

examples/orszag-tang/mhd.cxx

@@ -47,7 +47,7 @@ class MHD : public PhysicsModel {
     B.covariant = false; // evolve contravariant components
     bout_solve(B, "B");
 
-    Coordinates *coord = mesh->coordinates();
+    Coordinates *coord = mesh->getCoordinates();
     output.write("dx[0,0] = %e, dy[0,0] = %e, dz = %e\n", coord->dx(0, 0),
                  coord->dy(0, 0), coord->dz);
 

examples/tokamak-2fluid/2fluid.cxx

@@ -101,7 +101,7 @@ class TwoFluid : public PhysicsModel {
     Field2D I; // Shear factor 
 
     // Get the coordinate system
-    coord = mesh->coordinates();
+    coord = mesh->getCoordinates();
 
     output.write("Solving 6-variable 2-fluid equations\n");
 

include/bout/coordinates.hxx

@@ -123,8 +123,8 @@ public:
   const Field3D Div_par(const Field3D &f, CELL_LOC outloc=CELL_DEFAULT, DIFF_METHOD method=DIFF_DEFAULT);
 
   // Second derivative along magnetic field
-  const Field2D Grad2_par2(const Field2D &f, CELL_LOC outloc=CELL_DEFAULT);
-  const Field3D Grad2_par2(const Field3D &f, CELL_LOC outloc=CELL_DEFAULT);
+  const Field2D Grad2_par2(const Field2D &f, CELL_LOC outloc=CELL_DEFAULT, DIFF_METHOD method=DIFF_DEFAULT);
+  const Field3D Grad2_par2(const Field3D &f, CELL_LOC outloc=CELL_DEFAULT, DIFF_METHOD method=DIFF_DEFAULT);
 
   // Perpendicular Laplacian operator, using only X-Z derivatives
   // NOTE: This might be better bundled with the Laplacian inversion code

include/bout/mesh.hxx

@@ -430,20 +430,25 @@ class Mesh {
   bool IncIntShear; ///< Include integrated shear (if shifting X)
 
   /// Coordinate system
-  Coordinates *coordinates(const CELL_LOC location = CELL_CENTRE) {
+  Coordinates *getCoordinates(const CELL_LOC location = CELL_CENTRE) {
+    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 if (location == CELL_DEFAULT) {
-      throw BoutException("Ambiguous location 'CELL_DEFAULT' passed to mesh::coordinates");
     } 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.insert(std::pair<CELL_LOC, std::shared_ptr<Coordinates> >(location, createDefaultCoordinates(location)));
+      coords_map.emplace(location, createDefaultCoordinates(location));
       return coords_map[location].get();
     }
   }
 
+  Coordinates *DEPRECATED(coordinates(const CELL_LOC location = CELL_CENTRE)) {
+    return getCoordinates(location);
+  }
+
   // First derivatives in index space
   // Implemented in src/mesh/index_derivs.hxx
 
@@ -650,12 +655,12 @@ class Mesh {
   /// Get the named region from the region_map for the data iterator
   ///
   /// Throws if region_name not found
-  Region<> &getRegion(const std::string &region_name){
+  const Region<> &getRegion(const std::string &region_name) const{
     return getRegion3D(region_name);
   }
-  Region<Ind3D> &getRegion3D(const std::string &region_name);
-  Region<Ind2D> &getRegion2D(const std::string &region_name);
-  Region<IndPerp> &getRegionPerp(const std::string &region_name);
+  const Region<Ind3D> &getRegion3D(const std::string &region_name) const;
+  const Region<Ind2D> &getRegion2D(const std::string &region_name) const;
+  const Region<IndPerp> &getRegionPerp(const std::string &region_name) const;
 
   /// Add a new region to the region_map for the data iterator
   ///

include/bout/region.hxx

@@ -483,8 +483,10 @@ public:
   /// Note that if the indices are altered using these iterators, the
   /// blocks may become out of sync and will need to manually updated
   typename RegionIndices::iterator begin() { return std::begin(indices); };
+  typename RegionIndices::const_iterator begin() const { return std::begin(indices); };
   typename RegionIndices::const_iterator cbegin() const { return indices.cbegin(); };
   typename RegionIndices::iterator end() { return std::end(indices); };
+  typename RegionIndices::const_iterator end() const { return std::end(indices); };
   typename RegionIndices::const_iterator cend() const { return indices.cend(); };
 
   const ContiguousBlocks &getBlocks() const { return blocks; };

include/difops.hxx

@@ -173,7 +173,7 @@ const Field3D Div_par(const Field3D &f, const Field3D &v);
  *
  * Note: For parallel Laplacian use LaplacePar
  */
-const Field2D Grad2_par2(const Field2D &f, CELL_LOC outloc=CELL_DEFAULT);
+const Field2D Grad2_par2(const Field2D &f, CELL_LOC outloc=CELL_DEFAULT, DIFF_METHOD method=DIFF_DEFAULT);
 
 /*!
  * second parallel derivative
@@ -186,7 +186,7 @@ const Field2D Grad2_par2(const Field2D &f, CELL_LOC outloc=CELL_DEFAULT);
  * @param[in] f The field to be differentiated
  * @param[in] outloc The cell location of the result 
  */
-const Field3D Grad2_par2(const Field3D &f, CELL_LOC outloc=CELL_DEFAULT);
+const Field3D Grad2_par2(const Field3D &f, CELL_LOC outloc=CELL_DEFAULT, DIFF_METHOD method=DIFF_DEFAULT);
 
 /*!
  * Parallel derivatives, converting between cell-centred and lower cell boundary

include/field2d.hxx

@@ -154,6 +154,10 @@ class Field2D : public Field, public FieldData {
    */
   const IndexRange DEPRECATED(region(REGION rgn)) const override;
 
+  /// Return a Region<Ind2D> reference to use to iterate over this field
+  const Region<Ind2D>& getRegion(REGION region) const;  
+  const Region<Ind2D>& getRegion(const std::string &region_name) const;
+
   BoutReal& operator[](const Ind2D &d) {
     return data[d.ind];
   }

include/field3d.hxx

@@ -328,6 +328,10 @@ class Field3D : public Field, public FieldData {
    */
   const IndexRange DEPRECATED(region2D(REGION rgn)) const;
 
+  /// Return a Region<Ind3D> reference to use to iterate over this field
+  const Region<Ind3D>& getRegion(REGION region) const;  
+  const Region<Ind3D>& getRegion(const std::string &region_name) const;
+
   /*!
    * Direct data access using DataIterator object.
    * This uses operator(x,y,z) so checks will only be

include/field_data.hxx

@@ -62,12 +62,23 @@ class FieldVisitor;
 */
 class FieldData {
 public:
-  FieldData();
+  FieldData(Mesh *datamesh = nullptr)
+      : fielddatamesh(datamesh != nullptr ? datamesh : mesh), boundaryIsCopy(false),
+        boundaryIsSet(true) {}
+
   virtual ~FieldData();
 
   // Visitor pattern support
   virtual void accept(FieldVisitor &v) = 0;
-
+
+  Mesh *getDataMesh() const {
+    if (fielddatamesh != nullptr) {
+      return fielddatamesh;
+    } else {
+      return mesh;
+    }
+  }
+
   // Defines interface which must be implemented
   virtual bool isReal() const = 0; ///< Returns true if field consists of BoutReal values
   virtual bool is3D() const = 0;   ///< True if variable is 3D
@@ -92,6 +103,7 @@ public:
   FieldGeneratorPtr getBndryGenerator(BndryLoc location);
 
 protected:
+  Mesh* fielddatamesh;
   vector<BoundaryOp *> bndry_op; ///< Boundary conditions
   bool boundaryIsCopy;           ///< True if bndry_op is a copy
   bool boundaryIsSet;            ///< Set to true when setBoundary called

include/fieldperp.hxx

@@ -91,6 +91,10 @@ class FieldPerp : public Field {
 
   const IndexRange DEPRECATED(region(REGION rgn)) const override;
 
+  /// Return a Region<IndPerp> reference to use to iterate over this field
+  const Region<IndPerp>& getRegion(REGION region) const;  
+  const Region<IndPerp>& getRegion(const std::string &region_name) const;
+
   /*!
    * Direct data access using DataIterator indexing
    */