diff --git a/docs/CMakeLists.txt b/docs/CMakeLists.txt index 26ba1cf23db6bebdd9bc7160535b4defbcad4ffb..97febbed0fda483e3b61ab5de6ba1dbb72f1eca6 100644 --- a/docs/CMakeLists.txt +++ b/docs/CMakeLists.txt @@ -14,12 +14,11 @@ if(WITH_DOCUMENTATION) # --------------------- Doxygen setup --------------------- find_package(Doxygen) - find_file(DOXY name doxypy.py PATH ENV{PATH}) - if(DOXY-NOTFOUND) - message(STATUS "Warning, doxypy seems to be missing on your system. You may not be able to properly generate the documentation.") - endif() + # Output path for doxygen documentation - set(DOXYGEN_OUTPUT ${DOC_ROOT_DIR}/doxygen) + if(NOT DOXYGEN_OUTPUT) + set(DOXYGEN_OUTPUT ${DOC_ROOT_DIR}/doxygen) + endif() # doxygen tags are required for doxygen-sphinx link. file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/build/tags) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/build/html) @@ -43,7 +42,7 @@ if(WITH_DOCUMENTATION) set(GENERATE_HTML YES) set(GENERATE_XML YES) - configure_file(${CMAKE_CURRENT_SOURCE_DIR}/config/hysop.doxyfile.in ${CMAKE_CURRENT_BINARY_DIR}/config/hysop.doxyfile) + configure_file(${CMAKE_CURRENT_SOURCE_DIR}/config/hysop.doxyfile.in ${DOXY_CONFIG}) # create a target to build doxygen documentation # i.e. to run 'make doxygen' @@ -149,4 +148,4 @@ if(WITH_DOCUMENTATION) # add_dependencies(html doxygen) # add_dependencies(html apidoc) -endif() \ No newline at end of file +endif() diff --git a/docs/config/hysop.doxyfile.in b/docs/config/hysop.doxyfile.in index 78907963763fed73166db0e9e9b68f613fa52774..f38f4ce8f3cf093ed21e8e9c1958ef44a3746cf9 100644 --- a/docs/config/hysop.doxyfile.in +++ b/docs/config/hysop.doxyfile.in @@ -1,11 +1,11 @@ -# Doxyfile 1.8.5 +# Doxyfile 1.8.11 # This file describes the settings to be used by the documentation system # doxygen (www.doxygen.org) for a project. -# +# # All text after a double hash (##) is considered a comment and is placed in # front of the TAG it is preceding. -# +# # All text after a single hash (#) is considered a comment and will be ignored. # The format is: # TAG = value [value, ...] @@ -46,21 +46,21 @@ PROJECT_NUMBER = @PACKAGE_VERSION@ PROJECT_BRIEF = "Particle Methods simulation on hybrid architectures" -# With the PROJECT_LOGO tag one can specify an logo or icon that is included in -# the documentation. The maximum height of the logo should not exceed 55 pixels -# and the maximum width should not exceed 200 pixels. Doxygen will copy the logo -# to the output directory. +# With the PROJECT_LOGO tag one can specify a logo or an icon that is included +# in the documentation. The maximum height of the logo should not exceed 55 +# pixels and the maximum width should not exceed 200 pixels. Doxygen will copy +# the logo to the output directory. -PROJECT_LOGO = +PROJECT_LOGO = # The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute) path # into which the generated documentation will be written. If a relative path is # entered, it will be relative to the location where doxygen was started. If # left blank the current directory will be used. -OUTPUT_DIRECTORY = @CMAKE_BINARY_DIR@/DoxygenGeneratedDoc +OUTPUT_DIRECTORY = @DOXYGEN_OUTPUT@ -# If the CREATE_SUBDIRS tag is set to YES, then doxygen will create 4096 sub- +# If the CREATE_SUBDIRS tag is set to YES then doxygen will create 4096 sub- # directories (in 2 levels) under the output directory of each output format and # will distribute the generated files over these directories. Enabling this # option can be useful when feeding doxygen a huge amount of source files, where @@ -70,29 +70,39 @@ OUTPUT_DIRECTORY = @CMAKE_BINARY_DIR@/DoxygenGeneratedDoc CREATE_SUBDIRS = YES +# If the ALLOW_UNICODE_NAMES tag is set to YES, doxygen will allow non-ASCII +# characters to appear in the names of generated files. If set to NO, non-ASCII +# characters will be escaped, for example _xE3_x81_x84 will be used for Unicode +# U+3044. +# The default value is: NO. + +ALLOW_UNICODE_NAMES = NO + # The OUTPUT_LANGUAGE tag is used to specify the language in which all # documentation generated by doxygen is written. Doxygen will use this # information to generate all constant output in the proper language. -# Possible values are: Afrikaans, Arabic, Brazilian, Catalan, Chinese, Chinese- -# Traditional, Croatian, Czech, Danish, Dutch, English, Esperanto, Farsi, -# Finnish, French, German, Greek, Hungarian, Italian, Japanese, Japanese-en, -# Korean, Korean-en, Latvian, Norwegian, Macedonian, Persian, Polish, -# Portuguese, Romanian, Russian, Serbian, Slovak, Slovene, Spanish, Swedish, -# Turkish, Ukrainian and Vietnamese. +# Possible values are: Afrikaans, Arabic, Armenian, Brazilian, Catalan, Chinese, +# Chinese-Traditional, Croatian, Czech, Danish, Dutch, English (United States), +# Esperanto, Farsi (Persian), Finnish, French, German, Greek, Hungarian, +# Indonesian, Italian, Japanese, Japanese-en (Japanese with English messages), +# Korean, Korean-en (Korean with English messages), Latvian, Lithuanian, +# Macedonian, Norwegian, Persian (Farsi), Polish, Portuguese, Romanian, Russian, +# Serbian, Serbian-Cyrillic, Slovak, Slovene, Spanish, Swedish, Turkish, +# Ukrainian and Vietnamese. # The default value is: English. OUTPUT_LANGUAGE = English -# If the BRIEF_MEMBER_DESC tag is set to YES doxygen will include brief member +# If the BRIEF_MEMBER_DESC tag is set to YES, doxygen will include brief member # descriptions after the members that are listed in the file and class # documentation (similar to Javadoc). Set to NO to disable this. # The default value is: YES. BRIEF_MEMBER_DESC = YES -# If the REPEAT_BRIEF tag is set to YES doxygen will prepend the brief +# If the REPEAT_BRIEF tag is set to YES, doxygen will prepend the brief # description of a member or function before the detailed description -# +# # Note: If both HIDE_UNDOC_MEMBERS and BRIEF_MEMBER_DESC are set to NO, the # brief descriptions will be completely suppressed. # The default value is: YES. @@ -108,7 +118,7 @@ REPEAT_BRIEF = YES # the entity):The $name class, The $name widget, The $name file, is, provides, # specifies, contains, represents, a, an and the. -ABBREVIATE_BRIEF = +ABBREVIATE_BRIEF = # If the ALWAYS_DETAILED_SEC and REPEAT_BRIEF tags are both set to YES then # doxygen will generate a detailed section even if there is only a brief @@ -125,7 +135,7 @@ ALWAYS_DETAILED_SEC = NO INLINE_INHERITED_MEMB = NO -# If the FULL_PATH_NAMES tag is set to YES doxygen will prepend the full path +# If the FULL_PATH_NAMES tag is set to YES, doxygen will prepend the full path # before files name in the file list and in the header files. If set to NO the # shortest path that makes the file name unique will be used # The default value is: YES. @@ -137,12 +147,12 @@ FULL_PATH_NAMES = NO # part of the path. The tag can be used to show relative paths in the file list. # If left blank the directory from which doxygen is run is used as the path to # strip. -# +# # Note that you can specify absolute paths here, but also relative paths, which # will be relative from the directory where doxygen is started. # This tag requires that the tag FULL_PATH_NAMES is set to YES. -STRIP_FROM_PATH = +STRIP_FROM_PATH = # The STRIP_FROM_INC_PATH tag can be used to strip a user-defined part of the # path mentioned in the documentation of a class, which tells the reader which @@ -151,7 +161,7 @@ STRIP_FROM_PATH = # specify the list of include paths that are normally passed to the compiler # using the -I flag. -STRIP_FROM_INC_PATH = +STRIP_FROM_INC_PATH = # If the SHORT_NAMES tag is set to YES, doxygen will generate much shorter (but # less readable) file names. This can be useful is your file systems doesn't @@ -182,7 +192,7 @@ QT_AUTOBRIEF = NO # a brief description. This used to be the default behavior. The new default is # to treat a multi-line C++ comment block as a detailed description. Set this # tag to YES if you prefer the old behavior instead. -# +# # Note that setting this tag to YES also means that rational rose comments are # not recognized any more. # The default value is: NO. @@ -195,9 +205,9 @@ MULTILINE_CPP_IS_BRIEF = NO INHERIT_DOCS = YES -# If the SEPARATE_MEMBER_PAGES tag is set to YES, then doxygen will produce a -# new page for each member. If set to NO, the documentation of a member will be -# part of the file/class/namespace that contains it. +# If the SEPARATE_MEMBER_PAGES tag is set to YES then doxygen will produce a new +# page for each member. If set to NO, the documentation of a member will be part +# of the file/class/namespace that contains it. # The default value is: NO. SEPARATE_MEMBER_PAGES = NO @@ -218,13 +228,13 @@ TAB_SIZE = 8 # "Side Effects:". You can put \n's in the value part of an alias to insert # newlines. -ALIASES = +ALIASES = # This tag can be used to specify a number of word-keyword mappings (TCL only). # A mapping has the form "name=value". For example adding "class=itcl::class" # will allow you to use the command class in the itcl::class meaning. -TCL_SUBST = +TCL_SUBST = # Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C sources # only. Doxygen will then generate output that is more tailored for C. For @@ -259,16 +269,19 @@ OPTIMIZE_OUTPUT_VHDL = NO # extension. Doxygen has a built-in mapping, but you can override or extend it # using this tag. The format is ext=language, where ext is a file extension, and # language is one of the parsers supported by doxygen: IDL, Java, Javascript, -# C#, C, C++, D, PHP, Objective-C, Python, Fortran, VHDL. For instance to make -# doxygen treat .inc files as Fortran files (default is PHP), and .f files as C -# (default is Fortran), use: inc=Fortran f=C. -# -# Note For files without extension you can use no_extension as a placeholder. -# +# C#, C, C++, D, PHP, Objective-C, Python, Fortran (fixed format Fortran: +# FortranFixed, free formatted Fortran: FortranFree, unknown formatted Fortran: +# Fortran. In the later case the parser tries to guess whether the code is fixed +# or free formatted code, this is the default for Fortran type files), VHDL. For +# instance to make doxygen treat .inc files as Fortran files (default is PHP), +# and .f files as C (default is Fortran), use: inc=Fortran f=C. +# +# Note: For files without extension you can use no_extension as a placeholder. +# # Note that for custom extensions you also need to set FILE_PATTERNS otherwise # the files are not read by doxygen. -EXTENSION_MAPPING = +EXTENSION_MAPPING = # If the MARKDOWN_SUPPORT tag is enabled then doxygen pre-processes all comments # according to the Markdown format, which allows for more readable @@ -282,8 +295,8 @@ MARKDOWN_SUPPORT = NO # When enabled doxygen tries to link words that correspond to documented # classes, or namespaces to their corresponding documentation. Such a link can -# be prevented in individual cases by by putting a % sign in front of the word -# or globally by setting AUTOLINK_SUPPORT to NO. +# be prevented in individual cases by putting a % sign in front of the word or +# globally by setting AUTOLINK_SUPPORT to NO. # The default value is: YES. AUTOLINK_SUPPORT = YES @@ -323,13 +336,20 @@ SIP_SUPPORT = NO IDL_PROPERTY_SUPPORT = YES # If member grouping is used in the documentation and the DISTRIBUTE_GROUP_DOC -# tag is set to YES, then doxygen will reuse the documentation of the first +# tag is set to YES then doxygen will reuse the documentation of the first # member in the group (if any) for the other members of the group. By default # all members of a group must be documented explicitly. # The default value is: NO. DISTRIBUTE_GROUP_DOC = NO +# If one adds a struct or class to a group and this option is enabled, then also +# any nested class or struct is added to the same group. By default this option +# is disabled and one has to add nested compounds explicitly via \ingroup. +# The default value is: NO. + +GROUP_NESTED_COMPOUNDS = NO + # Set the SUBGROUPING tag to YES to allow class member groups of the same type # (for instance a group of public functions) to be put as a subgroup of that # type (e.g. under the Public Functions section). Set it to NO to prevent @@ -343,7 +363,7 @@ SUBGROUPING = YES # are shown inside the group in which they are included (e.g. using \ingroup) # instead of on a separate page (for HTML and Man pages) or section (for LaTeX # and RTF). -# +# # Note that this feature does not work in combination with # SEPARATE_MEMBER_PAGES. # The default value is: NO. @@ -388,7 +408,7 @@ LOOKUP_CACHE_SIZE = 0 # Build related configuration options #--------------------------------------------------------------------------- -# If the EXTRACT_ALL tag is set to YES doxygen will assume all entities in +# If the EXTRACT_ALL tag is set to YES, doxygen will assume all entities in # documentation are documented, even if no documentation was available. Private # class members and static file members will be hidden unless the # EXTRACT_PRIVATE respectively EXTRACT_STATIC tags are set to YES. @@ -398,35 +418,35 @@ LOOKUP_CACHE_SIZE = 0 EXTRACT_ALL = NO -# If the EXTRACT_PRIVATE tag is set to YES all private members of a class will +# If the EXTRACT_PRIVATE tag is set to YES, all private members of a class will # be included in the documentation. # The default value is: NO. EXTRACT_PRIVATE = NO -# If the EXTRACT_PACKAGE tag is set to YES all members with package or internal +# If the EXTRACT_PACKAGE tag is set to YES, all members with package or internal # scope will be included in the documentation. # The default value is: NO. EXTRACT_PACKAGE = YES -# If the EXTRACT_STATIC tag is set to YES all static members of a file will be +# If the EXTRACT_STATIC tag is set to YES, all static members of a file will be # included in the documentation. # The default value is: NO. EXTRACT_STATIC = YES -# If the EXTRACT_LOCAL_CLASSES tag is set to YES classes (and structs) defined -# locally in source files will be included in the documentation. If set to NO +# If the EXTRACT_LOCAL_CLASSES tag is set to YES, classes (and structs) defined +# locally in source files will be included in the documentation. If set to NO, # only classes defined in header files are included. Does not have any effect # for Java sources. # The default value is: YES. EXTRACT_LOCAL_CLASSES = YES -# This flag is only useful for Objective-C code. When set to YES local methods, +# This flag is only useful for Objective-C code. If set to YES, local methods, # which are defined in the implementation section but not in the interface are -# included in the documentation. If set to NO only methods in the interface are +# included in the documentation. If set to NO, only methods in the interface are # included. # The default value is: NO. @@ -451,21 +471,21 @@ HIDE_UNDOC_MEMBERS = YES # If the HIDE_UNDOC_CLASSES tag is set to YES, doxygen will hide all # undocumented classes that are normally visible in the class hierarchy. If set -# to NO these classes will be included in the various overviews. This option has -# no effect if EXTRACT_ALL is enabled. +# to NO, these classes will be included in the various overviews. This option +# has no effect if EXTRACT_ALL is enabled. # The default value is: NO. HIDE_UNDOC_CLASSES = YES # If the HIDE_FRIEND_COMPOUNDS tag is set to YES, doxygen will hide all friend -# (class|struct|union) declarations. If set to NO these declarations will be +# (class|struct|union) declarations. If set to NO, these declarations will be # included in the documentation. # The default value is: NO. HIDE_FRIEND_COMPOUNDS = NO # If the HIDE_IN_BODY_DOCS tag is set to YES, doxygen will hide any -# documentation blocks found inside the body of a function. If set to NO these +# documentation blocks found inside the body of a function. If set to NO, these # blocks will be appended to the function's detailed documentation block. # The default value is: NO. @@ -479,27 +499,41 @@ HIDE_IN_BODY_DOCS = NO INTERNAL_DOCS = NO # If the CASE_SENSE_NAMES tag is set to NO then doxygen will only generate file -# names in lower-case letters. If set to YES upper-case letters are also +# names in lower-case letters. If set to YES, upper-case letters are also # allowed. This is useful if you have classes or files whose names only differ # in case and if your file system supports case sensitive file names. Windows # and Mac users are advised to set this option to NO. # The default value is: system dependent. -CASE_SENSE_NAMES = NO +CASE_SENSE_NAMES = YES # If the HIDE_SCOPE_NAMES tag is set to NO then doxygen will show members with -# their full class and namespace scopes in the documentation. If set to YES the +# their full class and namespace scopes in the documentation. If set to YES, the # scope will be hidden. # The default value is: NO. HIDE_SCOPE_NAMES = NO +# If the HIDE_COMPOUND_REFERENCE tag is set to NO (default) then doxygen will +# append additional text to a page's title, such as Class Reference. If set to +# YES the compound reference will be hidden. +# The default value is: NO. + +HIDE_COMPOUND_REFERENCE= NO + # If the SHOW_INCLUDE_FILES tag is set to YES then doxygen will put a list of # the files that are included by a file in the documentation of that file. # The default value is: YES. SHOW_INCLUDE_FILES = YES +# If the SHOW_GROUPED_MEMB_INC tag is set to YES then Doxygen will add for each +# grouped member an include statement to the documentation, telling the reader +# which file to include in order to use the member. +# The default value is: NO. + +SHOW_GROUPED_MEMB_INC = NO + # If the FORCE_LOCAL_INCLUDES tag is set to YES then doxygen will list include # files with double quotes in the documentation rather than with sharp brackets. # The default value is: NO. @@ -514,14 +548,15 @@ INLINE_INFO = YES # If the SORT_MEMBER_DOCS tag is set to YES then doxygen will sort the # (detailed) documentation of file and class members alphabetically by member -# name. If set to NO the members will appear in declaration order. +# name. If set to NO, the members will appear in declaration order. # The default value is: YES. SORT_MEMBER_DOCS = NO # If the SORT_BRIEF_DOCS tag is set to YES then doxygen will sort the brief # descriptions of file, namespace and class members alphabetically by member -# name. If set to NO the members will appear in declaration order. +# name. If set to NO, the members will appear in declaration order. Note that +# this will also influence the order of the classes in the class list. # The default value is: NO. SORT_BRIEF_DOCS = NO @@ -565,27 +600,25 @@ SORT_BY_SCOPE_NAME = YES STRICT_PROTO_MATCHING = NO -# The GENERATE_TODOLIST tag can be used to enable ( YES) or disable ( NO) the -# todo list. This list is created by putting \todo commands in the -# documentation. +# The GENERATE_TODOLIST tag can be used to enable (YES) or disable (NO) the todo +# list. This list is created by putting \todo commands in the documentation. # The default value is: YES. GENERATE_TODOLIST = YES -# The GENERATE_TESTLIST tag can be used to enable ( YES) or disable ( NO) the -# test list. This list is created by putting \test commands in the -# documentation. +# The GENERATE_TESTLIST tag can be used to enable (YES) or disable (NO) the test +# list. This list is created by putting \test commands in the documentation. # The default value is: YES. GENERATE_TESTLIST = YES -# The GENERATE_BUGLIST tag can be used to enable ( YES) or disable ( NO) the bug +# The GENERATE_BUGLIST tag can be used to enable (YES) or disable (NO) the bug # list. This list is created by putting \bug commands in the documentation. # The default value is: YES. GENERATE_BUGLIST = YES -# The GENERATE_DEPRECATEDLIST tag can be used to enable ( YES) or disable ( NO) +# The GENERATE_DEPRECATEDLIST tag can be used to enable (YES) or disable (NO) # the deprecated list. This list is created by putting \deprecated commands in # the documentation. # The default value is: YES. @@ -596,7 +629,7 @@ GENERATE_DEPRECATEDLIST= YES # sections, marked by \if <section_label> ... \endif and \cond <section_label> # ... \endcond blocks. -ENABLED_SECTIONS = +ENABLED_SECTIONS = # The MAX_INITIALIZER_LINES tag determines the maximum number of lines that the # initial value of a variable or macro / define can have for it to appear in the @@ -610,8 +643,8 @@ ENABLED_SECTIONS = MAX_INITIALIZER_LINES = 29 # Set the SHOW_USED_FILES tag to NO to disable the list of files generated at -# the bottom of the documentation of classes and structs. If set to YES the list -# will mention the files that were used to generate the documentation. +# the bottom of the documentation of classes and structs. If set to YES, the +# list will mention the files that were used to generate the documentation. # The default value is: YES. SHOW_USED_FILES = YES @@ -638,7 +671,7 @@ SHOW_NAMESPACES = NO # by doxygen. Whatever the program writes to standard output is used as the file # version. For an example see the documentation. -FILE_VERSION_FILTER = +FILE_VERSION_FILTER = # The LAYOUT_FILE tag can be used to specify a layout file which will be parsed # by doxygen. The layout file controls the global structure of the generated @@ -646,12 +679,12 @@ FILE_VERSION_FILTER = # that represents doxygen's defaults, run doxygen with the -l option. You can # optionally specify a file name after the option, if omitted DoxygenLayout.xml # will be used as the name of the layout file. -# +# # Note that if you run doxygen from a directory containing a file called # DoxygenLayout.xml, doxygen will parse it automatically even if the LAYOUT_FILE # tag is left empty. -LAYOUT_FILE = +LAYOUT_FILE = # The CITE_BIB_FILES tag can be used to specify one or more bib files containing # the reference definitions. This must be a list of .bib files. The .bib @@ -659,10 +692,9 @@ LAYOUT_FILE = # to be installed. See also http://en.wikipedia.org/wiki/BibTeX for more info. # For LaTeX the style of the bibliography can be controlled using # LATEX_BIB_STYLE. To use this feature you need bibtex and perl available in the -# search path. Do not use file names with spaces, bibtex cannot handle them. See -# also \cite for info how to create references. +# search path. See also \cite for info how to create references. -CITE_BIB_FILES = +CITE_BIB_FILES = #--------------------------------------------------------------------------- # Configuration options related to warning and progress messages @@ -676,15 +708,15 @@ CITE_BIB_FILES = QUIET = NO # The WARNINGS tag can be used to turn on/off the warning messages that are -# generated to standard error ( stderr) by doxygen. If WARNINGS is set to YES +# generated to standard error (stderr) by doxygen. If WARNINGS is set to YES # this implies that the warnings are on. -# +# # Tip: Turn warnings on while writing the documentation. # The default value is: YES. WARNINGS = YES -# If the WARN_IF_UNDOCUMENTED tag is set to YES, then doxygen will generate +# If the WARN_IF_UNDOCUMENTED tag is set to YES then doxygen will generate # warnings for undocumented members. If EXTRACT_ALL is set to YES then this flag # will automatically be disabled. # The default value is: YES. @@ -701,12 +733,18 @@ WARN_IF_DOC_ERROR = YES # This WARN_NO_PARAMDOC option can be enabled to get warnings for functions that # are documented, but have no documentation for their parameters or return -# value. If set to NO doxygen will only warn about wrong or incomplete parameter -# documentation, but not about the absence of documentation. +# value. If set to NO, doxygen will only warn about wrong or incomplete +# parameter documentation, but not about the absence of documentation. # The default value is: NO. WARN_NO_PARAMDOC = YES +# If the WARN_AS_ERROR tag is set to YES then doxygen will immediately stop when +# a warning is encountered. +# The default value is: NO. + +WARN_AS_ERROR = NO + # The WARN_FORMAT tag determines the format of the warning messages that doxygen # can produce. The string should contain the $file, $line, and $text tags, which # will be replaced by the file and line number from which the warning originated @@ -721,7 +759,7 @@ WARN_FORMAT = "$file:$line: $text" # messages should be written. If left blank the output is written to standard # error (stderr). -WARN_LOGFILE = +WARN_LOGFILE = #--------------------------------------------------------------------------- # Configuration options related to the input files @@ -730,7 +768,7 @@ WARN_LOGFILE = # The INPUT tag is used to specify the files and/or directories that contain # documented source files. You may enter file names like myfile.cpp or # directories like /usr/src/myproject. Separate the files or directories with -# spaces. +# spaces. See also FILE_PATTERNS and EXTENSION_MAPPING # Note: If this tag is empty the current directory is searched. INPUT = @CMAKE_SOURCE_DIR@/hysop \ @@ -748,12 +786,17 @@ INPUT_ENCODING = UTF-8 # If the value of the INPUT tag contains directories, you can use the # FILE_PATTERNS tag to specify one or more wildcard patterns (like *.cpp and -# *.h) to filter out the source-files in the directories. If left blank the -# following patterns are tested:*.c, *.cc, *.cxx, *.cpp, *.c++, *.java, *.ii, -# *.ixx, *.ipp, *.i++, *.inl, *.idl, *.ddl, *.odl, *.h, *.hh, *.hxx, *.hpp, -# *.h++, *.cs, *.d, *.php, *.php4, *.php5, *.phtml, *.inc, *.m, *.markdown, -# *.md, *.mm, *.dox, *.py, *.f90, *.f, *.for, *.tcl, *.vhd, *.vhdl, *.ucf, -# *.qsf, *.as and *.js. +# *.h) to filter out the source-files in the directories. +# +# Note that for custom extensions or not directly supported extensions you also +# need to set EXTENSION_MAPPING for the extension otherwise the files are not +# read by doxygen. +# +# If left blank the following patterns are tested:*.c, *.cc, *.cxx, *.cpp, +# *.c++, *.java, *.ii, *.ixx, *.ipp, *.i++, *.inl, *.idl, *.ddl, *.odl, *.h, +# *.hh, *.hxx, *.hpp, *.h++, *.cs, *.d, *.php, *.php4, *.php5, *.phtml, *.inc, +# *.m, *.markdown, *.md, *.mm, *.dox, *.py, *.pyw, *.f90, *.f, *.for, *.tcl, +# *.vhd, *.vhdl, *.ucf, *.qsf, *.as and *.js. FILE_PATTERNS = *.doxygen \ *.py \ @@ -769,11 +812,11 @@ RECURSIVE = YES # The EXCLUDE tag can be used to specify files and/or directories that should be # excluded from the INPUT source files. This way you can easily exclude a # subdirectory from a directory tree whose root is specified with the INPUT tag. -# +# # Note that relative paths are relative to the directory from which doxygen is # run. -EXCLUDE = +EXCLUDE = # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or # directories that are symbolic links (a Unix file system feature) are excluded @@ -785,7 +828,7 @@ EXCLUDE_SYMLINKS = NO # If the value of the INPUT tag contains directories, you can use the # EXCLUDE_PATTERNS tag to specify one or more wildcard patterns to exclude # certain files from those directories. -# +# # Note that the wildcards are matched against the file with absolute path, so to # exclude all test directories for example use the pattern */test/* @@ -797,17 +840,17 @@ EXCLUDE_PATTERNS = */.svn/* \ # output. The symbol name can be a fully qualified name, a word, or if the # wildcard * is used, a substring. Examples: ANamespace, AClass, # AClass::ANamespace, ANamespace::*Test -# +# # Note that the wildcards are matched against the file with absolute path, so to # exclude all test directories use the pattern */test/* -EXCLUDE_SYMBOLS = +EXCLUDE_SYMBOLS = # The EXAMPLE_PATH tag can be used to specify one or more files or directories # that contain example code fragments that are included (see the \include # command). -EXAMPLE_PATH = +EXAMPLE_PATH = # If the value of the EXAMPLE_PATH tag contains directories, you can use the # EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp and @@ -827,24 +870,28 @@ EXAMPLE_RECURSIVE = NO # that contain images that are to be included in the documentation (see the # \image command). -IMAGE_PATH = +IMAGE_PATH = # The INPUT_FILTER tag can be used to specify a program that doxygen should # invoke to filter for each input file. Doxygen will invoke the filter program # by executing (via popen()) the command: -# +# # <filter> <input-file> -# +# # where <filter> is the value of the INPUT_FILTER tag, and <input-file> is the # name of an input file. Doxygen will then use the output that the filter # program writes to standard output. If FILTER_PATTERNS is specified, this tag # will be ignored. -# +# # Note that the filter must not add or remove lines; it is applied before the # code is scanned, but not when the output code is generated. If lines are added # or removed, the anchors will not be placed correctly. +# +# Note that for custom extensions or not directly supported extensions you also +# need to set EXTENSION_MAPPING for the extension otherwise the files are not +# properly processed by doxygen. -INPUT_FILTER = +INPUT_FILTER = # The FILTER_PATTERNS tag can be used to specify filters on a per file pattern # basis. Doxygen will compare the file name with each pattern and apply the @@ -852,11 +899,15 @@ INPUT_FILTER = # (like *.cpp=my_cpp_filter). See INPUT_FILTER for further information on how # filters are used. If the FILTER_PATTERNS tag is empty or if none of the # patterns match the file name, INPUT_FILTER is applied. +# +# Note that for custom extensions or not directly supported extensions you also +# need to set EXTENSION_MAPPING for the extension otherwise the files are not +# properly processed by doxygen. FILTER_PATTERNS = *.py=doxypy.py # If the FILTER_SOURCE_FILES tag is set to YES, the input filter (if set using -# INPUT_FILTER ) will also be used to filter the input files that are used for +# INPUT_FILTER) will also be used to filter the input files that are used for # producing the source files to browse (i.e. when SOURCE_BROWSER is set to YES). # The default value is: NO. @@ -868,14 +919,14 @@ FILTER_SOURCE_FILES = YES # *.ext= (so without naming a filter). # This tag requires that the tag FILTER_SOURCE_FILES is set to YES. -FILTER_SOURCE_PATTERNS = +FILTER_SOURCE_PATTERNS = # If the USE_MDFILE_AS_MAINPAGE tag refers to the name of a markdown file that # is part of the input, its contents will be placed on the main page # (index.html). This can be useful if you have a project on for instance GitHub # and want to reuse the introduction page also for the doxygen output. -USE_MDFILE_AS_MAINPAGE = +USE_MDFILE_AS_MAINPAGE = #--------------------------------------------------------------------------- # Configuration options related to source browsing @@ -883,7 +934,7 @@ USE_MDFILE_AS_MAINPAGE = # If the SOURCE_BROWSER tag is set to YES then a list of source files will be # generated. Documented entities will be cross-referenced with these sources. -# +# # Note: To get rid of all source code in the generated output, make sure that # also VERBATIM_HEADERS is set to NO. # The default value is: NO. @@ -916,7 +967,7 @@ REFERENCED_BY_RELATION = NO REFERENCES_RELATION = NO # If the REFERENCES_LINK_SOURCE tag is set to YES and SOURCE_BROWSER tag is set -# to YES, then the hyperlinks from functions in REFERENCES_RELATION and +# to YES then the hyperlinks from functions in REFERENCES_RELATION and # REFERENCED_BY_RELATION lists will link to the source code. Otherwise they will # link to the documentation. # The default value is: YES. @@ -938,16 +989,16 @@ SOURCE_TOOLTIPS = YES # source browser. The htags tool is part of GNU's global source tagging system # (see http://www.gnu.org/software/global/global.html). You will need version # 4.8.6 or higher. -# +# # To use it do the following: # - Install the latest version of global # - Enable SOURCE_BROWSER and USE_HTAGS in the config file # - Make sure the INPUT points to the root of the source tree # - Run doxygen as normal -# +# # Doxygen will invoke htags (and that will in turn invoke gtags), so these # tools must be available from the command line (i.e. in the search path). -# +# # The result: instead of the source browser generated by doxygen, the links to # source code will now point to the output of htags. # The default value is: NO. @@ -963,25 +1014,6 @@ USE_HTAGS = NO VERBATIM_HEADERS = YES -# If the CLANG_ASSISTED_PARSING tag is set to YES, then doxygen will use the -# clang parser (see: http://clang.llvm.org/) for more acurate parsing at the -# cost of reduced performance. This can be particularly helpful with template -# rich C++ code for which doxygen's built-in parser lacks the necessary type -# information. -# Note: The availability of this option depends on whether or not doxygen was -# compiled with the --with-libclang option. -# The default value is: NO. - -CLANG_ASSISTED_PARSING = NO - -# If clang assisted parsing is enabled you can provide the compiler with command -# line options that you would normally use when invoking the compiler. Note that -# the include paths will already be set by doxygen for the files and directories -# specified with INPUT and INCLUDE_PATH. -# This tag requires that the tag CLANG_ASSISTED_PARSING is set to YES. - -CLANG_OPTIONS = - #--------------------------------------------------------------------------- # Configuration options related to the alphabetical class index #--------------------------------------------------------------------------- @@ -1006,13 +1038,13 @@ COLS_IN_ALPHA_INDEX = 5 # while generating the index headers. # This tag requires that the tag ALPHABETICAL_INDEX is set to YES. -IGNORE_PREFIX = +IGNORE_PREFIX = #--------------------------------------------------------------------------- # Configuration options related to the HTML output #--------------------------------------------------------------------------- -# If the GENERATE_HTML tag is set to YES doxygen will generate HTML output +# If the GENERATE_HTML tag is set to YES, doxygen will generate HTML output # The default value is: YES. GENERATE_HTML = YES @@ -1035,7 +1067,7 @@ HTML_FILE_EXTENSION = .html # The HTML_HEADER tag can be used to specify a user-defined HTML header file for # each generated HTML page. If the tag is left blank doxygen will generate a # standard header. -# +# # To get valid HTML the header file that includes any scripts and style sheets # that doxygen needs, which is dependent on the configuration options used (e.g. # the setting GENERATE_TREEVIEW). It is highly recommended to start with a @@ -1050,7 +1082,7 @@ HTML_FILE_EXTENSION = .html # of the possible markers and block names see the documentation. # This tag requires that the tag GENERATE_HTML is set to YES. -HTML_HEADER = +HTML_HEADER = # The HTML_FOOTER tag can be used to specify a user-defined HTML footer for each # generated HTML page. If the tag is left blank doxygen will generate a standard @@ -1060,7 +1092,7 @@ HTML_HEADER = # that doxygen normally uses. # This tag requires that the tag GENERATE_HTML is set to YES. -HTML_FOOTER = +HTML_FOOTER = # The HTML_STYLESHEET tag can be used to specify a user-defined cascading style # sheet that is used by each HTML page. It can be used to fine-tune the look of @@ -1072,18 +1104,20 @@ HTML_FOOTER = # obsolete. # This tag requires that the tag GENERATE_HTML is set to YES. -HTML_STYLESHEET = +HTML_STYLESHEET = -# The HTML_EXTRA_STYLESHEET tag can be used to specify an additional user- -# defined cascading style sheet that is included after the standard style sheets +# The HTML_EXTRA_STYLESHEET tag can be used to specify additional user-defined +# cascading style sheets that are included after the standard style sheets # created by doxygen. Using this option one can overrule certain style aspects. # This is preferred over using HTML_STYLESHEET since it does not replace the -# standard style sheet and is therefor more robust against future updates. -# Doxygen will copy the style sheet file to the output directory. For an example -# see the documentation. +# standard style sheet and is therefore more robust against future updates. +# Doxygen will copy the style sheet files to the output directory. +# Note: The order of the extra style sheet files is of importance (e.g. the last +# style sheet in the list overrules the setting of the previous ones in the +# list). For an example see the documentation. # This tag requires that the tag GENERATE_HTML is set to YES. -HTML_EXTRA_STYLESHEET = +HTML_EXTRA_STYLESHEET = # The HTML_EXTRA_FILES tag can be used to specify one or more extra images or # other source files which should be copied to the HTML output directory. Note @@ -1093,10 +1127,10 @@ HTML_EXTRA_STYLESHEET = # files will be copied as-is; there are no commands or markers available. # This tag requires that the tag GENERATE_HTML is set to YES. -HTML_EXTRA_FILES = +HTML_EXTRA_FILES = # The HTML_COLORSTYLE_HUE tag controls the color of the HTML output. Doxygen -# will adjust the colors in the stylesheet and background images according to +# will adjust the colors in the style sheet and background images according to # this color. Hue is specified as an angle on a colorwheel, see # http://en.wikipedia.org/wiki/Hue for more information. For instance the value # 0 represents red, 60 is yellow, 120 is green, 180 is cyan, 240 is blue, 300 @@ -1127,8 +1161,9 @@ HTML_COLORSTYLE_GAMMA = 124 # If the HTML_TIMESTAMP tag is set to YES then the footer of each generated HTML # page will contain the date and time when the page was generated. Setting this -# to NO can help when comparing the output of multiple runs. -# The default value is: YES. +# to YES can help to show when doxygen was last run and thus if the +# documentation is up to date. +# The default value is: NO. # This tag requires that the tag GENERATE_HTML is set to YES. HTML_TIMESTAMP = YES @@ -1203,7 +1238,7 @@ DOCSET_PUBLISHER_NAME = Publisher # index.hhp is a project file that can be read by Microsoft's HTML Help Workshop # (see: http://www.microsoft.com/en-us/download/details.aspx?id=21138) on # Windows. -# +# # The HTML Help Workshop contains a compiler that can convert all HTML output # generated by doxygen into a single compiled HTML file (.chm). Compiled HTML # files are now used as the Windows 98 help format, and will replace the old @@ -1221,31 +1256,32 @@ GENERATE_HTMLHELP = NO # written to the html output directory. # This tag requires that the tag GENERATE_HTMLHELP is set to YES. -CHM_FILE = +CHM_FILE = # The HHC_LOCATION tag can be used to specify the location (absolute path -# including file name) of the HTML help compiler ( hhc.exe). If non-empty +# including file name) of the HTML help compiler (hhc.exe). If non-empty, # doxygen will try to run the HTML help compiler on the generated index.hhp. # The file has to be specified with full path. # This tag requires that the tag GENERATE_HTMLHELP is set to YES. -HHC_LOCATION = +HHC_LOCATION = -# The GENERATE_CHI flag controls if a separate .chi index file is generated ( -# YES) or that it should be included in the master .chm file ( NO). +# The GENERATE_CHI flag controls if a separate .chi index file is generated +# (YES) or that it should be included in the master .chm file (NO). # The default value is: NO. # This tag requires that the tag GENERATE_HTMLHELP is set to YES. GENERATE_CHI = NO -# The CHM_INDEX_ENCODING is used to encode HtmlHelp index ( hhk), content ( hhc) +# The CHM_INDEX_ENCODING is used to encode HtmlHelp index (hhk), content (hhc) # and project file content. # This tag requires that the tag GENERATE_HTMLHELP is set to YES. -CHM_INDEX_ENCODING = +CHM_INDEX_ENCODING = -# The BINARY_TOC flag controls whether a binary table of contents is generated ( -# YES) or a normal table of contents ( NO) in the .chm file. +# The BINARY_TOC flag controls whether a binary table of contents is generated +# (YES) or a normal table of contents (NO) in the .chm file. Furthermore it +# enables the Previous and Next buttons. # The default value is: NO. # This tag requires that the tag GENERATE_HTMLHELP is set to YES. @@ -1272,7 +1308,7 @@ GENERATE_QHP = NO # the HTML output folder. # This tag requires that the tag GENERATE_QHP is set to YES. -QCH_FILE = +QCH_FILE = # The QHP_NAMESPACE tag specifies the namespace to use when generating Qt Help # Project output. For more information please see Qt Help Project / Namespace @@ -1297,7 +1333,7 @@ QHP_VIRTUAL_FOLDER = doc # filters). # This tag requires that the tag GENERATE_QHP is set to YES. -QHP_CUST_FILTER_NAME = +QHP_CUST_FILTER_NAME = # The QHP_CUST_FILTER_ATTRS tag specifies the list of the attributes of the # custom filter to add. For more information please see Qt Help Project / Custom @@ -1305,21 +1341,21 @@ QHP_CUST_FILTER_NAME = # filters). # This tag requires that the tag GENERATE_QHP is set to YES. -QHP_CUST_FILTER_ATTRS = +QHP_CUST_FILTER_ATTRS = # The QHP_SECT_FILTER_ATTRS tag specifies the list of the attributes this # project's filter section matches. Qt Help Project / Filter Attributes (see: # http://qt-project.org/doc/qt-4.8/qthelpproject.html#filter-attributes). # This tag requires that the tag GENERATE_QHP is set to YES. -QHP_SECT_FILTER_ATTRS = +QHP_SECT_FILTER_ATTRS = # The QHG_LOCATION tag can be used to specify the location of Qt's # qhelpgenerator. If non-empty doxygen will try to run qhelpgenerator on the # generated .qhp file. # This tag requires that the tag GENERATE_QHP is set to YES. -QHG_LOCATION = +QHG_LOCATION = # If the GENERATE_ECLIPSEHELP tag is set to YES, additional index files will be # generated, together with the HTML files, they form an Eclipse help plugin. To @@ -1358,7 +1394,7 @@ DISABLE_INDEX = NO # index structure (just like the one that is generated for HTML Help). For this # to work a browser that supports JavaScript, DHTML, CSS and frames is required # (i.e. any modern browser). Windows users are probably better off using the -# HTML help feature. Via custom stylesheets (see HTML_EXTRA_STYLESHEET) one can +# HTML help feature. Via custom style sheets (see HTML_EXTRA_STYLESHEET) one can # further fine-tune the look of the index. As an example, the default style # sheet generated by doxygen has an example that shows how to put an image at # the root of the tree instead of the PROJECT_NAME. Since the tree basically has @@ -1371,7 +1407,7 @@ GENERATE_TREEVIEW = NO # The ENUM_VALUES_PER_LINE tag can be used to set the number of enum values that # doxygen will group on one line in the generated HTML documentation. -# +# # Note that a value of 0 will completely suppress the enum values from appearing # in the overview section. # Minimum value: 0, maximum value: 20, default value: 4. @@ -1386,7 +1422,7 @@ ENUM_VALUES_PER_LINE = 4 TREEVIEW_WIDTH = 250 -# When the EXT_LINKS_IN_WINDOW option is set to YES doxygen will open links to +# If the EXT_LINKS_IN_WINDOW option is set to YES, doxygen will open links to # external symbols imported via tag files in a separate window. # The default value is: NO. # This tag requires that the tag GENERATE_HTML is set to YES. @@ -1405,7 +1441,7 @@ FORMULA_FONTSIZE = 10 # Use the FORMULA_TRANPARENT tag to determine whether or not the images # generated for formulas are transparent PNGs. Transparent PNGs are not # supported properly for IE 6.0, but are supported on all modern browsers. -# +# # Note that when changing this option you need to delete any form_*.png files in # the HTML output directory before the changes have effect. # The default value is: YES. @@ -1415,7 +1451,7 @@ FORMULA_TRANSPARENT = YES # Enable the USE_MATHJAX option to render LaTeX formulas using MathJax (see # http://www.mathjax.org) which uses client side Javascript for the rendering -# instead of using prerendered bitmaps. Use this if you do not have LaTeX +# instead of using pre-rendered bitmaps. Use this if you do not have LaTeX # installed or if you want to formulas look prettier in the HTML output. When # enabled you may also need to install MathJax separately and configure the path # to it using the MATHJAX_RELPATH option. @@ -1452,7 +1488,7 @@ MATHJAX_RELPATH = http://www.mathjax.org/mathjax # MATHJAX_EXTENSIONS = TeX/AMSmath TeX/AMSsymbols # This tag requires that the tag USE_MATHJAX is set to YES. -MATHJAX_EXTENSIONS = +MATHJAX_EXTENSIONS = # The MATHJAX_CODEFILE tag can be used to specify a file with javascript pieces # of code that will be used on startup of the MathJax code. See the MathJax site @@ -1460,7 +1496,7 @@ MATHJAX_EXTENSIONS = # example see the documentation. # This tag requires that the tag USE_MATHJAX is set to YES. -MATHJAX_CODEFILE = +MATHJAX_CODEFILE = # When the SEARCHENGINE tag is enabled doxygen will generate a search box for # the HTML output. The underlying search engine uses javascript and DHTML and @@ -1485,11 +1521,11 @@ SEARCHENGINE = YES # When the SERVER_BASED_SEARCH tag is enabled the search engine will be # implemented using a web server instead of a web client using Javascript. There -# are two flavours of web server based searching depending on the -# EXTERNAL_SEARCH setting. When disabled, doxygen will generate a PHP script for -# searching and an index file used by the script. When EXTERNAL_SEARCH is -# enabled the indexing and searching needs to be provided by external tools. See -# the section "External Indexing and Searching" for details. +# are two flavors of web server based searching depending on the EXTERNAL_SEARCH +# setting. When disabled, doxygen will generate a PHP script for searching and +# an index file used by the script. When EXTERNAL_SEARCH is enabled the indexing +# and searching needs to be provided by external tools. See the section +# "External Indexing and Searching" for details. # The default value is: NO. # This tag requires that the tag SEARCHENGINE is set to YES. @@ -1500,11 +1536,11 @@ SERVER_BASED_SEARCH = NO # which needs to be processed by an external indexer. Doxygen will invoke an # external search engine pointed to by the SEARCHENGINE_URL option to obtain the # search results. -# -# Doxygen ships with an example indexer ( doxyindexer) and search engine +# +# Doxygen ships with an example indexer (doxyindexer) and search engine # (doxysearch.cgi) which are based on the open source search engine library # Xapian (see: http://xapian.org/). -# +# # See the section "External Indexing and Searching" for details. # The default value is: NO. # This tag requires that the tag SEARCHENGINE is set to YES. @@ -1513,14 +1549,14 @@ EXTERNAL_SEARCH = NO # The SEARCHENGINE_URL should point to a search engine hosted by a web server # which will return the search results when EXTERNAL_SEARCH is enabled. -# -# Doxygen ships with an example indexer ( doxyindexer) and search engine +# +# Doxygen ships with an example indexer (doxyindexer) and search engine # (doxysearch.cgi) which are based on the open source search engine library # Xapian (see: http://xapian.org/). See the section "External Indexing and # Searching" for details. # This tag requires that the tag SEARCHENGINE is set to YES. -SEARCHENGINE_URL = +SEARCHENGINE_URL = # When SERVER_BASED_SEARCH and EXTERNAL_SEARCH are both enabled the unindexed # search data is written to a file for indexing by an external tool. With the @@ -1536,7 +1572,7 @@ SEARCHDATA_FILE = searchdata.xml # projects and redirect the results back to the right project. # This tag requires that the tag SEARCHENGINE is set to YES. -EXTERNAL_SEARCH_ID = +EXTERNAL_SEARCH_ID = # The EXTRA_SEARCH_MAPPINGS tag can be used to enable searching through doxygen # projects other than the one defined by this configuration file, but that are @@ -1546,13 +1582,13 @@ EXTERNAL_SEARCH_ID = # EXTRA_SEARCH_MAPPINGS = tagname1=loc1 tagname2=loc2 ... # This tag requires that the tag SEARCHENGINE is set to YES. -EXTRA_SEARCH_MAPPINGS = +EXTRA_SEARCH_MAPPINGS = #--------------------------------------------------------------------------- # Configuration options related to the LaTeX output #--------------------------------------------------------------------------- -# If the GENERATE_LATEX tag is set to YES doxygen will generate LaTeX output. +# If the GENERATE_LATEX tag is set to YES, doxygen will generate LaTeX output. # The default value is: YES. GENERATE_LATEX = NO @@ -1567,7 +1603,7 @@ LATEX_OUTPUT = latex # The LATEX_CMD_NAME tag can be used to specify the LaTeX command name to be # invoked. -# +# # Note that when enabling USE_PDFLATEX this option is only used for generating # bitmaps for formulas in the HTML output, but not in the Makefile that is # written to the output directory. @@ -1583,7 +1619,7 @@ LATEX_CMD_NAME = latex MAKEINDEX_CMD_NAME = makeindex -# If the COMPACT_LATEX tag is set to YES doxygen generates more compact LaTeX +# If the COMPACT_LATEX tag is set to YES, doxygen generates more compact LaTeX # documents. This may be useful for small projects and may help to save some # trees in general. # The default value is: NO. @@ -1601,38 +1637,54 @@ COMPACT_LATEX = YES PAPER_TYPE = a4 # The EXTRA_PACKAGES tag can be used to specify one or more LaTeX package names -# that should be included in the LaTeX output. To get the times font for -# instance you can specify -# EXTRA_PACKAGES=times +# that should be included in the LaTeX output. The package can be specified just +# by its name or with the correct syntax as to be used with the LaTeX +# \usepackage command. To get the times font for instance you can specify : +# EXTRA_PACKAGES=times or EXTRA_PACKAGES={times} +# To use the option intlimits with the amsmath package you can specify: +# EXTRA_PACKAGES=[intlimits]{amsmath} # If left blank no extra packages will be included. # This tag requires that the tag GENERATE_LATEX is set to YES. -EXTRA_PACKAGES = +EXTRA_PACKAGES = # The LATEX_HEADER tag can be used to specify a personal LaTeX header for the # generated LaTeX document. The header should contain everything until the first # chapter. If it is left blank doxygen will generate a standard header. See # section "Doxygen usage" for information on how to let doxygen write the # default header to a separate file. -# +# # Note: Only use a user-defined header if you know what you are doing! The # following commands have a special meaning inside the header: $title, -# $datetime, $date, $doxygenversion, $projectname, $projectnumber. Doxygen will -# replace them by respectively the title of the page, the current date and time, -# only the current date, the version number of doxygen, the project name (see -# PROJECT_NAME), or the project number (see PROJECT_NUMBER). +# $datetime, $date, $doxygenversion, $projectname, $projectnumber, +# $projectbrief, $projectlogo. Doxygen will replace $title with the empty +# string, for the replacement values of the other commands the user is referred +# to HTML_HEADER. # This tag requires that the tag GENERATE_LATEX is set to YES. -LATEX_HEADER = +LATEX_HEADER = # The LATEX_FOOTER tag can be used to specify a personal LaTeX footer for the # generated LaTeX document. The footer should contain everything after the last -# chapter. If it is left blank doxygen will generate a standard footer. -# +# chapter. If it is left blank doxygen will generate a standard footer. See +# LATEX_HEADER for more information on how to generate a default footer and what +# special commands can be used inside the footer. +# # Note: Only use a user-defined footer if you know what you are doing! # This tag requires that the tag GENERATE_LATEX is set to YES. -LATEX_FOOTER = +LATEX_FOOTER = + +# The LATEX_EXTRA_STYLESHEET tag can be used to specify additional user-defined +# LaTeX style sheets that are included after the standard style sheets created +# by doxygen. Using this option one can overrule certain style aspects. Doxygen +# will copy the style sheet files to the output directory. +# Note: The order of the extra style sheet files is of importance (e.g. the last +# style sheet in the list overrules the setting of the previous ones in the +# list). +# This tag requires that the tag GENERATE_LATEX is set to YES. + +LATEX_EXTRA_STYLESHEET = # The LATEX_EXTRA_FILES tag can be used to specify one or more extra images or # other source files which should be copied to the LATEX_OUTPUT output @@ -1640,7 +1692,7 @@ LATEX_FOOTER = # markers available. # This tag requires that the tag GENERATE_LATEX is set to YES. -LATEX_EXTRA_FILES = +LATEX_EXTRA_FILES = # If the PDF_HYPERLINKS tag is set to YES, the LaTeX that is generated is # prepared for conversion to PDF (using ps2pdf or pdflatex). The PDF file will @@ -1651,8 +1703,8 @@ LATEX_EXTRA_FILES = PDF_HYPERLINKS = YES -# If the LATEX_PDFLATEX tag is set to YES, doxygen will use pdflatex to generate -# the PDF file directly from the LaTeX files. Set this option to YES to get a +# If the USE_PDFLATEX tag is set to YES, doxygen will use pdflatex to generate +# the PDF file directly from the LaTeX files. Set this option to YES, to get a # higher quality PDF documentation. # The default value is: YES. # This tag requires that the tag GENERATE_LATEX is set to YES. @@ -1677,7 +1729,7 @@ LATEX_HIDE_INDICES = NO # If the LATEX_SOURCE_CODE tag is set to YES then doxygen will include source # code with syntax highlighting in the LaTeX output. -# +# # Note that which sources are shown also depends on other settings such as # SOURCE_BROWSER. # The default value is: NO. @@ -1693,11 +1745,19 @@ LATEX_SOURCE_CODE = NO LATEX_BIB_STYLE = plain +# If the LATEX_TIMESTAMP tag is set to YES then the footer of each generated +# page will contain the date and time when the page was generated. Setting this +# to NO can help when comparing the output of multiple runs. +# The default value is: NO. +# This tag requires that the tag GENERATE_LATEX is set to YES. + +LATEX_TIMESTAMP = NO + #--------------------------------------------------------------------------- # Configuration options related to the RTF output #--------------------------------------------------------------------------- -# If the GENERATE_RTF tag is set to YES doxygen will generate RTF output. The +# If the GENERATE_RTF tag is set to YES, doxygen will generate RTF output. The # RTF output is optimized for Word 97 and may not look too pretty with other RTF # readers/editors. # The default value is: NO. @@ -1712,7 +1772,7 @@ GENERATE_RTF = NO RTF_OUTPUT = rtf -# If the COMPACT_RTF tag is set to YES doxygen generates more compact RTF +# If the COMPACT_RTF tag is set to YES, doxygen generates more compact RTF # documents. This may be useful for small projects and may help to save some # trees in general. # The default value is: NO. @@ -1725,7 +1785,7 @@ COMPACT_RTF = NO # output) instead of page references. This makes the output suitable for online # browsing using Word or some other Word compatible readers that support those # fields. -# +# # Note: WordPad (write) and others do not support links. # The default value is: NO. # This tag requires that the tag GENERATE_RTF is set to YES. @@ -1735,25 +1795,35 @@ RTF_HYPERLINKS = NO # Load stylesheet definitions from file. Syntax is similar to doxygen's config # file, i.e. a series of assignments. You only have to provide replacements, # missing definitions are set to their default value. -# +# # See also section "Doxygen usage" for information on how to generate the # default style sheet that doxygen normally uses. # This tag requires that the tag GENERATE_RTF is set to YES. -RTF_STYLESHEET_FILE = +RTF_STYLESHEET_FILE = # Set optional variables used in the generation of an RTF document. Syntax is # similar to doxygen's config file. A template extensions file can be generated # using doxygen -e rtf extensionFile. # This tag requires that the tag GENERATE_RTF is set to YES. -RTF_EXTENSIONS_FILE = +RTF_EXTENSIONS_FILE = + +# If the RTF_SOURCE_CODE tag is set to YES then doxygen will include source code +# with syntax highlighting in the RTF output. +# +# Note that which sources are shown also depends on other settings such as +# SOURCE_BROWSER. +# The default value is: NO. +# This tag requires that the tag GENERATE_RTF is set to YES. + +RTF_SOURCE_CODE = NO #--------------------------------------------------------------------------- # Configuration options related to the man page output #--------------------------------------------------------------------------- -# If the GENERATE_MAN tag is set to YES doxygen will generate man pages for +# If the GENERATE_MAN tag is set to YES, doxygen will generate man pages for # classes and files. # The default value is: NO. @@ -1777,6 +1847,13 @@ MAN_OUTPUT = man MAN_EXTENSION = .3 +# The MAN_SUBDIR tag determines the name of the directory created within +# MAN_OUTPUT in which the man pages are placed. If defaults to man followed by +# MAN_EXTENSION with the initial . removed. +# This tag requires that the tag GENERATE_MAN is set to YES. + +MAN_SUBDIR = + # If the MAN_LINKS tag is set to YES and doxygen generates man output, then it # will generate one additional man file for each entity documented in the real # man page(s). These additional files only source the real man page, but without @@ -1790,7 +1867,7 @@ MAN_LINKS = NO # Configuration options related to the XML output #--------------------------------------------------------------------------- -# If the GENERATE_XML tag is set to YES doxygen will generate an XML file that +# If the GENERATE_XML tag is set to YES, doxygen will generate an XML file that # captures the structure of the code including all documentation. # The default value is: NO. @@ -1804,19 +1881,7 @@ GENERATE_XML = NO XML_OUTPUT = xml -# The XML_SCHEMA tag can be used to specify a XML schema, which can be used by a -# validating XML parser to check the syntax of the XML files. -# This tag requires that the tag GENERATE_XML is set to YES. - -XML_SCHEMA = - -# The XML_DTD tag can be used to specify a XML DTD, which can be used by a -# validating XML parser to check the syntax of the XML files. -# This tag requires that the tag GENERATE_XML is set to YES. - -XML_DTD = - -# If the XML_PROGRAMLISTING tag is set to YES doxygen will dump the program +# If the XML_PROGRAMLISTING tag is set to YES, doxygen will dump the program # listings (including syntax highlighting and cross-referencing information) to # the XML output. Note that enabling this will significantly increase the size # of the XML output. @@ -1829,7 +1894,7 @@ XML_PROGRAMLISTING = YES # Configuration options related to the DOCBOOK output #--------------------------------------------------------------------------- -# If the GENERATE_DOCBOOK tag is set to YES doxygen will generate Docbook files +# If the GENERATE_DOCBOOK tag is set to YES, doxygen will generate Docbook files # that can be used to generate PDF. # The default value is: NO. @@ -1843,14 +1908,23 @@ GENERATE_DOCBOOK = NO DOCBOOK_OUTPUT = docbook +# If the DOCBOOK_PROGRAMLISTING tag is set to YES, doxygen will include the +# program listings (including syntax highlighting and cross-referencing +# information) to the DOCBOOK output. Note that enabling this will significantly +# increase the size of the DOCBOOK output. +# The default value is: NO. +# This tag requires that the tag GENERATE_DOCBOOK is set to YES. + +DOCBOOK_PROGRAMLISTING = NO + #--------------------------------------------------------------------------- # Configuration options for the AutoGen Definitions output #--------------------------------------------------------------------------- -# If the GENERATE_AUTOGEN_DEF tag is set to YES doxygen will generate an AutoGen -# Definitions (see http://autogen.sf.net) file that captures the structure of -# the code including all documentation. Note that this feature is still -# experimental and incomplete at the moment. +# If the GENERATE_AUTOGEN_DEF tag is set to YES, doxygen will generate an +# AutoGen Definitions (see http://autogen.sf.net) file that captures the +# structure of the code including all documentation. Note that this feature is +# still experimental and incomplete at the moment. # The default value is: NO. GENERATE_AUTOGEN_DEF = NO @@ -1859,15 +1933,15 @@ GENERATE_AUTOGEN_DEF = NO # Configuration options related to the Perl module output #--------------------------------------------------------------------------- -# If the GENERATE_PERLMOD tag is set to YES doxygen will generate a Perl module +# If the GENERATE_PERLMOD tag is set to YES, doxygen will generate a Perl module # file that captures the structure of the code including all documentation. -# +# # Note that this feature is still experimental and incomplete at the moment. # The default value is: NO. GENERATE_PERLMOD = NO -# If the PERLMOD_LATEX tag is set to YES doxygen will generate the necessary +# If the PERLMOD_LATEX tag is set to YES, doxygen will generate the necessary # Makefile rules, Perl scripts and LaTeX code to be able to generate PDF and DVI # output from the Perl module output. # The default value is: NO. @@ -1875,9 +1949,9 @@ GENERATE_PERLMOD = NO PERLMOD_LATEX = NO -# If the PERLMOD_PRETTY tag is set to YES the Perl module output will be nicely +# If the PERLMOD_PRETTY tag is set to YES, the Perl module output will be nicely # formatted so it can be parsed by a human reader. This is useful if you want to -# understand what is going on. On the other hand, if this tag is set to NO the +# understand what is going on. On the other hand, if this tag is set to NO, the # size of the Perl module output will be much smaller and Perl will parse it # just the same. # The default value is: YES. @@ -1891,20 +1965,20 @@ PERLMOD_PRETTY = YES # overwrite each other's variables. # This tag requires that the tag GENERATE_PERLMOD is set to YES. -PERLMOD_MAKEVAR_PREFIX = +PERLMOD_MAKEVAR_PREFIX = #--------------------------------------------------------------------------- # Configuration options related to the preprocessor #--------------------------------------------------------------------------- -# If the ENABLE_PREPROCESSING tag is set to YES doxygen will evaluate all +# If the ENABLE_PREPROCESSING tag is set to YES, doxygen will evaluate all # C-preprocessor directives found in the sources and include files. # The default value is: YES. ENABLE_PREPROCESSING = YES -# If the MACRO_EXPANSION tag is set to YES doxygen will expand all macro names -# in the source code. If set to NO only conditional compilation will be +# If the MACRO_EXPANSION tag is set to YES, doxygen will expand all macro names +# in the source code. If set to NO, only conditional compilation will be # performed. Macro expansion can be done in a controlled way by setting # EXPAND_ONLY_PREDEF to YES. # The default value is: NO. @@ -1920,7 +1994,7 @@ MACRO_EXPANSION = NO EXPAND_ONLY_PREDEF = NO -# If the SEARCH_INCLUDES tag is set to YES the includes files in the +# If the SEARCH_INCLUDES tag is set to YES, the include files in the # INCLUDE_PATH will be searched if a #include is found. # The default value is: YES. # This tag requires that the tag ENABLE_PREPROCESSING is set to YES. @@ -1932,7 +2006,7 @@ SEARCH_INCLUDES = YES # preprocessor. # This tag requires that the tag SEARCH_INCLUDES is set to YES. -INCLUDE_PATH = +INCLUDE_PATH = # You can use the INCLUDE_FILE_PATTERNS tag to specify one or more wildcard # patterns (like *.h and *.hpp) to filter out the header-files in the @@ -1940,7 +2014,7 @@ INCLUDE_PATH = # used. # This tag requires that the tag ENABLE_PREPROCESSING is set to YES. -INCLUDE_FILE_PATTERNS = +INCLUDE_FILE_PATTERNS = # The PREDEFINED tag can be used to specify one or more macro names that are # defined before the preprocessor is started (similar to the -D option of e.g. @@ -1950,7 +2024,7 @@ INCLUDE_FILE_PATTERNS = # recursively expanded use the := operator instead of the = operator. # This tag requires that the tag ENABLE_PREPROCESSING is set to YES. -PREDEFINED = +PREDEFINED = # If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then this # tag can be used to specify a list of macro names that should be expanded. The @@ -1959,12 +2033,12 @@ PREDEFINED = # definition found in the source code. # This tag requires that the tag ENABLE_PREPROCESSING is set to YES. -EXPAND_AS_DEFINED = +EXPAND_AS_DEFINED = # If the SKIP_FUNCTION_MACROS tag is set to YES then doxygen's preprocessor will -# remove all refrences to function-like macros that are alone on a line, have an -# all uppercase name, and do not end with a semicolon. Such function macros are -# typically used for boiler-plate code, and will confuse the parser if not +# remove all references to function-like macros that are alone on a line, have +# an all uppercase name, and do not end with a semicolon. Such function macros +# are typically used for boiler-plate code, and will confuse the parser if not # removed. # The default value is: YES. # This tag requires that the tag ENABLE_PREPROCESSING is set to YES. @@ -1984,32 +2058,33 @@ SKIP_FUNCTION_MACROS = YES # where loc1 and loc2 can be relative or absolute paths or URLs. See the # section "Linking to external documentation" for more information about the use # of tag files. -# Note: Each tag file must have an unique name (where the name does NOT include +# Note: Each tag file must have a unique name (where the name does NOT include # the path). If a tag file is not located in the directory in which doxygen is # run, you must also specify the path to the tagfile here. -TAGFILES = +TAGFILES = # When a file name is specified after GENERATE_TAGFILE, doxygen will create a # tag file that is based on the input files it reads. See section "Linking to # external documentation" for more information about the usage of tag files. -GENERATE_TAGFILE = +GENERATE_TAGFILE = -# If the ALLEXTERNALS tag is set to YES all external class will be listed in the -# class index. If set to NO only the inherited external classes will be listed. +# If the ALLEXTERNALS tag is set to YES, all external class will be listed in +# the class index. If set to NO, only the inherited external classes will be +# listed. # The default value is: NO. ALLEXTERNALS = NO -# If the EXTERNAL_GROUPS tag is set to YES all external groups will be listed in -# the modules index. If set to NO, only the current project's groups will be +# If the EXTERNAL_GROUPS tag is set to YES, all external groups will be listed +# in the modules index. If set to NO, only the current project's groups will be # listed. # The default value is: YES. EXTERNAL_GROUPS = YES -# If the EXTERNAL_PAGES tag is set to YES all external pages will be listed in +# If the EXTERNAL_PAGES tag is set to YES, all external pages will be listed in # the related pages index. If set to NO, only the current project's pages will # be listed. # The default value is: YES. @@ -2026,7 +2101,7 @@ PERL_PATH = /usr/bin/perl # Configuration options related to the dot tool #--------------------------------------------------------------------------- -# If the CLASS_DIAGRAMS tag is set to YES doxygen will generate a class diagram +# If the CLASS_DIAGRAMS tag is set to YES, doxygen will generate a class diagram # (in HTML and LaTeX) for classes with base or super classes. Setting the tag to # NO turns the diagrams off. Note that this option also works with HAVE_DOT # disabled, but it is recommended to install and use dot, since it yields more @@ -2042,9 +2117,16 @@ CLASS_DIAGRAMS = YES # the mscgen tool resides. If left empty the tool is assumed to be found in the # default search path. -MSCGEN_PATH = +MSCGEN_PATH = + +# You can include diagrams made with dia in doxygen documentation. Doxygen will +# then run dia to produce the diagram and insert it in the documentation. The +# DIA_PATH tag allows you to specify the directory where the dia binary resides. +# If left empty dia is assumed to be found in the default search path. + +DIA_PATH = -# If set to YES, the inheritance and collaboration graphs will hide inheritance +# If set to YES the inheritance and collaboration graphs will hide inheritance # and usage relations if the target is undocumented or is not a class. # The default value is: YES. @@ -2069,7 +2151,7 @@ HAVE_DOT = YES DOT_NUM_THREADS = 0 -# When you want a differently looking font n the dot files that doxygen +# When you want a differently looking font in the dot files that doxygen # generates you can specify the font name using DOT_FONTNAME. You need to make # sure dot is able to find the font, which can be done by putting it in a # standard location or by setting the DOTFONTPATH environment variable or by @@ -2091,7 +2173,7 @@ DOT_FONTSIZE = 10 # the path where dot can find it using this tag. # This tag requires that the tag HAVE_DOT is set to YES. -DOT_FONTPATH = +DOT_FONTPATH = # If the CLASS_GRAPH tag is set to YES then doxygen will generate a graph for # each documented class showing the direct and indirect inheritance relations. @@ -2117,7 +2199,7 @@ COLLABORATION_GRAPH = YES GROUP_GRAPHS = YES -# If the UML_LOOK tag is set to YES doxygen will generate inheritance and +# If the UML_LOOK tag is set to YES, doxygen will generate inheritance and # collaboration diagrams in a style similar to the OMG's Unified Modeling # Language. # The default value is: NO. @@ -2166,10 +2248,11 @@ INCLUDED_BY_GRAPH = YES # If the CALL_GRAPH tag is set to YES then doxygen will generate a call # dependency graph for every global function or class method. -# +# # Note that enabling this option will significantly increase the time of a run. # So in most cases it will be better to enable call graphs for selected -# functions only using the \callgraph command. +# functions only using the \callgraph command. Disabling a call graph can be +# accomplished by means of the command \hidecallgraph. # The default value is: NO. # This tag requires that the tag HAVE_DOT is set to YES. @@ -2177,10 +2260,11 @@ CALL_GRAPH = YES # If the CALLER_GRAPH tag is set to YES then doxygen will generate a caller # dependency graph for every global function or class method. -# +# # Note that enabling this option will significantly increase the time of a run. # So in most cases it will be better to enable caller graphs for selected -# functions only using the \callergraph command. +# functions only using the \callergraph command. Disabling a caller graph can be +# accomplished by means of the command \hidecallergraph. # The default value is: NO. # This tag requires that the tag HAVE_DOT is set to YES. @@ -2203,11 +2287,15 @@ GRAPHICAL_HIERARCHY = YES DIRECTORY_GRAPH = YES # The DOT_IMAGE_FORMAT tag can be used to set the image format of the images -# generated by dot. +# generated by dot. For an explanation of the image formats see the section +# output formats in the documentation of the dot tool (Graphviz (see: +# http://www.graphviz.org/)). # Note: If you choose svg you need to set HTML_FILE_EXTENSION to xhtml in order # to make the SVG files visible in IE 9+ (other browsers do not have this # requirement). -# Possible values are: png, jpg, gif and svg. +# Possible values are: png, jpg, gif, svg, png:gd, png:gd:gd, png:cairo, +# png:cairo:gd, png:cairo:cairo, png:cairo:gdiplus, png:gdiplus and +# png:gdiplus:gdiplus. # The default value is: png. # This tag requires that the tag HAVE_DOT is set to YES. @@ -2215,7 +2303,7 @@ DOT_IMAGE_FORMAT = png # If DOT_IMAGE_FORMAT is set to svg, then this option can be set to YES to # enable generation of interactive SVG images that allow zooming and panning. -# +# # Note that this requires a modern browser other than Internet Explorer. Tested # and working are Firefox, Chrome, Safari, and Opera. # Note: For IE 9+ you need to set HTML_FILE_EXTENSION to xhtml in order to make @@ -2236,13 +2324,32 @@ DOT_PATH = @DOXYGEN_DOT_PATH@ # command). # This tag requires that the tag HAVE_DOT is set to YES. -DOTFILE_DIRS = +DOTFILE_DIRS = # The MSCFILE_DIRS tag can be used to specify one or more directories that # contain msc files that are included in the documentation (see the \mscfile # command). -MSCFILE_DIRS = +MSCFILE_DIRS = + +# The DIAFILE_DIRS tag can be used to specify one or more directories that +# contain dia files that are included in the documentation (see the \diafile +# command). + +DIAFILE_DIRS = + +# When using plantuml, the PLANTUML_JAR_PATH tag should be used to specify the +# path where java can find the plantuml.jar file. If left blank, it is assumed +# PlantUML is not used or called during a preprocessing step. Doxygen will +# generate a warning when it encounters a \startuml command in this case and +# will not generate output for the diagram. + +PLANTUML_JAR_PATH = + +# When using plantuml, the specified paths are searched for files specified by +# the !include statement in a plantuml block. + +PLANTUML_INCLUDE_PATH = # The DOT_GRAPH_MAX_NODES tag can be used to set the maximum number of nodes # that will be shown in the graph. If the number of nodes in a graph becomes @@ -2271,7 +2378,7 @@ MAX_DOT_GRAPH_DEPTH = 0 # Set the DOT_TRANSPARENT tag to YES to generate images with a transparent # background. This is disabled by default, because dot on Windows does not seem # to support this out of the box. -# +# # Warning: Depending on the platform used, enabling this option may lead to # badly anti-aliased labels on the edges of a graph (i.e. they become hard to # read). @@ -2280,7 +2387,7 @@ MAX_DOT_GRAPH_DEPTH = 0 DOT_TRANSPARENT = NO -# Set the DOT_MULTI_TARGETS tag to YES allow dot to generate multiple output +# Set the DOT_MULTI_TARGETS tag to YES to allow dot to generate multiple output # files in one run (i.e. multiple -o and -T options on the command line). This # makes dot run faster, but since only newer versions of dot (>1.8.10) support # this, this feature is disabled by default. @@ -2297,7 +2404,7 @@ DOT_MULTI_TARGETS = NO GENERATE_LEGEND = YES -# If the DOT_CLEANUP tag is set to YES doxygen will remove the intermediate dot +# If the DOT_CLEANUP tag is set to YES, doxygen will remove the intermediate dot # files that are used to generate the various graphs. # The default value is: YES. # This tag requires that the tag HAVE_DOT is set to YES. diff --git a/docs/sphinx/conf.py.in b/docs/sphinx/conf.py.in index c6d9070893cec0f914d4809afb3758b7325d9e76..a1622861b790d827a166814924dc5ad2e250f6ac 100644 --- a/docs/sphinx/conf.py.in +++ b/docs/sphinx/conf.py.in @@ -26,6 +26,8 @@ # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. import sys +import sphinx_bootstrap_theme + extensions = [ 'sphinx.ext.autodoc', 'sphinx.ext.autosummary', @@ -138,29 +140,12 @@ html_theme = 'bootstrap' # Add any paths that contain custom themes here, relative to this directory. html_theme_path = ['@CMAKE_CURRENT_SOURCE_DIR@/sphinx/themes'] -import sphinx_bootstrap_theme html_theme_path += sphinx_bootstrap_theme.get_html_theme_path() # Theme options are theme-specific and customize the look and feel of a theme # further. For a list of options available for each theme, see the # documentation. html_theme_options = {'bootswatch_theme': "sandstone", } -# Activate the bootstrap theme if it's available -#try: -# import sphinx_bootstrap_theme -# html_theme = 'bootstrap' -# html_theme_path += sphinx_bootstrap_theme.get_html_theme_path() -# print "ooooooooooooooooooooo ", html_theme_path -# html_theme_options = { -# 'bootswatch_theme': "superhero", -# } -#except: -# print("""Warning: Module sphinx_bootstrap_theme not found, continuing -#anyway using old theme. To get the new theme, use,# -# -# pip install sphinx_bootstrap_theme""") - - # The name for this set of Sphinx documents. If None, it defaults to # "<project> v<release> documentation". # html_title = None diff --git a/docs/sphinx/index.rst.in b/docs/sphinx/index.rst.in index 1f86f101d59afb2a504e5898a9270a0cdbfc4127..a70e47b6ef165447fb4be63b39982cb68c75616d 100644 --- a/docs/sphinx/index.rst.in +++ b/docs/sphinx/index.rst.in @@ -13,7 +13,7 @@ HySoP software: hybrid simulation with particles getting_started/index users_guide/index examples/index - + reference/index license contacts diff --git a/docs/sphinx/users_guide/finite_differences.rst b/docs/sphinx/users_guide/finite_differences.rst index 81d277931d12716dff1f41adee43ac36cf2b6e3f..24671a8cc6e2c3a7e89594fcb795aaf8e8dc0631 100644 --- a/docs/sphinx/users_guide/finite_differences.rst +++ b/docs/sphinx/users_guide/finite_differences.rst @@ -5,17 +5,17 @@ Finite differences schemes -------------------------- -Differentiate some fields in one direction using finite differences. +To differentiate some field(s) in one direction using finite differences. -So, to compute +If tab is a numpy array representing the discrete values of a scalar field :math:`\rho` +on a grid, then to compute .. math:: \begin{eqnarray} result &=& \frac{\partial \rho}{\partial y} \end{eqnarray} -if tab is a numpy array representing the discrete values of the scalar field :math:`\rho` -on a grid, then the basic usage of such schemes is : +the basic usage of such schemes is : .. code:: @@ -25,10 +25,10 @@ on a grid, then the basic usage of such schemes is : dir = 1 # y direction result = scheme(tab, dir, result) -This will compute : +with : -result = diff(tab[indices], dir), diff depending on -the chosen scheme. +result[iout] = diff(tab[indices], dir), diff depending on +the chosen scheme. See examples below to understand properly indices and iout. Available schemes @@ -42,75 +42,96 @@ Available schemes A few important remarks ^^^^^^^^^^^^^^^^^^^^^^^ +* Mind that 'tab' arg is a numpy array not a field. Usually, + you will have to use some_field.data[i] as tab, some_field being a :class:`~hysop.fields.discrete.DiscreteField` and i the component number. * step is the space discretization step size in each direction, i.e. a list or numpy array - with d values, d being the dimension of the domain. A common case is :code:: + with d values, d being the dimension of the domain. A common case is code:: step = topo.mesh.space_step * indices represent the set of points on which the scheme must be applied. This is usually - a list of slices, for example, :code:: + a list of slices, for example, code:: - indices = topo.mesh.iCompute + indices = topo.mesh.compute_index * result must be a predefined numpy array of the 'right' size, here the same size/shape as tab. * To optimize memory usage and computational time, it's possible to reduce the size of - the output and/or to apply the scheme on a subset of the domain. All available possibilities - are summarized through the examples below. - -* the size of the ghost layer depends on the scheme but is not checked! You must ensure - that topo.ghosts() >= scheme.minimal_ghost_layer. + the output and/or to apply the scheme on a subset of the domain. + scheme.indices represents the indices (in topo.mesh) of points for which the derivative will be computed. + scheme.output_indices are positions in result where the resulting derivative will be saved. + All available possibilities are summarized through the examples below. + +* The size of the ghost layer depends on the scheme but is not checked! You must ensure + that topo.ghosts() >= scheme.ghosts_layer_size. -* for some schemes a work array may be provided during call. It must be a numpy array of - the same size of the result. It's shape is not really important, since during call +* For some schemes a work array may be provided during call. If so, it must be a numpy array of + the same size as result. It's shape is not really important, since during call work will be reshaped to be of the same shape as result. This allows us to provide 1D array that may be shared between different operators/methods whatever the internal required shapes are. * Notice that most of the time, finite-differences are defined as internal methods of operators and work arrays management, indices list or ghosts layers are set/checked internally. + Default case : apply scheme on all the domain with a full-size result ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +This is the standard case where result and input tab have the same shape. +Computation is performed on all grid points, excluding ghosts. + .. code:: from hysop.numerics.finite_differences import FDC2 + # --- Domain, discretization and topology --- box = Box(length=[1., 1.]) d2d = Discretization([33, 33], [1, 1]) topo = box.create_topology(d2d) + # --- Scalar field on the domain --- rho = Field(box, name='rho') # discretize and initialize rho rd = rho.randomize(topo) + # --- Grid parameters required for FD scheme --- # Get 'computational' points, i.e. grid points excluding ghosts - ic = topo.mesh.iCompute + ic = topo.mesh.compute_index # space step step = topo.mesh.space_step # field resolution on the grid defined by topo shape = topo.mesh.resolution + # create a new array to save results result = npw.zeros(shape) + # --- Create FD scheme --- scheme = FDC2(step, ic) - assert (topo.ghosts() >= scheme.minimal_ghost_layer).all() + assert (topo.ghosts() >= scheme.ghosts_layer_size).all() + # --- Apply FD scheme --- + # compute first derivative of scal on all points (except ghosts) + # of the 'topo' discretization of the domain result = scheme(rd.data[0], 1, result) In that case: * result.shape = (34, 34) -* scheme.indices = [slice(1,33), slice(1,33)] -* scheme.output_indices = [slice(0,32), slice(0,32)] +* scheme.indices = [slice(1,33), slice(1,33)] = topo.mesh.compute_index +* scheme.output_indices = scheme.indices +Notice that ghost points of result are not updated! Apply scheme on all the domain with a reduced result ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -If you do not want to allocate ghosts points for the result, then use: +Suppose that you want to compute a derivative on all the domain but +use non-standard shape array to save the result. This may be useful +if you do not want to allocate memory for ghost points in result. .. code:: - + + # Create a 'reduced-shape' result + # result shape == shape of topo.mesh without ghosts shape = np.asarray(topo.mesh.resolution).copy() shape -= 2 * topo.ghosts() shape = tuple(shape) result = npw.zeros(shape) - scheme = FDC2(step, ic, indices_out=True) - assert (topo.ghosts() >= scheme.minimal_ghost_layer).all() + # build and apply the scheme + scheme = FDC2(step, ic, output_indices=True) result = scheme(rd.data[0], 1, result) In that case: @@ -120,8 +141,8 @@ In that case: * scheme.output_indices = [slice(0,32), slice(0,32)] -Apply scheme on a subset of the domain with a full-size result -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Apply scheme on a subset of the domain +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. code:: @@ -131,14 +152,18 @@ Apply scheme on a subset of the domain with a full-size result sl = topo.domain.length * 0.5 orig = topo.domain.origin + 0.1 * topo.domain.length subbox = SubBox(parent=topo.domain, origin=orig, length=sl) + # then we get indices of points inside this subdomain, for + # the discretization 'topo' indices = subbox.discretize(topo)[0] + # and build a FD scheme for this set of indices scheme = FDC2(step, indices) result = npw.zeros_like(rd.data[0]) result = scheme(rd.data[0], 1, result) -In that case: +In that case derivative is computed only inside subdomain and +only the corresponding subpart of result is updated: -* result.shape = (34,34) +* result.shape = (34,34), i.e. same shape as the complete mesh of topo. * scheme.indices = [slice(4,21), slice(4,21)] * scheme.output_indices = scheme.indices @@ -146,6 +171,9 @@ In that case: Apply scheme on a subset of the domain with a reduced-size result ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Same as above but to limit memory print, a reduced-shape result +is used + .. code:: # We define a subset of the domain, @@ -155,7 +183,9 @@ Apply scheme on a subset of the domain with a reduced-size result orig = topo.domain.origin + 0.1 * topo.domain.length subbox = SubBox(parent=topo.domain, origin=orig, length=sl) indices = subbox.discretize(topo)[0] - scheme = FDC2(step, indices, indices_out=True) + # build and apply a FD scheme, notice the output_indices arg! + scheme = FDC2(step, indices, output_indices=True) + # result shape is computed using discretized subdomain shape shape = subbox.mesh[topo].resolution result = npw.zeros(shape) result = scheme(rd.data[0], 1, result) @@ -164,7 +194,7 @@ In that case: * result.shape = (17,17) * scheme.indices = [slice(4,21), slice(4,21)] -* scheme.output_indices = [slice(0,17), slice(0,17)]scheme.indices +* scheme.output_indices = [slice(0,17), slice(0,17)] Apply scheme on a subset of the domain with a predifined-size result @@ -183,17 +213,18 @@ the values of the computed derivative. subbox = SubBox(parent=topo.domain, origin=orig, length=sl) indices = subbox.discretize(topo)[0] shape = subbox.mesh[topo].resolution - result = npw.zeros_like(rd.data[0) + result = npw.zeros_like(rd.data[0]) + # output_indices are explicitely given iout = [slice(2, 19), slice(3, 20)] - scheme = FDC2(step, indices, indices_out=iout) + scheme = FDC2(step, indices, output_indices=iout) result = scheme(rd.data[0], 1, result) In that case: * result.shape = (17,17) * scheme.indices = [slice(4,21), slice(4,21)] -* scheme.output_indices = iout +* scheme.output_indices = iout = [slice(2, 19), slice(3, 20)] -This option is usefull when input tab is a work array, with a size different +This option is useful when input tab is a work array, with a size different from the topology resolution, and the result a field defined on the whole grid. diff --git a/docs/sphinx/users_guide/odesolvers.rst b/docs/sphinx/users_guide/odesolvers.rst index 5408afcb89936ef81d24eb3e1af033ba7e780af9..13909c7b213cb6817462bde64ef9f5112095660d 100644 --- a/docs/sphinx/users_guide/odesolvers.rst +++ b/docs/sphinx/users_guide/odesolvers.rst @@ -13,7 +13,7 @@ Consider the following problem: y(x, t_0) &=& y_0 -with y and rhs some vector fields. +with y and rhs some vector fields defined on a 1,2 or 3D domain. The following time integrators are implemented in hysop: @@ -25,10 +25,11 @@ Usage :: - # Build - ti = Euler(nb_components=2, topo=topo, f=func) + # Build the system with 'some_func' as right-hand side + ti = Euler(nb_components=2, topo=topo, f=some_func) - # Solve + # integrate between t and t+dt, with y representing y(x,t) on topo, + # save y(x,t+dt) in result: result = ti(t, y, dt, result) @@ -36,16 +37,15 @@ Usage To build a time integrator, it is necessary to define a python function to compute the right-hand side. This function must be of the following form:: - def some_name(t, y, work): + def some_func(t, y, work): work[0][...] = ... work[1][...] = ... ... return work -with y and work some lists of numpy arrays of length nb_components, nb_components being the number of -components of the right-hand side. y must not be modified by the function and work must contain the result of the computation. t is the time. - - - -to set the number of components of the rhs +* nb_components being the number of components of the right-hand side. +* y and work some lists of numpy arrays of length nb_components +* t is the time +* y must not be modified by the function +* work must contain the result of the computation diff --git a/hysop/__init__.py.in b/hysop/__init__.py.in index 4e2d010be9cb72165c2d4eb6350878e4af3bd4f3..664c72b139d5f14194c53ee683d06d9852307ddc 100755 --- a/hysop/__init__.py.in +++ b/hysop/__init__.py.in @@ -69,4 +69,3 @@ __DEFAULT_PLATFORM_ID__ = @OPENCL_DEFAULT_OPENCL_PLATFORM_ID@ __DEFAULT_DEVICE_ID__ = @OPENCL_DEFAULT_OPENCL_DEVICE_ID@ version = "1.0.0" - diff --git a/hysop/default_methods.py b/hysop/default_methods.py index 012ee19c356b4fafd4ec93b7479837ead049a21f..df41a328e367e5f79221d04f65ae1954c4646eb2 100644 --- a/hysop/default_methods.py +++ b/hysop/default_methods.py @@ -8,11 +8,13 @@ from hysop.numerics.odesolvers import RK2, RK3 from hysop.numerics.interpolation import Linear from hysop.numerics.remeshing import L2_1, Linear as Rmsh_Linear from hysop.numerics.finite_differences import FDC4, FDC2 +from hysop.problem.simulation import O2 + #from hysop.operator.discrete.stretching import Conservative ADVECTION = {TimeIntegrator: RK2, Interpolation: Linear, - Remesh: L2_1, Support: '', Splitting: 'o2', MultiScale: L2_1, + Remesh: L2_1, Support: '', Splitting: O2, MultiScale: L2_1, Precision: HYSOP_REAL} """Advection operators default setup """ diff --git a/hysop/fields/discrete.py b/hysop/fields/discrete.py index 264e1f5e41f0a79c4b2fba8fdf5f4d2275abe1ab..9a09128e48dfe617a60ef1e7cc344fd8ddc8184a 100755 --- a/hysop/fields/discrete.py +++ b/hysop/fields/discrete.py @@ -28,7 +28,7 @@ class DiscreteField(object): return object.__new__(cls, *args, **kw) @debug - def __init__(self, topology, is_vector=False, name=None, + def __init__(self, topology, name, is_vector=False, nb_components=None): """Creates a discrete field for a given topology. @@ -65,11 +65,8 @@ class DiscreteField(object): # Id (unique) for the field self.__id = self.__field_counter.next() # Field name. - if name is not None: - self.name = name - else: - self.name = 'unamed' - #: Field dimension. + self.name = name + # Field dimension. self.dimension = topology.domain.dimension # Field resolution. self.resolution = self.topology.mesh.resolution diff --git a/hysop/gpu/__init__.py b/hysop/gpu/__init__.py index a7389ea91dbd9d32e02696d5006599c686bc448d..197c5adb8decd131fff87771bb14404b9054d75a 100644 --- a/hysop/gpu/__init__.py +++ b/hysop/gpu/__init__.py @@ -28,16 +28,19 @@ see hysop.gpu.tools.parse_file import pyopencl import pyopencl.tools import pyopencl.array -## open cl underlying implementation +import os + cl = pyopencl -## PyOpenCL tools +"""open cl underlying implementation""" + clTools = pyopencl.tools -## PyOpenCL arrays +"""PyOpencl tools""" + clArray = pyopencl.array +"""PyOpenCL arrays""" -import os -## GPU deflault sources GPU_SRC = os.path.join(__path__[0], "cl_src", '') +"""Default path to cl kernels source files""" -## If use OpenCL profiling events to time computations CL_PROFILE = False +"""Boolean, true to enable OpenCL profiling events to time computations""" diff --git a/hysop/gpu/gpu_discrete.py b/hysop/gpu/gpu_discrete.py index d51a5fa23e4accd2a35507c6c9c0a3ebba80c93b..bc418c48a18470e4a1d697b140bf016116f7534f 100644 --- a/hysop/gpu/gpu_discrete.py +++ b/hysop/gpu/gpu_discrete.py @@ -1,7 +1,4 @@ -""" -@file gpu_discrete.py - -Contains class for discrete fields on GPU. +"""Discrete field defined on device (GPU) """ from hysop import __VERBOSE__ from hysop.constants import ORDER, np,\ @@ -57,57 +54,65 @@ toLayoutMgrFunc_2D = [ class GPUDiscreteField(DiscreteField): - """ - GPU Discrete vector field implementation. + """GPU Discrete vector field implementation. Allocates OpenCL device memory for the field. """ def __init__(self, cl_env, topology=None, is_vector=False, name="?", precision=HYSOP_REAL, layout=True, simple_layout=False): + """GPU Discrete vector field implementation. + Allocates OpenCL device memory for the field. + + Parameters + ---------- + + queue : OpenCL queue + topology : :class:`~hysop.mpi.topology.Cartesian`, optional + mpi topology and local meshes info + precision : np type, optional + Floating point precision, + default=:data:`~hysop.constants.HYSOP_REAL` + is_vector: boolean, optional + true if parent field is a vector field, default=False + name : string, optional + Field name + layour : boolean, optional + indicates if components are arranged in memory, default=True + i.e. all components are considered in the same way. + simple_layout : boolean, optional + indicates if in the Z direction, layout is ZYX (simple) or ZXY. """ - Constructor. - @param queue : OpenCL queue - @param precision : Floating point precision - @param parent : Continuous field. - @param topology : Topology informations - @param name : Field name - @param idFromParent : Index in the parent's discrete fields - @param layout : Boolean indicating if components are arranged in memory - Defaut : all components are considered in the same way. - @param simple_layout : Boolean indicating if in the Z direction, - layout is ZYX (simple) or ZXY. - @see hysop.fields.vector.VectorField.__init__ - """ + # init base class super(GPUDiscreteField, self).__init__(topology, is_vector, name) - ## OpenCL environment + # OpenCL environment self.cl_env = cl_env - ## Precision for the field + # Precision for the field self.precision = precision - ## Memory used + # Memory used self.mem_size = 0 ## Initialization OpenCL kernel as KernelLauncher self.init_kernel = None self._isReleased = False ## OpenCL Buffer pointer self.gpu_data = [None] * self.nb_components - ## Is the device allocations are performed + # True if device allocations have been done, + # (self.allocate call) self.gpu_allocated = False ## OpenCL Events list modifying this field self.events = [] - # Get the process number involved in this field discretisation - # By default, all mpi process are take, otherwise, user create and - # gives his own topologies. + # Get the ids of processes involved in the field discretisation. + # Default = all, otherwise, get info from input topology if given. if topology is None: from hysop.mpi import main_rank self._rank = main_rank else: self._rank = topology.rank - ## Data layout is direction dependant + # Data layout is direction dependant self.layout = layout - ## Layout for the Z direction + # Layout for the Z direction self.simple_layout = simple_layout - ## Layout and shape managers + # Layout and shape managers if self.domain.dimension == 3: if self.simple_layout: self._shapeFunc = shapeFunc_3D @@ -124,12 +129,12 @@ class GPUDiscreteField(DiscreteField): self.profiler += FProfiler("Transfer_toHost") self.profiler += FProfiler("Transfer_toDevice") - ## Transfer size counter (to device) + # Transfer size counter (to device) self.to_dev_size = 0. - ## Transfer size counter (to host) + # Transfer size counter (to host) self.to_host_size = 0. - ## Temporary cpu buffer to change data layout between cpu ang gpu + # Temporary cpu buffer to change data layout between cpu ang gpu self.host_data_pinned = [None, ] * self.nb_components def allocate(self): @@ -137,9 +142,12 @@ class GPUDiscreteField(DiscreteField): if not self.gpu_allocated: evt = [None, ] * self.nb_components for d in xrange(self.nb_components): + # convert data to required precision self.data[d] = np.asarray(self.data[d], dtype=self.precision, order=ORDER) + # create on-device buffer self.gpu_data[d] = self.cl_env.global_allocation(self.data[d]) + # update memory counter self.mem_size += self.gpu_data[d].size self.host_data_pinned[d], evt[d] = cl.enqueue_map_buffer( self.cl_env.queue, diff --git a/hysop/gpu/gpu_kernel.py b/hysop/gpu/gpu_kernel.py index 3508f316f69b905aaf2d7f918544b9c6b10deee8..fa968cc486f2ed9936895be300a2ba3454dd82f8 100644 --- a/hysop/gpu/gpu_kernel.py +++ b/hysop/gpu/gpu_kernel.py @@ -1,25 +1,29 @@ """ -@file gpu_kernel.py """ from hysop.constants import debug, S_DIR from hysop import __VERBOSE__ from hysop.gpu import cl, CL_PROFILE from hysop.tools.profiler import FProfiler + class KernelListLauncher(object): - """ - OpenCL kernel list launcher. + """OpenCL kernel list launcher. Manage launching of OpenCL kernels as a list. """ @debug def __init__(self, kernel, queue, gsize, lsize=None): - """ - Create a kernel list launcher. - @param kernel : kernel list. - @param queue : OpenCL command queue. - @param gsize : OpenCL global size index. - @param lsize : OpenCL local size index. + """Create a kernel list launcher. + + Parameters + ---------- + + kernel : list + queue : OpenCL command queue + gsize : int + OpenCL global size index. + lsize : int, optional + OpenCL local size index. """ ## OpenCL Kernel list self.kernel = kernel diff --git a/hysop/gpu/gpu_operator.py b/hysop/gpu/gpu_operator.py index 523bb6c99ff34acce24c40e7bf12656d4a5c02b1..d0e4f001134f24396aa1e4eede2c070a1dc7dd76 100644 --- a/hysop/gpu/gpu_operator.py +++ b/hysop/gpu/gpu_operator.py @@ -1,75 +1,82 @@ -""" -@file gpu_operator.py +"""Abstract class providing a common interface to all +discrete operators working on GPU. + +* :class:`~hysop.gpu_operator.GPUOperator` is an abstract class + used to provide a common interface to all discrete operators + working on GPU. + See for example :class:`~hysop.gpu.gpu_diffusion.GPUDiffusion` or + :class:`~hysop.gpu.gpu_particle_advection.GPUParticleAdvection`. -Discrete operator for GPU architecture. """ from abc import ABCMeta from hysop.constants import HYSOP_REAL, S_DIR -from hysop.methods_keys import Precision from hysop.gpu.tools import get_opencl_environment class GPUOperator(object): - """ - Abstract class for GPU operators. - - In practice, discrete operators must inherit from the classic - version and this abstract layer. + """Abstract class for discrete operators working on GPU. """ __metaclass__ = ABCMeta - def __init__(self, platform_id=None, device_id=None, device_type=None, - user_src=None, **kwds): + def __init__(self, topology, direction, platform_id=None, device_id=None, + device_type=None, gpu_precision=HYSOP_REAL): """ - Create the common attributes of all GPU discrete operators. - All in this interface is independant of a discrete operator. - - @param platform_id : OpenCL platform id (default = 0). - @param device_id : OpenCL device id (default = 0). - @param device_type : OpenCL device type (default = 'gpu'). - @param user_src : User OpenCL sources. + Parameters + ---------- + topology : :class:`~hysop.mpi.topology.Cartesian` + the topology on which this operator works + direction : int + leading direction of work (for example advection dir) + platform_id : int, optional + OpenCL id, default = 0. + device_id : int, optional + OpenCL id, default = 0. + device_type : string, optional. + OpenCL selected device, default = 'gpu'. + gpu_precision : numpy.dtype, optional + floatting precision for gpu, default=HYSOP_REAL """ - ## real type precision on GPU - self.gpu_precision = HYSOP_REAL - if 'method' in kwds and Precision in kwds['method'].keys(): - self.gpu_precision = kwds['method'][Precision] + super(GPUOperator, self).__init__() + # real type precision on GPU + self.gpu_precision = gpu_precision # Initialize opencl environment - comm_ref = self.variables[0].topology.comm - # use mpi_param comm instead? self.cl_env = get_opencl_environment( platform_id=platform_id, device_id=device_id, device_type=device_type, precision=self.gpu_precision, - comm=comm_ref) + comm=topology.comm) # Functions to get the appropriate vectors for the current direction - self.dim = self.domain.dimension self._reorderVect = lambda v: v - if self.dim == 2 and self.direction == 1: + dimension = topology.domain.dimension + if dimension == 2 and direction == 1: self._reorderVect = lambda v: (v[1], v[0]) - if self.dim == 3 and self.direction == 1: + if dimension == 3 and direction == 1: self._reorderVect = lambda v: (v[1], v[0], v[2]) - if self.dim == 3 and self.direction == 2: + if dimension == 3 and direction == 2: self._reorderVect = lambda v: (v[2], v[0], v[1]) # Size constants for local mesh size self._size_constants = "" self._kernel_cfg = \ - self.cl_env.kernels_config[self.dim][self.gpu_precision] + self.cl_env.kernels_config[dimension][self.gpu_precision] self._num_locMem = None - ## Global memory allocated on gpu by this operator + # Global memory allocated on gpu by this operator self.size_global_alloc = 0 - ## Local memory allocated on gpu by this operator + # Local memory allocated on gpu by this operator self.size_local_alloc = 0 def _append_size_constants(self, values, prefix='NB', suffix=None): - """ - Append to the string containing the constants for building kernels. - @param values : values to add - @param prefix : prefix of variables names - @param suffix : suffix of variables names + """Append to the string containing the constants for building kernels. + + Parameters + ---------- + values : list + prefix : string, optional + suffix : list of strings, optional + directions, default = `hysop.constants.S_DIR`. """ if suffix is None: suffix = S_DIR diff --git a/hysop/gpu/gpu_particle_advection.py b/hysop/gpu/gpu_particle_advection.py index 3b065ee716a03b9224cbb48a1e0eba2b83d5eb3b..8c2e736968caae213471ca268f179b8bc5e06636 100644 --- a/hysop/gpu/gpu_particle_advection.py +++ b/hysop/gpu/gpu_particle_advection.py @@ -1,11 +1,11 @@ """Discrete advection for GPU """ -from abc import abstractmethod +from abc import abstractmethod, ABCMeta from hysop import __VERBOSE__ from hysop.constants import np, debug, S_DIR, HYSOP_REAL from hysop.methods_keys import TimeIntegrator, Remesh, \ - Support, Splitting, MultiScale -from hysop.numerics.integrators.euler import Euler + Support, Splitting, MultiScale, Precision +from hysop.numerics.odesolvers import Euler from hysop.operator.discrete.particle_advection import ParticleAdvection from hysop.operator.discrete.discrete import get_extra_args_from_method from hysop.gpu import cl @@ -16,6 +16,7 @@ from hysop.gpu.gpu_discrete import GPUDiscreteField from hysop.gpu.gpu_operator import GPUOperator from hysop.tools.profiler import profile from hysop.numerics.update_ghosts import UpdateGhostsFull +from hysop.tools.misc import WorkSpaceTools class GPUParticleAdvection(ParticleAdvection, GPUOperator): @@ -48,128 +49,78 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): # init base class (i.e. ParticleAdvection) super(GPUParticleAdvection, self).__init__(**kwds) + + # mpi comm is set to advected field topology communicator. self.fields_topo = self.fields_on_grid[0].topology - self.velocity_topo = self.velocity.topology self._comm = self.fields_topo.comm self._comm_size = self._comm.Get_size() self._comm_rank = self._comm.Get_rank() + assert self._comm == self.velocity.topology.comm - user_src = get_extra_args_from_method(self, 'user_src', None) - - # init the second base class (the previous super only call - # the fist __init__ method found in the order of - # [ParticleAdvection, GPUOperator], i.e. ParticleAdvection.__init__) - # the GPUOperator.__init__ must be explicitely called. - # see http://stackoverflow.com/questions/3277367/how-does-pythons-super-work-with-multiple-inheritance # init second base class (i.e GPUOperator) + if Precision in self.method: + gpu_precision = self.method[Precision] + else: + gpu_precision = None GPUOperator.__init__( - self, + self, topology=self.velocity.topology, + direction=self.direction, platform_id=get_extra_args_from_method(self, 'platform_id', None), device_id=get_extra_args_from_method(self, 'device_id', None), device_type=get_extra_args_from_method(self, 'device_type', None), - user_src=user_src, - **kwds) + gpu_precision=gpu_precision) - ## Work arrays for fields on particles (cpu) - self.fields_on_part = None - - # The default is one kernel for all operations - self._is2kernel = False + # Choose between one kernel for all operations (default) ... + self._use_2_kernels = False if self.method[Support].find('gpu_2k') >= 0: - # different kernels for advection and remesh - self._is2kernel = True - - self._isMultiScale = False - if MultiScale in self.method: - if self.method[MultiScale] is not None: - self._isMultiScale = True + # ... or two different kernels for advection and remesh + self._use_2_kernels = True self._synchronize = None - if self._isMultiScale: + self._isMultiScale = False + if MultiScale in self.method and self.method[MultiScale] is not None: + self._isMultiScale = True self._synchronize = UpdateGhostsFull( self.velocity.topology, self.velocity.nb_components) # Compute resolutions for kernels for each direction. - ## Resolution of the local mesh but reoganized redarding - ## splitting direction: - ## direction X : XYZ - ## direction Y : YXZ - ## direction Z : ZYX in parallel, ZXY in sequentiel. - self.resol_dir = npw.dim_ones((self.dim,)) - self.v_resol_dir = npw.dim_ones((self.dim,)) + # Resolution of the local mesh but reoganized regarding + # splitting direction: + # direction X : XYZ + # direction Y : YXZ + # direction Z : ZYX in parallel, ZXY in sequentiel. + self.resol_dir = npw.dim_ones((self._dim,)) + self.v_resol_dir = npw.dim_ones((self._dim,)) shape = self.fields_topo.mesh.resolution - v_shape = self.velocity_topo.mesh.resolution + v_shape = self.velocity.topology.mesh.resolution # Local mesh resolution - resol = shape.copy() - self.resol_dir[:self.dim] = self._reorderVect(shape) - v_resol = v_shape.copy() - self.v_resol_dir[:self.dim] = self._reorderVect(v_shape) - - self._append_size_constants(resol) - self._append_size_constants(v_resol, prefix='V_NB') - self._append_size_constants( - [self.velocity_topo.ghosts()[self.direction]], - prefix='V_GHOSTS_NB', suffix=['']) - enum = ['I', 'II', 'III'] - self._append_size_constants( - self._reorderVect(['NB' + d for d in S_DIR[:self.dim]]), - prefix='NB_', suffix=enum[:self.dim]) - self._append_size_constants( - self._reorderVect(['V_NB' + d for d in S_DIR[:self.dim]]), - prefix='V_NB_', suffix=enum[:self.dim]) - - fields_topo = self.fields_topo - # Coordinates of the local origin - self._coord_min = npw.ones(4, dtype=self.gpu_precision) - self._coord_min[:self.dim] = fields_topo.mesh.origin + self.resol_dir[:self._dim] = self._reorderVect(shape) + self.v_resol_dir[:self._dim] = self._reorderVect(v_shape) + self._init_size_constants() - # Space step for fields - self._mesh_size = npw.ones(4, dtype=self.gpu_precision) - self._mesh_size[:self.dim] = self._reorderVect( - self.fields_topo.mesh.space_step) - - # Space step for velocity - self._v_mesh_size = npw.ones(4, dtype=self.gpu_precision) - self._v_mesh_size[:self.dim] = self._reorderVect( - self.velocity_topo.mesh.space_step) - - self._mesh_info = npw.ones((12, )) - self._mesh_info[:4] = self._mesh_size - self._mesh_info[4:8] = self._v_mesh_size - self._mesh_info[8] = self._coord_min[self.direction] - self._mesh_info[9] = 1. / self._mesh_size[0] - self._mesh_info[10] = 1. / self._v_mesh_size[0] - self._cl_mesh_info = cl.Buffer(self.cl_env.ctx, cl.mem_flags.READ_ONLY, - size=self._mesh_info.nbytes) - cl.enqueue_write_buffer(self.cl_env.queue, - self._cl_mesh_info, self._mesh_info).wait() - - assert self._coord_min.dtype == self.gpu_precision - assert self._mesh_size.dtype == self.gpu_precision - assert self._v_mesh_size.dtype == self.gpu_precision - - ## opencl kernels build options - self.build_options = "" + # Init mesh info + self._cl_mesh_info = None + self._init_cl_mesh_info() # user defined opencl sources self.prg = None - self._collect_usr_cl_src(user_src) + # get extra (opencl kernels) files from method, if any + self._collect_usr_cl_src() # Set copy kernel - self.copy = None - self._collect_kernels_cl_src_copy() + self.copy = self._collect_kernels_cl_src_copy() # Set transposition kernels self.transpose_xy, self.transpose_xy_r = None, None self.transpose_xz, self.transpose_xz_r = None, None self._collect_kernels_cl_src_transpositions_xy() - if self.dim == 3: + if self._dim == 3: self._collect_kernels_cl_src_transpositions_xz() # Set advection and remesh kernels self.num_advec, self.num_remesh = None, None self.num_advec_and_remesh = None - if self._is2kernel: + if self._use_2_kernels: self._collect_kernels_cl_src_2k() self._compute = self._compute_2k else: @@ -193,12 +144,76 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): # Particle initialisation OpenCL events for each field: self._init_events = {self.fields_on_grid[0]: []} - @abstractmethod - def globalMemoryUsagePreview(self, v_shape, shape): + def _init_size_constants(self): + """Fill size_constants attribute in + (_size_constant belongs to gpu_operator) """ - @param[in] v_shape: shape of the discretization of the velocity - @param[in] shape: shape of the discretization of the advected fields - @return size of the required memory + v_shape = self.velocity.topology.mesh.resolution + f_shape = self.fields_topo.mesh.resolution + self._append_size_constants(f_shape) + self._append_size_constants(v_shape, prefix='V_NB') + self._append_size_constants( + [self.velocity.topology.ghosts()[self.direction]], + prefix='V_GHOSTS_NB', suffix=['']) + enum = ['I', 'II', 'III'] + self._append_size_constants( + self._reorderVect(['NB' + d for d in S_DIR[:self._dim]]), + prefix='NB_', suffix=enum[:self._dim]) + self._append_size_constants( + self._reorderVect(['V_NB' + d for d in S_DIR[:self._dim]]), + prefix='V_NB_', suffix=enum[:self._dim]) + + def _init_cl_mesh_info(self): + """Collect mesh info from fields and velocity, + set and send opencl buffer (self._cl_mesh_info) to device. + """ + # Space step for fields + mesh_size = npw.ones(4, dtype=self.gpu_precision) + mesh_size[:self._dim] = self._reorderVect( + self.fields_topo.mesh.space_step) + + # Space step for velocity + self._v_mesh_size = npw.ones(4, dtype=self.gpu_precision) + self._v_mesh_size[:self._dim] = self._reorderVect( + self.velocity.topology.mesh.space_step) + + mesh_info = npw.ones((12, ), dtype=self.gpu_precision) + mesh_info[:4] = mesh_size + mesh_info[4:8] = self._v_mesh_size + # Coordinate of the local origin in advection dir. + mesh_info[8] = self.fields_topo.mesh.origin[self.direction] + mesh_info[9] = 1. / mesh_size[0] + mesh_info[10] = 1. / self._v_mesh_size[0] + assert mesh_size.dtype == self.gpu_precision + # returns an opencl buffer + self._cl_mesh_info = cl.Buffer(self.cl_env.ctx, cl.mem_flags.READ_ONLY, + size=mesh_info.nbytes) + cl.enqueue_write_buffer(self.cl_env.queue, + self._cl_mesh_info, mesh_info).wait() + + def _set_work_arrays(self, rwork=None, iwork=None): + + topo = self.fields_on_grid[0].topology + # Find number and shape of required work arrays + # For GPU version, no need of numerics works + # Shape of reference comes from fields, not from velocity + rwork_length = np.sum([f.nb_components for f in self.fields_on_grid]) + if self.method[Support].find('gpu_2k') >= 0: + rwork_length += 1 # work array for positions + + # check and/or allocate work arrays according to properties above + subshape = tuple(topo.mesh.resolution) + self._rwork = WorkSpaceTools.check_work_array(rwork_length, subshape, + rwork, HYSOP_REAL) + + @abstractmethod + def global_memory_usage(self, v_shape, shape): + """Returns an estimation of memory usage + + Parameters + ---------- + v_shape, shape : tuples + shapes of the velocity and advected fields """ pass @@ -210,13 +225,13 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): # 3D: X(dt/2), Y(dt/2), Z(dt), Y(dt/2), X(dt/2) # 2D: X(dt/2), Y(dt), X(dt/2) if self.method[Splitting] == 'o2': - if self.dim == 2: + if self._dim == 2: self.exec_list = [ [self._init_copy, self._compute], # X(dt/2) [self._init_transpose_xy, self._compute], # Y(dt) [self._init_transpose_xy, self._compute] # X(dt/2) ] - elif self.dim == 3: + elif self._dim == 3: self.exec_list = [ [self._init_copy, self._compute], # X(dt/2) [self._init_transpose_xy, self._compute], # Y(dt/2) @@ -228,14 +243,14 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): # Splitting Strang 2nd order (fullHalf): # X(dt/2), Y(dt/2), Z(dt/2), Z(dt/2), Y(dt/2), X(dt/2) elif self.method[Splitting] == 'o2_FullHalf': - if self.dim == 2: + if self._dim == 2: self.exec_list = [ [self._init_copy, self._compute], # X(dt/2) [self._init_transpose_xy, self._compute], # Y(dt) [self._init_copy, self._compute], # Y(dt) [self._init_transpose_xy, self._compute] # X(dt/2) ] - elif self.dim == 3: + elif self._dim == 3: self.exec_list = [ [self._init_copy, self._compute], # X(dt/2) [self._init_transpose_xy, self._compute], # Y(dt/2) @@ -253,8 +268,8 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): raise ValueError('Splitting type not yet implemeted on GPU: ' + self.method[Splitting]) - def globalMemoryUsagePreview(self, v_shape, shape): - if self._is2kernel: + def global_memory_usage(self, v_shape, shape): + if self._use_2_kernels: r = (self.velocity.nb_components * v_shape.prod() + (2 * self.fields_on_grid[0].nb_components + 1) * shape.prod()) else: @@ -285,19 +300,18 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): self.size_global_alloc += self.fields_on_grid[0].mem_size # Fields on particles - self.fields_on_part = {} start = 0 for f in self.fields_on_grid: for i in xrange(start, start + f.nb_components): - if type(self._rwork[i]) is np.ndarray: + if isinstance(self._rwork[i], np.ndarray): self._rwork[i] = \ self.cl_env.global_allocation(self._rwork[i]) self.fields_on_part[f] = self._rwork[start: start + f.nb_components] start += f.nb_components - if self._is2kernel: - ## Particles position - if type(self._rwork[start]) is np.ndarray: + if self._use_2_kernels: + # Particles position + if isinstance(self._rwork[start], np.ndarray): self._rwork[start] = \ self.cl_env.global_allocation(self._rwork[start]) self.part_position = self._rwork[start:start + 1] @@ -332,7 +346,6 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): """ Compile OpenCL sources for copy kernel. """ - # build_options = self.build_options # # copy settings # src, t_dim, b_rows, vec, f_space = self._kernel_cfg['copy'] # while t_dim > self.resol_dir[0] or (self.resol_dir[0] % t_dim) > 0: @@ -349,50 +362,51 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): # vec) # self.copy = KernelLauncher(prg.copy, # self.cl_env.queue, gwi, lwi) - self.copy = KernelLauncher(cl.enqueue_copy, - self.cl_env.queue) + return KernelLauncher(cl.enqueue_copy, self.cl_env.queue) def _collect_kernels_cl_src_transpositions_xy(self): """Compile OpenCL sources for transpositions kernel. - @remark : Transpositions kernels are launched as initialization. - Arrays are taken to their destination layout (for initialize in Y - directions, either we came from X or Z but shapes are the Y ones). + + Notes + ----- + + * Transpositions kernels are launched at initialization. This routine sets transpose_xy and transpose_xy_r. """ resol = self.fields_topo.mesh.resolution - resol_tmp = npw.zeros_like(resol) # XY transposition settings is_XY_needed = self.direction == 1 or self.direction == 0 if is_XY_needed: - resol_tmp[...] = resol[...] + resol_tmp = resol.copy() if self.direction == 1: # (XY -> YX) resol_tmp[0] = resol[1] resol_tmp[1] = resol[0] ocl_cte = " -D NB_I=NB_Y -D NB_II=NB_X -D NB_III=NB_Z" - elif self.direction == 0: # (YX -> XY) only for sequential + else: + # self.direction == 0: # (YX -> XY) only for sequential ocl_cte = " -D NB_I=NB_X -D NB_II=NB_Y -D NB_III=NB_Z" self.transpose_xy = self._make_transpose_xy(resol_tmp, ocl_cte) - # is_XY_r_needed = self.direction == 1 and self._comm_size > 1 - # if is_XY_r_needed: - # # Reversed XY transposition settings (YX -> XY), only in parallel - # resol_tmp[...] = resol[...] - # ocl_cte = " -D NB_I=NB_X -D NB_II=NB_Y -D NB_III=NB_Z" + def _make_transpose_xy(self, resolution, ocl_cte): + """Perform xy transposition - # self.transpose_xy_r = self._make_transpose_xy(resol_tmp, ocl_cte) - - def _make_transpose_xy(self, resol_tmp, ocl_cte): - - build_options = self.build_options + self._size_constants + Parameters + ---------- + resolution : numpy array + required shape for? + ocl_cte: string + compilation options + """ + build_options = self._size_constants src, t_dim, b_rows, is_padding, vec, f_space = \ self._kernel_cfg['transpose_xy'] - while t_dim > resol_tmp[0] or t_dim > resol_tmp[1] or \ - (resol_tmp[0] % t_dim) > 0 or (resol_tmp[1] % t_dim) > 0: + while t_dim > resolution[0] or t_dim > resolution[1] or \ + (resolution[0] % t_dim) > 0 or (resolution[1] % t_dim) > 0: t_dim /= 2 - gwi, lwi, blocs_nb = f_space(resol_tmp, t_dim, b_rows, vec) + gwi, lwi, blocs_nb = f_space(resolution, t_dim, b_rows, vec) if is_padding: build_options += " -D PADDING_XY=1" @@ -441,7 +455,7 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): def _make_transpose_xz(self, resol_tmp, ocl_cte): - build_options = self.build_options + self._size_constants + build_options = self._size_constants src, t_dim, b_rows, b_deph, is_padding, vec, f_space = \ self._kernel_cfg['transpose_xz'] @@ -465,25 +479,26 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): vec) return KernelLauncher(prg.transpose_xz, self.cl_env.queue, gwi, lwi) - def _collect_usr_cl_src(self, usr_src): - """ - Build user sources. + def _collect_usr_cl_src(self): + """Build user sources. """ - if usr_src is not None: - build_options = self.build_options + self._size_constants + # get extra (opencl kernels) files from method, if any + user_src = get_extra_args_from_method(self, 'user_src', None) + if user_src is not None: + build_options = self._size_constants workItemNb, gwi, lwi = self.cl_env.get_work_items(self.resol_dir) v_workItemNb, gwi, lwi = self.cl_env.get_work_items(self.v_resol_dir) build_options += " -D WI_NB=" + str(workItemNb) build_options += " -D V_WI_NB=" + str(v_workItemNb) - self.prg = self.cl_env.build_src(usr_src, build_options, 1) + self.prg = self.cl_env.build_src(user_src, build_options, 1) def _collect_kernels_cl_src_1k(self): """ Compile OpenCL sources for advection and remeshing kernel. """ - build_options = self.build_options + self._size_constants + build_options = self._size_constants src, is_noBC, vec, f_space = self._kernel_cfg['advec_and_remesh'] gwi, lwi = f_space(self.resol_dir, vec) WINb = lwi[0] @@ -516,7 +531,7 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): Compile OpenCL sources for advection and remeshing kernel. """ # Advection - build_options = self.build_options + self._size_constants + build_options = self._size_constants src, is_noBC, vec, f_space = self._kernel_cfg['advec'] gwi, lwi = f_space(self.resol_dir, vec) WINb = lwi[0] @@ -554,7 +569,7 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): prg.advection_kernel, self.cl_env.queue, gwi, lwi) # remeshing - build_options = self.build_options + self._size_constants + build_options = self._size_constants src, is_noBC, vec, f_space = self._kernel_cfg['remesh'] gwi, lwi = f_space(self.resol_dir, vec) WINb = lwi[0] @@ -574,7 +589,7 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): @debug @profile - def apply(self, simulation, dtCoeff, split_id, old_dir): + def apply(self, simulation, dt_coeff, split_id, old_dir): """ Apply operator along specified splitting direction. @param t : Current time @@ -583,15 +598,15 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): @param split_id : Splitting step id """ # If first direction of advection, wait for work gpu fields - # It avoid wait_for lists to increase indelinitely + # It avoid wait_for lists to increase indefinitely # In practice, all events are terminated so wait() resets events list if split_id == 0: for v in self.fields_on_grid + [self.velocity]: v.clean_events() for exe in self.exec_list[split_id]: - exe(simulation, dtCoeff, split_id, old_dir) + exe(simulation, dt_coeff, split_id, old_dir) - def _init_copy(self, simulation, dtCoeff, split_id, old_dir): + def _init_copy(self, simulation, dt_coeff, split_id, old_dir): wait_evt = self.fields_on_grid[0].events for g, p in zip(self.fields_on_grid[0].gpu_data, self.fields_on_part[self.fields_on_grid[0]]): @@ -599,7 +614,7 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): #evt = self.copy(g, p, wait_for=wait_evt) self._init_events[self.fields_on_grid[0]].append(evt) - # def _init_copy_r(self, simulation, dtCoeff, split_id, old_dir): + # def _init_copy_r(self, simulation, dt_coeff, split_id, old_dir): # wait_evt = self.fields_on_grid[0].events # for g, p in zip(self.fields_on_grid[0].gpu_data, # self.fields_on_part[self.fields_on_grid[0]]): @@ -607,36 +622,36 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): # #evt = self.copy(p, g, wait_for=wait_evt) # self._init_events[self.fields_on_grid[0]].append(evt) - def _init_transpose_xy(self, simulation, dtCoeff, split_id, old_dir): + def _init_transpose_xy(self, simulation, dt_coeff, split_id, old_dir): wait_evt = self.fields_on_grid[0].events for g, p in zip(self.fields_on_grid[0].gpu_data, self.fields_on_part[self.fields_on_grid[0]]): evt = self.transpose_xy(g, p, wait_for=wait_evt) self._init_events[self.fields_on_grid[0]].append(evt) - # def _init_transpose_xy_r(self, simulation, dtCoeff, split_id, old_dir): + # def _init_transpose_xy_r(self, simulation, dt_coeff, split_id, old_dir): # wait_evt = self.fields_on_grid[0].events # for g, p in zip(self.fields_on_grid[0].gpu_data, # self.fields_on_part[self.fields_on_grid[0]]): # evt = self.transpose_xy_r(p, g, wait_for=wait_evt) # self._init_events[self.fields_on_grid[0]].append(evt) - def _init_transpose_xz(self, simulation, dtCoeff, split_id, old_dir): + def _init_transpose_xz(self, simulation, dt_coeff, split_id, old_dir): wait_evt = self.fields_on_grid[0].events for g, p in zip(self.fields_on_grid[0].gpu_data, self.fields_on_part[self.fields_on_grid[0]]): evt = self.transpose_xz(g, p, wait_for=wait_evt) self._init_events[self.fields_on_grid[0]].append(evt) - # def _init_transpose_xz_r(self, simulation, dtCoeff, split_id, old_dir): + # def _init_transpose_xz_r(self, simulation, dt_coeff, split_id, old_dir): # wait_evt = self.fields_on_grid[0].events # for g, p in zip(self.fields_on_grid[0].gpu_data, # self.fields_on_part[self.fields_on_grid[0]]): # evt = self.transpose_xz_r(p, g, wait_for=wait_evt) # self._init_events[self.fields_on_grid[0]].append(evt) - def _compute_advec_euler_simpleechelle(self, simulation, dtCoeff, split_id, old_dir): - dt = simulation.time_step * dtCoeff + def _compute_advec_euler_simpleechelle(self, simulation, dt_coeff, split_id, old_dir): + dt = simulation.time_step * dt_coeff wait_evts = self.velocity.events + \ self._init_events[self.fields_on_grid[0]] # Advection @@ -648,8 +663,8 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): wait_for=wait_evts) self._init_events[self.fields_on_grid[0]].append(evt) - def _compute_advec_simpleechelle(self, simulation, dtCoeff, split_id, old_dir): - dt = simulation.time_step * dtCoeff + def _compute_advec_simpleechelle(self, simulation, dt_coeff, split_id, old_dir): + dt = simulation.time_step * dt_coeff wait_evts = self.velocity.events + \ self._init_events[self.fields_on_grid[0]] # Advection @@ -661,8 +676,8 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): wait_for=wait_evts) self._init_events[self.fields_on_grid[0]].append(evt) - def _compute_advec_multiechelle(self, simulation, dtCoeff, split_id, old_dir): - dt = simulation.time_step * dtCoeff + def _compute_advec_multiechelle(self, simulation, dt_coeff, split_id, old_dir): + dt = simulation.time_step * dt_coeff wait_evts = self.velocity.events + \ self._init_events[self.fields_on_grid[0]] # Advection @@ -676,8 +691,8 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): wait_for=wait_evts) self._init_events[self.fields_on_grid[0]].append(evt) - def _compute_2k(self, simulation, dtCoeff, split_id, old_dir): - self._compute_advec(simulation, dtCoeff, split_id, old_dir) + def _compute_2k(self, simulation, dt_coeff, split_id, old_dir): + self._compute_advec(simulation, dt_coeff, split_id, old_dir) wait_evts = self._init_events[self.fields_on_grid[0]] + \ self.fields_on_grid[0].events nbc = self.fields_on_grid[0].nb_components @@ -691,11 +706,11 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): self.fields_on_grid[0].events.append(evt) self._init_events[self.fields_on_grid[0]] = [] - def _compute_1k_multiechelle(self, simulation, dtCoeff, split_id, old_dir): + def _compute_1k_multiechelle(self, simulation, dt_coeff, split_id, old_dir): if split_id==0 and self._synchronize is not None: self._synchronize(self.velocity.data) self.velocity.toDevice() - dt = simulation.time_step * dtCoeff + dt = simulation.time_step * dt_coeff wait_evts = self.velocity.events + \ self._init_events[self.fields_on_grid[0]] + \ self.fields_on_grid[0].events @@ -713,8 +728,8 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): self.fields_on_grid[0].events.append(evt) self._init_events[self.fields_on_grid[0]] = [] - def _compute_1k_simpleechelle(self, simulation, dtCoeff, split_id, old_dir): - dt = simulation.time_step * dtCoeff + def _compute_1k_simpleechelle(self, simulation, dt_coeff, split_id, old_dir): + dt = simulation.time_step * dt_coeff wait_evts = self.velocity.events + \ self._init_events[self.fields_on_grid[0]] + \ self.fields_on_grid[0].events @@ -730,8 +745,8 @@ class GPUParticleAdvection(ParticleAdvection, GPUOperator): self._init_events[self.fields_on_grid[0]] = [] - def _compute_1k_euler_simpleechelle(self, simulation, dtCoeff, split_id, old_dir): - dt = simulation.time_step * dtCoeff + def _compute_1k_euler_simpleechelle(self, simulation, dt_coeff, split_id, old_dir): + dt = simulation.time_step * dt_coeff wait_evts = self.velocity.events + \ self._init_events[self.fields_on_grid[0]] + \ self.fields_on_grid[0].events diff --git a/hysop/gpu/multi_gpu_particle_advection.py b/hysop/gpu/multi_gpu_particle_advection.py index 41023d6b1ee8fed939d7abdd38bca0d49bb99b33..7fb22bd425bcaed81bb9546cb80667db959d2cfb 100644 --- a/hysop/gpu/multi_gpu_particle_advection.py +++ b/hysop/gpu/multi_gpu_particle_advection.py @@ -8,11 +8,11 @@ from hysop.constants import np, debug, HYSOP_INTEGER, HYSOP_REAL, ORDER,\ HYSOP_MPI_REAL, SIZEOF_HYSOP_REAL from hysop.gpu.gpu_particle_advection import GPUParticleAdvection from hysop.operator.discrete.discrete import get_extra_args_from_method -from hysop.methods_keys import Support, TimeIntegrator, MultiScale, Remesh -from hysop.numerics.integrators.runge_kutta2 import RK2 +from hysop.methods_keys import TimeIntegrator, MultiScale, Remesh +from hysop.numerics.odesolvers import RK2 from hysop.numerics.remeshing import Linear as Linear_rmsh from hysop.gpu.gpu_kernel import KernelLauncher -from hysop.tools.profiler import FProfiler, profile +from hysop.tools.profiler import FProfiler from hysop.gpu import cl, CL_PROFILE from hysop.mpi import Wtime import hysop.tools.numpywrappers as npw @@ -679,8 +679,8 @@ class MultiGPUParticleAdvection(GPUParticleAdvection): advec_gpu_time += (evt.profile.end - evt.profile.start) * 1e-9 self.profiler['comm_gpu_advec_set'] += advec_gpu_time - def _init_copy(self, simulation, dtCoeff, split_id, old_dir): - self._exchange_velocity_buffers(simulation.time_step * dtCoeff) + def _init_copy(self, simulation, dt_coeff, split_id, old_dir): + self._exchange_velocity_buffers(simulation.time_step * dt_coeff) wait_evt = self.fields_on_grid[0].events for g, p in zip(self.fields_on_grid[0].gpu_data, self.fields_on_part[self.fields_on_grid[0]]): @@ -688,24 +688,25 @@ class MultiGPUParticleAdvection(GPUParticleAdvection): #evt = self.copy(g, p, wait_for=wait_evt) self._init_events[self.fields_on_grid[0]].append(evt) - def _init_transpose_xy(self, simulation, dtCoeff, split_id, old_dir): - self._exchange_velocity_buffers(simulation.time_step * dtCoeff) + def _init_transpose_xy(self, simulation, dt_coeff, split_id, old_dir): + self._exchange_velocity_buffers(simulation.time_step * dt_coeff) wait_evt = self.fields_on_grid[0].events for g, p in zip(self.fields_on_grid[0].gpu_data, self.fields_on_part[self.fields_on_grid[0]]): evt = self.transpose_xy(g, p, wait_for=wait_evt) self._init_events[self.fields_on_grid[0]].append(evt) - def _init_transpose_xz(self, simulation, dtCoeff, split_id, old_dir): - self._exchange_velocity_buffers(simulation.time_step * dtCoeff) + def _init_transpose_xz(self, simulation, dt_coeff, split_id, old_dir): + self._exchange_velocity_buffers(simulation.time_step * dt_coeff) wait_evt = self.fields_on_grid[0].events for g, p in zip(self.fields_on_grid[0].gpu_data, self.fields_on_part[self.fields_on_grid[0]]): evt = self.transpose_xz(g, p, wait_for=wait_evt) self._init_events[self.fields_on_grid[0]].append(evt) - def _compute_advec_comm(self, simulation, dtCoeff, split_id, old_dir): - dt = simulation.time_step * dtCoeff + def _compute_advec_comm(self, simulation, dt_coeff, split_id, old_dir): + dt = simulation.time_step * dt_coeff + self._todevice_velocity_buffers() wait_evts = self.velocity.events + self._evt_l_v + self._evt_r_v + \ self._init_events[self.fields_on_grid[0]] @@ -829,10 +830,10 @@ class MultiGPUParticleAdvection(GPUParticleAdvection): self._cl_mesh_info, wait_for=wait_list) - def _compute_1c_comm(self, simulation, dtCoeff, split_id, old_dir): - dt = simulation.time_step * dtCoeff + def _compute_1c_comm(self, simulation, dt_coeff, split_id, old_dir): + dt = simulation.time_step * dt_coeff if self._is2kernel: - self._compute_advec_comm(simulation, dtCoeff, split_id, old_dir) + self._compute_advec_comm(simulation, dt_coeff, split_id, old_dir) else: self._todevice_velocity_buffers() wait_evts = self.velocity.events + \ diff --git a/hysop/gpu/tests/test_advection_nullVelocity.py b/hysop/gpu/tests/test_advection_nullVelocity.py deleted file mode 100644 index 5a26b98c2a6c57a56161730c74b51f25f1aef94e..0000000000000000000000000000000000000000 --- a/hysop/gpu/tests/test_advection_nullVelocity.py +++ /dev/null @@ -1,992 +0,0 @@ -""" -@file hysop.gpu.tests.test_advection_nullVelocity -Testing advection kernels with a null velocity. Basic functionnal test. -""" -from hysop.domain.box import Box -from hysop.fields.continuous import Field -from hysop.operator.advection import Advection -from hysop.constants import np, HYSOP_REAL -from hysop.problem.simulation import Simulation -from hysop.methods_keys import TimeIntegrator, Interpolation, Remesh, \ - Support, Splitting, Precision -from hysop.numerics.integrators.runge_kutta2 import RK2 -from hysop.numerics.interpolation import Linear -from hysop.numerics.remeshing import L2_1, L4_2, L6_3, M8Prime -from hysop.tools.parameters import Discretization -import hysop.tools.numpywrappers as npw - - -def setup_2D(): - box = Box(length=[1., 1.], origin=[0., 0.]) - scal = Field(domain=box, name='Scalar') - velo = Field(domain=box, name='Velocity', is_vector=True) - return scal, velo - - -def setup_3D(): - box = Box(length=[1., 1., 1.], origin=[0., 0., 0.]) - scal = Field(domain=box, name='Scalar') - velo = Field(domain=box, name='Velocity', is_vector=True) - return scal, velo - - -def assertion_2D(scal, velo, advec): - advec.discretize() - advec.setup() - - scal_d = scal.discreteFields.values()[0] - velo_d = velo.discreteFields.values()[0] - scal_d.data[0][...] = npw.asrealarray(np.random.random(scal_d.data[0].shape)) - velo_d.data[0][...] = npw.zeros_like(scal_d.data[0]) - velo_d.data[1][...] = npw.zeros_like(scal_d.data[0]) - scal_init = scal_d.data[0].copy() - scal_d.toDevice() - velo_d.toDevice() - scal_d.wait() - velo_d.wait() - - advec.apply(Simulation(start=0., end=0.01, nb_iter=1)) - - scal_d.toHost() - scal_d.wait() - advec.finalize() - return np.allclose(scal_init, scal_d.data[0]) - - -def assertion_2D_withPython(scal, velo, advec, advec_py): - advec.discretize() - advec_py.discretize() - advec.setup() - advec_py.setup() - - scal_d = scal.discreteFields.values()[0] - velo_d = velo.discreteFields.values()[0] - scal_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) - velo_d.data[0][...] = npw.zeros(scal_d.data[0].shape) - velo_d.data[1][...] = npw.zeros(scal_d.data[0].shape) - scal_d.toDevice() - velo_d.toDevice() - scal_d.wait() - velo_d.wait() - - advec_py.apply(Simulation(start=0., end=0.01, nb_iter=1)) - advec.apply(Simulation(start=0., end=0.01, nb_iter=1)) - - py_res = scal_d.data[0].copy() - scal_d.toHost() - scal_d.wait() - - advec.finalize() - print py_res, scal_d.data[0] - print py_res - scal_d.data[0] - return np.allclose(py_res, scal_d.data[0]) - - -def assertion_3D(scal, velo, advec): - advec.discretize() - advec.setup() - - scal_d = scal.discreteFields.values()[0] - velo_d = velo.discreteFields.values()[0] - scal_d.data[0][...] = npw.asrealarray(np.random.random(scal_d.data[0].shape)) - velo_d.data[0][...] = npw.zeros_like(scal_d.data[0]) - velo_d.data[1][...] = npw.zeros_like(scal_d.data[0]) - velo_d.data[2][...] = npw.zeros_like(scal_d.data[0]) - scal_init = scal_d.data[0].copy() - scal_d.toDevice() - velo_d.toDevice() - scal_d.wait() - velo_d.wait() - - advec.apply(Simulation(start=0., end=0.01, nb_iter=1)) - - scal_d.toHost() - scal_d.wait() - - advec.finalize() - return np.allclose(scal_init, scal_d.data[0]) - - -def assertion_3D_withPython(scal, velo, advec, advec_py): - advec.discretize() - advec_py.discretize() - advec.setup() - advec_py.setup() - - scal_d = scal.discreteFields.values()[0] - velo_d = velo.discreteFields.values()[0] - scal_d.data[0][...] = npw.asrealarray(np.random.random(scal_d.data[0].shape)) - velo_d.data[0][...] = npw.zeros_like(scal_d.data[0]) - velo_d.data[1][...] = npw.zeros_like(scal_d.data[0]) - velo_d.data[2][...] = npw.zeros_like(scal_d.data[0]) - scal_d.toDevice() - velo_d.toDevice() - scal_d.wait() - velo_d.wait() - - advec_py.apply(Simulation(start=0., end=0.01, nb_iter=1)) - advec.apply(Simulation(start=0., end=0.01, nb_iter=1)) - - py_res = scal_d.data[0].copy() - scal_d.toHost() - scal_d.wait() - - advec.finalize() - return np.allclose(py_res, scal_d.data[0]) - - -d2d = Discretization([33,33]) -d3d = Discretization([17, 17, 17]) - -# M6 tests -def test_2D_m6_1k(): - """ - Testing M6 remeshing formula in 2D, 1 kernel, - o2 splitting. - """ - scal, velo = setup_2D() - - advec = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: 'gpu_1k', - Splitting: 'o2', - Precision: HYSOP_REAL} - ) - advec_py = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: '', - Splitting: 'o2'}, - ) - assert assertion_2D_withPython(scal, velo, advec, advec_py) - #assert assertion_2D(scal, velo, advec) - - -def test_2D_m6_2k(): - """ - Testing M6 remeshing formula in 2D, 2 kernel, - o2 splitting. - """ - scal, velo = setup_2D() - - advec = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: 'gpu_2k', - Splitting: 'o2', - Precision: HYSOP_REAL} - ) - assert assertion_2D(scal, velo, advec) - - -def test_2D_m6_1k_sFH(): - """ - Testing M6 remeshing formula in 2D, 1 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_2D() - - advec = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: 'gpu_1k', - Splitting: 'o2_FullHalf', - Precision: HYSOP_REAL} - ) - assert assertion_2D(scal, velo, advec) - - -def test_2D_m6_2k_sFH(): - """ - Testing M6 remeshing formula in 2D, 2 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_2D() - - advec = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: 'gpu_2k', - Splitting: 'o2_FullHalf', - Precision: HYSOP_REAL} - ) - assert assertion_2D(scal, velo, advec) - - -def test_3D_m6_1k(): - """ - Testing M6 remeshing formula in 3D, 1 kernel, - o2 splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: 'gpu_1k', - Splitting: 'o2', - Precision: HYSOP_REAL} - ) - assert assertion_3D(scal, velo, advec) - - -def test_3D_m6_2k(): - """ - Testing M6 remeshing formula in 3D, 2 kernel, - o2 splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: 'gpu_2k', - Splitting: 'o2', - Precision: HYSOP_REAL} - ) - assert assertion_3D(scal, velo, advec) - - -def test_3D_m6_1k_sFH(): - """ - Testing M6 remeshing formula in 3D, 1 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: 'gpu_1k', - Splitting: 'o2_FullHalf', - Precision: HYSOP_REAL} - ) - assert assertion_3D(scal, velo, advec) - - -def test_3D_m6_2k_sFH(): - """ - Testing M6 remeshing formula in 3D, 2 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: 'gpu_2k', - Splitting: 'o2_FullHalf', - Precision: HYSOP_REAL}, - ) - assert assertion_3D(scal, velo, advec) - - -# M4 testing -def test_2D_m4_1k(): - """ - Testing M4 remeshing formula in 2D, 1 kernel, - o2 splitting. - """ - scal, velo = setup_2D() - - advec = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: 'gpu_1k', - Splitting: 'o2', - Precision: HYSOP_REAL} - ) - advec_py = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: '', - Splitting: 'o2'} - ) - assert assertion_2D_withPython(scal, velo, advec, advec_py) - - -def test_2D_m4_2k(): - """ - Testing M4 remeshing formula in 2D, 2 kernel, - o2 splitting. - """ - scal, velo = setup_2D() - - advec = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: 'gpu_2k', - Splitting: 'o2', - Precision: HYSOP_REAL - } - ) - advec_py = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: '', - Splitting: 'o2'}, - ) - assert assertion_2D_withPython(scal, velo, advec, advec_py) - - -def test_2D_m4_1k_sFH(): - """ - Testing M4 remeshing formula in 2D, 1 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_2D() - - advec = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: 'gpu_1k', - Splitting: 'o2_FullHalf', - Precision: HYSOP_REAL}, - ) - advec_py = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: '', - Splitting: 'o2_FullHalf'} - ) - assert assertion_2D_withPython(scal, velo, advec, advec_py) - - -def test_2D_m4_2k_sFH(): - """ - Testing M4 remeshing formula in 2D, 2 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_2D() - - advec = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: 'gpu_2k', - Splitting: 'o2_FullHalf'} - ) - advec_py = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: '', - Splitting: 'o2_FullHalf'} - ) - assert assertion_2D_withPython(scal, velo, advec, advec_py) - - -def test_3D_m4_1k(): - """ - Testing M4 remeshing formula in 3D, 1 kernel, - o2 splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: 'gpu_1k', - Splitting: 'o2'} - ) - advec_py = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: '', - Splitting: 'o2'}, - ) - assert assertion_3D_withPython(scal, velo, advec, advec_py) - - -def test_3D_m4_2k(): - """ - Testing M4 remeshing formula in 3D, 2 kernel, - o2 splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: 'gpu_2k', - Splitting: 'o2'} - ) - advec_py = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: '', - Splitting: 'o2'}, - ) - assert assertion_3D_withPython(scal, velo, advec, advec_py) - - -def test_3D_m4_1k_sFH(): - """ - Testing M4 remeshing formula in 3D, 1 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: 'gpu_1k', - Splitting: 'o2_FullHalf'}) - advec_py = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: '', - Splitting: 'o2_FullHalf'} - ) - assert assertion_3D_withPython(scal, velo, advec, advec_py) - - -def test_3D_m4_2k_sFH(): - """ - Testing M4 remeshing formula in 3D, 2 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: 'gpu_2k', - Splitting: 'o2_FullHalf'} - ) - advec_py = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Support: '', - Splitting: 'o2_FullHalf'} - ) - assert assertion_3D_withPython(scal, velo, advec, advec_py) - - -# M8 testing -def test_2D_m8_1k(): - """ - Testing M8 remeshing formula in 2D, 1 kernel, - o2 splitting. - """ - scal, velo = setup_2D() - - advec = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: 'gpu_1k', - Splitting: 'o2'} - ) - advec_py = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: '', - Splitting: 'o2'}, - ) - assert assertion_2D_withPython(scal, velo, advec, advec_py) - - -def test_2D_m8_2k(): - """ - Testing M8 remeshing formula in 2D, 2 kernel, - o2 splitting. - """ - scal, velo = setup_2D() - - advec = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: 'gpu_2k', - Splitting: 'o2'} - ) - advec_py = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: '', - Splitting: 'o2'}, - ) - assert assertion_2D_withPython(scal, velo, advec, advec_py) - - -def test_2D_m8_1k_sFH(): - """ - Testing M8 remeshing formula in 2D, 1 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_2D() - - advec = Advection(velo, scal, discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: 'gpu_1k', - Splitting: 'o2_FullHalf'} - ) - advec_py = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: '', - Splitting: 'o2_FullHalf'} - ) - assert assertion_2D_withPython(scal, velo, advec, advec_py) - - -def test_2D_m8_2k_sFH(): - """ - Testing M8 remeshing formula in 2D, 2 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_2D() - - advec = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: 'gpu_2k', - Splitting: 'o2_FullHalf'} - ) - advec_py = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: '', - Splitting: 'o2_FullHalf'} - ) - assert assertion_2D_withPython(scal, velo, advec, advec_py) - - -def test_3D_m8_1k(): - """ - Testing M8 remeshing formula in 3D, 1 kernel, - o2 splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: 'gpu_1k', - Splitting: 'o2'} - ) - advec_py = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: '', - Splitting: 'o2'}, - ) - assert assertion_3D_withPython(scal, velo, advec, advec_py) - - -def test_3D_m8_2k(): - """ - Testing M8 remeshing formula in 3D, 2 kernel, - o2 splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: 'gpu_2k', - Splitting: 'o2'} - ) - advec_py = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: '', - Splitting: 'o2'}, - ) - assert assertion_3D_withPython(scal, velo, advec, advec_py) - - -def test_3D_m8_1k_sFH(): - """ - Testing M8 remeshing formula in 3D, 1 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: 'gpu_1k', - Splitting: 'o2_FullHalf'} - ) - advec_py = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: '', - Splitting: 'o2_FullHalf'} - ) - assert assertion_3D_withPython(scal, velo, advec, advec_py) - - -def test_3D_m8_2k_sFH(): - """ - Testing M8 remeshing formula in 3D, 2 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: 'gpu_2k', - Splitting: 'o2_FullHalf'} - ) - advec_py = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: '', - Splitting: 'o2_FullHalf'} - ) - assert assertion_3D_withPython(scal, velo, advec, advec_py) - - -def test_2D_l6_2k(): - """ - Testing Lambda6star remeshing formula in 2D, 2 kernel, - o2 splitting. - """ - scal, velo = setup_2D() - - advec = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L6_3, - Support: 'gpu_1k', - Splitting: 'o2'} - ) - assert assertion_2D(scal, velo, advec) - - -def test_2D_l6_1k_sFH(): - """ - Testing Lambda6star remeshing formula in 2D, 1 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_2D() - - advec = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L6_3, - Support: 'gpu_1k', - Splitting: 'o2_FullHalf'} - ) - assert assertion_2D(scal, velo, advec) - - -def test_2D_l6_2k_sFH(): - """ - Testing Lambda6star remeshing formula in 2D, 2 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_2D() - - advec = Advection(velo, scal,discretization=d2d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L6_3, - Support: 'gpu_2k', - Splitting: 'o2_FullHalf'} - ) - assert assertion_2D(scal, velo, advec) - - -def test_3D_l6_1k(): - """ - Testing Lambda6star remeshing formula in 3D, 1 kernel, - o2 splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L6_3, - Support: 'gpu_1k', - Splitting: 'o2'} - ) - assert assertion_3D(scal, velo, advec) - - -def test_3D_l6_2k(): - """ - Testing Lambda6star remeshing formula in 3D, 2 kernel, - o2 splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L6_3, - Support: 'gpu_2k', - Splitting: 'o2'} - ) - assert assertion_3D(scal, velo, advec) - - -def test_3D_l6_1k_sFH(): - """ - Testing Lambda6star remeshing formula in 3D, 1 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L6_3, - Support: 'gpu_1k', - Splitting: 'o2_FullHalf'} - ) - assert assertion_3D(scal, velo, advec) - - -def test_3D_l6_2k_sFH(): - """ - Testing Lambda6star remeshing formula in 3D, 2 kernel, - o2_FullHalf splitting. - """ - scal, velo = setup_3D() - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L6_3, - Support: 'gpu_2k', - Splitting: 'o2_FullHalf'} - ) - assert assertion_3D(scal, velo, advec) - - -def test_rectangular_domain2D(): - box = Box(length=[1., 1.], origin=[0., 0.]) - scal = Field(domain=box, name='Scalar') - velo = Field(domain=box, name='Velocity', is_vector=True) - advec = Advection(velo, scal, discretization=Discretization([65, 33]), - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: 'gpu_1k', - Splitting: 'o2'} - ) - advec_py = Advection(velo, scal, discretization=Discretization([65, 33]), - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: '', - Splitting: 'o2'}, - ) - advec.discretize() - advec_py.discretize() - advec.setup() - advec_py.setup() - - scal_d = scal.discreteFields.values()[0] - velo_d = velo.discreteFields.values()[0] - scal_d.data[0][...] = npw.asrealarray(np.random.random(scal_d.data[0].shape)) - velo_d.data[0][...] = npw.zeros_like(scal_d.data[0]) - velo_d.data[1][...] = npw.zeros_like(scal_d.data[0]) - scal_init = scal_d.data[0].copy() - scal_d.toDevice() - velo_d.toDevice() - - advec.apply(Simulation(start=0., end=0.01, nb_iter=1)) - advec_py.apply(Simulation(start=0., end=0.01, nb_iter=1)) - scal_py_res = scal_d.data[0].copy() - - scal_d.toHost() - - advec.finalize() - assert np.allclose(scal_init, scal_d.data[0]) - assert np.allclose(scal_init, scal_py_res) - - -def test_rectangular_domain3D(): - box = Box(length=[1., 1., 1.], origin=[0., 0., 0.]) - scal = Field(domain=box, name='Scalar') - velo = Field(domain=box, name='Velocity', is_vector=True) - - advec = Advection(velo, scal, discretization=Discretization([65, 33, 17]), - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: 'gpu_1k', - Splitting: 'o2'} - ) - advec_py = Advection(velo, scal, discretization=Discretization([65, 33, 17]), - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: '', - Splitting: 'o2'}, - ) - advec.discretize() - advec_py.discretize() - advec.setup() - advec_py.setup() - - scal_d = scal.discreteFields.values()[0] - velo_d = velo.discreteFields.values()[0] - scal_d.data[0][...] = npw.asrealarray(np.random.random(scal_d.data[0].shape)) - velo_d.data[0][...] = npw.zeros_like(scal_d.data[0]) - velo_d.data[1][...] = npw.zeros_like(scal_d.data[0]) - velo_d.data[2][...] = npw.zeros_like(scal_d.data[0]) - scal_init = scal_d.data[0].copy() - scal_d.toDevice() - velo_d.toDevice() - - advec.apply(Simulation(start=0., end=0.01, nb_iter=1)) - advec_py.apply(Simulation(start=0., end=0.01, nb_iter=1)) - scal_py_res = scal_d.data[0].copy() - - scal_d.toHost() - - advec.finalize() - assert np.allclose(scal_init, scal_d.data[0]) - assert np.allclose(scal_init, scal_py_res) - - -def test_2D_vector(): - box = Box(length=[1., 1.], origin=[0., 0.]) - scal = Field(domain=box, name='Scalar', is_vector=True) - velo = Field(domain=box, name='Velocity', is_vector=True) - - advec = Advection(velo, scal, discretization=Discretization([129, 129]), - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: 'gpu_1k', - Splitting: 'o2'} - ) - advec_py = Advection(velo, scal, discretization=Discretization([129, 129]), - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: '', - Splitting: 'o2'}, - ) - advec.discretize() - advec_py.discretize() - advec.setup() - advec_py.setup() - - scal_d = scal.discreteFields.values()[0] - velo_d = velo.discreteFields.values()[0] - scal_d.data[0][...] = npw.asrealarray(np.random.random(scal_d.data[0].shape)) - scal_d.data[1][...] = npw.asrealarray(np.random.random(scal_d.data[0].shape)) - velo_d.data[0][...] = npw.zeros_like(scal_d.data[0]) - velo_d.data[1][...] = npw.zeros_like(scal_d.data[0]) - scal_init_X = scal_d.data[0].copy() - scal_init_Y = scal_d.data[1].copy() - scal_d.toDevice() - velo_d.toDevice() - - advec.apply(Simulation(start=0., end=0.01, nb_iter=1)) - advec_py.apply(Simulation(start=0., end=0.01, nb_iter=1)) - scal_py_res_X = scal_d.data[0].copy() - scal_py_res_Y = scal_d.data[1].copy() - - scal_d.toHost() - - advec.finalize() - assert np.allclose(scal_init_X, scal_d.data[0]) - assert np.allclose(scal_init_Y, scal_d.data[1]) - assert np.allclose(scal_init_X, scal_py_res_X) - assert np.allclose(scal_init_Y, scal_py_res_Y) - - -def test_3D_vector(): - box = Box(length=[1., 1., 1.], origin=[0., 0., 0.]) - scal = Field(domain=box, name='Scalar', is_vector=True) - velo = Field(domain=box, name='Velocity', is_vector=True) - - advec = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: 'gpu_1k', - Splitting: 'o2'} - ) - advec_py = Advection(velo, scal,discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Support: '', - Splitting: 'o2'}, - ) - advec.discretize() - advec_py.discretize() - advec.setup() - advec_py.setup() - - scal_d = scal.discreteFields.values()[0] - velo_d = velo.discreteFields.values()[0] - scal_d.data[0][...] = npw.asrealarray(np.random.random(scal_d.data[0].shape)) - scal_d.data[1][...] = npw.asrealarray(np.random.random(scal_d.data[0].shape)) - scal_d.data[2][...] = npw.asrealarray(np.random.random(scal_d.data[0].shape)) - velo_d.data[0][...] = npw.zeros_like(scal_d.data[0]) - velo_d.data[1][...] = npw.zeros_like(scal_d.data[0]) - velo_d.data[2][...] = npw.zeros_like(scal_d.data[0]) - scal_init_X = scal_d.data[0].copy() - scal_init_Y = scal_d.data[1].copy() - scal_init_Z = scal_d.data[2].copy() - scal_d.toDevice() - velo_d.toDevice() - - advec.apply(Simulation(start=0., end=0.01, nb_iter=1)) - advec_py.apply(Simulation(start=0., end=0.01, nb_iter=1)) - scal_py_res_X = scal_d.data[0].copy() - scal_py_res_Y = scal_d.data[1].copy() - scal_py_res_Z = scal_d.data[2].copy() - - scal_d.toHost() - - advec.finalize() - assert np.allclose(scal_init_X, scal_d.data[0]) - assert np.allclose(scal_init_Y, scal_d.data[1]) - assert np.allclose(scal_init_Z, scal_d.data[2]) - assert np.allclose(scal_init_X, scal_py_res_X) - assert np.allclose(scal_init_Y, scal_py_res_Y) - assert np.allclose(scal_init_Z, scal_py_res_Z) diff --git a/hysop/gpu/tests/test_gpu_advection_null_velocity.py b/hysop/gpu/tests/test_gpu_advection_null_velocity.py new file mode 100644 index 0000000000000000000000000000000000000000..ce9d2cacf2bd8ab4e5a7db8cde4922e07b37f5eb --- /dev/null +++ b/hysop/gpu/tests/test_gpu_advection_null_velocity.py @@ -0,0 +1,436 @@ +"""Testing advection kernels with null velocity. +""" +from hysop.domain.box import Box +from hysop.fields.continuous import Field +from hysop.operator.advection import Advection +from hysop.constants import np, HYSOP_REAL +from hysop.problem.simulation import Simulation +from hysop.methods_keys import TimeIntegrator, Interpolation, Remesh, \ + Support, Splitting, Precision +from hysop.numerics.odesolvers import RK2 +from hysop.numerics.interpolation import Linear +from hysop.numerics.remeshing import L2_1, L4_2, L6_3, M8Prime +from hysop.tools.parameters import Discretization +import hysop.tools.numpywrappers as npw +from hysop.problem.simulation import O2, O2FULLHALF + + +d2d = Discretization([33, 33]) +d3d = Discretization([17, 17, 17]) +m1k = {TimeIntegrator: RK2, Interpolation: Linear, + Remesh: L4_2, Support: 'gpu_1k', Splitting: O2, + Precision: HYSOP_REAL} +m2k = {TimeIntegrator: RK2, Interpolation: Linear, + Remesh: L4_2, Support: 'gpu_2k', Splitting: O2, + Precision: HYSOP_REAL} + + +def run_advection(discr, vector_field, method=None): + """Create advection operator, ref operator + and fields, run scales and python advection, + compare results. + + Parameters + ---------- + dicr : :class:`~hysop.tools.parameters.Discretization` + chosen discretization for operator --> set domain dimension + vector_field: bool + True to advect a vector field, else scalar field. + method : dictionnary + Set scales remeshing type. + """ + dimension = len(discr.resolution) + box = Box(length=[1., ] * dimension, origin=[0., ] * dimension) + scal = Field(domain=box, name='Scalar', is_vector=vector_field) + scal_ref = Field(domain=box, name='Scalar_ref', is_vector=vector_field) + velo = Field(domain=box, name='Velocity', + formula=lambda x, y, z, t: (0., 0., 0.), is_vector=True, + vectorize_formula=True) + # build gpu advection + if method is None: + method = m1k + + advec = Advection(velocity=velo, + advected_fields=scal, + discretization=discr, + method=method) + advec.discretize() + advec.setup() + + # Get and randomize discrete fields + topo = advec.advected_fields_topology() + scal_d = scal.randomize(topo) + assert (velo.norm(topo) == 0).all() + ic = topo.mesh.compute_index + # copy data for reference + scal_ref.copy(scal, topo) + scal_ref_d = scal_ref.discretize(topo) + topo_velo = advec.velocity_topology() + velo_d = velo.discretize(topo_velo) + # transfer data to gpu + scal_d.toDevice() + velo_d.toDevice() + scal_d.wait() + velo_d.wait() + + advec.apply(Simulation(start=0., end=0.1, nb_iter=1)) + # Get data back + scal_d.toHost() + scal_d.wait() + advec.finalize() + for d in xrange(len(scal_d.data)): + assert np.allclose(scal_ref_d.data[d][ic], scal_d.data[d][ic]) + + +def test_2D_m6_1k(): + """Test M6 remeshing formula in 2D, 1 kernel, + o2 splitting. + """ + run_advection(d2d, False, m1k) + + +# def test_2D_m6_2k(): +# """Test M6 remeshing formula in 2D, 2 kernels, +# o2 splitting. +# """ +# run_advection(d2d, False, m2k) + + +# def test_2D_m6_1k_sFH(): +# """Test M6 remeshing formula in 2D, 1 kernel, +# o2 full-half splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L4_2, Support: 'gpu_1k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d2d, False, meth) + + +# def test_2D_m6_2k_sFH(): +# """Test M6 remeshing formula in 2D, 2 kernels, +# o2 full-half splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L4_2, Support: 'gpu_2k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d2d, False, meth) + + +# def test_3D_m6_1k(): +# """Test M6 remeshing formula in 3D, 1 kernel, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L4_2, Support: 'gpu_1k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_3D_m6_2k(): +# """Test M6 remeshing formula in 3D, 2 kernels, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L4_2, Support: 'gpu_2k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_3D_m6_1k_sFH(): +# """Test M6 remeshing formula in 3D, 1 kernel, +# o2 fh splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L4_2, Support: 'gpu_1k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_3D_m6_2k_sFH(): +# """Test M6 remeshing formula in 3D, 2 kernels, +# o2 fh splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L4_2, Support: 'gpu_2k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_2D_m4_1k(): +# """Test M4 remeshing formula in 2D, 1 kernel, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L2_1, Support: 'gpu_1k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(d2d, False, meth) + + +# def test_2D_m4_2k(): +# """Test M4 remeshing formula in 2D, 2 kernels, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L2_1, Support: 'gpu_2k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(d2d, False, meth) + + +# def test_2D_m4_1k_sFH(): +# """Test M4 remeshing formula in 2D, 1 kernel, +# o2 full-half splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L2_1, Support: 'gpu_1k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d2d, False, meth) + + +# def test_2D_m4_2k_sFH(): +# """Test M4 remeshing formula in 2D, 2 kernels, +# o2 full-half splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L2_1, Support: 'gpu_2k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d2d, False, meth) + + +# def test_3D_m4_1k(): +# """Test M4 remeshing formula in 3D, 1 kernel, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L2_1, Support: 'gpu_1k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_3D_m4_2k(): +# """Test M4 remeshing formula in 3D, 2 kernels, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L2_1, Support: 'gpu_2k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_3D_m4_1k_sFH(): +# """Test M4 remeshing formula in 3D, 1 kernel, +# o2 fh splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L2_1, Support: 'gpu_1k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_3D_m4_2k_sFH(): +# """Test M4 remeshing formula in 3D, 2 kernels, +# o2 fh splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L2_1, Support: 'gpu_2k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_2D_m8_1k(): +# """Test M8 remeshing formula in 2D, 1 kernel, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: M8Prime, Support: 'gpu_1k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(d2d, False, meth) + + +# def test_2D_m8_2k(): +# """Test M8 remeshing formula in 2D, 2 kernels, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: M8Prime, Support: 'gpu_2k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(d2d, False, meth) + + +# def test_2D_m8_1k_sFH(): +# """Test M8 remeshing formula in 2D, 1 kernel, +# o2 full-half splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: M8Prime, Support: 'gpu_1k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d2d, False, meth) + + +# def test_2D_m8_2k_sFH(): +# """Test M8 remeshing formula in 2D, 2 kernels, +# o2 full-half splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: M8Prime, Support: 'gpu_2k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d2d, False, meth) + + +# def test_3D_m8_1k(): +# """Test M8 remeshing formula in 3D, 1 kernel, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: M8Prime, Support: 'gpu_1k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_3D_m8_2k(): +# """Test M8 remeshing formula in 3D, 2 kernels, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: M8Prime, Support: 'gpu_2k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_3D_m8_1k_sFH(): +# """Test M8 remeshing formula in 3D, 1 kernel, +# o2 fh splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: M8Prime, Support: 'gpu_1k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_3D_m8_2k_sFH(): +# """Test M8 remeshing formula in 3D, 2 kernels, +# o2 fh splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: M8Prime, Support: 'gpu_2k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_2D_l6_1k(): +# """Test L6 remeshing formula in 2D, 1 kernel, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L6_3, Support: 'gpu_1k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(d2d, False, meth) + + +# def test_2D_l6_2k(): +# """Test L6 remeshing formula in 2D, 2 kernels, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L6_3, Support: 'gpu_2k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(d2d, False, meth) + + +# def test_2D_l6_1k_sFH(): +# """Test L6 remeshing formula in 2D, 1 kernel, +# o2 full-half splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L6_3, Support: 'gpu_1k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d2d, False, meth) + + +# def test_2D_l6_2k_sFH(): +# """Test L6 remeshing formula in 2D, 2 kernels, +# o2 full-half splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L6_3, Support: 'gpu_2k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d2d, False, meth) + + +# def test_3D_l6_1k(): +# """Test L6 remeshing formula in 3D, 1 kernel, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L6_3, Support: 'gpu_1k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_3D_l6_2k(): +# """Test L6 remeshing formula in 3D, 2 kernels, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L6_3, Support: 'gpu_2k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_3D_l6_1k_sFH(): +# """Test L6 remeshing formula in 3D, 1 kernel, +# o2 fh splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L6_3, Support: 'gpu_1k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_3D_l6_2k_sFH(): +# """Test L6 remeshing formula in 3D, 2 kernels, +# o2 fh splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L6_3, Support: 'gpu_2k', Splitting: O2FULLHALF, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) + + +# def test_rectangular_domain2D(): +# """Test remeshing formula in 2D, with different resolutions in each dir, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L4_2, Support: 'gpu_1k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(Discretization([65, 33]), False, meth) + + +# def test_rectangular_domain3D(): +# """Test remeshing formula in 3D, with different resolutions in each dir, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L4_2, Support: 'gpu_1k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(Discretization([65, 33, 17]), False, meth) + + +# def test_vector_2D(): +# """Test remeshing formula in 2D, advection of a vector field, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L4_2, Support: 'gpu_1k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(Discretization([129, 129]), True, meth) + + +# def test_vector_3D(): +# """Test remeshing formula in 3D, advection of a vector field, +# o2 splitting. +# """ +# meth = {TimeIntegrator: RK2, Interpolation: Linear, +# Remesh: L4_2, Support: 'gpu_1k', Splitting: O2, +# Precision: HYSOP_REAL} +# run_advection(d3d, False, meth) diff --git a/hysop/gpu/tests/test_advection_randomVelocity.py b/hysop/gpu/tests/test_gpu_advection_random_velocity.py similarity index 100% rename from hysop/gpu/tests/test_advection_randomVelocity.py rename to hysop/gpu/tests/test_gpu_advection_random_velocity.py diff --git a/hysop/gpu/tests/test_opencl_environment.py b/hysop/gpu/tests/test_opencl_environment.py index 0d7d0d6d720b779e3507b85ede1c20d0a2ebfe51..d4da68f0c0b7e26d0f713e6ef40a4b676936d149 100644 --- a/hysop/gpu/tests/test_opencl_environment.py +++ b/hysop/gpu/tests/test_opencl_environment.py @@ -1,6 +1,46 @@ -from hysop.constants import np -from hysop.gpu.tools import get_opencl_environment +"""Test hysop implementation of OpenCL basic functionnalities""" +import numpy as np +from hysop.gpu.tools import get_opencl_environment, explore, OpenCLEnvironment +from hysop.constants import HYSOP_REAL FLOAT_GPU = np.float32 +from hysop.mpi import main_comm +import pyopencl as cl + + +def test_opencl_env_default(): + """Test gpu tools. Just call get_opencl_environment and explore functions. + """ + # Display devices info. + explore() + # Create default opencl env. + cl_env = get_opencl_environment() + assert isinstance(cl_env, OpenCLEnvironment) + assert cl_env.device is not None + assert cl_env.ctx is not None + assert cl_env.queue is not None + assert cl_env.precision == HYSOP_REAL + + +def test_opencl_env(): + """Test gpu tools. Just call get_opencl_environment and explore functions. + """ + # Display devices info. + explore() + # Create default opencl env. + comm = main_comm.Dup() + plt = cl.get_platforms()[-1] + device = plt.get_devices()[-1] + nb_platforms = len(cl.get_platforms()) + nb_devices = len(plt.get_devices()) + cl_env = get_opencl_environment(platform_id=nb_platforms - 1, + device_id=nb_devices - 1, + precision=FLOAT_GPU, comm=comm) + assert isinstance(cl_env, OpenCLEnvironment) + assert cl_env.platform == plt + assert cl_env.device == device + assert cl_env.ctx is not None + assert cl_env.queue is not None + assert cl_env.precision == FLOAT_GPU def test_queue_unique_creation(): diff --git a/hysop/gpu/tools.py b/hysop/gpu/tools.py index d0ac1424741aadc95066e485ed9edec0eb364a7f..ff13a59efc51c07016b66cda0d2858b491b18828 100644 --- a/hysop/gpu/tools.py +++ b/hysop/gpu/tools.py @@ -1,16 +1,30 @@ -"""Misc. tools for gpu/OpenCL management. +"""Classes and tools used to handle openCL interface. + + +* :class:`~hysop.gpu.tools.OpenCLEnvironment`: + object handling opencl platform, device ... info. +* :func:`~hysop.gpu.tools.get_opengl_shared_environment`: + build or get an OpenCL environment with openGL properties. +* :func:`~hysop.gpu.tools.get_opencl_environment`: + build or get an OpenCL environment. +* :func:`~hysop.gpu.tools.explore` + explore system and display platform, devices, memory ... info. + + """ +from __future__ import print_function from hysop import __VERBOSE__, __DEFAULT_PLATFORM_ID__, __DEFAULT_DEVICE_ID__ from hysop.constants import np, HYSOP_REAL -from hysop.gpu import cl, GPU_SRC, CL_PROFILE +from hysop.gpu import cl, CL_PROFILE, GPU_SRC import re +import os from hysop.mpi import MPI FLOAT_GPU, DOUBLE_GPU = np.float32, np.float64 -# Global variable handling the OpenCL Environment instance __cl_env = None +"""Global variable handling the OpenCL Environment instance """ class OpenCLEnvironment(object): @@ -37,7 +51,7 @@ class OpenCLEnvironment(object): Default=False. comm : mpi communicator, optional Communicator which handles the OpenCL env. - Default = hysop.mpi.main_comm + Default = :data:`~hysop.mpi.main_comm` """ self._platform_id = platform_id @@ -81,7 +95,8 @@ class OpenCLEnvironment(object): self.default_build_opts = "" if CL_PROFILE and self.device.vendor.find('NVIDIA') >= 0: self.default_build_opts += " -cl-nv-verbose" - self.default_build_opts += " -Werror" + self._get_precision_opts() + #self.default_build_opts += "-Werror" + self._get_precision_opts() + self.default_build_opts += self._get_precision_opts() # Kernels configuration dictionary if self.device.name == "Cayman": @@ -90,7 +105,7 @@ class OpenCLEnvironment(object): self.device.name == "Tesla K20Xm": from hysop.gpu.config_k20m import kernels_config as kernel_cfg else: - print "/!\\ Get a defautl kernels config for", self.device.name + print("/!\\ Get a defautl kernels config for", self.device.name) from hysop.gpu.config_default import kernels_config as kernel_cfg self.kernels_config = kernel_cfg self._locMem_Buffers = {} @@ -116,13 +131,13 @@ class OpenCLEnvironment(object): """ platform_changed, device_changed = False, False if not platform_id == self._platform_id: - print ("platform changed") + print("platform changed") self._platform_id = platform_id self.platform = self._get_platform(platform_id) platform_changed = True - if platform_changed or not (device_id is self._device_id - and device_type == self._device_type): - print ("device changed") + if platform_changed or not (device_id is self._device_id and + device_type == self._device_type): + print("device changed") self._device_id = device_id self._device_type = device_type self.device = self._get_device(self.platform, @@ -132,21 +147,23 @@ class OpenCLEnvironment(object): if platform_changed or device_changed or \ (not self._gl_sharing and gl_sharing is not self._gl_sharing): if self._gl_sharing and not gl_sharing: - print ("Warning: Loosing Gl shared context.") + print("Warning: Loosing Gl shared context.") self._gl_sharing = gl_sharing self.ctx = self._get_context(self.device, gl_sharing) self.queue = self._get_queue(self.ctx) if self.precision is not precision and precision is not None: if self.precision is not None: - print ("Warning, GPU precision is overrided from",) - print (self.precision, 'to', precision) + print("Warning, GPU precision is overrided from",) + print(self.precision, 'to', precision) self.precision = precision self.default_build_opts = "" if CL_PROFILE and self.device.vendor.find('NVIDIA') >= 0: self.default_build_opts += " -cl-nv-verbose" - self.default_build_opts += "-Werror" + self._get_precision_opts() + #self.default_build_opts += "-Werror" + self._get_precision_opts() + self.default_build_opts += self._get_precision_opts() - def _get_platform(self, platform_id): + @staticmethod + def _get_platform(platform_id): """Returns an OpenCL platform :param platform_id : OpenCL platform id @@ -155,17 +172,19 @@ class OpenCLEnvironment(object): # OpenCL platform platform = cl.get_platforms()[platform_id] except IndexError: - print (" Incorrect platform_id :", platform_id, ".",) - print (" Only ", len(cl.get_platforms()), " available.",) - print (" Getting default platform. ") - platform = cl.get_platforms()[0] + plist = cl.get_platforms() + platform = plist[0] + print(" Incorrect platform_id :", platform_id, ".",) + print(" Only ", len(plist), " available.",) + print(" --> getting default platform ", platform.name) if __VERBOSE__: - print (" Platform ") - print (" - Name :", platform.name) - print (" - Version :", platform.version) + print(" Platform ") + print(" - Name :", platform.name) + print(" - Version :", platform.version) return platform - def _get_device(self, platform, device_id, device_type): + @staticmethod + def _get_device(platform, device_id, device_type): """Returns an OpenCL device Parameters @@ -183,42 +202,42 @@ class OpenCLEnvironment(object): display = False try: if device_type is not None: - device = platform.get_devices( - eval("cl.device_type." + str(device_type.upper())) - )[device_id] + device_type_id = cl.device_type.__getattribute__( + cl.device_type, str(device_type.upper())) + device = platform.get_devices(device_type_id)[device_id] else: device = platform.get_devices()[device_id] except cl.RuntimeError as e: - print ("RuntimeError:", e) + print("RuntimeError:", e) device = cl.create_some_context().devices[0] display = True except AttributeError as e: - print ("AttributeError:", e) + print("AttributeError:", e) device = cl.create_some_context().devices[0] display = True except IndexError: - print (" Incorrect device_id :", device_id, ".",) - print (" Only ", len(platform.get_devices()), " available.",) + print(" Incorrect device_id :", device_id, ".",) + print(" Only ", len(platform.get_devices()), " available.",) if device_type is not None: - print (" Getting first device of type " + - str(device_type.upper())) + device_type = str(device_type.upper()) + print(" Getting first device of type " + device_type) else: - print (" Getting first device of the platform") + print(" Getting first device of the platform") device = platform.get_devices()[0] display = True if device_type is not None: assert device_type.upper() == cl.device_type.to_string(device.type) if display or __VERBOSE__: - print (" Device") - print (" - id :", device_id) - print (" - Name :",) - print (device.name) - print (" - Type :",) - print cl.device_type.to_string(device.type) + print(" Device") + print(" - id :", device_id) + print(" - Name :",) + print(device.name) + print(" - Type :",) + print(cl.device_type.to_string(device.type)) print (" - C Version :",) print (device.opencl_c_version) print (" - Global mem size :",) - print device.global_mem_size / (1024 ** 3), "GB" + print(device.global_mem_size / (1024 ** 3), "GB") return device def _get_context(self, device, gl_sharing): @@ -249,12 +268,13 @@ class OpenCLEnvironment(object): else: ctx = cl.Context([device]) if __VERBOSE__: - print " Context:" + print(" Context:") if props is not None: - print " - properties :", props + print(" - properties :", props) return ctx - def _get_queue(self, ctx): + @staticmethod + def _get_queue(ctx): """Returns OpenCL queue from context :param ctx : OpenCL context @@ -267,10 +287,10 @@ class OpenCLEnvironment(object): else: queue = cl.CommandQueue(ctx) if __VERBOSE__: - print " Queue" + print(" Queue") if props is not None: - print " - properties :", props - print "===" + print(" - properties :", props) + print("===") return queue def create_other_queue(self): @@ -278,7 +298,8 @@ class OpenCLEnvironment(object): """ return self._get_queue(self.ctx) - def get_work_items(self, resolution, vector_width=1): + @staticmethod + def get_work_items(resolution, vector_width=1): """Set the optimal work-item number and OpenCL space index. Parameters @@ -309,12 +330,12 @@ class OpenCLEnvironment(object): # Change work-item regarding problem size if resolution[0] % workItemNumber > 0: if len(resolution) == 3: - print "Warning : GPU best performances obtained for", - print "problem sizes multiples of 64" + print("Warning : GPU best performances obtained for",) + print("problem sizes multiples of 64") else: - print "Warning : GPU best performances obtained for", - print "problem sizes multiples of 256" - while(resolution[0] % workItemNumber > 0): + print("Warning : GPU best performances obtained for",) + print("problem sizes multiples of 256") + while resolution[0] % workItemNumber > 0: workItemNumber = workItemNumber / 2 # Change work-item regarding vector_width if workItemNumber * vector_width > resolution[0]: @@ -334,13 +355,11 @@ class OpenCLEnvironment(object): return workItemNumber, gwi, lwi def _get_precision_opts(self): - """Check if device is capable regarding given precision - and returns build options regarding the required precision. + """Check if device is capable to work with given precision + and returns build options considering this precision """ opts = "" # Precision supported - if __VERBOSE__: - print " Precision capability ", fp32_rounding_flag = True if self.precision is FLOAT_GPU: opts += " -cl-single-precision-constant" @@ -350,7 +369,8 @@ class OpenCLEnvironment(object): raise ValueError("Double Precision is not supported by device") prec = "double" if __VERBOSE__: - print "for " + prec + " Precision: " + print(" Precision capability ",) + print("for " + prec + " Precision: ") for v in ['DENORM', 'INF_NAN', 'ROUND_TO_NEAREST', 'ROUND_TO_ZERO', 'ROUND_TO_INF', 'FMA', 'CORRECTLY_ROUNDED_DIVIDE_SQRT', 'SOFT_FLOAT']: @@ -359,7 +379,7 @@ class OpenCLEnvironment(object): ' cl.device_fp_config.' + v + ') == cl.device_fp_config.' + v): if __VERBOSE__: - print v + print(v) else: if v is 'CORRECTLY_ROUNDED_DIVIDE_SQRT': fp32_rounding_flag = False @@ -367,28 +387,26 @@ class OpenCLEnvironment(object): if v is 'CORRECTLY_ROUNDED_DIVIDE_SQRT': fp32_rounding_flag = False if __VERBOSE__: - print v, 'is not supported in OpenCL C 1.2.\n', - print ' Exception catched : ', ae + print(v, 'is not supported in OpenCL C 1.2.\n',) + print(' Exception catched : ', ae) if fp32_rounding_flag: opts += " -cl-fp32-correctly-rounded-divide-sqrt" return opts - def build_src(self, files, options="", vector_width=4, - nb_remesh_components=1): + def _create_cl_program(self, file_list, vector_width=4, + nb_remesh_components=1): """Build OpenCL sources Parameters ---------- - files : file or list of files - sources - options : string - Compiler options to use for buildind - vector_width : int - OpenCL vector type width - nb_remesh_components : int - number of remeshed components + file_list : list of strings + user defined files names + vector_width : int, optional + OpenCL vector type width, default=4 + nb_remesh_components : int, optional + number of remeshed components, default=1 - Returns OpenCL binaries + Returns OpenCL kernel Parse the sources to handle single and double precision. """ @@ -396,28 +414,22 @@ class OpenCLEnvironment(object): if cl.device_type.to_string(self.device.type) == 'GPU' and \ self.precision is DOUBLE_GPU: gpu_src += "#pragma OPENCL EXTENSION cl_khr_fp64: enable \n" - if isinstance(files, list): - file_list = files - else: - file_list = [files] - if __VERBOSE__: - print "=== Kernel sources compiling ===" - for sf in file_list: - print " - ", sf for sf in file_list: + # search and open cl file. try: f = open(sf, 'r') except IOError as ioe: if ioe.errno == 2: + # path to cl files inside hysop.gpu package f = open(GPU_SRC + sf, 'r') else: raise ioe gpu_src += "".join( self.parse_file(f, vector_width, nb_remesh_components)) f.close() - #print gpu_src - for k in self.macros: - gpu_src = gpu_src.replace(k, str(self.macros[k])) + if self.macros is not None: + for k in self.macros: + gpu_src = gpu_src.replace(k, str(self.macros[k])) if self.precision is FLOAT_GPU: # Rexexp to add 'f' suffix to float constants # Match 1.2, 1.234, 1.2e3, 1.2E-05 @@ -427,41 +439,81 @@ class OpenCLEnvironment(object): float_replace.sub(r'\g<float>f', gpu_src)) else: prg = cl.Program(self.ctx, gpu_src.replace('float', 'double')) - # OpenCL program + return prg + + def build_src(self, files, options=None, vector_width=4, + nb_remesh_components=1): + """Build OpenCL sources + + Parameters + ---------- + files : string or list of strings + user defined files names + options : string, optional + Compiler options, default=None + vector_width : int, optional + OpenCL vector type width, default=4 + nb_remesh_components : int, optional + number of remeshed components, default=1 + + Returns OpenCL binaries + + Parse the sources to handle single and double precision. + """ + if isinstance(files, list): + file_list = files + else: + file_list = [files] + if __VERBOSE__: + print("=== Kernel sources compiling ===") + for sf in file_list: + print(" - ", sf) + + # --- create kernel from cl files --- + prg = self._create_cl_program(file_list, vector_width, + nb_remesh_components) + # --- build kernel --- + if options is None: + options = "" + # --- Build kernel --- try: build = prg.build(self.default_build_opts + options) except Exception, e: - print "Build files : " + print("Build files : ") for sf in file_list: - print " - ", sf - print "Build options : ", self.default_build_opts + options - print "Vectorization : ", vector_width + print(" - ", sf) + print("Build options : ", self.default_build_opts + options) + print("Vectorization : ", vector_width) raise e + + # display post-build info if __VERBOSE__: #print options - print "Build options : ", - print build.get_build_info( - self.device, cl.program_build_info.OPTIONS) - print "Compiler status : ", - print build.get_build_info( - self.device, cl.program_build_info.STATUS) - print "Compiler log : ", - print build.get_build_info(self.device, - cl.program_build_info.LOG) - print "===\n" + print("Build options : ",) + print(build.get_build_info( + self.device, cl.program_build_info.OPTIONS)) + print("Compiler status : ",) + print(build.get_build_info( + self.device, cl.program_build_info.STATUS)) + print("Compiler log : ",) + print(build.get_build_info(self.device, + cl.program_build_info.LOG)) + print("===\n") elif CL_PROFILE: - print "Build files: " + str(file_list) - print "With build options: " + self.default_build_opts + options - print "Compiler output : " + build.get_build_info( - self.device, cl.program_build_info.LOG) + print("Build files: " + str(file_list)) + print("With build options: " + self.default_build_opts + options) + print("Compiler output : " + build.get_build_info( + self.device, cl.program_build_info.LOG)) return build - def parse_file(self, f, n=8, nb_remesh_components=1): + @staticmethod + def parse_file(f, n=8, nb_remesh_components=1): """Parse a file containing OpenCL sources. Parameters ---------- - f : file + f : string + file name n : int, optional vector width, default=8 nb_remesh_components : int @@ -471,49 +523,51 @@ class OpenCLEnvironment(object): ------- string, the parsed sources. - - <code>__N__</code> is expanded as an integer corresponding to a - vector with - - <code>__NN__</code>, instruction is duplicated to operate on each - vector component: - - if line ends with '<code>;</code>', the whole instruciton is - duplicated. - - if line ends with '<code>,</code>' and contains - '<code>(float__N__)(</code>', the float element is duplicated - - Remeshing fields components are expanded as follows : - All code between '<code>__RCOMPONENT_S__</code>' and - '<code>__RCOMPONENT_E__</code>' flags are duplicated n times with n + Notes + ----- + * __N__ is expanded as an integer corresponding to vector width. + * __NN__ instruction is duplicated to operate on each vector component: + + * if line ends with ';', the whole instruciton is + duplicated. + * if line ends with ',' and contains + '(float__N__)(', the float element is duplicated + + * Remeshing fields components are expanded as follows : + All code between '__RCOMPONENT_S__' and + '__RCOMPONENT_E__' flags are duplicated n times with n the number of components to compute. In this duplicated code, the - flag '<code>__ID__</code>' is replaced by index of a range of lenght - the number of components. A flag '<code>__RCOMPONENT_S__P__</code>' - may be used and the duplicated elements are separated by ',' (for - function parameters expanding). - - Examples with a 4-width vector:\n - \code - float__N__ x; -> float4 x; - - x.s__NN__ = 1.0f; -> x.s0 = 1.0f; - x.s1 = 1.0f; - x.s2 = 1.0f; - x.s3 = 1.0f; - - x = (int__N__)(\__NN__, -> x = (int4)(0, - ); 1, - 2, - 3, - ); - \endcode - - Examples with a 2 components expansion:\n - __RCOMP_P __global const float* var__ID__, - -> __global const float* var0,__global const float* var1, - - __RCOMP_I var__ID__[i] = 0.0; - -> var0[i] = 0.0;var1[i] = 0.0; - - aFunction(__RCOMP_P var__ID__, __RCOMP_P other__ID__); - -> aFunction(var0, var1, other0, other1); - \endcode + flag '__ID__' is replaced by index of a range of lenght + the number of components. A flag '__RCOMPONENT_S__P__' + may be used and the duplicated elements are separated with ',' + (for function parameters expanding). + + Examples with a 4-width vector code:: + + float__N__ x; -> float4 x; + + x.s__NN__ = 1.0f; -> x.s0 = 1.0f; + x.s1 = 1.0f; + x.s2 = 1.0f; + x.s3 = 1.0f; + + x = (int__N__)(__NN__, -> x = (int4)(0, + ); 1, + 2, + 3, + ); + + Examples with a 2 components expansion code:: + + __RCOMP_P __global const float* var__ID__, + -> __global const float* var0,__global const float* var1, + + __RCOMP_I var__ID__[i] = 0.0; + -> var0[i] = 0.0;var1[i] = 0.0; + + aFunction(__RCOMP_P var__ID__, __RCOMP_P other__ID__); + -> aFunction(var0, var1, other0, other1); + """ src = "" # replacement for floatN elements @@ -560,14 +614,23 @@ class OpenCLEnvironment(object): return src def global_allocation(self, array): + """Allocate and returns an opencl buffer + + Parameters + ---------- + array : numpy array + source buffer, on host + """ + # create an opencl buffer from input array clBuff = cl.Buffer(self.ctx, cl.mem_flags.ALLOC_HOST_PTR, size=array.nbytes) # Touch the buffer on device to performs the allocation - # Transfers a single element in device (the precision no matters here) + # Transfers a single element in device (the precision does not matter) e = np.zeros((1,), dtype=np.float64) cl.enqueue_copy(self.queue, clBuff, e, buffer_origin=(0, 0, 0), host_origin=(0, 0, 0), region=(e.nbytes,)).wait() + # update memory counter self.available_mem -= clBuff.size return clBuff @@ -620,7 +683,7 @@ def get_opengl_shared_environment(platform_id=None, device_id=None, device_type=None, precision=HYSOP_REAL, comm=None): - """Build or get an OpenCL environment. + """Build or get an OpenCL environment with openGL properties. Parameters ---------- @@ -640,7 +703,7 @@ def get_opengl_shared_environment(platform_id=None, Returns ------- :class:`~hysop.gpu.tools.OpenCLEnvironment` - handling OpenCL platform, device, context and queue + object handling OpenCL platform, device, context and queue The context is obtained with gl-shared properties depending on the OS. """ @@ -681,8 +744,9 @@ def get_opencl_environment(platform_id=None, Returns ------- + :class:`~hysop.gpu.tools.OpenCLEnvironment` - handling OpenCL platform, device, context and queue + object handling OpenCL platform, device, context and queue """ if platform_id is None: @@ -700,8 +764,8 @@ def get_opencl_environment(platform_id=None, def explore(): - """Print environment details""" - print "OpenCL exploration : " + """Scan system and print OpenCL environment details""" + print("OpenCL exploration : ") platforms = cl.get_platforms() platforms_info = ["name", "version", "vendor", "profile", "extensions"] devices_info = ["name", @@ -731,11 +795,11 @@ def explore(): "preferred_vector_width_float", "preferred_vector_width_int"] for pltfm in platforms: - print "Platform:", pltfm.name + print("Platform:", pltfm.name) for pltfm_info in platforms_info: - print " |-", pltfm_info, ':', eval("pltfm." + pltfm_info) + print(" |-", pltfm_info, ':', eval("pltfm." + pltfm_info)) devices = pltfm.get_devices() for dvc in devices: - print " |- Device:", dvc.name + print(" |- Device:", dvc.name) for dvc_info in devices_info: - print " |-", dvc_info, ':', eval("dvc." + dvc_info) + print(" |-", dvc_info, ':', eval("dvc." + dvc_info)) diff --git a/hysop/methods_keys.py b/hysop/methods_keys.py index ec4601ed3d02e5d1432c25dd7f519beeee090b69..c00043a2d22ac3daa33a59cba803cd42b4512da8 100644 --- a/hysop/methods_keys.py +++ b/hysop/methods_keys.py @@ -26,8 +26,6 @@ Formulation = 44444 ## Space discretisation method ## (see for example hysop.numerics.finite_differences) SpaceDiscretisation = 55555 -## Scales method parameters -Scales = 66666 ## Splitting method Splitting = 77777 ## Device (either GPU or CPU) diff --git a/hysop/numerics/differential_operations.py b/hysop/numerics/differential_operations.py index c6ba5544573d8f50095bf5f44746e8035ad8a285..1f050dc8a84fb6fbb88895b609b40a0bf587db3c 100755 --- a/hysop/numerics/differential_operations.py +++ b/hysop/numerics/differential_operations.py @@ -308,7 +308,7 @@ class DivRhoV(DifferentialOperation): iout = None fd_scheme = self.method(topo.mesh.space_step, self._wk_indices, - indices_out=iout) + output_indices=iout) self._fd_work = self.method(topo.mesh.space_step, self._wk_indices) msg = 'Ghost layer is too small for the chosen FD scheme.' diff --git a/hysop/numerics/finite_differences.py b/hysop/numerics/finite_differences.py index dbf51487dd75011a5000718a51668e13639fe095..805682ca3fac11f7b13caf8a81bf74e9891d465e 100644 --- a/hysop/numerics/finite_differences.py +++ b/hysop/numerics/finite_differences.py @@ -69,6 +69,9 @@ class FiniteDifference(object): or result += the derivative of tab according to dir: >> scheme.compute_and_add(tab, dir, result) + + >>> from hysop import Box + >>> domain = Box() Notes FP : Compute method is much more performant than compute and add @@ -84,7 +87,7 @@ class FiniteDifference(object): def __new__(cls, *args, **kw): return object.__new__(cls, *args, **kw) - def __init__(self, step, indices, indices_out=None): + def __init__(self, step, indices, output_indices=None): """ Parameters @@ -95,9 +98,9 @@ class FiniteDifference(object): Represents the local mesh on which finite-differences will be applied, like compute_index in :class:`~hysop.domain.mesh.Mesh`. - indices_out : list of slices, optional + output_indices : list of slices, optional set which indices of result will be modified. By default, - indices_out = indices + output_indices = indices Attributes ---------- @@ -118,24 +121,24 @@ class FiniteDifference(object): self._coeff = None # List of slices representing the mesh on which fd scheme is applied self.indices = indices - if indices_out is None or indices_out is False: + if output_indices is None or output_indices is False: # This is the default case. It means that # result shape is unknown. It will # be the same as input tab shape during call # and input/output indices are the same: # we compute result[indices] = scheme(tab[indices]) self.output_indices = self.indices - elif isinstance(indices_out, list): + elif isinstance(output_indices, list): # In that case, result shape remains unknown but # the position where result is written is not the # same as input indices. # We compute: - # result[indices_out] = scheme(tab[indices]) - # Warning : result shape/indices_out compliance + # result[output_indices] = scheme(tab[indices]) + # Warning : result shape/output_indices compliance # is not checked! - self.output_indices = indices_out + self.output_indices = output_indices - elif indices_out is True: + elif output_indices is True: # Here, output shape is fixed to the minimal required # shape, computed with indices. # So, we compute: diff --git a/hysop/numerics/interpolation.py b/hysop/numerics/interpolation.py index 167efede0c16a21a25eb28f2258304e2ae4a19c2..d8099c7d5a8edeae58a5f4c3ba82f0fb143cfbeb 100644 --- a/hysop/numerics/interpolation.py +++ b/hysop/numerics/interpolation.py @@ -5,7 +5,7 @@ * :class:`~interpolation.Linear` """ -from hysop.constants import HYSOP_INTEGER, EPS, XDIR, YDIR, ZDIR +from hysop.constants import HYSOP_INTEGER, EPS from hysop.tools.misc import WorkSpaceTools import numpy as np import hysop.tools.numpywrappers as npw @@ -257,26 +257,37 @@ class Interpolation(object): work[0][ic][ilist] = work[1][ic][ilist] return True + def __str__(self): + """Print method + """ + disp = self.__class__.__name__ + ' interpolation scheme,' + disp += ' from grid of shape ' + str(self.topo_source.mesh.resolution) + disp += ' to grid of shape ' + str(self.topo_target.mesh.resolution) + return disp + class Linear(Interpolation): """Linear interpolation of a field""" ghosts_layer_size = 1 - def __init__(self, tab, **kwds): + def __init__(self, interpolated_field, **kwds): """ Parameters ---------- - tab : numpy array + interpolated_field : numpy array data to be interpolated kwds : base class parameters Notes ----- - * topo argument corresponds to the topology of field tab + * topo argument corresponds to the topology of field interpolated_field + * it would be probably better to pass interpolated_field + during __call__ but interpolation is used as rhs in odesolvers and + call function must fit with rhs arg list. """ - self.tab = tab + self.interpolated_field = interpolated_field super(Linear, self).__init__(**kwds) @staticmethod @@ -299,7 +310,7 @@ class Linear(Interpolation): return {'rwork': [noghosts_shape, ], 'iwork': [noghosts_shape, ]} def __call__(self, t, positions, result): - """Find 'self.tab' values at nodes positions using linear interpolation. + """interpolate values at nodes positions using linear interpolation. Parameters ---------- @@ -323,13 +334,16 @@ class Linear(Interpolation): assert isinstance(positions, list) assert isinstance(result, list) assert result[0].shape == positions[0].shape - assert self._assert_target_in_source(positions) + assert (self.interpolated_field.shape == + self.topo_source.mesh.resolution).all() + #assert self._assert_target_in_source(positions) # linear interpolation: - # input = self.tab (field values at points xi), + # input = interpolated_field (field values at points xi), # pos, a vector of p points - # output = res(pos), interpolation of tab at points pos - # res(pos_p) = (pos_p - xi)/h . tab_i+1 + (x_i+1 - pos_p)/h . tab_i - # = iy.tab_i+1 + (1 - iy).tab_i + # output = res(pos), interpolation of interpolated_field at points pos + # intf =interpolated_field + # res(pos_p) = (pos_p - xi)/h . intf_i+1 + (x_i+1 - pos_p)/h . intf_i + # = iy.intf_i+1 + (1 - iy).intf_i # where iy = (pos_p - x0)/h - i, i = floor((pos_p - x0)/h) # number of points to interpolate @@ -345,10 +359,10 @@ class Linear(Interpolation): left_points_indices = self._compute_iy(positions[0], wk) assert npw.arrays_share_data(left_points_indices[self.direction], self._iwork[0]) - wk[1][self._target_indices] = self.tab[left_points_indices] *\ - (1. - wk[0][...]) + wk[1][self._target_indices] = \ + self.interpolated_field[left_points_indices] * (1. - wk[0][...]) # compute 'right points' indices left_points_indices[self.direction] += 1 wk[1][self._target_indices] += \ - self.tab[left_points_indices] * wk[0][...] + self.interpolated_field[left_points_indices] * wk[0][...] return result diff --git a/hysop/numerics/odesolvers.py b/hysop/numerics/odesolvers.py index 1ef90454caee2dccefc7ae143851999d217505ac..d54eaebc0319ff1efec2842676cff402159dd83b 100755 --- a/hysop/numerics/odesolvers.py +++ b/hysop/numerics/odesolvers.py @@ -2,12 +2,13 @@ .. currentmodule hysop.numerics -* :class:`~odesolvers.Euler`, +* :class:`~hysop.numerics.odesolvers.Euler`, * :class:`~odesolvers.RK2`, * :class:`~odesolvers.RK3`, * :class:`~odesolvers.RK4`, * :class:`~odesolvers.ODESolver` (abstract base class). +See :ref:`odesolvers` """ @@ -294,9 +295,10 @@ class RK2(ODESolver): Note : since result may be equal to y, it can not be used as a temporary workspace. """ + ic = self.in_indices for i in xrange(self._nb_components): - cond = [y[i] is self.rwork[j] for j in xrange(len(self.rwork))] - assert cond.count(True) is 0 + for j in xrange(len(self.rwork)): + assert not npw.arrays_share_data(y[i], self.rwork[j]) assert len(result) == self._nb_components work0 = self.rwork[:self._nb_components] yn = self.rwork[self._nb_components:2 * self._nb_components] @@ -307,7 +309,7 @@ class RK2(ODESolver): # yn for i in xrange(self._nb_components): - np.add(y[i], work0[i] * 0.5 * dt, yn[i]) + np.add(y[i][ic], work0[i] * 0.5 * dt, yn[i]) # Update ghosts self._synchronize.apply(yn) @@ -318,7 +320,7 @@ class RK2(ODESolver): #np.multiply(work0[i], dt, work0[i]) work0[i] *= dt # result = y + work0 - np.add(work0[i], y[i], result[i]) + np.add(work0[i], y[i][ic], result[i]) return result diff --git a/hysop/numerics/tests/test_finite_differences.py b/hysop/numerics/tests/test_finite_differences.py index b098eb7b03b12cb43fd210fe543a9e8795af8a60..584df6a81fd2c8bb576b5fe80f27e5c03bedf954 100644 --- a/hysop/numerics/tests/test_finite_differences.py +++ b/hysop/numerics/tests/test_finite_differences.py @@ -116,9 +116,9 @@ def check_fd_schemes_reduce_output(discr, formulas): ind = topo.mesh.compute_index hh = topo.mesh.space_step sc = { - FDC2(hh, ind, indices_out=True): [refd, 0, 2], - FDC4(hh, ind, indices_out=True): [refd, 0, 4], - FD2C2(hh, ind, indices_out=True): [ref2d, 1, 2] + FDC2(hh, ind, output_indices=True): [refd, 0, 2], + FDC4(hh, ind, output_indices=True): [refd, 0, 4], + FD2C2(hh, ind, output_indices=True): [ref2d, 1, 2] } shape = np.asarray(topo.mesh.resolution).copy() shape -= 2 * g @@ -136,9 +136,9 @@ def check_fd_schemes_reduce_all(discr, formulas): subbox = SubBox(parent=topo.domain, origin=orig, length=sl) ind = subbox.discretize(topo)[0] sc = { - FDC2(hh, ind, indices_out=True): [refd, 0, 2], - FDC4(hh, ind, indices_out=True): [refd, 0, 4], - FD2C2(hh, ind, indices_out=True): [ref2d, 1, 2] + FDC2(hh, ind, output_indices=True): [refd, 0, 2], + FDC4(hh, ind, output_indices=True): [refd, 0, 4], + FD2C2(hh, ind, output_indices=True): [ref2d, 1, 2] } shape = subbox.mesh[topo].resolution result = npw.zeros(shape) @@ -159,9 +159,9 @@ def check_fd_schemes_setiout(discr, formulas): result = npw.zeros(shape) iout = [slice(3, shape[i] - 3) for i in xrange(len(shape))] sc = { - FDC2(hh, ind, indices_out=iout): [refd, 0, 2], - FDC4(hh, ind, indices_out=iout): [refd, 0, 4], - FD2C2(hh, ind, indices_out=iout): [ref2d, 1, 2] + FDC2(hh, ind, output_indices=iout): [refd, 0, 2], + FDC4(hh, ind, output_indices=iout): [refd, 0, 4], + FD2C2(hh, ind, output_indices=iout): [ref2d, 1, 2] } run_schemes(sc, f1d, result, ind, hh, iout) diff --git a/hysop/numerics/tests/test_interpolation.py b/hysop/numerics/tests/test_interpolation.py new file mode 100644 index 0000000000000000000000000000000000000000..51a9b4d3a16d0021408c7337e54ffa972b139de2 --- /dev/null +++ b/hysop/numerics/tests/test_interpolation.py @@ -0,0 +1,123 @@ +"""Test interpolation schemes +""" +from hysop.numerics.interpolation import Linear +import hysop.tools.numpywrappers as npw +from hysop import Box, Discretization +import numpy as np +cos = np.cos + +disc = [None, ] * 3 +g = 1 +nx = 33 +ny = 69 +nz = 129 +disc[0] = Discretization([nx, ], [g, ]) +disc[1] = Discretization([nx, ny], [g, g]) +disc[2] = Discretization([nx, ny, nz], [g, g, g]) +doms = [None, ] * 3 +topos = [None, ] * 3 +rworks = [None, ] * 3 +iworks = [None, ] * 3 +for case in xrange(3): + dim = case + 1 + doms[case] = Box(origin=[0.] * dim, length=[4. * np.pi, ] * dim) + topos[case] = doms[case].create_topology(disc[case]) + refsize = (np.prod(topos[case].mesh.compute_resolution) + 324, ) + iworks[case] = [npw.int_zeros(refsize)] + rworks[case] = [npw.zeros(refsize)] +dx = topos[2].mesh.space_step + + +def run_interp(dimension, velo=0., with_work=False): + """Interpolate a field and compare with reference + (cosine function), for a given domain dimension + """ + topo = topos[dimension - 1] + rwork, iwork = None, None + if with_work: + rwork = rworks[dimension - 1] + iwork = iworks[dimension - 1] + ind = topo.mesh.compute_index + field = npw.zeros(topo.mesh.resolution) + for direction in xrange(topo.domain.dimension): + # lref = doms[dimension - 1].length[direction] + h = topo.mesh.space_step[direction] + x = topo.mesh.coords[direction] + field[...] = np.cos(x) + # select a set of points which are advected randomly and interpolated + # --> 1 point over two from the original mesh + shape = list(topo.mesh.resolution) + coords = topo.mesh.coords[direction] + ic = [slice(0, coords.shape[i]) for i in xrange(topo.domain.dimension)] + ic[direction] = slice(0, topo.mesh.resolution[direction] - 3, 2) + shape[direction] = topo.mesh.coords[direction][ic].shape[direction] + positions = npw.random(shape) * velo + positions[...] += topo.mesh.coords[direction][ic] + interp = Linear(field, direction=direction, topo_source=topo, + rwork=rwork, iwork=iwork) + result = [npw.zeros_like(positions)] + result = interp(0., [positions], result) + reference = npw.zeros_like(result[0]) + reference[...] = np.cos(positions) + tol = 0.5 * h ** 2 + assert np.allclose(result[0][ind], reference[ind], rtol=tol) + + +def test_linear_interpolation_1d_velocity_null(): + """Linear interpolation on 1d domain""" + run_interp(1) + + +def test_linear_interpolation_1d(): + """Linear interpolation on 1d domain""" + run_interp(1, 0.5 * dx[0]) + + +def test_linear_interpolation_2d_velocity_null(): + """Linear interpolation on 2d domain""" + run_interp(2) + + +def test_linear_interpolation_2d(): + """Linear interpolation on 2d domain""" + run_interp(2, 0.5 * dx[1]) + + +def test_linear_interpolation_3d_velocity_null(): + """Linear interpolation on 3d domain""" + run_interp(3) + + +def test_linear_interpolation_3d(): + """Linear interpolation on 3d domain""" + run_interp(3, 0.5 * dx[2]) + + +def test_linear_interpolation_1d_wk_velocity_null(): + """Linear interpolation on 1d domain, external work arrays""" + run_interp(1, 0., True) + + +def test_linear_interpolation_1d_wk(): + """Linear interpolation on 1d domain, external work arrays""" + run_interp(1, 0.5 * dx[0], True) + + +def test_linear_interpolation_2d_wk_velocity_null(): + """Linear interpolation on 2d domain, external work arrays""" + run_interp(2, 0., True) + + +def test_linear_interpolation_2d_wk(): + """Linear interpolation on 2d domain, external work arrays""" + run_interp(2, 0.5 * dx[1], True) + + +def test_linear_interpolation_3d_wk_velocity_null(): + """Linear interpolation on 3d domain, external work arrays""" + run_interp(3, 0., True) + + +def test_linear_interpolation_3d_wk(): + """Linear interpolation on 3d domain, external work arrays""" + run_interp(3, 0.5 * dx[2], True) diff --git a/hysop/numerics/tests/test_remesh.py b/hysop/numerics/tests/test_remesh.py new file mode 100644 index 0000000000000000000000000000000000000000..8ff48ed3f5bbda67628a95201156a848af09f593 --- /dev/null +++ b/hysop/numerics/tests/test_remesh.py @@ -0,0 +1,223 @@ +"""Test remesh schemes +""" +from hysop.numerics.remeshing import Linear, L2_1, Remeshing +import hysop.tools.numpywrappers as npw +from hysop.constants import EPS +from hysop import Box, Discretization +import numpy as np +cos = np.cos + +disc = [None, ] * 3 +g = 2 +nx = 65 +ny = 49 +nz = 57 +disc[0] = Discretization([nx, ], [g, ]) +disc[1] = Discretization([nx, ny], [g, g]) +disc[2] = Discretization([nx, ny, nz], [g, g, g]) +doms = [None, ] * 3 +topos = [None, ] * 3 +rworks = [None, ] * 3 +iworks = [None, ] * 3 +for case in xrange(3): + dim = case + 1 + doms[case] = Box(origin=[0.] * dim, length=[4. * np.pi, ] * dim) + topos[case] = doms[case].create_topology(disc[case]) + refsize = (np.prod(topos[case].mesh.resolution) + 324, ) + iworks[case] = [npw.int_zeros(refsize)] + rworks[case] = [npw.zeros(refsize), npw.zeros(refsize)] +dx = [0., 0. , 0.] # topos[2].mesh.space_step + + +def run_remesh(kernel, dimension, velo=0., with_work=False): + """Interpolate a field and compare with reference + (cosine function), for a given domain dimension + """ + topo = topos[dimension - 1] + rwork, iwork = None, None + if with_work: + rwork = rworks[dimension - 1] + iwork = iworks[dimension - 1] + ind = topo.mesh.compute_index + lenfield = 2 + fieldp = [None, ] * lenfield + result = [None, ] * lenfield + reference = [None, ] * lenfield + for direction in xrange(topo.domain.dimension): + h = topo.mesh.space_step[direction] + x = topo.mesh.coords[direction] + reference = [npw.zeros(topo.mesh.resolution) + for _ in xrange(lenfield)] + reference[0][...] = np.cos(x) + if lenfield > 1: + reference[1][...] = np.sin(x) + shape = topo.mesh.resolution.copy() + shape[direction] -= 2 * g + # note : this should work for any number of points + # in direction, as soon as work array is large enough + # + EPS to be sure that particles are (almost) on grid + # but on the right of the grid point when velo=0 + positions = npw.random(shape) * velo + positions[...] += 10 * EPS + positions[...] += topo.mesh.compute_coords[direction] + for i in xrange(lenfield): + fieldp[i] = npw.zeros_like(positions) + result[i] = npw.zeros_like(reference[i]) + fieldp[0][...] = np.cos(positions) + if lenfield > 1: + fieldp[1][...] = np.sin(positions) + remesh = Remeshing(kernel, topo_source=topo, direction=direction, + rwork=rwork, iwork=iwork) + result = remesh(positions, fieldp, result) + tol = 0.5 * h ** 2 + #plt.plot(x.flat[2:-2], result[ind], 'o-', x.flat[2:-2], reference[ind], '+-') + #plt.show() + for i in xrange(lenfield): + assert np.allclose(result[i][ind], reference[i][ind], rtol=tol) + + +def test_linear_remesh_1d_velocity_null(): + """Linear remesh on 1d domain""" + run_remesh(Linear, 1) + + +def test_linear_remesh_1d(): + """Linear remesh on 1d domain""" + run_remesh(Linear, 1, 0.3 * dx[0]) + + +def test_linear_remesh_2d_velocity_null(): + """Linear remesh on 2d domain""" + run_remesh(Linear, 2) + + +def test_linear_remesh_2d(): + """Linear remesh on 2d domain""" + run_remesh(Linear, 2, 0.3 * dx[1]) + + +def test_linear_remesh_3d_velocity_null(): + """Linear remesh on 3d domain""" + run_remesh(Linear, 3) + + +def test_linear_remesh_3d(): + """Linear remesh on 3d domain""" + run_remesh(Linear, 3, 0.3 * dx[2]) + + +def test_linear_remesh_1d_wk_velocity_null(): + """Linear remesh on 1d domain, external work arrays""" + run_remesh(Linear, 1, 0., True) + + +def test_linear_remesh_1d_wk(): + """Linear remesh on 1d domain, external work arrays""" + run_remesh(Linear, 1, 0.3 * dx[0], True) + + +def test_linear_remesh_2d_wk_velocity_null(): + """Linear remesh on 2d domain, external work arrays""" + run_remesh(Linear, 2, 0., True) + + +def test_linear_remesh_2d_wk(): + """Linear remesh on 2d domain, external work arrays""" + run_remesh(Linear, 2, 0.3 * dx[1], True) + + +def test_linear_remesh_3d_wk_velocity_null(): + """Linear remesh on 3d domain, external work arrays""" + run_remesh(Linear, 3, 0., True) + + +def test_linear_remesh_3d_wk(): + """Linear remesh on 3d domain, external work arrays""" + run_remesh(Linear, 3, 0.3 * dx[2], True) + + +def test_l21_remesh_1d_velocity_null(): + """L2_1 remesh on 1d domain""" + run_remesh(L2_1, 1) + + +def test_l21_remesh_1d(): + """L2_1 remesh on 1d domain""" + run_remesh(L2_1, 1, 0.3 * dx[0]) + + +def test_l21_remesh_2d_velocity_null(): + """L2_1 remesh on 2d domain""" + run_remesh(L2_1, 2) + + +def test_l21_remesh_2d(): + """L2_1 remesh on 2d domain""" + run_remesh(L2_1, 2, 0.3 * dx[1]) + + +def test_l21_remesh_3d_velocity_null(): + """L2_1 remesh on 3d domain""" + run_remesh(L2_1, 3) + + +def test_l21_remesh_3d(): + """L2_1 remesh on 3d domain""" + run_remesh(L2_1, 3, 0.3 * dx[2]) + + +def test_l21_remesh_1d_wk_velocity_null(): + """L2_1 remesh on 1d domain, external work arrays""" + run_remesh(L2_1, 1, 0., True) + + +def test_l21_remesh_1d_wk(): + """L2_1 remesh on 1d domain, external work arrays""" + run_remesh(L2_1, 1, 0.3 * dx[0], True) + + +def test_l21_remesh_2d_wk_velocity_null(): + """L2_1 remesh on 2d domain, external work arrays""" + run_remesh(L2_1, 2, 0., True) + + +def test_l21_remesh_2d_wk(): + """L2_1 remesh on 2d domain, external work arrays""" + run_remesh(L2_1, 2, 0.3 * dx[1], True) + + +def test_l21_remesh_3d_wk_velocity_null(): + """L2_1 remesh on 3d domain, external work arrays""" + run_remesh(L2_1, 3, 0., True) + + +def test_l21_remesh_3d_wk(): + """L2_1 remesh on 3d domain, external work arrays""" + run_remesh(L2_1, 3, 0.3 * dx[2], True) + + +if __name__ == "__main__": + test_linear_remesh_1d() + test_linear_remesh_1d_velocity_null() + test_linear_remesh_2d() + test_linear_remesh_2d_velocity_null() + test_linear_remesh_3d() + test_linear_remesh_3d_velocity_null() + test_linear_remesh_1d_wk() + test_linear_remesh_1d_wk_velocity_null() + test_linear_remesh_2d_wk() + test_linear_remesh_2d_wk_velocity_null() + test_linear_remesh_3d_wk() + test_linear_remesh_3d_wk_velocity_null() + test_l21_remesh_1d() + test_l21_remesh_1d_velocity_null() + test_l21_remesh_2d() + test_l21_remesh_2d_velocity_null() + test_l21_remesh_3d() + test_l21_remesh_3d_velocity_null() + test_l21_remesh_1d_wk() + test_l21_remesh_1d_wk_velocity_null() + test_l21_remesh_2d_wk() + test_l21_remesh_2d_wk_velocity_null() + test_l21_remesh_3d_wk() + test_l21_remesh_3d_wk_velocity_null() diff --git a/hysop/operator/advection.py b/hysop/operator/advection.py index 85883c51a635a2432e4c56834ed25ae5b2f13271..fdb7b99cec0353afb646e0a9a5f11ba6edb8c8dc 100644 --- a/hysop/operator/advection.py +++ b/hysop/operator/advection.py @@ -2,143 +2,409 @@ See :ref:`advection` in user guide. """ +from __future__ import print_function +from abc import ABCMeta, abstractmethod from hysop.constants import debug, S_DIR, ZDIR +from hysop import __SCALES_ENABLED__, Field, __GPU_ENABLED__ from hysop.operator.computational import Computational -from hysop.methods_keys import Scales, TimeIntegrator, Interpolation,\ +from hysop.methods_keys import TimeIntegrator, Interpolation,\ Remesh, Support, Splitting, MultiScale -from hysop.numerics.remeshing import L2_1 +from hysop.numerics.remeshing import Remeshing, L2_1 from hysop.operator.continuous import opsetup, opapply -from hysop.operator.advection_dir import AdvectionDir import hysop.default_methods as default from hysop.tools.parameters import Discretization from hysop.mpi.topology import Cartesian import hysop.tools.numpywrappers as npw - - -class Advection(Computational): - """ - Advection of a field, - \f{eqnarray*} - X = Op(X,velocity) - \f} for : - \f{eqnarray*} - \frac{\partial{X}}{\partial{t}} + velocity.\nabla X = 0 - \f} - Note : we assume incompressible flow. - - Computations are performed within a dimensional splitting as folows: - - 2nd order: - - X-dir, half time step - - Y-dir, half time step - - Z-dir, full time step - - Y-dir, half time step - - X-dir, half time step - - 2nd order full half-steps: - - X-dir, half time step - - Y-dir, half time step - - Z-dir, half time step - - Z-dir, half time step - - Y-dir, half time step - - X-dir, half time step - - 1st order g: - - X-dir, half time step - - Y-dir, half time step - - Z-dir, half time step - +from hysop.problem.simulation import O2, SplittingParameters +from hysop.operator.discrete.particle_advection import ParticleAdvection +import numpy as np +if __GPU_ENABLED__: + from hysop.gpu.gpu_particle_advection import GPUParticleAdvection + from hysop.gpu.multi_gpu_particle_advection import MultiGPUParticleAdvection +if __SCALES_ENABLED__: + from hysop.f2hysop import scales2py as scales + + +class AdvectionBase(Computational): + """Abstract interface to advection operators """ + __metaclass__ = ABCMeta @debug - def __init__(self, velocity, advected_fields=None, **kwds): + def __init__(self, velocity, advected_fields=None, + discretization_fields=None, **kwds): """ - Advection of a set of fields for a given velocity. + Parameters + ---------- + velocity : :class:`hysop.field.continuous.Field` + the velocity field + advected_field: a list of :class:`hysop.field.continuous.Field`, + optional + fields to be advected + discretization_fields : :class:`hysop.mpi.topology.Cartesian` + or :class:`tools.parameters.Discretization` + Defined the data/mpi distribution for advected_fields. + Default = same as velocity + kwds : base class arguments. + + Notes + ----- + * one single discretization for all advected fields + * velocity discretization may be different + * to limit advection to one or two directions, set + directions=[0, 2] (advection in x and z dirs). - @param velocity : velocity field used for advection - @param advectedFields : the list of fields to be advected. - It may be a single field (no list). """ - ## Transport velocity - self.velocity = velocity - if 'variables' in kwds: - kw = kwds.copy() - kw['variables'] = kwds['variables'].copy() - # In that case, variables must contains only the advected fields - # with their discretization param. - # Velocity must always be given outside variables, with its - # own discretization. - assert advected_fields is None, 'too many input arguments.' - self.advected_fields = kwds['variables'].keys() - kw['variables'][self.velocity] = kwds['discretization'] - kw.pop('discretization') - super(Advection, self).__init__(**kw) - + # --- check inputs for fields (advected and velocity) --- + msg = 'discretization argument is required in Advection.' + assert 'discretization' in kwds, msg + msg = 'Wrong input for advection operator:' + msg += ' "variables" attribute is not allowed. Check manual.' + assert 'variables' not in kwds, msg + # advected_fields : + # == field or [field1, field2, ...] ==> same discretization for all, + # equal to velocity discretization if discretization_fields is None + self.advected_fields = [] + if advected_fields is not None: + if not isinstance(advected_fields, list): + assert isinstance(advected_fields, Field) + advected_fields = [advected_fields] + self.advected_fields = advected_fields + # velocity field + if discretization_fields is not None: + dfields = discretization_fields else: - v = [self.velocity] - if isinstance(advected_fields, list): - self.advected_fields = advected_fields - else: - self.advected_fields = [advected_fields] - v += self.advected_fields - super(Advection, self).__init__(variables=v, **kwds) + dfields = kwds['discretization'] + variables = {velocity: kwds['discretization']} + for field in advected_fields: + variables[field] = dfields + kwds.pop('discretization') + super(AdvectionBase, self).__init__(variables=variables, **kwds) + # velocity field + self.velocity = velocity + assert velocity.nb_components == self.domain.dimension - # Set default method, if required + # --- Set default method, if required --- if self.method is None: self.method = default.ADVECTION - self.output = self.advected_fields - self.input = [var for var in self.variables] - + # --- Set name --- vars_str = "_(" for vv in self.advected_fields: vars_str += vv.name + "," self.name += vars_str[0:-1] + ')' self.config = {} + self.input = [var for var in self.variables] + self.output = [var for var in self.variables + if var is not self.velocity] + + def get_work_properties(self): + return {'rwork': None, 'iwork': None} + + @opsetup + def setup(self, rwork=None, iwork=None): + if self._is_uptodate: + return + # select discretization of the advected fields + advected_discrete_fields = [self.discreteFields[v] + for v in self.variables + if v is not self.velocity] + toporef = advected_discrete_fields[0].topology + msg = 'All advected fields must have the same topology.' + for f in advected_discrete_fields: + assert f.topology == toporef, msg + + # call specific setup (see derived class) + self._setup(advected_discrete_fields, rwork, iwork) + self._is_uptodate = True + + def advected_fields_topology(self): + """Returns the topology used for advected fields + """ + msg = 'Operator must be discretized' + assert self._is_discretized, msg + return self.discreteFields[self.advected_fields[0]].topology + + def velocity_topology(self): + """Returns the topology used for velocity + """ + msg = 'Operator must be discretized' + assert self._is_discretized, msg + return self.discreteFields[self.velocity].topology + + @abstractmethod + def _setup(self, advected_discrete_fields, rwork=None, iwork=None): + """Setup specific to advection implementation (see derived class) + """ + - # Find which solver is used for advection, - # among Scales, pure-python and GPU-like. - # Check also operator-splitting type. +class Advection(AdvectionBase): + """Python or GPU advection + """ + @debug + def __init__(self, directions=None, **kwds): + """ + Parameters + ---------- + directions : list of int, optional + direction(s) in which advection must be performed. Default = all. + kwds : base class arguments. + + """ + super(Advection, self).__init__(**kwds) + self.extra_args = {} ### TO BE REVIEWED ### + if directions is None: + directions = [i for i in xrange(self.domain.dimension)] + assert isinstance(directions, list) + self.directions = directions + assert TimeIntegrator in self.method.keys() + assert Interpolation in self.method.keys() + assert Remesh in self.method.keys() + if Splitting not in self.method.keys(): + self.method[Splitting] = O2 + self.discrete_op = [None] * self.domain.dimension + + # Fields on particles + self.particle_fields = None + + # Positions of the particles + self.particle_positions = None + + # number of components of the rhs (time integration) + self._rhs_size = 1 + + # check if advection is done on gpu + if Support in self.method: + self._gpu_advection = self.method[Support].find('gpu') >= 0 + else: + self._gpu_advection = False + self.splitting = None + self._old_dir = 0 - # ---- Scales advection ---- - if Scales in self.method.keys(): - self._is_scales = True - if not self.domain.dimension == 3: - raise ValueError("Scales Advection not implemented in 2D.") - # Default splitting = Strang - if Splitting not in self.method.keys(): - self.method[Splitting] = 'strang' + def discretize(self): + """Discretisation (create topologies and discretize fields) + """ + if self._is_discretized: + return + cutdir = [False] * self.domain.dimension + cutdir[-1] = True + # call standard discretization (from computational base class) + self._standard_discretize(cutdir=cutdir) + # check if velocity and fields are defined with the same resolution + if self._single_topo: + self.method[MultiScale] = None + else: + # set a default value for interpolation method + if MultiScale not in self.method or\ + self.method[MultiScale] is None: + self.method[MultiScale] = L2_1 + mscale = self.method[MultiScale] + min_ghosts = mscale.ghosts_layer_size + # get ghost layer size used for velocity field + ghosts_v = self.variables[self.velocity].ghosts() + # and check if it fits with interpolation method + msg = 'Ghost layer required for velocity. Size min = ' + msg += str(min_ghosts) + " (" + str(ghosts_v) + " given)" + assert (ghosts_v >= min_ghosts).all(), msg + self._is_discretized = True - self._my_setup = self._setup_scales - self.advec_dir = None + def get_work_properties(self): + assert self._is_discretized + # Shape of reference comes from fields, not from velocity + advected_discrete_fields = [self.discreteFields[v] + for v in self.variables + if v is not self.velocity] + topo = advected_discrete_fields[0].topology + # Find number and shape of required work arrays + if not self._gpu_advection: + # -- pure python advection -- + # work array shape depends on the time integrator + # interpolation scheme and remeshing scheme + ti_work = self.method[TimeIntegrator].get_work_properties( + self._rhs_size, topo) + ti_rwork_length = len(ti_work['rwork']) + iw_prop = self.method[Interpolation].get_work_properties(topo) + rw_prop = Remeshing.get_work_properties(topo) + interp_iwork_length = len(iw_prop['iwork']) + interp_rwork_length = len(iw_prop['rwork']) + remesh_iwork_length = len(rw_prop['iwork']) + remesh_rwork_length = len(rw_prop['rwork']) + iwork_length = max(interp_iwork_length, remesh_iwork_length) + rwork_length = max(ti_rwork_length + interp_rwork_length, + remesh_rwork_length) else: - # ---- Python or GPU advection ---- - self._is_scales = False - assert TimeIntegrator in self.method.keys() - assert Interpolation in self.method.keys() - assert Remesh in self.method.keys() - assert Support in self.method.keys() - if Splitting not in self.method.keys(): - self.method[Splitting] = 'o2' - dimension = self.domain.dimension - self.advec_dir = [None] * dimension - name = vars_str[0:-1] + ')' - if 'variables' in kwds: - for i in xrange(self.domain.dimension): - self.advec_dir[i] = AdvectionDir( - self.velocity, direction=i, - name_suffix=name, **kwds) + # -- GPU advection -- + # no work array + iwork_length, rwork_length = 0, 0 + + # buffers for fields on particles + rwork_length += np.sum([f.nb_components for f in self.advected_fields]) + if not self._gpu_advection or \ + self.method[Support].find('gpu_2k') >= 0: + rwork_length += 1 # work array for positions + memsize = np.prod(topo.mesh.resolution) + return {'rwork': [(memsize,)] * rwork_length, + 'iwork': [(memsize,)] * iwork_length} + + def _setup(self, advected_discrete_fields, rwork=None, iwork=None): + for i in self.directions: + ref_shape = advected_discrete_fields[0].topology.shape[i] + if self._gpu_advection: + if ref_shape == 1: + DA = GPUParticleAdvection + else: + DA = MultiGPUParticleAdvection else: - for i in xrange(self.domain.dimension): - self.advec_dir[i] = AdvectionDir( - self.velocity, direction=i, - advected_fields=self.advected_fields, - name_suffix=name, **kwds) + DA = ParticleAdvection - self._my_setup = self._setup_python - self.apply = self._apply_python + self.discrete_op[i] = DA( + velocity=self.discreteFields[self.velocity], + fields_on_grid=advected_discrete_fields, + direction=i, method=self.method, + rwork=rwork, iwork=iwork) + #**self.extra_args) - self._old_dir = 0 - self.splitting = [] + if i == 0: + # work arrays can be shared between directions. + rwork = self.discrete_op[i]._rwork + iwork = self.discrete_op[i]._iwork + + # set splitting parameters (depends on method) + self.splitting = self.method[Splitting](self.domain.dimension) + assert isinstance(self.splitting, SplittingParameters) + + # configure gpu + if self._gpu_advection: + self._configure_gpu() + + def _configure_gpu(self): + """Setup for gpu related things (Memory ...) + """ + splitting_nbSteps = len(self.splitting) + for d in xrange(self.domain.dimension): + dOp = self.discrete_op[d] + assert len(dOp.exec_list) == splitting_nbSteps, \ + "Discrete operator execution " + \ + "list and splitting steps sizes must be equal " + \ + str(len(dOp.exec_list)) + " != " + \ + str(splitting_nbSteps) + s = "" + device_id = self.discrete_op[0].cl_env._device_id + gpu_comm = self.discrete_op[0].cl_env.gpu_comm + gpu_rank = gpu_comm.Get_rank() + if gpu_rank == 0: + s += "=== OpenCL buffers allocated" + s += " on Device:{0} ===\n".format(device_id) + s += "Global memory used:\n" + total_gmem = 0 + for d in xrange(self.domain.dimension): + g_mem_d = 0 + # allocate all variables in advec_dir + for df in self.discrete_op[d].variables: + if not df.gpu_allocated: + df.allocate() + g_mem_df = gpu_comm.allreduce(df.mem_size) + g_mem_d += g_mem_df + if gpu_rank == 0: + s += " Advection" + S_DIR[d] + ": {0:9d}".format(g_mem_d) + s += "Bytes ({0:5d} MB)\n".format(g_mem_d / (1024 ** 2)) + total_gmem += g_mem_d + if gpu_rank == 0: + s += " Total : {0:9d}".format(total_gmem) + s += "Bytes ({0:5d} MB)\n".format(total_gmem / (1024 ** 2)) + s += "Local memory used:\n" + total_lmem = 0 + for d in xrange(self.domain.dimension): + l_mem_d = gpu_comm.allreduce( + self.discrete_op[d].size_local_alloc) + if gpu_rank == 0: + s += " Advection" + S_DIR[d] + ": {0:9d}".format(l_mem_d) + s += "Bytes ({0:5d} MB)\n".format(l_mem_d / (1024 ** 2)) + total_lmem += l_mem_d + if gpu_rank == 0: + s += " Total : {0:9d}".format(total_lmem) + "Bytes" + print(s) + + @debug + @opapply + def apply(self, simulation=None): + """Redefinition of apply for advection --> dimensional splitting. + + Parameters + ---------- + simulation : `:class::~hysop.problem.simulation.Simulation` + + """ + assert simulation is not None + for split_id, split in enumerate(self.splitting()): + simulation.set_split(split_id, split) + self.discrete_op[simulation.current_dir].apply(simulation) + #, split[1], split_id, self._old_dir) + simulation.next_split() + + @debug + def finalize(self): + """Memory cleaning. + """ + for i in self.directions: + self.discrete_op[i].finalize() + + def get_profiling_info(self): + if self._is_uptodate: + for d in self.directions: + self.profiler += self.discrete_op[d].profiler + + def __str__(self): + """ + Common printings for operators + """ + short_name = str(self.__class__).rpartition('.')[-1][0:-2] + for i in self.directions: + if self.discrete_op[i] is not None: + s = str(self.discrete_op[i]) + else: + s = short_name + " operator. Not discretised." + return s + "\n" + + +class ScalesAdvection(AdvectionBase): + """Advection based on fortran (scales) interface + """ + @debug + def __init__(self, cutdir=None, **kwds): + """ + Parameters + ---------- + cutdir : list of bool + cutdir[d] = True to distribute data in direction d. + + kwds : base class arguments. + + Notes + ----- + * In scales cutdir[0] must be False, i.e. data are not + distributed in the first dir, which corresponds to contiguous + memory layout (fortran). + * Scales assumes that all subdomains have the same + local resolution. + * No ghosts points allowed in scales (at the moment) + + """ + super(ScalesAdvection, self).__init__(**kwds) + assert Remesh in self.method.keys() + msg = 'Scales is not available for your configuration.' + assert __SCALES_ENABLED__, msg + msg = 'Scales Advection is only implemented in 3D.' + assert self.domain.dimension == 3, msg + # Default splitting = Strang + if Splitting not in self.method.keys(): + self.method[Splitting] = 'strang' + if cutdir is None: + cutdir = [False, ] * self.domain.dimension + cutdir[-1] = True + self.cutdir = cutdir def scales_parameters(self): """ @@ -149,10 +415,10 @@ class Advection(Computational): for o in ['p_O2', 'p_O4', 'p_L2', 'p_M4', 'p_M6', 'p_M8', 'p_44', 'p_64', 'p_66', 'p_84']: - if self.method[Scales].find(o) >= 0: + if self.method[Remesh].find(o) >= 0: order = o if order is None: - print ('Unknown advection method, turn to default (p_M6).') + print('Unknown advection method, turn to default (p_M6).') order = 'p_M6' # - Extract splitting form self.method (default strang) - @@ -174,47 +440,9 @@ class Advection(Computational): """ if self._is_discretized: return - - if self._is_scales: - self._scales_discretize() - else: - self._no_scales_discretize() - advected_discrete_fields = [self.discreteFields[f] - for f in self.advected_fields] - toporef = advected_discrete_fields[0].topology - msg = 'All advected fields must have the same topology.' - for f in advected_discrete_fields: - assert f.topology == toporef, msg - - if self._single_topo: - self.method[MultiScale] = None - - @debug - def _scales_discretize(self): - """ - Discretization (create topologies and discretize fields) - when using SCALES fortran routines (3d only, list of vector - and/or scalar) - - 'p_O2' : order 4 method, corrected to allow large CFL number, - untagged particles - - 'p_O4' : order 4 method, corrected to allow large CFL number, - untagged particles - - 'p_L2' : limited and corrected lambda 2 - - 'p_M4' : Lambda_2,1 (=M'4) 4 point formula - - 'p_M6' (default) : Lambda_4,2 (=M'6) 6 point formula - - 'p_M8' : M8prime formula - - 'p_44' : Lambda_4,4 formula - - 'p_64' : Lambda_6,4 formula - - 'p_66' : Lambda_6,6 formula - - 'p_84' : Lambda_8,4 formula - """ - assert self._is_scales - # - Extract order form self.method (default p_M6) - - order, splitting = self.scales_parameters() - # Check if topos need to be created build_topos = self._check_variables() - from hysop.f2hysop import scales2py as scales + order, splitting = self.scales_parameters() # Scales, single resolution if self._single_topo: @@ -279,11 +507,23 @@ class Advection(Computational): # All topos are built, we can discretize fields. self._discretize_vars() + advected_discrete_fields = [self.discreteFields[f] + for f in self.advected_fields] + toporef = advected_discrete_fields[0].topology + msg = 'All advected fields must have the same topology.' + for f in advected_discrete_fields: + assert f.topology == toporef, msg + + if self._single_topo: + self.method[MultiScale] = None + def _create_scales_topo(self, d3d, order, splitting): - from hysop.f2hysop import scales2py as scales + """set specific MPI layout for scales. + """ comm = self._mpis.comm topodims = [1, 1, comm.Get_size()] - msg = 'Wrong type for parameter discretization (at init).' + str(self._discretization) + msg = 'Wrong type for parameter discretization (at init).' + msg += str(self._discretization) assert isinstance(d3d, Discretization), msg nbcells = d3d.resolution - 1 scalesres, global_start = \ @@ -299,7 +539,8 @@ class Advection(Computational): discretization=d3d, cdir=ZDIR) def _check_scales_topo(self, toporef, order, splitting): - from hysop.f2hysop import scales2py as scales + """Check if input topo fits with scales requirements + """ # In that case, self._discretization must be # a Cartesian object, used for all fields. # We use it to initialize scales solver @@ -317,228 +558,23 @@ class Advection(Computational): assert isinstance(topo, Cartesian), str(topo) #assert (topo.shape == topodims).all(), \ # str(topo.shape) + ' != ' + str(topodims) - assert not self._single_topo or (topo.mesh.resolution == scalesres).all(), \ + assert not self._single_topo or \ + (topo.mesh.resolution == scalesres).all(), \ str(topo.mesh.resolution) + ' != ' + str(scalesres) - assert not self._single_topo or (topo.mesh.start() == global_start).all(), \ + assert not self._single_topo or \ + (topo.mesh.start() == global_start).all(), \ str(topo.mesh.start()) + ' != ' + str(global_start) - def _no_scales_discretize(self): - """ - GPU or pure-python advection - """ - if not self._is_discretized: - for i in xrange(self.domain.dimension): - self.advec_dir[i].discretize() - self.discreteFields = self.advec_dir[0].discreteFields - self._single_topo = self.advec_dir[0]._single_topo - self._is_discretized = True - - def get_work_properties(self): - """ - Return the length of working arrays lists required - for the discrete operator. - @return shapes, shape of the arrays: - shapes['rwork'] == list of shapes for real arrays, - shapes['iwork'] == list of shapes for int arrays. - len(shapes['...'] gives the number of required arrays. - """ - if self._is_scales: - return {'rwork': None, 'iwork': None} - else: - if not self.advec_dir[0]._is_discretized: - msg = 'The operator must be discretized ' - msg += 'before any call to this function.' - raise RuntimeError(msg) - return self.advec_dir[0].get_work_properties() - - @opsetup - def setup(self, rwork=None, iwork=None): + def _setup(self, advected_discrete_fields, rwork=None, iwork=None): # Check resolutions to set multiscale case, if required. if not self._single_topo and MultiScale not in self.method: self.method[MultiScale] = L2_1 - if not self._is_uptodate: - self._my_setup(rwork, iwork) - - def _setup_scales(self, rwork=None, iwork=None): - - advected_discrete_fields = [self.discreteFields[f] - for f in self.advected_fields] # - Create the discrete_op from the # list of discrete fields - from hysop.operator.discrete.scales_advection import \ - ScalesAdvection - self.discrete_op = ScalesAdvection( + ScalesAdvection as SD + self.discrete_op = SD( self.discreteFields[self.velocity], advected_discrete_fields, method=self.method, rwork=rwork, iwork=iwork, **self.config) - self._is_uptodate = True - - def _setup_advec_dir(self, rwork=None, iwork=None): - """ - Local allocation of work arrays, - common to advec_dir operators and setup for those - operators - """ - wk_p = self.advec_dir[0].get_work_properties() - wk_length = len(wk_p['rwork']) - if rwork is None: - rwork = [] - for i in xrange(wk_length): - memshape = wk_p['rwork'][i] - rwork.append(npw.zeros(memshape)) - else: - assert len(rwork) == wk_length - for wk, refshape in zip(rwork, wk_p['rwork']): - assert wk.shape == refshape - wk_length = len(wk_p['iwork']) - if iwork is None: - iwork = [] - for i in xrange(wk_length): - memshape = wk_p['iwork'][i] - iwork.append(npw.int_zeros(memshape)) - else: - assert len(iwork) == wk_length - for wk, refshape in zip(iwork, wk_p['iwork']): - assert wk.shape == refshape - # Work arrays are common between all directions - # of advection. - for i in xrange(self.domain.dimension): - self.advec_dir[i].setup(rwork, iwork) - - def _setup_python(self, rwork=None, iwork=None): - - # setup for advection in each direction - self._setup_advec_dir(rwork, iwork) - # set splitting parameters (depends on method) - self._configure_splitting() - - # configure gpu - if self.method[Support].find('gpu') >= 0: - self._configure_gpu() - - self._is_uptodate = True - - def _configure_splitting(self): - dimension = self.domain.dimension - if self.method[Splitting] == 'o2_FullHalf': - ## Half timestep in all directions - [self.splitting.append((i, 0.5)) - for i in xrange(dimension)] - [self.splitting.append((dimension - 1 - i, 0.5)) - for i in xrange(dimension)] - elif self.method[Splitting] == 'o1': - [self.splitting.append((i, 1.)) for i in xrange(dimension)] - elif self.method[Splitting] == 'x_only': - self.splitting.append((0, 1.)) - elif self.method[Splitting] == 'y_only': - self.splitting.append((1, 1.)) - elif self.method[Splitting] == 'z_only': - self.splitting.append((2, 1.)) - elif self.method[Splitting] == 'o2': - ## Half timestep in all directions but last - [self.splitting.append((i, 0.5)) - for i in xrange(dimension - 1)] - self.splitting.append((dimension - 1, 1.)) - [self.splitting.append((dimension - 2 - i, 0.5)) - for i in xrange(dimension - 1)] - else: - raise ValueError('Unknown splitting configuration:' + - self.method[Splitting]) - - def _configure_gpu(self): - splitting_nbSteps = len(self.splitting) - for d in xrange(self.domain.dimension): - dOp = self.advec_dir[d].discrete_op - assert len(dOp.exec_list) == splitting_nbSteps, \ - "Discrete operator execution " + \ - "list and splitting steps sizes must be equal " + \ - str(len(dOp.exec_list)) + " != " + \ - str(splitting_nbSteps) - s = "" - device_id = self.advec_dir[0].discrete_op.cl_env._device_id - gpu_comm = self.advec_dir[0].discrete_op.cl_env.gpu_comm - gpu_rank = gpu_comm.Get_rank() - if gpu_rank == 0: - s += "=== OpenCL buffers allocated" - s += " on Device:{0} ===\n".format(device_id) - s += "Global memory used:\n" - total_gmem = 0 - for d in xrange(self.domain.dimension): - g_mem_d = 0 - # allocate all variables in advec_dir - for df in self.advec_dir[d].discrete_op.variables: - if not df.gpu_allocated: - df.allocate() - g_mem_df = gpu_comm.allreduce(df.mem_size) - g_mem_d += g_mem_df - if gpu_rank == 0: - s += " Advection" + S_DIR[d] + ": {0:9d}".format(g_mem_d) - s += "Bytes ({0:5d} MB)\n".format(g_mem_d / (1024 ** 2)) - total_gmem += g_mem_d - if gpu_rank == 0: - s += " Total : {0:9d}".format(total_gmem) - s += "Bytes ({0:5d} MB)\n".format(total_gmem / (1024 ** 2)) - s += "Local memory used:\n" - total_lmem = 0 - for d in xrange(self.domain.dimension): - l_mem_d = gpu_comm.allreduce( - self.advec_dir[d].discrete_op.size_local_alloc) - if gpu_rank == 0: - s += " Advection" + S_DIR[d] + ": {0:9d}".format(l_mem_d) - s += "Bytes ({0:5d} MB)\n".format(l_mem_d / (1024 ** 2)) - total_lmem += l_mem_d - if gpu_rank == 0: - s += " Total : {0:9d}".format(total_lmem) + "Bytes" - print s - - @debug - @opapply - def _apply_python(self, simulation=None): - """ - Apply this operator to its variables. - @param simulation : object that describes the simulation - parameters (time, time step, iteration number ...), see - hysop.problem.simulation.Simulation for details. - - Redefinition for advection. Applying a dimensional splitting. - """ - assert simulation is not None - for split_id, split in enumerate(self.splitting): - self.advec_dir[split[0]].apply( - simulation, split[1], split_id, self._old_dir) - self._old_dir = split[0] - - @debug - def finalize(self): - """ - Memory cleaning. - """ - if self._is_scales: - Computational.finalize(self) - else: - for dop in self.advec_dir: - dop.finalize() - - def get_profiling_info(self): - if self._is_uptodate: - if self._is_scales: - self.profiler += self.discrete_op.profiler - else: - for dop in self.advec_dir: - self.profiler += dop.profiler - - def __str__(self): - """ - Common printings for operators - """ - shortName = str(self.__class__).rpartition('.')[-1][0:-2] - if self._is_scales: - super(Advection, self).__str__() - else: - for i in xrange(self.domain.dimension): - if self.advec_dir[i].discrete_op is not None: - s = str(self.advec_dir[i].discrete_op) - else: - s = shortName + " operator. Not discretised." - return s + "\n" diff --git a/hysop/operator/advection_dir.py b/hysop/operator/advection_dir.py index 513949dd19803f4fadcf51705b545677e46a923c..dc563c92415b147e80b13be67f64a957401d4fe9 100644 --- a/hysop/operator/advection_dir.py +++ b/hysop/operator/advection_dir.py @@ -1,162 +1,72 @@ -""" -@file advectionDir.py +"""Advection of a field in a given direction + +See :ref:`advection` in user guide. -Advection of a field in a single direction. """ from hysop.constants import debug, S_DIR -from hysop.methods_keys import Support, MultiScale, \ - TimeIntegrator, Interpolation -from hysop.numerics.remeshing import L2_1, L4_2, L4_4, Remeshing, Linear +from hysop.methods_keys import Support, TimeIntegrator, Interpolation +from hysop.numerics.remeshing import Remeshing from hysop.operator.computational import Computational -# To get default method values for advection: -import hysop.default_methods as default import numpy as np from hysop.operator.continuous import opsetup, opapply +from hysop.operator.advection import Advection class AdvectionDir(Computational): - """ - Advection of a field, - \f{eqnarray*} - X = Op(X,velocity) - \f} for : - \f{eqnarray*} - \frac{\partial{X}}{\partial{t}} + velocity.\nabla X = 0 - \f} - Note : we assume incompressible flow. - - Computations are performed in a given direction. + """Advection of a scalar or vector field in a given direction, + assuming incompressible flow. """ @debug - def __init__(self, velocity, direction, advected_fields=None, - name_suffix='', **kwds): - ## Get the other arguments to pass to the discrete operators - self._kwds = kwds.copy() - self._kwds.pop('discretization') - ## Transport velocity - self.velocity = velocity - if 'variables' in kwds: - self._kwds.pop('variables') - kw = kwds.copy() - kw['variables'] = kwds['variables'].copy() - # In that case, variables must contains only the advected fields - # with their discretization param. - # Velocity must always be given outside variables, with its - # own discretization. - assert advected_fields is None, 'too many input arguments.' - self.advected_fields = kwds['variables'].keys() - kw['variables'][self.velocity] = kwds['discretization'] - kw.pop('discretization') - super(AdvectionDir, self).__init__(**kw) - else: - v = [self.velocity] - if isinstance(advected_fields, list): - self.advected_fields = advected_fields - else: - self.advected_fields = [advected_fields] - v += self.advected_fields - super(AdvectionDir, self).__init__(variables=v, **kwds) - - # Set default method, if required - if self.method is None: - self.method = default.ADVECTION - self.output = self.advected_fields - self.input = [var for var in self.variables] - - from hysop.methods_keys import Splitting - if Splitting not in self.method.keys(): - self.method[Splitting] = 'o2' - self.name += name_suffix + S_DIR[direction] - - ## direction to advect + def __init__(self, parent, direction, **kwds): + """" + Parameters + ---------- + parent : :class:`hysop.operator.advection.Advection` + main advection operator + direction: int + direction of advection + kwds : base class arguments. + """ + assert isinstance(parent, Advection) + self.parent = parent + # advected_fields : + self.advected_fields = parent.advected_fields + # velocity field + self.velocity = parent.velocity + self._kwds = kwds + super(AdvectionDir, self).__init__(variables=parent.variables, **kwds) + + self.method = parent.method + self.input = parent.input + self.output = parent.output + self.name = parent.name + S_DIR[direction] + + # direction to advect self.direction = direction - ## Fields on particles + # Fields on particles self.particle_fields = None - ## Positions of the particles + # Positions of the particles self.particle_positions = None + # number of components of the rhs (time integration) + self._rhs_size = 1 + + # check if advection is done on gpu + self._gpu_advection = self.method[Support].find('gpu') >= 0 + @debug def discretize(self): if self._is_discretized: return - build_topos = self._check_variables() - - # Check if multiscale is available - if not self._single_topo: - if self.method[Support].find('gpu') < 0: - raise ValueError("Multiscale advection is not yet supported " - "in pure Python, use Scales or GPU.") - - ## Default topology cutdir for parallel advection - cutdir = [False] * self.domain.dimension - cutdir[-1] = True - - if self._single_topo: - # One topo for all fields ... - self.method[MultiScale] = None - if build_topos: - topo = self.domain.create_topology( - discretization=self._discretization, cutdir=cutdir) - for v in self.variables: - self.variables[v] = topo - else: - # Topo is already built, just check if it is 1D - topo = self.variables.values()[0] - msg = str(topo.cutdir) + ' != ' + str(cutdir) - assert (topo.cutdir == cutdir).all(), msg - - else: - # ... or one topo for each field. - for v in self.variables: - if build_topos[v]: - topo = self.domain.create_topology( - discretization=self.variables[v], cutdir=cutdir) - self.variables[v] = topo - build_topos[v] = False - # compute velocity minimal ghost layer size - self._check_ghost_layer(build_topos) - - # All topos are built, we can discretize fields. - self._discretize_vars() - + # everything is done in parent ... + self.variables = self.parent.variables + self.discreteFields = self.parent.discreteFields self._is_discretized = True - def _check_ghost_layer(self, build_topos): - """ - Only meaningful if fields have different resolutions. - Check/set interpolation method for multiscale and - set ghost layer size, if required. - """ - # Set method to default if unknown - if MultiScale not in self.method.keys(): - self.method[MultiScale] = L2_1 - - mscale = self.method[MultiScale] - if mscale == Linear: - min_ghosts = 1 - elif mscale == L2_1: - min_ghosts = 2 - elif mscale == L4_2 or mscale == L4_4: - min_ghosts = 3 - else: - raise ValueError("Unknown multiscale method") - - # Topo or resolution associated with velocity - discr_v = self.variables[self.velocity] - if build_topos[self.velocity]: - # discr_v = Discretization - ghosts_v = discr_v.ghosts - else: - # discr_v = Cartesian - ghosts_v = discr_v.ghosts() - msg = 'Ghost layer required for velocity. Size min = ' - msg += str(min_ghosts) + " (" + str(ghosts_v) + " given)" - assert (ghosts_v >= min_ghosts).all(), msg - @opsetup def setup(self, rwork=None, iwork=None): # select discretization of the advected fields @@ -164,7 +74,7 @@ class AdvectionDir(Computational): for v in self.variables if v is not self.velocity] # GPU advection ... - if self.method[Support].find('gpu') >= 0: + if self._gpu_advection: topo_shape = advected_discrete_fields[0].topology.shape if topo_shape[self.direction] == 1: from hysop.gpu.gpu_particle_advection \ @@ -187,28 +97,19 @@ class AdvectionDir(Computational): self._is_uptodate = True def get_work_properties(self): - """ - Work vector for advection in one dir : - - [ interp , part_positions, fields_on_particles] - interp part is used also for remesh and time-integrator. - """ - assert self._is_discretized - + super(AdvectionDir, self).get_work_properties() # Shape of reference comes from fields, not from velocity advected_discrete_fields = [self.discreteFields[v] for v in self.variables if v is not self.velocity] topo = advected_discrete_fields[0].topology - # number of components of the rhs (time integration) - rhs_size = 1 - - if self.method[Support].find('gpu') < 0: + # Find number and shape of required work arrays + if not self._gpu_advection: # -- pure python advection -- # work array shape depends on the time integrator # interpolation scheme and remeshing scheme ti_work = self.method[TimeIntegrator].get_work_properties( - rhs_size, topo) + self._rhs_size, topo) ti_rwork_length = len(ti_work['rwork']) iw_prop = self.method[Interpolation].get_work_properties(topo) rw_prop = Remeshing.get_work_properties(topo) @@ -224,23 +125,31 @@ class AdvectionDir(Computational): # no work array iwork_length, rwork_length = 0, 0 - memshape = advected_discrete_fields[0].topology.mesh.resolution + # buffers for fields on particles rwork_length += np.sum([f.nb_components for f in self.advected_fields]) - if self.method[Support].find('gpu') < 0 or \ + if not self._gpu_advection or \ self.method[Support].find('gpu_2k') >= 0: - rwork_length += 1 # work array for position + rwork_length += 1 # work array for positions memsize = np.prod(topo.mesh.resolution) return {'rwork': [(memsize,)] * rwork_length, 'iwork': [(memsize,)] * iwork_length} @debug @opapply - def apply(self, simulation=None, dtCoeff=1.0, split_id=0, old_dir=0): + def apply(self, simulation=None, dt_coeff=1.0, split_id=0, old_dir=0): """ - Apply this operator to its variables. - @param simulation : object that describes the simulation - parameters (time, time step, iteration number ...), see - hysop.problem.simulation.Simulation for details. + + Parameters + ---------- + + simulation : `:class::~hysop.problem.simulation.Simulation` + dt_coeff : double + split_id : int, optional + old_dir : int, optional """ + if not self._single_topo and not self._gpu_advection: + raise ValueError("Multiscale advection is not yet supported " + "in pure Python, use Scales or GPU.") + self.discrete_op.apply(simulation, - dtCoeff, split_id, old_dir) + dt_coeff, split_id, old_dir) diff --git a/hysop/operator/computational.py b/hysop/operator/computational.py index 186e27489824717b9d3cb027ccf950c3b19a0675..1f55c96095b9cdeca8df2314367aa560426b6d0a 100755 --- a/hysop/operator/computational.py +++ b/hysop/operator/computational.py @@ -199,7 +199,7 @@ class Computational(Operator): return build_topos - def _standard_discretize(self, min_ghosts=0): + def _standard_discretize(self, min_ghosts=0, cutdir=None): """ This functions provides a standard way to discretize the operator, but some operators may need a specific discretization process. @@ -211,21 +211,22 @@ class Computational(Operator): if self._single_topo: # One topo for all fields ... if build_topos: - topo = self._build_topo(self._discretization, min_ghosts) + topo = self._build_topo(self._discretization, min_ghosts, + cutdir) for v in self.variables: self.variables[v] = topo else: - # Topo is already built, just check its ghosts + # Topo is already built, just check its ghosts and cutdir msg = 'The proposed ghost layer is not large enough.' ghosts = self.variables.values()[0].mesh.discretization.ghosts assert (ghosts >= min_ghosts).all(), msg - + assert (self.variables.values()[0].cutdir == cutdir).all() else: # ... or one topo for each field. for v in self.variables: if build_topos[v]: self.variables[v] = self._build_topo(self.variables[v], - min_ghosts) + min_ghosts, cutdir) build_topos[v] = False else: assert (self.variables[v].ghosts() >= min_ghosts).all() @@ -233,16 +234,16 @@ class Computational(Operator): # All topos are built, we can discretize fields. self._discretize_vars() - def _build_topo(self, discretization, min_ghosts): - """Build an mpi topology and its mesh from a given + def _build_topo(self, discretization, min_ghosts, cutdir=None): + """Build a mpi topology and its mesh from a given discretization. - """ # Reset ghosts if necessary ghosts = discretization.ghosts ghosts[ghosts < min_ghosts] = min_ghosts # build a topology from the given discretization - return self.domain.create_topology(discretization) + return self.domain.create_topology(discretization=discretization, + cutdir=cutdir) def _fftw_discretize(self): """ diff --git a/hysop/operator/discrete/discrete.py b/hysop/operator/discrete/discrete.py index 94e5418a36211662fa301e44a68349ff8d862da7..966d78865b014956c87a1863868911407f2c6224 100755 --- a/hysop/operator/discrete/discrete.py +++ b/hysop/operator/discrete/discrete.py @@ -20,8 +20,7 @@ class DiscreteOperator(object): @debug @abstractmethod - def __init__(self, variables, rwork=None, iwork=None, method=None, - mpi_params=None): + def __init__(self, variables, rwork=None, iwork=None, method=None): """ Parameters ----------- @@ -34,9 +33,6 @@ class DiscreteOperator(object): method : dictionnary, optional internal solver parameters (discretisation ...). If None, use default method from hysop.default_method.py. - mpi_params : :class:`hysop.tools.parameters.MPIParams - parameters to set mpi context. See Notes below. - Attributes ---------- @@ -46,13 +42,6 @@ class DiscreteOperator(object): input : fields used as input (i.e. read-only) output : fields used as in/out (i.e. modified during apply call) - Notes - ----- - - Methods : to be done ... - - MPIParams : to be done ... - """ if isinstance(variables, list): # variables @@ -166,8 +155,7 @@ class DiscreteOperator(object): def get_extra_args_from_method(op, key, default_value): - """ - Returns the given extra arguments dictionary from method attribute. + """Returns the given extra arguments dictionary from method attributes. Parameters ----------- @@ -179,6 +167,16 @@ def get_extra_args_from_method(op, key, default_value): default value when ExtraArgs is not in op.method or key is not in op.method[ExtraArgs] + Usage + ----- + + .. code:: + + method = {ExtraArgs: {'device_id': 2, user_src: ['./ker.cl']} + op = SomeOp(..., method=method) + val = get_extra_args_from_method(op, device_id, 6) + # set val to 2. If device_id or ExtraArgs does not exist, set val to 8. + """ from hysop.methods_keys import ExtraArgs try: diff --git a/hysop/operator/discrete/particle_advection.py b/hysop/operator/discrete/particle_advection.py index 9d81aa17de957ea7e1b0cc5752f318806cb1f2e7..a8fb2be0945138cc3073879311e48c9aa5631266 100644 --- a/hysop/operator/discrete/particle_advection.py +++ b/hysop/operator/discrete/particle_advection.py @@ -1,10 +1,10 @@ """Advection solver, particular method, pure-python version """ -from hysop.constants import debug, WITH_GUESS, HYSOP_REAL, HYSOP_INTEGER -from hysop.methods_keys import TimeIntegrator, Interpolation, Remesh, Support +from hysop.constants import debug, WITH_GUESS, HYSOP_REAL, HYSOP_INTEGER, EPS +from hysop.methods_keys import TimeIntegrator, Interpolation, Remesh from hysop.operator.discrete.discrete import DiscreteOperator -import hysop.default_methods as default +from hysop.default_methods import ADVECTION import numpy as np from hysop.numerics.remeshing import Remeshing from hysop.tools.profiler import profile @@ -47,7 +47,7 @@ class ParticleAdvection(DiscreteOperator): variables += self.fields_on_grid if 'method' not in kwds: - kwds['method'] = default.ADVECTION + kwds['method'] = ADVECTION # number of components of the RHS (see time integrator) self._rhs_size = 1 # dictionnary of internal buffers used @@ -80,9 +80,9 @@ class ParticleAdvection(DiscreteOperator): # on grid point to add a particle. # If None, particles everywhere. self._threshold = None - #self._init_particles_on_grid = self._init_particles_everywhere - #if self._threshold is not None: - # self._init_particles_on_grid = self._init_particles_with_threshold + self._init_particles_on_grid = self._init_particles_everywhere + if self._threshold is not None: + self._init_particles_on_grid = self._init_particles_with_threshold def _configure_numerical_methods(self): """Function to set the numerical method for python operator and link them @@ -116,7 +116,7 @@ class ParticleAdvection(DiscreteOperator): self._rhs_size, rwork=self._rw_integ, f=num_interpolate, topo=topo_fields, - optim=WITH_GUESS, + #optim=WITH_GUESS, indices=ic) # -- remesh -- @@ -128,139 +128,113 @@ class ParticleAdvection(DiscreteOperator): def _set_work_arrays(self, rwork=None, iwork=None): - if self.method[Support].find('gpu') < 0: - # Find number and shape of required work arrays - - # topologies for advected field(s) grid - topo_fields = self.fields_on_grid[0].topology - # -- work arrays properties for interpolation -- - # work arrays shape depends on the targeted grid, i.e. topo_fields. - interp_wp = self.method[Interpolation].get_work_properties(topo_fields) - interp_iwork_length = len(interp_wp['iwork']) - interp_rwork_length = len(interp_wp['rwork']) - - # -- work arrays properties for time-integrator -- - # depends on topo_fields. - ti_wp = self.method[TimeIntegrator].get_work_properties( - self._rhs_size, topo_fields) - ti_rwork_length = len(ti_wp['rwork']) - # -- work arrays properties for remeshing scheme -- - # depends on topo_fields. - remesh_wp = Remeshing.get_work_properties(topo_fields) - remesh_iwork_length = len(remesh_wp['iwork']) - remesh_rwork_length = len(remesh_wp['rwork']) - - # -- Find out the number of required work arrays -- - # Interpolation and odesolver work arrays must be different - # but after interpolation/advection, work arrays can be reused - # by remeshing. So the number of required work arrays is - # equal to - # max(nb works_for_time_integrator + nb work_for_interpolation, - # nb works_for_remesh) - rwork_length = max(remesh_rwork_length, - ti_rwork_length + interp_rwork_length) - iwork_length = max(interp_iwork_length, remesh_iwork_length) - - # -- Finally, particles positions and fields on particles - # values will also be in work array. - ppos_index = rwork_length - # max number of components for advected fields - nbc = max([f.nb_components for f in self.fields_on_grid]) - rwork_length += nbc - rwork_length += 1 # work array for positions - - # -- rwork organisation -- - - # rw = [ | rw_interp | rw_time_int | part_pos | fields_on_part |] - # | rw_remesh | ... |] - - # work array shape depends on the time integrator - # interpolation scheme and remeshing scheme - - # check and/or allocate work arrays according to properties above - # reshape stuff is done inside numerical methods. - # Here we just allocate flat memory. - subshape_field = tuple(topo_fields.mesh.resolution) - subsize_interp_ti = interp_wp['rwork'] + ti_wp['rwork'] - subsize_remesh = remesh_wp['rwork'] - # list of shapes for rwork arrays (!flat arrays!) - subsize_rw = max(subsize_interp_ti, subsize_remesh) - # list of shapes for iwork arrays (!flat arrays!) - subsize_iw = max(interp_wp['iwork'], remesh_wp['iwork']) - # one array for particles positions - # and arrays for fields on particles - for i in xrange(nbc + 1): - subsize_rw.append((np.prod(subshape_field), )) - - # Build iwork and rwork lists - self._rwork = WorkSpaceTools.check_work_array(rwork_length, - subsize_rw, - rwork, HYSOP_REAL) - self._iwork = WorkSpaceTools.check_work_array(iwork_length, - subsize_iw, - iwork, HYSOP_INTEGER) - - # -- set connections for remesh, ti and interp methods -- - # interpolation - self._rw_interp = [self._rwork[i] - for i in xrange(interp_rwork_length)] - self._iw_interp = [self._iwork[i] - for i in xrange(interp_iwork_length)] - assert np.asarray([npw.arrays_share_data(self._rw_interp[i], - self._rwork[i]) - for i in xrange(interp_rwork_length)]).all() - - # time integration - self._rw_integ = [self._rwork[i] - for i in xrange(interp_rwork_length, - interp_rwork_length + - ti_rwork_length)] - # remesh - self._rw_remesh = [self._rwork[i] - for i in xrange(remesh_rwork_length)] - self._iw_remesh = [self._iwork[i] - for i in xrange(remesh_iwork_length)] - - # --- links to work for particles positions buffers --- - self._work_part_pos = ppos_index - self.part_positions = WorkSpaceTools.check_work_array( - 1, subshape_field, [self._rwork[ppos_index]]) - assert npw.arrays_share_data(self.part_positions[0], - self._rwork[ppos_index]) - - # --- links to buffers for fields on particles --- - ppos_index += 1 - self.fields_on_part = WorkSpaceTools.check_work_array( - nbc, subshape_field, - [self._rwork[i] for i in xrange(ppos_index, ppos_index + nbc)]) - assert np.asarray( - [npw.arrays_share_data(self.fields_on_part[i], - self._rwork[ppos_index + i]) - for i in xrange(nbc)]).all() + # Find number and shape of required work arrays - else: - # --- GPU setup --- - topo = self.fields_on_grid[0].topology - # Find number and shape of required work arrays - # For GPU version, no need of numerics works - # Shape of reference comes from fields, not from velocity - rwork_length = np.sum([f.nb_components - for f in self.fields_on_grid]) - if self.method[Support].find('gpu_2k') >= 0: - rwork_length += 1 # work array for positions - - # check and/or allocate work arrays according to properties above - subshape = tuple(topo.mesh.resolution) - self._rwork = WorkSpaceTools.check_work_array(rwork_length, - subshape, - rwork, HYSOP_REAL) - # seems that iwork is not required for GPU, this must be checked - self._iwork = WorkSpaceTools.check_work_array(0, subshape, - iwork, HYSOP_INTEGER) + # topologies for advected field(s) grid + topo_fields = self.fields_on_grid[0].topology + # -- work arrays properties for interpolation -- + # work arrays shape depends on the targeted grid, i.e. topo_fields. + interp_wp = self.method[Interpolation].get_work_properties(topo_fields) + interp_iwork_length = len(interp_wp['iwork']) + interp_rwork_length = len(interp_wp['rwork']) + + # -- work arrays properties for time-integrator -- + # depends on topo_fields. + ti_wp = self.method[TimeIntegrator].get_work_properties( + self._rhs_size, topo_fields) + ti_rwork_length = len(ti_wp['rwork']) + # -- work arrays properties for remeshing scheme -- + # depends on topo_fields. + remesh_wp = Remeshing.get_work_properties(topo_fields) + remesh_iwork_length = len(remesh_wp['iwork']) + remesh_rwork_length = len(remesh_wp['rwork']) + + # -- Find out the number of required work arrays -- + # Interpolation and odesolver work arrays must be different + # but after interpolation/advection, work arrays can be reused + # by remeshing. So the number of required work arrays is + # equal to + # max(nb works_for_time_integrator + nb work_for_interpolation, + # nb works_for_remesh) + rwork_length = max(remesh_rwork_length, + ti_rwork_length + interp_rwork_length) + iwork_length = max(interp_iwork_length, remesh_iwork_length) + + # -- Finally, particles positions and fields on particles + # values will also be in work array. + ppos_index = rwork_length + # max number of components for advected fields + nbc = max([f.nb_components for f in self.fields_on_grid]) + rwork_length += nbc + rwork_length += 1 # work array for positions + + # -- rwork organisation -- + + # rw = [ | rw_interp | rw_time_int | part_pos | fields_on_part |] + # | rw_remesh | ... |] + + # work array shape depends on the time integrator + # interpolation scheme and remeshing scheme + + # check and/or allocate work arrays according to properties above + # reshape stuff is done inside numerical methods. + # Here we just allocate flat memory. + subshape_field = tuple(topo_fields.mesh.compute_resolution) + subsize_interp_ti = interp_wp['rwork'] + ti_wp['rwork'] + subsize_remesh = remesh_wp['rwork'] + # list of shapes for rwork arrays (!flat arrays!) + subsize_rw = max(subsize_interp_ti, subsize_remesh) + # list of shapes for iwork arrays (!flat arrays!) + subsize_iw = max(interp_wp['iwork'], remesh_wp['iwork']) + # one array for particles positions + # and arrays for fields on particles + for i in xrange(nbc + 1): + subsize_rw.append((np.prod(subshape_field), )) + + # Build iwork and rwork lists + self._rwork = WorkSpaceTools.check_work_array(rwork_length, subsize_rw, + rwork, HYSOP_REAL) + self._iwork = WorkSpaceTools.check_work_array(iwork_length, subsize_iw, + iwork, HYSOP_INTEGER) + + # -- set connections for remesh, ti and interp methods -- + # interpolation + self._rw_interp = [self._rwork[i] for i in xrange(interp_rwork_length)] + self._iw_interp = [self._iwork[i] for i in xrange(interp_iwork_length)] + assert np.asarray([npw.arrays_share_data(self._rw_interp[i], + self._rwork[i]) + for i in xrange(interp_rwork_length)]).all() + + # time integration + self._rw_integ = [self._rwork[i] + for i in xrange(interp_rwork_length, + interp_rwork_length + + ti_rwork_length)] + # remesh + self._rw_remesh = [self._rwork[i] for i in xrange(remesh_rwork_length)] + self._iw_remesh = [self._iwork[i] for i in xrange(remesh_iwork_length)] + + # --- links to work for particles positions buffers --- + self._work_part_pos = ppos_index + self.part_positions = WorkSpaceTools.check_work_array( + 1, subshape_field, [self._rwork[ppos_index]]) + assert npw.arrays_share_data(self.part_positions[0], + self._rwork[ppos_index]) + + # --- links to buffers for fields on particles --- + ppos_index += 1 + self.fields_on_part = WorkSpaceTools.check_work_array( + nbc, subshape_field, + [self._rwork[i] for i in xrange(ppos_index, ppos_index + nbc)]) + assert np.asarray( + [npw.arrays_share_data(self.fields_on_part[i], + self._rwork[ppos_index + i]) + for i in xrange(nbc)]).all() @debug @profile - def apply(self, simulation=None, dt_coeff=1., split_id=0, old_dir=0): + def apply(self, simulation=None): + #, dt_coeff=1., split_id=0, old_dir=0): """ Advection algorithm: - initialize particles and fields with their values on the grid. @@ -272,12 +246,10 @@ class ParticleAdvection(DiscreteOperator): """ assert simulation is not None, \ 'Simulation parameter is missing.' + t, dt = simulation.time, simulation.sub_step - t, dt = simulation.time, simulation.time_step * dt_coeff - - # Initialize particles on the grid - toporef = self.fields_on_grid[0].topology - self.part_positions[0][...] = toporef.mesh.coords[self.direction] + # -- Initialize particles on the grid -- + self._init_particles_on_grid() # -- Advect particles -- # RHS of odesolver must be computed and save in odesolver._rwork[0] @@ -285,7 +257,7 @@ class ParticleAdvection(DiscreteOperator): # - ghost-synchro of the rhs (if needed) is done by the odesolver # - interpolation of velocity field from grid to particles is done # 'inside' the ode solver (rhs = interpolation scheme) - self.num_advec.rwork[0][...] = self.velocity.data[self.direction][...] + #self.num_advec.rwork[0][...] = self.velocity.data[self.direction][...] self.part_positions = self.num_advec( t, self.part_positions, dt, result=self.part_positions) @@ -297,3 +269,37 @@ class ParticleAdvection(DiscreteOperator): # Remesh all components of the field fg = self.num_remesh(ppos=self.part_positions, pscal=self.fields_on_part, result=fg) + + def _init_particles_everywhere(self): + """Initialize particles on grid points + """ + # indices of computation points (i.e. no ghosts) + topo = self.fields_on_grid[0].topology + ic = topo.mesh.compute_index + self.part_positions[0][...] = \ + topo.mesh.coords[self.direction][ic[self.direction]] + 10 * EPS + # remark : we shift particles with EPS to avoid + # the case where, because of floating point precision + # xp (particle position) < xg (grid position), + # with xp - xg very small. If this happens, + # left point of xp will be xg - 1 rather than xg + # which is not what is expected. + + def _init_particles_with_threshold(self): + """"Initialize particles on grid points where + advected field value is greater than a given threshold + """ + raise AttributeError("Not yet implemented") + # first field is used as reference to check threshold + ref_field = self.fields_on_grid[0].data + topo = self.fields_on_grid[0].topology + ic = topo.mesh.compute_index + self._rwork[0][...] = 0. + for i in xrange(len(ref_field)): + self._rwork[0][...] += ref_field[i][ic] ** 2 + self._rwork[0][...] = np.sqrt(self._rwork[0]) + ipart = np.where(self._rwork < self._threshold) + part_shape = ref_field[0][ipart].shape + self.part_positions = WorkSpaceTools.check_work_array( + 1, part_shape, self._rwork[self._work_part_pos]) + self.part_positions[0][...] = topo.mesh.coords[self.direction][ipart] diff --git a/hysop/operator/tests/test_advec_scales.py b/hysop/operator/tests/test_advec_scales.py index eabb4b92519b3496f2d831d64508c56261b50255..1fa5c36efc284e900234e229cae3560cdd0fbe2e 100755 --- a/hysop/operator/tests/test_advec_scales.py +++ b/hysop/operator/tests/test_advec_scales.py @@ -1,484 +1,219 @@ """Testing Scales advection operator. """ +from __future__ import print_function import numpy as np -from hysop.methods_keys import Scales, TimeIntegrator, Interpolation,\ +from hysop.methods_keys import TimeIntegrator, Interpolation,\ Remesh, Support, Splitting from hysop.methods import RK2, L2_1, L4_2, M8Prime, Linear from hysop.domain.box import Box from hysop.fields.continuous import Field -from hysop.operator.advection import Advection +from hysop.operator.advection import Advection, ScalesAdvection +from hysop.problem.simulation import O2FULLHALF from hysop.problem.simulation import Simulation -import hysop.tools.numpywrappers as npw from hysop.testsenv import scales_failed from hysop.tools.parameters import Discretization d3d = Discretization([17, 17, 17]) d3d_g = Discretization([17, 17, 17], [2, 2, 2]) - -@scales_failed -def test_nullVelocity_m4(): - """Basic test with random velocity. Using M4prime""" +m1 = {TimeIntegrator: RK2, Interpolation: Linear, Remesh: L2_1, + Splitting: O2FULLHALF, Support: ''} +m2 = {TimeIntegrator: RK2, Interpolation: Linear, Remesh: L4_2, + Splitting: O2FULLHALF, Support: ''} +m3 = {TimeIntegrator: RK2, Interpolation: Linear, Remesh: M8Prime, + Splitting: O2FULLHALF, Support: ''} + + +def run_advection(vector_field, method, mref=None, random_velocity=False): + """Create advection operator, ref operator + and fields, run scales and python advection, + compare results. + + Parameters + ---------- + vector_field: bool + True to advect a vector field, else scalar field. + method : dictionnary + Set scales remeshing type. + mref : dictionnary, optional + method for pure python advection. Default = m1 + random_velocity: boolean + If True, randomize velocity values, else set it to 0. + """ box = Box(length=[1., 1., 1.], origin=[0., 0., 0.]) - scal = Field(domain=box, name='Scalar') - scal_ref = Field(domain=box, name='Scalar_ref') + scal = Field(domain=box, name='Scalar', is_vector=vector_field) + scal_ref = Field(domain=box, name='Scalar_ref', is_vector=vector_field) velo = Field(domain=box, name='Velocity', formula=lambda x, y, z, t: (0., 0., 0.), is_vector=True, vectorize_formula=True) - advec = Advection(velo, scal, discretization=d3d, method={Scales: 'p_M4'}) - advec_py = Advection(velo, scal_ref, discretization=d3d_g, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Splitting: 'o2_FullHalf', - Support: ''}) + advec = ScalesAdvection(velocity=velo, + advected_fields=scal, + discretization=d3d, + method=method) + if mref is None: + mref = m1 + advec_py = Advection(velocity=velo, + advected_fields=scal_ref, discretization=d3d_g, + method=mref) advec.discretize() advec_py.discretize() advec.setup() advec_py.setup() - scal_d = scal.discreteFields.values()[0] - scal_ref_d = scal_ref.discreteFields.values()[0] - topo = scal_d.topology - toporef = scal_ref_d.topology + # Get and randomize discrete fields + topo = advec.advected_fields_topology() + topo_ref = advec_py.advected_fields_topology() + + if random_velocity: + topo_velo = advec.velocity_topology() + vd = velo.randomize(topo_velo) + for d in xrange(velo.nb_components): + vd[d] /= 2 * topo_velo.mesh.resolution[d] + else: + assert (velo.norm(topo) == 0).all() ic = topo.mesh.compute_index - icref = toporef.mesh.compute_index - scal_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) - scal_ref_d.data[0][ic] = scal_d.data[0][icref] - assert (velo.norm(topo) == 0).all() + icref = topo_ref.mesh.compute_index + print(icref) + print(ic) + scal_d = scal.randomize(topo) + scal_ref_d = scal.discretize(topo_ref) + for d in xrange(len(scal_d.data)): + scal_ref_d[d][icref] = scal_d[d][ic] + assert np.allclose(scal_ref_d.data[d][icref], scal_d.data[d][ic]) advec.apply(Simulation(start=0., end=0.1, nb_iter=1)) advec_py.apply(Simulation(start=0., end=0.1, nb_iter=1)) - print (np.max(np.abs(scal_ref_d.data[0][icref] - scal_d.data[0][ic]))) - assert np.allclose(scal_ref_d.data[0][icref], scal_d.data[0][ic]) - - -@scales_failed -def test_nullVelocity_vec_m4(): - """Basic test with random velocity and vector field. Using M4prime""" + # compare values on grid (excluding ghosts) + for d in xrange(len(scal_d.data)): + assert np.allclose(scal_ref_d.data[d][icref], scal_d.data[d][ic]) + + +def run_advection_2(vector_field, method): + """Create advection operator, + fields, run and check results. + + Parameters + ---------- + vector_field: bool + True to advect a vector field, else scalar field. + method : dictionnary + Set scales remeshing type. + """ box = Box(length=[1., 1., 1.], origin=[0., 0., 0.]) - scal = Field(domain=box, name='Scalar', is_vector=True) - scal_ref = Field(domain=box, name='Scalar_ref', is_vector=True) + scal = Field(domain=box, name='Scalar', is_vector=vector_field) + scal_ref = Field(domain=box, name='Scalar_ref', is_vector=vector_field) velo = Field(domain=box, name='Velocity', formula=lambda x, y, z, t: (0., 0., 0.), is_vector=True, vectorize_formula=True) - advec = Advection(velo, scal, discretization=d3d, method={Scales: 'p_M4'}) - advec_py = Advection(velo, scal_ref, discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Splitting: 'o2_FullHalf', - Support: ''} - ) + advec = ScalesAdvection(velocity=velo, + advected_fields=scal, discretization=d3d, + method=method) advec.discretize() - advec_py.discretize() advec.setup() - advec_py.setup() - - scal_d = scal.discreteFields.values()[0] - scal_ref_d = scal_ref.discreteFields.values()[0] - - scal_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) - scal_d.data[1][...] = npw.asrealarray( - np.random.random(scal_d.data[1].shape)) - scal_d.data[2][...] = npw.asrealarray( - np.random.random(scal_d.data[2].shape)) - scal_ref_d.data[0][...] = scal_d.data[0][...] - scal_ref_d.data[1][...] = scal_d.data[1][...] - scal_ref_d.data[2][...] = scal_d.data[2][...] - + # Get and randomize discrete fields + topo = advec.advected_fields_topology() + assert (velo.norm(topo) == 0).all() + ic = topo.mesh.compute_index + scal_d = scal.randomize(topo) + scal_ref.copy(scal, topo) + scal_ref_d = scal.discretize(topo) + for d in xrange(len(scal_d.data)): + assert np.allclose(scal_ref_d.data[d][ic], scal_d.data[d][ic]) advec.apply(Simulation(start=0., end=0.1, nb_iter=1)) - advec_py.apply(Simulation(start=0., end=0.1, nb_iter=1)) - - assert np.allclose(scal_ref_d.data[0], scal_d.data[0]) - assert np.allclose(scal_ref_d.data[1], scal_d.data[1]) - assert np.allclose(scal_ref_d.data[2], scal_d.data[2]) + # compare values on grid (excluding ghosts) + for d in xrange(len(scal_d.data)): + assert np.allclose(scal_ref_d.data[d][ic], scal_d.data[d][ic]) @scales_failed -def test_nullVelocity_m6(): - """Basic test with null velocity. Using M6prime""" - box = Box(length=[1., 1., 1.], origin=[0., 0., 0.]) - scal = Field(domain=box, name='Scalar') - velo = Field(domain=box, name='Velocity', - formula=lambda x, y, z, t: (0., 0., 0.), is_vector=True, - vectorize_formula=True) - advec = Advection(velo, scal, discretization=d3d, method={Scales: 'p_M6'} - ) - advec.discretize() - advec.setup() - - scal_d = scal.discreteFields.values()[0] - scal_d.data[0][...] = np.asarray(np.random.random(scal_d.data[0].shape)) - scal_init = npw.copy(scal_d.data[0]) - - advec.apply(Simulation(start=0., end=0.1, nb_iter=1)) - - assert np.allclose(scal_init, scal_d.data[0]) +def test_null_velocity_m4(): + """Scalar field advection, null velocity, M4prime remesh + """ + method = {Remesh: 'p_M4'} + run_advection_2(False, method) @scales_failed -def test_nullVelocity_vec_m6(): - """Basic test with null velocity and vector field. Using M6prime""" - box = Box(length=[1., 1., 1.], origin=[0., 0., 0.]) - scal = Field(domain=box, name='Scalar', is_vector=True) - velo = Field(domain=box, name='Velocity', - formula=lambda x, y, z, t: (0., 0., 0.), is_vector=True, - vectorize_formula=True) - advec = Advection(velo, scal, discretization=d3d, method={Scales: 'p_M6'} - ) - advec.discretize() - advec.setup() - - scal_d = scal.discreteFields.values()[0] - scal_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) - scal_d.data[1][...] = npw.asrealarray( - np.random.random(scal_d.data[1].shape)) - scal_d.data[2][...] = npw.asrealarray( - np.random.random(scal_d.data[2].shape)) - scal_init0 = npw.copy(scal_d.data[0]) - scal_init1 = npw.copy(scal_d.data[1]) - scal_init2 = npw.copy(scal_d.data[2]) - - advec.apply(Simulation(start=0., end=0.1, nb_iter=1)) - - assert np.allclose(scal_init0, scal_d.data[0]) - assert np.allclose(scal_init1, scal_d.data[1]) - assert np.allclose(scal_init2, scal_d.data[2]) +def test_null_velocity_vec_m4(): + """Vector field advection, null velocity, M4prime remesh + """ + method = {Remesh: 'p_M4'} + run_advection_2(True, method) @scales_failed -def test_nullVelocity_m8(): - """Basic test with random velocity. Using M4prime""" - box = Box(length=[1., 1., 1.], origin=[0., 0., 0.]) - scal = Field(domain=box, name='Scalar') - scal_ref = Field(domain=box, name='Scalar_ref') - velo = Field(domain=box, name='Velocity', - formula=lambda x, y, z, t: (0., 0., 0.), is_vector=True, - vectorize_formula=True) - advec = Advection(velo, scal, discretization=d3d, method={Scales: 'p_M8'} - ) - advec_py = Advection(velo, scal_ref, discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Splitting: 'o2_FullHalf', - Support: ''} - ) - advec.discretize() - advec_py.discretize() - advec.setup() - advec_py.setup() - topo = advec.discreteFields[velo].topology - scal_d = scal.discreteFields[topo] - scal_ref_d = scal_ref.discreteFields[topo] - - scal_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) - scal_ref_d.data[0][...] = scal_d.data[0] - - advec.apply(Simulation(start=0., end=0.1, nb_iter=1)) - advec_py.apply(Simulation(start=0., end=0.1, nb_iter=1)) - assert np.allclose(scal_ref_d.data[0], scal_d.data[0]) +def test_null_velocity_m6(): + """Scalar field advection, null velocity, M6prime remesh + """ + method = {Remesh: 'p_M6'} + run_advection_2(False, method) @scales_failed -def test_nullVelocity_vec_m8(): - """Basic test with random velocity and vector field. Using M4prime""" - box = Box(length=[1., 1., 1.], origin=[0., 0., 0.]) - scal = Field(domain=box, name='Scalar', is_vector=True) - scal_ref = Field(domain=box, name='Scalar_ref', is_vector=True) - velo = Field(domain=box, name='Velocity', - formula=lambda x, y, z, t: (0., 0., 0.), is_vector=True, - vectorize_formula=True) - advec = Advection(velo, scal, discretization=d3d, method={Scales: 'p_M8'} - ) - advec_py = Advection(velo, scal_ref, discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Splitting: 'o2_FullHalf', - Support: ''} - ) - advec.discretize() - advec_py.discretize() - advec.setup() - advec_py.setup() - scal_d = scal.discreteFields.values()[0] - scal_ref_d = scal_ref.discreteFields.values()[0] - topo = scal_d.topology - assert (velo.norm(topo) == 0).all() - for i in xrange(box.dimension): - scal_d.data[i][...] = \ - npw.asrealarray(np.random.random(scal_d.data[i].shape)) - scal_ref_d.data[i][...] = scal_d.data[i][...] - - advec.apply(Simulation(start=0., end=0.075, nb_iter=1)) - advec_py.apply(Simulation(start=0., end=0.075, nb_iter=1)) - - for i in xrange(box.dimension): - assert np.allclose(scal_ref_d.data[i], scal_d.data[i]) +def test_null_velocity_vec_m6(): + """Vector field advection, null velocity, M6prime remesh + """ + method = {Remesh: 'p_M6'} + run_advection_2(True, method) @scales_failed -def _randomVelocity_m4(): - """Basic test with random velocity. Using M4prime""" - box = Box(length=[1., 1., 1.], origin=[0., 0., 0.]) - scal = Field(domain=box, name='Scalar') - scal_ref = Field(domain=box, name='Scalar_ref') - velo = Field(domain=box, name='Velocity', is_vector=True) - advec = Advection(velo, scal, discretization=d3d, method={Scales: 'p_M4'} - ) - advec_py = Advection(velo, scal_ref, discretization=d3d, method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Splitting: 'o2_FullHalf', - Support: ''} - ) - advec.discretize() - advec_py.discretize() - advec.setup() - advec_py.setup() - - scal_d = scal.discreteFields.values()[0] - scal_ref_d = scal_ref.discreteFields.values()[0] - velo_d = velo.discreteFields.values()[0] - - scal_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) - scal_ref_d.data[0][...] = scal_d.data[0][...] - velo_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) / (2. * scal_d.resolution[0]) - velo_d.data[1][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) / (2. * scal_d.resolution[1]) - velo_d.data[2][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) / (2. * scal_d.resolution[1]) - - advec.apply(Simulation(start=0., end=0.1, nb_iter=1)) - advec_py.apply(Simulation(start=0., end=0.1, nb_iter=1)) - - assert np.allclose(scal_ref_d.data[0], scal_d.data[0]) +def test_null_velocity_m8(): + """Scalar field advection, null velocity, M6prime remesh + """ + method = {Remesh: 'p_M8'} + run_advection_2(False, method) @scales_failed -def _randomVelocity_vec_m4(): - """Basic test with random velocity vector Field. Using M4prime""" - box = Box(length=[1., 1., 1.], origin=[0., 0., 0.]) - scal = Field(domain=box, name='Scalar', is_vector=True) - scal_ref = Field(domain=box, name='Scalar_ref', is_vector=True) - velo = Field(domain=box, name='Velocity', is_vector=True) - advec = Advection(velo, scal, discretization=d3d, method={Scales: 'p_M4'} - ) - advec_py = Advection(velo, scal_ref, discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L2_1, - Splitting: 'o2_FullHalf', - Support: ''} - ) - advec.discretize() - advec_py.discretize() - advec.setup() - advec_py.setup() - - scal_d = scal.discreteFields.values()[0] - scal_ref_d = scal_ref.discreteFields.values()[0] - velo_d = velo.discreteFields.values()[0] - - scal_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) - scal_d.data[1][...] = npw.asrealarray( - np.random.random(scal_d.data[1].shape)) - scal_d.data[2][...] = npw.asrealarray( - np.random.random(scal_d.data[2].shape)) - scal_ref_d.data[0][...] = scal_d.data[0][...] - scal_ref_d.data[1][...] = scal_d.data[1][...] - scal_ref_d.data[2][...] = scal_d.data[2][...] - velo_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) / (2. * scal_d.resolution[0]) - velo_d.data[1][...] = npw.asrealarray( - np.random.random(scal_d.data[1].shape)) / (2. * scal_d.resolution[1]) - velo_d.data[2][...] = npw.asrealarray( - np.random.random(scal_d.data[2].shape)) / (2. * scal_d.resolution[1]) - - advec.apply(Simulation(start=0., end=0.1, nb_iter=1)) - advec_py.apply(Simulation(start=0., end=0.1, nb_iter=1)) - - assert np.allclose(scal_ref_d.data[0], scal_d.data[0]) - assert np.allclose(scal_ref_d.data[1], scal_d.data[1]) - assert np.allclose(scal_ref_d.data[2], scal_d.data[2]) +def test_null_velocity_vec_m8(): + """Vector field advection, null velocity, M6prime remesh + """ + method = {Remesh: 'p_M8'} + run_advection_2(True, method) @scales_failed -def test_randomVelocity_m6(): - """Basic test with random velocity. Using M6prime""" - box = Box(length=[1., 1., 1.], origin=[0., 0., 0.]) - scal = Field(domain=box, name='Scalar') - scal_ref = Field(domain=box, name='Scalar_ref') - velo = Field(domain=box, name='Velocity', is_vector=True) - advec = Advection(velo, scal, discretization=d3d, method={Scales: 'p_M6'}) - advec_py = Advection(velo, scal_ref, discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Splitting: 'o2_FullHalf', - Support: ''} - ) - advec.discretize() - advec_py.discretize() - advec.setup() - advec_py.setup() - scal_d = scal.discreteFields.values()[0] - scal_ref_d = scal_ref.discreteFields.values()[0] - velo_d = velo.discreteFields.values()[0] - - scal_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) - scal_ref_d.data[0][...] = scal_d.data[0][...] - velo_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) / (2. * scal_d.resolution[0]) - velo_d.data[1][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) / (2. * scal_d.resolution[1]) - velo_d.data[2][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) / (2. * scal_d.resolution[1]) - advec.apply(Simulation(start=0., end=0.1, nb_iter=1)) - advec_py.apply(Simulation(start=0., end=0.1, nb_iter=1)) - assert np.allclose(scal_ref_d.data[0], scal_d.data[0]) +def test_random_velocity_m4(): + """Scalar field advection, random velocity, M4prime remesh + """ + method = {Remesh: 'p_M4'} + run_advection(False, method, True) @scales_failed -def test_randomVelocity_vec_m6(): - """Basic test with random velocity vector Field. Using M6prime""" - box = Box(length=[1., 1., 1.], origin=[0., 0., 0.]) - scal = Field(domain=box, name='Scalar', is_vector=True) - scal_ref = Field(domain=box, name='Scalar_ref', is_vector=True) - velo = Field(domain=box, name='Velocity', is_vector=True) - advec = Advection(velo, scal, discretization=d3d, method={Scales: 'p_M6'} - ) - advec_py = Advection(velo, scal_ref, discretization=d3d, method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: L4_2, - Splitting: 'o2_FullHalf', - Support: ''} - ) - advec.discretize() - advec_py.discretize() - advec.setup() - advec_py.setup() - - scal_d = scal.discreteFields.values()[0] - scal_ref_d = scal_ref.discreteFields.values()[0] - velo_d = velo.discreteFields.values()[0] - - scal_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) - scal_d.data[1][...] = npw.asrealarray( - np.random.random(scal_d.data[1].shape)) - scal_d.data[2][...] = npw.asrealarray( - np.random.random(scal_d.data[2].shape)) - scal_ref_d.data[0][...] = scal_d.data[0][...] - scal_ref_d.data[1][...] = scal_d.data[1][...] - scal_ref_d.data[2][...] = scal_d.data[2][...] - velo_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) / (2. * scal_d.resolution[0]) - velo_d.data[1][...] = npw.asrealarray( - np.random.random(scal_d.data[1].shape)) / (2. * scal_d.resolution[1]) - velo_d.data[2][...] = npw.asrealarray( - np.random.random(scal_d.data[2].shape)) / (2. * scal_d.resolution[1]) - - advec.apply(Simulation(start=0., end=0.1, nb_iter=1)) - advec_py.apply(Simulation(start=0., end=0.1, nb_iter=1)) - - assert np.allclose(scal_ref_d.data[0], scal_d.data[0]) - assert np.allclose(scal_ref_d.data[1], scal_d.data[1]) - assert np.allclose(scal_ref_d.data[2], scal_d.data[2]) +def test_random_velocity_vec_m4(): + """Vector field advection, random velocity, M4prime remesh + """ + method = {Remesh: 'p_M4'} + run_advection(True, method, True) @scales_failed -def test_randomVelocity_m8(): - """Basic test with random velocity. Using M8prime""" - box = Box(length=[1., 1., 1.], origin=[0., 0., 0.]) - scal = Field(domain=box, name='Scalar') - scal_ref = Field(domain=box, name='Scalar_ref') - velo = Field(domain=box, name='Velocity', is_vector=True) - advec = Advection(velo, scal, discretization=d3d, method={Scales: 'p_M8'} - ) - advec_py = Advection(velo, scal_ref, discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Splitting: 'o2_FullHalf', - Support: ''} - ) - advec.discretize() - advec_py.discretize() - advec.setup() - advec_py.setup() +def test_random_velocity_m6(): + """Scalar field advection, random velocity, M6prime remesh + """ + method = {Remesh: 'p_M6'} + run_advection(False, method, m2, True) - scal_d = scal.discreteFields.values()[0] - scal_ref_d = scal_ref.discreteFields.values()[0] - velo_d = velo.discreteFields.values()[0] - scal_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) - scal_ref_d.data[0][...] = scal_d.data[0][...] - velo_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) / (2. * scal_d.resolution[0]) - velo_d.data[1][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) / (2. * scal_d.resolution[1]) - velo_d.data[2][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) / (2. * scal_d.resolution[1]) +@scales_failed +def test_random_velocity_vec_m6(): + """Vector field advection, random velocity, M6prime remesh + """ + method = {Remesh: 'p_M6'} + run_advection(True, method, m2, True) - advec.apply(Simulation(start=0., end=0.1, nb_iter=1)) - advec_py.apply(Simulation(start=0., end=0.1, nb_iter=1)) - assert np.allclose(scal_ref_d.data[0], scal_d.data[0], atol=1e-07) +@scales_failed +def test_random_velocity_m8(): + """Scalar field advection, random velocity, M8prime remesh + """ + method = {Remesh: 'p_M8'} + run_advection(False, method, m3, True) @scales_failed -def test_randomVelocity_vec_m8(): - """Basic test with random velocity vector Field. Using M8prime""" - box = Box(length=[1., 1., 1.], origin=[0., 0., 0.]) - scal = Field(domain=box, name='Scalar', is_vector=True) - scal_ref = Field(domain=box, name='Scalar_ref', is_vector=True) - velo = Field(domain=box, name='Velocity', is_vector=True) - advec = Advection(velo, scal, discretization=d3d, method={Scales: 'p_M8'} - ) - advec_py = Advection(velo, scal_ref, discretization=d3d, - method={TimeIntegrator: RK2, - Interpolation: Linear, - Remesh: M8Prime, - Support: '', - Splitting: 'o2_FullHalf'} - ) - advec.discretize() - advec_py.discretize() - advec.setup() - advec_py.setup() +def test_random_velocity_vec_m8(): + """Vector field advection, random velocity, M8prime remesh + """ + method = {Remesh: 'p_M8'} + run_advection(True, method, m3, True) - scal_d = scal.discreteFields.values()[0] - scal_ref_d = scal_ref.discreteFields.values()[0] - velo_d = velo.discreteFields.values()[0] - - scal_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) - scal_d.data[1][...] = npw.asrealarray( - np.random.random(scal_d.data[1].shape)) - scal_d.data[2][...] = npw.asrealarray( - np.random.random(scal_d.data[2].shape)) - scal_ref_d.data[0][...] = scal_d.data[0][...] - scal_ref_d.data[1][...] = scal_d.data[1][...] - scal_ref_d.data[2][...] = scal_d.data[2][...] - velo_d.data[0][...] = npw.asrealarray( - np.random.random(scal_d.data[0].shape)) / (2. * scal_d.resolution[0]) - velo_d.data[1][...] = npw.asrealarray( - np.random.random(scal_d.data[1].shape)) / (2. * scal_d.resolution[1]) - velo_d.data[2][...] = npw.asrealarray( - np.random.random(scal_d.data[2].shape)) / (2. * scal_d.resolution[1]) - - advec.apply(Simulation(start=0., end=0.01, nb_iter=1)) - advec_py.apply(Simulation(start=0., end=0.01, nb_iter=1)) - - assert np.allclose(scal_ref_d.data[0], scal_d.data[0], atol=1e-07) - assert np.allclose(scal_ref_d.data[1], scal_d.data[1], atol=1e-07) - assert np.allclose(scal_ref_d.data[2], scal_d.data[2], atol=1e-07) diff --git a/hysop/problem/simulation.py b/hysop/problem/simulation.py index dea2196c18b2018fe72589ec9bb5811631353a7c..61f52e75ead5b72b10b143f2f0fbaa2f81fca778 100644 --- a/hysop/problem/simulation.py +++ b/hysop/problem/simulation.py @@ -218,3 +218,79 @@ class Simulation(object): self.current_iteration = -1 +class SplittingParameters(object): + """Configure splitting + + Time step is splitted into 'n' parts, + and for id in [0, n-1], self.splitting[id] = [dir, dt], where + dir is the concerned direction and dt the time step coefficient. + + n, dt, dir depend on the chosen scheme. + """ + __metaclass__ = ABCMeta + + def __init__(self, dimension): + self.splitting = self._conf(dimension) + + def __call__(self): + """Get splitting params""" + return self.splitting + + @abstractmethod + def _conf(self, dim=3): + """set parameters + """ + return None + + +class O2FULLHALF(SplittingParameters): + """timestep/2 in each directions""" + + def _conf(self, dimension=3): + params = [(i, 0.5) for i in xrange(dimension)] + for i in xrange(dimension): + params.append((dimension - 1 - i, 0.5)) + return params + + +class O1(SplittingParameters): + """ + """ + + def _conf(self, dimension=3): + return [((i, 1.)) for i in xrange(dimension)] + + +class O2(SplittingParameters): + """ + """ + + def _conf(self, dimension=3): + # Half timestep in all directions but last + splitting = [(i, 0.5) for i in xrange(dimension - 1)] + splitting.append((dimension - 1, 1.)) + for i in xrange(dimension - 1): + splitting.append((dimension - 2 - i, 0.5)) + return splitting + + +class x_only(SplittingParameters): + """ """ + + def _conf(self, dimension=3): + assert dimension > 0 + return [(0, 1.)] + +class y_only(SplittingParameters): + """ """ + + def _conf(self, dimension=3): + assert dimension > 1 + return [(1, 1.)] + +class z_only(SplittingParameters): + """ """ + + def _conf(self, dimension=3): + assert dimension > 2 + return [(2, 1.)]