From 48e5628d4a59ee2ecfed0f4315115a6ef6e27f77 Mon Sep 17 00:00:00 2001 From: Rao Garimella Date: Tue, 12 Jan 2021 10:51:54 -0700 Subject: [PATCH] Release 1.2.10 (#12) * Build system fixes * Fix warnings and errors * remove cinch; use googletest directly * fixes for building with LAPACKE; Also, use QR decomp if not LAPACK * Improve coding for least squares computation * Matrix inversion options and fixes * Initial work for non-cartesian coordinate systems * Efficiency improvements * Changes to swarm class (easier to construct and use) * MPI ghost update manager --- CMakeLists.txt | 29 +- cmake/FindTHRUST.cmake | 30 +- cmake/wontonConfig.cmake.in | 14 +- .../{axis_hpc.yml => axis_hpc_install.yml} | 6 +- jenkins/axis_hpc_nightly.yml | 21 + jenkins/axis_hpc_pr.yml | 4 +- jenkins/axis_varan.yml | 18 - jenkins/axis_varan_install.yml | 5 +- jenkins/axis_varan_nightly.yml | 21 + jenkins/axis_varan_pr.yml | 6 +- jenkins/build_docs.sh | 3 + jenkins/build_matrix_entry.sh | 134 ---- jenkins/build_matrix_entry_hpc.sh | 173 +++-- ...x_entry.sh => build_matrix_entry_varan.sh} | 149 ++-- wonton/CMakeLists.txt | 7 +- wonton/distributed/CMakeLists.txt | 57 ++ wonton/distributed/mpi_ghost_manager.h | 731 ++++++++++++++++++ .../test/test_mpi_ghost_manager.cc | 482 ++++++++++++ wonton/intersect/simple_intersect_for_tests.h | 117 +++ wonton/mesh/AuxMeshTopology.h | 73 +- wonton/mesh/CMakeLists.txt | 8 + wonton/mesh/simple/simple_mesh.h | 133 +++- wonton/mesh/simple/simple_mesh_wrapper.h | 11 +- wonton/mesh/test/my_mesh_wrapper.h | 76 ++ wonton/mesh/test/test_mesh_api.cc | 148 ++++ wonton/support/CoordinateSystem.h | 4 +- wonton/support/Matrix.cc | 14 +- wonton/support/Matrix.h | 32 +- wonton/support/Polytope.h | 94 +-- wonton/support/wonton.h | 11 +- wonton/swarm/swarm.h | 50 ++ wonton/swarm/swarm_state.h | 20 + wonton/swarm/test/test_swarm_state.cc | 6 +- 33 files changed, 2301 insertions(+), 386 deletions(-) rename jenkins/{axis_hpc.yml => axis_hpc_install.yml} (69%) create mode 100644 jenkins/axis_hpc_nightly.yml delete mode 100644 jenkins/axis_varan.yml create mode 100644 jenkins/axis_varan_nightly.yml delete mode 100755 jenkins/build_matrix_entry.sh rename jenkins/{install_matrix_entry.sh => build_matrix_entry_varan.sh} (52%) create mode 100644 wonton/distributed/CMakeLists.txt create mode 100644 wonton/distributed/mpi_ghost_manager.h create mode 100644 wonton/distributed/test/test_mpi_ghost_manager.cc create mode 100644 wonton/intersect/simple_intersect_for_tests.h create mode 100644 wonton/mesh/test/my_mesh_wrapper.h create mode 100644 wonton/mesh/test/test_mesh_api.cc diff --git a/CMakeLists.txt b/CMakeLists.txt index f4e539e..0770a51 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -33,7 +33,7 @@ endif () set(WONTON_VERSION_MAJOR 1) set(WONTON_VERSION_MINOR 2) -set(WONTON_VERSION_PATCH 2) +set(WONTON_VERSION_PATCH 8) # Top level target @@ -117,6 +117,12 @@ endif() #------------------------------------------------------------------------------# set(WONTON_ENABLE_Jali False CACHE BOOL "Jali Interface enabled?") + +# If MPI is enabled,allow the user to set Jali. If no MPI, throw an error. +if (WONTON_ENABLE_Jali AND NOT WONTON_ENABLE_MPI) + message(FATAL_ERROR "Jali cannot be enabled without mpi. Please set WONTON_ENABLE_MPI=True") +endif() + if (WONTON_ENABLE_Jali AND WONTON_ENABLE_MPI AND NOT Jali_LIBRARIES) # Look for the Jali package @@ -148,28 +154,11 @@ endif () set(WONTON_ENABLE_THRUST False CACHE BOOL "Is the Thrust library being used?") if (WONTON_ENABLE_THRUST) - - # Allow for swapping backends - set(THRUST_HOST_BACKEND "THRUST_HOST_SYSTEM_CPP" CACHE STRING "Thrust host backend") + set(THRUST_HOST_BACKEND "THRUST_HOST_SYSTEM_OMP" CACHE STRING "Thrust host backend") set(THRUST_DEVICE_BACKEND "THRUST_DEVICE_SYSTEM_OMP" CACHE STRING "Thrust device backend") - if ((${THRUST_HOST_BACKEND} STREQUAL "THRUST_HOST_SYSTEM_OMP") OR - (${THRUST_DEVICE_BACKEND} STREQUAL "THRUST_DEVICE_SYSTEM_OMP")) - list(APPEND _components OpenMP) - endif () - if (${THRUST_DEVICE_BACKEND} STREQUAL "THRUST_DEVICE_SYSTEM_CUDA") - list(APPEND _components CUDA) - endif () - - find_package(THRUST COMPONENTS ${_components} REQUIRED MODULE) - + find_package(THRUST COMPONENTS OpenMP REQUIRED MODULE) message(STATUS "Enabling compilation with Thrust") - message(STATUS "Using THRUST_ROOT=${THRUST_ROOT}") - message(STATUS "THRUST_INCLUDE_DIRS ${THRUST_INCLUDE_DIRS}") - - message(STATUS "Using ${THRUST_HOST_BACKEND} as Thrust host backend") - message(STATUS "Using ${THRUST_DEVICE_BACKEND} as Thrust device backend") - else () #----------------------------------------------------------------------------- diff --git a/cmake/FindTHRUST.cmake b/cmake/FindTHRUST.cmake index 05dcf55..b1541ac 100644 --- a/cmake/FindTHRUST.cmake +++ b/cmake/FindTHRUST.cmake @@ -11,11 +11,6 @@ # THRUST_COMPONENTS - Which backends were found (OpenMP, TBB, CUDA) #------------------------------------------------------------------------------# -# Create an INTERFACE library - -set(THRUST_LIBRARIES THRUST::THRUST) -add_library(${THRUST_LIBRARIES} INTERFACE IMPORTED) - #------------------------------------------------------------------------------# # Find the header file #------------------------------------------------------------------------------# @@ -29,7 +24,6 @@ else () endif () endif () -target_include_directories(${THRUST_LIBRARIES} INTERFACE ${THRUST_INCLUDE_DIRS}) #------------------------------------------------------------------------------# # Find the backend components @@ -38,25 +32,39 @@ target_include_directories(${THRUST_LIBRARIES} INTERFACE ${THRUST_INCLUDE_DIRS}) set(THRUST_COMPONENTS "") foreach (_component IN LISTS THRUST_FIND_COMPONENTS) if (${_component} STREQUAL "OpenMP") - find_package(OpenMP) + find_package(OpenMP COMPONENTS CXX) if (OpenMP_FOUND) - target_link_libraries(${THRUST_LIBRARIES} INTERFACE OpenMP::OpenMP_CXX) list(APPEND THRUST_COMPONENTS OpenMP) endif () elseif (${_component} STREQUAL "CUDA") find_package(CUDA) if (CUDA_FOUND) - target_link_libraries(${THRUST_LIBRARIES} INTERFACE "${CUDA_LIBRARIES}") list(APPEND THRUST_COMPONENTS CUDA) endif () - elseif () - message(FATAL_ERROR "Unknown component ${_component) requested for THRUST") + else () + message(FATAL_ERROR "Unknown component ${_component} requested for THRUST") endif () endforeach () message(STATUS "THRUST components found: ${THRUST_COMPONENTS}") +# Create an INTERFACE library + +set(THRUST_LIBRARIES THRUST::THRUST) +if (NOT TARGET ${THRUST_LIBRARIES}) + add_library(${THRUST_LIBRARIES} INTERFACE IMPORTED) + target_include_directories(${THRUST_LIBRARIES} INTERFACE ${THRUST_INCLUDE_DIRS}) + foreach (_component IN LISTS THRUST_COMPONENTS) + if (${_component} STREQUAL "OpenMP" AND OpenMP_FOUND) + target_link_libraries(${THRUST_LIBRARIES} INTERFACE OpenMP::OpenMP_CXX) + endif () + if (${_component} STREQUAL "CUDA" AND CUDA_FOUND) + target_link_libraries(${THRUST_LIBRARIES} INTERFACE "${CUDA_LIBRARIES}") + endif () + endforeach () +endif () + #------------------------------------------------------------------------------# # Set standard args stuff #------------------------------------------------------------------------------# diff --git a/cmake/wontonConfig.cmake.in b/cmake/wontonConfig.cmake.in index b97e907..65b3001 100644 --- a/cmake/wontonConfig.cmake.in +++ b/cmake/wontonConfig.cmake.in @@ -107,7 +107,19 @@ if (WONTON_ENABLE_LAPACKE) endif () if (WONTON_ENABLE_THRUST) - find_dependency(THRUST COMPONENTS ${THRUST_COMPONENTS}) + set(THRUST_LIBRARIES THRUST::THRUST) + + # First try discovery through a config file + find_package(THRUST NAMES Thrust CONFIG) + if (THRUST_FOUND) + string(REPLACE "THRUST_HOST_SYSTEM_" "" HOST_SYSTEM ${THRUST_HOST_BACKEND}) + string(REPLACE "THRUST_DEVICE_SYSTEM_" "" DEVICE_SYSTEM ${THRUST_DEVICE_BACKEND}) + if (NOT TARGET ${THRUST_LIBRARIES}) + thrust_create_target(${THRUST_LIBRARIES} HOST ${HOST_SYSTEM} DEVICE ${DEVICE_SYSTEM}) + endif () + else () + find_dependency(THRUST COMPONENTS ${THRUST_COMPONENTS} REQUIRED MODULE) + endif () endif () if (WONTON_ENABLE_Kokkos) diff --git a/jenkins/axis_hpc.yml b/jenkins/axis_hpc_install.yml similarity index 69% rename from jenkins/axis_hpc.yml rename to jenkins/axis_hpc_install.yml index 45113be..06897c5 100644 --- a/jenkins/axis_hpc.yml +++ b/jenkins/axis_hpc_install.yml @@ -3,7 +3,7 @@ COMPILER: - gcc6 - gcc7 -BUILD_TYPE: +CONFIG_TYPE: - base - debug - serial @@ -12,7 +12,7 @@ BUILD_TYPE: exclude: - COMPILER: gcc6 - BUILD_TYPE: readme + CONFIG_TYPE: readme - COMPILER: gcc7 - BUILD_TYPE: readme + CONFIG_TYPE: readme diff --git a/jenkins/axis_hpc_nightly.yml b/jenkins/axis_hpc_nightly.yml new file mode 100644 index 0000000..8b1290e --- /dev/null +++ b/jenkins/axis_hpc_nightly.yml @@ -0,0 +1,21 @@ +# axis_hpc_nightly.yml (some other combinations will get weeded out inside the script) +COMPILER: + - intel18 + - gcc6 + - gcc7 + +CONFIG_TYPE: + - base + - debug + - serial + - readme + - thrust + - kokkos + +exclude: + - COMPILER: gcc6 + CONFIG_TYPE: readme + - COMPILER: gcc7 + CONFIG_TYPE: readme + - COMPILER: intel18 + CONFIG_TYPE: kokkos diff --git a/jenkins/axis_hpc_pr.yml b/jenkins/axis_hpc_pr.yml index 8d692a8..6fc2448 100644 --- a/jenkins/axis_hpc_pr.yml +++ b/jenkins/axis_hpc_pr.yml @@ -2,7 +2,7 @@ COMPILER: - intel18 - gcc7 -BUILD_TYPE: +CONFIG_TYPE: - base - debug - serial @@ -11,5 +11,5 @@ BUILD_TYPE: exclude: - COMPILER: gcc7 - BUILD_TYPE: readme + CONFIG_TYPE: readme diff --git a/jenkins/axis_varan.yml b/jenkins/axis_varan.yml deleted file mode 100644 index ccb9a39..0000000 --- a/jenkins/axis_varan.yml +++ /dev/null @@ -1,18 +0,0 @@ -# axis_varan.yml -COMPILER: - - intel18 - - gcc6 - - gcc7 - -BUILD_TYPE: - - base - - debug - - serial - - readme - - thrust - -exclude: - - COMPILER: gcc6 - BUILD_TYPE: readme - - COMPILER: gcc7 - BUILD_TYPE: readme diff --git a/jenkins/axis_varan_install.yml b/jenkins/axis_varan_install.yml index 9c54eb7..700e2c8 100644 --- a/jenkins/axis_varan_install.yml +++ b/jenkins/axis_varan_install.yml @@ -4,13 +4,14 @@ COMPILER: - gcc6 - gcc7 -BUILD_TYPE: +CONFIG_TYPE: - base + - debug - serial - thrust - kokkos exclude: - COMPILER: intel18 - BUILD_TYPE: kokkos + CONFIG_TYPE: kokkos diff --git a/jenkins/axis_varan_nightly.yml b/jenkins/axis_varan_nightly.yml new file mode 100644 index 0000000..8a7d0a2 --- /dev/null +++ b/jenkins/axis_varan_nightly.yml @@ -0,0 +1,21 @@ +# axis_varan.yml (some other combinations will get weeded out inside the script) +COMPILER: + - intel18 + - gcc6 + - gcc7 + +CONFIG_TYPE: + - base + - debug + - serial + - readme + - thrust + - kokkos + +exclude: + - COMPILER: gcc6 + CONFIG_TYPE: readme + - COMPILER: gcc7 + CONFIG_TYPE: readme + - COMPILER: intel18 + CONFIG_TYPE: kokkos diff --git a/jenkins/axis_varan_pr.yml b/jenkins/axis_varan_pr.yml index 45113be..06897c5 100644 --- a/jenkins/axis_varan_pr.yml +++ b/jenkins/axis_varan_pr.yml @@ -3,7 +3,7 @@ COMPILER: - gcc6 - gcc7 -BUILD_TYPE: +CONFIG_TYPE: - base - debug - serial @@ -12,7 +12,7 @@ BUILD_TYPE: exclude: - COMPILER: gcc6 - BUILD_TYPE: readme + CONFIG_TYPE: readme - COMPILER: gcc7 - BUILD_TYPE: readme + CONFIG_TYPE: readme diff --git a/jenkins/build_docs.sh b/jenkins/build_docs.sh index c729013..9cf2166 100755 --- a/jenkins/build_docs.sh +++ b/jenkins/build_docs.sh @@ -9,6 +9,9 @@ set -e # Echo each command set -x +# Set umask so installations have rwx permissions for the group +umask 007 + JALI_VERSION=1.0.0 openmpi_version=2.1.2 diff --git a/jenkins/build_matrix_entry.sh b/jenkins/build_matrix_entry.sh deleted file mode 100755 index b27a5cd..0000000 --- a/jenkins/build_matrix_entry.sh +++ /dev/null @@ -1,134 +0,0 @@ -#!/usr/bin/env bash -# This script is executed on Jenkins using -# -# $WORKSPACE/jenkins/build_matrix_entry.sh -# -# The exit code determines if the test succeeded or failed. -# Note that the environment variable WORKSPACE must be set (Jenkins -# will do this automatically). - -# Exit on error -set -e -# Echo each command -set -x - -compiler=$1 -build_type=$2 - -echo "inside build_matrix entry" - -# special case for README builds -if [[ $build_type == "readme" ]]; then - python2 $WORKSPACE/jenkins/parseREADME.py $WORKSPACE/README.md $WORKSPACE - exit -fi -# special case for README builds -if [[ $build_type == "readme" ]]; then - - # Put a couple of settings in place to generate test output even if - # the README doesn't ask for it. - export CTEST_OUTPUT_ON_FAILURE=1 - CACHE_OPTIONS="-D ENABLE_JENKINS_OUTPUT=True" - sed "s/^ *cmake/& $CACHE_OPTIONS/g" $WORKSPACE/README.md >$WORKSPACE/README.md.1 - python2 $WORKSPACE/jenkins/parseREADME.py \ - $WORKSPACE/README.md.1 \ - $WORKSPACE \ - varan - exit - -fi - -# set modules and install paths - -jali_version=1.1.4 -lapack_version=3.8.0 - -export NGC=/usr/local/codes/ngc -ngc_include_dir=$NGC/private/include - -thrust_dir=${ngc_include_dir} - - -# compiler-specific settings -if [[ $compiler == "intel18" ]]; then - compiler_version=18.0.1 - cxxmodule=intel/${compiler_version} - compiler_suffix="-intel-${compiler_version}" - - openmpi_version=2.1.2 - mpi_module=openmpi/2.1.2 - mpi_suffix="-openmpi-${openmpi_version}" - -elif [[ $compiler =~ "gcc" ]]; then - - if [[ $compiler == "gcc6" ]]; then - compiler_version=6.4.0 - elif [[ $compiler == "gcc7" ]]; then - compiler_version=7.3.0 - fi - cxxmodule=gcc/${compiler_version} - compiler_suffix="-gcc-${compiler_version}" - - openmpi_version=2.1.2 - mpi_module=openmpi/${openmpi_version} - mpi_suffix="-openmpi-${openmpi_version}" - -fi -jali_install_dir=$NGC/private/jali/${jali_version}${compiler_suffix}${mpi_suffix} -lapacke_dir=$NGC/private/lapack/${lapack_version}-patched${compiler_suffix} - - -# build-type-specific settings -cmake_build_type=Release -mpi_flags="-D WONTON_ENABLE_MPI=True" -extra_flags= -thrust_flags= -jali_flags="-D WONTON_ENABLE_Jali=True -D Jali_ROOT:FILEPATH=$jali_install_dir" -lapacke_flags="-D WONTON_ENABLE_LAPACKE=True -D LAPACKE_ROOT:FILEPATH=$lapacke_dir" -if [[ $build_type == "debug" ]]; then - cmake_build_type=Debug -elif [[ $build_type == "serial" ]]; then - mpi_flags= - jali_flags= # jali is not available in serial -elif [[ $build_type == "thrust" ]]; then - thrust_flags="-D WONTON_ENABLE_THRUST=True -DTHRUST_ROOT:FILEPATH=${thrust_dir}" -fi - -flecsi_flags="-D WONTON_ENABLE_FleCSI:BOOL=False" -if [[ $compiler == "gcc6" && $build_type != "serial" ]]; then - flecsi_install_dir=$NGC/private/flecsi/374b56b-gcc-6.4.0 - flecsisp_install_dir=$NGC/private/flecsi-sp/e78c594-gcc-6.4.0 - flecsi_flags="-D WONTON_ENABLE_FleCSI:BOOL=True -D FleCSI_ROOT:PATH=$flecsi_install_dir -D FleCSISP_ROOT:PATH=$flecsisp_install_dir" -fi - -export SHELL=/bin/sh - -export MODULEPATH="" -. /opt/local/packages/Modules/default/init/sh -module load $cxxmodule -if [[ -n "$mpi_flags" ]]; then - module load ${mpi_module} -fi -module load cmake/3.14.0 # 3.13 or higher is required - -echo $WORKSPACE -cd $WORKSPACE - -mkdir build -cd build - -cmake \ - -D CMAKE_BUILD_TYPE=$cmake_build_type \ - -D CMAKE_CXX_FLAGS="-Wall -Werror" \ - -D CMAKE_PREFIX_PATH=$ngc_include_dir \ - -D ENABLE_UNIT_TESTS=True \ - -D ENABLE_APP_TESTS=True \ - -D ENABLE_JENKINS_OUTPUT=True \ - $mpi_flags \ - $jali_flags \ - $flecsi_flags \ - $lapacke_flags \ - $thrust_flags \ - .. -make -j2 -ctest -j2 --output-on-failure diff --git a/jenkins/build_matrix_entry_hpc.sh b/jenkins/build_matrix_entry_hpc.sh index acd8725..45fd354 100755 --- a/jenkins/build_matrix_entry_hpc.sh +++ b/jenkins/build_matrix_entry_hpc.sh @@ -1,24 +1,62 @@ #!/usr/bin/env bash # This script is executed on Jenkins using # -# $WORKSPACE/jenkins/build_matrix_entry.sh +# $WORKSPACE/jenkins/build_matrix_entry_hpc.sh BUILD_TYPE # -# The exit code determines if the test succeeded or failed. -# Note that the environment variable WORKSPACE must be set (Jenkins +# BUILD_TYPE - pr, nightly, install +# +# if VER is absent, the HEAD of the master branch will be taken +# (except for kokkos config_type which takes the HEAD of the kokkos +# branch). If BUILD_TYPE is 'install' and VER is specified, it will +# install it to /install_prefix/wonton/${VER}-blah-blah; if VER is not +# specified, it will install to /install_prefix/wonton/dev-blah-blah +# +# Note that the following environment variables must be set (Jenkins # will do this automatically). +# +# WORKSPACE - where the code is checked out +# CONFIG_TYPE - base, debug, serial, readme, thrust, kokkos +# COMPILER - intel, gcc6, gcc7 +# BRANCH_NAME - master, kokkos +# +# The exit code determines if the test succeeded or failed. # Exit on error set -e # Echo each command set -x -compiler=$1 -build_type=$2 +# Set umask so installations have rwx permissions for the group +umask 007 + +BUILD_TYPE=$1 +version=$2 +if [[ $version == "" ]]; then + version=dev +fi + + +# Don't build kokkos config from master branch and general configs +# from kokkos branch +if [[ $BRANCH_NAME == "master" && $CONFIG_TYPE == "kokkos" ]]; then + exit +fi +if [[ $BRANCH_NAME == "kokkos" && $CONFIG_TYPE != "kokkos" ]]; then + exit +fi + + +# versions of dependencies +jali_version=1.1.4 +thrust_version=1.8.3 +kokkos_version=3.2.00-openmp +lapack_version=3.8.0 + +echo "inside build_matrix on PLATFORM=$PLATFORM with BUILD_TYPE=$BUILD_TYPE $CONFIG_TYPE=$CONFIG_TYPE COMPILER=$COMPILER" -echo "inside build_matrix entry" # special case for README builds -if [[ $build_type == "readme" ]]; then +if [[ $BUILD_TYPE != "install" && $CONFIG_TYPE == "readme" ]]; then # Put a couple of settings in place to generate test output even if # the README doesn't ask for it. @@ -35,68 +73,94 @@ fi # set modules and install paths -jali_version=1.1.4 -lapack_version=3.8.0 - export NGC=/usr/projects/ngc ngc_include_dir=$NGC/private/include -thrust_dir=${ngc_include_dir} - # compiler-specific settings -if [[ $compiler == "intel18" ]]; then - compiler_version=18.0.5 - cxxmodule=intel/${compiler_version} - compiler_suffix="-intel-${compiler_version}" +if [[ $COMPILER =~ "intel" ]]; then - openmpi_version=2.1.2 - mpi_module=openmpi/2.1.2 - mpi_suffix="-openmpi-${openmpi_version}" + compiler_version=18.0.5 + cxxmodule=intel/${compiler_version} + compiler_suffix="-intel-${compiler_version}" -elif [[ $compiler =~ "gcc" ]]; then + openmpi_version=2.1.2 + mpi_module=openmpi/2.1.2 + mpi_suffix="-openmpi-${openmpi_version}" + +elif [[ $COMPILER =~ "gcc" ]]; then + + openmpi_version=2.1.2 + if [[ $COMPILER == "gcc6" ]]; then + compiler_version=6.4.0 + elif [[ $COMPILER == "gcc7" ]]; then + compiler_version=7.4.0 + fi + cxxmodule=gcc/${compiler_version} + compiler_suffix="-gcc-${compiler_version}" - if [[ $compiler == "gcc6" ]]; then - compiler_version=6.4.0 - elif [[ $compiler == "gcc7" ]]; then - compiler_version=7.4.0 - fi - cxxmodule=gcc/${compiler_version} - compiler_suffix="-gcc-${compiler_version}" - - openmpi_version=2.1.2 - mpi_module=openmpi/${openmpi_version} - mpi_suffix="-openmpi-${openmpi_version}" - + mpi_module=openmpi/${openmpi_version} + mpi_suffix="-openmpi-${openmpi_version}" + fi + +# Jali jali_install_dir=$NGC/private/jali/${jali_version}${compiler_suffix}${mpi_suffix} +jali_flags="-D WONTON_ENABLE_Jali:BOOL=True -D Jali_ROOT:FILEPATH=$jali_install_dir" + +# LAPACKE lapacke_dir=$NGC/private/lapack/${lapack_version}-patched${compiler_suffix} +lapacke_flags="-D WONTON_ENABLE_LAPACKE:BOOL=True -D LAPACKE_ROOT:FILEPATH=$lapacke_dir" +# Flecsi - Not yet installed on Snow +flecsi_flags="-D WONTON_ENABLE_FleCSI:BOOL=False" +# if [[ $COMPILER == "gcc6" ]]; then +# flecsi_install_dir=$NGC/private/flecsi/374b56b-gcc-6.4.0 +# flecsisp_install_dir=$NGC/private/flecsi-sp/e78c594-gcc-6.4.0 +# flecsi_flags="-D WONTON_ENABLE_FleCSI:BOOL=True -D FleCSI_ROOT:PATH=$flecsi_install_dir -D FleCSISP_ROOT:PATH=$flecsisp_install_dir" +# fi -# build-type-specific settings -cmake_build_type=Release -mpi_flags="-D WONTON_ENABLE_MPI=True" -extra_flags= +# THRUST thrust_flags= -jali_flags="-D WONTON_ENABLE_Jali=True -D Jali_ROOT:FILEPATH=$jali_install_dir" -lapacke_flags="-D WONTON_ENABLE_LAPACKE=True -D LAPACKE_ROOT:FILEPATH=$lapacke_dir" -if [[ $build_type == "debug" ]]; then +thrust_suffix= +thrust_install_dir=$NGC/private/thrust/${thrust_version} +if [[ $CONFIG_TYPE == "thrust" ]]; then + thrust_flags="-D WONTON_ENABLE_THRUST=True -D THRUST_ROOT=${thrust_install_dir}" + thrust_suffix="-thrust" +fi + +# Kokkos +kokkos_flags= +kokkos_suffix= +kokkos_install_dir=$NGC/private/kokkos/${kokkos_version}${compiler_suffix} +if [[ $CONFIG_TYPE == "kokkos" ]]; then + kokkos_flags="-D WONTON_ENABLE_Kokkos=True -D Kokkos_ROOT=${kokkos_install_dir}" + kokkos_suffix="-kokkos" +fi + +# MPI or not +mpi_flags="-D WONTON_ENABLE_MPI=True" +if [[ $CONFIG_TYPE == "serial" ]]; then + mpi_flags="-D WONTON_ENABLE_MPI=False" + mpi_suffix= + jali_flags= + flecsi_flags= +fi + +# Debug or Optimized build +cmake_build_type=Release +debug_suffix= +if [[ $CONFIG_TYPE == "debug" ]]; then cmake_build_type=Debug -elif [[ $build_type == "serial" ]]; then - mpi_flags= - jali_flags= # jali is not available in serial -elif [[ $build_type == "thrust" ]]; then - thrust_flags="-D WONTON_ENABLE_THRUST=True -DTHRUST_ROOT:FILEPATH=${thrust_dir}" + debug_suffix="-debug" fi -# We don't seem to have a suitable (read as "ancient") flecsi install on HPC -#if [[ $compiler == "gcc6" && $build_type != "serial" ]]; then -# flecsi_install_dir=$NGC/private/flecsi/374b56b-gcc-6.4.0 -# flecsisp_install_dir=$NGC/private/flecsi-sp/e78c594-gcc-6.4.0 -# flecsi_flags="-D WONTON_ENABLE_FleCSI:BOOL=True -D FleCSI_ROOT:PATH=$flecsi_install_dir -D FleCSISP_ROOT:PATH=$flecsisp_install_dir" -#fi +# Build up an install dir name +wonton_install_dir=$NGC/private/wonton/${version}${compiler_suffix}${mpi_suffix}${thrust_suffix}${kokkos_suffix}${debug_suffix} + export SHELL=/bin/sh + #Rely on default user environment to load modules; these scripts can be found in /etc/profile.d/ as of 8/25/20 the scripts that set up modules on snow are /etc/profile.d/z00_lmod.sh; /etc/profile.d/00-modulepath.sh; /etc/profile.d/z01-modules.lanl.sh; module load $cxxmodule if [[ -n "$mpi_flags" ]]; then @@ -104,15 +168,17 @@ if [[ -n "$mpi_flags" ]]; then fi module load cmake/3.14.0 # 3.13 or higher is required -echo $WORKSPACE +echo "JENKINS WORKSPACE = $WORKSPACE" cd $WORKSPACE +rm -fr build mkdir build cd build cmake \ -D CMAKE_BUILD_TYPE=$cmake_build_type \ -D CMAKE_CXX_FLAGS="-Wall -Werror" \ + -D CMAKE_INSTALL_PREFIX=$wonton_install_dir \ -D CMAKE_PREFIX_PATH=$ngc_include_dir \ -D ENABLE_UNIT_TESTS=True \ -D ENABLE_APP_TESTS=True \ @@ -122,6 +188,11 @@ cmake \ $flecsi_flags \ $lapacke_flags \ $thrust_flags \ + $kokkos_flags \ .. make -j36 ctest -j36 --output-on-failure + +if [[ $BUILD_TYPE == "install" ]]; then + make install +fi diff --git a/jenkins/install_matrix_entry.sh b/jenkins/build_matrix_entry_varan.sh similarity index 52% rename from jenkins/install_matrix_entry.sh rename to jenkins/build_matrix_entry_varan.sh index 188951e..ddaab6b 100755 --- a/jenkins/install_matrix_entry.sh +++ b/jenkins/build_matrix_entry_varan.sh @@ -1,37 +1,83 @@ #!/usr/bin/env bash # This script is executed on Jenkins using # -# $WORKSPACE/jenkins/build_matrix_entry.sh +# $WORKSPACE/jenkins/build_matrix_entry_varan.sh BUILD_TYPE # -# The exit code determines if the test succeeded or failed. -# Note that the environment variable WORKSPACE must be set (Jenkins +# BUILD_TYPE - pr, nightly, install +# +# if VER is absent, the HEAD of the master branch will be taken +# (except for kokkos config_type which takes the HEAD of the kokkos +# branch). If BUILD_TYPE is 'install' it will install it to +# /install_prefix/wonton/dev-blah-blah +# +# Note that the following environment variables must be set (Jenkins # will do this automatically). +# +# WORKSPACE - where the code is checked out +# CONFIG_TYPE - base, debug, serial, readme, thrust, kokkos +# COMPILER - intel, gcc6, gcc7 +# BRANCH_NAME - master, kokkos +# +# The exit code determines if the test succeeded or failed. # Exit on error set -e # Echo each command set -x -compiler=$1 -build_type=$2 +# Set umask so installations have rwx permissions for the group +umask 007 -echo "inside install_matrix entry" +BUILD_TYPE=$1 +version=$2 +if [[ $version == "" ]]; then + version=dev +fi -# set modules and install paths +# Don't build kokkos config from master branch and general configs +# from kokkos branch +if [[ $BRANCH_NAME == "master" && $CONFIG_TYPE == "kokkos" ]]; then + exit +fi +if [[ $BRANCH_NAME == "kokkos" && $CONFIG_TYPE != "kokkos" ]]; then + exit +fi + + +# versions of dependencies jali_version=1.1.4 +thrust_version=1.8.3 +kokkos_version=3.1.01 lapack_version=3.8.0 +echo "inside build_matrix on PLATFORM=$PLATFORM with BUILD_TYPE=$BUILD_TYPE $CONFIG_TYPE=$CONFIG_TYPE COMPILER=$COMPILER" + + +# special case for README builds +if [[ $BUILD_TYPE != "install" && $CONFIG_TYPE == "readme" ]]; then + + # Put a couple of settings in place to generate test output even if + # the README doesn't ask for it. + export CTEST_OUTPUT_ON_FAILURE=1 + CACHE_OPTIONS="-D ENABLE_JENKINS_OUTPUT=True" + sed "s/^ *cmake/& $CACHE_OPTIONS/g" $WORKSPACE/README.md >$WORKSPACE/README.md.1 + python2 $WORKSPACE/jenkins/parseREADME.py \ + $WORKSPACE/README.md.1 \ + $WORKSPACE \ + varan + exit + +fi + +# set modules and install paths + export NGC=/usr/local/codes/ngc ngc_include_dir=$NGC/private/include -wonton_install_dir=$NGC/private/wonton -wonton_version=$3 - -thrust_dir=${ngc_include_dir} # compiler-specific settings -if [[ $compiler =~ "intel" ]]; then +if [[ $COMPILER =~ "intel" ]]; then compiler_version=18.0.1 cxxmodule=intel/${compiler_version} @@ -41,18 +87,14 @@ if [[ $compiler =~ "intel" ]]; then mpi_module=openmpi/2.1.2 mpi_suffix="-openmpi-${openmpi_version}" -elif [[ $compiler =~ "gcc" ]]; then +elif [[ $COMPILER =~ "gcc" ]]; then openmpi_version=2.1.2 - if [[ $compiler == "gcc6" ]]; then + if [[ $COMPILER == "gcc6" ]]; then compiler_version=6.4.0 - elif [[ $compiler == "gcc7" ]]; then + elif [[ $COMPILER == "gcc7" ]]; then compiler_version=7.3.0 - elif [[ $compiler == "gcc8" ]]; then - compiler_ersion=8.2.0 - openmpi_version=3.1.3 - fi - + fi cxxmodule=gcc/${compiler_version} compiler_suffix="-gcc-${compiler_version}" @@ -61,48 +103,60 @@ elif [[ $compiler =~ "gcc" ]]; then fi +# Jali jali_install_dir=$NGC/private/jali/${jali_version}${compiler_suffix}${mpi_suffix} -jali_flags="-D WONTON_ENABLE_Jali:BOOL=True -D Jali_ROOT:PATH=$jali_install_dir" +jali_flags="-D WONTON_ENABLE_Jali:BOOL=True -D Jali_ROOT:FILEPATH=$jali_install_dir" +# LAPACKE lapacke_dir=$NGC/private/lapack/${lapack_version}-patched${compiler_suffix} -lapacke_flags="-D WONTON_ENABLE_LAPACKE:BOOL=True -D LAPACKE_ROOT:PATH=$lapacke_dir" +lapacke_flags="-D WONTON_ENABLE_LAPACKE:BOOL=True -D LAPACKE_ROOT:FILEPATH=$lapacke_dir" +# Flecsi flecsi_flags="-D WONTON_ENABLE_FleCSI:BOOL=False" -if [[ $compiler == "gcc6" ]]; then +if [[ $COMPILER == "gcc6" ]]; then flecsi_install_dir=$NGC/private/flecsi/374b56b-gcc-6.4.0 flecsisp_install_dir=$NGC/private/flecsi-sp/e78c594-gcc-6.4.0 flecsi_flags="-D WONTON_ENABLE_FleCSI:BOOL=True -D FleCSI_ROOT:PATH=$flecsi_install_dir -D FleCSISP_ROOT:PATH=$flecsisp_install_dir" fi - -mpi_flags="-D WONTON_ENABLE_MPI=True" -if [[ $build_type == "serial" ]]; then - mpi_flags="-D WONTON_ENABLE_MPI=False" - mpi_suffix= - jali_flags= - flecsi_flags= -fi - - +# THRUST thrust_flags= thrust_suffix= -thrust_version=1.8.3 thrust_install_dir=$NGC/private/thrust/${thrust_version} -if [[ $build_type == "thrust" ]]; then +if [[ $CONFIG_TYPE == "thrust" ]]; then thrust_flags="-D WONTON_ENABLE_THRUST=True -D THRUST_ROOT=${thrust_install_dir}" thrust_suffix="-thrust" fi +# Kokkos kokkos_flags= kokkos_suffix= -kokkos_version=3.1.01 kokkos_install_dir=$NGC/private/kokkos/${kokkos_version}${compiler_suffix} -if [[ $build_type == "kokkos" ]]; then +if [[ $CONFIG_TYPE == "kokkos" ]]; then kokkos_flags="-D WONTON_ENABLE_Kokkos=True -D Kokkos_ROOT=${kokkos_install_dir}" kokkos_suffix="-kokkos" fi -wonton_install_dir=$wonton_install_dir/${wonton_version}${compiler_suffix}${mpi_suffix}${thrust_suffix}${kokkos_suffix} +# MPI or not +mpi_flags="-D WONTON_ENABLE_MPI=True" +if [[ $CONFIG_TYPE == "serial" ]]; then + mpi_flags="-D WONTON_ENABLE_MPI=False" + mpi_suffix= + jali_flags= + flecsi_flags= +fi + +# Debug or Optimized build +cmake_build_type=Release +debug_suffix= +if [[ $CONFIG_TYPE == "debug" ]]; then + cmake_build_type=Debug + debug_suffix="-debug" +fi + +# Build up an install dir name +wonton_install_dir=$NGC/private/wonton/${version}${compiler_suffix}${mpi_suffix}${thrust_suffix}${kokkos_suffix}${debug_suffix} + export SHELL=/bin/sh @@ -116,21 +170,11 @@ module load cmake/3.14.0 # 3.13 or higher is required module load git -if [[ $BRANCH_NAME == "master" && build_type == "kokkos" ]]; then - exit -fi - -if [[ $BRANCH_NAME == "kokkos" && build_type != "kokkos" ]]; then - exit -fi - -cmake_build_type=Release - echo "JENKINS WORKSPACE = $WORKSPACE" cd $WORKSPACE rm -fr build -mkdir -p build +mkdir build cd build cmake \ @@ -149,5 +193,8 @@ cmake \ $kokkos_flags \ .. make -j2 -ctest --output-on-failure -make install +ctest -j2 --output-on-failure + +if [[ $BUILD_TYPE == "install" ]]; then + make install +fi diff --git a/wonton/CMakeLists.txt b/wonton/CMakeLists.txt index e7cb28e..7506664 100644 --- a/wonton/CMakeLists.txt +++ b/wonton/CMakeLists.txt @@ -5,7 +5,12 @@ Please see the license file at the root of this repository, or at: ]] add_subdirectory(support) + +if (WONTON_ENABLE_MPI) + add_subdirectory(distributed) +endif () + add_subdirectory(intersect) add_subdirectory(mesh) add_subdirectory(state) -add_subdirectory(swarm) \ No newline at end of file +add_subdirectory(swarm) diff --git a/wonton/distributed/CMakeLists.txt b/wonton/distributed/CMakeLists.txt new file mode 100644 index 0000000..8105d73 --- /dev/null +++ b/wonton/distributed/CMakeLists.txt @@ -0,0 +1,57 @@ +#[[ +This file is part of the Ristra Wonton project. +Please see the license file at the root of this repository, or at: + https://github.com/laristra/wonton/blob/master/LICENSE +]] +#-----------------------------------------------------------------------------~# +project(wonton_distributed) + +add_library(wonton_distributed INTERFACE) + +target_include_directories(wonton_distributed INTERFACE + $ + $ + $) + +set(wonton_distributed_HEADERS + mpi_ghost_manager.h + ) + +# Not yet allowed for INTERFACE libraries +# +# set_target_properties(wonton_distributed PROPERTIES +# PUBLIC_HEADER ${wonton_distributed_HEADERS} +# ) +# +# Directly install files instead +install(FILES ${wonton_distributed_HEADERS} DESTINATION include/wonton/distributed) + +target_link_libraries(wonton_distributed INTERFACE wonton_mesh wonton_state) + +target_link_libraries(wonton INTERFACE wonton_distributed) + +install(TARGETS wonton_distributed + EXPORT wonton_LIBRARIES + ARCHIVE DESTINATION lib + LIBRARY DESTINATION lib + RUNTIME DESTINATION bin + PUBLIC_HEADER DESTINATION include/wonton/distributed + INCLUDES DESTINATION include/wonton/distributed + ) + +# Unit tests + +if (ENABLE_UNIT_TESTS) + + if (WONTON_ENABLE_Jali) + wonton_add_unittest(test_mpi_ghost_manager + SOURCES test/test_mpi_ghost_manager.cc + LIBRARIES wonton_distributed + POLICY MPI + THREADS 4 + ) + endif () + +endif () + + diff --git a/wonton/distributed/mpi_ghost_manager.h b/wonton/distributed/mpi_ghost_manager.h new file mode 100644 index 0000000..c8cef3f --- /dev/null +++ b/wonton/distributed/mpi_ghost_manager.h @@ -0,0 +1,731 @@ +/* + * This file is part of the Ristra Wonton project. + * Please see the license file at the root of this repository, or at: + * https://github.com/laristra/wonton/blob/master/LICENSE + */ + +#ifndef WONTON_MPI_GHOST_MANAGER_H +#define WONTON_MPI_GHOST_MANAGER_H + +#include +#include +#include "wonton/support/wonton.h" +#include "wonton/support/Point.h" + + +#ifdef WONTON_ENABLE_MPI +namespace Wonton { + +/** + * @brief A communication manager to fill values on ghost cells. + * + * @tparam Mesh: the mesh type. + * @tparam State: its state type. + * @tparam entity: the entity kind. + */ +template +class MPI_GhostManager { + +public: + /** + * @brief Create and initialize the manager. + * + * @param in_mesh: the given mesh. + * @param in_state: its state. + * @param in_comm: the MPI communicator. + */ + MPI_GhostManager(Mesh const& in_mesh, State& in_state, MPI_Comm in_comm) + : mesh(in_mesh), + state(in_state), + comm(in_comm) + { cache_comm_matrices(); } + + /** + * @brief Delete the manager. + * + */ + ~MPI_GhostManager() = default; + + /** + * @brief Send owned entities values and receive ghost entities values. + * + * @param field: the field to be remapped. + * @param cache: whether to cache values or not for this field. + */ + void update_ghost_values(std::string const& field, bool cache = false) { + + auto field_type = state.field_type(entity, field); + bool multimat = (field_type == Field_type::MULTIMATERIAL_FIELD); + + std::type_info const& data_type = state.get_data_type(field); + + if (multimat) { + assert(entity == Wonton::CELL); + for (int m = 1; m < num_mats; ++m) { + if (data_type == typeid(double)) + update_material_field(field, m, cache); + else if (data_type == typeid(Vector<1>)) + update_material_field>(field, m, cache); + else if (data_type == typeid(Vector<2>)) + update_material_field>(field, m, cache); + else if (data_type == typeid(Vector<3>)) + update_material_field>(field, m, cache); + else if (data_type == typeid(Point<1>)) + update_material_field>(field, m, cache); + else if (data_type == typeid(Point<2>)) + update_material_field>(field, m, cache); + else if (data_type == typeid(Point<3>)) + update_material_field>(field, m, cache); + else + throw std::runtime_error("mpi_ghost_manager.h: incorrect material field data type"); + } + } else { // mesh field + if (data_type == typeid(double)) + update_mesh_field(field, cache); + else if (data_type == typeid(Vector<1>)) + update_mesh_field>(field, cache); + else if (data_type == typeid(Vector<2>)) + update_mesh_field>(field, cache); + else if (data_type == typeid(Vector<3>)) + update_mesh_field>(field, cache); + else if (data_type == typeid(Point<1>)) + update_mesh_field>(field, cache); + else if (data_type == typeid(Point<2>)) + update_mesh_field>(field, cache); + else if (data_type == typeid(Point<3>)) + update_mesh_field>(field, cache); + else + throw std::runtime_error("mpi_ghost_manager.h: incorrect material field data type"); + } + } + + /** @brief Update ghosts of compact array of material cell values + * + * @param[in|out] material_values: compact material cell data + * @param[in] m: material ID + * @param[in] cache: whether to cache values or not for this field. + * + * Here 'material_values' is a compact array of material cell values, + * sized to num_owned + num_ghost but with only the owned values + * filled in. The ghost values are filled in by this routine. + * Note that the material ID offset is zero unlike in 'update_ghost_values'. + */ + template + void update_material_ghost_values(T* material_values, int m, bool cache = false) { + + // convert to mesh cell-centric field + int const num_mesh_cells = mesh.num_owned_cells() + mesh.num_ghost_cells(); + std::vector mesh_values(num_mesh_cells); + + std::vector material_cells; + state.mat_get_cells(m, &material_cells); // gives ALL cells OWNED+GHOST + int const num_material_cells = material_cells.size(); + + for (int i = 0; i < num_material_cells; i++) { + mesh_values[material_cells[i]] = material_values[i]; + } + + update_ghost_values(mesh_values.data(), m + 1, cache); + + // convert back to material cell-centric field + for (auto&& c : material_cells) { + int const j = state.cell_index_in_material(c, m); + material_values[j] = mesh_values[c]; + } + } + + + /** @brief Update ghosts of compact array of material cell values + * + * @param meshdata [IN/OUT}: data on mesh entities + * @param cache: whether to cache values or not for this field. + * + * Sized to nowned+nghost entities but with only owned entities filled in + */ + template + void update_mesh_ghost_values(T* values, bool cache = false) { + update_ghost_values(values, 0, cache); + } + + + /** + * @brief Retrieve the list of owned entities to send to each rank. + * + * @param m: material index. + * @return matrix of communication. + */ + std::vector> const& owned_matrix(int m = 0) const { return send.matrix[m]; } + + /** + * @brief Retrieve the list of ghost entities to receive from each rank. + * + * @param m: material index. + * @return matrix of communication. + */ + std::vector> const& ghost_matrix(int m = 0) const { return take.matrix[m]; } + + /** + * @brief Retrieve values of owned entities sent to each rank. + * + * @return matrix of values. + */ + std::vector> const& owned_values(int m = 0) const { return send.values[m]; } + + /** + * @brief Retrieve values of ghost entities received from each rank. + * + * @return matrix of values. + */ + std::vector> const& ghost_values(int m = 0) const { return take.values[m]; } + +private: + /** + * @brief Build and store communication matrices to identify + * which owned entities data should be sent to each rank, + * and which ghost entities data should be received from + * each rank for each material. + */ + void cache_comm_matrices() { + + // step 0: preprocessing + MPI_Comm_rank(comm, &rank); + MPI_Comm_size(comm, &num_ranks); + + num_mats = state.num_materials() + 1; + send.matrix.resize(num_mats); + send.count.resize(num_mats); + take.matrix.resize(num_mats); + take.count.resize(num_mats); + + int const num_entities = mesh.num_entities(entity, Wonton::ALL); + + for (int i = 0; i < num_entities; ++i) { + int const& gid = mesh.get_global_id(i, entity); + auto& exchange = is_ghost(i) ? take : send; + exchange.lookup[gid] = i; + } + + for (int m = 0; m < num_mats; ++m) { + + send.matrix[m].resize(num_ranks); + take.matrix[m].resize(num_ranks); + send.count[m].resize(num_ranks); + take.count[m].resize(num_ranks); + + std::vector entities; + int num_ghosts[num_ranks]; + int offsets[num_ranks]; + std::vector received; + std::vector gids[num_ranks]; + + // step 1: cache local ghost entities list + if (m > 0) /* multi-material */{ + assert(entity == Wonton::CELL && "only for cell-centered fields"); + std::vector cells; + state.mat_get_cells(m - 1, &cells); + for (int const& i : cells) { + if (is_ghost(i)) { + int const& gid = mesh.get_global_id(i, entity); + entities.emplace_back(gid); + } + } + } else { + for (auto&& id : take.lookup) { entities.emplace_back(id.first); } + } + + + // step 2: gather number of ghosts for all ranks and deduce offsets + int local_num_ghosts = entities.size(); + MPI_Allgather(&local_num_ghosts, 1, MPI_INT, num_ghosts, 1, MPI_INT, comm); + + int const total_ghosts = std::accumulate(num_ghosts, num_ghosts + num_ranks, 0); + received.resize(total_ghosts); + + int index = 0; + for (int i = 0; i < num_ranks; ++i) { + offsets[i] = index; + index += num_ghosts[i]; + } + + // step 3: gather all ghost entities on all ranks + MPI_Allgatherv(entities.data(), num_ghosts[rank], MPI_INT, received.data(), num_ghosts, offsets, MPI_INT, comm); + + // step 4: check received ghost cells, build send matrix and send entities count. + int num_sent[num_ranks]; + std::vector requests; + + for (int i = 0; i < num_ranks; ++i) { + if (i != rank) { + int const& start = offsets[i]; + int const& extent = (i < num_ranks - 1 ? offsets[i+1] : total_ghosts); + for (int j = start; j < extent; ++j) { + int const& gid = received[j]; + if (send.lookup.count(gid)) { + send.matrix[m][i].emplace_back(send.lookup[gid]); + gids[i].emplace_back(gid); + } + } + + MPI_Request request; + num_sent[i] = send.matrix[m][i].size(); + MPI_Isend(num_sent + i, 1, MPI_INT, i, 0, comm, &request); + requests.emplace_back(request); + } + } + +#if !defined(NDEBUG) && defined(VERBOSE_OUTPUT) + for (int i = 0; i < num_ranks; ++i) { + if (i != rank) { + std::cout << "[" << rank << "] -> [" << i << "]: ("; + int const num_to_send = send.matrix[m][i].size(); + for (int j = 0; j < num_to_send; ++j) { + std::cout << mesh.get_global_id(send.matrix[m][i][j], entity); + if (j < num_to_send - 1) { + std::cout << ", "; + } + } + std::cout << ")" << std::endl; + } + } +#endif + + // step 5: receive ghost count per rank + for (int i = 0; i < num_ranks; ++i) { + if (i != rank) { + MPI_Request request; + MPI_Irecv(take.count[m].data() + i, 1, MPI_INT, i, 0, comm, &request); + requests.emplace_back(request); + } + } + + // step 6: send give entities, receive expected ghosts. + for (int i = 0; i < num_ranks; ++i) { + if (i != rank and num_sent[i]) { + MPI_Request request; + MPI_Isend(gids[i].data(), num_sent[i], MPI_INT, i, 1, comm, &request); + requests.emplace_back(request); + } + } + + MPI_Waitall(requests.size(), requests.data(), status); + requests.clear(); + + for (int i = 0; i < num_ranks; ++i) { + if (i != rank and take.count[m][i]) { + MPI_Request request; + take.matrix[m][i].resize(take.count[m][i]); + MPI_Irecv(take.matrix[m][i].data(), take.count[m][i], MPI_INT, i, 1, comm, &request); + requests.emplace_back(request); + } + } + + MPI_Waitall(requests.size(), requests.data(), status); + +#if !defined(NDEBUG) && defined(VERBOSE_OUTPUT) + for (int i = 0; i < num_ranks; ++i) { + if (i != rank) { + std::cout << "[" << rank << "] <- [" << i << "]: ("; + for (int j = 0; j < take.count[m][i]; ++j) { + std::cout << take.matrix[m][i][j]; + if (j < take.count[m][i] - 1) { + std::cout << ", "; + } + } + std::cout << ")" << std::endl; + } + } +#endif + +#if !defined(NDEBUG) + // verification + std::vector received_ids; + for (int i = 0; i < num_ranks; ++i) { + int nents = take.matrix[m][i].size(); + for (int j = 0; j < nents; j++) { + int id = take.matrix[m][i][j]; + if (std::find(received_ids.begin(), received_ids.end(), id) == + received_ids.end()) + received_ids.push_back(id); + } + } + assert(received_ids.size() == static_cast(num_ghosts[rank])); +#endif + MPI_Barrier(comm); + } + } + + /** + * @brief Check if the given entity is a ghost one. + * + * @param i: the index of the entity. + * @return true if a ghost entity, false otherwise. + */ + bool is_ghost(int i) const { + switch (entity) { + case Wonton::CELL: return mesh.cell_get_type(i) == Wonton::PARALLEL_GHOST; + case Wonton::NODE: return mesh.node_get_type(i) == Wonton::PARALLEL_GHOST; + default: return false; + } + } + + /** + * @brief Update the given material field values on ghost cells. + * + * @tparam T: the material field data type. + * @param field: the material field name. + * @param m: the material number. + * @param cache: whether to cache values or not. + */ + template + void update_material_field(std::string const& field, int m, bool cache = false) { + T* values; + state.mat_get_celldata(field, m - 1, &values); + update_material_ghost_values(values, m - 1, cache); + } + + /** + * @brief Update the given mesh field values on ghost cells. + * + * @tparam T: the mesh field data type. + * @param field: the mesh field name. + * @param cache: whether to cache values or not. + */ + template + void update_mesh_field(std::string const& field, bool cache = false) { + T* values; + state.mesh_get_data(entity, field, &values); + update_mesh_ghost_values(values, cache); + } + + /** + * @brief Send owned values and receive ghost values. + * + * @param data: field values array (size = num_owned_in_mesh + num_ghost_in_mesh) + * @param m: material index. + * @param cache: whether to cache values or not. + * + * Note that the data must not be the compact array of material cell + * values; rather it has to be the full mesh cell values array + * + * NOTE: if the inner routine does not know anything about materials + * per se and is just populating mesh sized array, then the comm + * matrices may not need to be material aware. The only savings we + * are getting is the size of the buffers we are communicating but + * not the number of communications. So, we can consider + * communicating the full mesh sized buffer for materials as well + * and doing the map from mesh cell values to mat cell values + * outside the lowest level update_values call using the state + * manager. + */ + void update_ghost_values(double* data, int m, bool cache) { + + // skip if no material data for this rank + if (data == nullptr) { return; } + + if (cache) { + if (send.values.empty() or take.values.empty()) { + send.values.resize(num_mats); + take.values.resize(num_mats); + for (int i = 0; i < num_mats; ++i) { + send.values[i].resize(num_ranks); + take.values[i].resize(num_ranks); + } + } + } + + std::vector buffer[num_ranks]; + std::vector requests; + + // step 1: send field values + for (int i = 0; i < num_ranks; ++i) { + if (i != rank and not send.matrix[m][i].empty()) { + buffer[i].clear(); + for (auto&& j : send.matrix[m][i]) { // 'j' local entity index + assert(not is_ghost(j)); + buffer[i].emplace_back(data[j]); + } + MPI_Request request; + MPI_Isend(buffer[i].data(), buffer[i].size(), MPI_DOUBLE, i, 0, comm, &request); + requests.emplace_back(request); + + if (cache) { + send.values[m][i].resize(buffer[i].size()); + std::copy(buffer[i].begin(), buffer[i].end(), send.values[m][i].begin()); + } + } + } + + // step 2: receive field data + for (int i = 0; i < num_ranks; ++i) { + if (i != rank and take.count[m][i]) { + buffer[i].resize(take.count[m][i]); + MPI_Request request; + MPI_Irecv(buffer[i].data(), take.count[m][i], MPI_DOUBLE, i, 0, comm, &request); + requests.emplace_back(request); + } + } + + // step 3: update state + MPI_Waitall(requests.size(), requests.data(), status); + + for (int i = 0; i < num_ranks; ++i) { + if (i != rank and take.count[m][i]) { + for (int j = 0; j < take.count[m][i]; ++j) { + int const& gid = take.matrix[m][i][j]; + int const& lid = take.lookup[gid]; + data[lid] = buffer[i][j]; + } + if (cache) { + take.values[m][i].resize(take.count[m][i]); + std::copy(buffer[i].begin(), buffer[i].end(), take.values[m][i].begin()); + } + buffer[i].clear(); + } + } + } + + /** + * @brief Update vector field values on ghost cells. + * + * @tparam dim: number of components of each vector. + * @param field: field Vector array (size = num_owned_in_mesh + num_ghost_in_mesh) + * @param m: the material index. + * @param cache: whether to cache values or not (for tests). + * + * Note that the data must not be the compact array of material cell + * values; rather it has to be the full mesh cell values array + + */ + template + void update_ghost_values(Vector* field, int m, bool cache = false) { + + static_assert(0 < dim and dim < 4, "invalid dimension"); + + if (cache) { + if (send.values.empty() or take.values.empty()) { + send.values.resize(num_mats); + take.values.resize(num_mats); + for (int i = 0; i < num_mats; ++i) { + send.values[i].resize(num_ranks); + take.values[i].resize(num_ranks); + } + } + } + + // skip if no material data for this rank + if (field == nullptr) { return; } + + // create a MPI contiguous type for serialization + MPI_Datatype MPI_Vector; + MPI_Type_contiguous(dim, MPI_DOUBLE, &MPI_Vector); + MPI_Type_commit(&MPI_Vector); + + std::vector buffer[num_ranks]; // stride = dim + std::vector requests; + + // step 1: send owned values + for (int i = 0; i < num_ranks; ++i) { + if (i != rank and not send.matrix[m][i].empty()) { + buffer[i].clear(); + send.count[m][i] = send.matrix[m][i].size(); + buffer[i].reserve(dim * send.count[m][i]); + + for (auto&& j : send.matrix[m][i]) { // 'j' local entity index + assert(not is_ghost(j)); + for (int d = 0; d < dim; ++d) { + buffer[i].emplace_back(field[j][d]); + } + } + + MPI_Request request; + MPI_Isend(buffer[i].data(), send.count[m][i], MPI_Vector, i, 0, comm, &request); + requests.emplace_back(request); + + if (cache) { + send.values[m][i].resize(buffer[i].size()); + std::copy(buffer[i].begin(), buffer[i].end(), send.values[m][i].begin()); + } + } + } + + // step 2: receive ghost values + for (int i = 0; i < num_ranks; ++i) { + if (i != rank and take.count[m][i]) { + buffer[i].resize(dim * take.count[m][i]); + MPI_Request request; + MPI_Irecv(buffer[i].data(), take.count[m][i], MPI_Vector, i, 0, comm, &request); + requests.emplace_back(request); + } + } + + // step 3: update vector field + MPI_Waitall(requests.size(), requests.data(), status); + + for (int i = 0; i < num_ranks; ++i) { + if (i != rank and take.count[m][i]) { + for (int j = 0; j < take.count[m][i]; ++j) { + int const& gid = take.matrix[m][i][j]; + int const& lid = take.lookup[gid]; + for (int d = 0; d < dim; ++d) { + field[lid][d] = buffer[i][j * dim + d]; + } + } + + if (cache) { + take.values[m][i].resize(dim * take.count[m][i]); + std::copy(buffer[i].begin(), buffer[i].end(), take.values[m][i].begin()); + } + buffer[i].clear(); + } + } + + // While it seems that a type can be committed repeatedly - what + // happens if someone tries to do a ghost exchange of Vector<2> + // and Vector<3> - say during a 3D to 2D remap. MPI_Vector will + // refer to one thing while the caller might assume something + // else. For this reason, free the user type + MPI_Type_free(&MPI_Vector); + } + + + /** + * @brief Update Point field values on ghost cells (e.g. material centroids) + * + * @tparam dim: number of components of each Point. + * @param field: field Point array (size = num_owned_in_mesh + num_ghost_in_mesh) + * @param m: the material index. + * @param cache: whether to cache values or not (for tests). + * + * Note that the data must not be the compact array of material cell + * values; rather it has to be the full mesh cell values array + + */ + template + void update_ghost_values(Point* field, int m, bool cache = false) { + + static_assert(0 < dim and dim < 4, "invalid dimension"); + + if (cache) { + if (send.values.empty() or take.values.empty()) { + send.values.resize(num_mats); + take.values.resize(num_mats); + for (int i = 0; i < num_mats; ++i) { + send.values[i].resize(num_ranks); + take.values[i].resize(num_ranks); + } + } + } + + // skip if no material data for this rank + if (field == nullptr) { return; } + + // create a MPI contiguous type for serialization + MPI_Datatype MPI_Point; + MPI_Type_contiguous(dim, MPI_DOUBLE, &MPI_Point); + MPI_Type_commit(&MPI_Point); + + std::vector buffer[num_ranks]; // stride = dim + std::vector requests; + + // step 1: send owned values + for (int i = 0; i < num_ranks; ++i) { + if (i != rank and not send.matrix[m][i].empty()) { + buffer[i].clear(); + send.count[m][i] = send.matrix[m][i].size(); + buffer[i].reserve(dim * send.count[m][i]); + + for (auto&& j : send.matrix[m][i]) { // 'j' local entity index + assert(not is_ghost(j)); + for (int d = 0; d < dim; ++d) { + buffer[i].emplace_back(field[j][d]); + } + } + + MPI_Request request; + MPI_Isend(buffer[i].data(), send.count[m][i], MPI_Point, i, 0, comm, &request); + requests.emplace_back(request); + + if (cache) { + send.values[m][i].resize(buffer[i].size()); + std::copy(buffer[i].begin(), buffer[i].end(), send.values[m][i].begin()); + } + } + } + + // step 2: receive ghost values + for (int i = 0; i < num_ranks; ++i) { + if (i != rank and take.count[m][i]) { + buffer[i].resize(dim * take.count[m][i]); + MPI_Request request; + MPI_Irecv(buffer[i].data(), take.count[m][i], MPI_Point, i, 0, comm, &request); + requests.emplace_back(request); + } + } + + // step 3: update vector field + MPI_Waitall(requests.size(), requests.data(), status); + + for (int i = 0; i < num_ranks; ++i) { + if (i != rank and take.count[m][i]) { + for (int j = 0; j < take.count[m][i]; ++j) { + int const& gid = take.matrix[m][i][j]; + int const& lid = take.lookup[gid]; + for (int d = 0; d < dim; ++d) { + field[lid][d] = buffer[i][j * dim + d]; + } + } + + if (cache) { + take.values[m][i].resize(dim * take.count[m][i]); + std::copy(buffer[i].begin(), buffer[i].end(), take.values[m][i].begin()); + } + buffer[i].clear(); + } + } + + // While it seems that a type can be committed repeatedly - what + // happens if someone tries to do a ghost exchange of Vector<2> + // and Vector<3> - say during a 3D to 2D remap. MPI_Vector will + // refer to one thing while the caller might assume something + // else. For this reason, free the user type + MPI_Type_free(&MPI_Point); + } + + /** + * @struct Data for MPI communications. + * + */ + struct Data { + /** lookup table for local indices */ + std::map lookup {}; + /** communication matrices per material */ + std::vector>> matrix {}; + /** exchanged entities count per material */ + std::vector> count {}; + /** cached values per rank for tests, stride = dim for vector fields */ + std::vector>> values {}; + }; + + /** mesh instance */ + Mesh const& mesh; + /** mesh state */ + State& state; + /** MPI communicator */ + MPI_Comm comm = MPI_COMM_NULL; + /** MPI status */ + MPI_Status* status = MPI_STATUS_IGNORE; + /** current MPI rank */ + int rank = 0; + /** number of ranks */ + int num_ranks = 1; + /** number of materials */ + int num_mats = 0; + /** data to send */ + Data send; + /** data to receive */ + Data take; +}; + +} // namespace Wonton +#endif // ifdef WONTON_ENABLE_MPI +#endif // ifndef WONTON_MPI_GHOST_MANAGER_H diff --git a/wonton/distributed/test/test_mpi_ghost_manager.cc b/wonton/distributed/test/test_mpi_ghost_manager.cc new file mode 100644 index 0000000..3474548 --- /dev/null +++ b/wonton/distributed/test/test_mpi_ghost_manager.cc @@ -0,0 +1,482 @@ +/* + * This file is part of the Ristra portage project. + * Please see the license file at the root of this repository, or at: + * https://github.com/laristra/portage/blob/master/LICENSE + */ + +/* NOTE: There are more comprehensive tests for this functionality in Portage */ + +#include "gtest/gtest.h" + +#include "wonton/support/wonton.h" +#include "wonton/distributed/mpi_ghost_manager.h" + +#include "wonton/mesh/jali/jali_mesh_wrapper.h" +#include "wonton/state/jali/jali_state_wrapper.h" + +#include "wonton/intersect/simple_intersect_for_tests.h" + +// Jali headers +#include "Mesh.hh" +#include "MeshFactory.hh" + +/** + * @brief Fixture class for ghost manager tests. + * + */ +class GhostManagerTest : public ::testing::Test { + +protected: + using CellGhostManager = Wonton::MPI_GhostManager; + using NodeGhostManager = Wonton::MPI_GhostManager; + + /** + * @brief Prepare the tests. + * + */ + GhostManagerTest() { + MPI_Comm_rank(comm, &rank); + MPI_Comm_size(comm, &num_ranks); + } + + /** + * @brief Verify send/take communication matrices. + * + * @param ghost_manager: the ghost manager. + * @param mesh: the mesh. + */ + template + void verify_communication_matrices(GhostManager const& ghost_manager, + Mesh const& mesh) const { + + std::vector> send[num_ranks]; /* owned cells that are ghost for some rank */ + std::vector> take[num_ranks]; /* ghost cells that are owned by some rank */ + std::vector num_sent[num_ranks]; /* number of owned cells sent by each rank */ + std::vector num_take[num_ranks]; /* number of ghost cells received from each rank */ + std::vector requests; /* list of asynchronous MPI requests */ + +// GhostManager ghost_manager(mesh, state, comm); + send[rank] = ghost_manager.owned_matrix(); + take[rank] = ghost_manager.ghost_matrix(); + + for (int i = 0; i < num_ranks; ++i) { + num_sent[rank].emplace_back(send[rank][i].size()); + num_take[rank].emplace_back(take[rank][i].size()); + // store GID for owned cells to ease verification + for (int& id : send[rank][i]) { + id = mesh.get_global_id(id, Wonton::CELL); + } + } + + if (rank > 0) { + MPI_Request request; + MPI_Isend(num_sent[rank].data(), num_ranks, MPI_INT, 0, rank, comm, &request); + requests.emplace_back(request); + MPI_Isend(num_take[rank].data(), num_ranks, MPI_INT, 0, rank, comm, &request); + requests.emplace_back(request); + } else { + for (int i = 1; i < num_ranks; ++i) { + MPI_Request request; + num_sent[i].resize(num_ranks, 0); + MPI_Irecv(num_sent[i].data(), num_ranks, MPI_INT, i, i, comm, &request); + requests.emplace_back(request); + num_take[i].resize(num_ranks, 0); + MPI_Irecv(num_take[i].data(), num_ranks, MPI_INT, i, i, comm, &request); + requests.emplace_back(request); + } + } + + MPI_Waitall(requests.size(), requests.data(), MPI_STATUS_IGNORE); + requests.clear(); + + if (rank > 0) { + int const offset = num_ranks * num_ranks; + for (int i = 0; i < num_ranks; ++i) { + if (i != rank) { + int tag = num_ranks * (rank - 1) + i; + MPI_Request request; + MPI_Isend(send[rank][i].data(), send[rank][i].size(), MPI_INT, 0, tag, comm, &request); + requests.emplace_back(request); + MPI_Isend(take[rank][i].data(), take[rank][i].size(), MPI_INT, 0, tag + offset, comm, &request); + requests.emplace_back(request); + } + } + } else { + // check entities count + for (int i = 0; i < num_ranks; ++i) { + for (int j = 0; j < num_ranks; ++j) { + ASSERT_EQ(num_sent[i][j], num_take[j][i]); + } + } + + int const offset = num_ranks * num_ranks; + for (int i = 1; i < num_ranks; ++i) { + send[i].resize(num_ranks); + take[i].resize(num_ranks); + for (int j = 0; j < num_ranks; ++j) { + if (i != j) { + int tag = num_ranks * (i - 1) + j; + MPI_Request request; + send[i][j].resize(num_sent[i][j]); + MPI_Irecv(send[i][j].data(), num_sent[i][j], MPI_INT, i, tag, comm, &request); + requests.emplace_back(request); + take[i][j].resize(num_take[i][j]); + MPI_Irecv(take[i][j].data(), num_take[i][j], MPI_INT, i, tag + offset, comm, &request); + requests.emplace_back(request); + } + } + } + } + + MPI_Waitall(requests.size(), requests.data(), MPI_STATUS_IGNORE); + + // check that we have a perfect matching on + // sent and received cells for each rank pair. + if (rank == 0) { + for (int i = 0; i < num_ranks; ++i) { + for (int j = 0; j < num_ranks; ++j) { + if (i == j) { continue; } + int const num_owned = send[i][j].size(); + int const num_ghost = take[j][i].size(); + ASSERT_EQ(num_owned, num_ghost); + for (int k = 0; k < num_owned; ++k) { + ASSERT_EQ(send[i][j][k], take[j][i][k]); + } + } + } + } + + MPI_Barrier(comm); + } + + /** + * @brief Verify updated ghost values. + * + * @param ghost_manager: the ghost manager + * @param rank: the current MPI rank. + * @param num_ranks: the number of MPI ranks. + * @param comm: the MPI communicator. + */ + template + void verify_values(GhostManager const& ghost_manager) const { + + // gather all sent/received values from all ranks to the master + std::vector> send[num_ranks]; /* values of owned cells that are ghost for some rank */ + std::vector> take[num_ranks]; /* values of ghost cells that are owned by some rank */ + std::vector num_sent[num_ranks]; /* number of owned cells sent by each rank */ + std::vector num_take[num_ranks]; /* number of ghost cells receive from each rank */ + std::vector requests; /* list of asynchronous MPI requests */ + + send[rank] = ghost_manager.owned_values(); + take[rank] = ghost_manager.ghost_values(); + + for (int i = 0; i < num_ranks; ++i) { + num_sent[rank].emplace_back(send[rank][i].size()); + num_take[rank].emplace_back(take[rank][i].size()); + } + + if (rank > 0) { + MPI_Request request; + MPI_Isend(num_sent[rank].data(), num_ranks, MPI_INT, 0, rank, comm, &request); + requests.emplace_back(request); + MPI_Isend(num_take[rank].data(), num_ranks, MPI_INT, 0, rank, comm, &request); + requests.emplace_back(request); + } else { + for (int i = 1; i < num_ranks; ++i) { + MPI_Request request; + num_sent[i].resize(num_ranks, 0); + MPI_Irecv(num_sent[i].data(), num_ranks, MPI_INT, i, i, comm, &request); + requests.emplace_back(request); + num_take[i].resize(num_ranks, 0); + MPI_Irecv(num_take[i].data(), num_ranks, MPI_INT, i, i, comm, &request); + requests.emplace_back(request); + } + } + + MPI_Waitall(requests.size(), requests.data(), MPI_STATUS_IGNORE); + requests.clear(); + + if (rank > 0) { + int const offset = num_ranks * num_ranks; + for (int i = 0; i < num_ranks; ++i) { + if (i != rank) { + int tag = num_ranks * (rank - 1) + i; + MPI_Request request; + MPI_Isend(send[rank][i].data(), send[rank][i].size(), MPI_DOUBLE, 0, tag, comm, &request); + requests.emplace_back(request); + MPI_Isend(take[rank][i].data(), take[rank][i].size(), MPI_DOUBLE, 0, tag + offset, comm, &request); + requests.emplace_back(request); + } + } + } else { + int const offset = num_ranks * num_ranks; + for (int i = 1; i < num_ranks; ++i) { + send[i].resize(num_ranks); + take[i].resize(num_ranks); + for (int j = 0; j < num_ranks; ++j) { + if (i != j) { + int tag = num_ranks * (i - 1) + j; + MPI_Request request; + send[i][j].resize(num_sent[i][j]); + MPI_Irecv(send[i][j].data(), num_sent[i][j], MPI_DOUBLE, i, tag, comm, &request); + requests.emplace_back(request); + take[i][j].resize(num_take[i][j]); + MPI_Irecv(take[i][j].data(), num_take[i][j], MPI_DOUBLE, i, tag + offset, comm, &request); + requests.emplace_back(request); + } + } + } + } + + MPI_Waitall(requests.size(), requests.data(), MPI_STATUS_IGNORE); + + if (rank == 0) { + for (int i = 0; i < num_ranks; ++i) { + for (int j = 0; j < num_ranks; ++j) { + if (i == j) { continue; } + int const num_owned = send[i][j].size(); + int const num_ghost = take[j][i].size(); + ASSERT_EQ(num_owned, num_ghost); + for (int k = 0; k < num_owned; ++k) { + ASSERT_DOUBLE_EQ(send[i][j][k], take[j][i][k]); + } + } + } + } + + MPI_Barrier(comm); + } + + /** current MPI rank */ + int rank = 0; + /** number of MPI ranks */ + int num_ranks = 1; + /** MPI communicator */ + MPI_Comm comm = MPI_COMM_WORLD; + /** number of owned cells */ + int num_owned_cells = 0; + /** number of ghost cells */ + int num_ghost_cells = 0; +}; + +TEST_F(GhostManagerTest, CommMatrices) { + + auto jali_mesh = Jali::MeshFactory(comm)(0.0, 0.0, 1.0, 1.0, 7, 6); + auto jali_state = Jali::State::create(jali_mesh); + Wonton::Jali_Mesh_Wrapper mesh(*jali_mesh); + Wonton::Jali_State_Wrapper state(*jali_state); + + CellGhostManager ghost_manager(mesh, state, comm); + + verify_communication_matrices(ghost_manager, mesh); +} + +TEST_F(GhostManagerTest, UpdateScalarField) { + + auto jali_mesh = Jali::MeshFactory(comm)(0.0, 0.0, 1.0, 1.0, 5, 5); + auto jali_state = Jali::State::create(jali_mesh); + Wonton::Jali_Mesh_Wrapper mesh(*jali_mesh); + Wonton::Jali_State_Wrapper state(*jali_state); + + int const num_owned_cells = mesh.num_owned_cells(); + int const num_ghost_cells = mesh.num_ghost_cells(); + double field[num_owned_cells + num_ghost_cells]; + + for (int i = 0; i < num_owned_cells; ++i) { + Wonton::Point<2> centroid; + mesh.cell_centroid(i, ¢roid); + field[i] = centroid[0] + 2 * centroid[1]; + } + + state.mesh_add_data(Wonton::CELL, "temperature", field); + + CellGhostManager ghost_manager(mesh, state, comm); + ghost_manager.update_ghost_values("temperature", true); + + verify_values(ghost_manager); +} + +TEST_F(GhostManagerTest, UpdateVectorField) { + + auto jali_mesh = Jali::MeshFactory(comm)(0.0, 0.0, 1.0, 1.0, 5, 5); + auto jali_state = Jali::State::create(jali_mesh); + Wonton::Jali_Mesh_Wrapper mesh(*jali_mesh); + Wonton::Jali_State_Wrapper state(*jali_state); + + int const num_owned_nodes = mesh.num_owned_nodes(); + int const num_ghost_nodes = mesh.num_ghost_nodes(); + Wonton::Vector<2> field[num_owned_nodes + num_ghost_nodes]; + + for (int i = 0; i < num_owned_nodes; ++i) { + Wonton::Point<2> coord; + mesh.node_get_coordinates(i, &coord); + field[i][0] = coord[0] + 2 * coord[1]; + field[i][1] = coord[1] + 2 * coord[0]; + } + + state.mesh_add_data>(Wonton::NODE, "velocity", field); + + NodeGhostManager ghost_manager(mesh, state, comm); + ghost_manager.update_ghost_values("velocity", true); + + verify_values(ghost_manager); +} + +TEST_F(GhostManagerTest, UpdateMMScalarField2D) { + + // 0,1 0.5,1 1,1 + // *-------------:------------* + // | : | + // | : 2 | + // | : mat2 | + // | : | + // | : | + // | 0 +............|1,0.5 + // | mat0 : | + // | : | + // | : mat1 | + // | : 1 | + // | : | + // *-------------:------------* + // 0,0 0.5,0 1,0 + + auto jali_mesh = Jali::MeshFactory(comm)(0.0, 0.0, 1.0, 1.0, 5, 5); + auto jali_state = Jali::State::create(jali_mesh); + Wonton::Jali_Mesh_Wrapper mesh(*jali_mesh); + Wonton::Jali_State_Wrapper state(*jali_state); + + int const num_owned_cells = mesh.num_owned_cells(); + int const num_ghost_cells = mesh.num_ghost_cells(); + int const num_cells = num_owned_cells + num_ghost_cells; + + int const num_mats = 3; + std::string const materials[num_mats] = {"mat0", "mat1", "mat2"}; + Wonton::Point<2> const lower_bound[num_mats] = {{0.0, 0.0}, {0.5, 0.0}, {0.5, 0.5}}; + Wonton::Point<2> const upper_bound[num_mats] = {{0.5, 1.0}, {1.0, 0.5}, {1.0, 1.0}}; + + std::vector mat_cells[num_mats]; + std::vector calc_volfracs[num_mats]; + std::vector> calc_centroids[num_mats]; + std::vector calc_vals[num_mats]; + std::vector> calc_vecs[num_mats]; + + // Calculate the cells that go into each material and material field + // values on those cells based on intersection with a known material + // geometry + + for (int c = 0; c < num_cells; c++) { + std::vector> coords; + mesh.cell_get_coordinates(c, &coords); + + double cellvol = mesh.cell_volume(c); + double const min_vol = 1.E-6; + + Wonton::Point<2> box[2]; + BOX_INTERSECT::bounding_box<2>(coords, box, box+1); + + std::vector moments; + for (int m = 0; m < num_mats; m++) { + bool intersected = + BOX_INTERSECT::intersect_boxes<2>(lower_bound[m], upper_bound[m], + box[0], box[1], &moments); + if (intersected and moments[0] > min_vol) { // non-trivial intersection + mat_cells[m].emplace_back(c); + + calc_volfracs[m].emplace_back(moments[0] / cellvol); + + Wonton::Vector<2> p(moments[1] / moments[0], moments[2] / moments[0]); + + calc_centroids[m].emplace_back(p); + calc_vals[m].emplace_back((m + 1) * (m + 1) * (p[0] + p[1])); + calc_vecs[m].emplace_back((m + 1) * p[0], (m + 2) * p[1]); + } + } + } + + // Create a multi-material state and initialize ONLY OWNED cell values + + std::vector state_volfracs[num_mats]; + std::vector> state_centroids[num_mats]; + std::vector state_vals[num_mats]; + std::vector> state_vecs[num_mats]; + + for (int m = 0; m < num_mats; m++) { + state.add_material(materials[m], mat_cells[m]); + int nmatcells = mat_cells[m].size(); + + state_volfracs[m].reserve(nmatcells); + state_centroids[m].reserve(nmatcells); + state_vals[m].reserve(nmatcells); + state_vecs[m].reserve(nmatcells); + for (int i = 0; i < nmatcells; i++) { + if (mesh.cell_get_type(mat_cells[m][i]) == Wonton::PARALLEL_OWNED) { + state_volfracs[m].push_back(calc_volfracs[m][i]); + state_centroids[m].push_back(calc_centroids[m][i]); + state_vals[m].push_back(calc_vals[m][i]); + state_vecs[m].push_back(calc_vecs[m][i]); + } + } + state.mat_add_celldata("matscalar", m, state_vals[m].data()); + state.mat_add_celldata("matvector", m, state_vecs[m].data()); + } + + // Create a ghost update manager + + CellGhostManager ghost_manager(mesh, state, comm); + + + // Update two sets of ghost values using the field name + ghost_manager.update_ghost_values("matscalar"); + ghost_manager.update_ghost_values("matvector"); + + + // Clear out the local copies of the multi-material arrays stored in + // the state manager and reinitialize them + for (int m = 0; m < num_mats; m++) { + int nmatcells = mat_cells[m].size(); + + state_vals[m].clear(); state_vals[m].resize(nmatcells); + double *temp1; + state.mat_get_celldata("matscalar", m, &temp1); + std::copy(temp1, temp1 + nmatcells, state_vals[m].begin()); + + state_vecs[m].clear(); state_vecs[m].resize(nmatcells); + Wonton::Vector<2> *temp2; + state.mat_get_celldata("matvector", m, &temp2); + std::copy(temp2, temp2 + nmatcells, state_vecs[m].begin()); + } + + // Now compare the state values of ALL cells with calculated values; + for (int m = 0; m < num_mats; m++) { + int nmatcells = mat_cells[m].size(); + for (int i = 0; i < nmatcells; i++) { + ASSERT_DOUBLE_EQ(state_vals[m][i], calc_vals[m][i]); + for (int d = 0; d < 2; d++) + ASSERT_DOUBLE_EQ(state_vecs[m][i][d], calc_vecs[m][i][d]); + } + } + + + // Update the other two multi-material arrays not stored in the + // state manager directly (try with and without caching even though + // we won't use caching here + for (int m = 0; m < num_mats; m++) { + ghost_manager.update_material_ghost_values(state_volfracs[m].data(), m); + ghost_manager.update_material_ghost_values(state_centroids[m].data(), m, true); + } + + // Now compare the state values of ALL cells with calculated values; + for (int m = 0; m < num_mats; m++) { + int nmatcells = mat_cells[m].size(); + for (int i = 0; i < nmatcells; i++) { + ASSERT_DOUBLE_EQ(state_volfracs[m][i], calc_volfracs[m][i]); + for (int d = 0; d < 2; d++) + ASSERT_DOUBLE_EQ(state_centroids[m][i][d], calc_centroids[m][i][d]); + } + } + +} diff --git a/wonton/intersect/simple_intersect_for_tests.h b/wonton/intersect/simple_intersect_for_tests.h new file mode 100644 index 0000000..49a078c --- /dev/null +++ b/wonton/intersect/simple_intersect_for_tests.h @@ -0,0 +1,117 @@ +/* +This file is part of the Ristra portage project. +Please see the license file at the root of this repository, or at: + https://github.com/laristra/portage/blob/master/LICENSE +*/ + +#include + +#include "wonton/support/wonton.h" +#include "wonton/support/Point.h" +#include "wonton/support/CoordinateSystem.h" + +namespace BOX_INTERSECT { + +template +void bounding_box(std::vector> coords, + Wonton::Point *pmin, Wonton::Point *pmax) { + *pmin = coords[0]; + *pmax = coords[0]; + for (auto pcoord : coords) { + for (int d = 0; d < D; d++) { + if (pcoord[d] < (*pmin)[d]) + (*pmin)[d] = pcoord[d]; + if (pcoord[d] > (*pmax)[d]) + (*pmax)[d] = pcoord[d]; + } + } +} + +template +bool intersect_boxes(Wonton::Point min1, Wonton::Point max1, + Wonton::Point min2, Wonton::Point max2, + std::vector *xsect_moments) { + + Wonton::Point intmin, intmax; + + for (int d = 0; d < D; ++d) { + // check for non-intersection in this dimension + + if (min1[d] > max2[d]) return false; + if (min2[d] > max1[d]) return false; + + // sort the min max vals in this dimension + double val[4]; + val[0] = min1[d]; val[1] = max1[d]; val[2] = min2[d]; val[3] = max2[d]; + for (int i = 0; i < 3; i++) + for (int j = i+1; j < 4; j++) + if (val[i] > val[j]) { + double tmp = val[i]; + val[i] = val[j]; + val[j] = tmp; + } + + // pick the middle two as the min max coordinates of intersection + // box in this dimension + + intmin[d] = val[1]; intmax[d] = val[2]; + } + + // Calculate the volume + + double vol = 1.0; + for (int d = 0; d < D; d++) + vol *= intmax[d]-intmin[d]; + + // Sanity check + + assert(vol >= 0.0); + + if (vol == 0.0) return false; + + // Calculate the centroid + + Wonton::Point centroid = (intmin+intmax)/2.0; + + // moments + + auto moments = centroid * vol; + CoordSys::modify_volume(vol, intmin, intmax); + CoordSys::modify_first_moments(moments, intmin, intmax); + + xsect_moments->clear(); + xsect_moments->push_back(vol); + for (int d = 0; d < D; d++) + xsect_moments->push_back(moments[d]); + + return true; +} + +template +void intersection_moments(std::vector> cell_xyz, + std::vector>> candidate_cells_xyz, + std::vector *xcells, + std::vector> *xwts) { + + int num_candidates = candidate_cells_xyz.size(); + + xwts->clear(); + + Wonton::Point cmin, cmax; + bounding_box(cell_xyz, &cmin, &cmax); + + for (int c = 0; c < num_candidates; ++c) { + Wonton::Point cmin2, cmax2; + bounding_box(candidate_cells_xyz[c], &cmin2, &cmax2); + + std::vector xsect_moments; + if (intersect_boxes(cmin, cmax, cmin2, cmax2, &xsect_moments)) { + xwts->push_back(xsect_moments); + xcells->push_back(c); + } + } +} + + +} // namespace BOX_INTERSECT + diff --git a/wonton/mesh/AuxMeshTopology.h b/wonton/mesh/AuxMeshTopology.h index 54e588f..ea689d5 100644 --- a/wonton/mesh/AuxMeshTopology.h +++ b/wonton/mesh/AuxMeshTopology.h @@ -94,9 +94,10 @@ void build_sides_3D(AuxMeshTopology& mesh); //! adjacency queries between these entities (In 2D, faces are the //! same as edges and in 1D, faces are the same as nodes). In //! particular, the basic mesh class is expected to support the -//! following methods to successfully instantiate this class: +//! following four groups of methods to successfully instantiate +//! this class: //! -//!~~~ +//!~~~ GROUP I: counts of mesh objects //! int space_dimension() const; // dimensionality of mesh points (1, 2, 3) //! //! int num_owned_cells() const; @@ -105,32 +106,36 @@ void build_sides_3D(AuxMeshTopology& mesh); //! int num_ghost_faces() const; //! int num_owned_nodes() const; //! int num_ghost_nodes() const; +//! +//! +//!~~~ GROUP II: information about mesh objects //! Wonton::Entity_type cell_get_type(int const cellid) const; //! Wonton::Entity_type node_get_type(int const nodeid) const; -//!~~~ //! -//! NOTE: Entity_type is Wonton::OWNED or Wonton::GHOST +//! NOTE: Entity_type is Wonton::PARALLEL_OWNED or Wonton::PARALLEL_GHOST //! -//!~~~ //! Wonton::Element_type cell_get_element_type(int const cellid) const; -//!~~~ //! -//! Can be Wonton::UNKNOWN_TOPOLOGY, Wonton::TRI, Wonton::QUAD, -//! Wonton::POLYGON, Wonton::TET, Wonton::PRISM, Wonton::PYRAMID, -//! Wonton::HEX, Wonton::POLYHEDRON +//! NOTE: Element_type can be Wonton::UNKNOWN_TOPOLOGY, Wonton::TRI, +//! Wonton::QUAD, Wonton::POLYGON, Wonton::TET, Wonton::PRISM, +//! Wonton::PYRAMID, Wonton::HEX, Wonton::POLYHEDRON //! +//! GID_t get_global_id(int const id, Entity_kind const kind) const; //! -//!~~~ +//! NOTE: GID_t == int64_t +//! Entity_kind can be Wonton::CELL, Wonton::SIDE, Wonton::WEDGE, +//! Wonton::FACE, Wonton::EDGE, Wonton::NODE, Wonton::CORNER +//! +//! +//!~~~ GROUP III: upward and downward connectivity //! void cell_get_faces_and_dirs(int const cellid, std::vector *cfaces, //! std::vector *cfdirs) const; -//!~~~ //! //! NOTE: The 'cfdirs' conveys the directions in which the faces are used by -//! the cell. If the natural normal of the face points out of the cell, its -//! direction should be returned as 1, if not, it should be returned as -1 -//! +//! the cell. If the natural normal of the face points out of the cell, +//! its direction should be returned as 1, if not, it should be returned +//! as -1. //! -//!~~~ //! void cell_get_nodes(int const cellid, std::vector *cnodes) const; //! //! void face_get_nodes(int const faceid, std::vector *fnodes) const; @@ -141,13 +146,19 @@ void build_sides_3D(AuxMeshTopology& mesh); //! void node_get_cells(int const nodeid, Wonton::Entity_type etype, //! std::vector *ncells) const; //! +//! +//!~~~ GROUP IV: mesh geometry //! template //! void node_get_coordinates(int const nodeid, Wonton::Point *pp) const; //! -//! GID_t get_global_id(int const id, Entity_kind const kind) const; -//! (where GID_t == int64_t) +//! NOTE: This is the only method that has to be templated on the space +//! dimension. If someone is implementing a mesh wrapper class to work +//! with a fixed space dimension (D=1 or 2 or 3), this method has to +//! be available for all dimensions. It can be empty or throw an exception +//! for a non-relevant dimensions. The other methods can be implemented +//! for a fixed space dimension. +//! //! -//!~~~ //! ******************************** NOTE *********************************** //! //! THIS IS AN INCOMPLETE CLASS DESIGNED TO BE USED IN A 'CRTP' (CURIOUSLY @@ -155,25 +166,29 @@ void build_sides_3D(AuxMeshTopology& mesh); //! PROVIDE A COMPLETE MESH CLASS //! (See https://en.m.wikipedia.org/wiki/Curiously_recurring_template_pattern) //! -//! So, If one is writing a mesh wrapper class called @c MY_MESH_WRAPPER, it is -//! declared like so +//! So, if one is writing a mesh wrapper class called @c MY_MESH_WRAPPER, it is +//! declared like so (see also example in portage/app/fixed_d_app) //! //!~~~ //! class MY_MESH_WRAPPER : public AuxMeshTopology -//! {......} +//! { +//! required methods(); // see four groups above +//! ... +//! optional methods(); // re-implementation of AuxMeshTopology API +//! ... +//! }; //!~~~ //! and it will automatically have the methods of the AuxMeshTopology class. //! -//! If MY_MESH_WRAPPER has equivalent classes for the ones in -//! AuxMeshTopology (possibly because they are more efficient), then -//! the ones from AuxMeshTopology are overridden +//! If MY_MESH_WRAPPER has methods for the ones in AuxMeshTopology (possibly +//! because they are more efficient), then the ones from AuxMeshTopology are +//! overridden. //! //! NOTE THAT THIS CLASS IS NOT DESIGNED TO EVER BE INSTANTIATED DIRECTLY //! //!*************************************************************************** - template class AuxMeshTopology { public: @@ -212,6 +227,9 @@ class AuxMeshTopology { // build_aux_entities() } + // ! Default coordinate system is Cartesian + + virtual Wonton::CoordSysType mesh_get_coordinate_system() const { return Wonton::CoordSysType::Cartesian; } // //! A method expected to be found in the BasicMesh class but defined // //! here as pure virtual to prevent this class from ever being @@ -2315,8 +2333,11 @@ void AuxMeshTopology::compute_cell_moments() { Polytope poly(cpoints, fnmap); int order = 1; - if (std::is_same::value) + if (std::is_same::value || + std::is_same::value) order = 2; + else if (std::is_same::value) + order = 3; auto moments = poly.moments(order); CoordSys::template shift_moments_list(moments); diff --git a/wonton/mesh/CMakeLists.txt b/wonton/mesh/CMakeLists.txt index e50d27e..b21179a 100644 --- a/wonton/mesh/CMakeLists.txt +++ b/wonton/mesh/CMakeLists.txt @@ -53,3 +53,11 @@ install(TARGETS wonton_mesh PUBLIC_HEADER DESTINATION include/wonton/mesh INCLUDES DESTINATION include/wonton/mesh ) + +if (ENABLE_UNIT_TESTS) + wonton_add_unittest(test_mesh_api + SOURCES test/test_mesh_api.cc + LIBRARIES wonton_mesh + POLICY SERIAL + ) +endif () diff --git a/wonton/mesh/simple/simple_mesh.h b/wonton/mesh/simple/simple_mesh.h index 419f28d..1f0a371 100644 --- a/wonton/mesh/simple/simple_mesh.h +++ b/wonton/mesh/simple/simple_mesh.h @@ -24,9 +24,9 @@ namespace Wonton { /*! @class Simple_Mesh "simple_mesh.h" - @brief A very light-weight, serial, 2D/3D Cartesian mesh. + @brief A very light-weight, serial, 1D/2D/3D Cartesian mesh. - A Simple_Mesh is a non-distributed (i.e. serial), 2D/3D, regular Cartesian mesh. + A Simple_Mesh is a non-distributed (i.e. serial), 1D/2D/3D, regular Cartesian mesh. The user need only specify the domain extents and the number of cells per direction, and the mesh class will create all of the connectivity information. The cells are created in row-first order. @@ -41,6 +41,10 @@ namespace Wonton { Likewise, each non-domain-boundary node is connected to 4 cells and each non-domain-boundary face is connected to two cells. + In a 1D mesh, each cell has 2 nodes and 2 faces. + Likewise, each non-domain-boundary node is connected to two cells and each + non-domain-boundary face is connected to two cells. + Nodes and faces on the domain boundary have less connections as there are no ghost mesh entities in Simple_Mesh. */ @@ -143,6 +147,49 @@ Simple_Mesh(double x0, double y0, double z0, faceids_owned_[i] = i; } + /*! + @brief Constructor for creating a serial, 1D Cartesian mesh. + @param[in] x0 The minimum coordinate of the domain. + @param[in] x1 The maximum coordinate of the domain. + @param[in] nx The number of _cells_. + + Connectivity information is automatically built from global IDs. + This mesh class has _zero_ ghost mesh entities. + */ + Simple_Mesh(double x0, double x1, int nx) : + nx_(nx), + x0_(x0), + x1_(x1) { + spacedim_ = 1; + num_cells_ = nx; + num_nodes_ = nx+1; + num_faces_ = nx_+1; + + nodes_per_face_ = 1; + nodes_per_cell_ = 2; + faces_per_cell_ = 2; + cells_per_node_aug_ = 3; + + // Construct the nodal coordinates from extents and number of nodes + build_node_coords_1d(); + + // Build cell <--> node, cell <--> face, face --> node adjacencies + build_cfn_adjacencies_1d(); + + // Build ownership information - no ghosts in Simple Mesh + nodeids_owned_.resize(num_nodes_); + for (int i(0); i < num_nodes_; ++i) + nodeids_owned_[i] = i; + + cellids_owned_.resize(num_cells_); + for (int i(0); i < num_cells_; ++i) + cellids_owned_[i] = i; + + faceids_owned_.resize(num_faces_); + for (int i(0); i < num_faces_; ++i) + faceids_owned_[i] = i; + } + /// Assignment operator (disabled). Simple_Mesh & operator=(const Simple_Mesh &) = delete; @@ -656,6 +703,77 @@ Simple_Mesh(double x0, double y0, double z0, } + /*! + @brief Constructs and stores the node coordinates from the extents and + number of cells per direction passed to the 2D constructor. + */ + void build_node_coords_1d() { + coordinates1d_.clear(); + + double hx = (x1_ - x0_)/nx_; + + for (int ix(0); ix <= nx_; ++ix) { + coordinates1d_.emplace_back(x0_+ix*hx); + } + } + + /* + @brief Builds the cell-face-node adjacency information for 1D mesh. + */ + void build_cfn_adjacencies_1d() { + // downward adjacencies + cell_to_node_.resize(nodes_per_cell_*num_cells_); + cell_to_face_.resize(faces_per_cell_*num_cells_); + cell_face_dirs_.resize(faces_per_cell_*num_cells_); + face_to_node_.resize(nodes_per_face_*num_faces_); + // upward adjacencies + node_to_cell_.resize(cells_per_node_aug_*num_nodes_); + face_to_cell_.resize(2*num_faces_, -1); + + // cell adjacencies + for (int ix = 0; ix < nx_; ++ix) { + auto thisCell = ix; + auto cstart = nodes_per_cell_ * thisCell; + auto fstart = faces_per_cell_ * thisCell; + + // downward connectivity: cell -> node + cell_to_node_[cstart] = ix; + cell_to_node_[cstart+1] = ix + 1; + + // upward connectivity: node -> cell + for (int iix = ix; iix <= ix+1; ++iix) { + auto cnstart = cells_per_node_aug_ * iix; + auto& count = node_to_cell_[cnstart]; + node_to_cell_[cnstart+count+1] = thisCell; + count++; + } + + // downward connectivity cell -> face + cell_to_face_[fstart] = ix; + cell_face_dirs_[fstart] = -1; + + cell_to_face_[fstart+1] = ix + 1; + cell_face_dirs_[fstart+1] = 1; + + // upward connectivity: face -> cell + int cfstart = 2 * ix; + face_to_cell_[cfstart+1] = thisCell; + + cfstart = 2 * (ix + 1); + face_to_cell_[cfstart] = thisCell; + } + + // fixing boundary effects + face_to_cell_[0] = face_to_cell_[1]; + face_to_cell_[1] = -1; + + // face adjacencies + for (int ix = 0; ix <= nx_; ++ix) { + auto nstart = nodes_per_face_ * ix; + face_to_node_[nstart] = ix; + } + } + /// @c Simple_Mesh can be 2D or 3D int spacedim_ ; @@ -668,6 +786,7 @@ Simple_Mesh(double x0, double y0, double z0, /// Node positions. std::vector> coordinates3d_; std::vector> coordinates2d_; + std::vector> coordinates1d_; /// Set in constructor for 2D quads and 3D hexes for now. int nodes_per_face_ ; @@ -721,7 +840,7 @@ Simple_Mesh(double x0, double y0, double z0, // Specializations /*! - @brief Get the 2D/3D coordinates of a specific node as @c Wonton::Point object. + @brief Get the 1D/2D/3D coordinates of a specific node as @c Wonton::Point object. @param[in] nodeid The ID of the node. @param[out] pp The @c Wonton::Point containing the coordinates for node @c nodeid. @@ -742,6 +861,14 @@ void Simple_Mesh::node_get_coordinates<2>(const ID nodeid, *pp = coordinates2d_[nodeid]; } +template<> +inline +void Simple_Mesh::node_get_coordinates<1>(const ID nodeid, + Point<1> *pp) const { + assert(spacedim_ == 1); + *pp = coordinates1d_[nodeid]; +} + template<> inline void Simple_Mesh::node_set_coordinates<3>(const ID nodeid, diff --git a/wonton/mesh/simple/simple_mesh_wrapper.h b/wonton/mesh/simple/simple_mesh_wrapper.h index 7586757..3d43c2c 100644 --- a/wonton/mesh/simple/simple_mesh_wrapper.h +++ b/wonton/mesh/simple/simple_mesh_wrapper.h @@ -59,9 +59,14 @@ class Simple_Mesh_Wrapper : public AuxMeshTopology { CoordSysType coord_sys = CoordSysType::Cartesian) : AuxMeshTopology( request_sides, request_wedges, request_corners), - mesh_(mesh) { + mesh_(mesh), + coord_sys_(coord_sys) { if (coord_sys == CoordSysType::CylindricalAxisymmetric) AuxMeshTopology::build_aux_entities(); + else if (coord_sys == CoordSysType::CylindricalRadial) + AuxMeshTopology::build_aux_entities(); + else if (coord_sys == CoordSysType::SphericalRadial) + AuxMeshTopology::build_aux_entities(); else AuxMeshTopology::build_aux_entities(); } // explicit constructor @@ -78,6 +83,9 @@ class Simple_Mesh_Wrapper : public AuxMeshTopology { ////////////////////////////////////////////////////////////////////// // The following methods are needed somewhere within AuxMeshTopology + /// mesh global information + virtual CoordSysType mesh_get_coordinate_system() const override { return coord_sys_; } + /// The spatial dimension of the mesh. int space_dimension() const { return mesh_.space_dimension(); @@ -227,6 +235,7 @@ class Simple_Mesh_Wrapper : public AuxMeshTopology { private: /// The mesh to wrap. Simple_Mesh const & mesh_; + const CoordSysType coord_sys_; }; // class Simple_Mesh_Wrapper } // namespace Wonton diff --git a/wonton/mesh/test/my_mesh_wrapper.h b/wonton/mesh/test/my_mesh_wrapper.h new file mode 100644 index 0000000..544101e --- /dev/null +++ b/wonton/mesh/test/my_mesh_wrapper.h @@ -0,0 +1,76 @@ +/* +This file is part of the Ristra Wonton project. +Please see the license file at the root of this repository, or at: + https://github.com/laristra/wonton/blob/master/LICENSE +*/ + +#ifndef PORTAGE_FIXED_D_MESH_WRAPPER_H_ +#define PORTAGE_FIXED_D_MESH_WRAPPER_H_ + +#include +#include + +#include "wonton/mesh/AuxMeshTopology.h" +#include "wonton/support/wonton.h" +#include "wonton/support/Point.h" + +/* +* This is compile-only test. It verifies that a fixed dimensional mesh +* wrapper is supported by the code. The actual return numbers do not +* matter but corresponds to a 2x2x2 cubic mesh. +*/ + +namespace Wonton { + +class FixedD_Mesh_Wrapper : public AuxMeshTopology { + public: + explicit FixedD_Mesh_Wrapper() + : AuxMeshTopology(true, true, true) { + AuxMeshTopology::build_aux_entities(); + } + ~FixedD_Mesh_Wrapper() = default; + + // Requid mesh API + // -- counts + int space_dimension() const { return 3; } + + int num_owned_cells() const { return 8; } + int num_owned_faces() const { return 36; } + int num_owned_sides() const { return 192; } + int num_owned_nodes() const { return 27; } + int num_owned_corners() const { return 64; } + + int num_ghost_cells() const { return 0; } + int num_ghost_faces() const { return 0; } + int num_ghost_sides() const { return 0; } + int num_ghost_nodes() const { return 0; } + + Entity_type cell_get_type(int c) const { return Entity_type::PARALLEL_OWNED; } + Entity_type node_get_type(int n) const { return Entity_type::PARALLEL_OWNED; } + Element_type cell_get_element_type(int c) const { return Element_type::HEX; } + + // -- connectivity information + void cell_get_faces_and_dirs(int c, std::vector* faces, std::vector* dirs) const { + faces->resize(6, 0); + dirs->resize(6, 0); + } + void cell_get_nodes(int c, std::vector* nodes) const { nodes->resize(8, 0); } + void face_get_nodes(int f, std::vector* nodes) const { nodes->resize(4, 0); } + void node_get_cells(int n, const Entity_type ptype, std::vector* cells) const { cells->resize(8, 0); } + + GID_t get_global_id(int id, const Entity_kind kind) const { return id; } + + // -- geometry + void cell_get_coordinates(int c, std::vector> *ccoords) const { ccoords->resize(8, Point<3>(0.0, 0.0, 0.0)); } + + // -- functions that are expected to work for all dimensions + template + void node_get_coordinates(int n, Point* pp) const { + if (D != 3) assert(false); + *pp = Wonton::Point(); + } +}; + +} // namespace Wonton + +#endif diff --git a/wonton/mesh/test/test_mesh_api.cc b/wonton/mesh/test/test_mesh_api.cc new file mode 100644 index 0000000..eac77ca --- /dev/null +++ b/wonton/mesh/test/test_mesh_api.cc @@ -0,0 +1,148 @@ +/* +This file is part of the Ristra portage project. +Please see the license file at the root of this repository, or at: + https://github.com/laristra/portage/blob/master/LICENSE +*/ + +#include +#include "gtest/gtest.h" + +// test includes +#include "my_mesh_wrapper.h" + + +TEST(Mesh_API, Fixed_D) { + Wonton::FixedD_Mesh_Wrapper mesh_wrapper; + + // test mesh + int ncells = mesh_wrapper.num_owned_cells(); + int nfaces = mesh_wrapper.num_owned_faces(); + int nsides = mesh_wrapper.num_owned_sides(); + int nwedges = mesh_wrapper.num_owned_wedges(); + int nnodes = mesh_wrapper.num_owned_nodes(); + int ncorners = mesh_wrapper.num_owned_corners(); + + // verify that there is no ghost data + int ncells_ghost = mesh_wrapper.num_ghost_cells(); + int nfaces_ghost = mesh_wrapper.num_ghost_faces(); + int nnodes_ghost = mesh_wrapper.num_ghost_nodes(); + int isum = ncells_ghost + nfaces_ghost + nnodes_ghost; + ASSERT_EQ(isum, 0); + + ASSERT_EQ(ncells, mesh_wrapper.num_entities(Wonton::CELL, Wonton::Entity_type::ALL)); + ASSERT_EQ(nfaces, mesh_wrapper.num_entities(Wonton::FACE, Wonton::Entity_type::ALL)); + ASSERT_EQ(nsides, mesh_wrapper.num_entities(Wonton::SIDE, Wonton::Entity_type::ALL)); + ASSERT_EQ(nnodes, mesh_wrapper.num_entities(Wonton::NODE, Wonton::Entity_type::ALL)); + + // cells: accumulate data into arrays and make small calculations to avoid any + // possibility for code optimization. + std::vector dirs; + Wonton::Point<3> xc, moment; + std::vector>> cell_coords(ncells); + std::vector> cell_cells_f(ncells), cell_cells_n(ncells); + std::vector> cell_sides(ncells), cell_wedges(ncells); + std::vector> cell_faces(ncells), cell_nodes(ncells), cell_corners(ncells); + + std::vector> fpoints; + std::vector, 4>> wpoints; + std::vector> facetpoints; + + for (int c = 0; c < ncells; ++c) { + mesh_wrapper.cell_get_coordinates(c, &(cell_coords[c])); + mesh_wrapper.cell_get_faces_and_dirs(c, &(cell_faces[c]), &dirs); + mesh_wrapper.cell_get_sides(c, &(cell_sides[c])); + mesh_wrapper.cell_get_wedges(c, &(cell_wedges[c])); + mesh_wrapper.cell_get_nodes(c, &(cell_nodes[c])); + mesh_wrapper.cell_get_corners(c, &(cell_corners[c])); + mesh_wrapper.cell_get_node_adj_cells(c, Wonton::Entity_type::ALL, &(cell_cells_n[c])); + mesh_wrapper.cell_get_face_adj_cells(c, Wonton::Entity_type::ALL, &(cell_cells_f[c])); + + double vol = mesh_wrapper.cell_volume(c); + mesh_wrapper.cell_centroid(c, &xc); + moment += xc * vol; + + mesh_wrapper.cell_get_facetization(c, &facetpoints, &fpoints); + mesh_wrapper.wedges_get_coordinates(c, &wpoints); + } + + // faces + std::vector> face_cells(nfaces), face_nodes(nfaces); + for (int f = 0; f < nfaces; ++f) { + mesh_wrapper.face_get_cells(f, Wonton::Entity_type::ALL, &(face_cells[f])); + mesh_wrapper.face_get_nodes(f, &(face_nodes[f])); + } + + // sides: verify that memory was initialized to zero + isum = 0; + double dsum(0.0); + std::vector, 4>> side_coords(nsides); + for (int s = 0; s < nsides; ++s) { + mesh_wrapper.side_get_coordinates(s, &(side_coords[s])); + + int n = mesh_wrapper.side_get_node(s, 0); + int c = mesh_wrapper.side_get_cell(s); + int f = mesh_wrapper.side_get_face(s); + int w = mesh_wrapper.side_get_wedge(s, 0); + int s2 = mesh_wrapper.side_get_opposite_side(s); + isum += n + c + f + w + s2; + + double vol = mesh_wrapper.side_volume(s); + dsum += vol; + } + ASSERT_GE(isum, 0); + ASSERT_LE(fabs(dsum), 1e-20); + + // wedges + isum = 0; + dsum = 0.0; + std::vector, 4>> wedge_coords(nwedges); + + for (int w = 0; w < nwedges; ++w) { + mesh_wrapper.wedge_get_coordinates(w, &(wedge_coords[w])); + + int c = mesh_wrapper.wedge_get_cell(w); + int s = mesh_wrapper.wedge_get_side(w); + int cn = mesh_wrapper.wedge_get_corner(w); + int w2 = mesh_wrapper.wedge_get_opposite_wedge(w); + int w3 = mesh_wrapper.wedge_get_adjacent_wedge(w); + int f = mesh_wrapper.wedge_get_face(w); + int n = mesh_wrapper.wedge_get_node(w); + isum += c + s + cn + w2 + w3 + f + n; + + double vol = mesh_wrapper.wedge_volume(w); + dsum += vol; + } + ASSERT_GE(isum, 0); + ASSERT_LE(fabs(dsum), 1e-20); + + // nodes + std::vector> node_coords(nnodes); + std::vector> node_cells(nnodes), node_wedges(nnodes); + std::vector> node_corners(nnodes), node_nodes(nnodes); + + for (int n = 0; n < nnodes; ++n) { + mesh_wrapper.node_get_coordinates(n, &(node_coords[n])); + mesh_wrapper.node_get_cells(n, Wonton::Entity_type::ALL, &(node_cells[n])); + mesh_wrapper.node_get_cell_adj_nodes(n, Wonton::Entity_type::ALL, &(node_nodes[n])); + mesh_wrapper.node_get_corners(n, Wonton::Entity_type::ALL, &(node_corners[n])); + mesh_wrapper.node_get_wedges(n, Wonton::Entity_type::ALL, &(node_wedges[n])); + } + + // corners + isum = 0; + dsum = 0.0; + std::vector> corner_wedges(ncorners); + + for (int cn = 0; cn < ncorners; ++cn) { + mesh_wrapper.corner_get_wedges(cn, &(corner_wedges[cn])); + int c = mesh_wrapper.corner_get_cell(cn); + int n = mesh_wrapper.corner_get_node(cn); + isum += c + n; + + double vol = mesh_wrapper.corner_volume(cn); + dsum += vol; + } + ASSERT_GE(isum, 0); + ASSERT_LE(fabs(dsum), 1e-20); +} + diff --git a/wonton/support/CoordinateSystem.h b/wonton/support/CoordinateSystem.h index ef9f463..af41fe7 100644 --- a/wonton/support/CoordinateSystem.h +++ b/wonton/support/CoordinateSystem.h @@ -79,7 +79,9 @@ namespace CoordinateSystem { // Typically, we template a class to add functionality. By majority vote, this // has not been done for curvilinier coordinate systems. enum class CoordSysType { Cartesian = 1, - CylindricalAxisymmetric }; + CylindricalRadial, + CylindricalAxisymmetric, + SphericalRadial }; // ============================================================================ /// Cartesian Coordinates diff --git a/wonton/support/Matrix.cc b/wonton/support/Matrix.cc index 986eb22..a98f988 100644 --- a/wonton/support/Matrix.cc +++ b/wonton/support/Matrix.cc @@ -161,6 +161,8 @@ Matrix Matrix::solve(Matrix const& B, std::string method, std::string &error) { if (error != "none") is_singular_ = 2; else is_singular_ = 1; } + else + is_singular_ = (info == 0) ? 1 : 2; } else if (method == "lapack-sysv") { // LAPACK symmetric matrix // check symmetric @@ -235,6 +237,9 @@ Matrix Matrix::solve(Matrix const& B, std::string method, std::string &error) { if (error != "none") is_singular_ = 2; else is_singular_ = 1; } + else + is_singular_ = (info == 0) ? 1 : 2; + } else if (method == "lapack-gesv") { // LAPACK general matrix @@ -308,6 +313,9 @@ Matrix Matrix::solve(Matrix const& B, std::string method, std::string &error) { if (error != "none") is_singular_ = 2; else is_singular_ = 1; } + else + is_singular_ = (info == 0) ? 1 : 2; + } else if (method == "lapack-sytr") { // LAPACK symmetric matrix @@ -364,7 +372,9 @@ Matrix Matrix::solve(Matrix const& B, std::string method, std::string &error) { if (error != "none") is_singular_ = 2; else is_singular_ = 1; } - + else + is_singular_ = (info == 0) ? 1 : 2; + // solve it if (not skip) { info = LAPACKE_dsytrs(matrix_layout,uplo,n,nrhs,a,lda,ipiv,b,ldb); @@ -379,6 +389,8 @@ Matrix Matrix::solve(Matrix const& B, std::string method, std::string &error) { if (error != "none") is_singular_ = 2; else is_singular_ = 1; } + else + is_singular_ = (info == 0) ? 1 : 2; } X = XT.transpose(); diff --git a/wonton/support/Matrix.h b/wonton/support/Matrix.h index ec15915..2b8ef31 100644 --- a/wonton/support/Matrix.h +++ b/wonton/support/Matrix.h @@ -375,6 +375,7 @@ inline Matrix operator*(const Vector& a, const Vector& b) { @param[in] A The system matrix @param[in] b The system right-hand side @param[out] x The solution vector + @return 0 if no errors were encounterted, 1 if A is singular Apparently generic template function definitions don't need "inline" keyword, only the specializations do @@ -382,7 +383,7 @@ inline Matrix operator*(const Vector& a, const Vector& b) { https://stackoverflow.com/questions/1759300/when-should-i-write-the-keyword-inline-for-a-function-method */ template -void solve(const Matrix& A, const Vector& b, Vector& x) { +int solve(const Matrix& A, const Vector& b, Vector& x) { int n = A.rows(); assert(n == A.columns()); assert(D == n); @@ -407,7 +408,9 @@ void solve(const Matrix& A, const Vector& b, Vector& x) { s += pow2(R[j][i]); } norm = sqrt(s + pow2(y[i])); - assert(std::fabs(norm > std::numeric_limits::epsilon())); + if (std::fabs(norm < std::numeric_limits::epsilon())) + return 1; + if (s == 0) continue; @@ -431,7 +434,8 @@ void solve(const Matrix& A, const Vector& b, Vector& x) { R = U*R; } // We need to confirm that the last element is nonzero - assert(std::fabs(R[n - 1][n - 1]) > std::numeric_limits::epsilon()); + if (std::fabs(R[n - 1][n - 1]) < std::numeric_limits::epsilon()) + return 1; Vector QTb = Q.transpose()*b; @@ -442,6 +446,7 @@ void solve(const Matrix& A, const Vector& b, Vector& x) { x[i] = QTb[i]/R[i][i]; } + return 0; } /*! @@ -449,12 +454,19 @@ void solve(const Matrix& A, const Vector& b, Vector& x) { @param[in] A The system matrix @param[in] b The system right-hand side @param[out] x The solution vector + @return 0 if no errors were encounterted, 1 if A is singular */ template<> inline -void solve<1>(const Matrix& A, const Vector<1>& b, Vector<1>& x) { - assert(std::fabs(A[0][0]) > std::numeric_limits::epsilon()); +int solve<1>(const Matrix& A, const Vector<1>& b, Vector<1>& x) { + assert(A.rows() == 1); + assert(A.rows() == A.columns()); + + if (std::fabs(A[0][0]) < std::numeric_limits::epsilon()) + return 1; + x[0] = b[0]/A[0][0]; + return 0; } /*! @@ -462,15 +474,21 @@ void solve<1>(const Matrix& A, const Vector<1>& b, Vector<1>& x) { @param[in] A The system matrix @param[in] b The system right-hand side @param[out] x The solution vector + @return 0 if no errors were encounterted, 1 if A is singular */ template<> inline -void solve<2>(const Matrix& A, const Vector<2>& b, Vector<2>& x) { +int solve<2>(const Matrix& A, const Vector<2>& b, Vector<2>& x) { + assert(A.rows() == 2); + assert(A.rows() == A.columns()); + double detA = A[0][0]*A[1][1] - A[0][1]*A[1][0]; - assert(std::fabs(detA) > std::numeric_limits::epsilon()); + if (std::fabs(detA) < std::numeric_limits::epsilon()) + return 1; x[0] = (A[1][1]*b[0] - A[0][1]*b[1])/detA; x[1] = (A[0][0]*b[1] - A[1][0]*b[0])/detA; + return 0; } } // namespace Wonton diff --git a/wonton/support/Polytope.h b/wonton/support/Polytope.h index 5db5b59..e7320bb 100644 --- a/wonton/support/Polytope.h +++ b/wonton/support/Polytope.h @@ -82,25 +82,6 @@ class Polytope { return face_vertices_[face_id]; } - /*! - @brief Coordinates of vertices of the polytope's face - in counter-clockwise order (i.e. so that the normal - to the face is pointing outside of the polytope) - @param face_id ID of the face of the polytope - @return Vector of coordinates of face's vertices - */ - std::vector< Point > face_points(int const face_id) const { - assert((face_id >= 0) && (unsigned(face_id) < face_vertices_.size())); - - int nvrts = static_cast(face_vertices_[face_id].size()); - std::vector< Point > fpoints; - fpoints.reserve(nvrts); - for (int ivrt = 0; ivrt < nvrts; ivrt++) - fpoints.push_back(vertex_points_[face_vertices_[face_id][ivrt]]); - - return fpoints; - } - /*! @brief Indices of vertices for all the faces of the polytope @return Vector of vectors of indices of faces' vertices @@ -303,7 +284,19 @@ Vector<3> Polytope<3>::face_normal(int face_id) const { template<> inline std::vector Polytope<1>::moments(int order) const { - throw std::runtime_error("Wonton::Polytope<1>::moments(int) not implemented"); + assert(order < 4); + int nmoments = order + 1; + std::vector poly_moments(nmoments, 0.0); + poly_moments[0] = vertex_points_[0][1] - vertex_points_[0][0]; + poly_moments[1] = (vertex_points_[0][1] + vertex_points_[0][0]) * poly_moments[0] / 2; + + if (order > 1) { + double b(vertex_points_[0][1]), a(vertex_points_[0][0]); + double b3(b * b * b), a3(a * a * a); + poly_moments[2] = (b3 - a3) / 3; + if (order > 2) poly_moments[3] = (b3 * b - a3 * a) / 4; + } + return poly_moments; } /*! @@ -359,39 +352,52 @@ std::vector Polytope<3>::moments(int order) const { std::vector poly_moments(4, 0.0); int nfaces = num_faces(); + Point<3> fcentroid; + Vector<3> vcp; + for (int iface = 0; iface < nfaces; iface++) { std::vector fmoments = face_moments(iface); if (fmoments[0] == 0.0) continue; - std::vector< Point<3> > face_pts = face_points(iface); - std::vector< std::vector > itri_pts; - - int nfvrts = face_vertices_[iface].size(); - if (nfvrts == 3) - itri_pts.push_back({0, 1, 2}); - else { - Point<3> fcentroid; - for (int ixyz = 0; ixyz < 3; ixyz++) - fcentroid[ixyz] = fmoments[ixyz + 1]/fmoments[0]; - face_pts.emplace_back(fcentroid); + // Coordinates of vertices of the polytope's face are + // in counter-clockwise order (i.e. so that the normal + // to the face is pointing outside of the polytope) + const auto& ids = face_vertices_[iface]; + int nfvrts = static_cast(ids.size()); - itri_pts.reserve(nfvrts); - for (int ivrt = 0; ivrt < nfvrts; ivrt++) - itri_pts.push_back({nfvrts, ivrt, (ivrt + 1)%nfvrts}); - } + if (nfvrts == 3) { + int v0 = ids[0], v1 = ids[1], v2 = ids[2]; + vcp = cross(vertex_points_[v1] - vertex_points_[v0], + vertex_points_[v2] - vertex_points_[v0]); + + poly_moments[0] += dot(vcp, vertex_points_[v0].asV()); + + for (int i = 0; i < 3; i++) { + double q0 = vertex_points_[v0][i] + vertex_points_[v1][i]; + double q1 = vertex_points_[v0][i] + vertex_points_[v2][i]; + double q2 = vertex_points_[v1][i] + vertex_points_[v2][i]; + poly_moments[i + 1] += vcp[i] * (q0 * q0 + q1 * q1 + q2 * q2); + } + } else { + for (int i = 0; i < 3; i++) + fcentroid[i] = fmoments[i + 1] / fmoments[0]; + + for (int n = 0; n < nfvrts; ++n) { + int v0 = ids[n]; + int v1 = ids[(n + 1)%nfvrts]; + + vcp = cross(vertex_points_[v0] - fcentroid, + vertex_points_[v1] - fcentroid); - for (auto&& point : itri_pts) { - Vector<3> vcp = cross(face_pts[point[1]] - face_pts[point[0]], - face_pts[point[2]] - face_pts[point[0]]); + poly_moments[0] += dot(vcp, vertex_points_[v0].asV()); - poly_moments[0] += dot(vcp, face_pts[point[0]].asV()); - for (int ixyz = 0; ixyz < 3; ixyz++) { - for (int ivrt = 0; ivrt < 3; ivrt++) { - int ifv = point[ivrt], isv = point[(ivrt + 1)%3]; - poly_moments[ixyz + 1] += vcp[ixyz]*pow(face_pts[ifv][ixyz] + - face_pts[isv][ixyz], 2); + for (int i = 0; i < 3; ++i) { + double q0 = vertex_points_[v0][i] + vertex_points_[v1][i]; + double q1 = vertex_points_[v0][i] + fcentroid[i]; + double q2 = vertex_points_[v1][i] + fcentroid[i]; + poly_moments[i + 1] += vcp[i] * (q0 * q0 + q1 * q1 + q2 * q2); } } } diff --git a/wonton/support/wonton.h b/wonton/support/wonton.h index 76eec46..a6c4474 100644 --- a/wonton/support/wonton.h +++ b/wonton/support/wonton.h @@ -238,7 +238,7 @@ struct SerialExecutor_type : Executor_type {}; // for RTTI #ifdef WONTON_ENABLE_MPI struct MPIExecutor_type : Executor_type { - MPIExecutor_type(MPI_Comm comm) : mpicomm(comm) {} + explicit MPIExecutor_type(MPI_Comm comm) : mpicomm(comm) {} MPI_Comm mpicomm = MPI_COMM_WORLD; }; #endif @@ -262,7 +262,8 @@ template inline OutputIterator transform(InputIterator first, InputIterator last, OutputIterator result, UnaryFunction op) { - return thrust::transform(first, last, result, op); + struct thrust::execution_policy exec; + return thrust::transform(exec, first, last, result, op); } template exec; + return thrust::transform(exec, first1, last1, first2, result, op); } template inline void for_each(InputIterator first, InputIterator last, UnaryFunction f) { - thrust::for_each(first, last, f); + struct thrust::execution_policy exec; + thrust::for_each(exec, first, last, f); } #else // no thrust diff --git a/wonton/swarm/swarm.h b/wonton/swarm/swarm.h index 58ec7b3..396fc73 100644 --- a/wonton/swarm/swarm.h +++ b/wonton/swarm/swarm.h @@ -35,6 +35,12 @@ template class Swarm { public: + /** + * @brief Create an empty swarm. + * + */ + Swarm() = default; + /** * @brief Create a particle field from a given list. * @@ -55,8 +61,52 @@ class Swarm { : points_(points), num_local_points_(points_.size()) {} + + /** + * @brief Initialize from a given list of points. + * + * @param points: the given list of points. + */ + inline void init(std::vector> const& points) { + num_local_points_ = points.size(); + points_.resize(num_local_points_); + std::copy(points.begin(), points.end(), points_.begin()); + } #endif + /** + * @brief Initialize from a given list of points. + * + * @param points: the given list of points. + */ + inline void init(Wonton::vector> const& points) { + num_local_points_ = points.size(); + points_.resize(num_local_points_); + std::copy(points.begin(), points.end(), points_.begin()); + } + + /** + * @brief Resize internal list of points. + * + * @param new_size: the new size + */ + inline void resize(int new_size) { + assert(new_size > 0); + num_local_points_ = new_size; + points_.resize(new_size); + } + + /** + * @brief Assign the given point at the given location. + * + * @param id: index of the point. + * @param p: coordinates of the point. + */ + inline void assign(int id, Wonton::Point const& p) { + assert(id >= 0 and id < num_local_points_); + points_[id] = p; + } + /** * @brief Generate a random or uniform particle field. * diff --git a/wonton/swarm/swarm_state.h b/wonton/swarm/swarm_state.h index af87eaf..8428f72 100644 --- a/wonton/swarm/swarm_state.h +++ b/wonton/swarm/swarm_state.h @@ -30,6 +30,12 @@ template class SwarmState { public: + /** + * @brief Create an empty state. + * + */ + SwarmState() = default; + /** * @brief Initialize from a reference swarm. * @@ -46,6 +52,20 @@ class SwarmState { */ explicit SwarmState(int size) : num_local_points_(size) {} + /** + * @brief Initialize from a reference swarm. + * + * @param swarm: the reference swarm. + */ + void init(Swarm const& swarm) { num_local_points_ = swarm.num_owned_particles(); } + + /** + * @brief Initialize with a field size. + * + * @param size: size of each field. + */ + void init(int size) { num_local_points_ = size; } + /** * @brief Create the swarm state from a given mesh state wrapper. * diff --git a/wonton/swarm/test/test_swarm_state.cc b/wonton/swarm/test/test_swarm_state.cc index b6cd87e..177fb11 100644 --- a/wonton/swarm/test/test_swarm_state.cc +++ b/wonton/swarm/test/test_swarm_state.cc @@ -26,6 +26,8 @@ TEST(SwarmState, basic) { int const num_points = 10; Wonton::vector> points(num_points); + Wonton::Swarm<3> swarm; + Wonton::SwarmState<3> state; for (int i = 0; i < num_points; i++) { // point coordinates are not always initialized in order @@ -38,8 +40,8 @@ TEST(SwarmState, basic) { } // create state - Wonton::Swarm<3> swarm(points); - Wonton::SwarmState<3> state(swarm); + swarm.init(points); + state.init(swarm); ASSERT_EQ(state.get_size(), num_points); // create state fields