diff --git a/.gitignore b/.gitignore index 934da1c0de..791afe8417 100644 --- a/.gitignore +++ b/.gitignore @@ -87,3 +87,7 @@ coverage/ /BOUT++-v*.tar.xz /CMakeCache.txt /CMakeFiles/cmake.check_cache + + +externalpackages/cpptrace/ +.vscode/ diff --git a/bin/bout-add-mod-path b/bin/bout-add-mod-path index 9faf1be3de..ee7ea66d25 100755 --- a/bin/bout-add-mod-path +++ b/bin/bout-add-mod-path @@ -104,8 +104,7 @@ def create_mod(modulepath, name, top, build): else: prereq = "" with open(filename, "w") as f: - f.write( - f"""#%Module 1.0 + f.write(f"""#%Module 1.0 # # BOUT++ module for use with 'environment-modules' package # Created by bout-add-mod-path v0.9 @@ -119,17 +118,14 @@ setenv BOUT_TOP {top} prepend-path PATH {top}/bin prepend-path PYTHONPATH {top}/tools/pylib prepend-path IDL_PATH +{top}/tools/idllib:'' -""" - ) +""") if build != top: - f.write( - f"""#%Module 1.0 + f.write(f"""#%Module 1.0 setenv BOUT_BUILD {build} prepend-path PATH {build}/bin prepend-path LD_LIBRARY_PATH {build}/lib prepend-path PYTHONPATH {build}/tools/pylib -""" - ) +""") print(f"created `{filename}`") diff --git a/bin/bout-pylib-cmd-to-bin b/bin/bout-pylib-cmd-to-bin index 8f88a5dbf4..8d6dbed7ae 100755 --- a/bin/bout-pylib-cmd-to-bin +++ b/bin/bout-pylib-cmd-to-bin @@ -126,8 +126,7 @@ def createwrapper(mod, func_name, func, name): out += end f.write(out) - fprint( - """#!/usr/bin/env python3 + fprint("""#!/usr/bin/env python3 # PYTHON_ARGCOMPLETE_OK import argparse @@ -136,8 +135,7 @@ try: except ImportError: argcomplete=None -""" - ) +""") doc = True para = False ret = False @@ -183,19 +181,16 @@ except ImportError: arg_help[curarg].append(esc(blas)) # Print functions that are needed if "str_to_slice" in arg_type.values(): - fprint( - """ + fprint(""" def str_to_slice(sstr): args=[] for s in sstr.split(','): args.append(int(s)) print(args) return slice(*args) -""" - ) +""") if "str_to_bool" in arg_type.values(): - fprint( - """ + fprint(""" def str_to_bool(sstr): low=sstr.lower() # no or false @@ -206,8 +201,7 @@ def str_to_bool(sstr): return True else: raise ArgumentTypeError("Cannot parse %s to bool type"%sstr) -""" - ) +""") # Create the parser fprint("parser = argparse.ArgumentParser(%s)" % (esc(docs))) spec = inspect.signature(func) @@ -247,24 +241,19 @@ def str_to_bool(sstr): pre = "\n " fprint(")") - fprint( - """ + fprint(""" if argcomplete: argcomplete.autocomplete(parser) -# late import for faster auto-complete""" - ) +# late import for faster auto-complete""") fprint("from %s import %s" % (mod, func_name)) - fprint( - """ + fprint(""" args = parser.parse_args() # Call the function %s, using command line arguments %s(**args.__dict__) -""" - % (func_name, func_name) - ) +""" % (func_name, func_name)) # alternative, but I think 0o755 is easier to read # import stat # os.chmod(filename,stat.S_IRWXU|stat.S_IRGRP|stat.S_IXGRP|stat.S_IROTH|stat.S_IXOTH) diff --git a/bin/bout-v5-factory-upgrader.py b/bin/bout-v5-factory-upgrader.py index 29fc07db30..ee24572a94 100755 --- a/bin/bout-v5-factory-upgrader.py +++ b/bin/bout-v5-factory-upgrader.py @@ -5,7 +5,6 @@ import difflib import re - # Dictionary of factory methods that may need updating factories = { "Interpolation": { @@ -62,9 +61,7 @@ def find_factory_calls(factory, source): \s*=\s* {factory_name}:: .*{create_method}.* - """.format( - **factory - ), + """.format(**factory), source, re.VERBOSE, ) @@ -75,9 +72,7 @@ def find_type_pointers(factory, source): r""" \b{type_name}\s*\*\s* # Type name and pointer ([\w_]+)\s*; # Variable name - """.format( - **factory - ), + """.format(**factory), source, re.VERBOSE, ) @@ -107,9 +102,7 @@ def fix_declarations(factory, variables, source): (.*?)(class\s*)? # optional "class" keyword \b({type_name})\s*\*\s* # Type-pointer ({variable_name})\s*; # Variable - """.format( - type_name=factory["type_name"], variable_name=variable - ), + """.format(type_name=factory["type_name"], variable_name=variable), r"\1std::unique_ptr<\3> \4{nullptr};", source, flags=re.VERBOSE, @@ -123,9 +116,7 @@ def fix_declarations(factory, variables, source): ({variable_name})\s* # Variable =\s* # Assignment from factory ({factory_name}::.*{create_method}.*); - """.format( - variable_name=variable, **factory - ), + """.format(variable_name=variable, **factory), r"\1auto \4 = \5;", source, flags=re.VERBOSE, @@ -139,9 +130,7 @@ def fix_declarations(factory, variables, source): ({variable_name})\s* # Variable =\s* # Assignment (0|nullptr|NULL); - """.format( - variable_name=variable, **factory - ), + """.format(variable_name=variable, **factory), r"\1std::unique_ptr<\2> \3{nullptr};", source, flags=re.VERBOSE, diff --git a/bin/bout-v5-format-upgrader.py b/bin/bout-v5-format-upgrader.py index 7c7d13ac7f..a534ca240a 100755 --- a/bin/bout-v5-format-upgrader.py +++ b/bin/bout-v5-format-upgrader.py @@ -5,7 +5,6 @@ import difflib import re - format_replacements = { "c": "c", "d": "d", diff --git a/bin/bout-v5-header-upgrader.py b/bin/bout-v5-header-upgrader.py index 49a8fbcbe4..77794ab920 100755 --- a/bin/bout-v5-header-upgrader.py +++ b/bin/bout-v5-header-upgrader.py @@ -9,7 +9,6 @@ from typing import List from subprocess import run - header_shim_sentinel = "// BOUT++ header shim" header_warning = f"""\ @@ -122,8 +121,7 @@ def create_patch(filename, original, modified): if __name__ == "__main__": parser = argparse.ArgumentParser( formatter_class=argparse.RawDescriptionHelpFormatter, - description=textwrap.dedent( - """\ + description=textwrap.dedent("""\ Fix deprecated header locations for BOUT++ v4 -> v5 All BOUT++ headers are now under ``include/bout`` and @@ -142,8 +140,7 @@ def create_patch(filename, original, modified): If you have staged changes, this tool will not work, so to avoid committing undesired or unrelated changes. - """ - ), + """), ) parser.add_argument( diff --git a/bin/bout-v5-input-file-upgrader.py b/bin/bout-v5-input-file-upgrader.py index e2940ff58a..ea979005d5 100755 --- a/bin/bout-v5-input-file-upgrader.py +++ b/bin/bout-v5-input-file-upgrader.py @@ -271,8 +271,7 @@ def possibly_apply_patch(patch, options_file, quiet=False, force=False): if __name__ == "__main__": parser = argparse.ArgumentParser( formatter_class=argparse.RawDescriptionHelpFormatter, - description=textwrap.dedent( - """\ + description=textwrap.dedent("""\ Fix input files for BOUT++ v5+ Please note that this will only fix input options in sections with @@ -300,8 +299,7 @@ def possibly_apply_patch(patch, options_file, quiet=False, force=False): Files that change in this way will have the "canonicalisation" patch presented first. If you choose not to apply this patch, the "upgrade - fixer" patch will still include it.""" - ), + fixer" patch will still include it."""), ) parser.add_argument("files", action="store", nargs="+", help="Input files") diff --git a/bin/bout-v5-macro-upgrader.py b/bin/bout-v5-macro-upgrader.py index 11b4926255..d644fed9e8 100755 --- a/bin/bout-v5-macro-upgrader.py +++ b/bin/bout-v5-macro-upgrader.py @@ -342,8 +342,7 @@ def create_patch(filename, original, modified): if __name__ == "__main__": parser = argparse.ArgumentParser( formatter_class=argparse.RawDescriptionHelpFormatter, - description=textwrap.dedent( - """\ + description=textwrap.dedent("""\ Fix macro defines for BOUT++ v4 -> v5 Please note that this is only slightly better than dumb text replacement. It @@ -359,8 +358,7 @@ def create_patch(filename, original, modified): still replace them in strings or comments. Please check the diff output carefully! - """ - ), + """), ) parser.add_argument("files", action="store", nargs="+", help="Input files") diff --git a/bin/bout-v5-physics-model-upgrader.py b/bin/bout-v5-physics-model-upgrader.py index 26fb8ef6e0..260d1d59ee 100755 --- a/bin/bout-v5-physics-model-upgrader.py +++ b/bin/bout-v5-physics-model-upgrader.py @@ -8,7 +8,6 @@ import textwrap import warnings - PHYSICS_MODEL_INCLUDE = '#include "bout/physicsmodel.hxx"' PHYSICS_MODEL_SKELETON = """ @@ -213,13 +212,11 @@ def fix_bout_constrain(source, error_on_warning): "\n ".join(["{}:{}".format(i, source_lines[i]) for i in line_range]) ) - message = textwrap.dedent( - """\ + message = textwrap.dedent("""\ Some uses of `bout_constrain` remain, but we could not automatically convert them to use `Solver::constraint`. Please fix them before continuing: - """ - ) + """) message += " " + "\n ".join(lines_context) if error_on_warning: @@ -389,8 +386,7 @@ def create_patch(filename, original, modified): if __name__ == "__main__": parser = argparse.ArgumentParser( formatter_class=argparse.RawDescriptionHelpFormatter, - description=textwrap.dedent( - """\ + description=textwrap.dedent("""\ Upgrade legacy physics models to use the PhysicsModel class This will do the bare minimum required to compile, and @@ -403,8 +399,7 @@ def create_patch(filename, original, modified): By default, this will use the file name stripped of file extensions as the name of the new class. Use '--name=' to give a different name. - """ - ), + """), ) parser.add_argument("files", action="store", nargs="+", help="Files to fix") diff --git a/bin/bout-v5-xzinterpolation-upgrader.py b/bin/bout-v5-xzinterpolation-upgrader.py index 37c79e0de8..e70c3c54ae 100755 --- a/bin/bout-v5-xzinterpolation-upgrader.py +++ b/bin/bout-v5-xzinterpolation-upgrader.py @@ -54,9 +54,7 @@ def fix_header_includes(old_header, new_header, source): (<|") ({header}) # Header name (>|") - """.format( - header=old_header - ), + """.format(header=old_header), r"\1\2{header}\4".format(header=new_header), source, flags=re.VERBOSE, @@ -67,9 +65,7 @@ def fix_interpolations(old_interpolation, new_interpolation, source): return re.sub( r""" \b{}\b - """.format( - old_interpolation - ), + """.format(old_interpolation), r"{}".format(new_interpolation), source, flags=re.VERBOSE, @@ -120,9 +116,7 @@ def fix_factories(old_factory, new_factory, source): return re.sub( r""" \b{}\b - """.format( - old_factory - ), + """.format(old_factory), r"{}".format(new_factory), source, flags=re.VERBOSE, diff --git a/bin/update_version_number_in_files.py b/bin/update_version_number_in_files.py index ec9a31bc32..3e9758a8c1 100755 --- a/bin/update_version_number_in_files.py +++ b/bin/update_version_number_in_files.py @@ -158,8 +158,7 @@ def create_patch(filename, original, modified): if __name__ == "__main__": parser = argparse.ArgumentParser( formatter_class=argparse.RawDescriptionHelpFormatter, - description=textwrap.dedent( - """\ + description=textwrap.dedent("""\ Update the software version number to the specified version, to be given in the form major.minor.patch, e.g. 5.10.3 @@ -172,8 +171,7 @@ def create_patch(filename, original, modified): the 'minor' version number of the provided version will be incremented by 1, e.g. 5.10.3 -> 5.11.3 - """ - ), + """), ) parser.add_argument( diff --git a/externalpackages/boutdata b/externalpackages/boutdata index 962d5a58b8..0e0626716f 160000 --- a/externalpackages/boutdata +++ b/externalpackages/boutdata @@ -1 +1 @@ -Subproject commit 962d5a58b8e6133a2977c3e44f1d4751a11d2b1c +Subproject commit 0e0626716f7fec01d318c427a96e596055540efe diff --git a/externalpackages/googletest b/externalpackages/googletest index 9156d4caac..a407966592 160000 --- a/externalpackages/googletest +++ b/externalpackages/googletest @@ -1 +1 @@ -Subproject commit 9156d4caac880b513264ecbe0aa4746a3fead3d7 +Subproject commit a407966592e8aab5c7b0757867a2f0585175fe80 diff --git a/include/bout/boutexception.hxx b/include/bout/boutexception.hxx index 238c88654c..324abcc0e6 100644 --- a/include/bout/boutexception.hxx +++ b/include/bout/boutexception.hxx @@ -18,9 +18,9 @@ class BoutException : public std::exception { public: BoutException(std::string msg); - template - BoutException(const S& format, const Args&... args) - : BoutException(fmt::format(format, args...)) {} + template + BoutException(fmt::format_string format, Args&&... args) + : BoutException(fmt::vformat(format, fmt::make_format_args(args...))) {} ~BoutException() override; @@ -45,17 +45,17 @@ private: class BoutRhsFail : public BoutException { public: BoutRhsFail(std::string message) : BoutException(std::move(message)) {} - template - BoutRhsFail(const S& format, const Args&... args) - : BoutRhsFail(fmt::format(format, args...)) {} + template + BoutRhsFail(fmt::format_string format, Args&&... args) + : BoutRhsFail(fmt::vformat(format, fmt::make_format_args(args...))) {} }; class BoutIterationFail : public BoutException { public: BoutIterationFail(std::string message) : BoutException(std::move(message)) {} - template - BoutIterationFail(const S& format, const Args&... args) - : BoutIterationFail(fmt::format(format, args...)) {} + template + BoutIterationFail(fmt::format_string format, Args&&... args) + : BoutIterationFail(fmt::vformat(format, fmt::make_format_args(args...))) {} }; #endif diff --git a/include/bout/field2d.hxx b/include/bout/field2d.hxx index 92658f1bbf..6bd2e49d03 100644 --- a/include/bout/field2d.hxx +++ b/include/bout/field2d.hxx @@ -254,7 +254,7 @@ public: // FieldData virtual functions - bool is3D() const override { return false; } + FieldType field_type() const override { return FieldType::field2d; } #if CHECK > 0 void doneComms() override { bndry_xin = bndry_xout = bndry_yup = bndry_ydown = true; } diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index a75e38df36..bd13ea3cfa 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -455,7 +455,7 @@ public: ///@} // FieldData virtual functions - bool is3D() const override { return true; } + FieldType field_type() const override { return FieldType::field3d; } #if CHECK > 0 void doneComms() override { bndry_xin = bndry_xout = bndry_yup = bndry_ydown = true; } diff --git a/include/bout/field_data.hxx b/include/bout/field_data.hxx index 185dcabf2d..8ddfbc8578 100644 --- a/include/bout/field_data.hxx +++ b/include/bout/field_data.hxx @@ -33,6 +33,7 @@ class FieldData; #include "bout/bout_types.hxx" #include "bout/unused.hxx" +#include #include #include #include @@ -71,9 +72,18 @@ public: /// Get variable location virtual CELL_LOC getLocation() const; + /// Enum to distinguish the different kinds of Fields + enum class FieldType : std::uint8_t { field3d, field2d, fieldperp }; + /// Is this an instance of `Field3D`, `Field2D`, or `FieldPerp`? + virtual FieldType field_type() const = 0; + // Defines interface which must be implemented /// True if variable is 3D - virtual bool is3D() const = 0; + [[deprecated("Use `field_type()` instead")]] + bool is3D() const { + return field_type() == FieldType::field3d; + } + /// Number of BoutReals in one element virtual int elementSize() const { return 1; } diff --git a/include/bout/fieldgroup.hxx b/include/bout/fieldgroup.hxx index 184766c6b8..a75f12db98 100644 --- a/include/bout/fieldgroup.hxx +++ b/include/bout/fieldgroup.hxx @@ -13,10 +13,10 @@ /// Group together fields for easier communication /// -/// Note: The FieldData class is used as a base class, -/// which is inherited by Field2D, Field3D, Vector2D and Vector3D -/// however Vector2D and Vector3D are stored by reference to their -/// components (x,y,z) as Field2D or Field3D objects. +/// Note: The `Field` class is used as a base class, +/// which is inherited by `Field2D`, `Field3D`, `FieldPerp`; +/// however `Vector2D` and `Vector3D` are stored by reference to their +/// components ``(x, y, z)`` as `Field2D` or `Field3D` objects. class FieldGroup { public: FieldGroup() = default; @@ -25,8 +25,8 @@ public: FieldGroup& operator=(const FieldGroup& other) = default; FieldGroup& operator=(FieldGroup&& other) = default; - /// Constructor with a single FieldData \p f - FieldGroup(FieldData& f) { fvec.push_back(&f); } + /// Constructor with a single Field \p f + FieldGroup(Field& f) { fvec.push_back(&f); } /// Constructor with a single Field3D \p f FieldGroup(Field3D& f) { @@ -56,7 +56,7 @@ public: } /// Variadic constructor. Allows an arbitrary number of - /// FieldData arguments + /// Field arguments /// /// The explicit keyword prevents FieldGroup being constructed with arbitrary /// types. In particular arguments to add() cannot be implicitly converted @@ -78,12 +78,12 @@ public: return *this; } - /// Add a FieldData \p f to the group. + /// Add a Field \p f to the group. /// /// A pointer to this field will be stored internally, /// so the lifetime of this variable should be longer /// than the lifetime of this group. - void add(FieldData& f) { fvec.push_back(&f); } + void add(Field& f) { fvec.push_back(&f); } // Add a 3D field \p f, which goes into both vectors. // @@ -126,7 +126,7 @@ public: /// treated as a special case. An arbitrary number of fields can be /// added. template - void add(FieldData& t, Ts&... ts) { + void add(Field& t, Ts&... ts) { add(t); // Add the first using functions above add(ts...); // Add the rest } @@ -165,16 +165,16 @@ public: } /// Iteration over all fields - using iterator = std::vector::iterator; + using iterator = std::vector::iterator; iterator begin() { return fvec.begin(); } iterator end() { return fvec.end(); } /// Const iteration over all fields - using const_iterator = std::vector::const_iterator; + using const_iterator = std::vector::const_iterator; const_iterator begin() const { return fvec.begin(); } const_iterator end() const { return fvec.end(); } - const std::vector& get() const { return fvec; } + const std::vector& get() const { return fvec; } /// Iteration over 3D fields const std::vector& field3d() const { return f3vec; } @@ -183,7 +183,7 @@ public: void makeUnique(); private: - std::vector fvec; // Vector of fields + std::vector fvec; // Vector of fields std::vector f3vec; // Vector of 3D fields }; diff --git a/include/bout/fieldperp.hxx b/include/bout/fieldperp.hxx index 6995308dbe..37e13e2363 100644 --- a/include/bout/fieldperp.hxx +++ b/include/bout/fieldperp.hxx @@ -286,7 +286,7 @@ public: */ int getNz() const override { return nz; }; - bool is3D() const override { return false; } + FieldType field_type() const override { return FieldType::fieldperp; } friend void swap(FieldPerp& first, FieldPerp& second) noexcept; diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index 02e2a23905..28b56bc48e 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -344,10 +344,15 @@ public: // non-local communications - virtual int getNXPE() = 0; ///< The number of processors in the X direction - virtual int getNYPE() = 0; ///< The number of processors in the Y direction - virtual int getXProcIndex() = 0; ///< This processor's index in X direction - virtual int getYProcIndex() = 0; ///< This processor's index in Y direction + virtual int getNXPE() const = 0; ///< The number of processors in the X direction + virtual int getNYPE() const = 0; ///< The number of processors in the Y direction + virtual int getNZPE() const = 0; ///< The number of processors in the Z direction + virtual int getXProcIndex() const = 0; ///< This processor's index in X direction + virtual int getYProcIndex() const = 0; ///< This processor's index in Y direction + virtual int getZProcIndex() const = 0; ///< This processor's index in Z direction + + /// The rank of the processor at (\p X, \p Y, \p Z) in the processor grid + virtual int getProcIndex(int X, int Y, int Z) const = 0; // X communications virtual bool firstX() @@ -393,6 +398,7 @@ public: } ///< Return communicator containing all processors in X virtual MPI_Comm getXcomm(int jy) const = 0; ///< Return X communicator virtual MPI_Comm getYcomm(int jx) const = 0; ///< Return Y communicator + virtual MPI_Comm getXZcomm() const = 0; ///< Communicator in X-Z /// Return pointer to the mesh's MPI Wrapper object MpiWrapper& getMpi() { return *mpi; } @@ -464,11 +470,11 @@ public: /// Is there a boundary on the lower guard cells in Y /// on any processor along the X direction? - bool hasBndryLowerY(); + virtual bool hasBndryLowerY() const = 0; /// Is there a boundary on the upper guard cells in Y /// on any processor along the X direction? - bool hasBndryUpperY(); + virtual bool hasBndryUpperY() const = 0; // Boundary regions /// Return a vector containing all the boundary regions on this processor @@ -500,8 +506,10 @@ public: virtual BoutReal GlobalX(int jx) const = 0; ///< Continuous X index between 0 and 1 virtual BoutReal GlobalY(int jy) const = 0; ///< Continuous Y index (0 -> 1) + virtual BoutReal GlobalZ(int jz) const = 0; ///< Continuous Z index (0 -> 1) virtual BoutReal GlobalX(BoutReal jx) const = 0; ///< Continuous X index between 0 and 1 virtual BoutReal GlobalY(BoutReal jy) const = 0; ///< Continuous Y index (0 -> 1) + virtual BoutReal GlobalZ(BoutReal jz) const = 0; ///< Continuous Z index (0 -> 1) ////////////////////////////////////////////////////////// @@ -802,8 +810,8 @@ protected: const std::vector readInts(const std::string& name, int n); /// Calculates the size of a message for a given x and y range - int msg_len(const std::vector& var_list, int xge, int xlt, int yge, - int ylt); + int msg_len(const std::vector& var_list, int xge, int xlt, int yge, + int ylt) const; /// Initialise derivatives void derivs_init(Options* options); diff --git a/include/bout/sys/gettext.hxx b/include/bout/sys/gettext.hxx index d67becd54d..d62cd66882 100644 --- a/include/bout/sys/gettext.hxx +++ b/include/bout/sys/gettext.hxx @@ -7,15 +7,47 @@ #if BOUT_HAS_GETTEXT -#include +#include // IWYU pragma: keep + #include #define GETTEXT_PACKAGE "libbout" -#define _(string) dgettext(GETTEXT_PACKAGE, string) +// If we have C++23, we can get fmt to do compile-time checks of our format +// strings, _and_ have gettext do runtime replacement +#if __cpp_if_consteval >= 202106L +constexpr const char* dgettext_wrap(const char* __domainname, const char* __msgid) __THROW + __attribute_format_arg__(2); + +constexpr const char* dgettext_wrap(const char* __domainname, const char* __msgid) { + if consteval { + return __msgid; + } + return dgettext(__domainname, __msgid); +} + +/// Gettext i18n macro for text containing fmt format specifiers +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) +#define _f(string) dgettext_wrap(GETTEXT_PACKAGE, string) #else +// We're pre-C++23, so all our i18n text must be fmt runtime formats +#include "fmt/base.h" // IWYU pragma: keep + +/// Gettext i18n macro for text containing fmt format specifiers +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) +#define _f(string) fmt::runtime(dgettext(GETTEXT_PACKAGE, string)) +#endif + +/// Gettext i18n macro for plain text that doesn't need formatting +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) +#define _(string) dgettext(GETTEXT_PACKAGE, string) + +#else // BOUT_HAS_GETTEXT +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) +#define _f(string) string +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) #define _(string) string #endif // BOUT_HAS_GETTEXT diff --git a/include/bout/vector2d.hxx b/include/bout/vector2d.hxx index bdc375e698..0317b210b6 100644 --- a/include/bout/vector2d.hxx +++ b/include/bout/vector2d.hxx @@ -142,7 +142,7 @@ public: CELL_LOC getLocation() const override; // FieldData virtual functions - bool is3D() const override { return false; } + FieldType field_type() const override { return FieldType::field2d; } int elementSize() const override { return 3; } /// Apply boundary condition to all fields diff --git a/include/bout/vector3d.hxx b/include/bout/vector3d.hxx index 0c71dcffa5..e230295e93 100644 --- a/include/bout/vector3d.hxx +++ b/include/bout/vector3d.hxx @@ -172,7 +172,7 @@ public: CELL_LOC getLocation() const override; // FieldData virtual functions - bool is3D() const override { return true; } + FieldType field_type() const override { return FieldType::field3d; } int elementSize() const override { return 3; } void applyBoundary(bool init = false) override; diff --git a/src/field/gen_fieldops.py b/src/field/gen_fieldops.py index 29631ff7aa..3e07d6fec4 100755 --- a/src/field/gen_fieldops.py +++ b/src/field/gen_fieldops.py @@ -11,7 +11,6 @@ """ - from __future__ import print_function from builtins import object diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index 574902ea7b..2cecc2354f 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -26,24 +26,45 @@ #include "boutmesh.hxx" +#include #include +#include #include #include #include #include #include #include +#include +#include +#include +#include +#include #include #include #include #include #include +#include +#include +#include #include #include +#include +#include + #include +#include +#include #include +#include +#include +#include #include +#include +#include +#include /// MPI type of BoutReal for communications #define PVEC_REAL_MPI_TYPE MPI_DOUBLE @@ -57,6 +78,8 @@ If you want the old setting, you have to specify mesh:symmetricGlobalY=false in << optionfile << "\n"; } OPTION(options, symmetricGlobalY, true); + OPTION(options, symmetricGlobalZ, false); // The default should be updated to true but + // this breaks backwards compatibility comm_x = MPI_COMM_NULL; comm_inner = MPI_COMM_NULL; @@ -84,6 +107,9 @@ BoutMesh::~BoutMesh() { if (comm_outer != MPI_COMM_NULL) { MPI_Comm_free(&comm_outer); } + if (comm_xz != MPI_COMM_NULL) { + MPI_Comm_free(&comm_xz); + } } BoutMesh::YDecompositionIndices @@ -152,82 +178,232 @@ void BoutMesh::setXDecompositionIndices(const XDecompositionIndices& indices) { ixseps2 = indices.ixseps2; } -namespace bout { -CheckMeshResult checkBoutMeshYDecomposition(int num_y_processors, int ny, - int num_y_guards, int jyseps1_1, - int jyseps2_1, int jyseps1_2, int jyseps2_2, - int ny_inner) { - const int num_local_y_points = ny / num_y_processors; +MeshTopology BoutMesh::getMeshTopology(int jyseps1_1_, int jyseps2_1_, //Returns MeshTopology that is type enum + int jyseps1_2_, int jyseps2_2_, + int ny_inner_, int ixseps1_, + int ixseps2_) { + + // Set member variables + jyseps1_1 = jyseps1_1_; + jyseps2_1 = jyseps2_1_; + jyseps1_2 = jyseps1_2_; + jyseps2_2 = jyseps2_2_; + ny_inner = ny_inner_; + ixseps1 = ixseps1_; + ixseps2 = ixseps2_; - // Check size of Y mesh if we've got multiple processors in Y - if (num_local_y_points < num_y_guards and num_y_processors != 1) { - return {false, - fmt::format(_("\t -> ny/NYPE ({:d}/{:d} = {:d}) must be >= MYG ({:d})\n"), ny, - num_y_processors, num_local_y_points, num_y_guards)}; - } - // Check branch cuts - if ((jyseps1_1 + 1) % num_local_y_points != 0) { - return {false, fmt::format(_("\t -> Leg region jyseps1_1+1 ({:d}) must be a " - "multiple of MYSUB ({:d})\n"), - jyseps1_1 + 1, num_local_y_points)}; - } + if (jyseps1_1 < 0 and jyseps2_2 >= ny - 1) { + return MeshTopology::closed_field_line; //write whole thing + } else if (jyseps2_1 == jyseps1_2) { + return MeshTopology::single_null; + } else if (ixseps1 == ixseps2) { + return MeshTopology::connected_double_null; + } + else{ + return MeshTopology::unconnected_double_null; + } +} - if (jyseps2_1 != jyseps1_2) { - // Double Null - if ((jyseps2_1 - jyseps1_1) % num_local_y_points != 0) { - return { - false, - fmt::format(_("\t -> Core region jyseps2_1-jyseps1_1 ({:d}-{:d} = {:d}) must " - "be a multiple of MYSUB ({:d})\n"), - jyseps2_1, jyseps1_1, jyseps2_1 - jyseps1_1, num_local_y_points)}; - } +namespace bout { + CheckMeshResult checkBoutMeshYDecomposition( + int num_y_processors, int ny, + int num_y_guards, + int jyseps1_1, int jyseps2_1, + int jyseps1_2, int jyseps2_2, + int ny_inner) { + + // Preserve legacy behaviour (Single / Double Null) + return checkBoutMeshYDecomposition( + num_y_processors, ny, num_y_guards, + jyseps1_1, jyseps2_1, + jyseps1_2, jyseps2_2, + ny_inner, + MeshTopology::unconnected_double_null); + } + + } // namespace bout + + namespace bout { + CheckMeshResult checkBoutMeshYDecomposition(int num_y_processors, int ny, + int num_y_guards, int jyseps1_1, + int jyseps2_1, int jyseps1_2, int jyseps2_2, + int ny_inner, MeshTopology mesh_topology) { + + const int num_local_y_points = ny / num_y_processors; + + // Check size of Y mesh if we've got multiple processors in Y + if (num_local_y_points < num_y_guards and num_y_processors != 1) { + return {false, + fmt::format(_f("\t -> ny/NYPE ({:d}/{:d} = {:d}) must be >= MYG ({:d})\n"), ny, + num_y_processors, num_local_y_points, num_y_guards)}; + } + // Check branch cuts + if ((jyseps1_1 + 1) % num_local_y_points != 0) { + return {false, fmt::format(_f("\t -> Leg region jyseps1_1+1 ({:d}) must be a " + "multiple of MYSUB ({:d})\n"), + jyseps1_1 + 1, num_local_y_points)}; + } + + if (mesh_topology == MeshTopology::unconnected_double_null || mesh_topology == MeshTopology::connected_double_null){ + if ((jyseps2_1 - jyseps1_1) % num_local_y_points != 0) { + return { + false, + fmt::format(_f("\t -> Core region jyseps2_1-jyseps1_1 ({:d}-{:d} = {:d}) must " + "be a multiple of MYSUB ({:d})\n"), + jyseps2_1, jyseps1_1, jyseps2_1 - jyseps1_1, num_local_y_points)}; + } - if ((jyseps2_2 - jyseps1_2) % num_local_y_points != 0) { - return { - false, - fmt::format(_("\t -> Core region jyseps2_2-jyseps1_2 ({:d}-{:d} = {:d}) must " - "be a multiple of MYSUB ({:d})\n"), - jyseps2_2, jyseps1_2, jyseps2_2 - jyseps1_2, num_local_y_points)}; - } + if ((jyseps2_2 - jyseps1_2) % num_local_y_points != 0) { + return { + false, + fmt::format(_f("\t -> Core region jyseps2_2-jyseps1_2 ({:d}-{:d} = {:d}) must " + "be a multiple of MYSUB ({:d})\n"), + jyseps2_2, jyseps1_2, jyseps2_2 - jyseps1_2, num_local_y_points)}; + } - // Check upper legs - if ((ny_inner - jyseps2_1 - 1) % num_local_y_points != 0) { - return { - false, - fmt::format(_("\t -> leg region ny_inner-jyseps2_1-1 ({:d}-{:d}-1 = {:d}) must " - "be a multiple of MYSUB ({:d})\n"), - ny_inner, jyseps2_1, ny_inner - jyseps2_1 - 1, num_local_y_points)}; - } - if ((jyseps1_2 - ny_inner + 1) % num_local_y_points != 0) { - return { - false, - fmt::format(_("\t -> leg region jyseps1_2-ny_inner+1 ({:d}-{:d}+1 = {:d}) must " - "be a multiple of MYSUB ({:d})\n"), - jyseps1_2, ny_inner, jyseps1_2 - ny_inner + 1, num_local_y_points)}; - } - } else { - // Single Null + // Check upper legs + if ((ny_inner - jyseps2_1 - 1) % num_local_y_points != 0) { + return { + false, + fmt::format(_f("\t -> leg region ny_inner-jyseps2_1-1 ({:d}-{:d}-1 = {:d}) must " + "be a multiple of MYSUB ({:d})\n"), + ny_inner, jyseps2_1, ny_inner - jyseps2_1 - 1, num_local_y_points)}; + } + if ((jyseps1_2 - ny_inner + 1) % num_local_y_points != 0) { + return { + false, + fmt::format(_f("\t -> leg region jyseps1_2-ny_inner+1 ({:d}-{:d}+1 = {:d}) must " + "be a multiple of MYSUB ({:d})\n"), + jyseps1_2, ny_inner, jyseps1_2 - ny_inner + 1, num_local_y_points)}; + } + } else if ((mesh_topology == MeshTopology::single_null) || (mesh_topology == MeshTopology::closed_field_line)){ + // Single Null or connected Double Null if ((jyseps2_2 - jyseps1_1) % num_local_y_points != 0) { return { false, - fmt::format(_("\t -> Core region jyseps2_2-jyseps1_1 ({:d}-{:d} = {:d}) must " - "be a multiple of MYSUB ({:d})\n"), + fmt::format(_f("\t -> Core region jyseps2_2-jyseps1_1 ({:d}-{:d} = {:d}) must " + "be a multiple of MYSUB ({:d})\n"), jyseps2_2, jyseps1_1, jyseps2_2 - jyseps1_1, num_local_y_points)}; } } - if ((ny - jyseps2_2 - 1) % num_local_y_points != 0) { + if ((ny - 1 - jyseps2_2) % num_local_y_points != 0) { return {false, fmt::format( - _("\t -> leg region ny-jyseps2_2-1 ({:d}-{:d}-1 = {:d}) must be a " - "multiple of MYSUB ({:d})\n"), - ny, jyseps2_2, ny - jyseps2_2 - 1, num_local_y_points)}; + _f("\t -> leg region ny-jyseps2_2-1 ({:d}-{:d}-1 = {:d}) must be a " + "multiple of MYSUB ({:d})\n"), + ny, jyseps2_2, ny - 1 - jyseps2_2, num_local_y_points)}; } return {true, ""}; + } + + + CheckMeshResult findValidProcessorNum(int ny, int nx, + int NPES, + int NYPE, + int NXPE) { + int best_nxpe = 0; + int best_nype = 0; + int best_npes = 0; + + for (int possible_nxpe = NXPE; possible_nxpe <= NPES; ++possible_nxpe) { + if (nx % possible_nxpe != 0) continue; + + for (int possible_nype = NYPE; + possible_nype <= NPES / possible_nxpe; + ++possible_nype) { + + if (possible_nype == 1 && possible_nxpe == 1 && NPES > 1) continue; // Skip single processor unless NPES=1 + if (ny % possible_nype != 0) continue; + if (NPES % (possible_nxpe * possible_nype) != 0) continue; + int possible_npes = possible_nxpe * possible_nype; + + if (possible_npes > best_npes) { + best_npes = possible_npes; + best_nxpe = possible_nxpe; + best_nype = possible_nype; + } + } + } + + if (best_npes > 0) { + return { + true, + fmt::format( + "\t -> Best processor decomposition found: " + "NPES={:d}, NXPE={:d}, and NYPE={:d}.", + best_npes, best_nxpe, best_nype) + }; + } + + return { + false, + fmt::format( + "\t -> No valid processor decomposition found for nx={:d}, ny={:d} with the number of given processors NPES = {:d}. " + "Try changing the number of points.", + nx, ny, NPES) + }; + } + + + CheckMeshResult findValidYDecomposition(int ny, + int num_y_processors, + int num_y_guards, + int jyseps1_1_start, + int jyseps2_1_start, + int jyseps1_2_start, + int jyseps2_2_start, + int ny_inner_start, + MeshTopology mesh_topology) { + + if (ny % num_y_processors != 0) { + return {false, fmt::format( + "\t ->ny ({:d}) must be divisible by NYPE ({:d}). Try changing the number of points or processors in Y.",ny, num_y_processors)}; + } + + for (int jyseps1_1 = jyseps1_1_start; jyseps1_1 < ny; ++jyseps1_1) { + + for (int jyseps2_1 = jyseps2_1_start; + jyseps2_1 < ny; ++jyseps2_1) { + + for (int jyseps1_2 = jyseps1_2_start; + jyseps1_2 < ny; ++jyseps1_2) { + + for (int ny_inner = ny_inner_start; + ny_inner < ny; ++ny_inner) { + + for (int jyseps2_2 = jyseps2_2_start; + jyseps2_2 < ny; ++jyseps2_2) { + + if (mesh_topology == MeshTopology::unconnected_double_null || mesh_topology == MeshTopology::connected_double_null){ + if (not (jyseps1_1 < jyseps2_1 && + jyseps2_1 < ny_inner && + ny_inner < jyseps1_2 && + jyseps1_2 < jyseps2_2)){ + continue; + }} + + auto result = bout::checkBoutMeshYDecomposition(num_y_processors, ny, num_y_guards, jyseps1_1, jyseps2_1, + jyseps1_2, jyseps2_2, ny_inner, mesh_topology); + + if (result.success) { + return {true, fmt::format( + "\t -> A valid decomposition in Y close to the one given in the grid would be: " + "jyseps1_1={:d}, jyseps2_1={:d}, jyseps1_2={:d}, jyseps2_2={:d}, ny_inner={:d}\n", + jyseps1_1, jyseps2_1, jyseps1_2, jyseps2_2, ny_inner)}; + } + } + } + } + } + } + + return {false, fmt::format( + "\t -> No valid Y decomposition found for ny = {:d} and NYPE = {:d}. Try changing the number of points or processors in Y.",ny, num_y_processors)}; + } } -} // namespace bout void BoutMesh::chooseProcessorSplit(Options& options) { // Possible issues: @@ -242,10 +418,16 @@ void BoutMesh::chooseProcessorSplit(Options& options) { .withDefault(1); if ((NPES % NXPE) != 0) { throw BoutException( - _("Number of processors ({:d}) not divisible by NPs in x direction ({:d})\n"), + _f("Number of processors ({:d}) not divisible by NPs in x direction ({:d})\n"), NPES, NXPE); } + if (nx % NXPE != 0) { + throw BoutException( + _f("Number of x points ({:d}) not divisible by NPs in x direction ({:d})\n"), nx, + NXPE); + } + NYPE = NPES / NXPE; } else { // NXPE not set, but NYPE is @@ -255,17 +437,26 @@ void BoutMesh::chooseProcessorSplit(Options& options) { .withDefault(1); if ((NPES % NYPE) != 0) { throw BoutException( - _("Number of processors ({:d}) not divisible by NPs in y direction ({:d})\n"), + _f("Number of processors ({:d}) not divisible by NPs in y direction ({:d})\n"), NPES, NYPE); } + if (ny % NYPE != 0) { + throw BoutException( + _f("Number of y points ({:d}) not divisible by NPs in y direction ({:d})\n"), nx, + NXPE); + } + NXPE = NPES / NYPE; } auto result = bout::checkBoutMeshYDecomposition(NYPE, ny, MYG, jyseps1_1, jyseps2_1, - jyseps1_2, jyseps2_2, ny_inner); + jyseps1_2, jyseps2_2, ny_inner, mesh_topology); if (not result.success) { + auto valid_y_decompostion= bout::findValidYDecomposition(ny, NYPE, MYG, jyseps1_1, jyseps2_1, + jyseps1_2, jyseps2_2, ny_inner, mesh_topology); + output_info.write(valid_y_decompostion.reason); throw BoutException(result.reason); } } @@ -278,19 +469,19 @@ void BoutMesh::findProcessorSplit() { // Results in square domains const BoutReal ideal = sqrt(MX * NPES / static_cast(ny)); - output_info.write(_("Finding value for NXPE (ideal = {:f})\n"), ideal); + output_info.write(_f("Finding value for NXPE (ideal = {:f})\n"), ideal); for (int i = 1; i <= NPES; i++) { // Loop over all possibilities if ((NPES % i == 0) && // Processors divide equally (MX % i == 0) && // Mesh in X divides equally (ny % (NPES / i) == 0)) { // Mesh in Y divides equally - output_info.write(_("\tCandidate value: {:d}\n"), i); + output_info.write(_f("\tCandidate value: {:d}\n"), i); const int nyp = NPES / i; auto result = bout::checkBoutMeshYDecomposition(nyp, ny, MYG, jyseps1_1, jyseps2_1, - jyseps1_2, jyseps2_2, ny_inner); + jyseps1_2, jyseps2_2, ny_inner, mesh_topology); if (not result.success) { output_info.write(result.reason); @@ -306,21 +497,21 @@ void BoutMesh::findProcessorSplit() { } if (NXPE < 1) { - throw BoutException(_("Could not find a valid value for NXPE. Try a different " + throw BoutException(_f("Could not find a valid value for NXPE. Try a different " "number of processors.")); } NYPE = NPES / NXPE; - output_progress.write(_("\tDomain split (NXPE={:d}, NYPE={:d}) into domains " - "(localNx={:d}, localNy={:d})\n"), + output_progress.write(_f("\tDomain split (NXPE={:d}, NYPE={:d}) into domains " + "(localNx={:d}, localNy={:d})\n"), NXPE, NYPE, MX / NXPE, ny / NYPE); } void BoutMesh::setDerivedGridSizes() { // Check that nx is large enough if (nx <= 2 * MXG) { - throw BoutException(_("Error: nx must be greater than 2 times MXG (2 * {:d})"), MXG); + throw BoutException(_f("Error: nx must be greater than 2 times MXG (2 * {:d})"), MXG); } GlobalNx = nx; @@ -345,7 +536,9 @@ void BoutMesh::setDerivedGridSizes() { MX = nx - 2 * MXG; MXSUB = MX / NXPE; if ((MX % NXPE) != 0) { - throw BoutException(_("Cannot split {:d} X points equally between {:d} processors\n"), + auto valid_process_num = bout::findValidProcessorNum(ny, nx, NPES, NXPE); + output_info.write(valid_process_num.reason); + throw BoutException(_f("Cannot split {:d} X points equally between {:d} processors\n"), MX, NXPE); } @@ -353,8 +546,10 @@ void BoutMesh::setDerivedGridSizes() { MY = ny; MYSUB = MY / NYPE; if ((MY % NYPE) != 0) { + auto valid_process_num = bout::findValidProcessorNum(ny, nx, NPES, NYPE); + output_info.write(valid_process_num.reason); throw BoutException( - _("\tERROR: Cannot split {:d} Y points equally between {:d} processors\n"), MY, + _f("\tERROR: Cannot split {:d} Y points equally between {:d} processors\n"), MY, NYPE); } @@ -362,7 +557,7 @@ void BoutMesh::setDerivedGridSizes() { MZSUB = MZ / NZPE; if ((MZ % NZPE) != 0) { throw BoutException( - _("\tERROR: Cannot split {:d} Z points equally between {:d} processors\n"), MZ, + _f("\tERROR: Cannot split {:d} Z points equally between {:d} processors\n"), MZ, NZPE); } @@ -431,7 +626,6 @@ void BoutMesh::setDerivedGridSizes() { } int BoutMesh::load() { - TRACE("BoutMesh::load()"); output_progress << _("Loading mesh") << endl; @@ -464,8 +658,8 @@ int BoutMesh::load() { if (!is_pow2(nz)) { // Should be a power of 2 for efficient FFTs output_warn.write( - _("WARNING: Number of toroidal points should be 2^n for efficient " - "FFT performance -- consider changing MZ ({:d}) if using FFTs\n"), + _f("WARNING: Number of toroidal points should be 2^n for efficient " + "FFT performance -- consider changing MZ ({:d}) if using FFTs\n"), nz); } } else { @@ -484,8 +678,19 @@ int BoutMesh::load() { } ASSERT0(MXG >= 0); - if (Mesh::get(MYG, "MYG") != 0) { - MYG = options["MYG"].doc("Number of guard cells on each side in Y").withDefault(2); + const bool meshHasMyg = Mesh::get(MYG, "MYG") == 0; + if (!meshHasMyg) { + MYG = 2; + } + int meshMyg = MYG; + + if (options.isSet("MYG") or (!meshHasMyg)) { + MYG = options["MYG"].doc("Number of guard cells on each side in Y").withDefault(MYG); + } + if (meshHasMyg && MYG != meshMyg) { + output_warn.write(_f("Options changed the number of y-guard cells. Grid has {} but " + "option specified {}! Continuing with {}"), + meshMyg, MYG, MYG); } ASSERT0(MYG >= 0); @@ -508,6 +713,12 @@ int BoutMesh::load() { Mesh::get(jyseps2_2, "jyseps2_2", ny - 1); Mesh::get(ny_inner, "ny_inner", jyseps2_1); + mesh_topology = getMeshTopology(jyseps1_1, jyseps2_1, + jyseps1_2, jyseps2_2, + ny_inner, ixseps1, ixseps2); + output_info << _("Detected mesh topology = ") + << toString(mesh_topology) << std::endl; + // Check inputs setYDecompositionIndices(jyseps1_1, jyseps2_1, jyseps1_2, jyseps2_2, ny_inner); @@ -517,9 +728,12 @@ int BoutMesh::load() { findProcessorSplit(); } - // Get X and Y processor indices + // Get X, Y, Z processor indices PE_YIND = MYPE / NXPE; PE_XIND = MYPE % NXPE; + PE_ZIND = 0; + + ASSERT2(MYPE == getProcIndex(PE_XIND, PE_YIND, PE_ZIND)); // Set the other grid sizes from nx, ny, nz setDerivedGridSizes(); @@ -545,6 +759,7 @@ int BoutMesh::load() { ///////////////////// TOPOLOGY ////////////////////////// /// Call topology to set layout of grid + /// Renamed variable topology because of this function topology(); TwistShift = options["twistshift"] @@ -607,18 +822,65 @@ int BoutMesh::load() { // Add boundary regions addBoundaryRegions(); + // Set cached values + { + int mybndry = static_cast(!(iterateBndryLowerY().isDone())); + int allbndry = 0; + mpi->MPI_Allreduce(&mybndry, &allbndry, 1, MPI_INT, MPI_BOR, getXcomm(yend)); + has_boundary_lower_y = static_cast(allbndry); + } + { + int mybndry = static_cast(!(iterateBndryUpperY().isDone())); + int allbndry = 0; + mpi->MPI_Allreduce(&mybndry, &allbndry, 1, MPI_INT, MPI_BOR, getXcomm(ystart)); + has_boundary_upper_y = static_cast(allbndry); + } + // Initialize default coordinates getCoordinates(); output_info.write(_("\tdone\n")); return 0; +} // unchanged, is not topology dependent + +namespace { +auto make_XZ_communicator(const BoutMesh& mesh, MPI_Group group_world) -> MPI_Comm { + std::vector ranks; + + const int yp = mesh.getYProcIndex(); + + // All processors with the same Y index + for (int xp = 0; xp < mesh.getNXPE(); ++xp) { + for (int zp = 0; zp < mesh.getNZPE(); ++zp) { + ranks.push_back(mesh.getProcIndex(xp, yp, zp)); + } + } + MPI_Group group{}; + if (MPI_Group_incl(group_world, static_cast(ranks.size()), ranks.data(), &group) + != MPI_SUCCESS) { + throw BoutException("Could not create X-Z communication group for ranks {}", + fmt::join(ranks, ", ")); + } + + MPI_Comm comm_xz{}; + if (MPI_Comm_create(BoutComm::get(), group, &comm_xz) != MPI_SUCCESS) { + throw BoutException("Could not create X-Z communicator for yp={} (xind={}, yind={}, " + "zind={}) ranks={}", + yp, mesh.getXProcIndex(), mesh.getYProcIndex(), + mesh.getZProcIndex(), fmt::join(ranks, ", ")); + } + + return comm_xz; } +} // namespace void BoutMesh::createCommunicators() { MPI_Group group_world{}; MPI_Comm_group(BoutComm::get(), &group_world); // Get the entire group + comm_xz = make_XZ_communicator(*this, group_world); + ////////////////////////////////////////////////////// /// Communicator in X @@ -667,10 +929,10 @@ void BoutMesh::createCommunicators() { MPI_Group group_tmp2{}; proc[2] = NXPE; // Stride in processor rank - // Outer SOL regions - if (jyseps1_2 == jyseps2_1) { - // Single-null. All processors with same PE_XIND + if (mesh_topology == MeshTopology::single_null || mesh_topology == MeshTopology::closed_field_line) { + // Single-null and CFL + //All processors with same PE_XIND TRACE("Creating Outer SOL communicators for Single Null operation"); for (int i = 0; i < NXPE; i++) { @@ -694,8 +956,10 @@ void BoutMesh::createCommunicators() { } MPI_Group_free(&group); } - } else { + + } else if (mesh_topology == MeshTopology::connected_double_null || mesh_topology == MeshTopology::unconnected_double_null) { // Double null + // Difference with UCD and CDN comes from a secondary inner SOL region (ixseps1 != ixseps2) TRACE("Creating Outer SOL communicators for Double Null operation"); for (int i = 0; i < NXPE; i++) { @@ -728,121 +992,18 @@ void BoutMesh::createCommunicators() { MPI_Group_free(&group); } } - - for (int i = 0; i < NXPE; i++) { - // Lower PF region - - if ((jyseps1_1 >= 0) || (jyseps2_2 + 1 < ny)) { - // A lower PF region exists - TRACE("Creating lower PF communicators for xp={:d}", i); - - output_debug << "Creating lower PF communicators for xp = " << i << endl; - - if (jyseps1_1 >= 0) { - proc[0] = PROC_NUM(i, 0); - proc[1] = PROC_NUM(i, YPROC(jyseps1_1)); - - output_debug << "PF1 " << proc[0] << ", " << proc[1] << endl; - - MPI_Group_range_incl(group_world, 1, &proc, &group_tmp1); - } else { - group_tmp1 = MPI_GROUP_EMPTY; - } - - if (jyseps2_2 + 1 < ny) { - proc[0] = PROC_NUM(i, YPROC(jyseps2_2 + 1)); - proc[1] = PROC_NUM(i, NYPE - 1); - - output_debug << "PF2 " << proc[0] << ", " << proc[1] << endl; - - MPI_Group_range_incl(group_world, 1, &proc, &group_tmp2); - } else { - group_tmp2 = MPI_GROUP_EMPTY; - } - - MPI_Group_union(group_tmp1, group_tmp2, &group); - MPI_Comm_create(BoutComm::get(), group, &comm_tmp); - if (comm_tmp != MPI_COMM_NULL) { - comm_inner = comm_tmp; - if (ixseps_lower == ixseps_outer) { - // Between the separatrices is still in the PF region - - output_debug << "-> Inner and middle\n"; - - comm_middle = comm_inner; - } else { - - output_debug << "-> Outer and middle\n"; - - comm_middle = comm_outer; - } - } - - output_debug << "Freeing\n"; - - MPI_Group_free(&group); - if (group_tmp1 != MPI_GROUP_EMPTY) { - MPI_Group_free(&group_tmp1); - } - if (group_tmp2 != MPI_GROUP_EMPTY) { - MPI_Group_free(&group_tmp2); - } - - output_debug << "done lower PF\n"; - } - - if (jyseps2_1 != jyseps1_2) { - // Upper PF region - // Note need to order processors so that a continuous surface is formed - TRACE("Creating upper PF communicators for xp={:d}", i); - - output_debug << "Creating upper PF communicators for xp = " << i << endl; - - proc[0] = PROC_NUM(i, YPROC(ny_inner)); - proc[1] = PROC_NUM(i, YPROC(jyseps1_2)); - - output_debug << "PF3 " << proc[0] << ", " << proc[1] << endl; - - MPI_Group_range_incl(group_world, 1, &proc, &group_tmp1); - proc[0] = PROC_NUM(i, YPROC(jyseps2_1 + 1)); - proc[1] = PROC_NUM(i, YPROC(ny_inner - 1)); - - output_debug << "PF4 " << proc[0] << ", " << proc[1] << endl; - - MPI_Group_range_incl(group_world, 1, &proc, &group_tmp2); - MPI_Group_union(group_tmp1, group_tmp2, &group); - MPI_Comm_create(BoutComm::get(), group, &comm_tmp); - if (comm_tmp != MPI_COMM_NULL) { - comm_inner = comm_tmp; - if (ixseps_upper == ixseps_outer) { - - output_debug << "-> Inner and middle\n"; - - comm_middle = comm_inner; - } else { - - output_debug << "-> Outer and middle\n"; - - comm_middle = comm_outer; - // MPI_Comm_dup(comm_outer, &comm_middle); - } - } - - output_debug << "Freeing\n"; - - MPI_Group_free(&group); - if (group_tmp1 != MPI_GROUP_EMPTY) { - MPI_Group_free(&group_tmp1); - } - if (group_tmp2 != MPI_GROUP_EMPTY) { - MPI_Group_free(&group_tmp2); - } - - output_debug << "done upper PF\n"; + else{ + std::string mesh_top = toString(mesh_topology); + throw BoutException(_f("Unsupported mesh topology {:s} for communicator creation\n"), mesh_top); } + for (int i = 0; i < NXPE; i++) { // Core region TRACE("Creating core communicators"); + //For CDN and UDN its the inner core region. + group_tmp1 = MPI_GROUP_EMPTY; + group_tmp2 = MPI_GROUP_EMPTY; + if (jyseps2_1 > jyseps1_1) { proc[0] = PROC_NUM(i, YPROC(jyseps1_1 + 1)); proc[1] = PROC_NUM(i, YPROC(jyseps2_1)); @@ -858,30 +1019,29 @@ void BoutMesh::createCommunicators() { group_tmp1 = MPI_GROUP_EMPTY; } - if (jyseps2_2 > jyseps1_2) { - proc[0] = PROC_NUM(i, YPROC(jyseps1_2 + 1)); - proc[1] = PROC_NUM(i, YPROC(jyseps2_2)); + //Only for CDN and UDN outer core region. Add check to ensure only created for these topologies. Maybe topology should go inside second if. + if (mesh_topology == MeshTopology::connected_double_null || mesh_topology == MeshTopology::unconnected_double_null){ + if (jyseps2_2 > jyseps1_2) { + proc[0] = PROC_NUM(i, YPROC(jyseps1_2 + 1)); + proc[1] = PROC_NUM(i, YPROC(jyseps2_2)); - output_debug << "CORE2 " << proc[0] << ", " << proc[1] << endl; + output_debug << "CORE2 " << proc[0] << ", " << proc[1] << endl; - if ((proc[0] < 0) || (proc[1] < 0)) { - group_tmp2 = MPI_GROUP_EMPTY; + if ((proc[0] < 0) || (proc[1] < 0)) { + group_tmp2 = MPI_GROUP_EMPTY; + } else { + MPI_Group_range_incl(group_world, 1, &proc, &group_tmp2); + } } else { - MPI_Group_range_incl(group_world, 1, &proc, &group_tmp2); + // no core region between jyseps1_2 and jyseps2_2 + group_tmp2 = MPI_GROUP_EMPTY; } - } else { - // no core region between jyseps1_2 and jyseps2_2 - group_tmp2 = MPI_GROUP_EMPTY; } MPI_Group_union(group_tmp1, group_tmp2, &group); MPI_Comm_create(BoutComm::get(), group, &comm_tmp); if (comm_tmp != MPI_COMM_NULL) { comm_inner = comm_tmp; - - if (ixseps_inner == ixseps_outer) { - MPI_Comm_dup(comm_inner, &comm_middle); - } } if (group_tmp1 != MPI_GROUP_EMPTY) { @@ -894,8 +1054,10 @@ void BoutMesh::createCommunicators() { } if (ixseps_inner == ixseps_outer) { - // Balanced null, so no middle - MPI_Comm_dup(comm_inner, &comm_middle); + // Balanced null: comm_middle = comm_inner (only for procs that are in the core group) + if (comm_inner != MPI_COMM_NULL) { + MPI_Comm_dup(comm_inner, &comm_middle); + } } else { // Need to handle unbalanced double-null case @@ -966,28 +1128,45 @@ void BoutMesh::createXBoundaries() { return; } - if (PE_XIND == 0) { - // Inner either core or PF + // Get a global index in this processor + const int yg = getGlobalYIndexNoBoundaries(MYG); - // Get a global index in this processor - const int yg = getGlobalYIndexNoBoundaries(MYG); + if (PE_XIND == 0) { + // Inner x face: either core or PF boundary. + // + // For CDN/UDN the y-range (jyseps1_2, jyseps2_2] is the outer core leg, + // so it gets a "core" boundary. + if ((mesh_topology == MeshTopology::connected_double_null) or (mesh_topology == MeshTopology::unconnected_double_null) + or (mesh_topology == MeshTopology::closed_field_line)){ + // CDN/UDN have two core legs; CFL is all core (no X-points). + // All three need both y-ranges checked. + const bool in_core = ((yg > jyseps1_1) and (yg <= jyseps2_1)) + or ((yg > jyseps1_2) and (yg <= jyseps2_2)); + + if (in_core) { + boundary.push_back(new BoundaryRegionXIn("core", ystart, yend, this)); + } else { + boundary.push_back(new BoundaryRegionXIn("pf", ystart, yend, this)); + } + } + else if (mesh_topology == MeshTopology::single_null){ + //SN has only one core region, but it goes from (jyseps1_1, jyseps2_2], ny_inner = jyseps1_2 = jyseps2_1 are not relevant for that case. + const bool in_core = ((yg > jyseps1_1) and (yg <= jyseps2_2)); - if (((yg > jyseps1_1) and (yg <= jyseps2_1)) - or ((yg > jyseps1_2) and (yg <= jyseps2_2))) { - // Core - boundary.push_back(new BoundaryRegionXIn("core", ystart, yend, this)); - } else { - // PF region - boundary.push_back(new BoundaryRegionXIn("pf", ystart, yend, this)); + if (in_core) { + boundary.push_back(new BoundaryRegionXIn("core", ystart, yend, this)); + } else { + boundary.push_back(new BoundaryRegionXIn("pf", ystart, yend, this)); + } } } if (PE_XIND == (NXPE - 1)) { - // Outer SOL - boundary.push_back(new BoundaryRegionXOut("sol", ystart, yend, this)); + boundary.push_back(new BoundaryRegionXOut("sol", ystart, yend, this)); } } +//Unchanged. Not topology dependent. void BoutMesh::createYBoundaries() { if (MYG <= 0) { return; @@ -1025,6 +1204,7 @@ void BoutMesh::createYBoundaries() { } } +// Unchanged. Not topology dependent. std::set BoutMesh::getPossibleBoundaries() const { // Result set: set so it automatically takes care of duplicates std::set all_boundaries{}; @@ -1349,7 +1529,6 @@ comm_handle BoutMesh::sendY(FieldGroup& g, comm_handle handle) { } int BoutMesh::wait(comm_handle handle) { - TRACE("BoutMesh::wait(comm_handle)"); if (handle == nullptr) { return 1; @@ -1507,13 +1686,24 @@ int BoutMesh::wait(comm_handle handle) { * Non-Local Communications ***************************************************************/ -int BoutMesh::getNXPE() { return NXPE; } +int BoutMesh::getNXPE() const { return NXPE; } + +int BoutMesh::getNYPE() const { return NYPE; } -int BoutMesh::getNYPE() { return NYPE; } +int BoutMesh::getNZPE() const { return NZPE; } -int BoutMesh::getXProcIndex() { return PE_XIND; } +int BoutMesh::getXProcIndex() const { return PE_XIND; } -int BoutMesh::getYProcIndex() { return PE_YIND; } +int BoutMesh::getYProcIndex() const { return PE_YIND; } + +int BoutMesh::getZProcIndex() const { return PE_ZIND; } + +int BoutMesh::getProcIndex(int X, int Y, int Z) const { + if (X < 0 || X >= NXPE || Y < 0 || Y >= NYPE || Z < 0 || Z >= NZPE) { + return -1; + } + return Y * NXPE * NZPE + X * NZPE + Z; +} /**************************************************************** * X COMMUNICATIONS @@ -1528,7 +1718,7 @@ bool BoutMesh::lastX() const { return PE_XIND == NXPE - 1; } int BoutMesh::sendXOut(BoutReal* buffer, int size, int tag) { Timer timer("comms"); - int proc {-1}; + int proc{-1}; if (PE_XIND == NXPE - 1) { if (periodicX) { // Wrap around to first processor in X @@ -1540,8 +1730,7 @@ int BoutMesh::sendXOut(BoutReal* buffer, int size, int tag) { proc = PROC_NUM(PE_XIND + 1, PE_YIND); } - mpi->MPI_Send(buffer, size, PVEC_REAL_MPI_TYPE, proc, tag, - BoutComm::get()); + mpi->MPI_Send(buffer, size, PVEC_REAL_MPI_TYPE, proc, tag, BoutComm::get()); return 0; } @@ -1549,7 +1738,7 @@ int BoutMesh::sendXOut(BoutReal* buffer, int size, int tag) { int BoutMesh::sendXIn(BoutReal* buffer, int size, int tag) { Timer timer("comms"); - int proc {-1}; + int proc{-1}; if (PE_XIND == 0) { if (periodicX) { // Wrap around to last processor in X @@ -1561,8 +1750,7 @@ int BoutMesh::sendXIn(BoutReal* buffer, int size, int tag) { proc = PROC_NUM(PE_XIND - 1, PE_YIND); } - mpi->MPI_Send(buffer, size, PVEC_REAL_MPI_TYPE, proc, tag, - BoutComm::get()); + mpi->MPI_Send(buffer, size, PVEC_REAL_MPI_TYPE, proc, tag, BoutComm::get()); return 0; } @@ -1570,7 +1758,7 @@ int BoutMesh::sendXIn(BoutReal* buffer, int size, int tag) { comm_handle BoutMesh::irecvXOut(BoutReal* buffer, int size, int tag) { Timer timer("comms"); - int proc {-1}; + int proc{-1}; if (PE_XIND == NXPE - 1) { if (periodicX) { // Wrap around to first processor in X @@ -1585,8 +1773,8 @@ comm_handle BoutMesh::irecvXOut(BoutReal* buffer, int size, int tag) { // Get a communications handle. Not fussy about size of arrays CommHandle* ch = get_handle(0, 0); - mpi->MPI_Irecv(buffer, size, PVEC_REAL_MPI_TYPE, proc, tag, - BoutComm::get(), ch->request); + mpi->MPI_Irecv(buffer, size, PVEC_REAL_MPI_TYPE, proc, tag, BoutComm::get(), + ch->request); ch->in_progress = true; @@ -1596,7 +1784,7 @@ comm_handle BoutMesh::irecvXOut(BoutReal* buffer, int size, int tag) { comm_handle BoutMesh::irecvXIn(BoutReal* buffer, int size, int tag) { Timer timer("comms"); - int proc {-1}; + int proc{-1}; if (PE_XIND == 0) { if (periodicX) { // Wrap around to last processor in X @@ -1611,8 +1799,8 @@ comm_handle BoutMesh::irecvXIn(BoutReal* buffer, int size, int tag) { // Get a communications handle. Not fussy about size of arrays CommHandle* ch = get_handle(0, 0); - mpi->MPI_Irecv(buffer, size, PVEC_REAL_MPI_TYPE, proc, tag, - BoutComm::get(), ch->request); + mpi->MPI_Irecv(buffer, size, PVEC_REAL_MPI_TYPE, proc, tag, BoutComm::get(), + ch->request); ch->in_progress = true; @@ -1631,33 +1819,53 @@ bool BoutMesh::lastY() const { return PE_YIND == NYPE - 1; } bool BoutMesh::firstY(int xpos) const { int xglobal = getGlobalXIndex(xpos); - int rank; + MPI_Comm comm; if (xglobal < ixseps_inner) { - MPI_Comm_rank(comm_inner, &rank); + comm = comm_inner; } else if (xglobal < ixseps_outer) { - MPI_Comm_rank(comm_middle, &rank); + comm = comm_middle; } else { - MPI_Comm_rank(comm_outer, &rank); + comm = comm_outer; + } + if (comm == MPI_COMM_NULL) { +// const bool lower_target = +// (xpos < DDATA_XSPLIT) ? (DDATA_INDEST < 0) : (DDATA_OUTDEST < 0); +// if (lower_target) { +// return true; +// } + return PE_YIND == 0; } + + int rank; + MPI_Comm_rank(comm, &rank); return rank == 0; } bool BoutMesh::lastY(int xpos) const { int xglobal = getGlobalXIndex(xpos); - int rank; - int size; + MPI_Comm comm; if (xglobal < ixseps_inner) { - MPI_Comm_size(comm_inner, &size); - MPI_Comm_rank(comm_inner, &rank); + comm = comm_inner; } else if (xglobal < ixseps_outer) { - MPI_Comm_size(comm_middle, &size); - MPI_Comm_rank(comm_middle, &rank); + comm = comm_middle; } else { - MPI_Comm_size(comm_outer, &size); - MPI_Comm_rank(comm_outer, &rank); + comm = comm_outer; + } + + if (comm == MPI_COMM_NULL) { +// const bool upper_target = +// (xpos < UDATA_XSPLIT) ? (UDATA_INDEST < 0) : (UDATA_OUTDEST < 0); +// if (upper_target) { +// return true; +// } + return PE_YIND == NYPE - 1; } + + int rank, size; + MPI_Comm_size(comm, &size); + MPI_Comm_rank(comm, &rank); return rank == size - 1; } @@ -1677,35 +1885,32 @@ int BoutMesh::PROC_NUM(int xind, int yind) const { return -1; } - return yind * NXPE + xind; + return (yind * NXPE) + xind; } -/// Returns the global X index given a local index -int BoutMesh::XGLOBAL(BoutReal xloc, BoutReal& xglo) const { - xglo = xloc + PE_XIND * MXSUB; - return static_cast(xglo); +BoutReal BoutMesh::getGlobalXIndex(BoutReal xloc) const { + return xloc + (PE_XIND * MXSUB); } -int BoutMesh::getGlobalXIndex(int xlocal) const { return xlocal + PE_XIND * MXSUB; } +int BoutMesh::getGlobalXIndex(int xlocal) const { return xlocal + (PE_XIND * MXSUB); } int BoutMesh::getGlobalXIndexNoBoundaries(int xlocal) const { - return xlocal + PE_XIND * MXSUB - MXG; + return xlocal + (PE_XIND * MXSUB) - MXG; } -int BoutMesh::getLocalXIndex(int xglobal) const { return xglobal - PE_XIND * MXSUB; } +int BoutMesh::getLocalXIndex(int xglobal) const { return xglobal - (PE_XIND * MXSUB); } int BoutMesh::getLocalXIndexNoBoundaries(int xglobal) const { - return xglobal - PE_XIND * MXSUB + MXG; + return xglobal - (PE_XIND * MXSUB) + MXG; } -int BoutMesh::YGLOBAL(BoutReal yloc, BoutReal& yglo) const { - yglo = yloc + PE_YIND * MYSUB - MYG; - return static_cast(yglo); +BoutReal BoutMesh::getGlobalYIndex(BoutReal yloc) const { + return yloc + (PE_YIND * MYSUB) - MYG; } int BoutMesh::getGlobalYIndex(int ylocal) const { - int yglobal = ylocal + PE_YIND * MYSUB; - if (jyseps1_2 > jyseps2_1 and PE_YIND * MYSUB + 2 * MYG + 1 > ny_inner) { + int yglobal = ylocal + (PE_YIND * MYSUB); + if (jyseps1_2 > jyseps2_1 and (PE_YIND * MYSUB) + (2 * MYG) + 1 > ny_inner) { // Double null, and we are past the upper target yglobal += 2 * MYG; } @@ -1713,12 +1918,12 @@ int BoutMesh::getGlobalYIndex(int ylocal) const { } int BoutMesh::getGlobalYIndexNoBoundaries(int ylocal) const { - return ylocal + PE_YIND * MYSUB - MYG; + return ylocal + (PE_YIND * MYSUB) - MYG; } int BoutMesh::getLocalYIndex(int yglobal) const { - int ylocal = yglobal - PE_YIND * MYSUB; - if (jyseps1_2 > jyseps2_1 and PE_YIND * MYSUB + 2 * MYG + 1 > ny_inner) { + int ylocal = yglobal - (PE_YIND * MYSUB); + if (jyseps1_2 > jyseps2_1 and (PE_YIND * MYSUB) + (2 * MYG) + 1 > ny_inner) { // Double null, and we are past the upper target ylocal -= 2 * MYG; } @@ -1726,19 +1931,25 @@ int BoutMesh::getLocalYIndex(int yglobal) const { } int BoutMesh::getLocalYIndexNoBoundaries(int yglobal) const { - return yglobal - PE_YIND * MYSUB + MYG; + return yglobal - (PE_YIND * MYSUB) + MYG; } -int BoutMesh::YGLOBAL(int yloc, int yproc) const { return yloc + yproc * MYSUB - MYG; } +int BoutMesh::YGLOBAL(int yloc, int yproc) const { return yloc + (yproc * MYSUB) - MYG; } -int BoutMesh::YLOCAL(int yglo, int yproc) const { return yglo - yproc * MYSUB + MYG; } +int BoutMesh::YLOCAL(int yglo, int yproc) const { return yglo - (yproc * MYSUB) + MYG; } -int BoutMesh::getGlobalZIndex(int zlocal) const { return zlocal; } +int BoutMesh::getGlobalZIndex(int zlocal) const { return zlocal + (PE_ZIND * MZSUB); } -int BoutMesh::getGlobalZIndexNoBoundaries(int zlocal) const { return zlocal; } +int BoutMesh::getGlobalZIndexNoBoundaries(int zlocal) const { + return zlocal + (PE_ZIND * MZSUB) - MZG; +} int BoutMesh::getLocalZIndex(int zglobal) const { return zglobal; } +BoutReal BoutMesh::getGlobalZIndex(BoutReal zloc) const { + return zloc + (PE_ZIND * MZSUB); +} + int BoutMesh::getLocalZIndexNoBoundaries(int zglobal) const { return zglobal; } int BoutMesh::YPROC(int yind) const { @@ -1806,17 +2017,19 @@ BoutMesh::BoutMesh(int input_nx, int input_ny, int input_nz, int mxg, int myg, i BoutMesh::BoutMesh(int input_nx, int input_ny, int input_nz, int mxg, int myg, int nxpe, int nype, int pe_xind, int pe_yind, bool symmetric_X, bool symmetric_Y, - bool periodicX_, int ixseps1_, int ixseps2_, int jyseps1_1_, + bool periodic_X_, int ixseps1_, int ixseps2_, int jyseps1_1_, int jyseps2_1_, int jyseps1_2_, int jyseps2_2_, int ny_inner_, bool create_regions) : nx(input_nx), ny(input_ny), nz(input_nz), NPES(nxpe * nype), - MYPE(nxpe * pe_yind + pe_xind), PE_YIND(pe_yind), NYPE(nype), NZPE(1), - ixseps1(ixseps1_), ixseps2(ixseps2_), symmetricGlobalX(symmetric_X), + MYPE((nxpe * pe_yind) + pe_xind), PE_XIND(pe_xind), NXPE(nxpe), PE_YIND(pe_yind), + NYPE(nype), ixseps1(ixseps1_), ixseps2(ixseps2_), symmetricGlobalX(symmetric_X), symmetricGlobalY(symmetric_Y), MXG(mxg), MYG(myg), MZG(0) { - NXPE = nxpe; - PE_XIND = pe_xind; - periodicX = periodicX_; + + periodicX = periodic_X_; setYDecompositionIndices(jyseps1_1_, jyseps2_1_, jyseps1_2_, jyseps2_2_, ny_inner_); + mesh_topology = getMeshTopology(jyseps1_1, jyseps2_1, + jyseps1_2, jyseps2_2, + ny_inner, ixseps1, ixseps2); setDerivedGridSizes(); topology(); if (create_regions) { @@ -2092,7 +2305,7 @@ void BoutMesh::topology() { true); // Twist-shift this connection set_connection(jyseps1_1, jyseps2_2 + 1, 0, ixseps1); // No twist-shift in PF region - } else { + } else if (mesh_topology == MeshTopology::connected_double_null || mesh_topology == MeshTopology::unconnected_double_null) { /*************** DOUBLE NULL OPERATION *******************/ /* UPPER LEGS: Do not have to be the same length as each other or lower legs, but do have to have an integer number @@ -2139,13 +2352,47 @@ void BoutMesh::topology() { // Add target plates at the top add_target(ny_inner - 1, 0, nx); - } - - if ((ixseps_inner > 0) + } + + // Additional limiters + // Each limiter needs 3 indices: A Y index, start and end X indices + int limiter_count = 0; + Mesh::get(limiter_count, "limiter_count", 0); + if (limiter_count > 0) { + std::vector limiter_yinds; + if (!source->get(this, limiter_yinds, "limiter_yinds", limiter_count)) { + throw BoutException("Couldn't read limiter_yinds vector of length {} from mesh", + limiter_count); + } + std::vector limiter_xstarts; + if (!source->get(this, limiter_xstarts, "limiter_xstarts", limiter_count)) { + throw BoutException("Couldn't read limiter_xstarts vector of length {} from mesh", + limiter_count); + } + std::vector limiter_xends; + if (!source->get(this, limiter_xends, "limiter_xends", limiter_count)) { + throw BoutException("Couldn't read limiter_xend vector of length {} from mesh", + limiter_count); + } + + for (int i = 0; i < limiter_count; ++i) { + const int yind = limiter_yinds[i]; + const int xstart = limiter_xstarts[i]; + const int xend = limiter_xends[i]; + output_info.write("Adding a limiter between y={} and {}. X indices {} to {}\n", + yind, yind + 1, xstart, xend); + add_target(yind, xstart, xend); + } + } + + if (mesh_topology == MeshTopology::unconnected_double_null || mesh_topology == MeshTopology::connected_double_null) { + //This is for DN topologies + if ((ixseps_inner > 0) && (((PE_YIND * MYSUB > jyseps1_1) && (PE_YIND * MYSUB <= jyseps2_1)) || ((PE_YIND * MYSUB > jyseps1_2) && (PE_YIND * MYSUB <= jyseps2_2)))) { - MYPE_IN_CORE = true; /* processor is in the core */ - } + MYPE_IN_CORE = true; /* processor is in the core */ + } + } if (DDATA_XSPLIT > LocalNx) { DDATA_XSPLIT = LocalNx; @@ -2311,72 +2558,112 @@ void BoutMesh::overlapHandleMemory(BoutMesh* yup, BoutMesh* ydown, BoutMesh* xin * Communication utilities ****************************************************************/ -int BoutMesh::pack_data(const std::vector& var_list, int xge, int xlt, - int yge, int ylt, BoutReal* buffer) { +int BoutMesh::pack_data(const std::vector& var_list, int xge, int xlt, int yge, + int ylt, BoutReal* buffer) const { + using FieldType = Field::FieldType; int len = 0; + const int zge = 0; + const int zlt = LocalNz; - /// Loop over variables for (const auto& var : var_list) { - if (var->is3D()) { - // 3D variable - auto* var3d_ref_ptr = dynamic_cast(var); + switch (var->field_type()) { + case FieldType::field3d: { + const auto* var3d_ref_ptr = dynamic_cast(var); ASSERT0(var3d_ref_ptr != nullptr); - auto& var3d_ref = *var3d_ref_ptr; + const auto& var3d_ref = *var3d_ref_ptr; ASSERT2(var3d_ref.isAllocated()); - for (int jx = xge; jx != xlt; jx++) { + for (int jx = xge; jx < xlt; jx++) { for (int jy = yge; jy < ylt; jy++) { - for (int jz = 0; jz < LocalNz; jz++, len++) { + for (int jz = zge; jz < zlt; jz++, len++) { buffer[len] = var3d_ref(jx, jy, jz); } } } - } else { - // 2D variable - auto* var2d_ref_ptr = dynamic_cast(var); + break; + } + case FieldType::field2d: { + const auto* var2d_ref_ptr = dynamic_cast(var); ASSERT0(var2d_ref_ptr != nullptr); - auto& var2d_ref = *var2d_ref_ptr; + const auto& var2d_ref = *var2d_ref_ptr; ASSERT2(var2d_ref.isAllocated()); - for (int jx = xge; jx != xlt; jx++) { + for (int jx = xge; jx < xlt; jx++) { for (int jy = yge; jy < ylt; jy++, len++) { buffer[len] = var2d_ref(jx, jy); } } + break; + } + case FieldType::fieldperp: { + const auto* varperp_ref_ptr = dynamic_cast(var); + ASSERT0(varperp_ref_ptr != nullptr); + const auto& varperp_ref = *varperp_ref_ptr; + ASSERT2(varperp_ref.isAllocated()); + for (int jx = xge; jx < xlt; jx++) { + for (int jz = zge; jz < zlt; jz++, len++) { + buffer[len] = varperp_ref(jx, jz); + } + } + break; + } } } - return (len); + return len; } -int BoutMesh::unpack_data(const std::vector& var_list, int xge, int xlt, - int yge, int ylt, BoutReal* buffer) { +int BoutMesh::unpack_data(const std::vector& var_list, int xge, int xlt, int yge, + int ylt, const BoutReal* buffer) const { + using FieldType = Field::FieldType; int len = 0; + const int zge = 0; + const int zlt = LocalNz; - /// Loop over variables for (const auto& var : var_list) { - if (var->is3D()) { - // 3D variable - auto& var3d_ref = *dynamic_cast(var); - for (int jx = xge; jx != xlt; jx++) { + switch (var->field_type()) { + case FieldType::field3d: { + auto* var3d_ref_ptr = dynamic_cast(var); + ASSERT0(var3d_ref_ptr != nullptr); + auto& var3d_ref = *var3d_ref_ptr; + ASSERT2(var3d_ref.isAllocated()); + for (int jx = xge; jx < xlt; jx++) { for (int jy = yge; jy < ylt; jy++) { - for (int jz = 0; jz < LocalNz; jz++, len++) { + for (int jz = zge; jz < zlt; jz++, len++) { var3d_ref(jx, jy, jz) = buffer[len]; } } } - } else { - // 2D variable - auto& var2d_ref = *dynamic_cast(var); - for (int jx = xge; jx != xlt; jx++) { + break; + } + case FieldType::field2d: { + auto* var2d_ref_ptr = dynamic_cast(var); + ASSERT0(var2d_ref_ptr != nullptr); + auto& var2d_ref = *var2d_ref_ptr; + ASSERT2(var2d_ref.isAllocated()); + for (int jx = xge; jx < xlt; jx++) { for (int jy = yge; jy < ylt; jy++, len++) { var2d_ref(jx, jy) = buffer[len]; } } + break; + } + case FieldType::fieldperp: { + auto* varperp_ref_ptr = dynamic_cast(var); + ASSERT0(varperp_ref_ptr != nullptr); + auto& varperp_ref = *varperp_ref_ptr; + ASSERT2(varperp_ref.isAllocated()); + for (int jx = xge; jx < xlt; jx++) { + for (int jz = zge; jz < zlt; jz++, len++) { + varperp_ref(jx, jz) = buffer[len]; + } + } + break; + } } } - return (len); + return len; } /**************************************************************** @@ -2438,41 +2725,42 @@ int BoutMesh::ySize(int xpos) const { int xglobal = getGlobalXIndex(xpos); int yglobal = getGlobalYIndexNoBoundaries(MYG); - if ((xglobal < ixseps_lower) && ((yglobal <= jyseps1_1) || (yglobal > jyseps2_2))) { - // Lower PF region - return (jyseps1_1 + 1) + (ny - jyseps2_2); + //Old divisions working for all other topologies. + if ((xglobal < ixseps_lower) && ((yglobal <= jyseps1_1) || (yglobal > jyseps2_2))) { + // Lower PF region + return (jyseps1_1 + 1) + (ny - jyseps2_2); - } else if ((xglobal < ixseps_upper) && (yglobal > jyseps2_1) - && (yglobal >= jyseps1_2)) { - // Upper PF region - return jyseps1_2 - jyseps2_1; + } else if ((xglobal < ixseps_upper) && (yglobal > jyseps2_1) + && (yglobal >= jyseps1_2)) { + // Upper PF region + return jyseps1_2 - jyseps2_1; - } else if (xglobal < ixseps_inner) { - // Core - return (jyseps2_1 - jyseps1_1) + (jyseps2_2 - jyseps1_2); + } else if (xglobal < ixseps_inner) { + return (jyseps2_1 - jyseps1_1) + (jyseps2_2 - jyseps1_2); - } else if (jyseps2_1 == jyseps1_2) { - // Single null, so in the SOL - return ny; + } else if (jyseps2_1 == jyseps1_2) { + // Single null, so in the SOL + return ny; - } else if ((xglobal >= ixseps_inner) && (xglobal < ixseps_outer)) { - // Intermediate SOL in DND + } else if ((xglobal >= ixseps_inner) && (xglobal < ixseps_outer)) { + // Intermediate SOL in DN - if (ixseps_lower < ixseps_upper) { - // Connects to lower divertor - return (jyseps2_1 + 1) + (ny - jyseps1_2); - } else { - // Connects to upper divertor - return jyseps2_2 - jyseps1_1; + if (ixseps_lower < ixseps_upper) { + // Connects to lower divertor + return (jyseps2_1 + 1) + (ny - jyseps1_2); + } else { + // Connects to upper divertor + return jyseps2_2 - jyseps1_1; + } + } else if (yglobal < ny_inner) { + // Inner SOL + return ny_inner; } - } else if (yglobal < ny_inner) { - // Inner SOL - return ny_inner; - } - // Outer SOL - return ny - ny_inner; + // Outer SOL + return ny - ny_inner; } +//Unchanged, not dependent on topology MPI_Comm BoutMesh::getYcomm(int xpos) const { int xglobal = getGlobalXIndex(xpos); @@ -2487,7 +2775,7 @@ MPI_Comm BoutMesh::getYcomm(int xpos) const { /**************************************************************** * Range iteration ****************************************************************/ - +//Not topology dependent void BoutMesh::addBoundaryRegions() { std::list all_boundaries; ///< Keep track of all boundary regions @@ -3070,47 +3358,6 @@ void BoutMesh::addBoundaryPar(std::shared_ptr bndry, par_boundary[static_cast(BoundaryParType::all)].push_back(bndry); } -Field3D BoutMesh::smoothSeparatrix(const Field3D& f) { - Field3D result{emptyFrom(f)}; - if ((ixseps_inner > 0) && (ixseps_inner < nx - 1)) { - if (XPROC(ixseps_inner) == PE_XIND) { - int x = getLocalXIndex(ixseps_inner); - for (int y = 0; y < LocalNy; y++) { - for (int z = 0; z < LocalNz; z++) { - result(x, y, z) = 0.5 * (f(x, y, z) + f(x - 1, y, z)); - } - } - } - if (XPROC(ixseps_inner - 1) == PE_XIND) { - int x = getLocalXIndex(ixseps_inner - 1); - for (int y = 0; y < LocalNy; y++) { - for (int z = 0; z < LocalNz; z++) { - result(x, y, z) = 0.5 * (f(x, y, z) + f(x + 1, y, z)); - } - } - } - } - if ((ixseps_outer > 0) && (ixseps_outer < nx - 1) && (ixseps_outer != ixseps_inner)) { - if (XPROC(ixseps_outer) == PE_XIND) { - int x = getLocalXIndex(ixseps_outer); - for (int y = 0; y < LocalNy; y++) { - for (int z = 0; z < LocalNz; z++) { - result(x, y, z) = 0.5 * (f(x, y, z) + f(x - 1, y, z)); - } - } - } - if (XPROC(ixseps_outer - 1) == PE_XIND) { - int x = getLocalXIndex(ixseps_outer - 1); - for (int y = 0; y < LocalNy; y++) { - for (int z = 0; z < LocalNz; z++) { - result(x, y, z) = 0.5 * (f(x, y, z) + f(x + 1, y, z)); - } - } - } - } - return result; -} - BoutReal BoutMesh::GlobalX(int jx) const { if (symmetricGlobalX) { // With this definition the boundary sits dx/2 away form the first/last inner points @@ -3122,8 +3369,7 @@ BoutReal BoutMesh::GlobalX(int jx) const { BoutReal BoutMesh::GlobalX(BoutReal jx) const { // Get global X index as a BoutReal - BoutReal xglo; - XGLOBAL(jx, xglo); + const BoutReal xglo = getGlobalXIndex(jx); if (symmetricGlobalX) { // With this definition the boundary sits dx/2 away form the first/last inner points @@ -3145,7 +3391,7 @@ BoutReal BoutMesh::GlobalY(int jy) const { yi -= jyseps1_1 + 0.5 + (jyseps1_2 - jyseps2_1); } return yi / nycore; - } + } int ly = getGlobalYIndexNoBoundaries(jy); // global poloidal index across subdomains int nycore = (jyseps2_1 - jyseps1_1) + (jyseps2_2 - jyseps1_2); @@ -3175,49 +3421,75 @@ BoutReal BoutMesh::GlobalY(int jy) const { } BoutReal BoutMesh::GlobalY(BoutReal jy) const { + + // Get global Y index as a BoutReal + //BoutReal yglo; + //YGLOBAL(jy, yglo); PREVIOUS BEHAVIOUR - // Get global Y index as a BoutReal - BoutReal yglo; - YGLOBAL(jy, yglo); - if (symmetricGlobalY) { - BoutReal yi = yglo; - int nycore = (jyseps2_1 - jyseps1_1) + (jyseps2_2 - jyseps1_2); + // Get global Y index as a BoutReal + BoutReal yglo = getGlobalYIndex(jy); - if (yi < ny_inner) { - // before upper target - yi -= jyseps1_1 + 0.5; - } else { - // after upper target - yi -= jyseps1_1 + 0.5 + (jyseps1_2 - jyseps2_1); + if (symmetricGlobalY) { + BoutReal yi = yglo; + int nycore = (jyseps2_1 - jyseps1_1) + (jyseps2_2 - jyseps1_2); + + if (yi < ny_inner) { + // before upper target + yi -= jyseps1_1 + 0.5; + } else { + // after upper target + yi -= jyseps1_1 + 0.5 + (jyseps1_2 - jyseps2_1); + } + return yi / nycore; } - return yi / nycore; - } - int nycore = (jyseps2_1 - jyseps1_1) + (jyseps2_2 - jyseps1_2); + int nycore = (jyseps2_1 - jyseps1_1) + (jyseps2_2 - jyseps1_2); - if (MYPE_IN_CORE) { - // Turn yglo into an index over the core cells onyglo - if (yglo <= jyseps2_1) { - yglo -= jyseps1_1 + 1; + if (MYPE_IN_CORE) { + // Turn yglo into an index over the core cells onyglo + if (yglo <= jyseps2_1) { + yglo -= jyseps1_1 + 1; + } else { + yglo -= jyseps1_1 + 1 + (jyseps1_2 - jyseps2_1); + } } else { - yglo -= jyseps1_1 + 1 + (jyseps1_2 - jyseps2_1); - } - } else { - // Not in core. Need to get the last "core" value - if (yglo <= jyseps1_1) { - // Inner lower leg - yglo = 0; - } else if ((yglo > jyseps2_1) && (yglo <= jyseps1_2)) { - // Upper legs - yglo = jyseps2_1 - jyseps1_1; - } else if (yglo > jyseps2_2) { - // Outer lower leg - yglo = nycore; + // Not in core. Need to get the last "core" value + if (yglo <= jyseps1_1) { + // Inner lower leg + yglo = 0; + } else if ((yglo > jyseps2_1) && (yglo <= jyseps1_2)) { + // Upper legs + yglo = jyseps2_1 - jyseps1_1; + } else if (yglo > jyseps2_2) { + // Outer lower leg + yglo = nycore; + } } + + return yglo / static_cast(nycore); +} + +BoutReal BoutMesh::GlobalZ(int jz) const { + if (symmetricGlobalZ) { + // With this definition the boundary sits dz/2 away form the first/last inner points + return (0.5 + getGlobalZIndexNoBoundaries(jz) - (nz - MZ) * 0.5) + / static_cast(MZ); } + return static_cast(getGlobalZIndexNoBoundaries(jz)) + / static_cast(MZ); +} + +BoutReal BoutMesh::GlobalZ(BoutReal jz) const { - return yglo / static_cast(nycore); + // Get global Z index as a BoutReal + const BoutReal zglo = getGlobalZIndex(jz); + + if (symmetricGlobalZ) { + // With this definition the boundary sits dz/2 away form the first/last inner points + return (0.5 + zglo - (nz - MZ) * 0.5) / static_cast(MZ); + } + return zglo / static_cast(MZ); } void BoutMesh::outputVars(Options& output_options) { @@ -3248,6 +3520,7 @@ void BoutMesh::outputVars(Options& output_options) { output_options["jyseps2_1"].force(jyseps2_1, "BoutMesh"); output_options["jyseps2_2"].force(jyseps2_2, "BoutMesh"); output_options["ny_inner"].force(ny_inner, "BoutMesh"); + output_options["mesh_topology"].force(mesh_topology, "BoutMesh"); getCoordinates()->outputVars(output_options); diff --git a/src/mesh/impls/bout/boutmesh.hxx b/src/mesh/impls/bout/boutmesh.hxx index cc674d401a..b1ff12d259 100644 --- a/src/mesh/impls/bout/boutmesh.hxx +++ b/src/mesh/impls/bout/boutmesh.hxx @@ -4,19 +4,29 @@ #include "mpi.h" +#include "bout/bout_types.hxx" #include "bout/unused.hxx" #include -#include #include #include #include #include +class Field; + /// Implementation of Mesh (mostly) compatible with BOUT /// /// Topology and communications compatible with BOUT /// conventions. + +BOUT_ENUM_CLASS(MeshTopology, + closed_field_line, // Closed field line + single_null, // Single null + unconnected_double_null, // Unconnected double null + connected_double_null); // Connected double null + + class BoutMesh : public Mesh { public: BoutMesh(GridDataSource* s, Options* options = nullptr); @@ -25,6 +35,10 @@ public: /// Read in the mesh from data sources int load() override; + MeshTopology getMeshTopology(int jyseps1_1_, int jyseps2_1_, int jyseps1_2_, + int jyseps2_2_, int ny_inner_, int ixseps1_, + int ixseps2_); + ///////////////////////////////////////////// // Communicate variables @@ -58,10 +72,13 @@ public: ///////////////////////////////////////////// // non-local communications - int getNXPE() override; ///< The number of processors in the X direction - int getNYPE() override; ///< The number of processors in the Y direction - int getXProcIndex() override; ///< This processor's index in X direction - int getYProcIndex() override; ///< This processor's index in Y direction + int getNXPE() const override; ///< The number of processors in the X direction + int getNYPE() const override; ///< The number of processors in the Y direction + int getNZPE() const override; ///< The number of processors in the Z direction + int getXProcIndex() const override; ///< This processor's index in X direction + int getYProcIndex() const override; ///< This processor's index in Y direction + int getZProcIndex() const override; ///< This processor's index in Z direction + int getProcIndex(int X, int Y, int Z) const override; ///////////////////////////////////////////// // X communications @@ -105,6 +122,7 @@ public: MPI_Comm getXcomm(int UNUSED(jy)) const override { return comm_x; } /// Return communicator containing all processors in Y MPI_Comm getYcomm(int xpos) const override; + MPI_Comm getXZcomm() const override { return comm_xz; } /// Is local X index \p jx periodic in Y? /// @@ -156,6 +174,9 @@ public: RangeIterator iterateBndryUpperInnerY() const override; RangeIterator iterateBndryUpperOuterY() const override; + bool hasBndryLowerY() const override { return has_boundary_lower_y; } + bool hasBndryUpperY() const override { return has_boundary_upper_y; } + // Boundary regions std::vector getBoundaries() override; std::vector> @@ -164,15 +185,15 @@ public: BoundaryParType type) override; std::set getPossibleBoundaries() const override; - Field3D smoothSeparatrix(const Field3D& f) override; - int getNx() const { return nx; } int getNy() const { return ny; } BoutReal GlobalX(int jx) const override; BoutReal GlobalY(int jy) const override; + BoutReal GlobalZ(int jz) const override; BoutReal GlobalX(BoutReal jx) const override; BoutReal GlobalY(BoutReal jy) const override; + BoutReal GlobalZ(BoutReal jz) const override; BoutReal getIxseps1() const { return ixseps1; } BoutReal getIxseps2() const { return ixseps2; } @@ -205,7 +226,7 @@ protected: /// `getPossibleBoundaries`. \p create_regions controls whether or /// not the various `Region`s are created on the new mesh BoutMesh(int input_nx, int input_ny, int input_nz, int mxg, int myg, int nxpe, int nype, - int pe_xind, int pe_yind, bool symmetric_X, bool symmetric_Y, bool periodic_X, + int pe_xind, int pe_yind, bool symmetric_X, bool symmetric_Y, bool periodic_X_, int ixseps1_, int ixseps2_, int jyseps1_1_, int jyseps2_1_, int jyseps1_2_, int jyseps2_2_, int ny_inner_, bool create_regions = true); @@ -275,6 +296,7 @@ protected: /// Create the various sub-communicators void createCommunicators(); + /// Create the boundary regions in X void createXBoundaries(); @@ -294,16 +316,24 @@ private: int NPES; ///< Number of processors int MYPE; ///< Rank of this processor + int PE_XIND; ///< X index of this processor + int NXPE; ///< Number of processors in the X direction + int PE_YIND; ///< Y index of this processor - int NYPE; // Number of processors in the Y direction + int NYPE; ///< Number of processors in the Y direction - int NZPE; + int PE_ZIND{0}; ///< Z index of this processor + int NZPE{1}; ///< Number of processors in the Z direction /// Is this processor in the core region? bool MYPE_IN_CORE{false}; - int XGLOBAL(BoutReal xloc, BoutReal& xglo) const; - int YGLOBAL(BoutReal yloc, BoutReal& yglo) const; + /// Returns the global X index given a local index + BoutReal getGlobalXIndex(BoutReal xloc) const; + /// Returns the global Y index given a local index + BoutReal getGlobalYIndex(BoutReal yloc) const; + /// Returns the global Z index given a local index + BoutReal getGlobalZIndex(BoutReal zloc) const; // Topology int ixseps1, ixseps2, jyseps1_1, jyseps2_1, jyseps1_2, jyseps2_2; @@ -312,6 +342,8 @@ private: std::vector ShiftAngle; ///< Angle for twist-shift location + MeshTopology mesh_topology; + protected: // These are protected so we can make them public in the test suite // for testing @@ -354,8 +386,9 @@ private: // Settings bool TwistShift; // Use a twist-shift condition in core? - bool symmetricGlobalX; ///< Use a symmetric definition in GlobalX() function - bool symmetricGlobalY; + bool symmetricGlobalX; ///< Use a symmetric definition in `GlobalX()` function + bool symmetricGlobalY; ///< Use a symmetric definition in `GlobalY()` function + bool symmetricGlobalZ{false}; ///< Use a symmetric definition in `GlobalZ()` function int zperiod; BoutReal ZMIN, ZMAX; // Range of the Z domain (in fractions of 2pi) @@ -394,12 +427,16 @@ protected: /// Adds 2D and 3D regions for boundaries void addBoundaryRegions(); + //void findValidProcessorNum(int ny, int nx); + private: std::vector boundary; // Vector of boundary regions std::array>, static_cast(BoundaryParType::SIZE)> par_boundary; // Vector of parallel boundary regions + bool has_boundary_lower_y{false}; + bool has_boundary_upper_y{false}; ////////////////////////////////////////////////// // Communications @@ -439,6 +476,8 @@ private: /// Communicator containing all processors in X MPI_Comm comm_x{MPI_COMM_NULL}; + /// Communicator for all processors in an XZ plane + MPI_Comm comm_xz{MPI_COMM_NULL}; ////////////////////////////////////////////////// // Surface communications @@ -460,12 +499,12 @@ private: void post_receiveY(CommHandle& ch); /// Take data from objects and put into a buffer - int pack_data(const std::vector& var_list, int xge, int xlt, int yge, - int ylt, BoutReal* buffer); + int pack_data(const std::vector& var_list, int xge, int xlt, int yge, int ylt, + BoutReal* buffer) const; /// Copy data from a buffer back into the fields + int unpack_data(const std::vector& var_list, int xge, int xlt, int yge, int ylt, + const BoutReal* buffer) const; - int unpack_data(const std::vector& var_list, int xge, int xlt, int yge, - int ylt, BoutReal* buffer); }; namespace { @@ -483,10 +522,28 @@ struct CheckMeshResult { /// Check that \p total_processors can be decomposed into \p /// num_y_processors in Y for the given `BoutMesh` topology parameters -CheckMeshResult checkBoutMeshYDecomposition(int num_y_processors, int ny, - int num_y_guards, int jyseps1_1, - int jyseps2_1, int jyseps1_2, int jyseps2_2, - int ny_inner); +CheckMeshResult checkBoutMeshYDecomposition( + int num_y_processors, int ny, + int num_y_guards, + int jyseps1_1, int jyseps2_1, + int jyseps1_2, int jyseps2_2, + int ny_inner); + +// New topology-aware +CheckMeshResult checkBoutMeshYDecomposition( + int num_y_processors, int ny, + int num_y_guards, + int jyseps1_1, int jyseps2_1, + int jyseps1_2, int jyseps2_2, + int ny_inner, + MeshTopology mesh_topology); + + CheckMeshResult findValidProcessorNum(int ny, int nx, int NPES, int NYPE = 1, int NXPE = 1); + + CheckMeshResult findValidYDecomposition(int ny, int NPES, int NYPE, + int jyseps1_1, int jyseps2_1, + int jyseps1_2, int jyseps2_2, + int ny_inner, MeshTopology mesh_topology); } // namespace bout #endif // BOUT_BOUTMESH_H diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx index cb6cb8410c..6d9bc9f681 100644 --- a/src/mesh/mesh.cxx +++ b/src/mesh/mesh.cxx @@ -377,16 +377,28 @@ void Mesh::communicate(FieldPerp& f) { wait(recv[1]); } -int Mesh::msg_len(const std::vector& var_list, int xge, int xlt, int yge, - int ylt) { +int Mesh::msg_len(const std::vector& var_list, int xge, int xlt, int yge, + int ylt) const { int len = 0; + using FieldType = Field::FieldType; + + const auto x_length = xlt - xge; + const auto y_length = ylt - yge; + const auto z_length = LocalNz; + /// Loop over variables for (const auto& var : var_list) { - if (var->is3D()) { - len += (xlt - xge) * (ylt - yge) * LocalNz * var->elementSize(); - } else { - len += (xlt - xge) * (ylt - yge) * var->elementSize(); + switch (var->field_type()) { + case FieldType::field3d: + len += x_length * y_length * z_length * var->elementSize(); + break; + case FieldType::field2d: + len += x_length * y_length * var->elementSize(); + break; + case FieldType::fieldperp: + len += x_length * z_length * var->elementSize(); + break; } } @@ -408,34 +420,6 @@ int Mesh::ySize(int jx) const { return all; } -bool Mesh::hasBndryLowerY() { - static bool calc = false, answer; - if (calc) { - return answer; // Already calculated - } - - int mybndry = static_cast(!(iterateBndryLowerY().isDone())); - int allbndry; - mpi->MPI_Allreduce(&mybndry, &allbndry, 1, MPI_INT, MPI_BOR, getXcomm(yend)); - answer = static_cast(allbndry); - calc = true; - return answer; -} - -bool Mesh::hasBndryUpperY() { - static bool calc = false, answer; - if (calc) { - return answer; // Already calculated - } - - int mybndry = static_cast(!(iterateBndryUpperY().isDone())); - int allbndry; - mpi->MPI_Allreduce(&mybndry, &allbndry, 1, MPI_INT, MPI_BOR, getXcomm(ystart)); - answer = static_cast(allbndry); - calc = true; - return answer; -} - int Mesh::localSize3D() { if (localNumCells3D < 0) { const int xs = firstX() ? xstart - 1 : xstart; diff --git a/tests/MMS/GBS/runtest-slab2d b/tests/MMS/GBS/runtest-slab2d index 97c7e11459..221c34e815 100755 --- a/tests/MMS/GBS/runtest-slab2d +++ b/tests/MMS/GBS/runtest-slab2d @@ -11,7 +11,6 @@ from boutdata.collect import collect from numpy import sqrt, max, abs, mean, array, log, concatenate - build_and_log("MMS test") diff --git a/tests/MMS/GBS/runtest-slab3d b/tests/MMS/GBS/runtest-slab3d index f193583e47..45acd53af7 100755 --- a/tests/MMS/GBS/runtest-slab3d +++ b/tests/MMS/GBS/runtest-slab3d @@ -14,7 +14,6 @@ from numpy import sqrt, max, abs, mean, array, log, concatenate import pickle - build_and_log("MMS test") # List of NX values to use diff --git a/tests/MMS/advection/runtest b/tests/MMS/advection/runtest index 1290a4c247..9a0a691287 100755 --- a/tests/MMS/advection/runtest +++ b/tests/MMS/advection/runtest @@ -20,7 +20,6 @@ import pickle import time - if __name__ == "__main__": build_and_log("MMS advection") diff --git a/tests/MMS/diffusion/runtest b/tests/MMS/diffusion/runtest index a25b3519a6..ef6b335f70 100755 --- a/tests/MMS/diffusion/runtest +++ b/tests/MMS/diffusion/runtest @@ -15,7 +15,6 @@ from boutdata.collect import collect from numpy import sqrt, max, abs, mean, array, log - build_and_log("MMS diffusion test") # List of NX values to use diff --git a/tests/MMS/diffusion2/Z/plot_error.py b/tests/MMS/diffusion2/Z/plot_error.py index 5c4e1164c1..4b4b666758 100644 --- a/tests/MMS/diffusion2/Z/plot_error.py +++ b/tests/MMS/diffusion2/Z/plot_error.py @@ -4,6 +4,7 @@ @author: yolen """ + import boutdata import matplotlib.pyplot as plt diff --git a/tests/MMS/diffusion2/runtest b/tests/MMS/diffusion2/runtest index 546be3b2c1..d1b8072fc8 100755 --- a/tests/MMS/diffusion2/runtest +++ b/tests/MMS/diffusion2/runtest @@ -20,7 +20,6 @@ from numpy import sqrt, max, abs, mean, array, log from os.path import join - build_and_log("MMS diffusion test") # List of input directories diff --git a/tests/MMS/fieldalign/runtest.broken b/tests/MMS/fieldalign/runtest.broken index e09fbeccd5..bfb2023c1d 100755 --- a/tests/MMS/fieldalign/runtest.broken +++ b/tests/MMS/fieldalign/runtest.broken @@ -18,7 +18,6 @@ from os.path import join import time - build_and_log("MMS test") # nxlist = [256, 128, 64, 32, 16, 8] # do in reverse order to save disk space diff --git a/tests/MMS/hw/runtest b/tests/MMS/hw/runtest index b6abadc3b9..62740fcbbd 100755 --- a/tests/MMS/hw/runtest +++ b/tests/MMS/hw/runtest @@ -15,7 +15,6 @@ from boutdata.collect import collect from numpy import sqrt, max, abs, mean, array, log, concatenate - build_and_log("MMS test") # List of NX values to use diff --git a/tests/MMS/laplace/runtest b/tests/MMS/laplace/runtest index 98eb62da87..80c8bf08af 100755 --- a/tests/MMS/laplace/runtest +++ b/tests/MMS/laplace/runtest @@ -14,7 +14,6 @@ from boutdata.collect import collect from numpy import sqrt, max, abs, mean, array, log, concatenate - build_and_log("MMS test") # List of NX values to use diff --git a/tests/MMS/spatial/diffusion/runtest b/tests/MMS/spatial/diffusion/runtest index 000cf8f09b..8d4e759d59 100755 --- a/tests/MMS/spatial/diffusion/runtest +++ b/tests/MMS/spatial/diffusion/runtest @@ -22,7 +22,6 @@ from os.path import join import matplotlib.pyplot as plt - build_and_log("MMS diffusion test") # List of input directories diff --git a/tests/MMS/wave-1d-y/runtest b/tests/MMS/wave-1d-y/runtest index 9a4de5a3bb..332d6864b0 100755 --- a/tests/MMS/wave-1d-y/runtest +++ b/tests/MMS/wave-1d-y/runtest @@ -20,7 +20,6 @@ from sys import stdout from numpy import sqrt, max, abs, mean, array, log, concatenate, pi - build_and_log("Making MMS wave test") # List of NX values to use diff --git a/tests/MMS/wave-1d/runtest b/tests/MMS/wave-1d/runtest index f13fb7bd11..1adb7c7d6d 100755 --- a/tests/MMS/wave-1d/runtest +++ b/tests/MMS/wave-1d/runtest @@ -15,7 +15,6 @@ from boutdata.collect import collect from numpy import sqrt, max, abs, mean, array, log, concatenate - build_and_log("Making MMS wave test") # List of NX values to use diff --git a/tests/integrated/test-cyclic/runtest b/tests/integrated/test-cyclic/runtest index 5e31da6239..c68f9b033e 100755 --- a/tests/integrated/test-cyclic/runtest +++ b/tests/integrated/test-cyclic/runtest @@ -17,7 +17,6 @@ from boututils.run_wrapper import build_and_log, shell, launch from boutdata.collect import collect from sys import stdout, exit - build_and_log("Cyclic Reduction test") flags = ["", "nsys=2", "nsys=5 periodic", "nsys=7 n=10"] diff --git a/tests/integrated/test-griddata/runtest b/tests/integrated/test-griddata/runtest index 1c57bf5be7..6aef955023 100755 --- a/tests/integrated/test-griddata/runtest +++ b/tests/integrated/test-griddata/runtest @@ -13,7 +13,6 @@ from boutdata.collect import collect import numpy as np from sys import stdout, exit - build_and_log("griddata test") for nproc in [1]: diff --git a/tests/integrated/test-gyro/runtest b/tests/integrated/test-gyro/runtest index 52ab0de4ab..7a0209fa57 100755 --- a/tests/integrated/test-gyro/runtest +++ b/tests/integrated/test-gyro/runtest @@ -25,7 +25,6 @@ from boutdata.collect import collect import numpy as np from sys import stdout, exit - build_and_log("Gyro-average inversion test") # Read benchmark values diff --git a/tests/integrated/test-invpar/runtest b/tests/integrated/test-invpar/runtest index de7f028528..53f6015fe9 100755 --- a/tests/integrated/test-invpar/runtest +++ b/tests/integrated/test-invpar/runtest @@ -12,7 +12,6 @@ from boututils.run_wrapper import build_and_log, shell, launch from boutdata.collect import collect from sys import stdout, exit - build_and_log("parallel inversion test") flags_src = [ diff --git a/tests/integrated/test-laplace/runtest b/tests/integrated/test-laplace/runtest index e54e46d0d8..d1e9c3d1c9 100755 --- a/tests/integrated/test-laplace/runtest +++ b/tests/integrated/test-laplace/runtest @@ -41,7 +41,6 @@ from boutdata.collect import collect import numpy as np from sys import stdout, exit - build_and_log("Laplacian inversion test") # Read benchmark values diff --git a/tests/integrated/test-multigrid_laplace/runtest b/tests/integrated/test-multigrid_laplace/runtest index 4a7455f80b..f325185c46 100755 --- a/tests/integrated/test-multigrid_laplace/runtest +++ b/tests/integrated/test-multigrid_laplace/runtest @@ -22,7 +22,6 @@ from boututils.run_wrapper import build_and_log, shell, launch_safe from boutdata.collect import collect from sys import exit - build_and_log("multigrid Laplacian inversion test") print("Running multigrid Laplacian inversion test") diff --git a/tests/integrated/test-multigrid_laplace/runtest_multiple_grids b/tests/integrated/test-multigrid_laplace/runtest_multiple_grids index 6817120b13..739d15f7d4 100755 --- a/tests/integrated/test-multigrid_laplace/runtest_multiple_grids +++ b/tests/integrated/test-multigrid_laplace/runtest_multiple_grids @@ -18,7 +18,6 @@ from boututils.run_wrapper import shell, build_and_log, launch_safe from boutdata.collect import collect from sys import exit - build_and_log("Multigrid Laplacian inversion test") print("Running multigrid Laplacian inversion test") diff --git a/tests/integrated/test-multigrid_laplace/runtest_unsheared b/tests/integrated/test-multigrid_laplace/runtest_unsheared index cda68f2167..d0dca3a808 100755 --- a/tests/integrated/test-multigrid_laplace/runtest_unsheared +++ b/tests/integrated/test-multigrid_laplace/runtest_unsheared @@ -18,7 +18,6 @@ from boututils.run_wrapper import shell, build_and_log, launch_safe from boutdata.collect import collect from sys import exit - build_and_log("Making multigrid Laplacian inversion test") print("Running multigrid Laplacian inversion test") diff --git a/tests/integrated/test-naulin-laplace/runtest b/tests/integrated/test-naulin-laplace/runtest index f972eab6cc..145257c23b 100755 --- a/tests/integrated/test-naulin-laplace/runtest +++ b/tests/integrated/test-naulin-laplace/runtest @@ -22,7 +22,6 @@ from boututils.run_wrapper import build_and_log, shell, launch_safe from boutdata.collect import collect from sys import exit - build_and_log("LaplaceNaulin inversion test") print("Running LaplaceNaulin inversion test") diff --git a/tests/integrated/test-naulin-laplace/runtest_multiple_grids b/tests/integrated/test-naulin-laplace/runtest_multiple_grids index c0281c3a4e..f56a80677e 100755 --- a/tests/integrated/test-naulin-laplace/runtest_multiple_grids +++ b/tests/integrated/test-naulin-laplace/runtest_multiple_grids @@ -18,7 +18,6 @@ from boututils.run_wrapper import shell, build_and_log, launch_safe from boutdata.collect import collect from sys import exit - build_and_log("Making LaplaceNaulin inversion test") print("Running LaplaceNaulin inversion test") diff --git a/tests/integrated/test-naulin-laplace/runtest_unsheared b/tests/integrated/test-naulin-laplace/runtest_unsheared index 8f47f33026..a845a9d92f 100755 --- a/tests/integrated/test-naulin-laplace/runtest_unsheared +++ b/tests/integrated/test-naulin-laplace/runtest_unsheared @@ -18,7 +18,6 @@ from boututils.run_wrapper import shell, build_and_log, launch_safe from boutdata.collect import collect from sys import exit - build_and_log("LaplaceNaulin inversion test") print("Running LaplaceNaulin inversion test") diff --git a/tests/integrated/test-petsc_laplace/runtest b/tests/integrated/test-petsc_laplace/runtest index ac248c4ce7..83e1006338 100755 --- a/tests/integrated/test-petsc_laplace/runtest +++ b/tests/integrated/test-petsc_laplace/runtest @@ -30,7 +30,6 @@ from boutdata.collect import collect # import numpy as np from sys import stdout, exit - build_and_log("PETSc Laplacian inversion test") print("Running PETSc Laplacian inversion test") diff --git a/tests/integrated/test-petsc_laplace_MAST-grid/runtest b/tests/integrated/test-petsc_laplace_MAST-grid/runtest index 5a4cbf875a..3be39949c3 100755 --- a/tests/integrated/test-petsc_laplace_MAST-grid/runtest +++ b/tests/integrated/test-petsc_laplace_MAST-grid/runtest @@ -28,7 +28,6 @@ from boututils.run_wrapper import build_and_log, shell, launch_safe from boutdata.collect import collect from sys import stdout, exit - build_and_log( "PETSc Laplacian inversion test with non-identity metric (taken from grid for MAST SOL)" ) diff --git a/tests/integrated/test-region-iterator/runtest b/tests/integrated/test-region-iterator/runtest index e5825285de..2277ef15e9 100755 --- a/tests/integrated/test-region-iterator/runtest +++ b/tests/integrated/test-region-iterator/runtest @@ -17,7 +17,6 @@ from boututils.run_wrapper import build_and_log, launch_safe from boutdata.collect import collect from sys import exit - build_and_log("Region Iterator test") flags = [""] diff --git a/tests/integrated/test-restarting/runtest b/tests/integrated/test-restarting/runtest index 262aa818b6..daf19ebac6 100755 --- a/tests/integrated/test-restarting/runtest +++ b/tests/integrated/test-restarting/runtest @@ -5,7 +5,6 @@ from boutdata.collect import collect import numpy as np from sys import stdout, exit - build_and_log("restart test") # Run once for 10 timesteps diff --git a/tests/integrated/test-smooth/runtest b/tests/integrated/test-smooth/runtest index bd7296341a..91ae793164 100755 --- a/tests/integrated/test-smooth/runtest +++ b/tests/integrated/test-smooth/runtest @@ -16,7 +16,6 @@ from boutdata.collect import collect import numpy as np from sys import stdout, exit - build_and_log("smoothing operator test") # Read benchmark values diff --git a/tests/integrated/test-squash/runtest b/tests/integrated/test-squash/runtest index 2f35e00589..8bf358ee8e 100755 --- a/tests/integrated/test-squash/runtest +++ b/tests/integrated/test-squash/runtest @@ -9,7 +9,6 @@ import argparse import re import os.path - # requires: all_tests # requires: netcdf # cores: 4 diff --git a/tests/integrated/test-yupdown/runtest b/tests/integrated/test-yupdown/runtest index 34fcd36496..13d7365b53 100755 --- a/tests/integrated/test-yupdown/runtest +++ b/tests/integrated/test-yupdown/runtest @@ -6,7 +6,6 @@ from sys import exit from numpy import max, abs - build_and_log("parallel slices test") failed = False diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx index e6f78f8767..1bc8cfe4c0 100644 --- a/tests/unit/fake_mesh.hxx +++ b/tests/unit/fake_mesh.hxx @@ -98,10 +98,16 @@ public: return nullptr; } int wait(comm_handle UNUSED(handle)) override { return 0; } - int getNXPE() override { return 1; } - int getNYPE() override { return 1; } - int getXProcIndex() override { return 1; } - int getYProcIndex() override { return 1; } + int getNXPE() const override { return 1; } + int getNYPE() const override { return 1; } + int getNZPE() const override { return 1; } + int getXProcIndex() const override { return 1; } + int getYProcIndex() const override { return 1; } + int getZProcIndex() const override { return 0; } + int getProcIndex([[maybe_unused]] int X, [[maybe_unused]] int Y, + [[maybe_unused]] int Z) const override { + return 0; + } bool firstX() const override { return true; } bool lastX() const override { return true; } int sendXOut(BoutReal* UNUSED(buffer), int UNUSED(size), int UNUSED(tag)) override { @@ -120,6 +126,7 @@ public: } MPI_Comm getXcomm(int UNUSED(jy)) const override { return BoutComm::get(); } MPI_Comm getYcomm(int UNUSED(jx)) const override { return BoutComm::get(); } + MPI_Comm getXZcomm() const override { return BoutComm::get(); } bool periodicY(int UNUSED(jx)) const override { return true; } bool periodicY(int UNUSED(jx), BoutReal& UNUSED(ts)) const override { return true; } int numberOfYBoundaries() const override { return 1; } @@ -143,6 +150,8 @@ public: RangeIterator iterateBndryLowerInnerY() const override { return RangeIterator(); } RangeIterator iterateBndryUpperOuterY() const override { return RangeIterator(); } RangeIterator iterateBndryUpperInnerY() const override { return RangeIterator(); } + bool hasBndryLowerY() const override { return false; } + bool hasBndryUpperY() const override { return false; } void addBoundary(BoundaryRegion* region) override { boundaries.push_back(region); } std::vector getBoundaries() override { return boundaries; } std::vector> @@ -151,8 +160,10 @@ public: } BoutReal GlobalX(int jx) const override { return jx; } BoutReal GlobalY(int jy) const override { return jy; } + BoutReal GlobalZ(int jz) const override { return jz; } BoutReal GlobalX(BoutReal jx) const override { return jx; } BoutReal GlobalY(BoutReal jy) const override { return jy; } + BoutReal GlobalZ(BoutReal jz) const override { return jz; } int getGlobalXIndex(int) const override { return 0; } int getGlobalXIndexNoBoundaries(int) const override { return 0; } int getGlobalYIndex(int y) const override { return y; } diff --git a/tests/unit/fake_parallel_mesh.hxx b/tests/unit/fake_parallel_mesh.hxx index 805dcb2a0a..f9b32c217b 100644 --- a/tests/unit/fake_parallel_mesh.hxx +++ b/tests/unit/fake_parallel_mesh.hxx @@ -163,7 +163,7 @@ public: /// (representing other processors) to see if a corresponding field /// has been created for them which can be used to communicate guard /// cells with. - void registerField(FieldData& f, int id) { + void registerField(Field& f, int id) { registeredFields.emplace(&f, id); registeredFieldIds.emplace(id, &f); } @@ -259,7 +259,7 @@ private: FakeParallelMesh *yUpMesh, *yDownMesh, *xInMesh, *xOutMesh; bool communicatingX = false, communicatingY = false; std::map registeredFields; - std::map registeredFieldIds; + std::map registeredFieldIds; std::map registeredFieldPerps; std::map registeredFieldPerpIds; std::unique_ptr mpiSmart; diff --git a/tests/unit/field/test_field.cxx b/tests/unit/field/test_field.cxx index 40d675e5b3..0e7b0d318f 100644 --- a/tests/unit/field/test_field.cxx +++ b/tests/unit/field/test_field.cxx @@ -24,7 +24,7 @@ class FieldSubClass : public Field { FieldSubClass(Mesh* localmesh, CELL_LOC location_in, DirectionTypes directions_in) : Field(localmesh, location_in, directions_in) {} - bool is3D() const override { return false; } + FieldType field_type() const override { return FieldType::field2d; } int size() const override { return 42; } }; } // namespace diff --git a/tests/unit/mesh/test_boutmesh.cxx b/tests/unit/mesh/test_boutmesh.cxx index fd2ba50d28..a35699ccb2 100644 --- a/tests/unit/mesh/test_boutmesh.cxx +++ b/tests/unit/mesh/test_boutmesh.cxx @@ -47,6 +47,7 @@ class BoutMeshExposer : public BoutMesh { using BoutMesh::XPROC; using BoutMesh::YDecompositionIndices; using BoutMesh::YPROC; + using BoutMesh::getMeshTopology; }; /// Minimal parameters need to construct a grid useful for testing @@ -239,24 +240,61 @@ BoutMeshParameters createDisconnectedDoubleNull(const BoutMeshGridInfo& grid) { grid.total_ny - grid.local_ny - 1, ny_inner}}; } +BoutMeshParameters createSnowflake(const BoutMeshGridInfo& grid) { + // Need at least 6 y-subdomains for a minimal snowflake + if (grid.nype < 6) { + throw BoutException( + "createSnowflake: Not enough processors for snowflake topology " + "(nype={}, needs at least 6)", + grid.nype); + } + + if ((grid.total_nx / 2) + 4 > grid.total_nx) { + throw BoutException( + "createSnowflake: Not enough points in x-direction " + "(need ixseps2 = ((nxpe * (local_nx - 2)) + 2) / 2 + 4 = {} to " + "be less than total_nx = (nxpe * (local_nx - 2)) + 2 = {}; nxpe={}, local_nx={}", + (grid.total_nx / 2) + 4, grid.total_nx, grid.nxpe, grid.local_nx); +} + + const int ny_inner = 4 * grid.local_ny; + // Separatrix indices + const int jyseps1_1 = grid.local_ny - 1; + const int jyseps2_1 = ny_inner - 2 * grid.local_ny - 1; + const int jyseps1_2 = ny_inner - grid.local_ny - 1; + const int jyseps2_2 = grid.total_ny - grid.local_ny - 1; + + return { + grid, + // X separatrices (same as standard snowflake assumption) + {grid.total_nx / 2, grid.total_nx / 2 + 4}, + // Y separatrices + ny_inner + {jyseps1_1, + jyseps2_1, + jyseps1_2, + jyseps2_2, + ny_inner} + }; +} + + //////////////////////////////////////////////////////////// // Start of tests -TEST(BoutMeshTest, NullOptionsCheck) { +struct BoutMeshTest : public ::testing::Test { + WithQuietOutput debug{output_debug}; WithQuietOutput info{output_info}; WithQuietOutput warn{output_warn}; + WithQuietOutput progress{output_progress}; +}; +TEST_F(BoutMeshTest, NullOptionsCheck) { EXPECT_NO_THROW(BoutMesh mesh(new FakeGridDataSource, nullptr)); } // Not a great test as it's not specific to the thing we want to test, // and can also take a whopping ~300ms! -TEST(BoutMeshTest, SingleCoreDecomposition) { - WithQuietOutput debug{output_debug}; - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - WithQuietOutput progress{output_progress}; - +TEST_F(BoutMeshTest, SingleCoreDecomposition) { Options options{}; options["ny"] = 1; options["nx"] = 4; @@ -271,6 +309,7 @@ TEST(BoutMeshTest, SingleCoreDecomposition) { bout::globals::mpi = nullptr; } + struct SetYDecompositionTestParameters { BoutMeshExposer::YDecompositionIndices input; BoutMeshExposer::YDecompositionIndices expected; @@ -327,13 +366,14 @@ TEST_P(BoutMeshSetYDecompositionTest, BasicTest) { EXPECT_EQ(mesh.numberOfXPoints, params.number_of_X_points); } -TEST(BoutMeshTest, SetYDecompositionIndicesJyseps22LowInconsistent) { - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, SetYDecompositionIndicesJyseps22LowInconsistent) { BoutMeshExposer mesh(1, 24, 1, 1, 1); EXPECT_THROW(mesh.setYDecompositionIndices({3, 7, 32, 8, 12}), BoutException); } +//New bit: + struct DecompositionTestParameters { int total_processors; int num_y_processors; @@ -343,20 +383,24 @@ struct DecompositionTestParameters { std::string expected_message; // Expect this fragment to be in the result.reason for bad // decompositions std::string name; + MeshTopology mesh_topology; // New: topology enum }; DecompositionTestParameters makeDecompositionTestParameters(const BoutMeshParameters& inputs, - const std::string& name) { + const std::string& name, + MeshTopology mesh_topology = MeshTopology::unconnected_double_null) { // default to unconnected_double_null return {inputs.grid.total_processors, inputs.grid.nype, inputs.grid.total_ny, inputs.grid.num_y_guards, inputs.y_indices, "", - name}; + name, + mesh_topology}; } + std::ostream& operator<<(std::ostream& out, const DecompositionTestParameters& value) { return out << fmt::format( "DecompositionTestParameters{{" @@ -369,10 +413,12 @@ std::ostream& operator<<(std::ostream& out, const DecompositionTestParameters& v "jyseps1_2 = {}, " "jyseps2_2 = {}, " "ny_inner = {}, " + "topology = {}, " "expected_message = {} }}", value.total_processors, value.num_y_processors, value.ny, value.num_y_guards, value.indices.jyseps1_1, value.indices.jyseps2_1, value.indices.jyseps1_2, - value.indices.jyseps2_2, value.indices.ny_inner, value.expected_message); + value.indices.jyseps2_2, value.indices.ny_inner, toString(value.mesh_topology), + value.expected_message); } std::string DecompositionTestParametersToString( @@ -380,6 +426,7 @@ std::string DecompositionTestParametersToString( return param.param.name; } + struct BoutMeshDecompositionTest : public testing::TestWithParam { virtual ~BoutMeshDecompositionTest() = default; @@ -388,36 +435,39 @@ struct BoutMeshDecompositionTest INSTANTIATE_TEST_SUITE_P( GoodDecompositions, BoutMeshDecompositionTest, ::testing::Values( - DecompositionTestParameters{1, 1, 1, 1, {-1, 0, 0, 0, 0}, "", "OnePoint"}, - DecompositionTestParameters{1, 1, 8, 1, {-1, 4, 4, 7, 4}, "", "EightPoints"}, + DecompositionTestParameters{1, 1, 1, 1, {-1, 0, 0, 0, 0}, "", "OnePoint", MeshTopology::single_null}, + DecompositionTestParameters{1, 1, 8, 1, {-1, 4, 4, 7, 4}, "", "EightPoints", MeshTopology::single_null}, DecompositionTestParameters{ - 2, 1, 8, 1, {-1, 4, 4, 7, 4}, "", "EightPointsTwoCores"}, + 2, 1, 8, 1, {-1, 4, 4, 7, 4}, "", "EightPointsTwoCores", MeshTopology::single_null}, DecompositionTestParameters{ - 2, 2, 8, 1, {-1, 4, 4, 7, 4}, "", "EightPointsTwoCoresNYPE2"}, + 2, 2, 8, 1, {-1, 4, 4, 7, 4}, "", "EightPointsTwoCoresNYPE2", MeshTopology::single_null}, // The following should basically all work by construction - makeDecompositionTestParameters(createCore({4, 4, 2, 2, 1, 1}), "Core"), - makeDecompositionTestParameters(createSOL({4, 4, 2, 2, 1, 1}), "SOL"), - makeDecompositionTestParameters(createLimiter({4, 4, 2, 2, 1, 1}), "Limiter"), - makeDecompositionTestParameters(createXPoint({4, 4, 2, 2, 1, 4}), "XPoint"), + makeDecompositionTestParameters(createCore({4, 4, 2, 2, 1, 1}), "Core", MeshTopology::single_null), + makeDecompositionTestParameters(createSOL({4, 4, 2, 2, 1, 1}), "SOL", MeshTopology::single_null), + makeDecompositionTestParameters(createLimiter({4, 4, 2, 2, 1, 1}), "Limiter", MeshTopology::single_null), + makeDecompositionTestParameters(createXPoint({4, 4, 2, 2, 1, 4}), "XPoint", MeshTopology::single_null), makeDecompositionTestParameters(createSingleNull({4, 4, 2, 2, 1, 3}), - "SingleNull"), + "SingleNull", MeshTopology::single_null), makeDecompositionTestParameters(createDoubleNull({4, 4, 2, 2, 1, 6}), - "DoubleNull"), + "DoubleNull", MeshTopology::connected_double_null), makeDecompositionTestParameters(createDisconnectedDoubleNull({12, 4, 2, 2, 1, 6}), - "DisconnectedDoubleNull")), + "DisconnectedDoubleNull", MeshTopology::unconnected_double_null)), DecompositionTestParametersToString); TEST_P(BoutMeshDecompositionTest, CheckYDecomposition) { const auto params = GetParam(); auto result = bout::checkBoutMeshYDecomposition( - params.num_y_processors, params.ny, 1, params.indices.jyseps1_1, - params.indices.jyseps2_1, params.indices.jyseps1_2, params.indices.jyseps2_2, - params.indices.ny_inner); + params.num_y_processors, params.ny, 1, + params.indices.jyseps1_1, params.indices.jyseps2_1, + params.indices.jyseps1_2, params.indices.jyseps2_2, + params.indices.ny_inner, + params.mesh_topology); // <- pass topology EXPECT_TRUE(result.success); EXPECT_TRUE(result.reason.empty()); } + using BadBoutMeshDecompositionTest = BoutMeshDecompositionTest; INSTANTIATE_TEST_SUITE_P( @@ -433,40 +483,115 @@ INSTANTIATE_TEST_SUITE_P( BadDoubleNull, BadBoutMeshDecompositionTest, ::testing::Values( DecompositionTestParameters{ - 1, 1, 4, 1, {3, 5, 6, 10, 0}, "Core region jyseps2_1", "CoreRegion1"}, + 1, 1, 4, 1, {3, 5, 6, 10, 0}, "Core region jyseps2_1", "CoreRegion1", MeshTopology::unconnected_double_null}, DecompositionTestParameters{ - 1, 1, 4, 1, {3, 7, 8, 11, 0}, "Core region jyseps2_2", "CoreRegion2"}, + 1, 1, 4, 1, {3, 7, 8, 11, 0}, "Core region jyseps2_2", "CoreRegion2", MeshTopology::unconnected_double_null}, DecompositionTestParameters{ - 1, 1, 4, 1, {3, 7, 8, 12, 11}, "leg region ny_inner", "UpperLeg1"}, + 1, 1, 4, 1, {3, 7, 8, 12, 11}, "leg region ny_inner", "UpperLeg1", MeshTopology::unconnected_double_null}, DecompositionTestParameters{ - 1, 1, 4, 1, {3, 7, 8, 12, 8}, "leg region jyseps1_2-ny_inner+1", "UpperLeg2"}, + 1, 1, 4, 1, {3, 7, 8, 12, 8}, "leg region jyseps1_2-ny_inner+1", "UpperLeg2", MeshTopology::unconnected_double_null}, DecompositionTestParameters{ - 1, 6, 25, 1, {3, 7, 15, 19, 12}, "leg region ny-jyseps2_2-1", "LegRegion"}), + 1, 6, 25, 1, {3, 7, 15, 19, 12}, "leg region ny-jyseps2_2-1", "LegRegion", MeshTopology::unconnected_double_null}), DecompositionTestParametersToString); INSTANTIATE_TEST_SUITE_P( BadSingleNull, BadBoutMeshDecompositionTest, ::testing::Values( DecompositionTestParameters{ - 1, 1, 4, 1, {3, 4, 4, 6, 0}, "Core region jyseps2_2-jyseps1_1", "CoreRegion"}, + 1, 1, 4, 1, {3, 4, 4, 6, 0}, "Core region jyseps2_2-jyseps1_1", "CoreRegion", MeshTopology::single_null}, DecompositionTestParameters{ - 1, 3, 13, 1, {3, 4, 4, 7, 0}, "leg region ny-jyseps2_2-1", "LegRegion"}), + 1, 3, 13, 1, {3, 4, 4, 7, 0}, "leg region ny-jyseps2_2-1", "LegRegion", MeshTopology::single_null}), DecompositionTestParametersToString); TEST_P(BadBoutMeshDecompositionTest, BadSingleCoreYDecomposition) { const auto params = GetParam(); auto result = bout::checkBoutMeshYDecomposition( - params.num_y_processors, params.ny, params.num_y_guards, params.indices.jyseps1_1, - params.indices.jyseps2_1, params.indices.jyseps1_2, params.indices.jyseps2_2, - params.indices.ny_inner); + params.num_y_processors, params.ny, params.num_y_guards, + params.indices.jyseps1_1, params.indices.jyseps2_1, + params.indices.jyseps1_2, params.indices.jyseps2_2, + params.indices.ny_inner, + params.mesh_topology); // <- pass topology using ::testing::HasSubstr; EXPECT_FALSE(result.success); + //Ask Peter about baddecomtest EXPECT_THAT(result.reason, HasSubstr(params.expected_message)); } -TEST(BoutMeshTest, ChooseProcessorSplitBadNXPE) { +TEST(BoutMeshDecompositionTest, InvalidYDecompositionBecuaseofTopologyUDN) { + int ny = 18; + int num_y_processors = 9; + int num_y_guards = 1; + + int jyseps1_1_start = 1; + int jyseps2_1_start = 1; + int jyseps1_2_start = 17; + int jyseps2_2_start = 1; + int ny_inner_start = 1; + + MeshTopology mesh_topology = MeshTopology::unconnected_double_null; + + auto result = bout::findValidYDecomposition(ny, num_y_processors, num_y_guards, + jyseps1_1_start, jyseps2_1_start, + jyseps1_2_start, jyseps2_2_start, + ny_inner_start, mesh_topology); + EXPECT_FALSE(result.success); +} + + +TEST(BoutMeshDecompositionTest, BasicValidProcessDecompositionDefaults) { + // 8x6 grid, up to 16 processors + auto result = bout::findValidProcessorNum(/*ny=*/8, /*nx=*/6, /*NPES=*/16); + using ::testing::HasSubstr; + EXPECT_TRUE(result.success); + EXPECT_THAT(result.reason, HasSubstr("NPES=16")); + EXPECT_THAT(result.reason, HasSubstr("NXPE=2")); + EXPECT_THAT(result.reason, HasSubstr("NYPE=8")); +} + +TEST(BoutMeshDecompositionTest, RespectsNXPE) { + int NXPE=2; + auto result = bout::findValidProcessorNum(/*ny=*/8, /*nx=*/8, /*NPES=*/16, + NXPE); + using ::testing::HasSubstr; + EXPECT_TRUE(result.success); + EXPECT_THAT(result.reason, HasSubstr("NPES=16")); + EXPECT_THAT(result.reason, HasSubstr("NXPE=2")); + EXPECT_THAT(result.reason, HasSubstr("NYPE=8")); +} + +TEST(BoutMeshDecompositionTest, RespectsNYPE) { + int NYPE=16; + auto result = bout::findValidProcessorNum(/*ny=*/16, /*nx=*/8, /*NPES=*/16, + NYPE); + using ::testing::HasSubstr; + EXPECT_TRUE(result.success); + EXPECT_THAT(result.reason, HasSubstr("NPES=16")); + EXPECT_THAT(result.reason, HasSubstr("NXPE=1")); + EXPECT_THAT(result.reason, HasSubstr("NYPE=16")); +} + + +TEST(BoutMeshDecompositionTest, NoValidDecomposition) { + // Prime sizes, limited processors + auto result = bout::findValidProcessorNum(/*ny=*/7, /*nx=*/8, /*NPES=*/5); + using ::testing::HasSubstr; + EXPECT_FALSE(result.success); + EXPECT_THAT(result.reason, HasSubstr("No valid processor decomposition found")); +} + +TEST(BoutMeshDecompositionTest, SingleProcessorOnly) { + auto result = bout::findValidProcessorNum(/*ny=*/10, /*nx=*/10, /*NPES=*/1); + using ::testing::HasSubstr; + EXPECT_TRUE(result.success); + EXPECT_THAT(result.reason, HasSubstr("NPES=1")); +} + + //End of new bit + +//End of the test +TEST_F(BoutMeshTest, ChooseProcessorSplitBadNXPE) { WithQuietOutput info{output_info}; Options options{{"NXPE", 3}}; @@ -475,7 +600,7 @@ TEST(BoutMeshTest, ChooseProcessorSplitBadNXPE) { EXPECT_THROW(mesh.chooseProcessorSplit(options), BoutException); } -TEST(BoutMeshTest, ChooseProcessorSplitBadNYPE) { +TEST_F(BoutMeshTest, ChooseProcessorSplitBadNYPETooManyYProcs) { WithQuietOutput info{output_info}; Options options{{"NYPE", 7}}; @@ -484,11 +609,28 @@ TEST(BoutMeshTest, ChooseProcessorSplitBadNYPE) { EXPECT_THROW(mesh.chooseProcessorSplit(options), BoutException); } -TEST(BoutMeshTest, ChooseProcessorSplitNXPE) { +TEST_F(BoutMeshTest, ChooseProcessorSplitBadNXPENotDivisibleByNYPE) { + WithQuietOutput info{output_info}; + Options options{{"NXPE", 5}}; + + BoutMeshExposer mesh(4, 24, 1, 1, 1, 8); + + EXPECT_THROW(mesh.chooseProcessorSplit(options), BoutException); +} + +TEST_F(BoutMeshTest, ChooseProcessorSplitBadNYPENotDivisibleByNYPE) { WithQuietOutput info{output_info}; + Options options{{"NYPE", 5}}; + + BoutMeshExposer mesh(5, 5, 1, 1, 1, 8); + + EXPECT_THROW(mesh.chooseProcessorSplit(options), BoutException); +} + +TEST_F(BoutMeshTest, ChooseProcessorSplitNXPE) { Options options{{"NXPE", 4}}; - BoutMeshExposer mesh(1, 24, 1, 1, 1, 8); + BoutMeshExposer mesh(4, 24, 1, 1, 1, 8); EXPECT_NO_THROW(mesh.chooseProcessorSplit(options)); @@ -496,8 +638,7 @@ TEST(BoutMeshTest, ChooseProcessorSplitNXPE) { EXPECT_EQ(mesh.getNYPE(), 2); } -TEST(BoutMeshTest, ChooseProcessorSplitBadNXPENotEnoughGuards) { - WithQuietOutput info{output_info}; +TEST_F(BoutMeshTest, ChooseProcessorSplitBadNXPENotEnoughGuards) { Options options{{"NXPE", 4}}; BoutMeshExposer mesh(1, 24, 1, 1, 13, 8); @@ -505,8 +646,7 @@ TEST(BoutMeshTest, ChooseProcessorSplitBadNXPENotEnoughGuards) { EXPECT_THROW(mesh.chooseProcessorSplit(options), BoutException); } -TEST(BoutMeshTest, ChooseProcessorSplitNYPE) { - WithQuietOutput info{output_info}; +TEST_F(BoutMeshTest, ChooseProcessorSplitNYPE) { Options options{{"NYPE", 4}}; BoutMeshExposer mesh(1, 24, 1, 1, 1, 8); @@ -517,6 +657,33 @@ TEST(BoutMeshTest, ChooseProcessorSplitNYPE) { EXPECT_EQ(mesh.getNYPE(), 4); } +TEST(getMeshTopologyTest, ReturnsCFLWhenNoXPoints) { + BoutMeshExposer mesh(8, 8, 1, 1, 1); + mesh.numberOfXPoints = 0; + EXPECT_EQ(mesh.getMeshTopology(-1, 2, 3, 10, 5, 6, 7), MeshTopology::closed_field_line); +} + +TEST(getMeshTopologyTest, ReturnsSNWhenOneXPoint) { + BoutMeshExposer mesh(8, 8, 1, 1, 1); + mesh.numberOfXPoints = 1; + EXPECT_EQ(mesh.getMeshTopology(1, 2, 2, 4, 5, 6, 7), MeshTopology::single_null); +} + + +TEST(getMeshTopologyTest, ReturnsUDNWhenTwoXPointsDifferentIndices) { + BoutMeshExposer mesh(8, 8, 1, 1, 1); + mesh.numberOfXPoints = 2; + // ny_inner not between jyseps1_2 and jyseps2_2 + EXPECT_EQ(mesh.getMeshTopology(0, 0, 10, 20, 25, 1, 2), MeshTopology::unconnected_double_null); +} + +TEST(getMeshTopologyTest, ReturnsCDNWhenTwoXPointsSameIndices) { + BoutMeshExposer mesh(8, 8, 1, 1, 1); + mesh.numberOfXPoints = 2; + // ny_inner not between jyseps1_2 and jyseps2_2 but ixseps1 == ixseps2 + EXPECT_EQ(mesh.getMeshTopology(0, 0, 10, 20, 25, 1, 1), MeshTopology::connected_double_null); +} + struct FindProcessorParameters { int total_processors; int nx; @@ -698,9 +865,7 @@ TEST_P(BoutMeshProcNumTest, ProcNum) { EXPECT_EQ(result, params.expected_result); } -TEST(BoutMeshTest, YProc) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, YProc) { // 2x2 processors, 3x3x1 (not including guards) on each processor BoutMeshExposer mesh(5, 3, 1, 2, 2, 0, 0); @@ -716,9 +881,7 @@ TEST(BoutMeshTest, YProc) { EXPECT_EQ(mesh.YPROC(7), -1); } -TEST(BoutMeshTest, XProc) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, XProc) { // 2x2 processors, 3x3x1 (not including guards) on each processor BoutMeshExposer mesh(5, 3, 1, 2, 2, 0, 0); @@ -735,9 +898,7 @@ TEST(BoutMeshTest, XProc) { // one example, so probably fine } -TEST(BoutMeshTest, GetGlobalXIndex) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, GetGlobalXIndex) { // 2x2 processors, 3x3x1 (not including guards) on each processor // Boundaries are included in the global index @@ -781,9 +942,7 @@ TEST(BoutMeshTest, GetGlobalXIndex) { EXPECT_EQ(mesh11.getGlobalXIndex(4), 7); } -TEST(BoutMeshTest, GetGlobalXIndexNoBoundaries) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, GetGlobalXIndexNoBoundaries) { // 2x2 processors, 3x3x1 (not including guards) on each processor // Global indices start counting from the first non-boundary point @@ -827,10 +986,7 @@ TEST(BoutMeshTest, GetGlobalXIndexNoBoundaries) { EXPECT_EQ(mesh11.getGlobalXIndexNoBoundaries(4), 6); } -TEST(BoutMeshTest, GlobalXIntSymmetricX) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, GlobalXIntSymmetricX) { BoutMeshExposer mesh01(4, 3, 1, 2, 2, 0, 1); EXPECT_EQ(mesh01.GlobalX(0), -0.125); EXPECT_EQ(mesh01.GlobalX(1), 0.125); @@ -839,10 +995,7 @@ TEST(BoutMeshTest, GlobalXIntSymmetricX) { EXPECT_EQ(mesh01.GlobalX(4), 0.875); } -TEST(BoutMeshTest, GlobalXIntAsymmetricX) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, GlobalXIntAsymmetricX) { BoutMeshExposer mesh01(4, 3, 1, 2, 2, 0, 1, false, false); EXPECT_EQ(mesh01.GlobalX(0), 0.); EXPECT_EQ(mesh01.GlobalX(1), 0.25); @@ -851,10 +1004,7 @@ TEST(BoutMeshTest, GlobalXIntAsymmetricX) { EXPECT_EQ(mesh01.GlobalX(4), 1.0); } -TEST(BoutMeshTest, GlobalXRealSymmetricX) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, GlobalXRealSymmetricX) { BoutMeshExposer mesh01(4, 3, 1, 2, 2, 0, 1); EXPECT_EQ(mesh01.GlobalX(0.5), 0.); EXPECT_EQ(mesh01.GlobalX(1.5), 0.25); @@ -863,10 +1013,7 @@ TEST(BoutMeshTest, GlobalXRealSymmetricX) { EXPECT_EQ(mesh01.GlobalX(4.5), 1.0); } -TEST(BoutMeshTest, GlobalXRealAsymmetricX) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, GlobalXRealAsymmetricX) { BoutMeshExposer mesh01(4, 3, 1, 2, 2, 0, 1, false, false); EXPECT_EQ(mesh01.GlobalX(0.5), 0.125); EXPECT_EQ(mesh01.GlobalX(1.5), 0.375); @@ -875,9 +1022,7 @@ TEST(BoutMeshTest, GlobalXRealAsymmetricX) { EXPECT_EQ(mesh01.GlobalX(4.5), 1.125); } -TEST(BoutMeshTest, GetLocalXIndex) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, GetLocalXIndex) { // 2x2 processors, 3x3x1 (not including guards) on each processor // Boundaries are included in the local index @@ -921,9 +1066,7 @@ TEST(BoutMeshTest, GetLocalXIndex) { EXPECT_EQ(mesh11.getLocalXIndex(7), 4); } -TEST(BoutMeshTest, GetLocalXIndexNoBoundaries) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, GetLocalXIndexNoBoundaries) { // 2x2 processors, 3x3x1 (not including guards) on each processor // Local indices start counting from the first non-boundary point @@ -967,9 +1110,7 @@ TEST(BoutMeshTest, GetLocalXIndexNoBoundaries) { EXPECT_EQ(mesh11.getLocalXIndexNoBoundaries(6), 4); } -TEST(BoutMeshTest, GetGlobalYIndexSingleNull) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, GetGlobalYIndexSingleNull) { // 2x2 processors, 3x3x1 (not including guards) on each processor // Boundaries are included in the global index @@ -1013,9 +1154,7 @@ TEST(BoutMeshTest, GetGlobalYIndexSingleNull) { EXPECT_EQ(mesh11.getGlobalYIndex(4), 7); } -TEST(BoutMeshTest, GetGlobalYIndexDoubleNull) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, GetGlobalYIndexDoubleNull) { // 2x2 processors, 3x3x1 (not including guards) on each processor // Boundaries are included in the global index @@ -1064,9 +1203,7 @@ TEST(BoutMeshTest, GetGlobalYIndexDoubleNull) { EXPECT_EQ(mesh11.getGlobalYIndex(4), 9); } -TEST(BoutMeshTest, GetGlobalYIndexNoBoundaries) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, GetGlobalYIndexNoBoundaries) { // 2x2 processors, 3x3x1 (not including guards) on each processor // Global indices start counting from the first non-boundary point @@ -1110,9 +1247,7 @@ TEST(BoutMeshTest, GetGlobalYIndexNoBoundaries) { EXPECT_EQ(mesh11.getGlobalYIndexNoBoundaries(4), 6); } -TEST(BoutMeshTest, GetLocalYIndexSingleNull) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, GetLocalYIndexSingleNull) { // 2x2 processors, 3x3x1 (not including guards) on each processor // Boundaries are included in the local index @@ -1156,9 +1291,7 @@ TEST(BoutMeshTest, GetLocalYIndexSingleNull) { EXPECT_EQ(mesh11.getLocalYIndex(7), 4); } -TEST(BoutMeshTest, GetLocalYIndexDoubleNull) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, GetLocalYIndexDoubleNull) { // 2x2 processors, 3x3x1 (not including guards) on each processor // Boundaries are included in the global index @@ -1207,9 +1340,7 @@ TEST(BoutMeshTest, GetLocalYIndexDoubleNull) { EXPECT_EQ(mesh11.getLocalYIndex(9), 4); } -TEST(BoutMeshTest, GetLocalYIndexNoBoundaries) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, GetLocalYIndexNoBoundaries) { // 2x2 processors, 3x3x1 (not including guards) on each processor // Local indices start counting from the first non-boundary point @@ -1253,10 +1384,7 @@ TEST(BoutMeshTest, GetLocalYIndexNoBoundaries) { EXPECT_EQ(mesh11.getLocalYIndexNoBoundaries(6), 4); } -TEST(BoutMeshTest, GlobalYIntSymmetricY) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, GlobalYIntSymmetricY) { BoutMeshExposer mesh_inner_pf(createDisconnectedDoubleNull({12, 4, 1, 1, 1, 6, 0, 0})); EXPECT_EQ(mesh_inner_pf.GlobalY(0), -0.5625); EXPECT_EQ(mesh_inner_pf.GlobalY(1), -0.4375); @@ -1284,10 +1412,7 @@ TEST(BoutMeshTest, GlobalYIntSymmetricY) { EXPECT_EQ(mesh_outer_pf.GlobalY(3), 1.3125); } -TEST(BoutMeshTest, GlobalYIntAsymmetricY) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, GlobalYIntAsymmetricY) { auto grid_inner_pf = createDisconnectedDoubleNull({12, 4, 1, 1, 1, 6, 0, 0}); grid_inner_pf.grid.symmetric_Y = false; BoutMeshExposer mesh_inner_pf(grid_inner_pf); @@ -1321,10 +1446,7 @@ TEST(BoutMeshTest, GlobalYIntAsymmetricY) { EXPECT_EQ(mesh_outer_pf.GlobalY(3), 1); } -TEST(BoutMeshTest, GlobalYRealSymmetricY) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, GlobalYRealSymmetricY) { BoutMeshExposer mesh_inner_pf(createDisconnectedDoubleNull({12, 4, 1, 1, 1, 6, 0, 0})); EXPECT_EQ(mesh_inner_pf.GlobalY(0.5), -0.5); EXPECT_EQ(mesh_inner_pf.GlobalY(1.5), -0.375); @@ -1352,10 +1474,7 @@ TEST(BoutMeshTest, GlobalYRealSymmetricY) { EXPECT_EQ(mesh_outer_pf.GlobalY(3.5), 1.375); } -TEST(BoutMeshTest, GlobalYRealAsymmetricY) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, GlobalYRealAsymmetricY) { auto grid_inner_pf = createDisconnectedDoubleNull({12, 4, 1, 1, 1, 6, 0, 0}); grid_inner_pf.grid.symmetric_Y = false; BoutMeshExposer mesh_inner_pf(grid_inner_pf); @@ -1389,9 +1508,7 @@ TEST(BoutMeshTest, GlobalYRealAsymmetricY) { EXPECT_EQ(mesh_outer_pf.GlobalY(3.5), 1); } -TEST(BoutMeshTest, GetGlobalZIndex) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, GetGlobalZIndex) { // 2x2 processors, 3x3x1 (not including guards) on each processor // Boundaries are included in the global index @@ -1427,9 +1544,7 @@ TEST(BoutMeshTest, GetGlobalZIndex) { EXPECT_EQ(mesh11.getGlobalZIndex(4), 4); } -TEST(BoutMeshTest, GetGlobalZIndexNoBoundaries) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, GetGlobalZIndexNoBoundaries) { // 2x2 processors, 3x3x1 (not including guards) on each processor BoutMeshExposer mesh00(5, 3, 4, 2, 2, 0, 0); @@ -1461,9 +1576,7 @@ TEST(BoutMeshTest, GetGlobalZIndexNoBoundaries) { EXPECT_EQ(mesh11.getGlobalZIndexNoBoundaries(4), 4); } -TEST(BoutMeshTest, GetLocalZIndex) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, GetLocalZIndex) { // 2x2 processors, 3x3x1 (not including guards) on each processor // Boundaries are included in the local index @@ -1507,9 +1620,7 @@ TEST(BoutMeshTest, GetLocalZIndex) { EXPECT_EQ(mesh11.getLocalZIndex(4), 4); } -TEST(BoutMeshTest, GetLocalZIndexNoBoundaries) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; +TEST_F(BoutMeshTest, GetLocalZIndexNoBoundaries) { // 2x2 processors, 3x3x1 (not including guards) on each processor // Local indices start counting from the first non-boundary point @@ -1553,10 +1664,7 @@ TEST(BoutMeshTest, GetLocalZIndexNoBoundaries) { EXPECT_EQ(mesh11.getLocalZIndexNoBoundaries(4), 4); } -TEST(BoutMeshTest, FirstX) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, FirstX) { BoutMeshExposer mesh00(5, 3, 4, 3, 3, 0, 0); EXPECT_TRUE(mesh00.firstX()); BoutMeshExposer mesh10(5, 3, 4, 3, 3, 1, 0); @@ -1577,10 +1685,7 @@ TEST(BoutMeshTest, FirstX) { EXPECT_FALSE(mesh22.firstX()); } -TEST(BoutMeshTest, LastX) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, LastX) { BoutMeshExposer mesh00(5, 3, 4, 3, 3, 0, 0); EXPECT_FALSE(mesh00.lastX()); BoutMeshExposer mesh10(5, 3, 4, 3, 3, 1, 0); @@ -1601,10 +1706,7 @@ TEST(BoutMeshTest, LastX) { EXPECT_TRUE(mesh22.lastX()); } -TEST(BoutMeshTest, FirstY) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, FirstY) { BoutMeshExposer mesh00(5, 3, 4, 3, 3, 0, 0); EXPECT_TRUE(mesh00.firstY()); BoutMeshExposer mesh10(5, 3, 4, 3, 3, 1, 0); @@ -1625,10 +1727,7 @@ TEST(BoutMeshTest, FirstY) { EXPECT_FALSE(mesh22.firstY()); } -TEST(BoutMeshTest, LastY) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, LastY) { BoutMeshExposer mesh00(5, 3, 4, 3, 3, 0, 0); EXPECT_FALSE(mesh00.lastY()); BoutMeshExposer mesh10(5, 3, 4, 3, 3, 1, 0); @@ -1666,8 +1765,7 @@ void checkRegionSizes(const BoutMeshExposer& mesh, std::array rgn_lower_ // These next few tests check both default_connections and the Region // creation, as these are quite tightly linked. -TEST(BoutMeshTest, DefaultConnectionsCore1x1) { - WithQuietOutput info{output_info}; +TEST_F(BoutMeshTest, DefaultConnectionsCore1x1) { // 5x3x1 grid on 1 processor, 1 boundary point. Boundaries should be // simple 1D rectangles, with 4 boundaries on this processor BoutMeshExposer mesh00(5, 3, 1, 1, 1, 0, 0, false); @@ -1685,9 +1783,7 @@ TEST(BoutMeshTest, DefaultConnectionsCore1x1) { checkRegionSizes(mesh00, {5, 0, 5}, {0, 5, 5}, {3, 3}); } -TEST(BoutMeshTest, TopologySOL2x2) { - WithQuietOutput info{output_info}; - +TEST_F(BoutMeshTest, TopologySOL2x2) { { SCOPED_TRACE("TopologySOL2x2, mesh00"); BoutMeshExposer mesh00(createSOL({3, 3, 1, 1, 2, 2, 0, 0})); @@ -1725,9 +1821,7 @@ TEST(BoutMeshTest, TopologySOL2x2) { } } -TEST(BoutMeshTest, TopologySOLPeriodicX2x2) { - WithQuietOutput info{output_info}; - +TEST_F(BoutMeshTest, TopologySOLPeriodicX2x2) { { SCOPED_TRACE("TopologySOLPeriodicX2x2, mesh00"); @@ -1766,9 +1860,7 @@ TEST(BoutMeshTest, TopologySOLPeriodicX2x2) { } } -TEST(BoutMeshTest, TopologySingleNull2x3) { - WithQuietOutput info{output_info}; - +TEST_F(BoutMeshTest, TopologySingleNull2x3) { { SCOPED_TRACE("TopologySingleNull2x3, mesh00"); BoutMeshExposer mesh00(createSingleNull({3, 3, 1, 1, 2, 3, 0, 0})); @@ -1824,9 +1916,7 @@ TEST(BoutMeshTest, TopologySingleNull2x3) { } } -TEST(BoutMeshTest, TopologyDisconnectedDoubleNull1x6) { - WithQuietOutput info{output_info}; - +TEST_F(BoutMeshTest, TopologyDisconnectedDoubleNull1x6) { { SCOPED_TRACE("TopologyDisconnectedDoubleNull1x6, mesh00"); // Inner lower leg BoutMeshExposer mesh00(createDisconnectedDoubleNull({12, 3, 1, 1, 1, 6, 0, 0})); @@ -1882,8 +1972,7 @@ TEST(BoutMeshTest, TopologyDisconnectedDoubleNull1x6) { } } -TEST(BoutMeshTest, SetDerivedGridSizes) { - WithQuietOutput info{output_info}; +TEST_F(BoutMeshTest, SetDerivedGridSizes) { BoutMeshGridInfo grid{12, 3, 1, 2, 3, 6, 2, 2}; BoutMeshExposer mesh(createDisconnectedDoubleNull(grid)); @@ -1911,8 +2000,7 @@ TEST(BoutMeshTest, SetDerivedGridSizes) { EXPECT_EQ(mesh.zend, 0); } -TEST(BoutMeshTest, CreateXBoundariesPeriodicX) { - WithQuietOutput info{output_info}; +TEST_F(BoutMeshTest, CreateXBoundariesPeriodicX) { // Periodic in X, so no boundaries BoutMeshExposer mesh(createDisconnectedDoubleNull({12, 3, 1, 1, 3, 6, 1, 0})); mesh.periodicX = true; @@ -1922,8 +2010,7 @@ TEST(BoutMeshTest, CreateXBoundariesPeriodicX) { EXPECT_TRUE(boundaries.empty()); } -TEST(BoutMeshTest, CreateXBoundariesNoGuards) { - WithQuietOutput info{output_info}; +TEST_F(BoutMeshTest, CreateXBoundariesNoGuards) { // No guards in X, so no boundaries BoutMeshExposer mesh(createDisconnectedDoubleNull({12, 3, 0, 1, 3, 6, 1, 0})); mesh.createXBoundaries(); @@ -1932,8 +2019,7 @@ TEST(BoutMeshTest, CreateXBoundariesNoGuards) { EXPECT_TRUE(boundaries.empty()); } -TEST(BoutMeshTest, CreateXBoundariesDoubleNullInsidePF) { - WithQuietOutput info{output_info}; +TEST_F(BoutMeshTest, CreateXBoundariesDoubleNullInsidePF) { // Three cores in X, inside core, one boundary BoutMeshExposer mesh_inside(createDisconnectedDoubleNull({12, 3, 1, 1, 3, 6, 0, 0})); mesh_inside.createXBoundaries(); @@ -1943,8 +2029,7 @@ TEST(BoutMeshTest, CreateXBoundariesDoubleNullInsidePF) { EXPECT_EQ(boundaries_inside[0]->label, "pf"); } -TEST(BoutMeshTest, CreateXBoundariesDoubleNullMiddlePF) { - WithQuietOutput info{output_info}; +TEST_F(BoutMeshTest, CreateXBoundariesDoubleNullMiddlePF) { // Three cores in X, middle core, so no boundaries BoutMeshExposer mesh_middle(createDisconnectedDoubleNull({12, 3, 1, 1, 3, 6, 1, 0})); mesh_middle.createXBoundaries(); @@ -1953,8 +2038,7 @@ TEST(BoutMeshTest, CreateXBoundariesDoubleNullMiddlePF) { EXPECT_TRUE(boundaries_middle.empty()); } -TEST(BoutMeshTest, CreateXBoundariesDoubleNullOutsidePF) { - WithQuietOutput info{output_info}; +TEST_F(BoutMeshTest, CreateXBoundariesDoubleNullOutsidePF) { // Three cores in X, outside core, one boundary BoutMeshExposer mesh_inside(createDisconnectedDoubleNull({12, 3, 1, 1, 3, 6, 0, 0})); mesh_inside.createXBoundaries(); @@ -1964,8 +2048,7 @@ TEST(BoutMeshTest, CreateXBoundariesDoubleNullOutsidePF) { EXPECT_EQ(boundaries_inside[0]->label, "pf"); } -TEST(BoutMeshTest, CreateXBoundariesDoubleNullInsideOutsideCore) { - WithQuietOutput info{output_info}; +TEST_F(BoutMeshTest, CreateXBoundariesDoubleNullInsideOutsideCore) { // One core in X, so we expect two boundaries BoutMeshExposer mesh(createDisconnectedDoubleNull({12, 3, 1, 1, 1, 6, 0, 1})); mesh.createXBoundaries(); @@ -1976,9 +2059,8 @@ TEST(BoutMeshTest, CreateXBoundariesDoubleNullInsideOutsideCore) { EXPECT_EQ(boundaries[1]->label, "sol"); } -TEST(BoutMeshTest, CreateYBoundariesNoGuards) { - WithQuietOutput info{output_info}; +TEST_F(BoutMeshTest, CreateYBoundariesNoGuards) { BoutMeshExposer mesh(createDisconnectedDoubleNull({12, 3, 1, 0, 1, 6, 0, 0})); mesh.createYBoundaries(); @@ -1986,10 +2068,7 @@ TEST(BoutMeshTest, CreateYBoundariesNoGuards) { EXPECT_TRUE(boundaries.empty()); } -TEST(BoutMeshTest, CreateYBoundariesClosedFieldLines) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, CreateYBoundariesClosedFieldLines) { BoutMeshExposer mesh(createCore({4, 4, 2, 2, 4, 4})); mesh.createYBoundaries(); @@ -1997,9 +2076,7 @@ TEST(BoutMeshTest, CreateYBoundariesClosedFieldLines) { EXPECT_TRUE(boundaries.empty()); } -TEST(BoutMeshTest, CreateYBoundariesInnerLower) { - WithQuietOutput info{output_info}; - +TEST_F(BoutMeshTest, CreateYBoundariesInnerLower) { BoutMeshExposer mesh(createDisconnectedDoubleNull({12, 3, 1, 1, 1, 6, 0, 0})); mesh.createYBoundaries(); @@ -2008,9 +2085,7 @@ TEST(BoutMeshTest, CreateYBoundariesInnerLower) { EXPECT_EQ(boundaries[0]->label, "lower_target"); } -TEST(BoutMeshTest, CreateYBoundariesInnerUpper) { - WithQuietOutput info{output_info}; - +TEST_F(BoutMeshTest, CreateYBoundariesInnerUpper) { BoutMeshExposer mesh(createDisconnectedDoubleNull({12, 3, 1, 1, 1, 6, 0, 2})); mesh.createYBoundaries(); @@ -2019,9 +2094,7 @@ TEST(BoutMeshTest, CreateYBoundariesInnerUpper) { EXPECT_EQ(boundaries[0]->label, "upper_target"); } -TEST(BoutMeshTest, CreateYBoundariesOuterUpper) { - WithQuietOutput info{output_info}; - +TEST_F(BoutMeshTest, CreateYBoundariesOuterUpper) { BoutMeshExposer mesh(createDisconnectedDoubleNull({12, 3, 1, 1, 1, 6, 0, 5})); mesh.createYBoundaries(); @@ -2030,9 +2103,7 @@ TEST(BoutMeshTest, CreateYBoundariesOuterUpper) { EXPECT_EQ(boundaries[0]->label, "upper_target"); } -TEST(BoutMeshTest, CreateYBoundariesOuterLower) { - WithQuietOutput info{output_info}; - +TEST_F(BoutMeshTest, CreateYBoundariesOuterLower) { BoutMeshExposer mesh(createDisconnectedDoubleNull({12, 3, 1, 1, 1, 6, 0, 3})); mesh.createYBoundaries(); @@ -2041,9 +2112,7 @@ TEST(BoutMeshTest, CreateYBoundariesOuterLower) { EXPECT_EQ(boundaries[0]->label, "lower_target"); } -TEST(BoutMestTest, PeriodicY) { - WithQuietOutput info{output_info}; - +TEST_F(BoutMeshTest, PeriodicY) { BoutMeshExposer mesh00(createDisconnectedDoubleNull({12, 3, 1, 1, 1, 6, 0, 0})); EXPECT_FALSE(mesh00.periodicY(2)); EXPECT_FALSE(mesh00.periodicY(10)); @@ -2053,9 +2122,7 @@ TEST(BoutMestTest, PeriodicY) { EXPECT_FALSE(mesh01.periodicY(10)); } -TEST(BoutMestTest, PeriodicYWithShiftAngle) { - WithQuietOutput info{output_info}; - +TEST_F(BoutMeshTest, PeriodicYWithShiftAngle) { const std::vector shift_angle = {-1., 11., 10., 9., 8., 7., 6., 5., 4., 3., 2., 1., 0., -1.}; @@ -2076,9 +2143,7 @@ TEST(BoutMestTest, PeriodicYWithShiftAngle) { EXPECT_EQ(twist_shift01, 0.); } -TEST(BoutMeshTest, NumberOfYBoundaries) { - WithQuietOutput info{output_info}; - +TEST_F(BoutMeshTest, NumberOfYBoundaries) { BoutMeshExposer mesh_SOL(createSOL({3, 3, 1, 1, 2, 2, 1, 1})); EXPECT_EQ(mesh_SOL.numberOfYBoundaries(), 1); @@ -2086,9 +2151,7 @@ TEST(BoutMeshTest, NumberOfYBoundaries) { EXPECT_EQ(mesh_DND.numberOfYBoundaries(), 2); } -TEST(BoutMeshTest, HasBranchCutLower) { - WithQuietOutput info{output_info}; - +TEST_F(BoutMeshTest, HasBranchCutLower) { BoutMeshExposer mesh_SOL(createSOL({3, 3, 1, 1, 2, 2, 1, 1})); EXPECT_EQ(mesh_SOL.hasBranchCutLower(2), std::make_pair(false, 0.)); @@ -2103,9 +2166,7 @@ TEST(BoutMeshTest, HasBranchCutLower) { EXPECT_EQ(mesh_DND04.hasBranchCutLower(2), std::make_pair(false, 0.)); } -TEST(BoutMeshTest, HasBranchCutUpper) { - WithQuietOutput info{output_info}; - +TEST_F(BoutMeshTest, HasBranchCutUpper) { BoutMeshExposer mesh_SOL(createSOL({3, 3, 1, 1, 2, 2, 1, 1})); EXPECT_EQ(mesh_SOL.hasBranchCutUpper(2), std::make_pair(false, 0.)); @@ -2120,10 +2181,7 @@ TEST(BoutMeshTest, HasBranchCutUpper) { EXPECT_EQ(mesh_DND04.hasBranchCutUpper(2), std::make_pair(true, 10.)); } -TEST(BoutMeshTest, GetPossibleBoundariesCore) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, GetPossibleBoundariesCore) { BoutMeshExposer mesh_core_1x1(createCore({12, 3, 1, 1, 1, 1, 0, 0})); BoutMeshExposer mesh_core_32x64(createCore({12, 3, 1, 1, 32, 64, 7, 4})); @@ -2133,10 +2191,7 @@ TEST(BoutMeshTest, GetPossibleBoundariesCore) { EXPECT_EQ(mesh_core_32x64.getPossibleBoundaries(), boundaries); } -TEST(BoutMeshTest, GetPossibleBoundariesCorePeriodicX) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, GetPossibleBoundariesCorePeriodicX) { BoutMeshExposer mesh_core_1x1(createCore({12, 3, 1, 1, 1, 1, 0, 0}), true); BoutMeshExposer mesh_core_32x64(createCore({12, 3, 1, 1, 32, 64, 7, 4}), true); @@ -2144,10 +2199,7 @@ TEST(BoutMeshTest, GetPossibleBoundariesCorePeriodicX) { EXPECT_TRUE(mesh_core_32x64.getPossibleBoundaries().empty()); } -TEST(BoutMeshTest, GetPossibleBoundariesDND) { - WithQuietOutput info{output_info}; - WithQuietOutput warn{output_warn}; - +TEST_F(BoutMeshTest, GetPossibleBoundariesDND) { BoutMeshExposer mesh_DND_1x6(createDisconnectedDoubleNull({12, 3, 1, 1, 1, 6, 0, 1})); BoutMeshExposer mesh_DND_32x64( createDisconnectedDoubleNull({12, 3, 1, 1, 32, 64, 0, 4})); diff --git a/tests/unit/mesh/test_mesh.cxx b/tests/unit/mesh/test_mesh.cxx index 4d0c9111dd..0cb8beae45 100644 --- a/tests/unit/mesh/test_mesh.cxx +++ b/tests/unit/mesh/test_mesh.cxx @@ -288,7 +288,7 @@ TEST_F(MeshTest, MsgLen) { Field2D f2D_1(0., &localmesh); Field2D f2D_2(0., &localmesh); - std::vector var_list{&f3D_1, &f2D_1, &f3D_2, &f2D_2}; + const std::vector var_list{&f3D_1, &f2D_1, &f3D_2, &f2D_2}; const int len = localmesh.msg_len(var_list, 0, nx, 0, ny); diff --git a/tools/README.md b/tools/README.md index da8ff6c037..3cb109ad99 100644 --- a/tools/README.md +++ b/tools/README.md @@ -2,19 +2,5 @@ pylib -- Library of Python tools. This is the mainly used interface for post-processing. -idllib --- -Library of IDL tools - -matlablib --- -Library of Matlab tools. - -eigensolver --- -Eigenvalue solver for ITG test case - -pdb2idl --- -Library for reading PDB files into IDL - +Other tools where removed with v5.3.0, as they where unmaintained. +You can find them in older versions, e.g. in [5.2.0](https://github.com/boutproject/BOUT-dev/releases/tag/v5.2.0) diff --git a/tools/archiving/README b/tools/archiving/README deleted file mode 100644 index 3ccdb594b3..0000000000 --- a/tools/archiving/README +++ /dev/null @@ -1,8 +0,0 @@ -Tools or scripts to help manage & archive BOUT++ data output - -mdsplus Tools to manage BOUT++ runs using MDSplus - -pdb2cdf Convert PDB files to netCDF - -sdctools Simulation Data Compression library - diff --git a/tools/archiving/dumpsample/Makefile b/tools/archiving/dumpsample/Makefile deleted file mode 100644 index a9f3bebd90..0000000000 --- a/tools/archiving/dumpsample/Makefile +++ /dev/null @@ -1,29 +0,0 @@ - -CC = c++ -LD = c++ - -CFLAGS = -Wall -g - -PDB_PATH=/hwdisks/home/bd512/local/ - -INCLUDE = -I$(PDB_PATH)/include -LIBS = -lm -L$(PDB_PATH)/lib -lpdb -lpml -lscore - -TARGET = pdbsample -OBJ = pdbsample.o - -.PHONY:all -all: $(TARGET) - -$(TARGET): $(OBJ) Makefile - $(LD) -o $(TARGET) $(OBJ) $(LIBS) - -$(OBJ): %.o: %.cpp Makefile - $(CC) $(CFLAGS) -c $< -o $@ $(INCLUDE) - -.PHONY:clean -clean: - rm -f $(OBJ) $(TARGET) - -.PHONY:force -force: clean all diff --git a/tools/archiving/dumpsample/README b/tools/archiving/dumpsample/README deleted file mode 100644 index 70b4b64da5..0000000000 --- a/tools/archiving/dumpsample/README +++ /dev/null @@ -1,7 +0,0 @@ - -Codes: - -pdbsample Samples a single PDB file -ncsample Sample a netCDF file - -dumpsample Perl script to sample BOUT++ outputs diff --git a/tools/archiving/dumpsample/dumpsample b/tools/archiving/dumpsample/dumpsample deleted file mode 100755 index 942fd9731e..0000000000 --- a/tools/archiving/dumpsample/dumpsample +++ /dev/null @@ -1,51 +0,0 @@ -#!/usr/bin/perl -# -# Reduce the time-resolution of BOUT++ output. -# - -use File::Basename; - -# Check command-line arguments -$numArgs = $#ARGV + 1; - -if( $numArgs != 3 ) { - print "Useage: \n"; - exit; -} - -$inputdir = $ARGV[0]; -$outputdir = $ARGV[1]; -$ts = $ARGV[2]; - -if ($inputdir eq $outputdir) { - # Input and output directories are the same - print "Input and output directories must be different\n"; - exit; -} - -# Check 3rd argument is an integer -if($ARGV[2] =~ /^\d+$/) { - print "Reducing time resolution by a factor of $ARGV[2]\n"; -}else { - print "3rd argument must be an integer\n"; - exit; -} - -@infiles = <$inputdir/BOUT.dmp.*>; - -foreach $file (@infiles) { - $f = $outputdir . "/" . basename($file); - - if($file =~ m/.pdb$/) { - # PDB file - print "PDB file: $file -> $f\n"; - system("pdbsample $file $f $ts") == 0 - or die "Could not execute pdbsample\n"; - }elsif( $file =~ m/(.nc|.cdl|.cdf|.ncdf)$/ ) { - # NetCDF file - print "NetCDF file\n"; - }else { - print $file . " has unknown file format\n"; - } - -} diff --git a/tools/archiving/dumpsample/pdbsample.cxx b/tools/archiving/dumpsample/pdbsample.cxx deleted file mode 100644 index 90ea69adb2..0000000000 --- a/tools/archiving/dumpsample/pdbsample.cxx +++ /dev/null @@ -1,232 +0,0 @@ -/******************************************************* - * PDB Sample - * - * - * Ben Dudson, University of York, July 2009 - *******************************************************/ - -#include -#include - -// The PDB library (C) -#include "pdb.h" - -int main(int argc, char** argv) { - if (argc < 4) { - fprintf(stderr, "Useage: %s \n", argv[0]); - return 1; - } - - int tstride; - if (sscanf(argv[3], "%d", &tstride) != 1) { - fprintf(stderr, "\tERROR: t stride must be an integer\n"); - return 1; - } - - // Open input file - PDBfile* in; - if ((in = PD_open(argv[1], "r")) == NULL) { - fprintf(stderr, "\tERROR: Could not open input file '%s'\n", argv[1]); - return 1; - } - - // Open output file - PDBfile* out; - if ((out = PD_open(argv[2], "w")) == NULL) { - fprintf(stderr, "\tERROR: Could not open output file '%s'\n", argv[2]); - return 1; - } - - // Get list of variables - - int nvars; - char** var_names = PD_ls(in, NULL, NULL, &nvars); - if ((var_names == (char**)NULL) || (nvars < 1)) { - fprintf(stderr, "\tERROR: No variables\n"); - return 1; - } - - // Go through the variables - char* varname; - for (int var = 0; var < nvars; var++) { - varname = var_names[var]; - - syment* ep; // PDB query types - dimdes* dims; - - // Query size of the variable - if ((ep = PD_query_entry(in, varname, NULL)) == (syment*)NULL) { - fprintf(stderr, "\tError querying variable %s\n", varname); - return 1; - } - dims = PD_entry_dimensions(ep); - int nd = 0; // Count number of dimensions - int varsize = 1; // Number of elements - long inds[12]; - while (dims != (dimdes*)NULL) { - long min, max; - min = dims->index_min; - max = dims->index_max; - - if (nd > 3) { - fprintf(stderr, "\tERROR: Can't handle variable '%s': more than 4D\n", varname); - return 2; - } - - inds[3 * nd] = min; - inds[3 * nd + 1] = max; - inds[3 * nd + 2] = 1L; - - varsize *= max - min + 1; - - nd++; - dims = dims->next; - } - - // Get variable type - char* typ; - typ = PD_entry_type(ep); - - if ((strcmp(varname, "t_array") == 0) && (nd == 1)) { - float* fdata = new float[varsize]; - - // Read the data from the PDB file - - inds[2] = tstride; - int nread; - if ((nread = PD_read_as_alt(in, varname, "float", fdata, inds)) == 0) { - fprintf(stderr, "\tWARNING: Could not read t_array. Ignoring\n"); - delete[] fdata; - continue; - } - - inds[0] = 0L; - inds[1] = nread - 1; - inds[2] = 1L; - - if (PD_write_alt(out, varname, "float", fdata, 1, inds) == FALSE) { - fprintf(stderr, "\tWARNING: Could not write '%s'. Ignoring\n", varname); - } - - delete[] fdata; - - } else if (nd == 4) { - // Reducing time resolution - - if (strcasecmp(typ, "integer") == 0) { - int* idata = new int[varsize]; - - // Read the data from the PDB file - inds[2] = tstride; - int nread; - if ((nread = PD_read_as_alt(in, varname, "integer", idata, inds)) == 0) { - fprintf(stderr, "\tWARNING: Could not read '%s'. Ignoring\n", varname); - delete[] idata; - continue; - } - - if (nread % (varsize / (inds[1] - inds[0] + 1)) != 0) { - fprintf(stderr, "ERROR: Accounting error: (%ld, %ld), %d, %d\n", inds[0], - inds[1], varsize, nread); - delete[] idata; - continue; - } - - nread = nread / (varsize / (inds[1] - inds[0] + 1)); - - inds[0] = 0L; - inds[1] = nread - 1; - inds[2] = 1L; - - if (PD_write_alt(out, varname, "integer", idata, 4, inds) == FALSE) { - fprintf(stderr, "\tWARNING: Could not write '%s'. Ignoring\n", varname); - } - - delete[] idata; - - } else if ((strcasecmp(typ, "float") == 0) || (strcasecmp(typ, "double") == 0)) { - // Convert doubles to floats - - float* fdata = new float[varsize]; - - // Read the data from the PDB file - inds[2] = tstride; - int nread; - if ((nread = PD_read_as_alt(in, varname, "float", fdata, inds)) == 0) { - fprintf(stderr, "\tWARNING: Could not read '%s'. Ignoring\n", varname); - delete[] fdata; - continue; - } - - if (nread % (varsize / (inds[1] - inds[0] + 1)) != 0) { - fprintf(stderr, "ERROR: Accounting error: (%ld, %ld), %d, %d\n", inds[0], - inds[1], varsize, nread); - delete[] fdata; - continue; - } - - nread = nread / (varsize / (inds[1] - inds[0] + 1)); - - inds[0] = 0L; - inds[1] = nread - 1; - inds[2] = 1L; - - if (PD_write_alt(out, varname, "float", fdata, 4, inds) == FALSE) { - fprintf(stderr, "\tWARNING: Could not write '%s'. Ignoring\n", varname); - } - - delete[] fdata; - - } else { - fprintf(stderr, "\tWARNING: '%s' has unrecognised type '%s'. Ignoring\n", varname, - typ); - } - } else { - // Just copy the data across - - if (strcasecmp(typ, "integer") == 0) { - - int* idata = new int[varsize]; - - // Read the data from the PDB file - if (PD_read_as(in, varname, "integer", idata) == 0) { - fprintf(stderr, "\t\tWARNING: Could not read variable. Ignoring\n"); - delete[] idata; - continue; - } - - if (PD_write_alt(out, varname, "integer", idata, nd, inds) == FALSE) { - fprintf(stderr, "\tWARNING: Could not write variable '%s'\n", varname); - } - - delete[] idata; - - } else if ((strcasecmp(typ, "float") == 0) || (strcasecmp(typ, "double") == 0)) { - // Convert doubles to floats - - float* fdata = new float[varsize]; - - // Read the data from the PDB file - if (PD_read_as(in, varname, "float", fdata) == 0) { - fprintf(stderr, "\tWARNING: Could not read variable '%s'. Ignoring\n", varname); - delete[] fdata; - continue; - } - - if (PD_write_alt(out, varname, "float", fdata, nd, inds) == FALSE) { - fprintf(stderr, "\tWARNING: Could not write variable '%s'\n", varname); - } - - delete[] fdata; - } else { - fprintf(stderr, "WARNING: '%s' has unrecognised type '%s'. Ignoring\n", varname, - typ); - } - } - } - - PD_close(in); - PD_close(out); - - return 0; -} diff --git a/tools/archiving/mdsplus/data/bout_mds_model.characteristics b/tools/archiving/mdsplus/data/bout_mds_model.characteristics deleted file mode 100644 index b3a5920b84..0000000000 Binary files a/tools/archiving/mdsplus/data/bout_mds_model.characteristics and /dev/null differ diff --git a/tools/archiving/mdsplus/data/bout_mds_model.datafile b/tools/archiving/mdsplus/data/bout_mds_model.datafile deleted file mode 100644 index ef462ba99f..0000000000 Binary files a/tools/archiving/mdsplus/data/bout_mds_model.datafile and /dev/null differ diff --git a/tools/archiving/mdsplus/data/bout_mds_model.tree b/tools/archiving/mdsplus/data/bout_mds_model.tree deleted file mode 100644 index f4171bb4c8..0000000000 Binary files a/tools/archiving/mdsplus/data/bout_mds_model.tree and /dev/null differ diff --git a/tools/archiving/mdsplus/idl/mds2inp.pro b/tools/archiving/mdsplus/idl/mds2inp.pro deleted file mode 100644 index 8f6fdc4575..0000000000 --- a/tools/archiving/mdsplus/idl/mds2inp.pro +++ /dev/null @@ -1,177 +0,0 @@ -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -; -; Reads settings from an MDSplus shot, outputting a BOUT.inp file -; Currently just prints the settings to stdout -; -; KEYWORDS -; -; server Name (and optionally the port) of MDS server -; Default 'localhost' -; tree MDS tree. Default 'bout_mds' -; shot Shot number. By default reads the latest -; -; Run server: mdsip -p 8000 -m -h $MDS_ROOT/etc/mdsip.hosts -; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; - -PRO print_settings, path, nremove, prefix, fd - - setting_nid=MDSVALUE('GETNCI("'+path+'","NID_NUMBER")',status=status) - IF NOT (status AND 1) THEN RETURN - - n = N_ELEMENTS(setting_nid) - ;PRINT, "Number of settings:", n - - nvars = 0 ; keep track of the number of variables - FOR i=0, n-1 DO BEGIN - str = mdsvalue('GETNCI($,"MINPATH")',setting_nid[i]) - ;PRINT, str - - ; look for members - pos = STRPOS(str, ':') - - IF pos NE -1 THEN BEGIN - ; this string is a leaf - - len = STRLEN(str) - - name = STRMID(str, pos+1, len - pos-1) - path = STRMID(str, nremove, pos - nremove) - - ; replace '.' with '/' - STRREPLACE, path, '.', '/' - - ;PRINT, "-> '" + path + "' - '"+name+"'" - - IF nvars EQ 0 THEN BEGIN - paths = [path] - names = [name] - nids = [setting_nid[i]] - ENDIF ELSE BEGIN - paths = [paths, path] - names = [names, name] - nids = [nids, setting_nid[i]] - ENDELSE - nvars = nvars + 1 - - ENDIF - - ENDFOR - - ; get the settings without section - - IF STRLEN(prefix) GT 0 THEN PRINTF, fd, '['+prefix+']' - - w = WHERE(paths EQ '', count) - IF count GT 0 THEN BEGIN - FOR i=0, count-1 DO BEGIN - ; get the value - p = MDSVALUE('GETNCI($, "PATH")', nids[i]) - - v = MDSVALUE(p, status=status, /quiet) - - IF status AND 1 THEN BEGIN - ; success - PRINTF, fd, names[i]+" = " + STRING(v) - ENDIF ELSE BEGIN - ; failed - no data - - PRINTF, fd, '# '+names[i]+ ' = ' - ENDELSE - ENDFOR - - w = WHERE(paths NE '', count) - IF count GT 0 THEN BEGIN - paths = paths[w] - names = names[w] - nids = nids[w] - ENDIF - nvars = count - ENDIF - - WHILE nvars GT 0 DO BEGIN - IF STRLEN(prefix) GT 0 THEN PRINTF, fd, '['+prefix+'/'+paths[0]+']' ELSE PRINTF, fd, '['+paths[0]+']' - - w = WHERE(paths EQ paths[0], count) - - FOR i=0, count-1 DO BEGIN - ; get the value - p = MDSVALUE('GETNCI($, "PATH")', nids[w[i]]) - v = MDSVALUE(p, status=status, /quiet) - - IF status AND 1 THEN BEGIN - ; success - PRINTF, fd, names[i]+" = " + STRING(v) - ENDIF ELSE BEGIN - ; failed - no data - - PRINTF, fd, '# '+names[i]+ ' = ' - ENDELSE - ENDFOR - - w = WHERE(paths NE paths[0], count) - IF count GT 0 THEN BEGIN - paths = paths[w] - names = names[w] - nids = nids[w] - ENDIF - nvars = count - ENDWHILE -END - -PRO mds2inp, server=server, tree=tree, shot=shot - - IF NOT KEYWORD_SET(server) THEN server = "localhost" - IF NOT KEYWORD_SET(tree) THEN tree = "bout_mds" - - mdsconnect, server, status = status - - IF NOT (status AND 1) THEN BEGIN - PRINT, "Could not connect" - RETURN - ENDIF - - IF NOT KEYWORD_SET(shot) THEN BEGIN - ; get the latest shot - shot = MDSVALUE('current_shot("'+tree+'")') - PRINT, "Latest shot = ", shot - ENDIF - - mdsopen, tree, shot, status=status - IF NOT (status AND 1) THEN BEGIN - PRINT, "Could not open tree" - mdsdisconnect - RETURN - ENDIF - - ; get the description and model - - desc = MDSVALUE(':DESC', status=status) - IF NOT (status AND 1) THEN desc = 'No description' - model = MDSVALUE(':MODEL', status=status) - IF NOT (status AND 1) THEN model = 'NO MODEL' - - ; NOTE: This method could be replaced with calls for CHILDREN_NIDS - ; but that doesn't seem to work (says no data for node) - - fd = -1 - - PRINTF, fd, "# BOUT++ input file. Generated from MDSplus data" - PRINTF, fd, "# Date: " + SYSTIME() - PRINTF, fd, "# Shot: " + STRTRIM(STRING(shot),2) - PRINTF, fd, "# Description: " + desc - PRINTF, fd, "# Model: " + model - PRINTF, fd, "" - - settings_path = "\\TOP.SETTINGS***" - print_settings, settings_path, 10, '', fd - - print_settings, "\\TOP."+model+".SETTINGS***", 10+STRLEN(model), model, fd - - - mdsclose - mdsdisconnect - -END - - diff --git a/tools/archiving/mdsplus/manual/Makefile b/tools/archiving/mdsplus/manual/Makefile deleted file mode 100644 index df1e184ecc..0000000000 --- a/tools/archiving/mdsplus/manual/Makefile +++ /dev/null @@ -1,8 +0,0 @@ - -TARGET = bout_mdsplus - -$(TARGET).pdf: $(TARGET).tex - latexmk -pdf $(TARGET).tex -interaction=batchmode - -clean: - latexmk -C diff --git a/tools/archiving/mdsplus/manual/bout_mdsplus.tex b/tools/archiving/mdsplus/manual/bout_mdsplus.tex deleted file mode 100644 index 17852e30f2..0000000000 --- a/tools/archiving/mdsplus/manual/bout_mdsplus.tex +++ /dev/null @@ -1,50 +0,0 @@ -% Manual describing the MDSplus interface to BOUT++ - -\documentclass[12pt]{article} -\usepackage[nofoot]{geometry} -\usepackage{graphicx} -\usepackage{fancyhdr} - -%% Modify margins -\addtolength{\oddsidemargin}{-.25in} -\addtolength{\evensidemargin}{-.25in} -\addtolength{\textwidth}{0.5in} -\addtolength{\textheight}{0.25in} -%% SET HEADERS AND FOOTERS - -\pagestyle{fancy} -\fancyfoot{} -\renewcommand{\sectionmark}[1]{ % Lower case Section marker style - \markright{\thesection.\ #1}} -\fancyhead[LE,RO]{\bfseries\thepage} % Page number (boldface) in left on even - % pages and right on odd pages -\renewcommand{\headrulewidth}{0.3pt} - -\newcommand{\code}[1]{\texttt{#1}} -\newcommand{\file}[1]{\texttt{\bf #1}} - -%% commands for boxes with important notes -\newlength{\notewidth} -\addtolength{\notewidth}{\textwidth} -\addtolength{\notewidth}{-3.\parindent} -\newcommand{\note}[1]{ -\fbox{ -\begin{minipage}{\notewidth} -{\bf NOTE}: #1 -\end{minipage} -}} - -\begin{document} - -\title{MDSplus interface for BOUT++} -\author{B.Dudson, University of York} - -\maketitle - -\tableofcontents - -\section{Introduction} - - - -\end{document} diff --git a/tools/archiving/pdb2cdf/Makefile b/tools/archiving/pdb2cdf/Makefile deleted file mode 100644 index d08edc127c..0000000000 --- a/tools/archiving/pdb2cdf/Makefile +++ /dev/null @@ -1,14 +0,0 @@ -BOUT_TOP = ../../.. -include $(BOUT_TOP)/make.config - -TARGET = pdb2cdf -OBJ = pdb2cdf.o - -.PHONY:all -all: $(TARGET) - -$(TARGET): $(OBJ) Makefile - $(LD) -o $(TARGET) $(OBJ) $(EXTRA_LIBS) - -$(OBJ): %.o: %.cxx Makefile - $(CC) $(CFLAGS) -c $< -o $@ $(EXTRA_INCS) diff --git a/tools/archiving/pdb2cdf/README b/tools/archiving/pdb2cdf/README deleted file mode 100644 index a8ce636d6a..0000000000 --- a/tools/archiving/pdb2cdf/README +++ /dev/null @@ -1,25 +0,0 @@ -PDB2CDF -======= - -Convert Portable Data Binary (PDB) files to netCDF format - -PDB: https://wci.llnl.gov/codes/pact/ -netCDF: http://www.unidata.ucar.edu/software/netcdf/ - -Handles up to 4-D integers, floats and doubles. Tries to assign the correct -dimensions to the variables using the following procedure: - -1. Find the variable with the largest number of dimensions -2. Labels the dimensions of this variable 'X', 'Y', 'Z', 'T' - (always in that order), and records the index ranges -3. For each variable, match index ranges to get the dimension, - using the first one which matches. Dimensions are assumed - to always be in the same order: XZ, YT, XY etc. but NOT ZX, TZ etc. -4. Variables which cannot be matched to dimensions will currently - cause the code to fall over. Should be improved later. - -The ordering of dimensions was chosen for BOUT++ files, but -apart from this, the code is not specific to BOUT++. - -B.Dudson, April 2009 -University of York diff --git a/tools/archiving/pdb2cdf/pdb2cdf.cxx b/tools/archiving/pdb2cdf/pdb2cdf.cxx deleted file mode 100644 index ac8bc85de6..0000000000 --- a/tools/archiving/pdb2cdf/pdb2cdf.cxx +++ /dev/null @@ -1,382 +0,0 @@ -/******************************************************* - * PDB2CDF - * - * Convert PDB files to netCDF - *******************************************************/ - -#include -#include - -// The PDB library (C) -#include "pdb.h" - -// netCDF C++ library -#include - -// Dimension -struct TDim { - char* label; // Label for the dimension - int size; - int min, max; // Minimum, maximum index - NcDim* nDim; // netCDF dimension -}; - -// List of dimensions. Handles up to 3 -static TDim dimlist3d[] = {{"x", 0}, {"y", 0}, {"z", 0}}; - -// Special case for 4D -static TDim dimlist4d[] = {{"t", 0}, {"x", 0}, {"y", 0}, {"z", 0}}; - -int main(int argc, char** argv) { - TDim* dimlist; - - if (argc < 2) { - fprintf(stderr, "Useage: %s file1 file2 ...\n", argv[0]); - return 1; - } - - for (int i = 1; i < argc; i++) { // Go through each PDB file - char* inname = argv[i]; - PDBfile* in; - - // Open input file - - if ((in = PD_open(inname, "r")) == NULL) { - fprintf(stderr, "ERROR: Could not open input file '%s'\n", inname); - return 1; - } - - // Get list of variables - - int nvars; - char** var_names = PD_ls(in, NULL, NULL, &nvars); - if ((var_names == (char**)NULL) || (nvars < 1)) { - fprintf(stderr, "ERROR: No variables\n"); - return 1; - } - - // Create a filename for the output - char* outname; - int len = strlen(inname); - if (len < 5) { - // Not long enough - just append .cdl - outname = new char[len + 5]; - snprintf(outname, len + 5, "%s%s", inname, ".cdl"); - } else { - // Replace '.pdb' by '.cdl' - outname = new char[len + 1]; - strncpy(outname, inname, len + 1); - strncpy(outname + len - 3, "cdl", 3); - } - - // Open output file, overwriting if it exists - NcFile dataFile(outname, NcFile::Replace); - - if (!dataFile.is_valid()) { - fprintf(stderr, "ERROR: Could not open output file '%s'\n", outname); - return 1; - } - - printf("Converting %s -> %s\n", inname, outname); - - // Find the variable with the largest number of dimensions - printf("\tAnalysing dimensions..."); - int maxdims = -1; // maximum number of variables - int maxdimvar = 0; // variable with the most dimensions - - syment* ep; // PDB query types - dimdes* dims; - - char* varname; - for (int var = 0; var < nvars; var++) { - varname = var_names[var]; - // Query size of the variable - - if ((ep = PD_query_entry(in, varname, NULL)) == (syment*)NULL) { - fprintf(stderr, "Error querying variable %s\n", varname); - return 1; - } - dims = PD_entry_dimensions(ep); - int nd = 0; /* Count number of dimensions */ - while (dims != (dimdes*)NULL) { - nd++; - dims = dims->next; - } - if (nd > maxdims) { - maxdims = nd; - maxdimvar = var; - } - } - - printf("%d dimensions\n", maxdims); - - if (maxdims < 4) { - dimlist = dimlist3d; - } else { - dimlist = dimlist4d; - } - - if (maxdims > 4) { - fprintf(stderr, "ERROR: Can only handle up to 4D variables\n"); - return 1; - } - - // Get the size of each dimension - varname = var_names[maxdimvar]; - if ((ep = PD_query_entry(in, varname, NULL)) == (syment*)NULL) { - fprintf(stderr, "Error querying variable %s\n", varname); - return 1; - } - dims = PD_entry_dimensions(ep); - - for (int d = 0; d < maxdims; d++) { - dimlist[d].min = dims->index_min; - dimlist[d].max = dims->index_max; - dimlist[d].size = dims->index_max - dims->index_min + 1; - - // Create a netCDF dimension - dimlist[d].nDim = dataFile.add_dim(dimlist[d].label, dimlist[d].size); - - printf("\t\t%s: %d -> %d (%d)\n", dimlist[d].label, dimlist[d].min, dimlist[d].max, - dimlist[d].size); - - dims = dims->next; - } - - // Go through each variable - for (int var = 0; var < nvars; var++) { - varname = var_names[var]; - if ((ep = PD_query_entry(in, varname, NULL)) == (syment*)NULL) { - fprintf(stderr, "Error querying variable %s\n", varname); - return 1; - } - - printf("\t%s", varname); - - // Get dimensions - int nd = 0; // Number of dimensions - int vardim[4]; // Indices in the dimlist array - int lastdim = -1; // Always assume indices have the same order - int varsize = 1; // Number of elements - - bool gotdims = true; - - dims = PD_entry_dimensions(ep); - while (dims != (dimdes*)NULL) { - int min, max; - min = dims->index_min; - max = dims->index_max; - - varsize *= max - min + 1; - - vardim[nd] = -1; - // Find which dimension this corresponds to - for (int d = lastdim + 1; d < maxdims; d++) { - if ((dimlist[d].min == min) && (dimlist[d].max == max)) { - if (lastdim == -1) { - printf("[%s", dimlist[d].label); - } else { - printf(",%s", dimlist[d].label); - } - vardim[nd] = d; - lastdim = d; - break; - } - } - if (vardim[nd] == -1) { - // Not an existing dimension. Should create a new dimension - fprintf(stderr, "ERROR: %s has an unrecognised %d dimension [%d -> %d]\n", - varname, nd + 1, min, max); - gotdims = false; - break; - } - - nd++; - dims = dims->next; - } - - if (!gotdims) { - continue; // Skip this variable - } - - if (lastdim != -1) { - printf("] (%d elements) ", varsize); - } - - // Now know number of dimensions nd, and a list of dimension indices vardim - - // Get variable type - char* typ; - typ = PD_entry_type(ep); - - printf(" Type: %s ", typ); - fflush(stdout); - - if (strcasecmp(typ, "integer") == 0) { - - int* idata = new int[varsize]; - - // Read the data from the PDB file - if (PD_read_as(in, varname, "integer", idata) == FALSE) { - fprintf(stderr, "\t\tWARNING: Could not read variable. Ignoring\n"); - continue; - } - switch (nd) { - case 0: { - // Add a 0-D variable to the netCDF file - NcVar* ncdata = dataFile.add_var(varname, ncInt); - - // Write data - ncdata->put(idata); - break; - } - case 1: { - NcVar* ncdata = dataFile.add_var(varname, ncInt, dimlist[vardim[0]].nDim); - ncdata->put(idata, varsize); - break; - } - case 2: { - NcVar* ncdata = dataFile.add_var(varname, ncInt, dimlist[vardim[0]].nDim, - dimlist[vardim[1]].nDim); - ncdata->put(idata, dimlist[vardim[0]].size, dimlist[vardim[1]].size); - break; - } - case 3: { - NcVar* ncdata = - dataFile.add_var(varname, ncInt, dimlist[vardim[0]].nDim, - dimlist[vardim[1]].nDim, dimlist[vardim[2]].nDim); - ncdata->put(idata, dimlist[vardim[0]].size, dimlist[vardim[1]].size, - dimlist[vardim[2]].size); - break; - } - case 4: { - NcVar* ncdata = dataFile.add_var( - varname, ncInt, dimlist[vardim[0]].nDim, dimlist[vardim[1]].nDim, - dimlist[vardim[2]].nDim, dimlist[vardim[3]].nDim); - ncdata->put(idata, dimlist[vardim[0]].size, dimlist[vardim[1]].size, - dimlist[vardim[2]].size, dimlist[vardim[3]].size); - break; - } - default: { - fprintf(stderr, "\t\tWARNING: Cannot yet handle %d-D integers. Ignoring\n", nd); - } - } - delete[] idata; - - ////////////////////////////////////////////////////////////////// - - } else if (strcasecmp(typ, "float") == 0) { - - float* fdata = new float[varsize]; - - // Read the data from the PDB file - if (PD_read_as(in, varname, "float", fdata) == FALSE) { - fprintf(stderr, "\t\tWARNING: Could not read variable. Ignoring\n"); - continue; - } - switch (nd) { - case 0: { - NcVar* ncdata = dataFile.add_var(varname, ncFloat); - - // Write data - ncdata->put(fdata); - break; - } - case 1: { - NcVar* ncdata = dataFile.add_var(varname, ncFloat, dimlist[vardim[0]].nDim); - ncdata->put(fdata, varsize); - break; - } - case 2: { - NcVar* ncdata = dataFile.add_var(varname, ncFloat, dimlist[vardim[0]].nDim, - dimlist[vardim[1]].nDim); - ncdata->put(fdata, dimlist[vardim[0]].size, dimlist[vardim[1]].size); - break; - } - case 3: { - NcVar* ncdata = - dataFile.add_var(varname, ncFloat, dimlist[vardim[0]].nDim, - dimlist[vardim[1]].nDim, dimlist[vardim[2]].nDim); - ncdata->put(fdata, dimlist[vardim[0]].size, dimlist[vardim[1]].size, - dimlist[vardim[2]].size); - break; - } - case 4: { - NcVar* ncdata = dataFile.add_var( - varname, ncFloat, dimlist[vardim[0]].nDim, dimlist[vardim[1]].nDim, - dimlist[vardim[2]].nDim, dimlist[vardim[3]].nDim); - ncdata->put(fdata, dimlist[vardim[0]].size, dimlist[vardim[1]].size, - dimlist[vardim[2]].size, dimlist[vardim[3]].size); - break; - } - default: { - fprintf(stderr, "\t\tWARNING: Cannot yet handle %d-D floats. Ignoring\n", nd); - } - } - delete[] fdata; - - ////////////////////////////////////////////////////////////////// - - } else if (strcasecmp(typ, "double") == 0) { - - double* ddata = new double[varsize]; - - // Read the data from the PDB file - if (PD_read_as(in, varname, "double", ddata) == FALSE) { - fprintf(stderr, "\t\tWARNING: Could not read variable. Ignoring\n"); - continue; - } - switch (nd) { - case 0: { - NcVar* ncdata = dataFile.add_var(varname, ncDouble); - - // Write data - ncdata->put(ddata); - break; - } - case 1: { - NcVar* ncdata = dataFile.add_var(varname, ncDouble, dimlist[vardim[0]].nDim); - ncdata->put(ddata, varsize); - break; - } - case 2: { - NcVar* ncdata = dataFile.add_var(varname, ncDouble, dimlist[vardim[0]].nDim, - dimlist[vardim[1]].nDim); - ncdata->put(ddata, dimlist[vardim[0]].size, dimlist[vardim[1]].size); - break; - } - case 3: { - NcVar* ncdata = - dataFile.add_var(varname, ncDouble, dimlist[vardim[0]].nDim, - dimlist[vardim[1]].nDim, dimlist[vardim[2]].nDim); - ncdata->put(ddata, dimlist[vardim[0]].size, dimlist[vardim[1]].size, - dimlist[vardim[2]].size); - break; - } - case 4: { - NcVar* ncdata = dataFile.add_var( - varname, ncDouble, dimlist[vardim[0]].nDim, dimlist[vardim[1]].nDim, - dimlist[vardim[2]].nDim, dimlist[vardim[3]].nDim); - ncdata->put(ddata, dimlist[vardim[0]].size, dimlist[vardim[1]].size, - dimlist[vardim[2]].size, dimlist[vardim[3]].size); - break; - } - default: { - fprintf(stderr, "\t\tWARNING: Cannot yet handle %d-D doubles. Ignoring\n", nd); - } - } - delete[] ddata; - } else { - fprintf(stderr, "\tWARNING: '%s' has unrecognised type '%s'. Ignoring\n", varname, - typ); - } - printf("\n"); - } - - delete[] outname; - - dataFile.close(); // Close the output file. Probably optional. - } - - return 0; -} diff --git a/tools/archiving/sdctools/LICENSE b/tools/archiving/sdctools/LICENSE deleted file mode 100644 index a44189fdca..0000000000 --- a/tools/archiving/sdctools/LICENSE +++ /dev/null @@ -1,18 +0,0 @@ -Copyright (c) 2007 B.Dudson, University of York - -Permission is hereby granted, free of charge, to any person obtaining -a copy of this software and associated documentation files (the "Software"), -to deal in the Software without restriction, including without limitation -the rights to use, copy, modify, merge, publish, distribute, sublicense, -and/or sell copies of the Software, and to permit persons to whom the Software -is furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, -INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A -PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT -HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION -OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE -SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/tools/archiving/sdctools/Makefile b/tools/archiving/sdctools/Makefile deleted file mode 100644 index d501067023..0000000000 --- a/tools/archiving/sdctools/Makefile +++ /dev/null @@ -1,35 +0,0 @@ - -SDC = ./sdclib/ -PACT = /hwdisks/home/bd512/local/ -IDL_DIR = /hwdisks/sfw/idl/idl6.3/idl/ - -CC = gcc -LD = gcc - -CFLAGS = -Wall -O - -INCLUDES = -I$(SDC) -I$(PACT)/include -LIBS = $(SDC)/sdclib.o -L$(PACT)/lib -lpdb -lpml -lscore -lm - - -.PHONY: all - -all: compress sdc2idl.so - -compress: compress.o - $(LD) -o compress compress.o $(LIBS) - -compress.o: compress.c - $(CC) $(CFLAGS) -c compress.c $(INCLUDES) - - -sdc2idl.so: sdc2idl.o - $(LD) -o sdc2idl.so sdc2idl.o $(SDC)/sdclib.so -shared -lm - -sdc2idl.o: sdc2idl.c - $(CC) -c $(CFLAGS) -fPIC -c sdc2idl.c -I$(SDC) -I$(IDL_DIR)/external - - -.PHONY: clean -clean: - rm *.o compress sdc2idl.so diff --git a/tools/archiving/sdctools/README b/tools/archiving/sdctools/README deleted file mode 100644 index af4c37dbb2..0000000000 --- a/tools/archiving/sdctools/README +++ /dev/null @@ -1,207 +0,0 @@ - - -Scientific Data Compression Library (SDCLib) Version 1.0 -Ben Dudson, University of York, August 2007 - -Introduction -============ - -This is a simple compression library, intended for the relatively -smoothly varying output from simulations such as BOUT++. -Uses polynomial interpolation in time to reduce number of stored -time-points - If the input is a polynomial in time, then the minimum -number of points needed will be stored, for example to store a -quadratic 3 points will be stored. -Compression is lossy - some accuracy is sacrificed in order to -reduce the size of the data, although the level of accuracy can -be specified and changed mid-way through the data. -This accuracy is guaranteed - the reconstruction is guaranteed -to be within the absolute and relative error bounds specified. - -NOTES: -1. This library is currently under development and probably still has some bugs. - Please don't use it to store anything you cannot replace. - -2. Currently doesn't take account of different integer/float sizes - (e.g. 32 vs 64-bit machines) apart from checking that variables - are the correct size for the archive. - If you attempt to read using 64-bit code an archive created by 32-bit code, - an error will be thrown. - -3. Archives cannot be altered once created. Appending data to files may happen - eventually, but modifying files is just too much of a pain. - -Any comments/suggestions/flames to Ben Dudson, bd512@york.ac.uk - - -Operation -========= - -Because this is a compressed format, searching for a given time-point would be difficult -unless the entire dataset were periodically output in one block. -As with a video file, these are called i-frames, and a table of their locations is stored -at the end of the file (like AVI). - -In between these i-frames, different parts of the domain will (probably) change at -different rates, and so the domain is split into several regions - when creating an -archive the number of regions must be specified. - -Creating an SDC archive ------------------------ - -SDC stores data in blocks. The idea is that a file format would be customised for a -given application, with it's own header describing the data. -After this would come the SDC header and then the data. - -First you must open a file - -FILE *fp -fp = fopen("somefile.sdc", "wb"); - -Write whatever header you want to this file, then call - -SDCfile* sdc_newfile(FILE *fp, int N, int order, int reset, int nregions); - -which returns NULL if an error occurs. -N is the number of data points in each time-slice -order is the maximum order of the interpolation (1 - linear, 2 - quadratic etc) -reset is the number of time-points between i-frames where all data is output (like a video, for faster searching) -nregions is the number of regions to split the data into. This doesn't have to be a divisor of N. - -This will write the SDC header to the file, allocate memory and prepare to write data. - -To specify the accuracy to be used, call - -void sdc_set_tol(SDCfile *f, float abstol, float reltol, float eta); - -which sets the absolute tolerance abstol, the relative tolerance reltol, -and a small quantity used in the relative tolerance eta. -The criterion for writing a block of data is to make sure the error e never satisfies: - -( |e| > reltol ) or ( |e| / (|data| + eta) > reltol ) - -Hence values below eta are ignored for the relative tolerance calculation. -These tolerances can be changed whilst writing the file to adapt the accuracy if needed. -Currently the default values for these are abstol = 1.0e-3, reltol = 1.0e-2, -and eta = 1.0e-10 which seem to give pretty decent results. - -To write a time-slice of the data, call - -int sdc_write(SDCfile *f, float *data); - -where data is a pointer to an array of N elements. This will return non-zero -if an error occurs. When all the data has been written, YOU MUST call: - -int sdc_close(SDCfile *f); - -which will write a final i-frame and generally tidy up. If you fail to do this, -terrible things will happen, and your data may be corrupted too. -After this, you can do what you want to the file, then close it - -fclose(fp); - -Reading an archive ------------------- - -First open the file, read whatever headers you wrote and then call - -SDCfile* sdc_open(FILE *fp); - -which will read the SDC data header. To read a time-slice, call - -int sdc_read(SDCfile *f, int t, float *data); - -where t is the time-index (counted from zero at the start of the file), -and data is a pointer to a block of memory. This function will return -non-zero if an error occurs. Reading is most efficient in order from -first to last time-point; reading in reverse order is the least efficient. - -Thats it! When finished reading, you can call - -int sdc_close(SDCfile *f); - -to free memory etc. As with the writing, you have to close the file itself separately. - -Tricks / Examples -================= - -Interleaved data ----------------- - -Several SDC data sets can be interleaved in a single file, for example with two sets - -FILE *fp -fp = fopen("somefile.sdc", "wb"); - -/* Create headers */ -SDCfile *f1, *f2; -f1 = sdc_newfile(fp, N1, order1, reset1, nregions1); -f2 = sdc_newfile(fp, N2, order2, reset2, nregions2); - -Can now write data to f1 or f2 in any order and then close both of them. -To read the data, you must read headers in the same order: - -f1 = sdc_open(fp); -f2 = sdc_open(fp); - -but data can then be read in any order from f1 or f2. - -Compress utility ----------------- - -The imaginatively named "compress" code takes a set of PDB files, selects a single variable -and concatenates all the data into a single output SDC file called "output.sdc". -The data can be any floating-point variable of 1,2,3 or 4 dimensions. The final index of the variable -is assumed to be time. The variable in each file needs to be the same size in each spatial dimension, but -not necessarily the same number of time-points. - -For example, to take all the collected BOUT data, and compress the density to an output file, - -compress ni_xyzt data*.pdb - -You will then be asked for the following parameters (decent values for full 4D 50x64x64x