Test difops

.travis.yml

@@ -52,7 +52,7 @@ matrix:
       - *default_env
       - CONFIGURE_OPTIONS='--enable-shared'
       - SCRIPT_FLAGS="-uim -t python -t shared"
-      - PIP_PACKAGES='netcdf4 cython'
+      - PIP_PACKAGES='netcdf4 cython sympy'
   - env:
       - *default_env
       - CONFIGURE_OPTIONS='--enable-openmp'

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])
 

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/griddata.hxx

@@ -63,6 +63,10 @@ public:
                    Direction dir = GridDataSource::X) = 0;
   virtual bool get(Mesh *m, vector<BoutReal> &var, const string &name, int len,
                    int offset = 0, Direction dir = GridDataSource::X) = 0;
+
+  /// Test if grid data source includes y-boundary guard cells.
+  /// Older grid files may not: this then requires some special handling.
+  virtual bool hasYGuards() { return true; }
 };
 
 /// Interface to grid data in a file
@@ -89,6 +93,14 @@ public:
   bool get(Mesh *m, vector<BoutReal> &var, const string &name, int len, int offset = 0,
            GridDataSource::Direction dir = GridDataSource::X) override;
 
+  // Grid files from hypnotoad don't have y-guard cells at all. But even if
+  // they did, they don't have guard cells for the y-boundaries of the core
+  // region because the grid structure means the grid cells adjacent to the
+  // first/last core cells in the input arrays in the file store data for the
+  // PF region. So until a GridFile supports input with different domains, it
+  // cannot contain *all* y-guard cells
+  bool hasYGuards() override { return false; }
+
 private:
   std::unique_ptr<DataFormat> file;
   string filename;

include/bout/mesh.hxx

@@ -151,10 +151,10 @@ class Mesh {
   /// @param[out] var   This will be set to the value. Will be allocated if needed
   /// @param[in] name   Name of the variable to read
   /// @param[in] def    The default value if not found
-  /// @param[in] communicate  Should the field be communicated to fill guard cells?
-  ///
+  /// @param[in] allow_communicate  Allow the field to be communicated if
+  ///                               necessary to fill guard cells
   /// @returns zero if successful, non-zero on failure
-  int get(Field3D &var, const string &name, BoutReal def=0.0, bool communicate=true);
+  int get(Field3D &var, const string &name, BoutReal def=0.0, bool allow_communicate=true);
 
   /// Get a Vector2D from the input source.
   /// If \p var is covariant then this gets three
@@ -650,12 +650,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

@@ -79,7 +79,8 @@ const Field3D Grad_parP(const Field3D &apar, const Field3D &f);
  * @param[in] f  The scalar field to be differentiated
  * 
  */
-const Field2D Vpar_Grad_par(const Field2D &v, const Field2D &f);
+const Field2D Vpar_Grad_par(const Field2D &v, const Field2D &f,
+			    CELL_LOC outloc=CELL_DEFAULT, DIFF_METHOD method=DIFF_DEFAULT);
 
 /*!
  * vpar times parallel derivative along unperturbed B-field (upwinding)
@@ -121,9 +122,12 @@ const Field3D Vpar_Grad_par(const Field3D &v, const Field3D &f, DIFF_METHOD meth
  * \f]
  * 
  * @param[in] f  The component of a vector along the magnetic field 
+ * @param[in] outloc  The cell location for the result. By default the same as \p f
+ * @param[in] method  The numerical method to use
  * 
  */
-const Field2D Div_par(const Field2D &f);
+const Field2D Div_par(const Field2D &f,
+    CELL_LOC outloc=CELL_DEFAULT, DIFF_METHOD method=DIFF_DEFAULT);
 
 /*!
  * parallel divergence operator
@@ -158,6 +162,8 @@ const Field3D Div_par(const Field3D &f, DIFF_METHOD method, CELL_LOC outloc = CE
 const Field3D Div_par_flux(const Field3D &v, const Field3D &f, 
 		      CELL_LOC outloc=CELL_DEFAULT, DIFF_METHOD method=DIFF_DEFAULT);
 const Field3D Div_par_flux(const Field3D &v, const Field3D &f, DIFF_METHOD method, CELL_LOC outloc = CELL_DEFAULT);
+const Field2D Div_par_flux(const Field2D &v, const Field2D &f,
+		      CELL_LOC outloc=CELL_DEFAULT, DIFF_METHOD method=DIFF_DEFAULT);
 
 // Divergence of a parallel flow: Div(f*v)
 // Both f and v are interpolated onto cell boundaries
@@ -173,7 +179,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 +192,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
@@ -214,10 +220,10 @@ const Field2D Div_par_CtoL(const Field2D &var);
  */
 const Field2D Div_par_K_Grad_par(BoutReal kY, const Field2D &f, CELL_LOC outloc=CELL_DEFAULT);
 const Field3D Div_par_K_Grad_par(BoutReal kY, const Field3D &f, CELL_LOC outloc=CELL_DEFAULT);
-const Field2D Div_par_K_Grad_par(const Field2D &kY, const Field2D &f, CELL_LOC outloc=CELL_DEFAULT);
-const Field3D Div_par_K_Grad_par(const Field2D &kY, const Field3D &f, CELL_LOC outloc=CELL_DEFAULT);
-const Field3D Div_par_K_Grad_par(const Field3D &kY, const Field2D &f, CELL_LOC outloc=CELL_DEFAULT);
-const Field3D Div_par_K_Grad_par(const Field3D &kY, const Field3D &f, CELL_LOC outloc=CELL_DEFAULT);
+const Field2D Div_par_K_Grad_par(const Field2D &kY, const Field2D &f, CELL_LOC outloc=CELL_DEFAULT, DIFF_METHOD method=DIFF_DEFAULT);
+const Field3D Div_par_K_Grad_par(const Field2D &kY, const Field3D &f, CELL_LOC outloc=CELL_DEFAULT, DIFF_METHOD method=DIFF_DEFAULT);
+const Field3D Div_par_K_Grad_par(const Field3D &kY, const Field2D &f, CELL_LOC outloc=CELL_DEFAULT, DIFF_METHOD method=DIFF_DEFAULT);
+const Field3D Div_par_K_Grad_par(const Field3D &kY, const Field3D &f, CELL_LOC outloc=CELL_DEFAULT, DIFF_METHOD method=DIFF_DEFAULT);
 
 /*!
  * Perpendicular Laplacian operator

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,20 @@ class FieldVisitor;
 */
 class FieldData {
 public:
-  FieldData();
+  FieldData(Mesh* m);
   virtual ~FieldData();
 
   // Visitor pattern support
   virtual void accept(FieldVisitor &v) = 0;
 
+  virtual Mesh * getDataMesh() const{
+    if (fielddatamesh){
+      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 +100,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
    */