diff --git a/.travis.yml b/.travis.yml index 70fd3f23ad..eadb4dc225 100644 --- a/.travis.yml +++ b/.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' diff --git a/CITATION.cff b/CITATION.cff index 345fd55d06..ab96c9484b 100644 --- a/CITATION.cff +++ b/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 diff --git a/bin/bout-squashoutput b/bin/bout-squashoutput index cfc6154710..cd3ccf5f92 100755 --- a/bin/bout-squashoutput +++ b/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__) diff --git a/configure b/configure index 81ee0d7f7a..dba763cd86 100755 --- a/configure +++ b/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 . # @@ -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\\" diff --git a/configure.ac b/configure.ac index 7cd03e2948..21dd288bbc 100644 --- a/configure.ac +++ b/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]) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index ba9522c6c3..e649cdfcc7 100644 --- a/include/bout/coordinates.hxx +++ b/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 diff --git a/include/bout/griddata.hxx b/include/bout/griddata.hxx index afc661303b..dd4d6e7907 100644 --- a/include/bout/griddata.hxx +++ b/include/bout/griddata.hxx @@ -63,6 +63,10 @@ public: Direction dir = GridDataSource::X) = 0; virtual bool get(Mesh *m, vector &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 &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 file; string filename; diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index f3682c9dbd..72d315ba17 100644 --- a/include/bout/mesh.hxx +++ b/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 ®ion_name){ + const Region<> &getRegion(const std::string ®ion_name) const{ return getRegion3D(region_name); } - Region &getRegion3D(const std::string ®ion_name); - Region &getRegion2D(const std::string ®ion_name); - Region &getRegionPerp(const std::string ®ion_name); + const Region &getRegion3D(const std::string ®ion_name) const; + const Region &getRegion2D(const std::string ®ion_name) const; + const Region &getRegionPerp(const std::string ®ion_name) const; /// Add a new region to the region_map for the data iterator /// diff --git a/include/bout/region.hxx b/include/bout/region.hxx index d11587ef98..d4671ebe25 100644 --- a/include/bout/region.hxx +++ b/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; }; diff --git a/include/difops.hxx b/include/difops.hxx index 029fade7eb..c7eeb4fd8f 100644 --- a/include/difops.hxx +++ b/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 diff --git a/include/field2d.hxx b/include/field2d.hxx index b261ca8456..51e14607ab 100644 --- a/include/field2d.hxx +++ b/include/field2d.hxx @@ -154,6 +154,10 @@ class Field2D : public Field, public FieldData { */ const IndexRange DEPRECATED(region(REGION rgn)) const override; + /// Return a Region reference to use to iterate over this field + const Region& getRegion(REGION region) const; + const Region& getRegion(const std::string ®ion_name) const; + BoutReal& operator[](const Ind2D &d) { return data[d.ind]; } diff --git a/include/field3d.hxx b/include/field3d.hxx index bd085f7c73..1025f935ab 100644 --- a/include/field3d.hxx +++ b/include/field3d.hxx @@ -328,6 +328,10 @@ class Field3D : public Field, public FieldData { */ const IndexRange DEPRECATED(region2D(REGION rgn)) const; + /// Return a Region reference to use to iterate over this field + const Region& getRegion(REGION region) const; + const Region& getRegion(const std::string ®ion_name) const; + /*! * Direct data access using DataIterator object. * This uses operator(x,y,z) so checks will only be diff --git a/include/field_data.hxx b/include/field_data.hxx index c4fadd4e95..d75e89c2ca 100644 --- a/include/field_data.hxx +++ b/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 bndry_op; ///< Boundary conditions bool boundaryIsCopy; ///< True if bndry_op is a copy bool boundaryIsSet; ///< Set to true when setBoundary called diff --git a/include/fieldperp.hxx b/include/fieldperp.hxx index 8fef875e8a..90fdcf03d8 100644 --- a/include/fieldperp.hxx +++ b/include/fieldperp.hxx @@ -91,6 +91,10 @@ class FieldPerp : public Field { const IndexRange DEPRECATED(region(REGION rgn)) const override; + /// Return a Region reference to use to iterate over this field + const Region& getRegion(REGION region) const; + const Region& getRegion(const std::string ®ion_name) const; + /*! * Direct data access using DataIterator indexing */ diff --git a/manual/doxygen/Doxyfile b/manual/doxygen/Doxyfile index d7067d0171..3271607d7a 100644 --- a/manual/doxygen/Doxyfile +++ b/manual/doxygen/Doxyfile @@ -38,7 +38,7 @@ PROJECT_NAME = BOUT++ # could be handy for archiving the generated documentation or if some version # control system is used. -PROJECT_NUMBER = 4.1.2 +PROJECT_NUMBER = 4.2.0 # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewer a diff --git a/manual/doxygen/Doxyfile_readthedocs b/manual/doxygen/Doxyfile_readthedocs index 4249ec2ebd..7f3dfd8a94 100644 --- a/manual/doxygen/Doxyfile_readthedocs +++ b/manual/doxygen/Doxyfile_readthedocs @@ -38,7 +38,7 @@ PROJECT_NAME = BOUT++ # could be handy for archiving the generated documentation or if some version # control system is used. -PROJECT_NUMBER = 4.1.2 +PROJECT_NUMBER = 4.2.0 # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewer a diff --git a/manual/sphinx/conf.py b/manual/sphinx/conf.py index 5e514e3283..9c39c2ac27 100755 --- a/manual/sphinx/conf.py +++ b/manual/sphinx/conf.py @@ -131,9 +131,9 @@ def __getattr__(cls, name): # built documents. # # The short X.Y version. -version = '4.1' +version = '4.2' # The full version, including alpha/beta/rc tags. -release = '4.1.2' +release = '4.2.0' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/manual/sphinx/developer_docs/data_types.rst b/manual/sphinx/developer_docs/data_types.rst index aec326d7f3..6bfffff328 100644 --- a/manual/sphinx/developer_docs/data_types.rst +++ b/manual/sphinx/developer_docs/data_types.rst @@ -265,7 +265,7 @@ to OpenMP parallelise or vectorise:: } If you wish to vectorise but can't use OpenMP then there is a serial -verion of the macro: +verion of the macro:: BoutReal max=0.; BOUT_FOR_SERIAL(i, region) { @@ -285,10 +285,10 @@ For loops inside parallel regions, there is ``BOUT_FOR_INNER``:: If a more general OpenMP directive is needed, there is ``BOUT_FOR_OMP``:: - BoutReal result=0.; - BOUT_FOR_OMP(i, region, parallel for reduction(max:result)) { - result = f[i] > result ? f[i] : result; - } + BoutReal result=0.; + BOUT_FOR_OMP(i, region, parallel for reduction(max:result)) { + result = f[i] > result ? f[i] : result; + } The iterator provides access to the x, y, z indices:: @@ -304,11 +304,11 @@ modulo operators are needed to calculate individual indices. To perform finite difference or similar operators, index offsets can be calculated:: - Field3D f = ...; - Field3D g(0.0); - BOUT_FOR(i, f.getMesh()->getRegion3D("RGN_NOBNDRY")) { - g[i] = f[i.xp()] - f[i.xm()]; - } + Field3D f = ...; + Field3D g(0.0); + BOUT_FOR(i, f.getMesh()->getRegion3D("RGN_NOBNDRY")) { + g[i] = f[i.xp()] - f[i.xm()]; + } The ``xp()`` function by default produces an offset of ``+1`` in ``X``, ``xm()`` an offset of ``-1`` in the ``X`` direction. These functions can also diff --git a/manual/sphinx/user_docs/advanced_install.rst b/manual/sphinx/user_docs/advanced_install.rst index ee6e80ecde..6468f31d47 100644 --- a/manual/sphinx/user_docs/advanced_install.rst +++ b/manual/sphinx/user_docs/advanced_install.rst @@ -53,7 +53,7 @@ control over how BOUT++ is built: - ``SUNDIALS_EXTRA_LIBS`` specifies additional libraries for linking to SUNDIALS, which are put at the end of the link command. - + It is possible to change flags for BOUT++ after running configure, by editing the ``make.config`` file. Note that this is not recommended, as e.g. PVODE will not be built with these flags. @@ -165,13 +165,14 @@ To compile for the SKL partition, configure with to enable AVX512 vectorization. -.. note:: As of 20/04/2018, an issue with the netcdf and netcdf-cxx4 modules - means that you will need to remove ``-lnetcdf`` from ``EXTRA_LIBS`` in - ``make.config`` after running ``./configure`` and before running - ``make``. ``-lnetcdf`` needs also to be removed from ``bin/bout-config`` - to allow a successful build of the python interface. Recreation of - ``boutcore.pyx`` needs to be manually triggered, if - ``boutcore.pyx`` has already been created. +.. note:: As of 20/04/2018, an issue with the netcdf and netcdf-cxx4 + modules means that you will need to remove ``-lnetcdf`` from + ``EXTRA_LIBS`` in ``make.config`` after running + ``./configure`` and before running ``make``. ``-lnetcdf`` + needs also to be removed from ``bin/bout-config`` to allow a + successful build of the python interface. Recreation of + ``boutcore.pyx`` needs to be manually triggered, if + ``boutcore.pyx`` has already been created. Ubgl ~~~~ @@ -258,12 +259,38 @@ appropriately. OpenMP ------ -BOUT++ can make use of Single-Instruction Multiple-Data (SIMD) -parallelism through OpenMP. To enable OpenMP, use the +BOUT++ can make use of OpenMP parallelism. To enable OpenMP, use the ``--enable-openmp`` flag to configure:: ./configure --enable-openmp +OpenMP can be used to parallelise in more directions than can be +achieved with MPI alone. For example, it is currently difficult to +parallelise in X using pure MPI if FCI is used, and impossible to +parallelise at all in Z with pure MPI. + +OpenMP is in a large number of places now, such that a decent speed-up +can be achieved with OpenMP alone. Hybrid parallelisation with both +MPI and OpenMP can lead to more significant speed-ups, but it +sometimes requires some fine tuning of numerical parameters in order +to achieve this. This greatly depends on the details not just of your +system, but also your particular problem. We have tried to choose +"sensible" defaults that will work well for the most common cases, but +this is not always possible. You may need to perform some testing +yourself to find e.g. the optimum split of OpenMP threads and MPI +ranks. + +One such parameter that can potentially have a significant effect (for +some problem sizes on some machines) is setting the OpenMP schedule +used in some of the OpenMP loops (specifically those using +`BOUT_FOR`). This can be set using:: + + ./configure --enable-openmp --with-openmp-schedule= + +with ```` being one of: ``static`` (the default), +``dynamic``, ``guided``, ``auto`` or ``runtime``. + + .. note:: If you want to use OpenMP with Clang, you will need Clang 3.7+, and either ``libomp`` or ``libiomp``. @@ -277,6 +304,14 @@ parallelism through OpenMP. To enable OpenMP, use the By default PVODE is built without OpenMP support. To enable this add ``--enable-pvode-openmp`` to the configure command. + +.. note:: + OpenMP will attempt to use all available threads by default. This + can cause oversubscription problems on certain systems. You can + limit the number of threads OpenMP uses with the + ``OMP_NUM_THREADS`` environment variable. See your system + documentation for more details. + .. _sec-sundials: SUNDIALS diff --git a/src/field/field2d.cxx b/src/field/field2d.cxx index c7859f1d46..ad0885f2ab 100644 --- a/src/field/field2d.cxx +++ b/src/field/field2d.cxx @@ -45,7 +45,8 @@ #include -Field2D::Field2D(Mesh *localmesh) : Field(localmesh), deriv(nullptr) { +Field2D::Field2D(Mesh *localmesh) : + Field(localmesh), FieldData(localmesh), deriv(nullptr) { boundaryIsSet = false; @@ -66,6 +67,7 @@ Field2D::Field2D(Mesh *localmesh) : Field(localmesh), deriv(nullptr) { } Field2D::Field2D(const Field2D& f) : Field(f.fieldmesh), // The mesh containing array sizes + FieldData(f.fieldmesh), data(f.data), // This handles references to the data array deriv(nullptr) { TRACE("Field2D(Field2D&)"); @@ -95,7 +97,9 @@ Field2D::Field2D(const Field2D& f) : Field(f.fieldmesh), // The mesh containing boundaryIsSet = false; } -Field2D::Field2D(BoutReal val, Mesh *localmesh) : Field(localmesh), deriv(nullptr) { +Field2D::Field2D(BoutReal val, Mesh *localmesh) : + Field(localmesh), FieldData(localmesh), deriv(nullptr) { + boundaryIsSet = false; nx = fieldmesh->LocalNx; @@ -179,6 +183,13 @@ const IndexRange Field2D::region(REGION rgn) const { }; } +const Region &Field2D::getRegion(REGION region) const { + return fieldmesh->getRegion2D(REGION_STRING(region)); +}; +const Region &Field2D::getRegion(const std::string ®ion_name) const { + return fieldmesh->getRegion2D(region_name); +}; + void Field2D::setLocation(CELL_LOC new_location) { if (getMesh()->StaggerGrids) { if (new_location == CELL_VSHIFT) { diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index c62900285a..7a141d2a7f 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -45,8 +45,8 @@ /// Constructor Field3D::Field3D(Mesh *localmesh) - : Field(localmesh), background(nullptr), deriv(nullptr), yup_field(nullptr), - ydown_field(nullptr) { + : Field(localmesh), FieldData(localmesh), background(nullptr), + deriv(nullptr), yup_field(nullptr), ydown_field(nullptr) { #ifdef TRACK name = ""; #endif @@ -71,6 +71,7 @@ Field3D::Field3D(Mesh *localmesh) /// later) Field3D::Field3D(const Field3D &f) : Field(f.fieldmesh), // The mesh containing array sizes + FieldData(f.fieldmesh), background(nullptr), data(f.data), // This handles references to the data array deriv(nullptr), yup_field(nullptr), ydown_field(nullptr) { @@ -100,8 +101,8 @@ Field3D::Field3D(const Field3D &f) } Field3D::Field3D(const Field2D &f) - : Field(f.getMesh()), background(nullptr), deriv(nullptr), yup_field(nullptr), - ydown_field(nullptr) { + : Field(f.getMesh()), FieldData(f.getMesh()), background(nullptr), + deriv(nullptr), yup_field(nullptr), ydown_field(nullptr) { TRACE("Field3D: Copy constructor from Field2D"); @@ -118,8 +119,8 @@ Field3D::Field3D(const Field2D &f) } Field3D::Field3D(const BoutReal val, Mesh *localmesh) - : Field(localmesh), background(nullptr), deriv(nullptr), yup_field(nullptr), - ydown_field(nullptr) { + : Field(localmesh), FieldData(localmesh), background(nullptr), + deriv(nullptr), yup_field(nullptr), ydown_field(nullptr) { TRACE("Field3D: Copy constructor from value"); @@ -355,6 +356,13 @@ const IndexRange Field3D::region2D(REGION rgn) const { }; } +const Region &Field3D::getRegion(REGION region) const { + return fieldmesh->getRegion3D(REGION_STRING(region)); +}; +const Region &Field3D::getRegion(const std::string ®ion_name) const { + return fieldmesh->getRegion3D(region_name); +}; + /////////////////// ASSIGNMENT //////////////////// Field3D & Field3D::operator=(const Field3D &rhs) { diff --git a/src/field/field_data.cxx b/src/field/field_data.cxx index 6ab19709c7..48c48e4847 100644 --- a/src/field/field_data.cxx +++ b/src/field/field_data.cxx @@ -6,8 +6,11 @@ #include #include "unused.hxx" -FieldData::FieldData() : boundaryIsCopy(false), boundaryIsSet(true) { - +FieldData::FieldData(Mesh* m) : + fielddatamesh(m), boundaryIsCopy(false), boundaryIsSet(true) { + if (fielddatamesh == nullptr) { + fielddatamesh = mesh; + } } FieldData::~FieldData() { @@ -24,7 +27,7 @@ void FieldData::setBoundary(const string &name) { output_info << "Setting boundary for variable " << name << endl; /// Loop over the mesh boundary regions - for(const auto& reg : mesh->getBoundaries()) { + for(const auto& reg : getDataMesh()->getBoundaries()) { BoundaryOp* op = static_cast(bfact->createFromOptions(name, reg)); if (op != nullptr) bndry_op.push_back(op); @@ -32,9 +35,9 @@ void FieldData::setBoundary(const string &name) { } /// Get the mesh boundary regions - vector par_reg = mesh->getBoundariesPar(); + vector par_reg = getDataMesh()->getBoundariesPar(); /// Loop over the mesh parallel boundary regions - for(const auto& reg : mesh->getBoundariesPar()) { + for(const auto& reg : getDataMesh()->getBoundariesPar()) { BoundaryOpPar* op = static_cast(bfact->createFromOptions(name, reg)); if (op != nullptr) bndry_op_par.push_back(op); @@ -47,7 +50,7 @@ void FieldData::setBoundary(const string &name) { void FieldData::setBoundary(const string &UNUSED(region), BoundaryOp *op) { /// Get the mesh boundary regions - vector reg = mesh->getBoundaries(); + vector reg = getDataMesh()->getBoundaries(); /// Find the region @@ -76,7 +79,7 @@ void FieldData::addBndryFunction(FuncPtr userfunc, BndryLoc location){ void FieldData::addBndryGenerator(FieldGeneratorPtr gen, BndryLoc location) { if(location == BNDRY_ALL){ - for(const auto& reg : mesh->getBoundaries()) { + for(const auto& reg : getDataMesh()->getBoundaries()) { bndry_generator[reg->location] = gen; } } else { diff --git a/src/field/fieldperp.cxx b/src/field/fieldperp.cxx index 856b6c0272..5bd7765e91 100644 --- a/src/field/fieldperp.cxx +++ b/src/field/fieldperp.cxx @@ -215,6 +215,13 @@ const IndexRange FieldPerp::region(REGION rgn) const { }; } +const Region &FieldPerp::getRegion(REGION region) const { + return fieldmesh->getRegionPerp(REGION_STRING(region)); +}; +const Region &FieldPerp::getRegion(const std::string ®ion_name) const { + return fieldmesh->getRegionPerp(region_name); +}; + //////////////// NON-MEMBER FUNCTIONS ////////////////// ////////////// NON-MEMBER OVERLOADED OPERATORS ////////////// diff --git a/src/field/vector2d.cxx b/src/field/vector2d.cxx index e94e26b512..e1aa7d50ec 100644 --- a/src/field/vector2d.cxx +++ b/src/field/vector2d.cxx @@ -36,10 +36,12 @@ #include Vector2D::Vector2D(Mesh *localmesh) - : x(localmesh), y(localmesh), z(localmesh), covariant(true), deriv(nullptr), location(CELL_CENTRE) {} + : FieldData(localmesh), x(localmesh), y(localmesh), z(localmesh), + covariant(true), deriv(nullptr), location(CELL_CENTRE) {} Vector2D::Vector2D(const Vector2D &f) - : x(f.x), y(f.y), z(f.z), covariant(f.covariant), deriv(nullptr), location(CELL_CENTRE) {} + : FieldData(f.fielddatamesh), x(f.x), y(f.y), z(f.z), covariant(f.covariant), + deriv(nullptr), location(CELL_CENTRE) {} Vector2D::~Vector2D() { if (deriv != nullptr) { @@ -147,6 +149,8 @@ Vector2D* Vector2D::timeDeriv() { /////////////////// ASSIGNMENT //////////////////// Vector2D & Vector2D::operator=(const Vector2D &rhs) { + fielddatamesh = rhs.fielddatamesh; + x = rhs.x; y = rhs.y; z = rhs.z; diff --git a/src/field/vector3d.cxx b/src/field/vector3d.cxx index c7df254d0e..82d737d177 100644 --- a/src/field/vector3d.cxx +++ b/src/field/vector3d.cxx @@ -37,10 +37,12 @@ #include Vector3D::Vector3D(Mesh *localmesh) - : x(localmesh), y(localmesh), z(localmesh), covariant(true), deriv(nullptr), location(CELL_CENTRE) {} + : FieldData(localmesh), x(localmesh), y(localmesh), z(localmesh), + covariant(true), deriv(nullptr), location(CELL_CENTRE) {} Vector3D::Vector3D(const Vector3D &f) - : x(f.x), y(f.y), z(f.z), covariant(f.covariant), deriv(nullptr), location(CELL_CENTRE) {} + : FieldData(f.fielddatamesh), x(f.x), y(f.y), z(f.z), covariant(f.covariant), + deriv(nullptr), location(CELL_CENTRE) {} Vector3D::~Vector3D() { if (deriv != nullptr) { @@ -148,6 +150,8 @@ Vector3D* Vector3D::timeDeriv() { /////////////////// ASSIGNMENT //////////////////// Vector3D & Vector3D::operator=(const Vector3D &rhs) { + fielddatamesh = rhs.fielddatamesh; + x = rhs.x; y = rhs.y; z = rhs.z; @@ -158,6 +162,8 @@ Vector3D & Vector3D::operator=(const Vector3D &rhs) { } Vector3D & Vector3D::operator=(const Vector2D &rhs) { + fielddatamesh = rhs.x.getMesh(); + x = rhs.x; y = rhs.y; z = rhs.z; diff --git a/src/mesh/boundary_standard.cxx b/src/mesh/boundary_standard.cxx index c82d8f7c0f..e0a69d1444 100644 --- a/src/mesh/boundary_standard.cxx +++ b/src/mesh/boundary_standard.cxx @@ -23,88 +23,92 @@ lead to an out of bounds access error later but we add it here to provide a more explanatory message. */ -void verifyNumPoints(BoundaryRegion *region, int ptsRequired) { - TRACE("Verifying number of points available for BC"); +namespace { + void verifyNumPoints(BoundaryRegion *region, int ptsRequired) { + TRACE("Verifying number of points available for BC"); #ifndef CHECK - return; //No checking so just return + return; //No checking so just return #else - int ptsAvailGlobal, ptsAvailLocal, ptsAvail; - string side, gridType; - - //Initialise var in case of no match and CHECK<=2 - ptsAvail = ptsRequired; //Ensures test passes without exception - - switch(region->location) { - case BNDRY_XIN: - case BNDRY_XOUT: { - side = "x"; - - //Here 2*mesh->xstart is the total number of guard/boundary cells - ptsAvailGlobal = mesh->GlobalNx - 2*mesh->xstart; - - //Work out how many processor local points we have excluding boundaries - //but including ghost/guard cells - ptsAvailLocal = mesh->LocalNx; - if(mesh->firstX()) ptsAvailLocal -= mesh->xstart; - if(mesh->lastX()) ptsAvailLocal -= mesh->xstart; - - //Now decide if it's a local or global limit, prefer global if a tie - if(ptsAvailGlobal <= ptsAvailLocal){ - ptsAvail = ptsAvailGlobal; - gridType = "global"; - }else{ - ptsAvail = ptsAvailLocal; - gridType = "local"; - } + Mesh* localmesh = region->localmesh; - break; - } - case BNDRY_YUP: - case BNDRY_YDOWN: { - side = "y"; + int ptsAvailGlobal, ptsAvailLocal, ptsAvail; + string side, gridType; + + //Initialise var in case of no match and CHECK<=2 + ptsAvail = ptsRequired; //Ensures test passes without exception - //Here 2*mesh->ystart is the total number of guard/boundary cells - ptsAvailGlobal = mesh->GlobalNy - 2*mesh->ystart; + switch(region->location) { + case BNDRY_XIN: + case BNDRY_XOUT: { + side = "x"; - //Work out how many processor local points we have excluding boundaries - //but including ghost/guard cells - ptsAvailLocal = mesh->LocalNy; - if(mesh->firstY()) ptsAvailLocal -= mesh->ystart; - if(mesh->lastY()) ptsAvailLocal -= mesh->ystart; + //Here 2*localmesh->xstart is the total number of guard/boundary cells + ptsAvailGlobal = localmesh->GlobalNx - 2*localmesh->xstart; - //Now decide if it's a local or global limit, prefer global if a tie - if(ptsAvailGlobal <= ptsAvailLocal){ - ptsAvail = ptsAvailGlobal; - gridType = "global"; - }else{ - ptsAvail = ptsAvailLocal; - gridType = "local"; + //Work out how many processor local points we have excluding boundaries + //but including ghost/guard cells + ptsAvailLocal = localmesh->LocalNx; + if(localmesh->firstX()) ptsAvailLocal -= localmesh->xstart; + if(localmesh->lastX()) ptsAvailLocal -= localmesh->xstart; + + //Now decide if it's a local or global limit, prefer global if a tie + if(ptsAvailGlobal <= ptsAvailLocal){ + ptsAvail = ptsAvailGlobal; + gridType = "global"; + }else{ + ptsAvail = ptsAvailLocal; + gridType = "local"; + } + + break; } + case BNDRY_YUP: + case BNDRY_YDOWN: { + side = "y"; - break; - } + //Here 2*localmesh->ystart is the total number of guard/boundary cells + ptsAvailGlobal = localmesh->GlobalNy - 2*localmesh->ystart; + + //Work out how many processor local points we have excluding boundaries + //but including ghost/guard cells + ptsAvailLocal = localmesh->LocalNy; + if(localmesh->firstY()) ptsAvailLocal -= localmesh->ystart; + if(localmesh->lastY()) ptsAvailLocal -= localmesh->ystart; + + //Now decide if it's a local or global limit, prefer global if a tie + if(ptsAvailGlobal <= ptsAvailLocal){ + ptsAvail = ptsAvailGlobal; + gridType = "global"; + }else{ + ptsAvail = ptsAvailLocal; + gridType = "local"; + } + + break; + } #if CHECK > 2 //Only fail on Unrecognised boundary for extreme checking - default : { - throw BoutException("Unrecognised boundary region (%s) for verifyNumPoints.",region->location); - } + default : { + throw BoutException("Unrecognised boundary region (%s) for verifyNumPoints.",region->location); + } #endif - } + } - //Now check we have enough points and if not throw an exception - if(ptsAvail < ptsRequired){ - throw BoutException("Too few %s grid points for %s boundary, have %d but need at least %d", - gridType.c_str(),side.c_str(),ptsAvail,ptsRequired); - } + //Now check we have enough points and if not throw an exception + if(ptsAvail < ptsRequired){ + throw BoutException("Too few %s grid points for %s boundary, have %d but need at least %d", + gridType.c_str(),side.c_str(),ptsAvail,ptsRequired); + } #endif + } } /////////////////////////////////////////////////////////////// BoundaryOp* BoundaryDirichlet::clone(BoundaryRegion *region, const list &args){ - verifyNumPoints(region,1); + verifyNumPoints(region, 1); std::shared_ptr newgen; if(!args.empty()) { @@ -122,6 +126,8 @@ void BoundaryDirichlet::apply(Field2D &f,BoutReal t) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* localmesh = f.getMesh(); + bndry->first(); // Decide which generator to use @@ -135,7 +141,7 @@ void BoundaryDirichlet::apply(Field2D &f,BoutReal t) { // Check for staggered grids CELL_LOC loc = f.getLocation(); - if(mesh->StaggerGrids && loc != CELL_CENTRE) { + if(localmesh->StaggerGrids && loc != CELL_CENTRE && loc != CELL_ZLOW) { // Staggered. Need to apply slightly differently if( loc == CELL_XLOW ) { @@ -146,9 +152,9 @@ void BoundaryDirichlet::apply(Field2D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { if(fg) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) - + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -167,9 +173,9 @@ void BoundaryDirichlet::apply(Field2D &f,BoutReal t) { // Inner x boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { if(fg) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) - + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -190,8 +196,8 @@ void BoundaryDirichlet::apply(Field2D &f,BoutReal t) { if(fg) { // x norm is shifted by half a grid point because it is staggered. // y norm is located half way between first grid cell and guard cell. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - 1) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - 1) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - bndry->by) ); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } f(bndry->x,bndry->y) = 2*val - f(bndry->x-bndry->bx, bndry->y-bndry->by); @@ -200,8 +206,8 @@ void BoundaryDirichlet::apply(Field2D &f,BoutReal t) { for(int i=1;iwidth;i++) { int xi = bndry->x ; int yi = bndry->y + i*bndry->by; - f(xi, yi) = 2*f(xi, yi - bndry->by) - f(xi, yi - 2*bndry->by); - } + f(xi, yi) = 2*f(xi, yi - bndry->by) - f(xi, yi - 2*bndry->by); + } } } } @@ -213,9 +219,9 @@ void BoundaryDirichlet::apply(Field2D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { if(fg) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) - + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + + localmesh->GlobalY(bndry->y - bndry->by) ); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -233,9 +239,9 @@ void BoundaryDirichlet::apply(Field2D &f,BoutReal t) { // Lower y boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { if(fg) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) - + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + + localmesh->GlobalY(bndry->y - bndry->by) ); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -255,10 +261,10 @@ void BoundaryDirichlet::apply(Field2D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { if(fg) { - // x norm is located half way between first grid cell and guard cell. - // y norm is shifted by half a grid point because it is staggered. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - 1) ); + // x norm is located half way between first grid cell and guard cell. + // y norm is shifted by half a grid point because it is staggered. + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - 1) ); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } f(bndry->x,bndry->y) = 2*val - f(bndry->x-bndry->bx, bndry->y-bndry->by); @@ -271,6 +277,8 @@ void BoundaryDirichlet::apply(Field2D &f,BoutReal t) { } } } + } else { + throw BoutException("Unhandled staggering in BoundaryDirichlet."); } } else { // Non-staggered, standard case @@ -279,11 +287,11 @@ void BoundaryDirichlet::apply(Field2D &f,BoutReal t) { if(fg) { // Calculate the X and Y normalised values half-way between the guard cell and grid cell - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) // In the guard cell - + mesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) // In the guard cell + + localmesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) // In the guard cell - + mesh->GlobalY(bndry->y - bndry->by) ); // the grid cell + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) // In the guard cell + + localmesh->GlobalY(bndry->y - bndry->by) ); // the grid cell val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -309,6 +317,8 @@ void BoundaryDirichlet::apply(Field3D &f,BoutReal t) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* localmesh = f.getMesh(); + bndry->first(); // Decide which generator to use @@ -321,7 +331,7 @@ void BoundaryDirichlet::apply(Field3D &f,BoutReal t) { // Check for staggered grids CELL_LOC loc = f.getLocation(); - if(mesh->StaggerGrids && loc != CELL_CENTRE) { + if(localmesh->StaggerGrids && loc != CELL_CENTRE && loc != CELL_ZLOW) { // Staggered. Need to apply slightly differently if( loc == CELL_XLOW ) { @@ -331,13 +341,13 @@ void BoundaryDirichlet::apply(Field3D &f,BoutReal t) { // Outer x boundary for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) - + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); } f(bndry->x,bndry->y, zk) = val; @@ -355,13 +365,13 @@ void BoundaryDirichlet::apply(Field3D &f,BoutReal t) { // Inner x boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) - + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); } f(bndry->x - bndry->bx,bndry->y, zk) = val; f(bndry->x,bndry->y, zk) = f(bndry->x - bndry->bx,bndry->y, zk); @@ -381,12 +391,12 @@ void BoundaryDirichlet::apply(Field3D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { // x norm is shifted by half a grid point because it is staggered. // y norm is located half way between first grid cell and guard cell. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - 1) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - 1) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - bndry->by) ); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); } f(bndry->x,bndry->y,zk) = 2*val - f(bndry->x-bndry->bx, bndry->y-bndry->by, zk); @@ -408,11 +418,11 @@ void BoundaryDirichlet::apply(Field3D &f,BoutReal t) { // Upper y boundary boundary for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - bndry->by) ); - for(int zk=0;zkLocalNz;zk++) { + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - bndry->by) ); + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); } f(bndry->x,bndry->y,zk) = val; @@ -430,12 +440,12 @@ void BoundaryDirichlet::apply(Field3D &f,BoutReal t) { // Lower y boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - bndry->by) ); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); } f(bndry->x,bndry->y - bndry->by, zk) = val; @@ -452,14 +462,14 @@ void BoundaryDirichlet::apply(Field3D &f,BoutReal t) { if(bndry->bx != 0){ // x boundaries for(; !bndry->isDone(); bndry->next1d()) { - // x norm is located half way between first grid cell and guard cell. - // y norm is shifted by half a grid point because it is staggered. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - 1) ); + // x norm is located half way between first grid cell and guard cell. + // y norm is shifted by half a grid point because it is staggered. + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - 1) ); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg) - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); f(bndry->x,bndry->y,zk) = 2*val - f(bndry->x-bndry->bx, bndry->y-bndry->by, zk); @@ -473,21 +483,23 @@ void BoundaryDirichlet::apply(Field3D &f,BoutReal t) { } } } + } else { + throw BoutException("Unhandled staggering in BoundaryDirichlet."); } } else { // Standard (non-staggered) case for(; !bndry->isDone(); bndry->next1d()) { // Calculate the X and Y normalised values half-way between the guard cell and grid cell - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) // In the guard cell - + mesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) // In the guard cell + + localmesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) // In the guard cell - + mesh->GlobalY(bndry->y - bndry->by) ); // the grid cell + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) // In the guard cell + + localmesh->GlobalY(bndry->y - bndry->by) ); // the grid cell - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); } f(bndry->x,bndry->y,zk) = 2*val - f(bndry->x-bndry->bx, bndry->y-bndry->by, zk); @@ -525,11 +537,11 @@ void BoundaryDirichlet::apply(Field3D &f,BoutReal t) { // Set any other guard cells using the values on the cells int xi = bndry->x + i*bndry->bx; int yi = bndry->y + i*bndry->by; - xnorm = mesh->GlobalX(xi); - ynorm = mesh->GlobalY(yi); - for(int zk=0;zkLocalNz;zk++) { + xnorm = localmesh->GlobalX(xi); + ynorm = localmesh->GlobalY(yi); + for(int zk=0;zkLocalNz;zk++) { if(fg) { - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); } f(xi, yi, zk) = val; } @@ -546,9 +558,10 @@ void BoundaryDirichlet::apply_ddt(Field2D &f) { } void BoundaryDirichlet::apply_ddt(Field3D &f) { + Mesh* localmesh = f.getMesh(); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) + for(int z=0;zLocalNz;z++) (*dt)(bndry->x,bndry->y,z) = 0.; // Set time derivative to zero } @@ -557,7 +570,7 @@ void BoundaryDirichlet::apply_ddt(Field3D &f) { // New implementation, accurate to higher order BoundaryOp* BoundaryDirichlet_O3::clone(BoundaryRegion *region, const list &args){ - verifyNumPoints(region,2); + verifyNumPoints(region, 2); std::shared_ptr newgen = nullptr; if(!args.empty()) { // First argument should be an expression @@ -574,6 +587,8 @@ void BoundaryDirichlet_O3::apply(Field2D &f,BoutReal t) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* localmesh = f.getMesh(); + bndry->first(); // Decide which generator to use @@ -587,7 +602,7 @@ void BoundaryDirichlet_O3::apply(Field2D &f,BoutReal t) { // Check for staggered grids CELL_LOC loc = f.getLocation(); - if(mesh->StaggerGrids && loc != CELL_CENTRE) { + if(localmesh->StaggerGrids && loc != CELL_CENTRE && loc != CELL_ZLOW) { // Staggered. Need to apply slightly differently if( loc == CELL_XLOW) { @@ -597,8 +612,8 @@ void BoundaryDirichlet_O3::apply(Field2D &f,BoutReal t) { // Outer x boundary for(; !bndry->isDone(); bndry->next1d()) { if(fg) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -616,8 +631,8 @@ void BoundaryDirichlet_O3::apply(Field2D &f,BoutReal t) { // Inner x boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { if(fg) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } f(bndry->x - bndry->bx,bndry->y) = val; @@ -634,9 +649,9 @@ void BoundaryDirichlet_O3::apply(Field2D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { if(fg) { // x norm is shifted by half a grid point because it is staggered. - // y norm is located half way between first grid cell and guard cell. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - 1) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - bndry->by) ); + // y norm is located half way between first grid cell and guard cell. + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - 1) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - bndry->by) ); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -661,8 +676,8 @@ void BoundaryDirichlet_O3::apply(Field2D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { if(fg) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - bndry->by) ); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -681,8 +696,8 @@ void BoundaryDirichlet_O3::apply(Field2D &f,BoutReal t) { // Lower y boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { if(fg) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - bndry->by) ); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -703,8 +718,8 @@ void BoundaryDirichlet_O3::apply(Field2D &f,BoutReal t) { if(fg) { // x norm is located half way between first grid cell and guard cell. // y norm is shifted by half a grid point because it is staggered. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - 1) ); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - 1) ); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -719,6 +734,8 @@ void BoundaryDirichlet_O3::apply(Field2D &f,BoutReal t) { } } } + } else { + throw BoutException("Unhandled staggering in BoundaryDirichlet_O3."); } } else { @@ -728,11 +745,11 @@ void BoundaryDirichlet_O3::apply(Field2D &f,BoutReal t) { if(fg) { // Calculate the X and Y normalised values half-way between the guard cell and grid cell - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) // In the guard cell - + mesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) // In the guard cell + + localmesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) // In the guard cell - + mesh->GlobalY(bndry->y - bndry->by) ); // the grid cell + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) // In the guard cell + + localmesh->GlobalY(bndry->y - bndry->by) ); // the grid cell val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -759,6 +776,8 @@ void BoundaryDirichlet_O3::apply(Field3D &f,BoutReal t) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* localmesh = f.getMesh(); + bndry->first(); // Decide which generator to use @@ -771,7 +790,7 @@ void BoundaryDirichlet_O3::apply(Field3D &f,BoutReal t) { // Check for staggered grids CELL_LOC loc = f.getLocation(); - if(mesh->StaggerGrids && loc != CELL_CENTRE) { + if(localmesh->StaggerGrids && loc != CELL_CENTRE && loc != CELL_ZLOW) { // Staggered. Need to apply slightly differently if( loc == CELL_XLOW ) { @@ -781,12 +800,12 @@ void BoundaryDirichlet_O3::apply(Field3D &f,BoutReal t) { // Outer x boundary for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); } f(bndry->x,bndry->y, zk) = val; @@ -804,12 +823,12 @@ void BoundaryDirichlet_O3::apply(Field3D &f,BoutReal t) { // Inner x boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); } f(bndry->x - bndry->bx,bndry->y, zk) = val; @@ -828,13 +847,13 @@ void BoundaryDirichlet_O3::apply(Field3D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { // x norm is shifted by half a grid point because it is staggered. - // y norm is located half way between first grid cell and guard cell. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - 1) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - bndry->by) ); + // y norm is located half way between first grid cell and guard cell. + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - 1) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - bndry->by) ); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg) - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); f(bndry->x,bndry->y,zk) = (8./3)*val - 2.*f(bndry->x-bndry->bx, bndry->y-bndry->by,zk) + f(bndry->x-2*bndry->bx, bndry->y-2*bndry->by,zk)/3.; @@ -856,11 +875,11 @@ void BoundaryDirichlet_O3::apply(Field3D &f,BoutReal t) { // Upper y boundary for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - bndry->by) ); - for(int zk=0;zkLocalNz;zk++) { + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - bndry->by) ); + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); } f(bndry->x,bndry->y,zk) = val; @@ -878,12 +897,12 @@ void BoundaryDirichlet_O3::apply(Field3D &f,BoutReal t) { // Lower y boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*(mesh->GlobalY(bndry->y)+ mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*(localmesh->GlobalY(bndry->y)+ localmesh->GlobalY(bndry->y - bndry->by) ); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); } f(bndry->x,bndry->y - bndry->by, zk) = val; @@ -902,12 +921,12 @@ void BoundaryDirichlet_O3::apply(Field3D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { // x norm is located half way between first grid cell and guard cell. // y norm is shifted by half a grid point because it is staggered. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - 1) ); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - 1) ); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg) - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); f(bndry->x,bndry->y,zk) = (8./3)*val - 2.*f(bndry->x-bndry->bx, bndry->y-bndry->by,zk) + f(bndry->x-2*bndry->bx, bndry->y-2*bndry->by,zk)/3.; @@ -921,21 +940,23 @@ void BoundaryDirichlet_O3::apply(Field3D &f,BoutReal t) { } } } + } else { + throw BoutException("Unhandled staggering in BoundaryDirichlet_O3."); } } else { // Standard (non-staggered) case for(; !bndry->isDone(); bndry->next1d()) { // Calculate the X and Y normalised values half-way between the guard cell and grid cell - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) // In the guard cell - + mesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) // In the guard cell + + localmesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) // In the guard cell - + mesh->GlobalY(bndry->y - bndry->by) ); // the grid cell + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) // In the guard cell + + localmesh->GlobalY(bndry->y - bndry->by) ); // the grid cell - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg) - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); f(bndry->x,bndry->y,zk) = (8./3)*val - 2.*f(bndry->x-bndry->bx, bndry->y-bndry->by,zk) + f(bndry->x-2*bndry->bx, bndry->y-2*bndry->by,zk)/3.; @@ -958,11 +979,11 @@ void BoundaryDirichlet_O3::apply_ddt(Field2D &f) { } void BoundaryDirichlet_O3::apply_ddt(Field3D &f) { + Mesh* localmesh = f.getMesh(); Field3D *dt = f.timeDeriv(); - bndry->first() ; for(bndry->first(); !bndry->isDone(); bndry->next()){ - for(int z=0;zLocalNz;z++){ + for(int z=0;zLocalNz;z++){ (*dt)(bndry->x,bndry->y,z) = 0.; // Set time derivative to zero } } @@ -972,7 +993,7 @@ void BoundaryDirichlet_O3::apply_ddt(Field3D &f) { // Extrapolate to calculate boundary cell to 4th-order BoundaryOp* BoundaryDirichlet_O4::clone(BoundaryRegion *region, const list &args){ - verifyNumPoints(region,3); + verifyNumPoints(region, 3); std::shared_ptr newgen = nullptr; if(!args.empty()) { // First argument should be an expression @@ -989,6 +1010,8 @@ void BoundaryDirichlet_O4::apply(Field2D &f,BoutReal t) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* localmesh = f.getMesh(); + bndry->first(); // Decide which generator to use @@ -1002,7 +1025,7 @@ void BoundaryDirichlet_O4::apply(Field2D &f,BoutReal t) { // Check for staggered grids CELL_LOC loc = f.getLocation(); - if(mesh->StaggerGrids && loc != CELL_CENTRE) { + if(localmesh->StaggerGrids && loc != CELL_CENTRE && loc != CELL_ZLOW) { // Staggered. Need to apply slightly differently if(loc == CELL_XLOW ) { @@ -1013,8 +1036,8 @@ void BoundaryDirichlet_O4::apply(Field2D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { if(fg) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } f(bndry->x,bndry->y) = val; @@ -1033,9 +1056,9 @@ void BoundaryDirichlet_O4::apply(Field2D &f,BoutReal t) { // Inner boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { if(fg) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) - + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -1057,8 +1080,8 @@ void BoundaryDirichlet_O4::apply(Field2D &f,BoutReal t) { if(fg) { // x norm is shifted by half a grid point because it is staggered. // y norm is located half way between first grid cell and guard cell. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - 1) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - 1) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - bndry->by) ); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -1082,9 +1105,9 @@ void BoundaryDirichlet_O4::apply(Field2D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { if(fg) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) - + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + + localmesh->GlobalY(bndry->y - bndry->by) ); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } f(bndry->x,bndry->y) = val; @@ -1102,9 +1125,9 @@ void BoundaryDirichlet_O4::apply(Field2D &f,BoutReal t) { // Inner y boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { if(fg) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) - + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + + localmesh->GlobalY(bndry->y - bndry->by) ); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -1128,8 +1151,8 @@ void BoundaryDirichlet_O4::apply(Field2D &f,BoutReal t) { if(fg) { // x norm is located half way between first grid cell and guard cell. // y norm is shifted by half a grid point because it is staggered. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - 1) ); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - 1) ); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -1145,6 +1168,8 @@ void BoundaryDirichlet_O4::apply(Field2D &f,BoutReal t) { } } } + } else { + throw BoutException("Unhandled staggering in BoundaryDirichlet_O4."); } } else { @@ -1154,11 +1179,11 @@ void BoundaryDirichlet_O4::apply(Field2D &f,BoutReal t) { if(fg) { // Calculate the X and Y normalised values half-way between the guard cell and grid cell - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) // In the guard cell - + mesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) // In the guard cell + + localmesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) // In the guard cell - + mesh->GlobalY(bndry->y - bndry->by) ); // the grid cell + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) // In the guard cell + + localmesh->GlobalY(bndry->y - bndry->by) ); // the grid cell val = fg->generate(xnorm,TWOPI*ynorm,0.0, t); } @@ -1186,6 +1211,8 @@ void BoundaryDirichlet_O4::apply(Field3D &f,BoutReal t) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* localmesh = f.getMesh(); + bndry->first(); // Decide which generator to use @@ -1198,7 +1225,7 @@ void BoundaryDirichlet_O4::apply(Field3D &f,BoutReal t) { // Check for staggered grids CELL_LOC loc = f.getLocation(); - if(mesh->StaggerGrids && loc != CELL_CENTRE) { + if(localmesh->StaggerGrids && loc != CELL_CENTRE && loc != CELL_ZLOW) { // Staggered. Need to apply slightly differently if( loc == CELL_XLOW ) { @@ -1208,13 +1235,13 @@ void BoundaryDirichlet_O4::apply(Field3D &f,BoutReal t) { // Outer x boundary for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) - + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); } f(bndry->x,bndry->y, zk) = val; @@ -1232,13 +1259,13 @@ void BoundaryDirichlet_O4::apply(Field3D &f,BoutReal t) { // Inner x boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) - + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg) - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); f(bndry->x - bndry->bx,bndry->y, zk) = val; @@ -1258,12 +1285,12 @@ void BoundaryDirichlet_O4::apply(Field3D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { // x norm is shifted by half a grid point because it is staggered. // y norm is located half way between first grid cell and guard cell. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - 1) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - 1) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - bndry->by) ); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg) { - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); } f(bndry->x,bndry->y,zk) = (16./5)*val - 3.*f(bndry->x-bndry->bx, bndry->y-bndry->by,zk) + f(bndry->x-2*bndry->bx, bndry->y-2*bndry->by,zk) - (1./5)*f(bndry->x-3*bndry->bx, bndry->y-3*bndry->by,zk); @@ -1285,12 +1312,12 @@ void BoundaryDirichlet_O4::apply(Field3D &f,BoutReal t) { // Outer y boundary for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) - + mesh->GlobalY(bndry->y - bndry->by) ); - for(int zk=0;zkLocalNz;zk++) { + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + + localmesh->GlobalY(bndry->y - bndry->by) ); + for(int zk=0;zkLocalNz;zk++) { if(fg) - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); f(bndry->x,bndry->y,zk) = val; @@ -1308,13 +1335,13 @@ void BoundaryDirichlet_O4::apply(Field3D &f,BoutReal t) { // Inner y boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) - + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + + localmesh->GlobalY(bndry->y - bndry->by) ); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg) - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); f(bndry->x,bndry->y - bndry->by, zk) = val; @@ -1334,12 +1361,12 @@ void BoundaryDirichlet_O4::apply(Field3D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { // x norm is located half way between first grid cell and guard cell. // y norm is shifted by half a grid point because it is staggered. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - 1) ); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - 1) ); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg) - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); f(bndry->x,bndry->y,zk) = (16./5)*val - 3.*f(bndry->x-bndry->bx, bndry->y-bndry->by,zk) + f(bndry->x-2*bndry->bx, bndry->y-2*bndry->by,zk) - (1./5)*f(bndry->x-3*bndry->bx, bndry->y-3*bndry->by,zk); @@ -1353,21 +1380,23 @@ void BoundaryDirichlet_O4::apply(Field3D &f,BoutReal t) { } } } + } else { + throw BoutException("Unhandled staggering in BoundaryDirichlet_O4."); } } else { // Standard (non-staggered) case for(; !bndry->isDone(); bndry->next1d()) { // Calculate the X and Y normalised values half-way between the guard cell and grid cell - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) // In the guard cell - + mesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) // In the guard cell + + localmesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) // In the guard cell - + mesh->GlobalY(bndry->y - bndry->by) ); // the grid cell + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) // In the guard cell + + localmesh->GlobalY(bndry->y - bndry->by) ); // the grid cell - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg) - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz), t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz), t); f(bndry->x,bndry->y,zk) = (16./5)*val - 3.*f(bndry->x-bndry->bx, bndry->y-bndry->by,zk) + f(bndry->x-2*bndry->bx, bndry->y-2*bndry->by,zk) - (1./5)*f(bndry->x-3*bndry->bx, bndry->y-3*bndry->by,zk); @@ -1390,9 +1419,10 @@ void BoundaryDirichlet_O4::apply_ddt(Field2D &f) { } void BoundaryDirichlet_O4::apply_ddt(Field3D &f) { + Mesh* localmesh = f.getMesh(); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) + for(int z=0;zLocalNz;z++) (*dt)(bndry->x,bndry->y,z) = 0.; // Set time derivative to zero } @@ -1401,7 +1431,7 @@ void BoundaryDirichlet_O4::apply_ddt(Field3D &f) { BoundaryOp* BoundaryDirichlet_2ndOrder::clone(BoundaryRegion *region, const list &args) { output << "WARNING: Use of boundary condition \"dirichlet_2ndorder\" is deprecated!\n"; output << " Consider using \"dirichlet\" instead\n"; - verifyNumPoints(region,2); + verifyNumPoints(region, 2); if(!args.empty()) { // First argument should be a value val = stringToReal(args.front()); @@ -1424,10 +1454,11 @@ void BoundaryDirichlet_2ndOrder::apply(Field2D &f) { } void BoundaryDirichlet_2ndOrder::apply(Field3D &f) { + Mesh* localmesh = f.getMesh(); // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used for(bndry->first(); !bndry->isDone(); bndry->next1d()) - for(int z=0;zLocalNz;z++) { + for(int z=0;zLocalNz;z++) { f(bndry->x,bndry->y,z) = 8./3.*val - 2.*f(bndry->x-bndry->bx,bndry->y-bndry->by,z) + 1./3.*f(bndry->x-2*bndry->bx,bndry->y-2*bndry->by,z); #ifdef BOUNDARY_CONDITIONS_UPGRADE_EXTRAPOLATE_FOR_2ND_ORDER f(bndry->x+bndry->bx,bndry->y+bndry->by,z) = 3.*f(bndry->x,bndry->y,z) - 3.*f(bndry->x-bndry->bx,bndry->y-bndry->by,z) + f(bndry->x-2*bndry->bx,bndry->y-2*bndry->by,z); @@ -1444,16 +1475,17 @@ void BoundaryDirichlet_2ndOrder::apply_ddt(Field2D &f) { } void BoundaryDirichlet_2ndOrder::apply_ddt(Field3D &f) { + Mesh* localmesh = f.getMesh(); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) + for(int z=0;zLocalNz;z++) (*dt)(bndry->x,bndry->y,z) = 0.; // Set time derivative to zero } /////////////////////////////////////////////////////////////// BoundaryOp* BoundaryDirichlet_4thOrder::clone(BoundaryRegion *region, const list &args) { - verifyNumPoints(region,4); + verifyNumPoints(region, 4); if(!args.empty()) { // First argument should be a value val = stringToReal(args.front()); @@ -1471,9 +1503,10 @@ void BoundaryDirichlet_4thOrder::apply(Field2D &f) { } void BoundaryDirichlet_4thOrder::apply(Field3D &f) { + Mesh* localmesh = f.getMesh(); // Set (at 4th order) the value at the mid-point between the guard cell and the grid cell to be val for(bndry->first(); !bndry->isDone(); bndry->next1d()) - for(int z=0;zLocalNz;z++) { + for(int z=0;zLocalNz;z++) { f(bndry->x,bndry->y,z) = 128./35.*val - 4.*f(bndry->x-bndry->bx,bndry->y-bndry->by,z) + 2.*f(bndry->x-2*bndry->bx,bndry->y-2*bndry->by,z) - 4./3.*f(bndry->x-3*bndry->bx,bndry->y-3*bndry->by,z) + 1./7.*f(bndry->x-4*bndry->bx,bndry->y-4*bndry->by,z); f(bndry->x+bndry->bx,bndry->y+bndry->by,z) = -128./5.*val + 9.*f(bndry->x,bndry->y,z) + 18.*f(bndry->x-bndry->bx,bndry->y-bndry->by,z) -4.*f(bndry->x-2*bndry->bx,bndry->y-2*bndry->by,z) + 3./5.*f(bndry->x-3*bndry->bx,bndry->y-3*bndry->by,z); } @@ -1486,16 +1519,17 @@ void BoundaryDirichlet_4thOrder::apply_ddt(Field2D &f) { } void BoundaryDirichlet_4thOrder::apply_ddt(Field3D &f) { + Mesh* localmesh = f.getMesh(); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) + for(int z=0;zLocalNz;z++) (*dt)(bndry->x,bndry->y,z) = 0.; // Set time derivative to zero } /////////////////////////////////////////////////////////////// BoundaryOp* BoundaryNeumann_NonOrthogonal::clone(BoundaryRegion *region, const list &args) { - verifyNumPoints(region,1); + verifyNumPoints(region, 1); if(!args.empty()) { output << "WARNING: arguments is set to BoundaryNeumann None Zero Gradient\n"; // First argument should be a value @@ -1506,9 +1540,10 @@ BoundaryOp* BoundaryNeumann_NonOrthogonal::clone(BoundaryRegion *region, const l } void BoundaryNeumann_NonOrthogonal::apply(Field2D &f) { + Mesh* localmesh = f.getMesh(); Coordinates *metric = f.getCoordinates(); // Calculate derivatives for metric use - mesh->communicate(f); + localmesh->communicate(f); Field2D dfdy = DDY(f); // Loop over all elements and set equal to the next point in for(bndry->first(); !bndry->isDone(); bndry->next1d()) { @@ -1546,9 +1581,10 @@ void BoundaryNeumann_NonOrthogonal::apply(Field2D &f) { } void BoundaryNeumann_NonOrthogonal::apply(Field3D &f) { + Mesh* localmesh = f.getMesh(); Coordinates *metric = f.getCoordinates(); // Calculate derivatives for metric use - mesh->communicate(f); + localmesh->communicate(f); Field3D dfdy = DDY(f); Field3D dfdz = DDZ(f); // Loop over all elements and set equal to the next point in @@ -1560,7 +1596,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D &f) { // Have to use derivatives at last gridpoint instead of derivatives on boundary layer // because derivative values don't exist in boundary region // NOTE: should be fixed to interpolate to boundary line - for(int z=0;zLocalNz;z++) { + for(int z=0;zLocalNz;z++) { BoutReal xshift = g12shift*dfdy(bndry->x-bndry->bx,bndry->y,z) + g13shift*dfdz(bndry->x-bndry->bx,bndry->y,z); if(bndry->bx != 0 && bndry->by == 0) { @@ -1594,7 +1630,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D &f) { BoundaryOp* BoundaryNeumann2::clone(BoundaryRegion *region, const list &args) { output << "WARNING: Use of boundary condition \"neumann2\" is deprecated!\n"; output << " Consider using \"neumann\" instead\n"; - verifyNumPoints(region,2); + verifyNumPoints(region, 2); if(!args.empty()) { output << "WARNING: Ignoring arguments to BoundaryNeumann2\n"; } @@ -1608,8 +1644,9 @@ void BoundaryNeumann2::apply(Field2D &f) { } void BoundaryNeumann2::apply(Field3D &f) { + Mesh* localmesh = f.getMesh(); for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) + for(int z=0;zLocalNz;z++) f(bndry->x, bndry->y, z) = (4.*f(bndry->x - bndry->bx, bndry->y - bndry->by, z) - f(bndry->x - 2*bndry->bx, bndry->y - 2*bndry->by, z))/3.; } @@ -1619,9 +1656,9 @@ BoundaryOp* BoundaryNeumann_2ndOrder::clone(BoundaryRegion *region, const listfirst(); !bndry->isDone(); bndry->next1d()) - for(int z=0;zLocalNz;z++) { + for(int z=0;zLocalNz;z++) { BoutReal delta = bndry->bx*metric->dx(bndry->x,bndry->y)+bndry->by*metric->dy(bndry->x,bndry->y); f(bndry->x,bndry->y,z) = f(bndry->x-bndry->bx,bndry->y-bndry->by,z) + val*delta; #ifdef BOUNDARY_CONDITIONS_UPGRADE_EXTRAPOLATE_FOR_2ND_ORDER @@ -1671,16 +1709,17 @@ void BoundaryNeumann_2ndOrder::apply_ddt(Field2D &f) { } void BoundaryNeumann_2ndOrder::apply_ddt(Field3D &f) { + Mesh* localmesh = f.getMesh(); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) + for(int z=0;zLocalNz;z++) (*dt)(bndry->x,bndry->y,z) = 0.; // Set time derivative to zero } /////////////////////////////////////////////////////////////// BoundaryOp* BoundaryNeumann::clone(BoundaryRegion *region, const list &args){ - verifyNumPoints(region,1); + verifyNumPoints(region, 1); std::shared_ptr newgen = nullptr; if(!args.empty()) { // First argument should be an expression @@ -1698,6 +1737,8 @@ void BoundaryNeumann::apply(Field2D &f,BoutReal t) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* localmesh = f.getMesh(); + Coordinates *metric = f.getCoordinates(); bndry->first(); @@ -1712,7 +1753,7 @@ void BoundaryNeumann::apply(Field2D &f,BoutReal t) { // Check for staggered grids CELL_LOC loc = f.getLocation(); - if(mesh->StaggerGrids && loc != CELL_CENTRE) { + if(localmesh->StaggerGrids && loc != CELL_CENTRE && loc != CELL_ZLOW) { // Staggered. Need to apply slightly differently // Use one-sided differencing. Cell is now on // the boundary, so use one-sided differencing @@ -1726,9 +1767,9 @@ void BoundaryNeumann::apply(Field2D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { if(fg) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) - + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t) * metric->dx(bndry->x, bndry->y); } @@ -1750,9 +1791,9 @@ void BoundaryNeumann::apply(Field2D &f,BoutReal t) { if(fg) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) - + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t) * metric->dx(bndry->x, bndry->y); } @@ -1778,8 +1819,8 @@ void BoundaryNeumann::apply(Field2D &f,BoutReal t) { if(fg) { // x norm is shifted by half a grid point because it is staggered. // y norm is located half way between first grid cell and guard cell. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - 1) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - 1) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - bndry->by) ); val = fg->generate(xnorm, TWOPI*ynorm, 0.0, t); } @@ -1799,9 +1840,9 @@ void BoundaryNeumann::apply(Field2D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { if(fg) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) - + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + + localmesh->GlobalY(bndry->y - bndry->by) ); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t) * metric->dx(bndry->x, bndry->y); } @@ -1823,9 +1864,9 @@ void BoundaryNeumann::apply(Field2D &f,BoutReal t) { if(fg) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) - + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + + localmesh->GlobalY(bndry->y - bndry->by) ); val = fg->generate(xnorm,TWOPI*ynorm,0.0, t) * metric->dx(bndry->x, bndry->y - bndry->by); } @@ -1849,8 +1890,8 @@ void BoundaryNeumann::apply(Field2D &f,BoutReal t) { if(fg) { // x norm is located half way between first grid cell and guard cell. // y norm is shifted by half a grid point because it is staggered. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - 1) ); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - 1) ); val = fg->generate(xnorm, TWOPI*ynorm, 0.0, t); } @@ -1861,6 +1902,8 @@ void BoundaryNeumann::apply(Field2D &f,BoutReal t) { } } } + } else { + throw BoutException("Unhandled staggering in BoundaryNeumann."); } } else { @@ -1871,11 +1914,11 @@ void BoundaryNeumann::apply(Field2D &f,BoutReal t) { if(fg) { // Calculate the X and Y normalised values half-way between the guard cell and grid cell - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) // In the guard cell - + mesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) // In the guard cell + + localmesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) // In the guard cell - + mesh->GlobalY(bndry->y - bndry->by) ); // the grid cell + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) // In the guard cell + + localmesh->GlobalY(bndry->y - bndry->by) ); // the grid cell val = fg->generate(xnorm, TWOPI*ynorm, 0.0, t); } @@ -1895,6 +1938,8 @@ void BoundaryNeumann::apply(Field3D &f) { void BoundaryNeumann::apply(Field3D &f,BoutReal t) { + Mesh* localmesh = f.getMesh(); + Coordinates *metric = f.getCoordinates(); bndry->first(); @@ -1909,7 +1954,7 @@ void BoundaryNeumann::apply(Field3D &f,BoutReal t) { // Check for staggered grids CELL_LOC loc = f.getLocation(); - if(mesh->StaggerGrids && loc != CELL_CENTRE) { + if(localmesh->StaggerGrids && loc != CELL_CENTRE && loc != CELL_ZLOW) { // Staggered. Need to apply slightly differently // Use one-sided differencing. Cell is now on // the boundary, so use one-sided differencing @@ -1920,13 +1965,13 @@ void BoundaryNeumann::apply(Field3D &f,BoutReal t) { if(bndry->bx > 0) { // Outer x boundary for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) - + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg) - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz),t) * metric->dx(bndry->x, bndry->y); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz),t) * metric->dx(bndry->x, bndry->y); f(bndry->x,bndry->y, zk) = (4.*f(bndry->x - bndry->bx, bndry->y,zk) - f(bndry->x - 2*bndry->bx, bndry->y,zk) + 2.*val)/3.; @@ -1946,14 +1991,14 @@ void BoundaryNeumann::apply(Field3D &f,BoutReal t) { // Inner x boundary for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) - + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = mesh->GlobalY(bndry->y); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = localmesh->GlobalY(bndry->y); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg) - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz),t) * metric->dx(bndry->x - bndry->bx, bndry->y); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz),t) * metric->dx(bndry->x - bndry->bx, bndry->y); f(bndry->x - bndry->bx,bndry->y, zk) = (4.*f(bndry->x - 2*bndry->bx, bndry->y,zk) - f(bndry->x - 3*bndry->bx, bndry->y,zk) - 2.*val)/3.; @@ -1973,14 +2018,14 @@ void BoundaryNeumann::apply(Field3D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { // x norm is shifted by half a grid point because it is staggered. // y norm is located half way between first grid cell and guard cell. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - 1) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - bndry->by) ); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - 1) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - bndry->by) ); BoutReal delta = bndry->bx*metric->dx(bndry->x,bndry->y)+bndry->by*metric->dy(bndry->x,bndry->y); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz),t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz),t); } f(bndry->x,bndry->y, zk) = f(bndry->x-bndry->bx, bndry->y-bndry->by, zk) + delta*val; if (bndry->width == 2){ @@ -1997,13 +2042,13 @@ void BoundaryNeumann::apply(Field3D &f,BoutReal t) { // Outer y boundary for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) - + mesh->GlobalY(bndry->y - bndry->by) ); - for(int zk=0;zkLocalNz;zk++) { + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + + localmesh->GlobalY(bndry->y - bndry->by) ); + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz),t) * metric->dy(bndry->x, bndry->y); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz),t) * metric->dy(bndry->x, bndry->y); } f(bndry->x,bndry->y,zk) = (4.*f(bndry->x, bndry->y - bndry->by,zk) - f(bndry->x, bndry->y - 2*bndry->by,zk) + 2.*val)/3.; @@ -2023,12 +2068,12 @@ void BoundaryNeumann::apply(Field3D &f,BoutReal t) { // Inner y boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { - BoutReal xnorm = mesh->GlobalX(bndry->x); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) - + mesh->GlobalY(bndry->y - bndry->by) ); - for(int zk=0;zkLocalNz;zk++) { + BoutReal xnorm = localmesh->GlobalX(bndry->x); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + + localmesh->GlobalY(bndry->y - bndry->by) ); + for(int zk=0;zkLocalNz;zk++) { if(fg) - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz),t) * metric->dy(bndry->x, bndry->y - bndry->by); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz),t) * metric->dy(bndry->x, bndry->y - bndry->by); f(bndry->x,bndry->y - bndry->by,zk) = (4.*f(bndry->x, bndry->y - 2*bndry->by,zk) - f(bndry->x, bndry->y - 3*bndry->by,zk) - 2.*val)/3.; @@ -2049,14 +2094,14 @@ void BoundaryNeumann::apply(Field3D &f,BoutReal t) { for(; !bndry->isDone(); bndry->next1d()) { // x norm is located half way between first grid cell and guard cell. // y norm is shifted by half a grid point because it is staggered. - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) + mesh->GlobalX(bndry->x - bndry->bx) ); - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) + mesh->GlobalY(bndry->y - 1) ); + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) + localmesh->GlobalX(bndry->x - bndry->bx) ); + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) + localmesh->GlobalY(bndry->y - 1) ); BoutReal delta = bndry->bx*metric->dx(bndry->x,bndry->y)+bndry->by*metric->dy(bndry->x,bndry->y); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz),t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz),t); } f(bndry->x,bndry->y, zk) = f(bndry->x-bndry->bx, bndry->y-bndry->by, zk) + delta*val; if (bndry->width == 2){ @@ -2065,22 +2110,24 @@ void BoundaryNeumann::apply(Field3D &f,BoutReal t) { } } } + } else { + throw BoutException("Unhandled staggering in BoundaryNeumann."); } } else { for(; !bndry->isDone(); bndry->next1d()) { // Calculate the X and Y normalised values half-way between the guard cell and grid cell - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) // In the guard cell - + mesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) // In the guard cell + + localmesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) // In the guard cell - + mesh->GlobalY(bndry->y - bndry->by) ); // the grid cell + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) // In the guard cell + + localmesh->GlobalY(bndry->y - bndry->by) ); // the grid cell BoutReal delta = bndry->bx*metric->dx(bndry->x,bndry->y)+bndry->by*metric->dy(bndry->x,bndry->y); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz),t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz),t); } f(bndry->x,bndry->y, zk) = f(bndry->x-bndry->bx, bndry->y-bndry->by, zk) + delta*val; if (bndry->width == 2){ @@ -2098,9 +2145,10 @@ void BoundaryNeumann::apply_ddt(Field2D &f) { } void BoundaryNeumann::apply_ddt(Field3D &f) { + Mesh* localmesh = f.getMesh(); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) + for(int z=0;zLocalNz;z++) (*dt)(bndry->x,bndry->y,z) = 0.; // Set time derivative to zero } @@ -2123,6 +2171,8 @@ void BoundaryNeumann_O4::apply(Field2D &f,BoutReal t) { // Set (at 4th order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* localmesh = f.getMesh(); + bndry->first(); // Decide which generator to use @@ -2134,7 +2184,7 @@ void BoundaryNeumann_O4::apply(Field2D &f,BoutReal t) { // Check for staggered grids CELL_LOC loc = f.getLocation(); - if(mesh->StaggerGrids && loc != CELL_CENTRE) { + if(localmesh->StaggerGrids && loc != CELL_CENTRE && loc != CELL_ZLOW) { throw BoutException("neumann_o4 not implemented with staggered grid yet"); } else { @@ -2147,11 +2197,11 @@ void BoundaryNeumann_O4::apply(Field2D &f,BoutReal t) { if(fg) { // Calculate the X and Y normalised values half-way between the guard cell and grid cell - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) // In the guard cell - + mesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) // In the guard cell + + localmesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) // In the guard cell - + mesh->GlobalY(bndry->y - bndry->by) ); // the grid cell + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) // In the guard cell + + localmesh->GlobalY(bndry->y - bndry->by) ); // the grid cell val = fg->generate(xnorm, TWOPI*ynorm, 0.0, t); } @@ -2177,6 +2227,8 @@ void BoundaryNeumann_O4::apply(Field3D &f) { } void BoundaryNeumann_O4::apply(Field3D &f,BoutReal t) { + Mesh* localmesh = f.getMesh(); + bndry->first(); // Decide which generator to use @@ -2188,24 +2240,24 @@ void BoundaryNeumann_O4::apply(Field3D &f,BoutReal t) { // Check for staggered grids CELL_LOC loc = f.getLocation(); - if(mesh->StaggerGrids && loc != CELL_CENTRE) { + if(localmesh->StaggerGrids && loc != CELL_CENTRE && loc != CELL_ZLOW) { throw BoutException("neumann_o4 not implemented with staggered grid yet"); } else { Coordinates *coords = f.getCoordinates(); for(; !bndry->isDone(); bndry->next1d()) { // Calculate the X and Y normalised values half-way between the guard cell and grid cell - BoutReal xnorm = 0.5*( mesh->GlobalX(bndry->x) // In the guard cell - + mesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell + BoutReal xnorm = 0.5*( localmesh->GlobalX(bndry->x) // In the guard cell + + localmesh->GlobalX(bndry->x - bndry->bx) ); // the grid cell - BoutReal ynorm = 0.5*( mesh->GlobalY(bndry->y) // In the guard cell - + mesh->GlobalY(bndry->y - bndry->by) ); // the grid cell + BoutReal ynorm = 0.5*( localmesh->GlobalY(bndry->y) // In the guard cell + + localmesh->GlobalY(bndry->y - bndry->by) ); // the grid cell BoutReal delta = bndry->bx*coords->dx(bndry->x,bndry->y)+bndry->by*coords->dy(bndry->x,bndry->y); - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { if(fg){ - val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(mesh->LocalNz),t); + val = fg->generate(xnorm,TWOPI*ynorm,TWOPI*zk/(localmesh->LocalNz),t); } f(bndry->x,bndry->y, zk) = 12.*delta*val/11. @@ -2232,16 +2284,17 @@ void BoundaryNeumann_O4::apply_ddt(Field2D &f) { } void BoundaryNeumann_O4::apply_ddt(Field3D &f) { + Mesh* localmesh = f.getMesh(); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) + for(int z=0;zLocalNz;z++) (*dt)(bndry->x,bndry->y,z) = 0.; // Set time derivative to zero } /////////////////////////////////////////////////////////////// BoundaryOp* BoundaryNeumann_4thOrder::clone(BoundaryRegion *region, const list &args) { - verifyNumPoints(region,4); + verifyNumPoints(region, 4); if(!args.empty()) { // First argument should be a value val = stringToReal(args.front()); @@ -2262,11 +2315,12 @@ void BoundaryNeumann_4thOrder::apply(Field2D &f) { } void BoundaryNeumann_4thOrder::apply(Field3D &f) { + Mesh* localmesh = f.getMesh(); Coordinates *metric = f.getCoordinates(); // Set (at 4th order) the gradient at the mid-point between the guard cell and the grid cell to be val // This sets the value of the co-ordinate derivative, i.e. DDX/DDY not Grad_par/Grad_perp.x for(bndry->first(); !bndry->isDone(); bndry->next1d()) - for(int z=0;zLocalNz;z++) { + for(int z=0;zLocalNz;z++) { BoutReal delta = -(bndry->bx*metric->dx(bndry->x,bndry->y)+bndry->by*metric->dy(bndry->x,bndry->y)); f(bndry->x,bndry->y,z) = 12.*delta/11.*val + 17./22.*f(bndry->x-bndry->bx,bndry->y-bndry->by,z) + 9./22.*f(bndry->x-2*bndry->bx,bndry->y-2*bndry->by,z) - 5./22.*f(bndry->x-3*bndry->bx,bndry->y-3*bndry->by,z) + 1./22.*f(bndry->x-4*bndry->bx,bndry->y-4*bndry->by,z); f(bndry->x+bndry->bx,bndry->y+bndry->by,z) = -24.*delta*val + 27.*f(bndry->x,bndry->y,z) - 27.*f(bndry->x-bndry->bx,bndry->y-bndry->by,z) + f(bndry->x-2*bndry->bx,bndry->y-2*bndry->by,z); // The f(bndry->x-4*bndry->bx,bndry->y-4*bndry->by,z) term vanishes, so that this sets to zero the 4th order central difference first derivative at the point half way between the guard cell and the grid cell @@ -2280,16 +2334,17 @@ void BoundaryNeumann_4thOrder::apply_ddt(Field2D &f) { } void BoundaryNeumann_4thOrder::apply_ddt(Field3D &f) { + Mesh* localmesh = f.getMesh(); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) + for(int z=0;zLocalNz;z++) (*dt)(bndry->x,bndry->y,z) = 0.; // Set time derivative to zero } /////////////////////////////////////////////////////////////// BoundaryOp* BoundaryNeumannPar::clone(BoundaryRegion *region, const list &args) { - verifyNumPoints(region,1); + verifyNumPoints(region, 1); if(!args.empty()) { output << "WARNING: Ignoring arguments to BoundaryNeumann2\n"; } @@ -2305,16 +2360,17 @@ void BoundaryNeumannPar::apply(Field2D &f) { } void BoundaryNeumannPar::apply(Field3D &f) { + Mesh* localmesh = f.getMesh(); Coordinates *metric = f.getCoordinates(); for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) + for(int z=0;zLocalNz;z++) f(bndry->x,bndry->y,z) = f(bndry->x - bndry->bx,bndry->y - bndry->by,z)*sqrt(metric->g_22(bndry->x, bndry->y)/metric->g_22(bndry->x - bndry->bx, bndry->y - bndry->by)); } /////////////////////////////////////////////////////////////// BoundaryOp* BoundaryRobin::clone(BoundaryRegion *region, const list &args) { - verifyNumPoints(region,1); + verifyNumPoints(region, 1); BoutReal a = 0.5, b = 1.0, g = 0.; list::const_iterator it = args.begin(); @@ -2358,16 +2414,17 @@ void BoundaryRobin::apply(Field2D &f) { } void BoundaryRobin::apply(Field3D &f) { + Mesh* localmesh = f.getMesh(); if(fabs(bval) < 1.e-12) { for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) + for(int z=0;zLocalNz;z++) f(bndry->x, bndry->y, z) = gval / aval; }else { BoutReal sign = 1.; if( (bndry->bx < 0) || (bndry->by < 0)) sign = -1.; for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) + for(int z=0;zLocalNz;z++) f(bndry->x, bndry->y, z) = f(bndry->x - bndry->bx, bndry->y - bndry->by, z) + sign*(gval - aval*f(bndry->x - bndry->bx, bndry->y - bndry->by, z) ) / bval; } } @@ -2381,15 +2438,16 @@ void BoundaryConstGradient::apply(Field2D &f){ } void BoundaryConstGradient::apply(Field3D &f) { + Mesh* localmesh = f.getMesh(); for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) + for(int z=0;zLocalNz;z++) f(bndry->x, bndry->y, z) = 2.*f(bndry->x - bndry->bx, bndry->y - bndry->by, z) - f(bndry->x - 2*bndry->bx,bndry->y - 2*bndry->by,z); } /////////////////////////////////////////////////////////////// BoundaryOp* BoundaryConstGradient::clone(BoundaryRegion *region, const list &args) { - verifyNumPoints(region,2); + verifyNumPoints(region, 2); if(!args.empty()) { output << "WARNING: Ignoring arguments to BoundaryConstGradient\n"; } @@ -2399,7 +2457,7 @@ BoundaryOp* BoundaryConstGradient::clone(BoundaryRegion *region, const list &args) { - verifyNumPoints(region,2); + verifyNumPoints(region, 2); if(!args.empty()) { output << "WARNING: Ignoring arguments to BoundaryZeroLaplace\n"; } @@ -2429,7 +2487,9 @@ void BoundaryZeroLaplace::apply(Field2D &f) { } void BoundaryZeroLaplace::apply(Field3D &f) { - int ncz = mesh->LocalNz; + Mesh* localmesh = f.getMesh(); + + int ncz = localmesh->LocalNz; Coordinates *metric = f.getCoordinates(); @@ -2452,8 +2512,8 @@ void BoundaryZeroLaplace::apply(Field3D &f) { int y = bndry->y; // Take FFT of last 2 points in domain - rfft(f(x - bx, y), mesh->LocalNz, c0.begin()); - rfft(f(x - 2 * bx, y), mesh->LocalNz, c1.begin()); + rfft(f(x - bx, y), localmesh->LocalNz, c0.begin()); + rfft(f(x - 2 * bx, y), localmesh->LocalNz, c1.begin()); c1[0] = c0[0] - c1[0]; // Only need gradient // Solve metric->g11*d2f/dx2 - metric->g33*kz^2f = 0 @@ -2472,7 +2532,7 @@ void BoundaryZeroLaplace::apply(Field3D &f) { c0[jz] *= exp(coef * kwave); // The decaying solution only } // Reverse FFT - irfft(c0.begin(), mesh->LocalNz, f(x, y)); + irfft(c0.begin(), localmesh->LocalNz, f(x, y)); bndry->nextX(); x = bndry->x; @@ -2485,7 +2545,7 @@ void BoundaryZeroLaplace::apply(Field3D &f) { BoundaryOp *BoundaryZeroLaplace2::clone(BoundaryRegion *region, const list &args) { - verifyNumPoints(region, 3); + verifyNumPoints(region, 3); if (!args.empty()) { output << "WARNING: Ignoring arguments to BoundaryZeroLaplace2\n"; } @@ -2519,7 +2579,9 @@ void BoundaryZeroLaplace2::apply(Field2D &f) { } void BoundaryZeroLaplace2::apply(Field3D &f) { - int ncz = mesh->LocalNz; + Mesh* localmesh = f.getMesh(); + + int ncz = localmesh->LocalNz; ASSERT0(ncz % 2 == 0); // Allocation assumes even number @@ -2572,7 +2634,7 @@ void BoundaryZeroLaplace2::apply(Field3D &f) { /////////////////////////////////////////////////////////////// BoundaryOp* BoundaryConstLaplace::clone(BoundaryRegion *region, const list &args) { - verifyNumPoints(region,2); + verifyNumPoints(region, 2); if(!args.empty()) { output << "WARNING: Ignoring arguments to BoundaryConstLaplace\n"; } @@ -2615,9 +2677,11 @@ void BoundaryConstLaplace::apply(Field3D &f) { throw BoutException("ERROR: Can't apply Zero Laplace condition to non-X boundaries\n"); } + Mesh* localmesh = f.getMesh(); + Coordinates *metric = f.getCoordinates(); - int ncz = mesh->LocalNz; + int ncz = localmesh->LocalNz; // Allocate memory Array c0(ncz/2 + 1), c1(ncz/2 + 1), c2(ncz/2 + 1); @@ -2686,21 +2750,23 @@ void BoundaryDivCurl::apply(Vector3D &var) { int jx, jy, jz, jzp, jzm; BoutReal tmp; - Coordinates *metric = mesh->coordinates(var.getLocation()); + Mesh* localmesh = var.x.getMesh(); + + Coordinates *metric = localmesh->coordinates(var.getLocation()); - int ncz = mesh->LocalNz; + int ncz = localmesh->LocalNz; if(bndry->location != BNDRY_XOUT) { throw BoutException("ERROR: DivCurl boundary only works for outer X currently\n"); } var.toCovariant(); - if(mesh->xstart > 2) { + if(localmesh->xstart > 2) { throw BoutException("Error: Div = Curl = 0 boundary condition doesn't work for MXG > 2. Sorry\n"); } - jx = mesh->xend+1; - for(jy=1;jyLocalNy-1;jy++) { + jx = localmesh->xend+1; + for(jy=1;jyLocalNy-1;jy++) { for(jz=0;jzdy(jx-1,jy-1) + metric->dy(jx-1,jy)); var.y(jx,jy,jz) = var.y(jx-2,jy,jz) + (metric->dx(jx-2,jy) + metric->dx(jx-1,jy)) * tmp; - if(mesh->xstart == 2) + if(localmesh->xstart == 2) // 4th order to get last point var.y(jx+1,jy,jz) = var.y(jx-3,jy,jz) + 4.*metric->dx(jx,jy)*tmp; @@ -2720,7 +2786,7 @@ void BoundaryDivCurl::apply(Vector3D &var) { tmp = (var.x(jx-1,jy,jzp) - var.x(jx-1,jy,jzm)) / (2.*metric->dz); var.z(jx,jy,jz) = var.z(jx-2,jy,jz) + (metric->dx(jx-2,jy) + metric->dx(jx-1,jy)) * tmp; - if(mesh->xstart == 2) + if(localmesh->xstart == 2) var.z(jx+1,jy,jz) = var.z(jx-3,jy,jz) + 4.*metric->dx(jx,jy)*tmp; // d/dx( Jmetric->g11 B_x ) = - d/dx( Jmetric->g12 B_y + Jmetric->g13 B_z) @@ -2739,7 +2805,7 @@ void BoundaryDivCurl::apply(Vector3D &var) { var.x(jx,jy,jz) = ( metric->J(jx-2,jy)*metric->g11(jx-2,jy)*var.x(jx-2,jy,jz) + (metric->dx(jx-2,jy) + metric->dx(jx-1,jy)) * tmp ) / metric->J(jx,jy)*metric->g11(jx,jy); - if(mesh->xstart == 2) + if(localmesh->xstart == 2) var.x(jx+1,jy,jz) = ( metric->J(jx-3,jy)*metric->g11(jx-3,jy)*var.x(jx-3,jy,jz) + 4.*metric->dx(jx,jy)*tmp ) / metric->J(jx+1,jy)*metric->g11(jx+1,jy); } @@ -2780,7 +2846,7 @@ void BoundaryFree::apply_ddt(Field3D &UNUSED(f)) { // 2nd order extrapolation: BoundaryOp* BoundaryFree_O2::clone(BoundaryRegion *region, const list &args){ - verifyNumPoints(region,2); + verifyNumPoints(region, 2); if(!args.empty()) { output << "WARNING: Ignoring arguments to BoundaryFree\n"; } @@ -2791,12 +2857,14 @@ void BoundaryFree_O2::apply(Field2D &f) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* localmesh = f.getMesh(); + bndry->first(); // Check for staggered grids CELL_LOC loc = f.getLocation(); - if(mesh->StaggerGrids && loc != CELL_CENTRE) { + if(localmesh->StaggerGrids && loc != CELL_CENTRE && loc != CELL_ZLOW) { // Staggered. Need to apply slightly differently if( loc == CELL_XLOW) { @@ -2869,6 +2937,8 @@ void BoundaryFree_O2::apply(Field2D &f) { } } } + } else { + throw BoutException("Unhandled staggering in BoundaryFree_O2."); } } else { @@ -2890,11 +2960,12 @@ void BoundaryFree_O2::apply(Field3D &f) { bndry->first(); + Mesh* localmesh = f.getMesh(); // Check for staggered grids CELL_LOC loc = f.getLocation(); - if(mesh->StaggerGrids && loc != CELL_CENTRE) { + if(localmesh->StaggerGrids && loc != CELL_CENTRE && loc != CELL_ZLOW) { // Staggered. Need to apply slightly differently if( loc == CELL_XLOW ) { @@ -2905,7 +2976,7 @@ void BoundaryFree_O2::apply(Field3D &f) { for(; !bndry->isDone(); bndry->next1d()) { - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { for(int i=0;iwidth;i++) { int xi = bndry->x + i*bndry->bx; int yi = bndry->y + i*bndry->by; @@ -2918,7 +2989,7 @@ void BoundaryFree_O2::apply(Field3D &f) { // Inner x boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { for(int i=-1;iwidth;i++) { int xi = bndry->x + i*bndry->bx; int yi = bndry->y + i*bndry->by; @@ -2932,7 +3003,7 @@ void BoundaryFree_O2::apply(Field3D &f) { for(; !bndry->isDone(); bndry->next1d()) { - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { for(int i=0;iwidth;i++) { int xi = bndry->x + i*bndry->bx; int yi = bndry->y + i*bndry->by; @@ -2948,7 +3019,7 @@ void BoundaryFree_O2::apply(Field3D &f) { if(bndry->by > 0) { // Upper y boundary for(; !bndry->isDone(); bndry->next1d()) { - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { for(int i=0;iwidth;i++) { int xi = bndry->x + i*bndry->bx; int yi = bndry->y + i*bndry->by; @@ -2961,7 +3032,7 @@ void BoundaryFree_O2::apply(Field3D &f) { // Lower y boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { for(int i=-1;iwidth;i++) { int xi = bndry->x + i*bndry->bx; int yi = bndry->y + i*bndry->by; @@ -2974,7 +3045,7 @@ void BoundaryFree_O2::apply(Field3D &f) { // x boundaries for(; !bndry->isDone(); bndry->next1d()) { - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { for(int i=0;iwidth;i++) { int xi = bndry->x + i*bndry->bx; int yi = bndry->y + i*bndry->by; @@ -2983,13 +3054,15 @@ void BoundaryFree_O2::apply(Field3D &f) { } } } + } else { + throw BoutException("Unhandled staggering in BoundaryFree_O2."); } } else { // Standard (non-staggered) case for(; !bndry->isDone(); bndry->next1d()) { - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { for(int i=0;iwidth;i++) { int xi = bndry->x + i*bndry->bx; int yi = bndry->y + i*bndry->by; @@ -3007,9 +3080,10 @@ void BoundaryFree_O2::apply_ddt(Field2D &f) { } void BoundaryFree_O2::apply_ddt(Field3D &f) { + Mesh* localmesh = f.getMesh(); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) + for(int z=0;zLocalNz;z++) (*dt)(bndry->x,bndry->y,z) = 0.; // Set time derivative to zero } @@ -3018,7 +3092,7 @@ void BoundaryFree_O2::apply_ddt(Field3D &f) { // Third order extrapolation: ////////////////////////////////// BoundaryOp* BoundaryFree_O3::clone(BoundaryRegion *region, const list &args){ - verifyNumPoints(region,3); + verifyNumPoints(region, 3); if(!args.empty()) { output << "WARNING: Ignoring arguments to BoundaryConstLaplace\n"; @@ -3028,12 +3102,14 @@ BoundaryOp* BoundaryFree_O3::clone(BoundaryRegion *region, const list &a void BoundaryFree_O3::apply(Field2D &f) { + Mesh* localmesh = f.getMesh(); + bndry->first(); // Check for staggered grids CELL_LOC loc = f.getLocation(); - if(mesh->StaggerGrids && loc != CELL_CENTRE) { + if(localmesh->StaggerGrids && loc != CELL_CENTRE && loc != CELL_ZLOW) { // Staggered. Need to apply slightly differently if( loc == CELL_XLOW) { @@ -3107,6 +3183,8 @@ void BoundaryFree_O3::apply(Field2D &f) { } } } + } else { + throw BoutException("Unhandled staggering in BoundaryFree_O3."); } } else { @@ -3126,13 +3204,14 @@ void BoundaryFree_O3::apply(Field2D &f) { void BoundaryFree_O3::apply(Field3D &f) { // Extrapolate from the last evolved simulation cells into the guard cells at 3rd order. - bndry->first(); + Mesh* localmesh = f.getMesh(); + bndry->first(); // Check for staggered grids CELL_LOC loc = f.getLocation(); - if(mesh->StaggerGrids && loc != CELL_CENTRE) { + if(localmesh->StaggerGrids && loc != CELL_CENTRE && loc != CELL_ZLOW) { // Staggered. Need to apply slightly differently if( loc == CELL_XLOW ) { @@ -3143,7 +3222,7 @@ void BoundaryFree_O3::apply(Field3D &f) { for(; !bndry->isDone(); bndry->next1d()) { - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { for(int i=0;iwidth;i++) { int xi = bndry->x + i*bndry->bx; int yi = bndry->y + i*bndry->by; @@ -3157,7 +3236,7 @@ void BoundaryFree_O3::apply(Field3D &f) { // Inner x boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { for(int i=-1;iwidth;i++) { int xi = bndry->x + i*bndry->bx; int yi = bndry->y + i*bndry->by; @@ -3172,7 +3251,7 @@ void BoundaryFree_O3::apply(Field3D &f) { for(; !bndry->isDone(); bndry->next1d()) { - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { for(int i=0;iwidth;i++) { int xi = bndry->x + i*bndry->bx; int yi = bndry->y + i*bndry->by; @@ -3189,7 +3268,7 @@ void BoundaryFree_O3::apply(Field3D &f) { if(bndry->by > 0) { // Upper y boundary for(; !bndry->isDone(); bndry->next1d()) { - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { for(int i=0;iwidth;i++) { int xi = bndry->x + i*bndry->bx; int yi = bndry->y + i*bndry->by; @@ -3203,7 +3282,7 @@ void BoundaryFree_O3::apply(Field3D &f) { // Lower y boundary. Set one point inwards for(; !bndry->isDone(); bndry->next1d()) { - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { for(int i=-1;iwidth;i++) { int xi = bndry->x + i*bndry->bx; int yi = bndry->y + i*bndry->by; @@ -3217,7 +3296,7 @@ void BoundaryFree_O3::apply(Field3D &f) { // x boundaries for(; !bndry->isDone(); bndry->next1d()) { - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { for(int i=0;iwidth;i++) { int xi = bndry->x + i*bndry->bx; int yi = bndry->y + i*bndry->by; @@ -3228,13 +3307,15 @@ void BoundaryFree_O3::apply(Field3D &f) { } } } + } else { + throw BoutException("Unhandled staggering in BoundaryFree_O3."); } } else { // Standard (non-staggered) case for(; !bndry->isDone(); bndry->next1d()) { - for(int zk=0;zkLocalNz;zk++) { + for(int zk=0;zkLocalNz;zk++) { for(int i=0;iwidth;i++) { int xi = bndry->x + i*bndry->bx; int yi = bndry->y + i*bndry->by; @@ -3253,9 +3334,10 @@ void BoundaryFree_O3::apply_ddt(Field2D &f) { } void BoundaryFree_O3::apply_ddt(Field3D &f) { + Mesh* localmesh = f.getMesh(); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) + for(int z=0;zLocalNz;z++) (*dt)(bndry->x,bndry->y,z) = 0.; // Set time derivative to zero } @@ -3304,13 +3386,15 @@ void BoundaryRelax::apply_ddt(Field2D &f) { void BoundaryRelax::apply_ddt(Field3D &f) { TRACE("BoundaryRelax::apply_ddt(Field3D)"); + Mesh* localmesh = f.getMesh(); + // Make a copy of f Field3D g = f; // NOTE: This is not very efficient... copying entire field // Apply the boundary to g op->apply(g); // Set time-derivatives for(bndry->first(); !bndry->isDone(); bndry->next()) - for(int z=0;zLocalNz;z++) { + for(int z=0;zLocalNz;z++) { ddt(f)(bndry->x, bndry->y, z) = r * (g(bndry->x, bndry->y, z) - f(bndry->x, bndry->y, z)); } } @@ -3377,14 +3461,16 @@ void BoundaryToFieldAligned::apply(Field2D &f, BoutReal t) { } void BoundaryToFieldAligned::apply(Field3D &f, BoutReal t) { + Mesh* localmesh = f.getMesh(); + //NOTE: This is not very efficient... updating entire field - f = mesh->fromFieldAligned(f); + f = localmesh->fromFieldAligned(f); // Apply the boundary to shifted field op->apply(f, t); //Shift back - f = mesh->toFieldAligned(f); + f = localmesh->toFieldAligned(f); //This is inefficient -- could instead use the shiftZ just in the bndry //but this is not portable to other parallel transforms -- we could instead @@ -3396,10 +3482,12 @@ void BoundaryToFieldAligned::apply_ddt(Field2D &f) { } void BoundaryToFieldAligned::apply_ddt(Field3D &f) { - f = mesh->fromFieldAligned(f); - ddt(f) = mesh->fromFieldAligned(ddt(f)); + Mesh* localmesh = f.getMesh(); + + f = localmesh->fromFieldAligned(f); + ddt(f) = localmesh->fromFieldAligned(ddt(f)); op->apply_ddt(f); - ddt(f) = mesh->toFieldAligned(ddt(f)); + ddt(f) = localmesh->toFieldAligned(ddt(f)); } @@ -3420,14 +3508,16 @@ void BoundaryFromFieldAligned::apply(Field2D &f, BoutReal t) { } void BoundaryFromFieldAligned::apply(Field3D &f, BoutReal t) { + Mesh* localmesh = f.getMesh(); + //NOTE: This is not very efficient... shifting entire field - f = mesh->toFieldAligned(f); + f = localmesh->toFieldAligned(f); // Apply the boundary to shifted field op->apply(f, t); //Shift back - f = mesh->fromFieldAligned(f); + f = localmesh->fromFieldAligned(f); //This is inefficient -- could instead use the shiftZ just in the bndry //but this is not portable to other parallel transforms -- we could instead @@ -3439,8 +3529,10 @@ void BoundaryFromFieldAligned::apply_ddt(Field2D &f) { } void BoundaryFromFieldAligned::apply_ddt(Field3D &f) { - f = mesh->toFieldAligned(f); - ddt(f) = mesh->toFieldAligned(ddt(f)); + Mesh* localmesh = f.getMesh(); + + f = localmesh->toFieldAligned(f); + ddt(f) = localmesh->toFieldAligned(ddt(f)); op->apply_ddt(f); - ddt(f) = mesh->fromFieldAligned(ddt(f)); + ddt(f) = localmesh->fromFieldAligned(ddt(f)); } diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 10dcd9b4df..a67e93dc54 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -28,23 +28,23 @@ Coordinates::Coordinates(Mesh *mesh) G3_23(mesh), G1(mesh), G2(mesh), G3(mesh), ShiftTorsion(mesh), IntShiftTorsion(mesh), localmesh(mesh), location(CELL_CENTRE) { - if (mesh->get(dx, "dx")) { + if (localmesh->get(dx, "dx")) { output_warn.write("\tWARNING: differencing quantity 'dx' not found. Set to 1.0\n"); dx = 1.0; } - if (mesh->periodicX) { - mesh->communicate(dx); + if (localmesh->periodicX) { + localmesh->communicate(dx); } - if (mesh->get(dy, "dy")) { + if (localmesh->get(dy, "dy")) { output_warn.write("\tWARNING: differencing quantity 'dy' not found. Set to 1.0\n"); dy = 1.0; } - nz = mesh->LocalNz; + nz = localmesh->LocalNz; - if (mesh->get(dz, "dz")) { + if (localmesh->get(dz, "dz")) { // Couldn't read dz from input int zperiod; BoutReal ZMIN, ZMAX; @@ -64,14 +64,14 @@ Coordinates::Coordinates(Mesh *mesh) } // Diagonal components of metric tensor g^{ij} (default to 1) - mesh->get(g11, "g11", 1.0); - mesh->get(g22, "g22", 1.0); - mesh->get(g33, "g33", 1.0); + localmesh->get(g11, "g11", 1.0); + localmesh->get(g22, "g22", 1.0); + localmesh->get(g33, "g33", 1.0); // Off-diagonal elements. Default to 0 - mesh->get(g12, "g12", 0.0); - mesh->get(g13, "g13", 0.0); - mesh->get(g23, "g23", 0.0); + localmesh->get(g12, "g12", 0.0); + localmesh->get(g13, "g13", 0.0); + localmesh->get(g23, "g23", 0.0); // Check input metrics if ((!finite(g11)) || (!finite(g22)) || (!finite(g33))) { @@ -86,19 +86,19 @@ Coordinates::Coordinates(Mesh *mesh) /// Find covariant metric components // Check if any of the components are present - if (mesh->sourceHasVar("g_11") or mesh->sourceHasVar("g_22") or - mesh->sourceHasVar("g_33") or mesh->sourceHasVar("g_12") or - mesh->sourceHasVar("g_13") or mesh->sourceHasVar("g_23")) { + if (localmesh->sourceHasVar("g_11") or localmesh->sourceHasVar("g_22") or + localmesh->sourceHasVar("g_33") or localmesh->sourceHasVar("g_12") or + localmesh->sourceHasVar("g_13") or localmesh->sourceHasVar("g_23")) { // Check that all components are present - if (mesh->sourceHasVar("g_11") and mesh->sourceHasVar("g_22") and - mesh->sourceHasVar("g_33") and mesh->sourceHasVar("g_12") and - mesh->sourceHasVar("g_13") and mesh->sourceHasVar("g_23")) { - mesh->get(g_11, "g_11"); - mesh->get(g_22, "g_22"); - mesh->get(g_33, "g_33"); - mesh->get(g_12, "g_12"); - mesh->get(g_13, "g_13"); - mesh->get(g_23, "g_23"); + if (localmesh->sourceHasVar("g_11") and localmesh->sourceHasVar("g_22") and + localmesh->sourceHasVar("g_33") and localmesh->sourceHasVar("g_12") and + localmesh->sourceHasVar("g_13") and localmesh->sourceHasVar("g_23")) { + localmesh->get(g_11, "g_11"); + localmesh->get(g_22, "g_22"); + localmesh->get(g_33, "g_33"); + localmesh->get(g_12, "g_12"); + localmesh->get(g_13, "g_13"); + localmesh->get(g_23, "g_23"); output_warn.write("\tWARNING! Covariant components of metric tensor set manually. " "Contravariant components NOT recalculated\n"); @@ -124,7 +124,7 @@ Coordinates::Coordinates(Mesh *mesh) // Attempt to read J from the grid file Field2D Jcalc = J; - if (mesh->get(J, "J")) { + if (localmesh->get(J, "J")) { output_warn.write("\tWARNING: Jacobian 'J' not found. Calculating from metric tensor\n"); J = Jcalc; } else { @@ -137,7 +137,7 @@ Coordinates::Coordinates(Mesh *mesh) // Attempt to read Bxy from the grid file Field2D Bcalc = Bxy; - if (mesh->get(Bxy, "Bxy")) { + if (localmesh->get(Bxy, "Bxy")) { output_warn.write("\tWARNING: Magnitude of B field 'Bxy' not found. Calculating from " "metric tensor\n"); Bxy = Bcalc; @@ -155,15 +155,15 @@ Coordinates::Coordinates(Mesh *mesh) throw BoutException("Differential geometry failed\n"); } - if (mesh->get(ShiftTorsion, "ShiftTorsion")) { + if (localmesh->get(ShiftTorsion, "ShiftTorsion")) { output_warn.write("\tWARNING: No Torsion specified for zShift. Derivatives may not be correct\n"); ShiftTorsion = 0.0; } ////////////////////////////////////////////////////// - if (mesh->IncIntShear) { - if (mesh->get(IntShiftTorsion, "IntShiftTorsion")) { + if (localmesh->IncIntShear) { + if (localmesh->get(IntShiftTorsion, "IntShiftTorsion")) { output_warn.write("\tWARNING: No Integrated torsion specified\n"); IntShiftTorsion = 0.0; } @@ -174,32 +174,73 @@ Coordinates::Coordinates(Mesh *mesh) namespace { /// Interpolate a Field2D to a new CELL_LOC with interp_to. /// Communicates to set internal guard cells. - /// Boundary guard cells are set equal to the nearest grid point (equivalent to - /// 2nd order accurate Neumann boundary condition). + /// Boundary guard cells are set by extrapolating from the grid, like + /// 'free_o3' boundary conditions /// Corner guard cells are set to BoutNaN - Field2D interpolateAndNeumann(const Field2D &f, CELL_LOC location) { + Field2D interpolateAndExtrapolate(const Field2D &f, CELL_LOC location) { 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(result); - // Copy nearest value into boundaries so that differential geometry terms can - // be interpolated if necessary - // Note: cannot use applyBoundary("neumann") here because applyBoundary() + // Extrapolate into boundaries 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 + // 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 (bndry->bx != 0) { - // If bx!=0 we are on an x-boundary, inner if bx>0 and outer if bx<0 - for(bndry->first(); !bndry->isDone(); bndry->next1d()) { - for (int i=0; ixstart; i++) - result(bndry->x+i*bndry->bx,bndry->y) = result(bndry->x+(i-1)*bndry->bx, bndry->y-bndry->by); + int extrap_start = 0; + 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 + if (extrap_start>0) { + // 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.; } - } - if (bndry->by != 0) { - // If by!=0 we are on a y-boundary, upper if by>0 and lower if by<0 - for(bndry->first(); !bndry->isDone(); bndry->next1d()) { - for (int i=0; iystart; i++) - result(bndry->x,bndry->y+i*bndry->by) = result(bndry->x-bndry->bx, bndry->y+(i-1)*bndry->by); + + // 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;iwidth;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); + } + } else { + // not enough grid points to extrapolate, set equal to last grid point + for(int i=extrap_start;iwidth;i++) { + result(bndry->x + i*bndry->bx, bndry->y + i*bndry->by) = result(bndry->x - bndry->bx, bndry->y - bndry->by); + } } } } @@ -229,22 +270,22 @@ Coordinates::Coordinates(Mesh *mesh, const CELL_LOC loc, const Coordinates* coor G3_23(mesh), G1(mesh), G2(mesh), G3(mesh), ShiftTorsion(mesh), IntShiftTorsion(mesh), localmesh(mesh), location(loc) { - dx = interpolateAndNeumann(coords_in->dx, location); - dy = interpolateAndNeumann(coords_in->dy, location); + dx = interpolateAndExtrapolate(coords_in->dx, location); + dy = interpolateAndExtrapolate(coords_in->dy, location); - nz = mesh->LocalNz; + nz = localmesh->LocalNz; dz = coords_in->dz; // Diagonal components of metric tensor g^{ij} - g11 = interpolateAndNeumann(coords_in->g11, location); - g22 = interpolateAndNeumann(coords_in->g22, location); - g33 = interpolateAndNeumann(coords_in->g33, location); + g11 = interpolateAndExtrapolate(coords_in->g11, location); + g22 = interpolateAndExtrapolate(coords_in->g22, location); + g33 = interpolateAndExtrapolate(coords_in->g33, location); // Off-diagonal elements. - g12 = interpolateAndNeumann(coords_in->g12, location); - g13 = interpolateAndNeumann(coords_in->g13, location); - g23 = interpolateAndNeumann(coords_in->g23, location); + g12 = interpolateAndExtrapolate(coords_in->g12, location); + g13 = interpolateAndExtrapolate(coords_in->g13, location); + g23 = interpolateAndExtrapolate(coords_in->g23, location); // Check input metrics if ((!finite(g11, RGN_NOBNDRY)) || (!finite(g22, RGN_NOBNDRY)) || (!finite(g33, RGN_NOBNDRY))) { @@ -273,12 +314,12 @@ Coordinates::Coordinates(Mesh *mesh, const CELL_LOC loc, const Coordinates* coor throw BoutException("Differential geometry failed\n"); } - ShiftTorsion = interpolateAndNeumann(coords_in->ShiftTorsion, location); + ShiftTorsion = interpolateAndExtrapolate(coords_in->ShiftTorsion, location); ////////////////////////////////////////////////////// - if (mesh->IncIntShear) { - IntShiftTorsion = interpolateAndNeumann(coords_in->IntShiftTorsion, location); + if (localmesh->IncIntShear) { + IntShiftTorsion = interpolateAndExtrapolate(coords_in->IntShiftTorsion, location); } } @@ -648,11 +689,11 @@ const Field2D Coordinates::DDZ(const Field2D &f, CELL_LOC loc, // Parallel gradient const Field2D Coordinates::Grad_par(const Field2D &var, CELL_LOC outloc, - DIFF_METHOD UNUSED(method)) { + DIFF_METHOD method) { TRACE("Coordinates::Grad_par( Field2D )"); - ASSERT1(location == outloc || outloc == CELL_DEFAULT); + ASSERT1(location == outloc || (outloc == CELL_DEFAULT && location == var.getLocation())); - return DDY(var) / sqrt(g_22); + return DDY(var, outloc, method) / sqrt(g_22); } const Field3D Coordinates::Grad_par(const Field3D &var, CELL_LOC outloc, @@ -669,9 +710,9 @@ const Field3D Coordinates::Grad_par(const Field3D &var, CELL_LOC outloc, const Field2D Coordinates::Vpar_Grad_par(const Field2D &v, const Field2D &f, CELL_LOC outloc, - DIFF_METHOD UNUSED(method)) { - ASSERT1(location == outloc || outloc == CELL_DEFAULT); - return VDDY(v, f) / sqrt(g_22); + DIFF_METHOD method) { + ASSERT1(location == outloc || (outloc == CELL_DEFAULT && location == f.getLocation())); + return VDDY(v, f, outloc, method) / sqrt(g_22); } const Field3D Coordinates::Vpar_Grad_par(const Field3D &v, const Field3D &f, CELL_LOC outloc, @@ -728,43 +769,37 @@ const Field3D Coordinates::Div_par(const Field3D &f, CELL_LOC outloc, // second parallel derivative (b dot Grad)(b dot Grad) // Note: For parallel Laplacian use Laplace_par -const Field2D Coordinates::Grad2_par2(const Field2D &f, CELL_LOC outloc) { +const Field2D Coordinates::Grad2_par2(const Field2D &f, CELL_LOC outloc, DIFF_METHOD method) { TRACE("Coordinates::Grad2_par2( Field2D )"); - ASSERT1(location == outloc || outloc == CELL_DEFAULT); + ASSERT1(location == outloc || (outloc == CELL_DEFAULT && location == f.getLocation())); Field2D sg = sqrt(g_22); - Field2D result = DDY(1. / sg, outloc) * DDY(f, outloc) / sg + D2DY2(f, outloc) / g_22; + Field2D result = DDY(1. / sg, outloc, method) * DDY(f, outloc, method) / sg + D2DY2(f, outloc, method) / g_22; return result; } -const Field3D Coordinates::Grad2_par2(const Field3D &f, CELL_LOC outloc) { +const Field3D Coordinates::Grad2_par2(const Field3D &f, CELL_LOC outloc, DIFF_METHOD method) { TRACE("Coordinates::Grad2_par2( Field3D )"); - ASSERT1(location == outloc || outloc == CELL_DEFAULT); + if (outloc == CELL_DEFAULT) { + outloc = f.getLocation(); + } + ASSERT1(location == outloc); Field2D sg(localmesh); Field3D result(localmesh), r2(localmesh); sg = sqrt(g_22); - sg = DDY(1. / sg) / sg; + sg = DDY(1. / sg, outloc, method) / sg; - if (outloc == CELL_DEFAULT) { - outloc = f.getLocation(); - } - if (sg.getLocation() != outloc) { - localmesh->communicate(sg); - sg = interp_to(sg, outloc); - } - - result = ::DDY(f, outloc); + result = ::DDY(f, outloc, method); - r2 = D2DY2(f, outloc) / interp_to(g_22, outloc); + r2 = D2DY2(f, outloc, method) / g_22; result = sg * result + r2; - ASSERT2(((outloc == CELL_DEFAULT) && (result.getLocation() == f.getLocation())) || - (result.getLocation() == outloc)); + ASSERT2(result.getLocation() == outloc); return result; } @@ -785,7 +820,10 @@ const Field2D Coordinates::Delp2(const Field2D &f, CELL_LOC outloc) { const Field3D Coordinates::Delp2(const Field3D &f, CELL_LOC outloc) { TRACE("Coordinates::Delp2( Field3D )"); - ASSERT1(location == outloc || outloc == CELL_DEFAULT); + if (outloc == CELL_DEFAULT) { + outloc = f.getLocation(); + } + ASSERT1(location == outloc); if (localmesh->GlobalNx == 1 && localmesh->GlobalNz == 1) { // copy mesh, location, etc @@ -793,7 +831,6 @@ const Field3D Coordinates::Delp2(const Field3D &f, CELL_LOC outloc) { } ASSERT2(localmesh->xstart > 0); // Need at least one guard cell - if (outloc == CELL_DEFAULT) outloc = f.getLocation(); ASSERT2(f.getLocation() == outloc); Field3D result(localmesh); @@ -855,7 +892,7 @@ const FieldPerp Coordinates::Delp2(const FieldPerp &f, CELL_LOC outloc) { if (outloc == CELL_DEFAULT) outloc = f.getLocation(); - ASSERT1(location == outloc || outloc == CELL_DEFAULT); + ASSERT1(location == outloc); ASSERT2(f.getLocation() == outloc); FieldPerp result(localmesh); diff --git a/src/mesh/data/gridfromfile.cxx b/src/mesh/data/gridfromfile.cxx index 2a2758525e..670a7e29c8 100644 --- a/src/mesh/data/gridfromfile.cxx +++ b/src/mesh/data/gridfromfile.cxx @@ -367,7 +367,6 @@ bool GridFile::get(Mesh *UNUSED(m), vector &var, const string &name, return true; } - ///////////////////////////////////////////////////////////// // Private routines diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 58625c77d6..2728c1afb2 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -167,15 +167,15 @@ const Field3D Grad_parP(const Field3D &apar, const Field3D &f) { * vparallel times the parallel derivative along unperturbed B-field *******************************************************************************/ -const Field2D Vpar_Grad_par(const Field2D &v, const Field2D &f, const CELL_LOC outloc) { - return f.getCoordinates(outloc)->Vpar_Grad_par(v, f, outloc); +const Field2D Vpar_Grad_par(const Field2D &v, const Field2D &f, CELL_LOC outloc, DIFF_METHOD method) { + return f.getCoordinates(outloc)->Vpar_Grad_par(v, f, outloc, method); } -const Field3D Vpar_Grad_par(const Field3D &v, const Field3D &f, CELL_LOC outloc, DIFF_METHOD method) { +const Field3D Vpar_Grad_par(const Field3D &v, const Field3D &f, const CELL_LOC outloc, const DIFF_METHOD method) { return f.getCoordinates(outloc)->Vpar_Grad_par(v, f, outloc, method); } -const Field3D Vpar_Grad_par(const Field3D &v, const Field3D &f, DIFF_METHOD method, CELL_LOC outloc) { +const Field3D Vpar_Grad_par(const Field3D &v, const Field3D &f, const DIFF_METHOD method, const CELL_LOC outloc) { return f.getCoordinates(outloc)->Vpar_Grad_par(v, f, outloc, method); } @@ -184,8 +184,8 @@ const Field3D Vpar_Grad_par(const Field3D &v, const Field3D &f, DIFF_METHOD meth * parallel divergence operator B \partial_{||} (F/B) *******************************************************************************/ -const Field2D Div_par(const Field2D &f, const CELL_LOC outloc) { - return f.getCoordinates(outloc)->Div_par(f, outloc); +const Field2D Div_par(const Field2D &f, CELL_LOC outloc, DIFF_METHOD method) { + return f.getCoordinates(outloc)->Div_par(f, outloc, method); } const Field3D Div_par(const Field3D &f, CELL_LOC outloc, DIFF_METHOD method) { @@ -246,6 +246,11 @@ const Field3D Div_par_flux(const Field3D &v, const Field3D &f, DIFF_METHOD metho return Div_par_flux(v,f, outloc, method); } +const Field2D Div_par_flux(const Field2D &v, const Field2D &f, CELL_LOC outloc, DIFF_METHOD method) { + Coordinates *metric = v.getCoordinates(outloc); + return metric->Bxy*FDDY(v, f/f.getCoordinates()->Bxy, outloc, method)/sqrt(metric->g_22); +} + /******************************************************************************* * Parallel derivatives converting between left and cell centred * NOTE: These are a quick hack to test if this works. The whole staggered grid @@ -342,58 +347,11 @@ const Field3D Vpar_Grad_par_LCtoC(const Field3D &v, const Field3D &f, REGION reg } } } - else if (vUseUpDown) { - // Only v has up/down fields - // f must shift to field aligned coordinates - Field3D f_fa = vMesh->toFieldAligned(f); - - BOUT_OMP(parallel) { - stencil fval, vval; - BOUT_FOR_INNER(i, vMesh->getRegion3D(region_str)) { - fval.mm = f_fa[i.ymm()]; - fval.m = f_fa[i.ym()]; - fval.c = f_fa[i]; - fval.p = f_fa[i.yp()]; - fval.pp = f_fa[i.ypp()]; - - vval.m = v.ydown()[i.ym()]; - vval.c = v[i]; - vval.p = v.yup()[i.yp()]; - - // Left side - result[i] = (vval.c >= 0.0) ? vval.c * fval.m : vval.c * fval.c; - // Right side - result[i] -= (vval.p >= 0.0) ? vval.p * fval.c : vval.p * fval.p; - } - } - } - else if (fUseUpDown) { - // Only f has up/down fields - // v must shift to field aligned coordinates - Field3D v_fa = vMesh->toFieldAligned(v); - - BOUT_OMP(parallel) { - stencil fval, vval; - BOUT_FOR_INNER(i, vMesh->getRegion3D(region_str)) { - fval.m = f.ydown()[i.ym()]; - fval.c = f[i]; - fval.p = f.yup()[i.yp()]; - - vval.mm = v_fa[i.ymm()]; - vval.m = v_fa[i.ym()]; - vval.c = v_fa[i]; - vval.p = v_fa[i.yp()]; - vval.pp = v_fa[i.ypp()]; - - // Left side - result[i] = (vval.c >= 0.0) ? vval.c * fval.m : vval.c * fval.c; - // Right side - result[i] -= (vval.p >= 0.0) ? vval.p * fval.c : vval.p * fval.p; - } - } - } else { // Both must shift to field aligned + // (even if one of v and f has yup/ydown fields, it doesn't make sense to + // multiply them with one in field-aligned and one in non-field-aligned + // coordinates) Field3D v_fa = vMesh->toFieldAligned(v); Field3D f_fa = vMesh->toFieldAligned(f); @@ -413,6 +371,8 @@ const Field3D Vpar_Grad_par_LCtoC(const Field3D &v, const Field3D &f, REGION reg // Right side result[i] -= (vval.p >= 0.0) ? vval.p * fval.c : vval.p * fval.p; } + + result = vMesh->fromFieldAligned(result); } } @@ -534,12 +494,12 @@ const Field3D Div_par_CtoL(const Field3D &var) { * Note: For parallel Laplacian use LaplacePar *******************************************************************************/ -const Field2D Grad2_par2(const Field2D &f, const CELL_LOC outloc) { - return f.getCoordinates(outloc)->Grad2_par2(f, outloc); +const Field2D Grad2_par2(const Field2D &f, CELL_LOC outloc, DIFF_METHOD method) { + return f.getCoordinates(outloc)->Grad2_par2(f, outloc, method); } -const Field3D Grad2_par2(const Field3D &f, const CELL_LOC outloc) { - return f.getCoordinates(outloc)->Grad2_par2(f, outloc); +const Field3D Grad2_par2(const Field3D &f, CELL_LOC outloc, DIFF_METHOD method) { + return f.getCoordinates(outloc)->Grad2_par2(f, outloc, method); } /******************************************************************************* @@ -547,28 +507,28 @@ const Field3D Grad2_par2(const Field3D &f, const CELL_LOC outloc) { * Parallel divergence of diffusive flux, K*Grad_par *******************************************************************************/ -const Field2D Div_par_K_Grad_par(BoutReal kY, const Field2D &f, const CELL_LOC outloc) { +const Field2D Div_par_K_Grad_par(BoutReal kY, const Field2D &f, CELL_LOC outloc) { return kY*Grad2_par2(f, outloc); } -const Field3D Div_par_K_Grad_par(BoutReal kY, const Field3D &f, const CELL_LOC outloc) { +const Field3D Div_par_K_Grad_par(BoutReal kY, const Field3D &f, CELL_LOC outloc) { return kY*Grad2_par2(f, outloc); } -const Field2D Div_par_K_Grad_par(const Field2D &kY, const Field2D &f, const CELL_LOC outloc) { - return interp_to(kY, outloc)*Grad2_par2(f, outloc) + Div_par(kY, outloc)*Grad_par(f, outloc); +const Field2D Div_par_K_Grad_par(const Field2D &kY, const Field2D &f, CELL_LOC outloc, DIFF_METHOD method) { + return interp_to(kY, outloc)*Grad2_par2(f, outloc, method) + Div_par(kY, outloc, method)*Grad_par(f, outloc, method); } -const Field3D Div_par_K_Grad_par(const Field2D &kY, const Field3D &f, const CELL_LOC outloc) { - return interp_to(kY, outloc)*Grad2_par2(f, outloc) + Div_par(kY, outloc)*Grad_par(f, outloc); +const Field3D Div_par_K_Grad_par(const Field2D &kY, const Field3D &f, CELL_LOC outloc, DIFF_METHOD method) { + return interp_to(kY, outloc)*Grad2_par2(f, outloc, method) + Div_par(kY, outloc, method)*Grad_par(f, outloc, method); } -const Field3D Div_par_K_Grad_par(const Field3D &kY, const Field2D &f, const CELL_LOC outloc) { - return interp_to(kY, outloc)*Grad2_par2(f, outloc) + Div_par(kY, outloc)*Grad_par(f, outloc); +const Field3D Div_par_K_Grad_par(const Field3D &kY, const Field2D &f, CELL_LOC outloc, DIFF_METHOD method) { + return interp_to(kY, outloc)*Grad2_par2(f, outloc, method) + Div_par(kY, outloc, method)*Grad_par(f, outloc, method); } -const Field3D Div_par_K_Grad_par(const Field3D &kY, const Field3D &f, const CELL_LOC outloc) { - return interp_to(kY, outloc)*Grad2_par2(f, outloc) + Div_par(kY, outloc)*Grad_par(f, outloc); +const Field3D Div_par_K_Grad_par(const Field3D &kY, const Field3D &f, CELL_LOC outloc, DIFF_METHOD method) { + return interp_to(kY, outloc)*Grad2_par2(f, outloc, method) + Div_par(kY, outloc, method)*Grad_par(f, outloc, method); } /******************************************************************************* @@ -576,15 +536,15 @@ const Field3D Div_par_K_Grad_par(const Field3D &kY, const Field3D &f, const CELL * perpendicular Laplacian operator *******************************************************************************/ -const Field2D Delp2(const Field2D &f, const CELL_LOC outloc) { +const Field2D Delp2(const Field2D &f, CELL_LOC outloc) { return f.getCoordinates(outloc)->Delp2(f, outloc); } -const Field3D Delp2(const Field3D &f, BoutReal UNUSED(zsmooth), const CELL_LOC outloc) { +const Field3D Delp2(const Field3D &f, BoutReal UNUSED(zsmooth), CELL_LOC outloc) { return f.getCoordinates(outloc)->Delp2(f, outloc); } -const FieldPerp Delp2(const FieldPerp &f, BoutReal UNUSED(zsmooth), const CELL_LOC outloc) { +const FieldPerp Delp2(const FieldPerp &f, BoutReal UNUSED(zsmooth), CELL_LOC outloc) { return f.getCoordinates(outloc)->Delp2(f, outloc); } @@ -595,11 +555,11 @@ const FieldPerp Delp2(const FieldPerp &f, BoutReal UNUSED(zsmooth), const CELL_L * Laplace_perp = Laplace - Laplace_par *******************************************************************************/ -const Field2D Laplace_perp(const Field2D &f, const CELL_LOC outloc) { +const Field2D Laplace_perp(const Field2D &f, CELL_LOC outloc) { return Laplace(f, outloc) - Laplace_par(f, outloc); } -const Field3D Laplace_perp(const Field3D &f, const CELL_LOC outloc) { +const Field3D Laplace_perp(const Field3D &f, CELL_LOC outloc) { return Laplace(f, outloc) - Laplace_par(f, outloc); } @@ -611,11 +571,11 @@ const Field3D Laplace_perp(const Field3D &f, const CELL_LOC outloc) { * *******************************************************************************/ -const Field2D Laplace_par(const Field2D &f, const CELL_LOC outloc) { +const Field2D Laplace_par(const Field2D &f, CELL_LOC outloc) { return f.getCoordinates(outloc)->Laplace_par(f, outloc); } -const Field3D Laplace_par(const Field3D &f, const CELL_LOC outloc) { +const Field3D Laplace_par(const Field3D &f, CELL_LOC outloc) { return f.getCoordinates(outloc)->Laplace_par(f, outloc); } @@ -624,11 +584,11 @@ const Field3D Laplace_par(const Field3D &f, const CELL_LOC outloc) { * Full Laplacian operator on scalar field *******************************************************************************/ -const Field2D Laplace(const Field2D &f, const CELL_LOC outloc) { +const Field2D Laplace(const Field2D &f, CELL_LOC outloc) { return f.getCoordinates(outloc)->Laplace(f, outloc); } -const Field3D Laplace(const Field3D &f, const CELL_LOC outloc) { +const Field3D Laplace(const Field3D &f, CELL_LOC outloc) { return f.getCoordinates(outloc)->Laplace(f, outloc); } @@ -775,8 +735,7 @@ const Field3D b0xGrad_dot_Grad(const Field3D &phi, const Field3D &A, CELL_LOC ou result.name = "b0xGrad_dot_Grad("+phi.name+","+A.name+")"; #endif - ASSERT2(((outloc == CELL_DEFAULT) && (result.getLocation() == A.getLocation())) || - (result.getLocation() == outloc)); + ASSERT2(result.getLocation() == outloc); return result; } @@ -786,52 +745,26 @@ const Field3D b0xGrad_dot_Grad(const Field3D &phi, const Field3D &A, CELL_LOC ou * Terms of form b0 x Grad(f) dot Grad(g) / B = [f, g] *******************************************************************************/ -/*! - * Calculate location of result - */ -// use anonymous namespace so this function is only available in this file -namespace { - CELL_LOC bracket_location(const CELL_LOC &f_loc, const CELL_LOC &g_loc, const CELL_LOC &outloc, Mesh* localmesh=mesh) { - if(!localmesh->StaggerGrids) - return CELL_CENTRE; - - if(outloc == CELL_DEFAULT){ - // Check that f and g are in the same location - if (f_loc != g_loc){ - throw BoutException("Bracket currently requires both fields to have the same cell location"); - }else { - return f_loc; // Location of result - } - } - - // Check that f, and g are in the same location as the specified output location - if(f_loc != g_loc || f_loc != outloc){ - throw BoutException("Bracket currently requires the location of both fields and the output locaton to be the same"); - } - - return outloc; // Location of result - } -} - const Field2D bracket(const Field2D &f, const Field2D &g, BRACKET_METHOD method, CELL_LOC outloc, Solver *UNUSED(solver)) { TRACE("bracket(Field2D, Field2D)"); ASSERT1(f.getMesh() == g.getMesh()); + if (outloc == CELL_DEFAULT) { + outloc = g.getLocation(); + } + ASSERT1(f.getLocation() == g.getLocation() && outloc == f.getLocation()) Field2D result(f.getMesh()); - // Sort out cell locations - CELL_LOC result_loc = bracket_location(f.getLocation(), g.getLocation(), outloc, f.getMesh()); - if( (method == BRACKET_SIMPLE) || (method == BRACKET_ARAKAWA)) { // Use a subset of terms for comparison to BOUT-06 result = 0.0; + result.setLocation(outloc); }else { // Use full expression with all terms - result = b0xGrad_dot_Grad(f, g) / f.getCoordinates(result_loc)->Bxy; + result = b0xGrad_dot_Grad(f, g, outloc) / f.getCoordinates(outloc)->Bxy; } - result.setLocation(result_loc); return result; } @@ -840,14 +773,16 @@ const Field3D bracket(const Field3D &f, const Field2D &g, BRACKET_METHOD method, TRACE("bracket(Field3D, Field2D)"); ASSERT1(f.getMesh() == g.getMesh()); + if (outloc == CELL_DEFAULT) { + outloc = g.getLocation(); + } + ASSERT1(f.getLocation() == g.getLocation() && outloc == f.getLocation()) Mesh *mesh = f.getMesh(); Field3D result(mesh); - CELL_LOC result_loc = bracket_location(f.getLocation(), g.getLocation(), outloc, f.getMesh()); - - Coordinates *metric = f.getCoordinates(result_loc); + Coordinates *metric = f.getCoordinates(outloc); switch(method) { case BRACKET_CTU: { @@ -858,6 +793,7 @@ const Field3D bracket(const Field3D &f, const Field2D &g, BRACKET_METHOD method, throw BoutException("CTU method requires access to the solver"); result.allocate(); + result.setLocation(outloc); int ncz = mesh->LocalNz; for(int x=mesh->xstart;x<=mesh->xend;x++) @@ -895,6 +831,7 @@ const Field3D bracket(const Field3D &f, const Field2D &g, BRACKET_METHOD method, // Arakawa scheme for perpendicular flow. Here as a test result.allocate(); + result.setLocation(outloc); const BoutReal fac = 1.0 / (12 * metric->dz); const int ncz = mesh->LocalNz; @@ -965,6 +902,7 @@ const Field3D bracket(const Field3D &f, const Field2D &g, BRACKET_METHOD method, } case BRACKET_ARAKAWA_OLD: { result.allocate(); + result.setLocation(outloc); const int ncz = mesh->LocalNz; const BoutReal partialFactor = 1.0/(12 * metric->dz); BOUT_OMP(parallel for) @@ -1004,15 +942,14 @@ const Field3D bracket(const Field3D &f, const Field2D &g, BRACKET_METHOD method, } case BRACKET_SIMPLE: { // Use a subset of terms for comparison to BOUT-06 - result = VDDX(DDZ(f), g); + result = VDDX(DDZ(f, outloc), g, outloc); break; } default: { // Use full expression with all terms - result = b0xGrad_dot_Grad(f, g) / metric->Bxy; + result = b0xGrad_dot_Grad(f, g, outloc) / metric->Bxy; } } - result.setLocation(result_loc); return result; } @@ -1021,13 +958,15 @@ const Field3D bracket(const Field2D &f, const Field3D &g, BRACKET_METHOD method, TRACE("bracket(Field2D, Field3D)"); ASSERT1(f.getMesh() == g.getMesh()); + if (outloc == CELL_DEFAULT) { + outloc = g.getLocation(); + } + ASSERT1(f.getLocation() == g.getLocation() && outloc == f.getLocation()) Mesh *mesh = f.getMesh(); Field3D result(mesh); - CELL_LOC result_loc = bracket_location(f.getLocation(), g.getLocation(), outloc, f.getMesh()); - switch(method) { case BRACKET_CTU: throw BoutException("Bracket method CTU is not yet implemented for [2d,3d] fields."); @@ -1038,16 +977,15 @@ const Field3D bracket(const Field2D &f, const Field3D &g, BRACKET_METHOD method, break; case BRACKET_SIMPLE: { // Use a subset of terms for comparison to BOUT-06 - result = VDDZ(-DDX(f), g); + result = VDDZ(-DDX(f, outloc), g, outloc); break; } default: { // Use full expression with all terms - Coordinates *metric = f.getCoordinates(result_loc); - result = b0xGrad_dot_Grad(f, g) / metric->Bxy; + Coordinates *metric = f.getCoordinates(outloc); + result = b0xGrad_dot_Grad(f, g, outloc) / metric->Bxy; } } - result.setLocation(result_loc) ; return result; } @@ -1057,18 +995,20 @@ const Field3D bracket(const Field3D &f, const Field3D &g, BRACKET_METHOD method, TRACE("Field3D, Field3D"); ASSERT1(f.getMesh() == g.getMesh()); + if (outloc == CELL_DEFAULT) { + outloc = g.getLocation(); + } + ASSERT1(f.getLocation() == g.getLocation() && outloc == f.getLocation()) Mesh *mesh = f.getMesh(); Field3D result(mesh); - CELL_LOC result_loc = bracket_location(f.getLocation(), g.getLocation(), outloc, f.getMesh()); - - Coordinates *metric = f.getCoordinates(result_loc); + Coordinates *metric = f.getCoordinates(outloc); if (mesh->GlobalNx == 1 || mesh->GlobalNz == 1) { result=0; - result.setLocation(result_loc); + result.setLocation(outloc); return result; } @@ -1084,6 +1024,7 @@ const Field3D bracket(const Field3D &f, const Field3D &g, BRACKET_METHOD method, BoutReal dt = solver->getCurrentTimestep(); result.allocate(); + result.setLocation(outloc); FieldPerp vx(mesh), vz(mesh); vx.allocate(); @@ -1174,6 +1115,7 @@ const Field3D bracket(const Field3D &f, const Field3D &g, BRACKET_METHOD method, // Arakawa scheme for perpendicular flow result.allocate(); + result.setLocation(outloc); const int ncz = mesh->LocalNz; const BoutReal partialFactor = 1.0/(12 * metric->dz); @@ -1267,6 +1209,7 @@ const Field3D bracket(const Field3D &f, const Field3D &g, BRACKET_METHOD method, // Arakawa scheme for perpendicular flow result.allocate(); + result.setLocation(outloc); const int ncz = mesh->LocalNz; const BoutReal partialFactor = 1.0 / (12 * metric->dz); @@ -1316,16 +1259,14 @@ const Field3D bracket(const Field3D &f, const Field3D &g, BRACKET_METHOD method, } case BRACKET_SIMPLE: { // Use a subset of terms for comparison to BOUT-06 - result = VDDX(DDZ(f), g) + VDDZ(-DDX(f), g); + result = VDDX(DDZ(f, outloc), g, outloc) + VDDZ(-DDX(f, outloc), g, outloc); break; } default: { // Use full expression with all terms - result = b0xGrad_dot_Grad(f, g) / metric->Bxy; + result = b0xGrad_dot_Grad(f, g, outloc) / metric->Bxy; } } - result.setLocation(result_loc) ; - return result; } diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx index d872b70ab4..f684cce292 100644 --- a/src/mesh/index_derivs.cxx +++ b/src/mesh/index_derivs.cxx @@ -662,6 +662,7 @@ const Field2D Mesh::applyXdiff(const Field2D &var, Mesh::deriv_func func, Field2D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (this->StaggerGrids && (outloc != inloc)) { // Staggered differencing @@ -752,8 +753,6 @@ const Field2D Mesh::applyXdiff(const Field2D &var, Mesh::deriv_func func, } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -787,6 +786,7 @@ const Field3D Mesh::applyXdiff(const Field3D &var, Mesh::deriv_func func, Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (this->StaggerGrids && (outloc != inloc)) { // Staggered differencing @@ -877,8 +877,6 @@ const Field3D Mesh::applyXdiff(const Field3D &var, Mesh::deriv_func func, } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -912,6 +910,7 @@ const Field2D Mesh::applyYdiff(const Field2D &var, Mesh::deriv_func func, CELL_L Field2D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (this->ystart > 1) { // More than one guard cell, so set pp and mm values @@ -943,8 +942,6 @@ const Field2D Mesh::applyYdiff(const Field2D &var, Mesh::deriv_func func, CELL_L } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_yup = result.bndry_ydown = false; @@ -977,6 +974,7 @@ const Field3D Mesh::applyYdiff(const Field3D &var, Mesh::deriv_func func, CELL_L Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (var.hasYupYdown() && ((&var.yup() != &var) || (&var.ydown() != &var))) { // Field "var" has distinct yup and ydown fields which @@ -1122,8 +1120,6 @@ const Field3D Mesh::applyYdiff(const Field3D &var, Mesh::deriv_func func, CELL_L result = this->fromFieldAligned(result); } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -1156,6 +1152,7 @@ const Field3D Mesh::applyZdiff(const Field3D &var, Mesh::deriv_func func, CELL_L Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); // Check that the input variable has data ASSERT1(var.isAllocated()); @@ -1174,8 +1171,6 @@ const Field3D Mesh::applyZdiff(const Field3D &var, Mesh::deriv_func func, CELL_L } } - result.setLocation(outloc); - return result; } @@ -1197,8 +1192,6 @@ const Field3D Mesh::indexDDX(const Field3D &f, CELL_LOC outloc, DIFF_METHOD meth ASSERT1(outloc == inloc || (outloc == CELL_CENTRE && inloc == CELL_XLOW) || (outloc == CELL_XLOW && inloc == CELL_CENTRE)); - Field3D result(this); - if (this->StaggerGrids && (outloc != inloc)) { // Shifting in X. Centre -> Xlow, or Xlow -> Centre @@ -1213,9 +1206,7 @@ const Field3D Mesh::indexDDX(const Field3D &f, CELL_LOC outloc, DIFF_METHOD meth throw BoutException("Cannot use FFT for X derivatives"); } - result = applyXdiff(f, func, outloc, region); - - return result; + return applyXdiff(f, func, outloc, region); } const Field2D Mesh::indexDDX(const Field2D &f, CELL_LOC outloc, @@ -1239,8 +1230,6 @@ const Field3D Mesh::indexDDY(const Field3D &f, CELL_LOC outloc, DIFF_METHOD meth ASSERT1(outloc == inloc || (outloc == CELL_CENTRE && inloc == CELL_YLOW) || (outloc == CELL_YLOW && inloc == CELL_CENTRE)); - Field3D result(this); - if (this->StaggerGrids && (outloc != inloc)) { // Shifting in Y. Centre -> Ylow, or Ylow -> Centre func = sfDDY; // Set default @@ -1254,9 +1243,7 @@ const Field3D Mesh::indexDDY(const Field3D &f, CELL_LOC outloc, DIFF_METHOD meth throw BoutException("Cannot use FFT for Y derivatives"); } - result = applyYdiff(f, func, outloc, region); - - return result; + return applyYdiff(f, func, outloc, region); } const Field2D Mesh::indexDDY(const Field2D &f, CELL_LOC outloc, @@ -1420,8 +1407,6 @@ const Field3D Mesh::indexD2DX2(const Field3D &f, CELL_LOC outloc, ASSERT1(this == f.getMesh()); - Field3D result(this); - if (StaggerGrids && (outloc != inloc)) { // Shifting in X. Centre -> Xlow, or Xlow -> Centre func = sfD2DX2; // Set default @@ -1435,9 +1420,7 @@ const Field3D Mesh::indexD2DX2(const Field3D &f, CELL_LOC outloc, throw BoutException("Cannot use FFT for X derivatives"); } - result = applyXdiff(f, func, outloc, region); - - return result; + return applyXdiff(f, func, outloc, region); } /*! @@ -1483,8 +1466,6 @@ const Field3D Mesh::indexD2DY2(const Field3D &f, CELL_LOC outloc, ASSERT1(outloc == inloc || (outloc == CELL_CENTRE && inloc == CELL_YLOW) || (outloc == CELL_YLOW && inloc == CELL_CENTRE)); - Field3D result(this); - if (StaggerGrids && (outloc != inloc)) { // Shifting in Y. Centre -> Ylow, or Ylow -> Centre func = sfD2DY2; // Set default @@ -1498,9 +1479,7 @@ const Field3D Mesh::indexD2DY2(const Field3D &f, CELL_LOC outloc, throw BoutException("Cannot use FFT for Y derivatives"); } - result = applyYdiff(f, func, outloc, region); - - return result; + return applyYdiff(f, func, outloc, region); } /*! @@ -1726,6 +1705,7 @@ const Field2D Mesh::indexVDDX(const Field2D &v, const Field2D &f, CELL_LOC outlo Field2D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (this->xstart > 1) { // Two or more guard cells @@ -1761,8 +1741,6 @@ const Field2D Mesh::indexVDDX(const Field2D &v, const Field2D &f, CELL_LOC outlo result.bndry_xin = result.bndry_xout = false; #endif - result.setLocation(outloc); - return result; } @@ -1776,6 +1754,7 @@ const Field3D Mesh::indexVDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo ASSERT1(this == v.getMesh()); ASSERT1(this == f.getMesh()); + CELL_LOC vloc = v.getLocation(); CELL_LOC inloc = f.getLocation(); // Input location if (outloc == CELL_DEFAULT) @@ -1787,6 +1766,7 @@ const Field3D Mesh::indexVDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); /// Convert REGION enum to a Region string identifier const auto region_str = REGION_STRING(region); @@ -1926,8 +1906,6 @@ const Field3D Mesh::indexVDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -1957,10 +1935,10 @@ const Field2D Mesh::indexVDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo Field2D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (this->LocalNy == 1){ result=0; - result.setLocation(outloc); return result; } @@ -2102,8 +2080,6 @@ const Field2D Mesh::indexVDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -2131,10 +2107,10 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); if (this->LocalNy == 1){ result=0; - result.setLocation(outloc); return result; } @@ -2158,12 +2134,8 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo func = lookupFunc(table, method); } - // There are four cases, corresponding to whether or not f and v - // have yup, ydown fields. - - // If vUseUpDown is true, field "v" has distinct yup and ydown fields which - // will be used to calculate a derivative along - // the magnetic field + // If *UseUpDown is true, field "*" has distinct yup and ydown fields which + // will be used to calculate a derivative along the magnetic field bool vUseUpDown = (v.hasYupYdown() && ((&v.yup() != &v) || (&v.ydown() != &v))); bool fUseUpDown = (f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f))); @@ -2195,6 +2167,9 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo } } else { // Both must shift to field aligned + // (even if one of v and f has yup/ydown fields, it doesn't make sense to + // multiply them with one in field-aligned and one in non-field-aligned + // coordinates) Field3D v_fa = this->toFieldAligned(v); Field3D f_fa = this->toFieldAligned(f); BOUT_OMP(parallel) { @@ -2208,7 +2183,7 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo fval.mm = f_fa[i.ymm()]; fval.m = f_fa[i.ym()]; - fval.c = f[i]; + fval.c = f_fa[i]; fval.p = f_fa[i.yp()]; fval.pp = f_fa[i.ypp()]; @@ -2225,6 +2200,8 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo result[i] = func(vval, fval); } } + + result = this->fromFieldAligned(result); } } else { // Non-staggered case @@ -2286,8 +2263,6 @@ const Field3D Mesh::indexVDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -2317,6 +2292,7 @@ const Field3D Mesh::indexVDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); /// Convert REGION enum to a Region string identifier const auto region_str = REGION_STRING(region); @@ -2387,8 +2363,6 @@ const Field3D Mesh::indexVDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -2426,6 +2400,7 @@ const Field2D Mesh::indexFDDX(const Field2D &v, const Field2D &f, CELL_LOC outlo Field2D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); ASSERT1(this == v.getMesh()); ASSERT1(this == f.getMesh()); @@ -2520,6 +2495,7 @@ const Field3D Mesh::indexFDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); /// Convert REGION enum to a Region string identifier const auto region_str = REGION_STRING(region); @@ -2651,8 +2627,6 @@ const Field3D Mesh::indexFDDX(const Field3D &v, const Field3D &f, CELL_LOC outlo } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = result.bndry_yup = result.bndry_ydown = false; @@ -2690,7 +2664,7 @@ const Field2D Mesh::indexFDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo Field2D result(this); result.allocate(); // Make sure data allocated - result.setLocation(f.getLocation()); + result.setLocation(outloc); /// Convert REGION enum to a Region string identifier const auto region_str = REGION_STRING(region); @@ -2734,8 +2708,6 @@ const Field2D Mesh::indexFDDY(const Field2D &v, const Field2D &f, CELL_LOC outlo } } - result.setLocation(outloc); - #if CHECK > 0 // Mark boundaries as invalid result.bndry_xin = result.bndry_xout = false; @@ -2787,13 +2759,10 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); - // There are four cases, corresponding to whether or not f and v - // have yup, ydown fields. - - // If vUseUpDown is true, field "v" has distinct yup and ydown fields which - // will be used to calculate a derivative along - // the magnetic field + // If *UseUpDown is true, field "*" has distinct yup and ydown fields which + // will be used to calculate a derivative along the magnetic field bool vUseUpDown = (v.hasYupYdown() && ((&v.yup() != &v) || (&v.ydown() != &v))); bool fUseUpDown = (f.hasYupYdown() && ((&f.yup() != &f) || (&f.ydown() != &f))); @@ -2830,6 +2799,9 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo } } else { // Both must shift to field aligned + // (even if one of v and f has yup/ydown fields, it doesn't make sense to + // multiply them with one in field-aligned and one in non-field-aligned + // coordinates) Field3D v_fa = this->toFieldAligned(v); Field3D f_fa = this->toFieldAligned(f); BOUT_OMP(parallel) { @@ -2862,9 +2834,9 @@ const Field3D Mesh::indexFDDY(const Field3D &v, const Field3D &f, CELL_LOC outlo result[i] = func(vval, fval); } } - } - result.setLocation(outloc); + result = this->fromFieldAligned(result); + } #if CHECK > 0 // Mark boundaries as invalid @@ -2916,6 +2888,7 @@ const Field3D Mesh::indexFDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo Field3D result(this); result.allocate(); // Make sure data allocated + result.setLocation(outloc); /// Convert REGION enum to a Region string identifier const auto region_str = REGION_STRING(region); @@ -2951,7 +2924,6 @@ const Field3D Mesh::indexFDDZ(const Field3D &v, const Field3D &f, CELL_LOC outlo result[i] = func(vval, fval); } } - result.setLocation(outloc); #if CHECK > 0 // Mark boundaries as invalid diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx index 34b1c7709f..65494bd908 100644 --- a/src/mesh/mesh.cxx +++ b/src/mesh/mesh.cxx @@ -77,8 +77,14 @@ int Mesh::get(Field2D &var, const string &name, BoutReal def) { if (source == nullptr or !source->get(this, var, name, def)) return 1; - // Communicate to get guard cell data - Mesh::communicate(var); + if(!source->hasYGuards()) { + // If the source does not include y-boundary guard cells then communications + // are needed to fill them. This provides backwards compatibility even though + // the communicated guard cell values may not be correct, because the + // integrated shear is not necessarily continuous when poloidal angle goes + // from 2pi->0 + Mesh::communicate(var); + } // Check that the data is valid checkData(var); @@ -86,7 +92,7 @@ int Mesh::get(Field2D &var, const string &name, BoutReal def) { return 0; } -int Mesh::get(Field3D &var, const string &name, BoutReal def, bool communicate) { +int Mesh::get(Field3D &var, const string &name, BoutReal def, bool allow_communicate) { TRACE("Loading 3D field: Mesh::get(Field3D, %s)", name.c_str()); // Ensure data allocated @@ -95,8 +101,12 @@ int Mesh::get(Field3D &var, const string &name, BoutReal def, bool communicate) if (source == nullptr or !source->get(this, var, name, def)) return 1; - // Communicate to get guard cell data - if(communicate) { + if(!source->hasYGuards() && allow_communicate) { + // If the source does not include y-boundary guard cells then communications + // are needed to fill them. This provides backwards compatibility even though + // the communicated guard cell values may not be correct, because the + // integrated shear is not necessarily continuous when poloidal angle goes + // from 2pi->0 Mesh::communicate(var); } @@ -335,24 +345,24 @@ std::shared_ptr Mesh::createDefaultCoordinates(const CELL_LOC locat } -Region<> & Mesh::getRegion3D(const std::string ®ion_name){ - auto found = regionMap3D.find(region_name); - if (found == end(regionMap3D)) { - throw BoutException("Couldn't find region %s in regionMap3D", region_name.c_str()); - } - return found->second; +const Region<> & Mesh::getRegion3D(const std::string ®ion_name) const { + const auto found = regionMap3D.find(region_name); + if (found == end(regionMap3D)) { + throw BoutException("Couldn't find region %s in regionMap3D", region_name.c_str()); + } + return found->second; } -Region & Mesh::getRegion2D(const std::string ®ion_name){ - auto found = regionMap2D.find(region_name); - if (found == end(regionMap2D)) { - throw BoutException("Couldn't find region %s in regionMap2D", region_name.c_str()); - } - return found->second; +const Region & Mesh::getRegion2D(const std::string ®ion_name) const { + const auto found = regionMap2D.find(region_name); + if (found == end(regionMap2D)) { + throw BoutException("Couldn't find region %s in regionMap2D", region_name.c_str()); + } + return found->second; } -Region &Mesh::getRegionPerp(const std::string ®ion_name) { - auto found = regionMapPerp.find(region_name); +const Region &Mesh::getRegionPerp(const std::string ®ion_name) const { + const auto found = regionMapPerp.find(region_name); if (found == end(regionMapPerp)) { throw BoutException("Couldn't find region %s in regionMapPerp", region_name.c_str()); } diff --git a/tests/MMS/derivatives3/data/BOUT.inp b/tests/MMS/derivatives3/data/BOUT.inp index 13db66011d..e203480ccf 100644 --- a/tests/MMS/derivatives3/data/BOUT.inp +++ b/tests/MMS/derivatives3/data/BOUT.inp @@ -9,20 +9,17 @@ ny=32 nz=n mxg=0 myg=2 -dy=1/ny [solver] rtol = 1e-14 atol = 1e-18 - [meshz] staggergrids=true nx=1 ny=1 nz=n n=2 -dz=2*Pi/nz MXG=0 MYG=0 @@ -32,18 +29,15 @@ nx=1 ny=n nz=1 n=6 -dy=2*Pi/ny MXG=0 MYG=2 - [meshx] staggergrids=true nx=n+2*mxg ny=1 nz=1 n=2 -dx=2*Pi/ny MXG=2 MYG=0 diff --git a/tests/MMS/derivatives3/runtest b/tests/MMS/derivatives3/runtest index c3f64624cb..21d14349bd 100755 --- a/tests/MMS/derivatives3/runtest +++ b/tests/MMS/derivatives3/runtest @@ -30,12 +30,23 @@ def runtests(functions,derivatives,directions,stag,msg): boutcore.setOption("meshD:nD".replace("D",direction) ,"%d"% (nz+ (2*guards if direction == "x" else 0)),force=True) boutcore.setOption("meshD:dD".replace("D",direction,) - ,"2*pi/(%d)"%(nz),force=True) + ,"2*pi/(%d)"%(nz) if direction == "z" else "2*pi*(1+.1*sin(2*pi*x)+.1*sin(y))/(%d)"%(nz),force=True) + boutcore.setOption("meshD:ddD:first".replace("D",direction,),diff,force=True) + boutcore.setOption("meshD:ddD:second".replace("D",direction,),diff,force=True) + if stag: + boutcore.setOption("meshD:ddD:firststag".replace("D",direction,),diff,force=True) + # not all methods available for secondstag, we only test "C2" + # so can skip setting here (get error otherwise) + #boutcore.setOption("meshD:ddD:secondstag".replace("D",direction,),diff,force=True) + else: + # method 'diff' may not be available for staggered derivatives, so set C2 explicitly + boutcore.setOption("meshD:ddD:firststag".replace("D",direction,),"C2",force=True) + boutcore.setOption("meshD:ddD:secondstag".replace("D",direction,),"C2",force=True) dirnfac=direction+"*"+fac mesh=boutcore.Mesh(section="mesh"+direction) f=boutcore.create3D(infunc.replace("%s",dirnfac),mesh ,outloc=inloc) - sim=diff_func(f,method=diff,outloc=outloc) + sim=diff_func(f,outloc=outloc) if sim.getLocation() != outloc: cent=['CENTRE','CENTER'] if outloc in cent and sim.getLocation() in cent: @@ -71,8 +82,8 @@ doPlot=False nzs=np.logspace(start,mmax,num=mmax-start+1,base=2) functions=[ - ["sin(%s)","cos(%s)"] , - ["cos(%s)", "-sin(%s)"] + ["sin(%s)","cos(%s)/(1+.1*sin(2*pi*x)+.1*sin(y))"] , + ["cos(%s)", "-sin(%s)/(1+.1*sin(2*pi*x)+.1*sin(y))"] ] derivatives=[ @@ -99,9 +110,10 @@ derivatives=[ runtests(functions,derivatives,directions,stag=True,msg="DD") +# NB need to include derivative of grid-spacing here functions=[ - ["sin(%s)","-sin(%s)"], - ["cos(%s)" , "-cos(%s)"] + ["cos(%s)" , "-cos(%s)/(1+.1*sin(2*pi*x)+.1*sin(y))^2 + .1*sin(%s)*cos(%s)/(1+.1*sin(2*pi*x)+.1*sin(y))^3"], + ["sin(%s)","-sin(%s)/(1+.1*sin(2*pi*x)+.1*sin(y))^2 - .1*(cos(%s)^2)/(1+.1*sin(2*pi*x)+.1*sin(y))^3"], ] derivatives=[ diff --git a/tests/MMS/difops/data/BOUT.inp b/tests/MMS/difops/data/BOUT.inp new file mode 100644 index 0000000000..2a4d87d691 --- /dev/null +++ b/tests/MMS/difops/data/BOUT.inp @@ -0,0 +1,16 @@ +# partial input file, further options are set using boutcore +# interface in runtest + +[mesh] +nx = 6 +ny = 2 +nz = 1 + +[testmesh:ddx] +upwind = U2 + +[testmesh:ddy] +upwind = U2 + +[testmesh:ddz] +upwind = U2 diff --git a/tests/MMS/difops/runtest b/tests/MMS/difops/runtest new file mode 100755 index 0000000000..7cfe03da38 --- /dev/null +++ b/tests/MMS/difops/runtest @@ -0,0 +1,9 @@ +#!/bin/bash + +#requires boutcore +#requires all_tests +#requires not make + +./runtest.py + +exit # exit with status of last command diff --git a/tests/MMS/difops/runtest.py b/tests/MMS/difops/runtest.py new file mode 100755 index 0000000000..5e58c244f0 --- /dev/null +++ b/tests/MMS/difops/runtest.py @@ -0,0 +1,538 @@ +#!/usr/bin/env python3 + +# MMS test for differential operators (that use the metric) + +import boutcore +from boutdata.mms_alternate import * + +import numpy +import sympy +from copy import copy +from sys import exit + +# get command line arguments +import argparse +parser = argparse.ArgumentParser() +parser.add_argument('--short', action='store_true', default=False) +args = parser.parse_args() +full_test = not args.short + +# geometry for simple circular tokamak +tokamak = SimpleTokamak() + +# rescale x and y coordinates so dx and dy are not constants +tokamak.set_scalex(1 + .1*sin(2*pi*metric.x+metric.y)) +tokamak.set_scaley(1 + .1*sin(2*pi*metric.x-metric.y)) +# re-calculate metric terms +tokamak.metric() + +def test_operator(ngrids, testfunc, dimensions, boutcore_operator, symbolic_operator, order, ftype, method, stagger, mesh_in=None): + + testfunc = copy(testfunc) # ensure we don't change global testfunc + error_list = [] + if full_test: + print('testing',boutcore_operator, ftype, stagger) + for n in ngrids: + if full_test: + print('n =',n) + if mesh_in is not None: + mesh = mesh_in + else: + # set options + boutcore.setOption('mxg', str(mxg), force=True) + boutcore.setOption('myg', str(myg), force=True) + # set up mesh input + if 'x' in dimensions: + boutcore.setOption('testmesh:nx', exprToStr(n+2*mxg), force=True) + else: + boutcore.setOption('testmesh:nx', exprToStr(default_n+2*mxg), force=True) + if 'y' in dimensions: + boutcore.setOption('testmesh:ny', exprToStr(n), force=True) + else: + boutcore.setOption('testmesh:ny', exprToStr(default_n), force=True) + if 'z' in dimensions: + boutcore.setOption('testmesh:nz', exprToStr(n), force=True) + else: + boutcore.setOption('testmesh:nz', exprToStr(default_n), force=True) + boutcore.setOption('testmesh:dx', exprToStr(metric.psiwidth*metric.scalex/n), force=True) + boutcore.setOption('testmesh:dy', exprToStr(2.*pi*metric.scaley/n), force=True) + boutcore.setOption('testmesh:dz', exprToStr(2.*pi/n), force=True) + boutcore.setOption('testmesh:g11', exprToStr(metric.g11), force=True) + boutcore.setOption('testmesh:g22', exprToStr(metric.g22), force=True) + boutcore.setOption('testmesh:g33', exprToStr(metric.g33), force=True) + boutcore.setOption('testmesh:g12', exprToStr(metric.g12), force=True) + boutcore.setOption('testmesh:g13', exprToStr(metric.g13), force=True) + boutcore.setOption('testmesh:g23', exprToStr(metric.g23), force=True) + boutcore.setOption('testmesh:g_11', exprToStr(metric.g_11), force=True) + boutcore.setOption('testmesh:g_22', exprToStr(metric.g_22), force=True) + boutcore.setOption('testmesh:g_33', exprToStr(metric.g_33), force=True) + boutcore.setOption('testmesh:g_12', exprToStr(metric.g_12), force=True) + boutcore.setOption('testmesh:g_13', exprToStr(metric.g_13), force=True) + boutcore.setOption('testmesh:g_23', exprToStr(metric.g_23), force=True) + if stagger is None: + boutcore.setOption('testmesh:staggergrids', str('false'), force=True) + else: + boutcore.setOption('testmesh:staggergrids', str('true'), force=True) + + # create new Mesh object + mesh = boutcore.Mesh(section='testmesh') + + if stagger is None: + inloc = 'CENTRE' + outloc = 'CENTRE' + else: + inloc = stagger[0] + outloc = stagger[1] + + # calculate result of differential operator using BOUT++ implementation + if ftype == '2D': + # cannot have z-dependence + testfunc = testfunc.replace('z', '0') + bout_input = boutcore.create2D(testfunc, mesh, outloc=inloc) + elif ftype == '3D': + bout_input = boutcore.create3D(testfunc, mesh, outloc=inloc) + else: + raise ValueError('Unexpected ftype argument '+str(ftype)) + if method is None: + bout_result = boutcore_operator(bout_input, outloc=outloc) + else: + bout_result = boutcore_operator(bout_input, outloc=outloc, method=method) + + # calculate result of differential operator symbolically, then convert to boutcore.Field3D/Field2D + analytic_input = sympy.sympify(testfunc) + analytic_func = symbolic_operator(analytic_input) + if ftype == '2D': + analytic_result = boutcore.create2D(exprToStr(analytic_func), mesh, outloc=outloc) + elif ftype == '3D': + analytic_result = boutcore.create3D(exprToStr(analytic_func), mesh, outloc=outloc) + else: + raise ValueError('Unexpected ftype argument '+str(ftype)) + + # calculate max error + error = bout_result - analytic_result # as Field3D/Field2D + error = error.get()[mxg:-mxg, myg:-myg] # numpy array, without guard cells + error_list.append(numpy.max(numpy.abs(error))) # max error + + logerrors = numpy.log(error_list[-2]/error_list[-1]) + logspacing = numpy.log(ngrids[-1]/ngrids[-2]) + convergence = logerrors/logspacing + + if order-.1 < convergence < order+.2: + return ['pass'] + else: + if plot_error: + from matplotlib import pyplot + pyplot.loglog(1./ngrids, error_list) + pyplot.show() + from boututils.showdata import showdata + showdata(error) + return [str(boutcore_operator)+' is not working for '+inloc+'->'+outloc+' '+str(ftype)+' '+str(method)+'. Expected '+str(order)+', got '+str(convergence)+'.'] + +def cycle_staggering(stagger_directions, base_dimensions, ngrids, testfunc, boutcore_operator, symbolic_operator, order, types, method=None): + """ + Loop over different parameters, calling test_operator for each + + Parameters + ---------- + + stagger_directions : str + Directions in which this operator can be staggered. String containing + 'x', 'y' or 'z' will test both permutations of staggering between + CENTRE and the corresponding direction. + base_dimensions : str + String containing any of 'x', 'y' and 'z': the grid must be refined in + the directions contained in this argument in order to converge. E.g. + Grad_par requires refinement only in the y-direction. + ngrids : numpy.array(int) + Array of grid sizes to use. All directions being refined are given the + same size. + testfunc : sympy expression + The input that will be given to the function being tested + boutcore_operator : function + function from boutcore to be tested + symbolic_operator : function + function using sympy to do the symbolic equivalent of boutcore_operator + order : int + expected order of convergence of boutcore_operator + types : list of str + list of strings giving the type of the argument to the operators + respectively. '2D' for Field2D and '3D' for Field3D. + """ + + # all derivatives at same inloc/outloc should work + dimensions_staggers = [(base_dimensions, None), # no staggering + (base_dimensions+'x', ('CENTRE', 'CENTRE')), # include all-centred, but with StaggerGrids=true + (base_dimensions+'x', ('XLOW', 'XLOW')), + (base_dimensions+'y', ('YLOW', 'YLOW')), + (base_dimensions+'z', ('ZLOW', 'ZLOW'))] + if 'x' in stagger_directions: + dimensions_staggers += [(base_dimensions+'x', ('CENTRE', 'XLOW')), + (base_dimensions+'x', ('XLOW', 'CENTRE'))] + if 'y' in stagger_directions: + dimensions_staggers += [(base_dimensions+'y', ('CENTRE', 'YLOW')), + (base_dimensions+'y', ('YLOW', 'CENTRE'))] + if 'z' in stagger_directions: + dimensions_staggers += [(base_dimensions+'z', ('CENTRE', 'ZLOW')), + (base_dimensions+'z', ('ZLOW', 'CENTRE'))] + result = [] + for ftype in types: + for dimensions, stagger in dimensions_staggers: + result += test_operator(ngrids, testfunc, dimensions, boutcore_operator, symbolic_operator, order, ftype, method, stagger) + + if test_throw: + # check that unsupported combinations of locations throw an exception + + locations = ['CENTRE', 'XLOW', 'YLOW', 'ZLOW'] + # first make a list of all permutations + fail_staggers = [(x,y) for x in locations for y in locations] + # now remove the ones have already tested + for dimensions, stagger in dimensions_staggers: + if stagger is not None: + index = fail_staggers.index(stagger) + del fail_staggers[index] + boutcore.setOption('failmesh:nx', '8', force=True) + boutcore.setOption('failmesh:ny', '4', force=True) + boutcore.setOption('failmesh:nz', '4', force=True) + boutcore.setOption('failmesh:staggergrids', 'true', force=True) + failmesh = boutcore.Mesh(section='failmesh') + for stagger in fail_staggers: + # check that an exception is throw for combinations of directions that we expect to fail + try: + test = test_operator(numpy.array([4, 8]), testfunc, dimensions, boutcore_operator, symbolic_operator, order, ftype, method, stagger, mesh_in=failmesh) + except RuntimeError: + result += ['pass'] + else: + result += ['Expected '+str(boutcore_operator)+' to throw for '+stagger[0]+'->'+stagger[1]+' '+' '+str(ftype)+' '+str(method)+' but it did not.'] + + return result + +def test_operator2(ngrids, testfunc1, testfunc2, dimensions, boutcore_operator, symbolic_operator, order, ftypes, method, stagger, mesh_in=None): + + testfunc1 = copy(testfunc1) # ensure we don't change global testfunc1 + testfunc2 = copy(testfunc2) # ensure we don't change global testfunc2 + error_list = [] + if full_test: + print('testing', boutcore_operator, ftypes, stagger) + for n in ngrids: + if full_test: + print('n =',n) + if mesh_in is not None: + mesh = mesh_in + else: + # set options + boutcore.setOption('mxg', str(mxg), force=True) + boutcore.setOption('myg', str(myg), force=True) + # set up mesh input + if 'x' in dimensions: + boutcore.setOption('testmesh:nx', exprToStr(n+2*mxg), force=True) + else: + boutcore.setOption('testmesh:nx', exprToStr(default_n+2*mxg), force=True) + if 'y' in dimensions: + boutcore.setOption('testmesh:ny', exprToStr(n), force=True) + else: + boutcore.setOption('testmesh:ny', exprToStr(default_n), force=True) + if 'z' in dimensions: + boutcore.setOption('testmesh:nz', exprToStr(n), force=True) + else: + boutcore.setOption('testmesh:nz', exprToStr(default_n), force=True) + boutcore.setOption('testmesh:dx', exprToStr(metric.psiwidth*metric.scalex/n), force=True) + boutcore.setOption('testmesh:dy', exprToStr(2.*pi*metric.scaley/n), force=True) + boutcore.setOption('testmesh:dz', exprToStr(2.*pi/n), force=True) + boutcore.setOption('testmesh:g11', exprToStr(metric.g11), force=True) + boutcore.setOption('testmesh:g22', exprToStr(metric.g22), force=True) + boutcore.setOption('testmesh:g33', exprToStr(metric.g33), force=True) + boutcore.setOption('testmesh:g12', exprToStr(metric.g12), force=True) + boutcore.setOption('testmesh:g13', exprToStr(metric.g13), force=True) + boutcore.setOption('testmesh:g23', exprToStr(metric.g23), force=True) + boutcore.setOption('testmesh:g_11', exprToStr(metric.g_11), force=True) + boutcore.setOption('testmesh:g_22', exprToStr(metric.g_22), force=True) + boutcore.setOption('testmesh:g_33', exprToStr(metric.g_33), force=True) + boutcore.setOption('testmesh:g_12', exprToStr(metric.g_12), force=True) + boutcore.setOption('testmesh:g_13', exprToStr(metric.g_13), force=True) + boutcore.setOption('testmesh:g_23', exprToStr(metric.g_23), force=True) + if stagger is None: + boutcore.setOption('testmesh:staggergrids', str('false'), force=True) + else: + boutcore.setOption('testmesh:staggergrids', str('true'), force=True) + + # create new Mesh object + mesh = boutcore.Mesh(section='testmesh') + + if stagger is None: + vloc = 'CENTRE' + inloc = 'CENTRE' + outloc = 'CENTRE' + else: + vloc = stagger[0] + inloc = stagger[1] + outloc = stagger[2] + + # calculate result of differential operator using BOUT++ implementation + if ftypes[0] == '2D': + # cannot have z-dependence + testfunc1 = testfunc1.replace('z', '0') + bout_input1 = boutcore.create2D(testfunc1, mesh, outloc=vloc) + elif ftypes[0] == '3D': + bout_input1 = boutcore.create3D(testfunc1, mesh, outloc=vloc) + else: + raise ValueError('Unexpected ftype argument '+str(ftypes[0])) + if ftypes[1] == '2D': + # cannot have z-dependence + testfunc2 = testfunc2.replace('z', '0') + bout_input2 = boutcore.create2D(testfunc2, mesh, outloc=inloc) + elif ftypes[1] == '3D': + bout_input2 = boutcore.create3D(testfunc2, mesh, outloc=inloc) + else: + raise ValueError('Unexpected ftype argument '+str(ftypes[1])) + if method is None: + bout_result = boutcore_operator(bout_input1, bout_input2, outloc=outloc) + else: + bout_result = boutcore_operator(bout_input1, bout_input2, outloc=outloc, method=method) + + # calculate result of differential operator symbolically, then convert to boutcore.Field3D/Field2D + analytic_input1 = sympy.sympify(testfunc1) + analytic_input2 = sympy.sympify(testfunc2) + analytic_func = symbolic_operator(analytic_input1, analytic_input2) + if ftypes[0] == '2D' and ftypes[1] == '2D': + analytic_result = boutcore.create2D(exprToStr(analytic_func), mesh, outloc=outloc) + else: + analytic_result = boutcore.create3D(exprToStr(analytic_func), mesh, outloc=outloc) + + # calculate max error + error = bout_result - analytic_result # as Field3D + error = error.get()[mxg:-mxg, myg:-myg] # numpy array, without guard cells + error_list.append(numpy.max(numpy.abs(error))) # max error + + logerrors = numpy.log(error_list[-2]/error_list[-1]) + logspacing = numpy.log(ngrids[-1]/ngrids[-2]) + convergence = logerrors/logspacing + + if order-.1 < convergence < order+.2: + return ['pass'] + else: + if plot_error: + from matplotlib import pyplot + pyplot.loglog(1./ngrids, error_list) + pyplot.show() + from boututils.showdata import showdata + showdata([error, bout_result.get()[mxg:-mxg, myg:-myg], analytic_result.get()[mxg:-mxg, myg:-myg]]) + return [str(boutcore_operator)+' is not working for '+inloc+'->'+outloc+' '+str(ftypes)+' '+str(method)+'. Expected '+str(order)+', got '+str(convergence)+'.'] + +def cycle_staggering2(stagger_directions, base_dimensions, ngrids, testfunc1, testfunc2, boutcore_operator, symbolic_operator, order, types, method=None): + """ + Loop over different parameters, calling test_operator2 for each + + Parameters + ---------- + + stagger_directions : str + Directions in which this operator can be staggered. String containing + 'xx', 'yy' or 'zz' will test all permutations of staggering between + CENTRE and the corresponding direction. String containing 'x', 'y' or + 'z' will keep second argument and outloc at the same location, and + stagger the first argument in the corresponding direction. + base_dimensions : str + String containing any of 'x', 'y' and 'z': the grid must be refined in + the directions contained in this argument in order to converge. E.g. + Grad_par requires refinement only in the y-direction. + ngrids : numpy.array(int) + Array of grid sizes to use. All directions being refined are given the + same size. + testfunc1, testfunc2 : sympy expression + The inputs that will be given to the function being tested + boutcore_operator : function + function from boutcore to be tested + symbolic_operator : function + function using sympy to do the symbolic equivalent of boutcore_operator + order : int + expected order of convergence of boutcore_operator + types : list of tuples of str + list of pairs (type1, type2) giving the type of first and second + arguments to the operators respectively. '2D' for Field2D and '3D' for + Field3D. + method : str + DIFF_METHOD to pass to boutcore_operator. If method=None then method + argument will not be passed, so the default will be used. + """ + + # all derivatives at same inloc/outloc should work + dimensions_staggers = [(base_dimensions, None), # no staggering + (base_dimensions+'x', ('CENTRE', 'CENTRE', 'CENTRE')), # include all-centred, but with StaggerGrids=true + (base_dimensions+'x', ('XLOW', 'XLOW', 'XLOW')), + (base_dimensions+'y', ('YLOW', 'YLOW', 'YLOW')), + (base_dimensions+'z', ('ZLOW', 'ZLOW', 'ZLOW'))] + if 'xx' in stagger_directions: + # for xx include all permutations of XLOW and CENTRE + dimensions_staggers += [(base_dimensions+'x', ('CENTRE', 'CENTRE', 'XLOW')), + (base_dimensions+'x', ('CENTRE', 'XLOW', 'CENTRE')), + (base_dimensions+'x', ('XLOW', 'CENTRE', 'CENTRE')), + (base_dimensions+'x', ('XLOW', 'XLOW', 'CENTRE')), + (base_dimensions+'x', ('XLOW', 'CENTRE', 'XLOW')), + (base_dimensions+'x', ('CENTRE', 'XLOW', 'XLOW'))] + elif 'x' in stagger_directions: + dimensions_staggers += [(base_dimensions+'x', ('CENTRE', 'XLOW', 'XLOW')), + (base_dimensions+'x', ('XLOW', 'CENTRE', 'CENTRE'))] + if 'yy' in stagger_directions: + # for yy include all permutations of YLOW and CENTRE + dimensions_staggers += [(base_dimensions+'y', ('CENTRE', 'CENTRE', 'YLOW')), + (base_dimensions+'y', ('CENTRE', 'YLOW', 'CENTRE')), + (base_dimensions+'y', ('YLOW', 'CENTRE', 'CENTRE')), + (base_dimensions+'y', ('YLOW', 'YLOW', 'CENTRE')), + (base_dimensions+'y', ('YLOW', 'CENTRE', 'YLOW')), + (base_dimensions+'y', ('CENTRE', 'YLOW', 'YLOW'))] + elif 'y' in stagger_directions: + dimensions_staggers += [(base_dimensions+'y', ('CENTRE', 'YLOW', 'YLOW')), + (base_dimensions+'y', ('YLOW', 'CENTRE', 'CENTRE'))] + if 'zz' in stagger_directions: + # for zz include all permutations of ZLOW and CENTRE + dimensions_staggers += [(base_dimensions+'z', ('CENTRE', 'CENTRE', 'ZLOW')), + (base_dimensions+'z', ('CENTRE', 'ZLOW', 'CENTRE')), + (base_dimensions+'z', ('ZLOW', 'CENTRE', 'CENTRE')), + (base_dimensions+'z', ('ZLOW', 'ZLOW', 'CENTRE')), + (base_dimensions+'z', ('ZLOW', 'CENTRE', 'ZLOW')), + (base_dimensions+'z', ('CENTRE', 'ZLOW', 'ZLOW'))] + elif 'z' in stagger_directions: + dimensions_staggers += [(base_dimensions+'z', ('CENTRE', 'ZLOW', 'ZLOW')), + (base_dimensions+'z', ('ZLOW', 'CENTRE', 'CENTRE'))] + result = [] + for ftypes in types: + for dimensions, stagger in dimensions_staggers: + result += test_operator2(ngrids, testfunc1, testfunc2, dimensions, boutcore_operator, symbolic_operator, order, ftypes, method, stagger) + + if test_throw: + # check that unsupported combinations of locations throw an exception + + locations = ['CENTRE', 'XLOW', 'YLOW', 'ZLOW'] + # first make a list of all permutations + fail_staggers = [(x,y,z) for x in locations for y in locations for z in locations] + # now remove the ones have already tested + for dimensions, stagger in dimensions_staggers: + if stagger is not None: + index = fail_staggers.index(stagger) + del fail_staggers[index] + boutcore.setOption('failmesh:nx', '8', force=True) + boutcore.setOption('failmesh:ny', '4', force=True) + boutcore.setOption('failmesh:nz', '4', force=True) + boutcore.setOption('failmesh:staggergrids', 'true', force=True) + failmesh = boutcore.Mesh(section='failmesh') + for stagger in fail_staggers: + # check that an exception is throw for combinations of directions that we expect to fail + try: + test = test_operator2(numpy.array([4, 8]), testfunc1, testfunc2, dimensions, boutcore_operator, symbolic_operator, order, ftypes, method, stagger, mesh_in=failmesh) + except RuntimeError: + result += ['pass'] + else: + result += ['Expected '+str(boutcore_operator)+' to throw for '+stagger[0]+','+stagger[1]+'->'+stagger[2]+' '+str(ftypes)+' '+str(method)+' but it did not.'] + + return result + +min_exponent = 6 +max_exponent = 7 +ngrids = numpy.logspace(min_exponent, max_exponent, num=max_exponent-min_exponent+1, base=2).astype(int) +default_n = 4 +mxg = 2 +myg = 2 +testfunc = 'cos(2*pi*x+y+z)' +testfunc2 = 'sin(4*pi*x+2*y+2*z)+cos(2*pi*x-z)' +order = 2 +plot_error = False +test_deriv_ops = full_test +tests_3d = full_test +test_throw = full_test + +if test_throw: + if boutcore.bout_CHECK < 1: + print('Warning: CHECK='+str(boutcore.bout_CHECK)+' so exceptions will ' + 'not be thrown for unsupported locations. Setting ' + 'test_throw=False...') + test_throw = False + +boutcore.init('-q -q -q -q') + +results = [] + +# single-argument operators +all_types = ('2D', '3D') # eventually, should be able to use this when Field2D operators support staggering properly +type_3d = ('3D',) +type_2d = ('2D',) +results += cycle_staggering('y', 'y', ngrids, testfunc, boutcore.Grad_par, Grad_par, order, type_3d) # staggering in y-direction allowed +results += cycle_staggering('', 'y', ngrids, testfunc, boutcore.Grad_par, Grad_par, order, type_2d) # no staggering allowed +results += cycle_staggering('y', 'y', ngrids, testfunc, boutcore.Div_par, Div_par, order, type_3d) # staggering in y-direction allowed +results += cycle_staggering('', 'y', ngrids, testfunc, boutcore.Div_par, Div_par, order, type_2d) # no staggering allowed +results += cycle_staggering('y', 'y', ngrids, testfunc, boutcore.Grad2_par2, Grad2_par2, order, type_3d) # staggering in y-direction allowed +results += cycle_staggering('', 'y', ngrids, testfunc, boutcore.Grad2_par2, Grad2_par2, order, type_2d) # no staggering allowed +if tests_3d: + results += cycle_staggering('', 'xyz', ngrids, testfunc, boutcore.Laplace, Laplace, order, all_types) # no staggering allowed +results += cycle_staggering('y', 'y', ngrids, testfunc, boutcore.Laplace_par, Laplace_par, order, type_3d) # staggering in y-direction allowed +results += cycle_staggering('', 'y', ngrids, testfunc, boutcore.Laplace_par, Laplace_par, order, type_2d) # no staggering allowed +# note Laplace_perp uses Laplace, so needs y-dimension refinement to converge +if tests_3d: + results += cycle_staggering('', 'xyz', ngrids, testfunc, boutcore.Laplace_perp, Laplace_perp, order, all_types) # no staggering allowed +# Delp2 uses the global mesh, which we can't reset, so can't test here +#results += cycle_staggering('x', 'xz', ngrids, testfunc, boutcore.Delp2, Delp2, order) + +# two-argument operators +all_types2 = [('2D', '2D'), ('3D', '3D')] # expand this to include mixed 2D/3D types at some point +types_2d = [('2D', '2D')] +types_3d = [('3D', '3D')] +results += cycle_staggering2('y', 'y', ngrids, testfunc, testfunc2, boutcore.Vpar_Grad_par, Vpar_Grad_par, order, types_3d) # some staggering in y-direction allowed +results += cycle_staggering2('y', 'y', ngrids, testfunc, testfunc2, boutcore.Vpar_Grad_par, Vpar_Grad_par, order, types_2d) # no staggering allowed +results += cycle_staggering2('yy', 'y', ngrids, testfunc, testfunc2, boutcore.Div_par_K_Grad_par, Div_par_K_Grad_par, order, types_3d) # any staggering in y-direction allowed +results += cycle_staggering2('', 'y', ngrids, testfunc, testfunc2, boutcore.Div_par_K_Grad_par, Div_par_K_Grad_par, order, types_2d) # no staggering allowed +results += cycle_staggering2('y', 'y', ngrids, testfunc, testfunc2, boutcore.Div_par_flux, lambda v,f: Div_par(v*f), 1, types_3d) # some staggering in y-direction allowed +results += cycle_staggering2('', 'y', ngrids, testfunc, testfunc2, boutcore.Div_par_flux, lambda v,f: Div_par(v*f), 1, types_2d) # no staggering allowed +results += cycle_staggering2('y', 'y', ngrids, testfunc, testfunc2, boutcore.Div_par_flux, lambda v,f: Div_par(v*f), order, types_3d, method='C2') # some staggering in y-direction allowed +results += cycle_staggering2('', 'y', ngrids, testfunc, testfunc2, boutcore.Div_par_flux, lambda v,f: Div_par(v*f), order, types_2d, method='C2') # no staggering allowed +# note bracket(Field2D, Field2D) is exactly zero, so doesn't make sense to MMS test +results += cycle_staggering2('', 'xz', ngrids, testfunc, testfunc2, boutcore.bracket, bracket, order, types_3d, method='BRACKET_ARAKAWA') # no staggering allowed +# Note BRACKET_STD version of bracket includes parallel derivatives, so needs +# y-dimension refinement to converge. +# Also it converges faster than 2nd order at 64->128, but approaches closer +# when resolution is increased. Allow test to pass anyway by increasing +# expected order +if tests_3d: + results += cycle_staggering2('', 'xyz', ngrids, testfunc, testfunc2, boutcore.bracket, lambda a,b: b0xGrad_dot_Grad(a,b)/metric.B, 2.6, types_3d, method='BRACKET_STD') # no staggering allowed + +if test_deriv_ops: + # test derivative operators + results += cycle_staggering('x', 'x', ngrids, testfunc, boutcore.DDX, DDX, order, type_3d) + results += cycle_staggering('', 'x', ngrids, testfunc, boutcore.DDX, DDX, order, type_2d) + results += cycle_staggering('y', 'y', ngrids, testfunc, boutcore.DDY, DDY, order, type_3d) + results += cycle_staggering('', 'y', ngrids, testfunc, boutcore.DDY, DDY, order, type_2d) + results += cycle_staggering('', 'z', ngrids, testfunc, boutcore.DDZ, DDZ, order, type_3d, method='C2') + results += cycle_staggering('x', 'x', ngrids, testfunc, boutcore.D2DX2, D2DX2, order, type_3d) + results += cycle_staggering('', 'x', ngrids, testfunc, boutcore.D2DX2, D2DX2, order, type_2d) + results += cycle_staggering('y', 'y', ngrids, testfunc, boutcore.D2DY2, D2DY2, order, type_3d) + results += cycle_staggering('', 'y', ngrids, testfunc, boutcore.D2DY2, D2DY2, order, type_2d) + results += cycle_staggering('', 'z', ngrids, testfunc, boutcore.D2DZ2, D2DZ2, order, type_3d, method='C2') + results += cycle_staggering('', 'x', ngrids, testfunc, boutcore.D4DX4, D4DX4, order, type_3d) + results += cycle_staggering('', 'x', ngrids, testfunc, boutcore.D4DX4, D4DX4, order, type_2d) + results += cycle_staggering('', 'y', ngrids, testfunc, boutcore.D4DY4, D4DY4, order, type_3d) + results += cycle_staggering('', 'y', ngrids, testfunc, boutcore.D4DY4, D4DY4, order, type_2d) + results += cycle_staggering('', 'z', ngrids, testfunc, boutcore.D4DZ4, D4DZ4, order, type_3d) # D4DZ4 is hard coded to use DIFF_C2 + results += cycle_staggering2('x', 'x', ngrids, testfunc, testfunc2, boutcore.VDDX, lambda v,f: v*DDX(f), order, types_3d) + results += cycle_staggering2('', 'x', ngrids, testfunc, testfunc2, boutcore.VDDX, lambda v,f: v*DDX(f), order, types_2d) + results += cycle_staggering2('y', 'y', ngrids, testfunc, testfunc2, boutcore.VDDY, lambda v,f: v*DDY(f), order, types_3d) + results += cycle_staggering2('y', 'y', ngrids, testfunc, testfunc2, boutcore.VDDY, lambda v,f: v*DDY(f), order, types_2d) + results += cycle_staggering2('z', 'z', ngrids, testfunc, testfunc2, boutcore.VDDZ, lambda v,f: v*DDZ(f), order, types_3d, method='C2') + results += cycle_staggering2('x', 'x', ngrids, testfunc, testfunc2, boutcore.FDDX, lambda v,f: DDX(v*f), 1, types_3d) + results += cycle_staggering2('', 'x', ngrids, testfunc, testfunc2, boutcore.FDDX, lambda v,f: DDX(v*f), 1, types_2d) + results += cycle_staggering2('y', 'y', ngrids, testfunc, testfunc2, boutcore.FDDY, lambda v,f: DDY(v*f), 1, types_3d) + results += cycle_staggering2('', 'y', ngrids, testfunc, testfunc2, boutcore.FDDY, lambda v,f: DDY(v*f), 1, types_2d) + results += cycle_staggering2('z', 'z', ngrids, testfunc, testfunc2, boutcore.FDDZ, lambda v,f: DDZ(v*f), 1, types_3d) + results += cycle_staggering('y', 'xy', ngrids, testfunc, boutcore.D2DXDY, D2DXDY, order, type_3d) + results += cycle_staggering('', 'xy', ngrids, testfunc, boutcore.D2DXDY, D2DXDY, order, type_2d) + results += cycle_staggering('', 'xz', ngrids, testfunc, boutcore.D2DXDZ, D2DXDZ, order, type_3d) + results += cycle_staggering('', 'yz', ngrids, testfunc, boutcore.D2DYDZ, D2DYDZ, order, type_3d) + +# check results of tests +fail = False +for result in results: + if result is not 'pass': + print(result) + fail = True +if fail: + exit(1) +else: + print('pass') + exit(0) diff --git a/tests/MMS/difops_short/data/BOUT.inp b/tests/MMS/difops_short/data/BOUT.inp new file mode 100644 index 0000000000..2a4d87d691 --- /dev/null +++ b/tests/MMS/difops_short/data/BOUT.inp @@ -0,0 +1,16 @@ +# partial input file, further options are set using boutcore +# interface in runtest + +[mesh] +nx = 6 +ny = 2 +nz = 1 + +[testmesh:ddx] +upwind = U2 + +[testmesh:ddy] +upwind = U2 + +[testmesh:ddz] +upwind = U2 diff --git a/tests/MMS/difops_short/runtest b/tests/MMS/difops_short/runtest new file mode 100755 index 0000000000..e6a1c47b36 --- /dev/null +++ b/tests/MMS/difops_short/runtest @@ -0,0 +1,8 @@ +#!/bin/bash + +#requires boutcore +#requires not make + +../difops/runtest.py --short + +exit # exit with status of last command diff --git a/tests/integrated/test-squash/.gitignore b/tests/integrated/test-squash/.gitignore new file mode 100644 index 0000000000..a12058acc2 --- /dev/null +++ b/tests/integrated/test-squash/.gitignore @@ -0,0 +1,2 @@ +*.nc +squash \ No newline at end of file diff --git a/tests/integrated/test-squash/data/BOUT.inp b/tests/integrated/test-squash/data/BOUT.inp new file mode 100644 index 0000000000..3253314f8d --- /dev/null +++ b/tests/integrated/test-squash/data/BOUT.inp @@ -0,0 +1,22 @@ +timestep = 1. +nout = 1 + +MZ = 1 + +[mesh] +MXG=1 +MYG=1 +nx = 4 +ny = 2 + +dx = 1. +dy = 1. + +[solver] + +[f2] +scale = 1. +function = 0. + +[f3] +function = 0. diff --git a/tests/integrated/test-squash/makefile b/tests/integrated/test-squash/makefile new file mode 100644 index 0000000000..cecff798ec --- /dev/null +++ b/tests/integrated/test-squash/makefile @@ -0,0 +1,6 @@ + +BOUT_TOP = ../../.. + +SOURCEC = squash.cxx + +include $(BOUT_TOP)/make.config diff --git a/tests/integrated/test-squash/runtest b/tests/integrated/test-squash/runtest new file mode 100755 index 0000000000..cc440597bf --- /dev/null +++ b/tests/integrated/test-squash/runtest @@ -0,0 +1,101 @@ +#!/usr/bin/env python3 + +from boututils.datafile import DataFile +import itertools +import time +import numpy as np +from boututils.run_wrapper import launch_safe, shell_safe + +#requires: all_tests +#requires: netcdf + + +class timer(object): + """Context manager for printing how long a command took + + """ + def __init__(self, msg): + self.msg = msg + + def __enter__(self): + self.start = time.time() + + def __exit__(self, exc_type, exc_value, traceback): + end = time.time() + print("{:12.8f}s {}".format(end - self.start, self.msg)) + + +def timed_shell_safe(cmd, *args, **kwargs): + """Wraps shell_safe in a timer + + """ + with timer(cmd): + shell_safe(cmd, *args, **kwargs) + + +def timed_launch_safe(cmd, *args, **kwargs): + """Wraps launch_safe in a timer + + """ + with timer(cmd): + launch_safe(cmd, *args, **kwargs) + + +def verify(f1, f2): + """Verifies that two BOUT++ files are identical + + """ + with timer("verify %s %s" % (f1, f2)): + d1 = DataFile(f1) + d2 = DataFile(f2) + for v in d1.keys(): + if d1[v].shape != d2[v].shape: + raise RuntimeError("shape mismatch in ", v, d1[v], d2[v]) + if v in ["MXSUB", "MYSUB", "NXPE", "NYPE", "iteration"]: + continue + if not np.allclose(d1[v], d2[v]): + err = "" + dimensions = [range(x) for x in d1[v].shape] + for i in itertools.product(*dimensions): + if d1[v][i] != d2[v][i]: + err += "{}: {} != {}\n".format(i, d1[v][i], d2[v][i]) + raise RuntimeError("data mismatch in ", v, err, d1[v], d2[v]) + + +timed_shell_safe("make") + +# Run once to get normal data +timed_shell_safe("./squash -q -q -q nout=2") +timed_shell_safe("mv data/BOUT.dmp.0.nc f1.nc") + +# Parallel test +timed_shell_safe("rm -f f2.nc") +timed_launch_safe("./squash -q -q -q nout=2", nproc=4, mthread=1) +timed_shell_safe("../../../bin/bout-squashoutput -qdcl 9 data --outputname ../f2.nc") + +verify("f1.nc", "f2.nc") + +# Parallel and in two pieces +timed_shell_safe("rm -f f2.nc") +timed_launch_safe("./squash -q -q -q", nproc=4, mthread=1) +timed_shell_safe("../../../bin/bout-squashoutput -qdcl 9 data --outputname ../f2.nc") +timed_launch_safe("./squash -q -q -q restart", nproc=4, mthread=1) +timed_shell_safe("../../../bin/bout-squashoutput -qdcal 9 data --outputname ../f2.nc") + +verify("f1.nc", "f2.nc") + +# Parallel and in two pieces without dump_on_restart +timed_shell_safe("rm -f f2.nc") +timed_launch_safe("./squash -q -q -q", nproc=4, mthread=1) +timed_shell_safe("../../../bin/bout-squashoutput -qdcl 9 data --outputname ../f2.nc") +timed_launch_safe("./squash -q -q -q restart dump_on_restart=false", nproc=4, mthread=1) +timed_shell_safe("../../../bin/bout-squashoutput -qdcal 9 data --outputname ../f2.nc") + +verify("f1.nc", "f2.nc") + +# Sequential test +timed_shell_safe("rm -f f2.nc") +timed_shell_safe("./squash -q -q -q nout=2") +timed_shell_safe("../../../bin/bout-squashoutput -qdcl 9 data --outputname ../f2.nc") + +verify("f1.nc", "f2.nc") diff --git a/tests/integrated/test-squash/squash.cxx b/tests/integrated/test-squash/squash.cxx new file mode 100644 index 0000000000..98528ba860 --- /dev/null +++ b/tests/integrated/test-squash/squash.cxx @@ -0,0 +1,30 @@ +/* + */ + +#include + +class SquashRun : public PhysicsModel { +protected: + // Initialisation + int init(bool restarting) { + solver->add(f2, "f2"); + solver->add(f3, "f3"); + return 0; + } + + // Calculate time-derivatives + int rhs(BoutReal t) { + ddt(f2) = 1; + ddt(f3) = -1; + f2.applyBoundary(); + f3.applyBoundary(); + return 0; + } + +private: + Field2D f2; + Field3D f3; +}; + +// Create a default main() +BOUTMAIN(SquashRun); diff --git a/tools/idllib/collect.pro b/tools/idllib/collect.pro index 6fcc78fb90..7c30b22bb5 100644 --- a/tools/idllib/collect.pro +++ b/tools/idllib/collect.pro @@ -20,6 +20,8 @@ FUNCTION collect, arg, xind=xind, yind=yind, zind=zind, tind=tind, $ path=path, var=var, t_array=t_array, use=use, old=old, $ quiet=quiet, debug=debug, prefix=prefix + MESSAGE, "This is currently broken for BOUT++ > v4.0.0. See issue #394" + IF NOT KEYWORD_SET(prefix) THEN prefix="BOUT.dmp" IF NOT KEYWORD_SET(debug) THEN BEGIN diff --git a/tools/mathematicalib/BoutCollect.m b/tools/mathematicalib/BoutCollect.m index b664812092..ae92fde2ed 100644 --- a/tools/mathematicalib/BoutCollect.m +++ b/tools/mathematicalib/BoutCollect.m @@ -7,6 +7,8 @@ {Xind,Yind,Zind,Tind,Path,Yguards,Info,Prefix, varnameissymbol,vars,position,dimensions,nxpe,nype,mxsub,mysub,mxg,myg,mz,tarray,files,nfiles,data,tempdata,ts,te,xs,xe,ys,ye,zs,ze,localx,localy,import,lxs,lxe,lys,lye}, + Throw["This is currently broken for BOUT++ > v4.0.0. See issue #394"] + Xind=OptionValue[xind]; Yind=OptionValue[yind]; Zind=OptionValue[zind]; diff --git a/tools/matlablib/import_data_netcdf.m b/tools/matlablib/import_data_netcdf.m index 31795f0d19..c3c3772840 100644 --- a/tools/matlablib/import_data_netcdf.m +++ b/tools/matlablib/import_data_netcdf.m @@ -21,6 +21,8 @@ % Last two variables important only for [X,Y,Z,T] format and any number % will be ok if we wish to plot [X,Y,Z] and [X,Y] type data. +error("This is currently broken for BOUT++ > v4.0.0. See issue #394") + % Check input arguments if ( nargin < 4 ) fprintf('\tBoth dump file path and variable name are requisite input arguments.\n'); diff --git a/tools/matlablib/import_dmp.m b/tools/matlablib/import_dmp.m index adadd93e8b..e5baa7f416 100644 --- a/tools/matlablib/import_dmp.m +++ b/tools/matlablib/import_dmp.m @@ -9,6 +9,8 @@ % % Coded by Minwoo Kim(Mar. 2012) +error("This is currently broken for BOUT++ > v4.0.0. See issue #394") + % Check input arguments if ( nargin < 2 ) fprintf('\tBoth dump file path and variable name are requisite input arguments.\n'); diff --git a/tools/octave/bcollect.m b/tools/octave/bcollect.m index 0005f4c726..d1128daf19 100644 --- a/tools/octave/bcollect.m +++ b/tools/octave/bcollect.m @@ -17,6 +17,7 @@ # Collect metadata from collection of data files function desc = bcollect(path) + error("This is currently broken for BOUT++ > v4.0.0. See issue #394") narg = nargin(); if (narg < 1) # No path specified, so use current directory diff --git a/tools/pylib/_boutcore_build/boutcore.pyx.in b/tools/pylib/_boutcore_build/boutcore.pyx.in index 6e20f61ae6..8f4b6caeda 100755 --- a/tools/pylib/_boutcore_build/boutcore.pyx.in +++ b/tools/pylib/_boutcore_build/boutcore.pyx.in @@ -60,6 +60,11 @@ cimport numpy as np cimport resolve_enum as benum from libc.stdlib cimport malloc, free import copy + +cdef extern from "helper.h": + int _bout_check + +bout_CHECK = _bout_check EOF # make a list of the format "nx, ny, nz" or something similar @@ -83,10 +88,12 @@ do fdd="f3d" ndim=3 dims=(x y z) + ftypeFromObj="f3dFromObj" elif [ $ftype = "Field2D" ]; then fdd="f2d" ndim=2 dims=(x y) + ftypeFromObj="f2dFromObj" else echo "Error, unimplemented ftype" exit 1 @@ -605,6 +612,37 @@ cat <fact).cobj.create${ndim}D(str_,0,0 + ,outloc_,time)) + + EOF done cat < f).cobj[0]) + if isinstance(f, Field3D): + fg.add(( f).cobj[0]) + elif isinstance(f, Field2D): + fg.add(( f).cobj[0]) + else: + raise ValueError("Cannot communicate type "+type(f)) self.cobj.communicate(fg[0]) del fg return self - @property - def coordinates(self): + def coordinates(self, loc=None): """ Get the Coordinates object of this mesh """ - if self._coords is None: - self._coords = coordsFromObj(self.cobj.coordinates()) - return self._coords + cdef benum.CELL_LOC loc_ + if loc is None: + loc_= benum.resolve_cell_loc('CENTRE') + if self._coords is None: + self._coords = coordsFromObj(self.cobj.coordinates(loc_)) + return self._coords + else: + loc_= benum.resolve_cell_loc(loc) + return coordsFromObj(self.cobj.coordinates(loc_)) + + def getXProcIndex(self): + """ + Returns the x-index of the processor + """ + return self.cobj.getXProcIndex() + + def getYProcIndex(self): + """ + Returns the y-index of the processor + """ + return self.cobj.getYProcIndex() cdef Coordinates coordsFromObj(c.Coordinates * obj): coords = Coordinates() @@ -777,7 +837,7 @@ cdef class Coordinates: def _setmembers(self): EOF -for f in "dx" "dy" "J" "Bxy" "g11" "g22" "g33" "g12" "g13" "g23" "g_11" "g_22" "g_33" "g_12" "g_13" "g_23" "G1_11" "G1_22" "G1_33" "G1_12" "G1_13" "G1_23" "G2_11" "G2_22" "G2_33" "G2_12" "G2_13" "G2_23" "G3_11" "G3_22" "G3_33" "G3_12" "G3_13" "G3_23" "G1" "G2" "G3" "ShiftTorsion" "IntShiftTorsion" +for f in "dx" "dy" "J" "Bxy" "g11" "g22" "g33" "g12" "g13" "g23" "g_11" "g_22" "g_33" "g_12" "g_13" "g_23" "G1_11" "G1_22" "G1_33" "G1_12" "G1_13" "G1_23" "G2_11" "G2_22" "G2_33" "G2_12" "G2_13" "G2_23" "G3_11" "G3_22" "G3_33" "G3_12" "G3_13" "G3_23" "G1" "G2" "G3" "ShiftTorsion" do echo " self.${f} = f2dFromObj(self.cobj.${f})" done @@ -787,6 +847,11 @@ do done cat <<"EOF" + # IntShiftTorsion is not always initialized. It should only be loaded + # if localmesh->IncIntShear is true, but localmesh is a private member + # of Coordinates which we cannot access here. + #self.IntShiftTorsion = f2dFromObj(self.cobj.IntShiftTorsion) + cdef class Laplacian: """ Laplacian inversion solver @@ -1097,13 +1162,50 @@ def checkInit(): EOF -f_desc_f="field : Field3D - The Field3D object of which to calculate the derivative" -f_desc_vf="field: Field3D - The Field3D object of which to calculate the derivative - velocity : Field3D - The Field3D object of which the field is advected" -fun () { +# note extra arguments to be passed to the function can be given as the third +# argument, but must start with a comma +dispatch_field_type1 () { +cat <${2}).cobj[0]${3})) + elif isinstance(${2}, Field2D): + return f2dFromObj(c.${1}((${2}).cobj[0]${3})) + else: + raise NotImplementedError("In ${1}: unexpected argument type '"+str(type(${2}))+"' - not supported (yet?).") +EOF +} + +dispatch_field_type_to_BoutReal () { +cat <${2}).cobj[0]${3}) + elif isinstance(${2}, Field2D): + return c.${1}((${2}).cobj[0]${3}) + else: + raise NotImplementedError("In ${1}: unexpected argument type '"+str(type(${2}))+"' - not supported (yet?).") +EOF +} + +# note extra arguments to be passed to the function can be given as the fourth +# argument, but must start with a comma +dispatch_field_type2 () { +cat <${2}).cobj[0], (${3}).cobj[0]${4})) + elif isinstance(${2}, Field2D) and isinstance(${3}, Field2D): + return f2dFromObj(c.${1}((${2}).cobj[0], (${3}).cobj[0]${4})) + else: + raise NotImplementedError("In ${1}: unexpected argument types '"+str(type(${2}))+" and "+str(type(${3}))+"' - not supported (yet?).") +EOF +} + +f_desc_f="field : Field3D Field2D + The Field3D/Field2D object of which to calculate the derivative" +f_desc_vf="field: Field3D Field2D + The Field3D/Field2D object of which to calculate the derivative + velocity : Field3D Field2D + The Field3D/Field2D object by which the field is advected" +fun1 () { cat <a).cobj[0])) - elif isinstance(a, Field2D): - return f2dFromObj(c.$fun((a).cobj[0])) - else: - raise NotImplementedError("In $fun: unexpected argument type '"+str(type(a))+"' - not supported (yet?).") +$(dispatch_field_type1 $fun a) EOF done -cat <<"EOF" - -def pow(Field3D a, exponent): - """ - Returns a**e where a is a Field3D and e is a number - - Parameters - ---------- - a : Field3D - The field for which to calculate the power - exponent : float - The exponent - - Returns - ------- - Field3D - The a**exponent - """ - return f3dFromObj(c.pow(a.cobj[0],float(exponent))) -def min(Field3D a): +for fun in min max +do +cat <fact).cobj.create3D(str_,0,0 - ,outloc_,time)) +$(dispatch_field_type1 pow a ", float(exponent)") -def interp_to(Field3D f3d,location): +def interp_to(f,location): """ - Interpolate a Field3D to a given location + Interpolate a Field3D/Field2D to a given location Parameters ---------- - f3d : Field3D + f3d : Field3D Field2D The field to interpolate location : string The location to which to interploate Returns ------- - Field3D + Field3D Field2D the interpolated field """ checkInit() cdef benum.CELL_LOC location_ = benum.resolve_cell_loc(location) - return f3dFromObj(c.interp_to(f3d.cobj[0],location_)) - +$(dispatch_field_type1 interp_to f ", location_") def setOption(name, value, source="PyInterface", force=False): """ Set an option in the global Options tree. Prefer - `Options.set` to avoid unexpected results if several Option + 'Options.set' to avoid unexpected results if several Option roots are avalaible. Parameters ---------- name : string the name of the value to be set. Can be relative, - e.g. `mesh:ddx:first`. + e.g. 'mesh:ddx:first'. value : string the value to be set source : string @@ -1407,7 +1486,7 @@ def setOption(name, value, source="PyInterface", force=False): track of where what was set. force : bool If a value is overwritten, an exception is - thrown. setting this to `True` avoids the exception. + thrown. setting this to 'True' avoids the exception. """ checkInit() root=Options('') @@ -1463,7 +1542,7 @@ cdef class Options: ---------- name : string the name of the value to be set. Can be relative, - e.g. `mesh:ddx:first`. + e.g. 'mesh:ddx:first'. value : string the value to be set source : string @@ -1471,7 +1550,7 @@ cdef class Options: track of where what was set. force : bool If a value is overwritten, an exception is - thrown. setting this to `True` avoids the exception. + thrown. setting this to 'True' avoids the exception. """ cdef c.Options * opt=self.cobj cdef c.string sec_ @@ -1493,7 +1572,7 @@ cdef class Options: ---------- name : string the name of the value to get. Can be relative, - e.g. `mesh:ddx:first`. + e.g. 'mesh:ddx:first'. default : bool, string or float Depending on the type of the default, different things will be returned. Supported types are bool, string or float diff --git a/tools/pylib/_boutcore_build/boutcpp.pxd.in b/tools/pylib/_boutcore_build/boutcpp.pxd.in index 49edc12332..42b9f62683 100644 --- a/tools/pylib/_boutcore_build/boutcpp.pxd.in +++ b/tools/pylib/_boutcore_build/boutcpp.pxd.in @@ -60,7 +60,7 @@ cdef extern from "bout/mesh.hxx": int ystart int LocalNx int LocalNy - Coordinates * coordinates() + Coordinates * coordinates(benum.CELL_LOC) cdef extern from "bout/coordinates.hxx": cppclass Coordinates: @@ -86,6 +86,7 @@ cdef extern from "bout/fieldgroup.hxx": cppclass FieldGroup: FieldGroup() void add(Field3D&) + void add(Field2D&) cdef extern from "invert_laplace.hxx": cppclass Laplacian: @staticmethod @@ -103,12 +104,24 @@ cdef extern from "invert_laplace.hxx": void setCoefEz(Field3D) cdef extern from "difops.hxx": - Field3D Div_par(Field3D, benum.CELL_LOC, benum.DIFF_METHOD) - Field3D Grad_par(Field3D, benum.CELL_LOC, benum.DIFF_METHOD) - Field3D Laplace(Field3D) - Field3D Vpar_Grad_par(Field3D, Field3D, benum.CELL_LOC, benum.DIFF_METHOD) - Field3D bracket(Field3D,Field3D, benum.BRACKET_METHOD, benum.CELL_LOC) - Field3D Delp2(Field3D,double) +EOF +for ftype in Field3D Field2D +do + cat < #include +#ifndef CHECK +int _bout_check = 0; +#else +int _bout_check = CHECK; +#endif + EOF for ftype in "Field3D" "Field2D" do diff --git a/tools/pylib/boutdata/mms_alternate.py b/tools/pylib/boutdata/mms_alternate.py new file mode 100644 index 0000000000..d89e333b84 --- /dev/null +++ b/tools/pylib/boutdata/mms_alternate.py @@ -0,0 +1,684 @@ +""" Functions for calculating sources for the + Method of Manufactured Solutions (MMS) + +""" +from __future__ import print_function +from __future__ import division +#from builtins import str +#from builtins import object + +from sympy import symbols, cos, sin, diff, sqrt, pi, simplify, trigsimp, Wild, integrate + +from numpy import arange, newaxis +#from numpy import arange, zeros + +global metric + +# Constants +qe = 1.602e-19 +Mp = 1.67262158e-27 +mu0 = 4.e-7*3.141592653589793 + +# Define symbols + +x = symbols('x') +y = symbols('y') +z = symbols('z') +t = symbols('t') + +class Metric(object): + def __init__(self): + # Create an identity metric + self.x = x + self.y = y + self.z = z + + self.g11 = self.g22 = self.g33 = 1.0 + self.g12 = self.g23 = self.g13 = 0.0 + + self.g_11 = self.g_22 = self.g_33 = 1.0 + self.g_12 = self.g_23 = self.g_13 = 0.0 + + self.J = 1.0 + self.B = 1.0 + + # coordinate transformations for 'mesh refinement' + # expressions for these must average to 1. + self.scalex = 1.0 + self.scaley = 1.0 + +metric = Metric() + +# Basic differencing +def ddt(f): + """Time derivative""" + return diff(f, t) + + +def DDX(f): + # psiwidth = dx/dx_in + return diff(f, metric.x)/metric.psiwidth/metric.scalex + +def DDY(f): + return diff(f, metric.y)/metric.scaley + +def DDZ(f): + return diff(f, metric.z)*metric.zperiod + + +def D2DX2(f): + return DDX(DDX(f)) + +def D2DY2(f): + return DDY(DDY(f)) + +def D2DZ2(f): + return DDZ(DDZ(f)) + + +# don't include derivatives of scalex/scaley, to match BOUT++ implementation +# where D4D*4 are used just for numerical 'hyperdiffusion' +def D4DX4(f): + return diff(f, metric.x, 4)/metric.psiwidth**4/metric.scalex**4 + +def D4DY4(f): + return diff(f, metric.y, 4)/metric.scaley**4 + +def D4DZ4(f): + return diff(f, metric.z, 4)*metric.zperiod**4 + + +def D2DXDY(f): + return DDX(DDY(f)) + +def D2DXDZ(f): + return DDX(DDZ(f)) + +def D2DYDZ(f): + return DDY(DDZ(f)) + +# Operators + +def bracket(f, g): + """ + Calculates [f,g] symbolically + """ + + dfdx = DDX(f) + dfdz = DDZ(f) + + dgdx = DDX(g) + dgdz = DDZ(g) + + return dfdz * dgdx - dfdx * dgdz + +def b0xGrad_dot_Grad(phi, A): + """ + Perpendicular advection operator, including + derivatives in y + + Note: If y derivatives are neglected, then this reduces + to bracket(f, g) * metric.B + (in a Clebsch coordinate system) + + """ + dpdx = DDX(phi) + dpdy = DDY(phi) + dpdz = DDZ(phi) + + vx = metric.g_22*dpdz - metric.g_23*dpdy; + vy = metric.g_23*dpdx - metric.g_12*dpdz; + vz = metric.g_12*dpdy - metric.g_22*dpdx; + + return (+ vx*DDX(A) + + vy*DDY(A) + + vz*DDZ(A) ) / (metric.J*sqrt(metric.g_22)) + +def Delp2(f, all_terms=True): + """ Laplacian in X-Z + + If all_terms is false then first derivative terms are neglected. + By default all_terms is true, but can be disabled + in the BOUT.inp file (laplace section) + + """ + d2fdx2 = D2DX2(f) + d2fdz2 = D2DZ2(f) + d2fdxdz = D2DXDZ(f) + + result = metric.g11*d2fdx2 + metric.g33*d2fdz2 + 2.*metric.g13*d2fdxdz + + if all_terms: + result += metric.G1 * DDX(f) + metric.G3 * DDZ(f) + + return result + +def Delp4(f): + d4fdx4 = D2DX2(D2DX2(f)) + d4fdz4 = D2DZ2(D2DZ2(f)) + + return d4fdx4 + d4fdz4 + +def Grad_par(f): + """The parallel gradient""" + return DDY(f) / sqrt(metric.g_22) + +def Grad2_par2(f): + """The parallel 2nd derivative""" + return Grad_par(Grad_par(f)) + +def Vpar_Grad_par(v, f): + """Parallel advection operator v*grad_||(f)""" + return v * Grad_par(f) + +def Div_par_K_Grad_par(K, f): + """Parallel diffusion operator Div(b K b.Grad(f))""" + return Div_par(K*Grad_par(f)) + +def Div_par(f): + ''' + Divergence of magnetic field aligned vector v = \bhat f + \nabla \cdot (\bhat f) = 1/J \partial_y (f/B) + = B Grad_par(f/B) + ''' + return metric.B*Grad_par(f/metric.B) + +def Laplace(f): + """The full Laplace operator""" + result = metric.G1*DDX(f) + metric.G2*DDY(f) + metric.G3*DDZ(f)\ + + metric.g11*D2DX2(f) + metric.g22*D2DY2(f) + metric.g33*D2DZ2(f)\ + + 2.0*(metric.g12*D2DXDY(f) + metric.g13*D2DXDZ(f) + metric.g23*D2DYDZ(f)) + + return result + +def Laplace_par(f): + """ + Div( b (b.Grad(f) ) ) = (1/J) d/dy ( J/g_22 * df/dy ) + """ + return DDY( metric.J/metric.g_22 * DDY(f) )/ metric.J + +def Laplace_perp(f): + """ + The perpendicular Laplace operator + + Laplace_perp = Laplace - Laplace_par + """ + return Laplace(f) - Laplace_par(f) + +# Convert expression to string + +def trySimplify(expr): + """ + Tries to simplify an expression + """ + try: + return simplify(expr) + except ValueError: + return expr + +def exprToStr(expr): + """ + Convert a sympy expression to a string for BOUT++ input + """ + + s = str(expr).replace("**", "^") # Replace exponent operator + + # Try to remove lots of 1.0*... + s = s.replace("(1.0*", "(") + s = s.replace(" 1.0*", " ") + + return s + +def exprMag(expr): + """ + Estimate the magnitude of an expression + + """ + + # Replace all sin, cos with 1 + any = Wild('a') # Wildcard + expr = expr.replace(sin(any), 1.0) + expr = expr.replace(cos(any), 1.0) + + # Pick maximum values of x,y,z + expr = expr.subs(x, 1.0) + expr = expr.subs(y, 2.*pi) + expr = expr.subs(z, 2.*pi) + + return expr.evalf() + +################################## + +class BaseTokamak(object): + """ + Virtual base class defining various useful functions for child *Tokamak classes + + Child class __init__ methods must define the member fields: + self.x + self.y + + # Major radius of axis + self.R + + self.dr + + # Minor radius + self.r + + # Safety factor + self.q + + # Toroidal angle of a field-line as function + # of poloidal angle y + self.zShift + + # Field-line pitch + self.nu + + # Coordinates of grid points + self.Rxy + self.Zxy + + # Poloidal arc length + self.hthe + + # Toroidal magnetic field + self.Btxy + + # Poloidal magnetic field + self.Bpxy + + # Total magnetic field + self.Bxy + + # dx = Bp * R * dr -- width of the box in psi space + self.psiwidth + + # Integrated shear + self.sinty + + # If true, set integrated shear to zero in metric (for shifted-metric schemes) + self.shifted + + # Extra expressions to add to grid file + self._extra + """ + + psiwidth = None + r0 = None + + def __init__(self): + raise ValueError("Error: An instance of BaseTokamak should never be created, use a derived class like SimpleTokamak or ShapedTokamak instead") + + def add(self, expr, name): + """ + Add an additional expression to be written to the grid files + + """ + self._extra[name] = expr + + def write(self, nx, ny, output, MXG=2): + """ + Outputs a tokamak shape to a grid file + + nx - Number of radial grid points, not including guard cells + ny - Number of poloidal (parallel) grid points + output - boututils.datafile object, e.g., an open netCDF file + MXG, Number of guard cells in the x-direction + """ + + ngx = nx + 2*MXG + ngy = ny + + # Create an x and y grid to evaluate expressions on + xarr = (arange(nx + 2*MXG) - MXG + 0.5) / nx + xarr = xarr[:,newaxis] + yarr = 2.*pi*arange(ny)/ny + yarr = yarr[newaxis,:] + + output.write("nx", ngx) + output.write("ny", ngy) + + dx = self.psiwidth * metric.scalex / nx + 0.*self.x + dy = 2.*pi * metric.scaley / ny + 0.*self.x + + for name, var in [ ("dx", dx), + ("dy", dy), + ("Rxy", self.Rxy), + ("Zxy", self.Zxy), + ("Btxy", self.Btxy), + ("Bpxy", self.Bpxy), + ("Bxy", self.Bxy), + ("hthe", self.hthe), + ("sinty", self.sinty), + ("zShift", self.zShift)]: + + varfunc = lambdify([x, y], var) + values = varfunc(xarr, yarr) + ## Note: This is slow, and could be improved using something like lambdify + #values = zeros([ngx, ngy]) + #for i, x in enumerate(xarr): + # for j, y in enumerate(yarr): + # values[i,j] = var.evalf(subs={self.x:x, self.y:y}) + + output.write(name, values) + + for name, var in list(self._extra.items()): + values = zeros([ngx, ngy]) + for i, x in enumerate(xarr): + for j, y in enumerate(yarr): + values[i,j] = var.evalf(subs={self.x:x, self.y:y}) + + output.write(name, values) + + shiftAngle = zeros(ngx) + for i, x in enumerate(xarr): + shiftAngle[i] = self.shiftAngle.evalf(subs={self.x:x}) + + output.write("ShiftAngle", shiftAngle) + + def metric(self): + """ + Calculates an analytic metric tensor + """ + + # Set symbols for x and y directions + metric.x = self.x + metric.y = self.y + + # Calculate metric tensor + sbp = 1. # sign of Bp + if (self.Bpxy.evalf(subs={x:0., y:0.}) < 0.): + sbp = -1. + + metric.g11 = (self.Rxy * self.Bpxy)**2 + metric.g22 = 1/self.hthe**2 + metric.g33 = self.sinty**2*metric.g11 + self.Bxy**2/metric.g11 + metric.g12 = 0*x + metric.g13 = -self.sinty*metric.g11 + metric.g23 = -sbp*self.Btxy / (self.hthe * self.Bpxy * self.Rxy) + + metric.g_11 = 1./metric.g11 + (self.sinty*self.Rxy)**2 + metric.g_22 = (self.Bxy * self.hthe / self.Bpxy)**2 + metric.g_33 = self.Rxy**2 + metric.g_12 = sbp*self.Btxy*self.hthe*self.sinty*self.Rxy / self.Bpxy + metric.g_13 = self.sinty*self.Rxy**2 + metric.g_23 = sbp*self.Btxy*self.hthe*self.Rxy / self.Bpxy + + metric.zShift = self.zShift + metric.shiftAngle = self.shiftAngle + + metric.J = self.hthe / self.Bpxy + metric.B = self.Bxy + + metric.psiwidth = self.psiwidth + metric.zperiod = self.zperiod + + # normalize lengths + if self.r0 is not None: + metric.g11 *= self.r0**2 + metric.g22 *= self.r0**2 + metric.g33 *= self.r0**2 + metric.g12 *= self.r0**2 + metric.g13 *= self.r0**2 + metric.g23 *= self.r0**2 + metric.g_11 /= self.r0**2 + metric.g_22 /= self.r0**2 + metric.g_33 /= self.r0**2 + metric.g_12 /= self.r0**2 + metric.g_13 /= self.r0**2 + metric.g_23 /= self.r0**2 + + self.metric_is_set = True + + # Christoffel symbols + metric.G1_11 = 0.5 * metric.g11 * DDX(metric.g_11) + metric.g12 * (DDX(metric.g_12) - 0.5 * DDY(metric.g_11)) + metric.g13 * (DDX(metric.g_13) - 0.5 * DDZ(metric.g_11)) + metric.G1_22 = metric.g11 * (DDY(metric.g_12) - 0.5 * DDX(metric.g_22)) + 0.5 * metric.g12 * DDY(metric.g_22) + metric.g13 * (DDY(metric.g_23) - 0.5 * DDZ(metric.g_22)) + metric.G1_33 = metric.g11 * (DDZ(metric.g_13) - 0.5 * DDX(metric.g_33)) + metric.g12 * (DDZ(metric.g_23) - 0.5 * DDY(metric.g_33)) + 0.5 * metric.g13 * DDZ(metric.g_33) + metric.G1_12 = 0.5 * metric.g11 * DDY(metric.g_11) + 0.5 * metric.g12 * DDX(metric.g_22) + 0.5 * metric.g13 * (DDY(metric.g_13) + DDX(metric.g_23) - DDZ(metric.g_12)) + metric.G1_13 = 0.5 * metric.g11 * DDZ(metric.g_11) + 0.5 * metric.g12 * (DDZ(metric.g_12) + DDX(metric.g_23) - DDY(metric.g_13)) + 0.5 * metric.g13 * DDX(metric.g_33) + metric.G1_23 = 0.5 * metric.g11 * (DDZ(metric.g_12) + DDY(metric.g_13) - DDX(metric.g_23)) + 0.5 * metric.g12 * (DDZ(metric.g_22) + DDY(metric.g_23) - DDY(metric.g_23)) + 0.5 * metric.g13 * DDY(metric.g_33) + + metric.G2_11 = 0.5 * metric.g12 * DDX(metric.g_11) + metric.g22 * (DDX(metric.g_12) - 0.5 * DDY(metric.g_11)) + metric.g23 * (DDX(metric.g_13) - 0.5 * DDZ(metric.g_11)) + metric.G2_22 = metric.g12 * (DDY(metric.g_12) - 0.5 * DDX(metric.g_22)) + 0.5 * metric.g22 * DDY(metric.g_22) + metric.g23 * (DDY(metric.g23) - 0.5 * DDZ(metric.g_22)) + metric.G2_33 = metric.g12 * (DDZ(metric.g_13) - 0.5 * DDX(metric.g_33)) + metric.g22 * (DDZ(metric.g_23) - 0.5 * DDY(metric.g_33)) + 0.5 * metric.g23 * DDZ(metric.g_33) + metric.G2_12 = 0.5 * metric.g12 * DDY(metric.g_11) + 0.5 * metric.g22 * DDX(metric.g_22) + 0.5 * metric.g23 * (DDY(metric.g_13) + DDX(metric.g_23) - DDZ(metric.g_12)) + metric.G2_13 = 0.5 * metric.g12 * (DDZ(metric.g_11) + DDX(metric.g_13) - DDX(metric.g_13)) + 0.5 * metric.g22 * (DDZ(metric.g_12) + DDX(metric.g_23) - DDY(metric.g_13)) + 0.5 * metric.g23 * DDX(metric.g_33) + metric.G2_23 = 0.5 * metric.g12 * (DDZ(metric.g_12) + DDY(metric.g_13) - DDX(metric.g_23)) + 0.5 * metric.g22 * DDZ(metric.g_22) + 0.5 * metric.g23 * DDY(metric.g_33) + + metric.G3_11 = 0.5 * metric.g13 * DDX(metric.g_11) + metric.g23 * (DDX(metric.g_12) - 0.5 * DDY(metric.g_11)) + metric.g33 * (DDX(metric.g_13) - 0.5 * DDZ(metric.g_11)) + metric.G3_22 = metric.g13 * (DDY(metric.g_12) - 0.5 * DDX(metric.g_22)) + 0.5 * metric.g23 * DDY(metric.g_22) + metric.g33 * (DDY(metric.g_23) - 0.5 * DDZ(metric.g_22)) + metric.G3_33 = metric.g13 * (DDZ(metric.g_13) - 0.5 * DDX(metric.g_33)) + metric.g23 * (DDZ(metric.g_23) - 0.5 * DDY(metric.g_33)) + 0.5 * metric.g33 * DDZ(metric.g_33) + metric.G3_12 = 0.5 * metric.g13 * DDY(metric.g_11) + 0.5 * metric.g23 * DDX(metric.g_22) + 0.5 * metric.g33 * (DDY(metric.g_13) + DDX(metric.g_23) - DDZ(metric.g_12)) + metric.G3_13 = 0.5 * metric.g13 * DDZ(metric.g_11) + 0.5 * metric.g23 * (DDZ(metric.g_12) + DDX(metric.g_23) - DDY(metric.g_13)) + 0.5 * metric.g33 * DDX(metric.g_33) + metric.G3_23 = 0.5 * metric.g13 * (DDZ(metric.g_12) + DDY(metric.g_13) - DDX(metric.g_23)) + 0.5 * metric.g23 * DDZ(metric.g_22) + 0.5 * metric.g33 * DDY(metric.g_33) + + metric.G1 = (DDX(metric.J * metric.g11) + DDY(metric.J * metric.g12) + DDZ(metric.J * metric.g13)) / metric.J + metric.G2 = (DDX(metric.J * metric.g12) + DDY(metric.J * metric.g22) + DDZ(metric.J * metric.g23)) / metric.J + metric.G3 = (DDX(metric.J * metric.g13) + DDY(metric.J * metric.g23) + DDZ(metric.J * metric.g33)) / metric.J + + def print_mesh(self): + """ + Prints the metrics to stdout to be copied to a BOUT.inp file + """ + if not self.metric_is_set: + raise ValueError("Error: metric has not been calculated yet, so cannot print") + + print("dx = "+exprToStr(metric.psiwidth*metric.scalex)+"/(nx-2*mxg)") + print("dy = 2.*pi*"+exprToStr(metric.scaley)+"/ny") + print("dz = 2.*pi/nz/"+exprToStr(metric.zperiod)) + print("g11 = "+exprToStr(metric.g11)) + print("g22 = "+exprToStr(metric.g22)) + print("g33 = "+exprToStr(metric.g33)) + print("g12 = "+exprToStr(metric.g12)) + print("g13 = "+exprToStr(metric.g13)) + print("g23 = "+exprToStr(metric.g23)) + print("g_11 = "+exprToStr(metric.g_11)) + print("g_22 = "+exprToStr(metric.g_22)) + print("g_33 = "+exprToStr(metric.g_33)) + print("g_12 = "+exprToStr(metric.g_12)) + print("g_13 = "+exprToStr(metric.g_13)) + print("g_23 = "+exprToStr(metric.g_23)) + print("J = "+exprToStr(metric.J)) + print("Bxy = "+exprToStr(metric.B)) + print("G1_11 = "+exprToStr(metric.G1_11)) + print("G1_22 = "+exprToStr(metric.G1_22)) + print("G1_33 = "+exprToStr(metric.G1_33)) + print("G1_12 = "+exprToStr(metric.G1_12)) + print("G1_13 = "+exprToStr(metric.G1_13)) + print("G1_23 = "+exprToStr(metric.G1_23)) + print("G2_11 = "+exprToStr(metric.G2_11)) + print("G2_22 = "+exprToStr(metric.G2_22)) + print("G2_33 = "+exprToStr(metric.G2_33)) + print("G2_12 = "+exprToStr(metric.G2_12)) + print("G2_13 = "+exprToStr(metric.G2_13)) + print("G2_23 = "+exprToStr(metric.G2_23)) + print("G3_11 = "+exprToStr(metric.G3_11)) + print("G3_22 = "+exprToStr(metric.G3_22)) + print("G3_33 = "+exprToStr(metric.G3_33)) + print("G3_12 = "+exprToStr(metric.G3_12)) + print("G3_13 = "+exprToStr(metric.G3_13)) + print("G3_23 = "+exprToStr(metric.G3_23)) + print("G1 = "+exprToStr(metric.G1)) + print("G2 = "+exprToStr(metric.G2)) + print("G3 = "+exprToStr(metric.G3)) + + print("Lx = "+exprToStr(metric.psiwidth)) + print("Lz = "+exprToStr((2.*pi*self.R/metric.zperiod).evalf())) + + def set_scalex(self, expr): + # check average of expr is 1. + average = integrate(expr, (x, 0, 1), (y, 0, 2*pi))/2/pi + if not average == 1: + raise ValueError("scalex must average to 1. Got "+str(average)) + else: + metric.scalex = expr + + def set_scaley(self, expr): + # check average of expr is 1. + average = integrate(expr, (x, 0, 1), (y, 0, 2*pi))/2/pi + if not average == 1: + raise ValueError("scaley must average to 1. Got "+str(average)) + else: + metric.scaley = expr + +################################## + +class SimpleTokamak(BaseTokamak): + """ + Simple tokamak + + NOTE: This is NOT an equilibrium calculation. The input + is intended solely for testing with MMS + """ + def __init__(self, R = 2, Bt = 1.0, eps = 0.1, dr=0.02, psiN0=0.5, zperiod=1, r0=1., q = lambda psiN:2+psiN**2, lower_legs_fraction=0., upper_legs_fraction=0., shifted=False): + """ + R - Major radius [metric] + + Bt - Toroidal field [T] + + eps - Inverse aspect ratio + + dr - Width of the radial region [metric] + + q(psiN) - A function which returns the safety factor + as a function of psiN in range [0,1] + + psiN0- Normalized poloidal flux of inner edge of grid + + r0 - Length-scale to normalize metric components (default 1m) + + lower_legs_fraction - what fraction of the poloidal grid is in the lower divertor legs (per leg, so lower_legs_fraction+upper_legs_fraction<=0.5) + + upper_legs_fraction - what fraction of the poloidal grid is in the upper divertor legs (per leg, so lower_legs_fraction+upper_legs_fraction<=0.5) + + Coordinates: + x - Radial, [0,psiwidth], x=psi-psi_inner so dx=dpsi but x=0 at the inner edge of the grid + y - Poloidal, [0,2pi]. Origin is at inboard midplane. + + + """ + + # Have we calculated metric components yet? + self.metric_is_set = False + + self.x = x + self.y = y + self.zperiod = zperiod + self.r0 = r0 + self.R = R + self.dr = dr + self.psiN0 = psiN0 + self.lower_legs_fraction = lower_legs_fraction + self.upper_legs_fraction = upper_legs_fraction + self.shifted = shifted + + # Minor radius + self.r = R * eps + + # Approximate poloidal field for radial width calculation + Bp0 = Bt * self.r*self.psiN0 / (q(self.psiN0) * self.R) + + # dpsi = Bp * R * dr -- width of the box in psi space + self.psiwidth = Bp0 * self.R * self.dr + self.psi0 = Bp0 * R * self.r # value of psi at 'separatrix' taken to be at r, psi=0 at magnetic axis + + # Get safety factor + self.q = q((x + self.psiN0*self.psi0)/self.psi0) + + # Toroidal angle of a field-line as function + # of poloidal angle y + self.zShift = self.q*(self.y-pi + eps * sin(y-pi)) + + # Toroidal shift of field line at branch cut where poloidal angle goes 2pi->0 + self.shiftAngle = 2*pi*self.q + + # Field-line pitch + self.nu = self.q*(1 + eps*cos(y-pi)) #diff(self.zShift, y) + + # Coordinates of grid points + self.Rxy = R - self.r*self.psiN0 * cos(y-pi) + self.Zxy = self.r*self.psiN0 * sin(y-pi) + + # Poloidal arc length + self.hthe = self.r*self.psiN0 + 0.*y + + # Toroidal magnetic field + self.Btxy = Bt * R / self.Rxy + + # Poloidal magnetic field + self.Bpxy = self.Btxy * self.hthe / (self.nu * self.Rxy) + + # Total magnetic field + self.Bxy = sqrt(self.Btxy**2 + self.Bpxy**2) + + # Integrated shear + if shifted: + self.sinty = 0*y + else: + self.sinty = diff(self.zShift, x) + + # Convert all "x" symbols from flux to [0,1] + xsub = metric.x * self.psiwidth + + self.q = self.q.subs(x, xsub) + self.zShift = self.zShift.subs(x, xsub) + self.shiftAngle = self.shiftAngle.subs(x, xsub) + self.nu = self.nu.subs(x, xsub) + self.Rxy = self.Rxy.subs(x, xsub) + self.Zxy = self.Zxy.subs(x, xsub) + self.hthe = self.hthe.subs(x, xsub) + self.Btxy = self.Btxy.subs(x, xsub) + self.Bpxy = self.Bpxy.subs(x, xsub) + self.Bxy = self.Bxy.subs(x, xsub) + self.sinty = self.sinty.subs(x, xsub) + + # calculating grid spacing needs grid szie: initialize as None + self.dx = None + self.dy = None + self.dz = None + + # calculating topology parameters needs grid size: initialize as None + self.ixseps1 = None + self.ixseps2 = None + self.jyseps1_1 = None + self.jyseps1_2 = None + self.jyseps2_1 = None + self.jyseps2_2 = None + + # Extra expressions to add to grid file + self._extra = {} + + # Calculate metric terms + self.metric() + + def setNs(self, nx=None, ny=None, nz=None, mxg=2): + """ + Pass in the grid sizes and hence calculate grid spacing and topology parameters + """ + self.nx = nx + self.ny = ny + self.nz = nz + + core_fraction = 1.-2.*self.lower_legs_fraction-2.*self.upper_legs_fraction + + # calculate grid spacings + self.dx = metric.psiwidth*metric.scalex/self.nx + self.dy = 2.*pi * metric.scaley / (core_fraction * self.ny) + self.dz = 2.*pi/self.nz + + # find position of separatrix (uniform grid in x), i.e. x-index for which psiN=1 + self.ixseps1 = int(self.nx * (1. - self.psiN0)*self.psi0/self.psiwidth)-1 + mxg # numbering for ixseps* includes x-guard cells + self.ixseps2 = self.nx+2*mxg # assume connected double-null if double-null: set ixseps2 to BoutMesh default + + if self.ixseps1 < 0 or self.ixseps1 >= self.nx: + # no separatrix, so no branch cuts: make everywhere 'core', like BoutMesh defaults + self.jyseps1_1 = -1 + self.jyseps1_2 = self.nx/2 + self.jyseps2_1 = self.jyseps1_2 + self.jyseps2_2 = self.ny-1 + else: + self.jyseps1_1 = int(self.lower_legs_fraction*self.ny)-1 + self.jyseps2_2 = int((1.-self.lower_legs_fraction)*self.ny)-1 + self.jyseps1_2 = int(self.ny//2 - self.upper_legs_fraction*self.ny)-1 + self.jyseps2_1 = int(self.ny//2 + self.upper_legs_fraction*self.ny)-1 diff --git a/tools/pylib/boutdata/squashoutput.py b/tools/pylib/boutdata/squashoutput.py index 0c8ac48029..52bcaa40d1 100644 --- a/tools/pylib/boutdata/squashoutput.py +++ b/tools/pylib/boutdata/squashoutput.py @@ -20,6 +20,11 @@ from boututils.boutarray import BoutArray import numpy import os +import gc +import tempfile +import shutil +import glob + def squashoutput(datadir=".", outputname="BOUT.dmp.nc", format="NETCDF4", tind=None, xind=None, yind=None, zind=None, singleprecision=False, compress=False, @@ -72,59 +77,68 @@ def squashoutput(datadir=".", outputname="BOUT.dmp.nc", format="NETCDF4", tind=N Delete the original files after squashing. """ - import gc - - fullpath = os.path.join(datadir,outputname) + fullpath = os.path.join(datadir, outputname) if append: - import tempfile - import shutil - import glob datadirnew = tempfile.mkdtemp(dir=datadir) - for f in glob.glob(datadir+"/BOUT.dmp.*.??"): + for f in glob.glob(datadir + "/BOUT.dmp.*.??"): if not quiet: - print("moving",f) - shutil.move(f,datadirnew) - oldfile=datadirnew+"/"+outputname - datadir=datadirnew + print("moving", f) + shutil.move(f, datadirnew) + oldfile = datadirnew + "/" + outputname + datadir = datadirnew if os.path.isfile(fullpath) and not append: - raise ValueError(fullpath+" already exists. Collect may try to read from this file, which is presumably not desired behaviour.") + raise ValueError( + fullpath + " already exists. Collect may try to read from this file, which is presumably not desired behaviour.") # useful object from BOUT pylib to access output data - outputs = BoutOutputs(datadir, info=False, xguards=True, yguards=True, tind=tind, xind=xind, yind=yind, zind=zind) + outputs = BoutOutputs(datadir, info=False, xguards=True, + yguards=True, tind=tind, xind=xind, yind=yind, zind=zind) outputvars = outputs.keys() # Read a value to cache the files outputs[outputvars[0]] if append: # move only after the file list is cached - shutil.move(fullpath,oldfile) + shutil.move(fullpath, oldfile) t_array_index = outputvars.index("t_array") outputvars.append(outputvars.pop(t_array_index)) - kwargs={} + kwargs = {} if compress: - kwargs['zlib']=True + kwargs['zlib'] = True if least_significant_digit is not None: - kwargs['least_significant_digit']=least_significant_digit + kwargs['least_significant_digit'] = least_significant_digit if complevel is not None: - kwargs['complevel']=complevel + kwargs['complevel'] = complevel if append: - old=DataFile(oldfile) + old = DataFile(oldfile) + # Check if dump on restart was enabled + # If so, we want to drop the duplicated entry + cropnew = 0 + if old['t_array'][-1] == outputs['t_array'][0]: + cropnew = 1 + # Make sure we don't end up with duplicated data: + for ot in old['t_array']: + if ot in outputs['t_array'][cropnew:]: + raise RuntimeError( + "For some reason t_array has some duplicated entries in the new and old file.") # Create single file for output and write data - with DataFile(fullpath,create=True,write=True,format=format, **kwargs) as f: + with DataFile(fullpath, create=True, write=True, format=format, **kwargs) as f: for varname in outputvars: if not quiet: print(varname) var = outputs[varname] if append: - dims=old.dimensions(varname) + dims = outputs.dimensions[varname] if 't' in dims: - varold=old[varname] - var=BoutArray(numpy.append(varold,var,axis=0),var.attributes) + var = var[cropnew:, ...] + varold = old[varname] + var = BoutArray(numpy.append( + varold, var, axis=0), var.attributes) if singleprecision: if not isinstance(var, int): @@ -133,15 +147,15 @@ def squashoutput(datadir=".", outputname="BOUT.dmp.nc", format="NETCDF4", tind=N f.write(varname, var) # Write changes, free memory f.sync() - var=None + var = None gc.collect() if delete: if append: os.remove(oldfile) - for f in glob.glob(datadir+"/BOUT.dmp.*.??"): + for f in glob.glob(datadir + "/BOUT.dmp.*.??"): if not quiet: - print("Deleting",f) + print("Deleting", f) os.remove(f) if append: os.rmdir(datadir)