Sampling support + testing + omnitrace namespace (#19)

* omnitrace namespace

* Kokkos + Lulesh example/tests

* Sampling support + more

- OMNITRACE_BUILD_TESTING option
- sampling support
- pthread_gotcha
- fixes to labels for mpi_gotcha, fork_gotcha, omnitrace_component
- tasking::block_signals, tasking::unblock_signals
- instrumentation mode option in omnitrace exe
- argument option groups in omnitrace exe
- categories in omnitrace settings
- remove TIMEMORY_ prefixed options

* Release workflow updates

* Updated settings printing

* Fixed defaults in README

* Tweak setting defaults in README

* CMake fixes

* cmake-format

* clang-format

* LULESH_USE_MPI OFF

* LULESH_USE_MPI fix

* timemory add_secondary fix

* timemory ambiguous internal namespace fix

* Update timemory submodule

* Handle output path/prefix in omnitrace

- updated timemory
- updated test environment

* sampling + papi fix

* Fix to sampling without PAPI

* Fix for using too many processors in CI

* formatting

* Updated CI

- minor cmake tweaks
- updated timemory submodule

* Updated CI

* Updated CI

* CI + timemory updates

- data race fixes

* CI updates + debug for sampling

* Sampling updates

- moved tasking::{block,unblock}_signals to sampling namespace
- improvements to sampling w.r.t. thread-locality

* Minimum OMNITRACE_THREAD_COUNT of 128

* Handle multiple dims in sampler data

* Configure libunwind support for timemory

* Improved safeguards for sampling

- updated CI
- lulesh runtime-instrument test tweak

* formatting

* CI updates + sampler updates + misc

- fixed stack-buffer-overflow in omnitrace (get_*file_line_info)
- test labels
- steady_clock instead of system_clock in sampler
- update dyninst submodule with upgradePlaceholder fix
- disable OMNITRACE_BUILD_TESTING by default

* Updated timemory submodule

- hidden visibility for timemory
- storage finalizers do not capture this

* Update timemory submodule

- component visibility updates

* Reworked header includes

- use <...> for timemory headers
- always include <library/defines.hpp>

* Rename some config options

* Update PTL submodule

* Update kokkos submodule

* Updated sampling

* Updated CI

* Reworked instrumentation exe

- lowered min-address-range threshold to 256
- extended whole function exclude

* CI fix + timemory submodule update

- TIMEMORY_VISIBLE on component base
- RelWithDebugInfo -> RelWithDebInfo
- Info output for parallel-overhead

* Sampling flags + transpose update + CI update

- disable critical trace for parallel-overhead in CI
- SA_RESTART only in sampler
- reworked transpose example to use fewer threads

* CI update

- removed ubuntu-focal-external-debug
- reduced data artifacts upload

* CI timeouts

- updated timemory submodule
- minor tweaks to omnitrace exe logging

* LICENSE updates (partial)

* CI Test stage timeout extension

* Docker and Packaging updates

* Miscellaneous fixes/tweaks

- gpu.hpp / gpu.cpp
- disable roctracer component if no devices
- re-enable InstrStackFrames by default
- disable sampling by default
- pthread_gotcha::m_enable_sampling is false by default
- timemory submodule update w/ sampler and pop(tid) updates
- fix minor bug in sampler logic
- CMake: OMNITRACE_USE_HIP option
- roctracer + timemory fix

* Replaced OMNITRACE_USE_ROCTRACER with OMNITRACE_USE_HIP where appropriate

* cmake format

* Sampler deadlock fixes

* Removed debug messages from sampler

* Fix for MPI detection + test tweaks + misc

* Sampler deadlock fixes + misc

- removed papi_tot_ins
- pthread_gotcha blocks signals globally until sampler is setup
- metadata specialization for sampling components
- OMNITRACE_INSTRUMENTATION_MODE -> OMNITRACE_MODE
- default sampling delay increased to 0.05 from 1.0e-6
- removed {block,unblock}_signals from critical_trace and ptl
    - no longer necessary to use
- sampling delay minimum is 1.0e-3
- OMNITRACE_BUILD_HIDDEN_VISIBILITY

* omnitrace-avail + libunwind update + restructure

- restructured omnitrace components
- build custom omnitrace-avail executable
- updated libunwind to avoid malloc in get_unw_backtrace

* Fix remaining reorganization issues

- removed some duplicate code
- fixed some trait specializations after implicit instatiation
- formatting

* ensure_storage fix + avail improvements

- fix ensure_storage when component not avail
- suppress irrelevant info in omnitrace-avail

* Delay settings initialization

- slight tweak to tests w/ MPI

* Disable OpenMPI testing w/ ubuntu-bionic

- MPI testing is hanging bc of network interface issue on system:

> [[20462,1],0]: A high-performance Open MPI point-to-point messaging module
> was unable to find any relevant network interfaces:
> Module: OpenFabrics (openib)
>   Host: fv-az19-371
> Another transport will be used instead, although this may result in
> lower performance.
> NOTE: You can disable this warning by setting the MCA parameter
> btl_base_warn_component_unused to 0.
Bu işleme şunda yer alıyor:
Jonathan R. Madsen
2022-01-24 20:49:17 -06:00
işlemeyi yapan: GitHub
ebeveyn 39cf760a4e
işleme 778af2a760
78 değiştirilmiş dosya ile 14296 ekleme ve 1071 silme
+4
Dosyayı Görüntüle
@@ -7,3 +7,7 @@ set(CMAKE_CXX_VISIBILITY_PRESET "default")
add_subdirectory(transpose)
add_subdirectory(parallel-overhead)
option(BUILD_SHARED_LIBS "Build dynamic libraries" ON)
add_subdirectory(lulesh)
+60
Dosyayı Görüntüle
@@ -0,0 +1,60 @@
cmake_minimum_required(VERSION 3.15 FATAL_ERROR)
project(lulesh LANGUAGES C CXX)
list(INSERT CMAKE_MODULE_PATH 0 ${PROJECT_SOURCE_DIR}/cmake/Modules)
add_subdirectory(external)
set(CMAKE_CXX_EXTENSIONS OFF)
if("${CMAKE_BUILD_TYPE}" STREQUAL "")
set(CMAKE_BUILD_TYPE
"RelWithDebInfo"
CACHE STRING "CMake build type" FORCE)
endif()
if(DEFINED OMNITRACE_USE_MPI)
option(LULESH_USE_MPI "Enable MPI" ${OMNITRACE_USE_MPI})
else()
option(LULESH_USE_MPI "Enable MPI" OFF)
endif()
add_library(lulesh-mpi INTERFACE)
if(LULESH_USE_MPI)
find_package(MPI REQUIRED)
target_compile_definitions(lulesh-mpi INTERFACE USE_MPI=1)
target_link_libraries(lulesh-mpi INTERFACE MPI::MPI_C MPI::MPI_CXX)
else()
target_compile_definitions(lulesh-mpi INTERFACE USE_MPI=0)
endif()
if(NOT TARGET Kokkos::kokkos)
find_package(Kokkos REQUIRED)
endif()
file(GLOB headers ${PROJECT_SOURCE_DIR}/*.h ${PROJECT_SOURCE_DIR}/*.hxx)
file(GLOB sources ${PROJECT_SOURCE_DIR}/*.cc)
add_executable(${PROJECT_NAME} ${sources} ${headers})
target_include_directories(${PROJECT_NAME} PRIVATE ${PROJECT_SOURCE_DIR}/includes)
target_link_libraries(${PROJECT_NAME} PRIVATE Kokkos::kokkos lulesh-mpi)
if(NOT CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME)
set_target_properties(${PROJECT_NAME} PROPERTIES RUNTIME_OUTPUT_DIRECTORY
${CMAKE_BINARY_DIR})
endif()
enable_testing()
if(LULESH_USE_MPI)
add_test(
NAME lulesh
COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 8
$<TARGET_FILE:${PROJECT_NAME}> -i 100 -s 20 -p
WORKING_DIRECTORY ${PROJECT_BINARY_DIR})
else()
add_test(
NAME lulesh
COMMAND $<TARGET_FILE:${PROJECT_NAME}> -i 100 -s 20 -p
WORKING_DIRECTORY ${PROJECT_BINARY_DIR})
endif()
+315
Dosyayı Görüntüle
@@ -0,0 +1,315 @@
# include guard
include_guard(DIRECTORY)
# MacroUtilities - useful macros and functions for generic tasks
#
include(CMakeDependentOption)
include(CMakeParseArguments)
# -----------------------------------------------------------------------
# function - capitalize - make a string capitalized (first letter is capital) usage:
# capitalize("SHARED" CShared) message(STATUS "-- CShared is \"${CShared}\"") $ -- CShared
# is "Shared"
function(CAPITALIZE str var)
# make string lower
string(TOLOWER "${str}" str)
string(SUBSTRING "${str}" 0 1 _first)
string(TOUPPER "${_first}" _first)
string(SUBSTRING "${str}" 1 -1 _remainder)
string(CONCAT str "${_first}" "${_remainder}")
set(${var}
"${str}"
PARENT_SCOPE)
endfunction()
# ----------------------------------------------------------------------------------------#
# macro CHECKOUT_GIT_SUBMODULE()
#
# Run "git submodule update" if a file in a submodule does not exist
#
# ARGS: RECURSIVE (option) -- add "--recursive" flag RELATIVE_PATH (one value) --
# typically the relative path to submodule from PROJECT_SOURCE_DIR WORKING_DIRECTORY (one
# value) -- (default: PROJECT_SOURCE_DIR) TEST_FILE (one value) -- file to check for
# (default: CMakeLists.txt) ADDITIONAL_CMDS (many value) -- any addition commands to pass
#
function(CHECKOUT_GIT_SUBMODULE)
# parse args
cmake_parse_arguments(
CHECKOUT "RECURSIVE"
"RELATIVE_PATH;WORKING_DIRECTORY;TEST_FILE;REPO_URL;REPO_BRANCH"
"ADDITIONAL_CMDS" ${ARGN})
if(NOT CHECKOUT_WORKING_DIRECTORY)
set(CHECKOUT_WORKING_DIRECTORY ${PROJECT_SOURCE_DIR})
endif()
if(NOT CHECKOUT_TEST_FILE)
set(CHECKOUT_TEST_FILE "Makefile")
endif()
# default assumption
if(NOT CHECKOUT_REPO_BRANCH)
set(CHECKOUT_REPO_BRANCH "master")
endif()
find_package(Git)
set(_DIR "${CHECKOUT_WORKING_DIRECTORY}/${CHECKOUT_RELATIVE_PATH}")
# ensure the (possibly empty) directory exists
if(NOT EXISTS "${_DIR}")
if(NOT CHECKOUT_REPO_URL)
message(FATAL_ERROR "submodule directory does not exist")
endif()
endif()
# if this file exists --> project has been checked out if not exists --> not been
# checked out
set(_TEST_FILE "${_DIR}/${CHECKOUT_TEST_FILE}")
# assuming a .gitmodules file exists
set(_SUBMODULE "${PROJECT_SOURCE_DIR}/.gitmodules")
set(_TEST_FILE_EXISTS OFF)
if(EXISTS "${_TEST_FILE}" AND NOT IS_DIRECTORY "${_TEST_FILE}")
set(_TEST_FILE_EXISTS ON)
endif()
if(_TEST_FILE_EXISTS)
return()
endif()
find_package(Git REQUIRED)
set(_SUBMODULE_EXISTS OFF)
if(EXISTS "${_SUBMODULE}" AND NOT IS_DIRECTORY "${_SUBMODULE}")
set(_SUBMODULE_EXISTS ON)
endif()
set(_HAS_REPO_URL OFF)
if(NOT "${CHECKOUT_REPO_URL}" STREQUAL "")
set(_HAS_REPO_URL ON)
endif()
# if the module has not been checked out
if(NOT _TEST_FILE_EXISTS AND _SUBMODULE_EXISTS)
# perform the checkout
execute_process(
COMMAND ${GIT_EXECUTABLE} submodule update --init ${_RECURSE}
${CHECKOUT_ADDITIONAL_CMDS} ${CHECKOUT_RELATIVE_PATH}
WORKING_DIRECTORY ${CHECKOUT_WORKING_DIRECTORY}
RESULT_VARIABLE RET)
# check the return code
if(RET GREATER 0)
set(_CMD "${GIT_EXECUTABLE} submodule update --init ${_RECURSE}
${CHECKOUT_ADDITIONAL_CMDS} ${CHECKOUT_RELATIVE_PATH}")
message(STATUS "function(CHECKOUT_GIT_SUBMODULE) failed.")
message(FATAL_ERROR "Command: \"${_CMD}\"")
else()
set(_TEST_FILE_EXISTS ON)
endif()
endif()
if(NOT _TEST_FILE_EXISTS AND _HAS_REPO_URL)
message(
STATUS "Checking out '${CHECKOUT_REPO_URL}' @ '${CHECKOUT_REPO_BRANCH}'...")
# remove the existing directory
if(EXISTS "${_DIR}")
execute_process(COMMAND ${CMAKE_COMMAND} -E remove_directory ${_DIR})
endif()
# perform the checkout
execute_process(
COMMAND
${GIT_EXECUTABLE} clone -b ${CHECKOUT_REPO_BRANCH}
${CHECKOUT_ADDITIONAL_CMDS} ${CHECKOUT_REPO_URL} ${CHECKOUT_RELATIVE_PATH}
WORKING_DIRECTORY ${CHECKOUT_WORKING_DIRECTORY}
RESULT_VARIABLE RET)
# perform the submodule update
if(CHECKOUT_RECURSIVE
AND EXISTS "${_DIR}"
AND IS_DIRECTORY "${_DIR}")
execute_process(
COMMAND ${GIT_EXECUTABLE} submodule update --init ${_RECURSE}
WORKING_DIRECTORY ${_DIR}
RESULT_VARIABLE RET)
endif()
# check the return code
if(RET GREATER 0)
set(_CMD
"${GIT_EXECUTABLE} clone -b ${CHECKOUT_REPO_BRANCH}
${CHECKOUT_ADDITIONAL_CMDS} ${CHECKOUT_REPO_URL} ${CHECKOUT_RELATIVE_PATH}"
)
message(STATUS "function(CHECKOUT_GIT_SUBMODULE) failed.")
message(FATAL_ERROR "Command: \"${_CMD}\"")
else()
set(_TEST_FILE_EXISTS ON)
endif()
endif()
if(NOT EXISTS "${_TEST_FILE}" OR NOT _TEST_FILE_EXISTS)
message(
FATAL_ERROR
"Error checking out submodule: '${CHECKOUT_RELATIVE_PATH}' to '${_DIR}'")
endif()
endfunction()
# ----------------------------------------------------------------------------------------#
# require variable
#
function(CHECK_REQUIRED VAR)
if(NOT DEFINED ${VAR} OR "${${VAR}}" STREQUAL "")
message(FATAL_ERROR "Variable '${VAR}' must be defined and not empty")
endif()
endfunction()
# -----------------------------------------------------------------------
# function add_feature(<NAME> <DOCSTRING>) Add a project feature, whose activation is
# specified by the existence of the variable <NAME>, to the list of enabled/disabled
# features, plus a docstring describing the feature
#
function(ADD_FEATURE _var _description)
set(EXTRA_DESC "")
foreach(currentArg ${ARGN})
if(NOT "${currentArg}" STREQUAL "${_var}" AND NOT "${currentArg}" STREQUAL
"${_description}")
set(EXTRA_DESC "${EXTA_DESC}${currentArg}")
endif()
endforeach()
set_property(GLOBAL APPEND PROPERTY ${PROJECT_NAME}_FEATURES ${_var})
set_property(GLOBAL PROPERTY ${_var}_DESCRIPTION "${_description}${EXTRA_DESC}")
if("CMAKE_DEFINE" IN_LIST ARGN)
set_property(GLOBAL APPEND PROPERTY ${PROJECT_NAME}_CMAKE_DEFINES
"${_var} @${_var}@")
endif()
endfunction()
# ----------------------------------------------------------------------------------------#
# function add_option(<OPTION_NAME> <DOCSRING> <DEFAULT_SETTING> [NO_FEATURE]) Add an
# option and add as a feature if NO_FEATURE is not provided
#
function(ADD_OPTION _NAME _MESSAGE _DEFAULT)
option(${_NAME} "${_MESSAGE}" ${_DEFAULT})
if("NO_FEATURE" IN_LIST ARGN)
mark_as_advanced(${_NAME})
else()
add_feature(${_NAME} "${_MESSAGE}")
endif()
if("ADVANCED" IN_LIST ARGN)
mark_as_advanced(${_NAME})
endif()
endfunction()
# ----------------------------------------------------------------------------------------#
# function print_enabled_features() Print enabled features plus their docstrings.
#
function(PRINT_ENABLED_FEATURES)
set(_basemsg "The following features are defined/enabled (+):")
set(_currentFeatureText "${_basemsg}")
get_property(_features GLOBAL PROPERTY ${PROJECT_NAME}_FEATURES)
if(NOT "${_features}" STREQUAL "")
list(REMOVE_DUPLICATES _features)
list(SORT _features)
endif()
foreach(_feature ${_features})
if(${_feature})
# add feature to text
set(_currentFeatureText "${_currentFeatureText}\n ${_feature}")
# get description
get_property(_desc GLOBAL PROPERTY ${_feature}_DESCRIPTION)
# print description, if not standard ON/OFF, print what is set to
if(_desc)
if(NOT "${${_feature}}" STREQUAL "ON" AND NOT "${${_feature}}" STREQUAL
"TRUE")
set(_currentFeatureText
"${_currentFeatureText}: ${_desc} -- [\"${${_feature}}\"]")
else()
string(REGEX REPLACE "^${PROJECT_NAME}_USE_" "" _feature_tmp
"${_feature}")
string(TOLOWER "${_feature_tmp}" _feature_tmp_l)
capitalize("${_feature_tmp}" _feature_tmp_c)
foreach(_var _feature _feature_tmp _feature_tmp_l _feature_tmp_c)
set(_ver "${${${_var}}_VERSION}")
if(NOT "${_ver}" STREQUAL "")
set(_desc "${_desc} -- [found version ${_ver}]")
break()
endif()
unset(_ver)
endforeach()
set(_currentFeatureText "${_currentFeatureText}: ${_desc}")
endif()
set(_desc NOTFOUND)
endif()
endif()
endforeach()
if(NOT "${_currentFeatureText}" STREQUAL "${_basemsg}")
message(STATUS "${_currentFeatureText}\n")
endif()
endfunction()
# ----------------------------------------------------------------------------------------#
# function print_disabled_features() Print disabled features plus their docstrings.
#
function(PRINT_DISABLED_FEATURES)
set(_basemsg "The following features are NOT defined/enabled (-):")
set(_currentFeatureText "${_basemsg}")
get_property(_features GLOBAL PROPERTY ${PROJECT_NAME}_FEATURES)
if(NOT "${_features}" STREQUAL "")
list(REMOVE_DUPLICATES _features)
list(SORT _features)
endif()
foreach(_feature ${_features})
if(NOT ${_feature})
set(_currentFeatureText "${_currentFeatureText}\n ${_feature}")
get_property(_desc GLOBAL PROPERTY ${_feature}_DESCRIPTION)
if(_desc)
set(_currentFeatureText "${_currentFeatureText}: ${_desc}")
set(_desc NOTFOUND)
endif(_desc)
endif()
endforeach(_feature)
if(NOT "${_currentFeatureText}" STREQUAL "${_basemsg}")
message(STATUS "${_currentFeatureText}\n")
endif()
endfunction()
# ----------------------------------------------------------------------------------------#
# function print_features() Print all features plus their docstrings.
#
function(PRINT_FEATURES)
message(STATUS "")
print_enabled_features()
print_disabled_features()
endfunction()
# ----------------------------------------------------------------------------------------#
# macro ADD_SUBPROJECT() Does a git submodule update + add_subdirectory
#
macro(ADD_SUBPROJECT PACKAGE_NAME)
# parse args
cmake_parse_arguments(PACKAGE "SUBMODULE" "DIRECTORY" "" ${ARGN})
if(NOT PACKAGE_DIRECTORY)
set(PACKAGE_DIRECTORY ${PACKAGE_NAME})
endif()
# if specified in options
if("${PACKAGE_NAME}" IN_LIST PROJECTS)
if(PACKAGE_SUBMODULE)
checkout_git_submodule(RECURSIVE RELATIVE_PATH ${PACKAGE_DIRECTORY})
endif()
if(NOT EXISTS "${PROJECT_SOURCE_DIR}/${PACKAGE_DIRECTORY}/CMakeLists.txt")
message(
STATUS
"Warning! '${PROJECT_SOURCE_DIR}/${PACKAGE_DIRECTORY}/CMakeLists.txt' does not exist!"
)
else()
add_subdirectory(${PACKAGE_DIRECTORY})
endif()
endif()
endmacro()
+28
Dosyayı Görüntüle
@@ -0,0 +1,28 @@
set(Kokkos_ENABLE_SERIAL
ON
CACHE BOOL "Enable Serial")
set(Kokkos_ENABLE_OPENMP
ON
CACHE BOOL "Enable OpenMP")
if(USE_CUDA)
set(Kokkos_ENABLE_CUDA
ON
CACHE BOOL "Enable CUDA")
set(Kokkos_ENABLE_CUDA_UVM
ON
CACHE BOOL "Enable CUDA UVM")
set(Kokkos_ENABLE_CUDA_LAMBDA
ON
CACHE BOOL "Enable CUDA UVM")
set(Kokkos_ENABLE_CUDA_CONSTEXPR
ON
CACHE BOOL "Enable CUDA UVM")
endif()
checkout_git_submodule(
RELATIVE_PATH external/kokkos WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} REPO_URL
https://github.com/kokkos/kokkos.git REPO_BRANCH develop)
set(CMAKE_SKIP_INSTALL_ALL_DEPENDENCY ON)
add_subdirectory(kokkos)
Sağlanmış Alt modül
+1
+127
Dosyayı Görüntüle
@@ -0,0 +1,127 @@
/*!
******************************************************************************
*
* \file
*
* \brief RAJA header file for simple class that can be used to
* time code sections.
*
* \author Rich Hornung, Center for Applied Scientific Computing, LLNL
* \author Jeff Keasler, Applications, Simulations And Quality, LLNL
*
******************************************************************************
*/
#ifndef RAJA_Timer_HXX
#define RAJA_Timer_HXX
#if defined(RAJA_USE_CYCLE)
# include "./cycle.h"
typedef ticks TimeType;
#elif defined(RAJA_USE_CLOCK)
# include <time.h>
typedef clock_t TimeType;
#elif defined(RAJA_USE_GETTIME)
# include <time.h>
typedef timespec TimeType;
#else
# error RAJA_TIMER_TYPE is undefined!
#endif
namespace RAJA
{
/*!
******************************************************************************
*
* \brief Simple timer class to time code sections.
*
******************************************************************************
*/
class Timer
{
public:
#if defined(RAJA_USE_CYCLE) || defined(RAJA_USE_CLOCK)
Timer()
: telapsed(0)
{
;
}
#endif
#if defined(RAJA_USE_GETTIME)
Timer()
: telapsed(0)
, stime_elapsed(0)
, nstime_elapsed(0)
{
;
}
#endif
#if defined(RAJA_USE_CYCLE)
void start() { tstart = getticks(); }
void stop()
{
tstop = getticks();
set_elapsed();
}
long double elapsed() { return static_cast<long double>(telapsed); }
#endif
#if defined(RAJA_USE_CLOCK)
void start() { tstart = clock(); }
void stop()
{
tstop = clock();
set_elapsed();
}
long double elapsed() { return static_cast<long double>(telapsed) / CLOCKS_PER_SEC; }
#endif
#if defined(RAJA_USE_GETTIME)
# if 0
void start() { clock_gettime(CLOCK_REALTIME, &tstart); }
void stop() { clock_gettime(CLOCK_REALTIME, &tstop); set_elapsed(); }
# else
void start() { clock_gettime(CLOCK_MONOTONIC, &tstart); }
void stop()
{
clock_gettime(CLOCK_MONOTONIC, &tstop);
set_elapsed();
}
# endif
long double elapsed() { return (stime_elapsed + nstime_elapsed); }
#endif
private:
TimeType tstart;
TimeType tstop;
long double telapsed;
#if defined(RAJA_USE_CYCLE) || defined(RAJA_USE_CLOCK)
void set_elapsed() { telapsed += (tstop - tstart); }
#elif defined(RAJA_USE_GETTIME)
long double stime_elapsed;
long double nstime_elapsed;
void set_elapsed()
{
stime_elapsed += static_cast<long double>(tstop.tv_sec - tstart.tv_sec);
nstime_elapsed +=
static_cast<long double>(tstop.tv_nsec - tstart.tv_nsec) / 1000000000.0;
}
#endif
};
} // namespace RAJA
#endif // closing endif for header file include guard
+545
Dosyayı Görüntüle
@@ -0,0 +1,545 @@
/*
* Copyright (c) 2003, 2007-8 Matteo Frigo
* Copyright (c) 2003, 2007-8 Massachusetts Institute of Technology
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*
*/
/* machine-dependent cycle counters code. Needs to be inlined. */
/***************************************************************************/
/* To use the cycle counters in your code, simply #include "cycle.h" (this
file), and then use the functions/macros:
ticks getticks(void);
ticks is an opaque typedef defined below, representing the current time.
You extract the elapsed time between two calls to gettick() via:
double elapsed(ticks t1, ticks t0);
which returns a double-precision variable in arbitrary units. You
are not expected to convert this into human units like seconds; it
is intended only for *comparisons* of time intervals.
(In order to use some of the OS-dependent timer routines like
Solaris' gethrtime, you need to paste the autoconf snippet below
into your configure.ac file and #include "config.h" before cycle.h,
or define the relevant macros manually if you are not using autoconf.)
*/
/***************************************************************************/
/* This file uses macros like HAVE_GETHRTIME that are assumed to be
defined according to whether the corresponding function/type/header
is available on your system. The necessary macros are most
conveniently defined if you are using GNU autoconf, via the tests:
dnl ---------------------------------------------------------------------
AC_C_INLINE
AC_HEADER_TIME
AC_CHECK_HEADERS([sys/time.h c_asm.h intrinsics.h mach/mach_time.h])
AC_CHECK_TYPE([hrtime_t],[AC_DEFINE(HAVE_HRTIME_T, 1, [Define to 1 if hrtime_t is
defined in <sys/time.h>])],,[#if HAVE_SYS_TIME_H #include <sys/time.h> #endif])
AC_CHECK_FUNCS([gethrtime read_real_time time_base_to_time clock_gettime
mach_absolute_time])
dnl Cray UNICOS _rtc() (real-time clock) intrinsic
AC_MSG_CHECKING([for _rtc intrinsic])
rtc_ok=yes
AC_TRY_LINK([#ifdef HAVE_INTRINSICS_H
#include <intrinsics.h>
#endif], [_rtc()], [AC_DEFINE(HAVE__RTC,1,[Define if you have the UNICOS _rtc()
intrinsic.])], [rtc_ok=no]) AC_MSG_RESULT($rtc_ok)
dnl ---------------------------------------------------------------------
*/
/***************************************************************************/
#if TIME_WITH_SYS_TIME
# include <sys/time.h>
# include <time.h>
#else
# if HAVE_SYS_TIME_H
# include <sys/time.h>
# else
# include <time.h>
# endif
#endif
#define INLINE_ELAPSED(INL) \
static INL double elapsed(ticks t1, ticks t0) { return (double) t1 - (double) t0; }
/*----------------------------------------------------------------*/
/* Solaris */
#if defined(HAVE_GETHRTIME) && defined(HAVE_HRTIME_T) && !defined(HAVE_TICK_COUNTER)
typedef hrtime_t ticks;
# define getticks gethrtime
INLINE_ELAPSED(inline)
# define HAVE_TICK_COUNTER
#endif
/*----------------------------------------------------------------*/
/* AIX v. 4+ routines to read the real-time clock or time-base register */
#if defined(HAVE_READ_REAL_TIME) && defined(HAVE_TIME_BASE_TO_TIME) && \
!defined(HAVE_TICK_COUNTER)
typedef timebasestruct_t ticks;
static __inline ticks
getticks(void)
{
ticks t;
read_real_time(&t, TIMEBASE_SZ);
return t;
}
static __inline double
elapsed(ticks t1, ticks t0) /* time in nanoseconds */
{
time_base_to_time(&t1, TIMEBASE_SZ);
time_base_to_time(&t0, TIMEBASE_SZ);
return (((double) t1.tb_high - (double) t0.tb_high) * 1.0e9 +
((double) t1.tb_low - (double) t0.tb_low));
}
# define HAVE_TICK_COUNTER
#endif
/*----------------------------------------------------------------*/
/*
* PowerPC ``cycle'' counter using the time base register.
*/
#if((((defined(__GNUC__) && (defined(__powerpc__) || defined(__ppc__))) || \
(defined(__MWERKS__) && defined(macintosh)))) || \
(defined(__IBM_GCC_ASM) && (defined(__powerpc__) || defined(__ppc__)))) && \
!defined(HAVE_TICK_COUNTER)
typedef unsigned long long ticks;
static __inline__ ticks
getticks(void)
{
unsigned int tbl, tbu0, tbu1;
do
{
__asm__ __volatile__("mftbu %0" : "=r"(tbu0));
__asm__ __volatile__("mftb %0" : "=r"(tbl));
__asm__ __volatile__("mftbu %0" : "=r"(tbu1));
} while(tbu0 != tbu1);
return (((unsigned long long) tbu0) << 32) | tbl;
}
INLINE_ELAPSED(__inline__)
# define HAVE_TICK_COUNTER
#endif
/* MacOS/Mach (Darwin) time-base register interface (unlike UpTime,
from Carbon, requires no additional libraries to be linked). */
#if defined(HAVE_MACH_ABSOLUTE_TIME) && defined(HAVE_MACH_MACH_TIME_H) && \
!defined(HAVE_TICK_COUNTER)
# include <mach/mach_time.h>
typedef uint64_t ticks;
# define getticks mach_absolute_time
INLINE_ELAPSED(__inline__)
# define HAVE_TICK_COUNTER
#endif
/*----------------------------------------------------------------*/
/*
* Pentium cycle counter
*/
#if(defined(__GNUC__) || defined(__ICC)) && defined(__i386__) && \
!defined(HAVE_TICK_COUNTER)
typedef unsigned long long ticks;
static __inline__ ticks
getticks(void)
{
ticks ret;
__asm__ __volatile__("rdtsc" : "=A"(ret));
/* no input, nothing else clobbered */
return ret;
}
INLINE_ELAPSED(__inline__)
# define HAVE_TICK_COUNTER
# define TIME_MIN 5000.0 /* unreliable pentium IV cycle counter */
#endif
/* Visual C++ -- thanks to Morten Nissov for his help with this */
#if _MSC_VER >= 1200 && _M_IX86 >= 500 && !defined(HAVE_TICK_COUNTER)
# include <windows.h>
typedef LARGE_INTEGER ticks;
# define RDTSC __asm __emit 0fh __asm __emit 031h /* hack for VC++ 5.0 */
static __inline ticks
getticks(void)
{
ticks retval;
__asm {
RDTSC
mov retval.HighPart, edx
mov retval.LowPart, eax
}
return retval;
}
static __inline double
elapsed(ticks t1, ticks t0)
{
return (double) t1.QuadPart - (double) t0.QuadPart;
}
# define HAVE_TICK_COUNTER
# define TIME_MIN 5000.0 /* unreliable pentium IV cycle counter */
#endif
/*----------------------------------------------------------------*/
/*
* X86-64 cycle counter
*/
#if(defined(__GNUC__) || defined(__ICC) || defined(__SUNPRO_C)) && \
defined(__x86_64__) && !defined(HAVE_TICK_COUNTER)
typedef unsigned long long ticks;
static __inline__ ticks
getticks(void)
{
unsigned a, d;
__asm__ volatile("rdtsc" : "=a"(a), "=d"(d));
return ((ticks) a) | (((ticks) d) << 32);
}
INLINE_ELAPSED(__inline__)
# define HAVE_TICK_COUNTER
#endif
/* PGI compiler, courtesy Cristiano Calonaci, Andrea Tarsi, & Roberto Gori.
NOTE: this code will fail to link unless you use the -Masmkeyword compiler
option (grrr). */
#if defined(__PGI) && defined(__x86_64__) && !defined(HAVE_TICK_COUNTER)
typedef unsigned long long ticks;
static ticks
getticks(void)
{
asm(" rdtsc; shl $0x20,%rdx; mov %eax,%eax; or %rdx,%rax; ");
}
INLINE_ELAPSED(__inline__)
# define HAVE_TICK_COUNTER
#endif
/* Visual C++, courtesy of Dirk Michaelis */
#if _MSC_VER >= 1400 && (defined(_M_AMD64) || defined(_M_X64)) && \
!defined(HAVE_TICK_COUNTER)
# include <intrin.h>
# pragma intrinsic(__rdtsc)
typedef unsigned __int64 ticks;
# define getticks __rdtsc
INLINE_ELAPSED(__inline)
# define HAVE_TICK_COUNTER
#endif
/*----------------------------------------------------------------*/
/*
* IA64 cycle counter
*/
/* intel's icc/ecc compiler */
#if(defined(__EDG_VERSION) || defined(__ECC)) && defined(__ia64__) && \
!defined(HAVE_TICK_COUNTER)
typedef unsigned long ticks;
# include <ia64intrin.h>
static __inline__ ticks
getticks(void)
{
return __getReg(_IA64_REG_AR_ITC);
}
INLINE_ELAPSED(__inline__)
# define HAVE_TICK_COUNTER
#endif
/* gcc */
#if defined(__GNUC__) && defined(__ia64__) && !defined(HAVE_TICK_COUNTER)
typedef unsigned long ticks;
static __inline__ ticks
getticks(void)
{
ticks ret;
__asm__ __volatile__("mov %0=ar.itc" : "=r"(ret));
return ret;
}
INLINE_ELAPSED(__inline__)
# define HAVE_TICK_COUNTER
#endif
/* HP/UX IA64 compiler, courtesy Teresa L. Johnson: */
#if defined(__hpux) && defined(__ia64) && !defined(HAVE_TICK_COUNTER)
# include <machine/sys/inline.h>
typedef unsigned long ticks;
static inline ticks
getticks(void)
{
ticks ret;
ret = _Asm_mov_from_ar(_AREG_ITC);
return ret;
}
INLINE_ELAPSED(inline)
# define HAVE_TICK_COUNTER
#endif
/* Microsoft Visual C++ */
#if defined(_MSC_VER) && defined(_M_IA64) && !defined(HAVE_TICK_COUNTER)
typedef unsigned __int64 ticks;
# ifdef __cplusplus
extern "C"
# endif
ticks
__getReg(int whichReg);
# pragma intrinsic(__getReg)
static __inline ticks
getticks(void)
{
volatile ticks temp;
temp = __getReg(3116);
return temp;
}
INLINE_ELAPSED(inline)
# define HAVE_TICK_COUNTER
#endif
/*----------------------------------------------------------------*/
/*
* PA-RISC cycle counter
*/
#if defined(__hppa__) || defined(__hppa) && !defined(HAVE_TICK_COUNTER)
typedef unsigned long ticks;
# ifdef __GNUC__
static __inline__ ticks
getticks(void)
{
ticks ret;
__asm__ __volatile__("mfctl 16, %0" : "=r"(ret));
/* no input, nothing else clobbered */
return ret;
}
# else
# include <machine/inline.h>
static inline unsigned long
getticks(void)
{
register ticks ret;
_MFCTL(16, ret);
return ret;
}
# endif
INLINE_ELAPSED(inline)
# define HAVE_TICK_COUNTER
#endif
/*----------------------------------------------------------------*/
/* S390, courtesy of James Treacy */
#if defined(__GNUC__) && defined(__s390__) && !defined(HAVE_TICK_COUNTER)
typedef unsigned long long ticks;
static __inline__ ticks
getticks(void)
{
ticks cycles;
__asm__("stck 0(%0)" : : "a"(&(cycles)) : "memory", "cc");
return cycles;
}
INLINE_ELAPSED(__inline__)
# define HAVE_TICK_COUNTER
#endif
/*----------------------------------------------------------------*/
#if defined(__GNUC__) && defined(__alpha__) && !defined(HAVE_TICK_COUNTER)
/*
* The 32-bit cycle counter on alpha overflows pretty quickly,
* unfortunately. A 1GHz machine overflows in 4 seconds.
*/
typedef unsigned int ticks;
static __inline__ ticks
getticks(void)
{
unsigned long cc;
__asm__ __volatile__("rpcc %0" : "=r"(cc));
return (cc & 0xFFFFFFFF);
}
INLINE_ELAPSED(__inline__)
# define HAVE_TICK_COUNTER
#endif
/*----------------------------------------------------------------*/
#if defined(__GNUC__) && defined(__sparc_v9__) && !defined(HAVE_TICK_COUNTER)
typedef unsigned long ticks;
static __inline__ ticks
getticks(void)
{
ticks ret;
__asm__ __volatile__("rd %%tick, %0" : "=r"(ret));
return ret;
}
INLINE_ELAPSED(__inline__)
# define HAVE_TICK_COUNTER
#endif
/*----------------------------------------------------------------*/
#if(defined(__DECC) || defined(__DECCXX)) && defined(__alpha) && \
defined(HAVE_C_ASM_H) && !defined(HAVE_TICK_COUNTER)
# include <c_asm.h>
typedef unsigned int ticks;
static __inline ticks
getticks(void)
{
unsigned long cc;
cc = asm("rpcc %v0");
return (cc & 0xFFFFFFFF);
}
INLINE_ELAPSED(__inline)
# define HAVE_TICK_COUNTER
#endif
/*----------------------------------------------------------------*/
/* SGI/Irix */
#if defined(HAVE_CLOCK_GETTIME) && defined(CLOCK_SGI_CYCLE) && !defined(HAVE_TICK_COUNTER)
typedef struct timespec ticks;
static inline ticks
getticks(void)
{
struct timespec t;
clock_gettime(CLOCK_SGI_CYCLE, &t);
return t;
}
static inline double
elapsed(ticks t1, ticks t0)
{
return ((double) t1.tv_sec - (double) t0.tv_sec) * 1.0E9 +
((double) t1.tv_nsec - (double) t0.tv_nsec);
}
# define HAVE_TICK_COUNTER
#endif
/*----------------------------------------------------------------*/
/* Cray UNICOS _rtc() intrinsic function */
#if defined(HAVE__RTC) && !defined(HAVE_TICK_COUNTER)
# ifdef HAVE_INTRINSICS_H
# include <intrinsics.h>
# endif
typedef long long ticks;
# define getticks _rtc
INLINE_ELAPSED(inline)
# define HAVE_TICK_COUNTER
#endif
/*----------------------------------------------------------------*/
/* MIPS ZBus */
#if HAVE_MIPS_ZBUS_TIMER
# if defined(__mips__) && !defined(HAVE_TICK_COUNTER)
# include <fcntl.h>
# include <sys/mman.h>
# include <unistd.h>
typedef uint64_t ticks;
static inline ticks
getticks(void)
{
static uint64_t* addr = 0;
if(addr == 0)
{
uint32_t rq_addr = 0x10030000;
int fd;
int pgsize;
pgsize = getpagesize();
fd = open("/dev/mem", O_RDONLY | O_SYNC, 0);
if(fd < 0)
{
perror("open");
return NULL;
}
addr = mmap(0, pgsize, PROT_READ, MAP_SHARED, fd, rq_addr);
close(fd);
if(addr == (uint64_t*) -1)
{
perror("mmap");
return NULL;
}
}
return *addr;
}
INLINE_ELAPSED(inline)
# define HAVE_TICK_COUNTER
# endif
#endif /* HAVE_MIPS_ZBUS_TIMER */
Dosya farkı çok büyük olduğundan ihmal edildi Fark Yükle
+886
Dosyayı Görüntüle
@@ -0,0 +1,886 @@
#include <math.h>
#if USE_MPI
# include <mpi.h>
#endif
#if _OPENMP
# include <omp.h>
#endif
#include "lulesh.h"
#include <cstdlib>
#include <limits.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
static KOKKOS_INLINE_FUNCTION Real_t
CalcElemVolume(const Real_t x0, const Real_t x1, const Real_t x2, const Real_t x3,
const Real_t x4, const Real_t x5, const Real_t x6, const Real_t x7,
const Real_t y0, const Real_t y1, const Real_t y2, const Real_t y3,
const Real_t y4, const Real_t y5, const Real_t y6, const Real_t y7,
const Real_t z0, const Real_t z1, const Real_t z2, const Real_t z3,
const Real_t z4, const Real_t z5, const Real_t z6, const Real_t z7)
{
Real_t twelveth = Real_t(1.0) / Real_t(12.0);
Real_t dx61 = x6 - x1;
Real_t dy61 = y6 - y1;
Real_t dz61 = z6 - z1;
Real_t dx70 = x7 - x0;
Real_t dy70 = y7 - y0;
Real_t dz70 = z7 - z0;
Real_t dx63 = x6 - x3;
Real_t dy63 = y6 - y3;
Real_t dz63 = z6 - z3;
Real_t dx20 = x2 - x0;
Real_t dy20 = y2 - y0;
Real_t dz20 = z2 - z0;
Real_t dx50 = x5 - x0;
Real_t dy50 = y5 - y0;
Real_t dz50 = z5 - z0;
Real_t dx64 = x6 - x4;
Real_t dy64 = y6 - y4;
Real_t dz64 = z6 - z4;
Real_t dx31 = x3 - x1;
Real_t dy31 = y3 - y1;
Real_t dz31 = z3 - z1;
Real_t dx72 = x7 - x2;
Real_t dy72 = y7 - y2;
Real_t dz72 = z7 - z2;
Real_t dx43 = x4 - x3;
Real_t dy43 = y4 - y3;
Real_t dz43 = z4 - z3;
Real_t dx57 = x5 - x7;
Real_t dy57 = y5 - y7;
Real_t dz57 = z5 - z7;
Real_t dx14 = x1 - x4;
Real_t dy14 = y1 - y4;
Real_t dz14 = z1 - z4;
Real_t dx25 = x2 - x5;
Real_t dy25 = y2 - y5;
Real_t dz25 = z2 - z5;
#define TRIPLE_PRODUCT(x1, y1, z1, x2, y2, z2, x3, y3, z3) \
((x1) * ((y2) * (z3) - (z2) * (y3)) + (x2) * ((z1) * (y3) - (y1) * (z3)) + \
(x3) * ((y1) * (z2) - (z1) * (y2)))
Real_t volume = TRIPLE_PRODUCT(dx31 + dx72, dx63, dx20, dy31 + dy72, dy63, dy20,
dz31 + dz72, dz63, dz20) +
TRIPLE_PRODUCT(dx43 + dx57, dx64, dx70, dy43 + dy57, dy64, dy70,
dz43 + dz57, dz64, dz70) +
TRIPLE_PRODUCT(dx14 + dx25, dx61, dx50, dy14 + dy25, dy61, dy50,
dz14 + dz25, dz61, dz50);
#undef TRIPLE_PRODUCT
volume *= twelveth;
return volume;
}
/******************************************/
KOKKOS_INLINE_FUNCTION
Real_t
CalcElemVolume(const Real_t x[8], const Real_t y[8], const Real_t z[8])
{
return CalcElemVolume(x[0], x[1], x[2], x[3], x[4], x[5], x[6], x[7], y[0], y[1],
y[2], y[3], y[4], y[5], y[6], y[7], z[0], z[1], z[2], z[3],
z[4], z[5], z[6], z[7]);
}
/////////////////////////////////////////////////////////////////////
Domain::Domain(Int_t numRanks, Index_t colLoc, Index_t rowLoc, Index_t planeLoc,
Index_t nx, int tp, int nr, int balance, Int_t cost)
: m_e_cut(Real_t(1.0e-7))
, m_p_cut(Real_t(1.0e-7))
, m_q_cut(Real_t(1.0e-7))
, m_v_cut(Real_t(1.0e-10))
, m_u_cut(Real_t(1.0e-7))
, m_hgcoef(Real_t(3.0))
, m_ss4o3(Real_t(4.0) / Real_t(3.0))
, m_qstop(Real_t(1.0e+12))
, m_monoq_max_slope(Real_t(1.0))
, m_monoq_limiter_mult(Real_t(2.0))
, m_qlc_monoq(Real_t(0.5))
, m_qqc_monoq(Real_t(2.0) / Real_t(3.0))
, m_qqc(Real_t(2.0))
, m_eosvmax(Real_t(1.0e+9))
, m_eosvmin(Real_t(1.0e-9))
, m_pmin(Real_t(0.))
, m_emin(Real_t(-1.0e+15))
, m_dvovmax(Real_t(0.1))
, m_refdens(Real_t(1.0))
,
//
// set pointers to (potentially) "new'd" arrays to null to
// simplify deallocation.
//
m_regNumList(0)
, m_nodeElemStart(0)
, m_nodeElemCornerList(0)
, m_regElemSize(0)
, m_regElemlist(0)
#if USE_MPI
, commDataSend(0)
, commDataRecv(0)
#endif
{
Index_t edgeElems = nx;
Index_t edgeNodes = edgeElems + 1;
this->cost() = cost;
m_tp = tp;
m_numRanks = numRanks;
///////////////////////////////
// Initialize Sedov Mesh
///////////////////////////////
// construct a uniform box for this processor
m_colLoc = colLoc;
m_rowLoc = rowLoc;
m_planeLoc = planeLoc;
m_sizeX = edgeElems;
m_sizeY = edgeElems;
m_sizeZ = edgeElems;
m_numElem = edgeElems * edgeElems * edgeElems;
m_numNode = edgeNodes * edgeNodes * edgeNodes;
m_regNumList = Allocate<Index_t>(numElem()); // material indexset
// Elem-centered
AllocateElemPersistent(numElem());
// Node-centered
AllocateNodePersistent(numNode());
SetupCommBuffers(edgeNodes);
// Basic Field Initialization
for(Index_t i = 0; i < numElem(); ++i)
{
e(i) = Real_t(0.0);
p(i) = Real_t(0.0);
q(i) = Real_t(0.0);
ss(i) = Real_t(0.0);
}
// Note - v initializes to 1.0, not 0.0!
for(Index_t i = 0; i < numElem(); ++i)
{
v(i) = Real_t(1.0);
}
for(Index_t i = 0; i < numNode(); ++i)
{
xd(i) = Real_t(0.0);
yd(i) = Real_t(0.0);
zd(i) = Real_t(0.0);
}
for(Index_t i = 0; i < numNode(); ++i)
{
xdd(i) = Real_t(0.0);
ydd(i) = Real_t(0.0);
zdd(i) = Real_t(0.0);
}
for(Index_t i = 0; i < numNode(); ++i)
{
nodalMass(i) = Real_t(0.0);
}
BuildMesh(nx, edgeNodes, edgeElems);
#if _OPENMP
SetupThreadSupportStructures();
#else
// These arrays are not used if we're not threaded
m_nodeElemStart = NULL;
m_nodeElemCornerList = NULL;
#endif
// Setup region index sets. For now, these are constant sized
// throughout the run, but could be changed every cycle to
// simulate effects of ALE on the lagrange solver
CreateRegionIndexSets(nr, balance);
// Setup symmetry nodesets
SetupSymmetryPlanes(edgeNodes);
// Setup element connectivities
SetupElementConnectivities(edgeElems);
// Setup symmetry planes and free surface boundary arrays
SetupBoundaryConditions(edgeElems);
// Setup defaults
// These can be changed (requires recompile) if you want to run
// with a fixed timestep, or to a different end time, but it's
// probably easier/better to just run a fixed number of timesteps
// using the -i flag in 2.x
dtfixed() = Real_t(-1.0e-6); // Negative means use courant condition
stoptime() = Real_t(1.0e-2); // *Real_t(edgeElems*tp/45.0) ;
// Initial conditions
deltatimemultlb() = Real_t(1.1);
deltatimemultub() = Real_t(1.2);
dtcourant() = Real_t(1.0e+20);
dthydro() = Real_t(1.0e+20);
dtmax() = Real_t(1.0e-2);
time() = Real_t(0.);
cycle() = Int_t(0);
// initialize field data
for(Index_t i = 0; i < numElem(); ++i)
{
Real_t x_local[8], y_local[8], z_local[8];
Index_t* elemToNode = nodelist(i);
for(Index_t lnode = 0; lnode < 8; ++lnode)
{
Index_t gnode = elemToNode[lnode];
x_local[lnode] = x(gnode);
y_local[lnode] = y(gnode);
z_local[lnode] = z(gnode);
}
// volume calculations
Real_t volume = CalcElemVolume(x_local, y_local, z_local);
volo(i) = volume;
elemMass(i) = volume;
for(Index_t j = 0; j < 8; ++j)
{
Index_t idx = elemToNode[j];
nodalMass(idx) += volume / Real_t(8.0);
}
}
// deposit initial energy
// An energy of 3.948746e+7 is correct for a problem with
// 45 zones along a side - we need to scale it
const Real_t ebase = Real_t(3.948746e+7);
Real_t scale = (nx * m_tp) / Real_t(45.0);
Real_t einit = ebase * scale * scale * scale;
if(m_rowLoc + m_colLoc + m_planeLoc == 0)
{
// Dump into the first zone (which we know is in the corner)
// of the domain that sits at the origin
e(0) = einit;
}
// set initial deltatime base on analytic CFL calculation
deltatime() = (Real_t(.5) * cbrt(volo(0))) / sqrt(Real_t(2.0) * einit);
} // End constructor
////////////////////////////////////////////////////////////////////////////////
Domain::~Domain()
{
/* Release(&m_regNumList);
Release(&m_nodeElemStart);
Release(&m_nodeElemCornerList);
Release(&m_regElemSize);
for (Index_t i=0 ; i<numReg() ; ++i) {
Release(&m_regElemlist[i]);
}
Release(&m_regElemlist);
#if USE_MPI
Release(&commDataSend);
Release(&commDataRecv);
#endif
*/
} // End destructor
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void
Domain::BuildMesh(Int_t nx, Int_t edgeNodes, Int_t edgeElems)
{
Index_t meshEdgeElems = m_tp * nx;
// initialize nodal coordinates
Index_t nidx = 0;
Real_t tz = Real_t(1.125) * Real_t(m_planeLoc * nx) / Real_t(meshEdgeElems);
for(Index_t plane = 0; plane < edgeNodes; ++plane)
{
Real_t ty = Real_t(1.125) * Real_t(m_rowLoc * nx) / Real_t(meshEdgeElems);
for(Index_t row = 0; row < edgeNodes; ++row)
{
Real_t tx = Real_t(1.125) * Real_t(m_colLoc * nx) / Real_t(meshEdgeElems);
for(Index_t col = 0; col < edgeNodes; ++col)
{
x(nidx) = tx;
y(nidx) = ty;
z(nidx) = tz;
++nidx;
// tx += ds ; // may accumulate roundoff...
tx = Real_t(1.125) * Real_t(m_colLoc * nx + col + 1) /
Real_t(meshEdgeElems);
}
// ty += ds ; // may accumulate roundoff...
ty = Real_t(1.125) * Real_t(m_rowLoc * nx + row + 1) / Real_t(meshEdgeElems);
}
// tz += ds ; // may accumulate roundoff...
tz = Real_t(1.125) * Real_t(m_planeLoc * nx + plane + 1) / Real_t(meshEdgeElems);
}
// embed hexehedral elements in nodal point lattice
Index_t zidx = 0;
nidx = 0;
for(Index_t plane = 0; plane < edgeElems; ++plane)
{
for(Index_t row = 0; row < edgeElems; ++row)
{
for(Index_t col = 0; col < edgeElems; ++col)
{
Index_t* localNode = nodelist(zidx);
localNode[0] = nidx;
localNode[1] = nidx + 1;
localNode[2] = nidx + edgeNodes + 1;
localNode[3] = nidx + edgeNodes;
localNode[4] = nidx + edgeNodes * edgeNodes;
localNode[5] = nidx + edgeNodes * edgeNodes + 1;
localNode[6] = nidx + edgeNodes * edgeNodes + edgeNodes + 1;
localNode[7] = nidx + edgeNodes * edgeNodes + edgeNodes;
++zidx;
++nidx;
}
++nidx;
}
nidx += edgeNodes;
}
}
////////////////////////////////////////////////////////////////////////////////
void
Domain::SetupThreadSupportStructures()
{
// set up node-centered indexing of elements
Index_t* nodeElemCount = Allocate<Index_t>(numNode());
for(Index_t i = 0; i < numNode(); ++i)
{
nodeElemCount[i] = 0;
}
for(Index_t i = 0; i < numElem(); ++i)
{
Index_t* nl = nodelist(i);
for(Index_t j = 0; j < 8; ++j)
{
++(nodeElemCount[nl[j]]);
}
}
m_nodeElemStart = Allocate<Index_t>(numNode() + 1);
m_nodeElemStart[0] = 0;
for(Index_t i = 1; i <= numNode(); ++i)
{
m_nodeElemStart[i] = m_nodeElemStart[i - 1] + nodeElemCount[i - 1];
}
m_nodeElemCornerList = Allocate<Index_t>(m_nodeElemStart[numNode()]);
for(Index_t i = 0; i < numNode(); ++i)
{
nodeElemCount[i] = 0;
}
for(Index_t i = 0; i < numElem(); ++i)
{
Index_t* nl = nodelist(i);
for(Index_t j = 0; j < 8; ++j)
{
Index_t m = nl[j];
Index_t k = i * 8 + j;
Index_t offset = m_nodeElemStart[m] + nodeElemCount[m];
m_nodeElemCornerList[offset] = k;
++(nodeElemCount[m]);
}
}
Index_t clSize = m_nodeElemStart[numNode()];
for(Index_t i = 0; i < clSize; ++i)
{
Index_t clv = m_nodeElemCornerList[i];
if((clv < 0) || (clv > numElem() * 8))
{
fprintf(
stderr,
"AllocateNodeElemIndexes(): nodeElemCornerList entry out of range!\n");
#if USE_MPI
MPI_Abort(MPI_COMM_WORLD, -1);
#else
exit(-1);
#endif
}
}
Release<Index_t>(&nodeElemCount);
}
////////////////////////////////////////////////////////////////////////////////
void
Domain::SetupCommBuffers(Int_t edgeNodes)
{
// allocate a buffer large enough for nodal ghost data
Index_t maxEdgeSize = MAX(this->sizeX(), MAX(this->sizeY(), this->sizeZ())) + 1;
m_maxPlaneSize = CACHE_ALIGN_REAL(maxEdgeSize * maxEdgeSize);
m_maxEdgeSize = CACHE_ALIGN_REAL(maxEdgeSize);
// assume communication to 6 neighbors by default
m_rowMin = (m_rowLoc == 0) ? 0 : 1;
m_rowMax = (m_rowLoc == m_tp - 1) ? 0 : 1;
m_colMin = (m_colLoc == 0) ? 0 : 1;
m_colMax = (m_colLoc == m_tp - 1) ? 0 : 1;
m_planeMin = (m_planeLoc == 0) ? 0 : 1;
m_planeMax = (m_planeLoc == m_tp - 1) ? 0 : 1;
#if USE_MPI
// account for face communication
Index_t comBufSize =
(m_rowMin + m_rowMax + m_colMin + m_colMax + m_planeMin + m_planeMax) *
m_maxPlaneSize * MAX_FIELDS_PER_MPI_COMM;
// account for edge communication
comBufSize +=
((m_rowMin & m_colMin) + (m_rowMin & m_planeMin) + (m_colMin & m_planeMin) +
(m_rowMax & m_colMax) + (m_rowMax & m_planeMax) + (m_colMax & m_planeMax) +
(m_rowMax & m_colMin) + (m_rowMin & m_planeMax) + (m_colMin & m_planeMax) +
(m_rowMin & m_colMax) + (m_rowMax & m_planeMin) + (m_colMax & m_planeMin)) *
m_maxEdgeSize * MAX_FIELDS_PER_MPI_COMM;
// account for corner communication
// factor of 16 is so each buffer has its own cache line
comBufSize +=
((m_rowMin & m_colMin & m_planeMin) + (m_rowMin & m_colMin & m_planeMax) +
(m_rowMin & m_colMax & m_planeMin) + (m_rowMin & m_colMax & m_planeMax) +
(m_rowMax & m_colMin & m_planeMin) + (m_rowMax & m_colMin & m_planeMax) +
(m_rowMax & m_colMax & m_planeMin) + (m_rowMax & m_colMax & m_planeMax)) *
CACHE_COHERENCE_PAD_REAL;
this->commDataSend = Allocate<Real_t>(comBufSize);
this->commDataRecv = Allocate<Real_t>(comBufSize);
// prevent floating point exceptions
memset(this->commDataSend, 0, comBufSize * sizeof(Real_t));
memset(this->commDataRecv, 0, comBufSize * sizeof(Real_t));
#endif
// Boundary nodesets
if(m_colLoc == 0)
m_symmX.resize(edgeNodes * edgeNodes);
if(m_rowLoc == 0)
m_symmY.resize(edgeNodes * edgeNodes);
if(m_planeLoc == 0)
m_symmZ.resize(edgeNodes * edgeNodes);
}
////////////////////////////////////////////////////////////////////////////////
void
Domain::CreateRegionIndexSets(Int_t nr, Int_t balance)
{
#if USE_MPI
Index_t myRank;
MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
srand(myRank);
#else
srand(0);
Index_t myRank = 0;
#endif
this->numReg() = nr;
m_regElemSize = Allocate<Index_t>(numReg());
m_regElemlist = Allocate<Index_t*>(numReg());
Index_t nextIndex = 0;
// if we only have one region just fill it
// Fill out the regNumList with material numbers, which are always
// the region index plus one
if(numReg() == 1)
{
while(nextIndex < numElem())
{
this->regNumList(nextIndex) = 1;
nextIndex++;
}
regElemSize(0) = 0;
}
// If we have more than one region distribute the elements.
else
{
Int_t regionNum;
Int_t regionVar;
Int_t lastReg = -1;
Int_t binSize;
Index_t elements;
Index_t runto = 0;
Int_t costDenominator = 0;
Int_t* regBinEnd = Allocate<Int_t>(numReg());
// Determine the relative weights of all the regions. This is based off the -b
// flag. Balance is the value passed into b.
for(Index_t i = 0; i < numReg(); ++i)
{
regElemSize(i) = 0;
costDenominator += pow((i + 1), balance); // Total sum of all regions weights
regBinEnd[i] =
costDenominator; // Chance of hitting a given region is (regBinEnd[i] -
// regBinEdn[i-1])/costDenominator
}
// Until all elements are assigned
while(nextIndex < numElem())
{
// pick the region
regionVar = rand() % costDenominator;
Index_t i = 0;
while(regionVar >= regBinEnd[i])
i++;
// rotate the regions based on MPI rank. Rotation is Rank % NumRegions this
// makes each domain have a different region with the highest representation
regionNum = ((i + myRank) % numReg()) + 1;
// make sure we don't pick the same region twice in a row
while(regionNum == lastReg)
{
regionVar = rand() % costDenominator;
i = 0;
while(regionVar >= regBinEnd[i])
i++;
regionNum = ((i + myRank) % numReg()) + 1;
}
// Pick the bin size of the region and determine the number of elements.
binSize = rand() % 1000;
if(binSize < 773)
{
elements = rand() % 15 + 1;
}
else if(binSize < 937)
{
elements = rand() % 16 + 16;
}
else if(binSize < 970)
{
elements = rand() % 32 + 32;
}
else if(binSize < 974)
{
elements = rand() % 64 + 64;
}
else if(binSize < 978)
{
elements = rand() % 128 + 128;
}
else if(binSize < 981)
{
elements = rand() % 256 + 256;
}
else
elements = rand() % 1537 + 512;
runto = elements + nextIndex;
// Store the elements. If we hit the end before we run out of elements then
// just stop.
while(nextIndex < runto && nextIndex < numElem())
{
this->regNumList(nextIndex) = regionNum;
nextIndex++;
}
lastReg = regionNum;
}
}
// Convert regNumList to region index sets
// First, count size of each region
for(Index_t i = 0; i < numElem(); ++i)
{
int r = this->regNumList(i) - 1; // region index == regnum-1
regElemSize(r)++;
}
// Second, allocate each region index set
for(Index_t i = 0; i < numReg(); ++i)
{
m_regElemlist[i] = Allocate<Int_t>(regElemSize(i));
regElemSize(i) = 0;
}
// Third, fill index sets
for(Index_t i = 0; i < numElem(); ++i)
{
Index_t r = regNumList(i) - 1; // region index == regnum-1
Index_t regndx = regElemSize(r)++; // Note increment
regElemlist(r, regndx) = i;
}
}
/////////////////////////////////////////////////////////////
void
Domain::SetupSymmetryPlanes(Int_t edgeNodes)
{
Index_t nidx = 0;
for(Index_t i = 0; i < edgeNodes; ++i)
{
Index_t planeInc = i * edgeNodes * edgeNodes;
Index_t rowInc = i * edgeNodes;
for(Index_t j = 0; j < edgeNodes; ++j)
{
if(m_planeLoc == 0)
{
m_symmZ[nidx] = rowInc + j;
}
if(m_rowLoc == 0)
{
m_symmY[nidx] = planeInc + j;
}
if(m_colLoc == 0)
{
m_symmX[nidx] = planeInc + j * edgeNodes;
}
++nidx;
}
}
}
/////////////////////////////////////////////////////////////
void
Domain::SetupElementConnectivities(Int_t edgeElems)
{
lxim(0) = 0;
for(Index_t i = 1; i < numElem(); ++i)
{
lxim(i) = i - 1;
lxip(i - 1) = i;
}
lxip(numElem() - 1) = numElem() - 1;
for(Index_t i = 0; i < edgeElems; ++i)
{
letam(i) = i;
letap(numElem() - edgeElems + i) = numElem() - edgeElems + i;
}
for(Index_t i = edgeElems; i < numElem(); ++i)
{
letam(i) = i - edgeElems;
letap(i - edgeElems) = i;
}
for(Index_t i = 0; i < edgeElems * edgeElems; ++i)
{
lzetam(i) = i;
lzetap(numElem() - edgeElems * edgeElems + i) =
numElem() - edgeElems * edgeElems + i;
}
for(Index_t i = edgeElems * edgeElems; i < numElem(); ++i)
{
lzetam(i) = i - edgeElems * edgeElems;
lzetap(i - edgeElems * edgeElems) = i;
}
}
/////////////////////////////////////////////////////////////
void
Domain::SetupBoundaryConditions(Int_t edgeElems)
{
Index_t ghostIdx[6]; // offsets to ghost locations
// set up boundary condition information
for(Index_t i = 0; i < numElem(); ++i)
{
elemBC(i) = Int_t(0);
}
for(Index_t i = 0; i < 6; ++i)
{
ghostIdx[i] = INT_MIN;
}
Int_t pidx = numElem();
if(m_planeMin != 0)
{
ghostIdx[0] = pidx;
pidx += sizeX() * sizeY();
}
if(m_planeMax != 0)
{
ghostIdx[1] = pidx;
pidx += sizeX() * sizeY();
}
if(m_rowMin != 0)
{
ghostIdx[2] = pidx;
pidx += sizeX() * sizeZ();
}
if(m_rowMax != 0)
{
ghostIdx[3] = pidx;
pidx += sizeX() * sizeZ();
}
if(m_colMin != 0)
{
ghostIdx[4] = pidx;
pidx += sizeY() * sizeZ();
}
if(m_colMax != 0)
{
ghostIdx[5] = pidx;
}
// symmetry plane or free surface BCs
for(Index_t i = 0; i < edgeElems; ++i)
{
Index_t planeInc = i * edgeElems * edgeElems;
Index_t rowInc = i * edgeElems;
for(Index_t j = 0; j < edgeElems; ++j)
{
if(m_planeLoc == 0)
{
elemBC(rowInc + j) |= ZETA_M_SYMM;
}
else
{
elemBC(rowInc + j) |= ZETA_M_COMM;
lzetam(rowInc + j) = ghostIdx[0] + rowInc + j;
}
if(m_planeLoc == m_tp - 1)
{
elemBC(rowInc + j + numElem() - edgeElems * edgeElems) |= ZETA_P_FREE;
}
else
{
elemBC(rowInc + j + numElem() - edgeElems * edgeElems) |= ZETA_P_COMM;
lzetap(rowInc + j + numElem() - edgeElems * edgeElems) =
ghostIdx[1] + rowInc + j;
}
if(m_rowLoc == 0)
{
elemBC(planeInc + j) |= ETA_M_SYMM;
}
else
{
elemBC(planeInc + j) |= ETA_M_COMM;
letam(planeInc + j) = ghostIdx[2] + rowInc + j;
}
if(m_rowLoc == m_tp - 1)
{
elemBC(planeInc + j + edgeElems * edgeElems - edgeElems) |= ETA_P_FREE;
}
else
{
elemBC(planeInc + j + edgeElems * edgeElems - edgeElems) |= ETA_P_COMM;
letap(planeInc + j + edgeElems * edgeElems - edgeElems) =
ghostIdx[3] + rowInc + j;
}
if(m_colLoc == 0)
{
elemBC(planeInc + j * edgeElems) |= XI_M_SYMM;
}
else
{
elemBC(planeInc + j * edgeElems) |= XI_M_COMM;
lxim(planeInc + j * edgeElems) = ghostIdx[4] + rowInc + j;
}
if(m_colLoc == m_tp - 1)
{
elemBC(planeInc + j * edgeElems + edgeElems - 1) |= XI_P_FREE;
}
else
{
elemBC(planeInc + j * edgeElems + edgeElems - 1) |= XI_P_COMM;
lxip(planeInc + j * edgeElems + edgeElems - 1) = ghostIdx[5] + rowInc + j;
}
}
}
}
///////////////////////////////////////////////////////////////////////////
void
InitMeshDecomp(Int_t numRanks, Int_t myRank, Int_t* col, Int_t* row, Int_t* plane,
Int_t* side)
{
Int_t testProcs;
Int_t dx, dy, dz;
Int_t myDom;
// Assume cube processor layout for now
testProcs = Int_t(cbrt(Real_t(numRanks)) + 0.5);
if(testProcs * testProcs * testProcs != numRanks)
{
printf("Num processors must be a cube of an integer (1, 8, 27, ...)\n");
#if USE_MPI
MPI_Abort(MPI_COMM_WORLD, -1);
#else
exit(-1);
#endif
}
if(sizeof(Real_t) != 4 && sizeof(Real_t) != 8)
{
printf("MPI operations only support float and double right now...\n");
#if USE_MPI
MPI_Abort(MPI_COMM_WORLD, -1);
#else
exit(-1);
#endif
}
if(MAX_FIELDS_PER_MPI_COMM > CACHE_COHERENCE_PAD_REAL)
{
printf("corner element comm buffers too small. Fix code.\n");
#if USE_MPI
MPI_Abort(MPI_COMM_WORLD, -1);
#else
exit(-1);
#endif
}
dx = testProcs;
dy = testProcs;
dz = testProcs;
// temporary test
if(dx * dy * dz != numRanks)
{
printf("error -- must have as many domains as procs\n");
#if USE_MPI
MPI_Abort(MPI_COMM_WORLD, -1);
#else
exit(-1);
#endif
}
Int_t remainder = dx * dy * dz % numRanks;
if(myRank < remainder)
{
myDom = myRank * (1 + (dx * dy * dz / numRanks));
}
else
{
myDom = remainder * (1 + (dx * dy * dz / numRanks)) +
(myRank - remainder) * (dx * dy * dz / numRanks);
}
*col = myDom % dx;
*row = (myDom / dx) % dy;
*plane = myDom / (dx * dy);
*side = testProcs;
return;
}
+273
Dosyayı Görüntüle
@@ -0,0 +1,273 @@
#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#if USE_MPI
# include <mpi.h>
#endif
#include "lulesh.h"
/* Helper function for converting strings to ints, with error checking */
int
StrToInt(const char* token, int* retVal)
{
const char* c;
char* endptr;
const int decimal_base = 10;
if(token == NULL)
return 0;
c = token;
*retVal = (int) strtol(c, &endptr, decimal_base);
if((endptr != c) && ((*endptr == ' ') || (*endptr == '\0')))
return 1;
else
return 0;
}
static void
PrintCommandLineOptions(char* execname, int myRank)
{
if(myRank == 0)
{
printf("Usage: %s [opts]\n", execname);
printf(" where [opts] is one or more of:\n");
printf(" -q : quiet mode - suppress all stdout\n");
printf(" -i <iterations> : number of cycles to run\n");
printf(" -s <size> : length of cube mesh along side\n");
printf(" -r <numregions> : Number of distinct regions (def: 11)\n");
printf(" -b <balance> : Load balance between regions of a domain (def: 1)\n");
printf(" -c <cost> : Extra cost of more expensive regions (def: 1)\n");
printf(" -f <numfiles> : Number of files to split viz dump into (def: "
"(np+10)/9)\n");
printf(" -p : Print out progress\n");
printf(
" -v : Output viz file (requires compiling with -DVIZ_MESH\n");
printf(" -h : This message\n");
printf("\n\n");
}
}
static void
ParseError(const char* message, int myRank)
{
if(myRank == 0)
{
printf("%s\n", message);
#if USE_MPI
MPI_Abort(MPI_COMM_WORLD, -1);
#else
exit(-1);
#endif
}
}
void
ParseCommandLineOptions(int argc, char* argv[], int myRank, struct cmdLineOpts* opts)
{
if(argc > 1)
{
int i = 1;
while(i < argc)
{
int ok;
/* -i <iterations> */
if(strcmp(argv[i], "-i") == 0)
{
if(i + 1 >= argc)
{
ParseError("Missing integer argument to -i", myRank);
}
ok = StrToInt(argv[i + 1], &(opts->its));
if(!ok)
{
ParseError("Parse Error on option -i integer value required after "
"argument\n",
myRank);
}
i += 2;
}
/* -s <size, sidelength> */
else if(strcmp(argv[i], "-s") == 0)
{
if(i + 1 >= argc)
{
ParseError("Missing integer argument to -s\n", myRank);
}
ok = StrToInt(argv[i + 1], &(opts->nx));
if(!ok)
{
ParseError("Parse Error on option -s integer value required after "
"argument\n",
myRank);
}
i += 2;
}
/* -r <numregions> */
else if(strcmp(argv[i], "-r") == 0)
{
if(i + 1 >= argc)
{
ParseError("Missing integer argument to -r\n", myRank);
}
ok = StrToInt(argv[i + 1], &(opts->numReg));
if(!ok)
{
ParseError("Parse Error on option -r integer value required after "
"argument\n",
myRank);
}
i += 2;
}
/* -f <numfilepieces> */
else if(strcmp(argv[i], "-f") == 0)
{
if(i + 1 >= argc)
{
ParseError("Missing integer argument to -f\n", myRank);
}
ok = StrToInt(argv[i + 1], &(opts->numFiles));
if(!ok)
{
ParseError("Parse Error on option -f integer value required after "
"argument\n",
myRank);
}
i += 2;
}
/* -p */
else if(strcmp(argv[i], "-p") == 0)
{
opts->showProg = 1;
i++;
}
/* -q */
else if(strcmp(argv[i], "-q") == 0)
{
opts->quiet = 1;
i++;
}
/* -q */
else if(strcmp(argv[i], "-a") == 0)
{
opts->do_atomic = 1;
i++;
}
else if(strcmp(argv[i], "-b") == 0)
{
if(i + 1 >= argc)
{
ParseError("Missing integer argument to -b\n", myRank);
}
ok = StrToInt(argv[i + 1], &(opts->balance));
if(!ok)
{
ParseError("Parse Error on option -b integer value required after "
"argument\n",
myRank);
}
i += 2;
}
else if(strcmp(argv[i], "-c") == 0)
{
if(i + 1 >= argc)
{
ParseError("Missing integer argument to -c\n", myRank);
}
ok = StrToInt(argv[i + 1], &(opts->cost));
if(!ok)
{
ParseError("Parse Error on option -c integer value required after "
"argument\n",
myRank);
}
i += 2;
}
/* -v */
else if(strcmp(argv[i], "-v") == 0)
{
#if VIZ_MESH
opts->viz = 1;
#else
ParseError("Use of -v requires compiling with -DVIZ_MESH\n", myRank);
#endif
i++;
}
/* -h */
else if(strcmp(argv[i], "-h") == 0)
{
PrintCommandLineOptions(argv[0], myRank);
#if USE_MPI
MPI_Abort(MPI_COMM_WORLD, 0);
#else
exit(0);
#endif
}
else
{
char msg[80];
PrintCommandLineOptions(argv[0], myRank);
sprintf(msg, "ERROR: Unknown command line argument: %s\n", argv[i]);
ParseError(msg, myRank);
}
}
}
}
/////////////////////////////////////////////////////////////////////
void
VerifyAndWriteFinalOutput(Real_t elapsed_time, Domain& locDom, Int_t nx, Int_t numRanks)
{
// GrindTime1 only takes a single domain into account, and is thus a good way to
// measure processor speed indepdendent of MPI parallelism. GrindTime2 takes into
// account speedups from MPI parallelism
Real_t grindTime1 = ((elapsed_time * 1e6) / locDom.cycle()) / (nx * nx * nx);
Real_t grindTime2 =
((elapsed_time * 1e6) / locDom.cycle()) / (nx * nx * nx * numRanks);
Index_t ElemId = 0;
printf("Run completed: \n");
printf(" Problem size = %i \n", nx);
printf(" MPI tasks = %i \n", numRanks);
printf(" Iteration count = %i \n", locDom.cycle());
printf(" Final Origin Energy = %12.6e \n", locDom.e(ElemId));
Real_t MaxAbsDiff = Real_t(0.0);
Real_t TotalAbsDiff = Real_t(0.0);
Real_t MaxRelDiff = Real_t(0.0);
for(Index_t j = 0; j < nx; ++j)
{
for(Index_t k = j + 1; k < nx; ++k)
{
Real_t AbsDiff = FABS(locDom.e(j * nx + k) - locDom.e(k * nx + j));
TotalAbsDiff += AbsDiff;
if(MaxAbsDiff < AbsDiff)
MaxAbsDiff = AbsDiff;
Real_t RelDiff = AbsDiff / locDom.e(k * nx + j);
if(MaxRelDiff < RelDiff)
MaxRelDiff = RelDiff;
}
}
// Quick symmetry check
printf(" Testing Plane 0 of Energy Array on rank 0:\n");
printf(" MaxAbsDiff = %12.6e\n", MaxAbsDiff);
printf(" TotalAbsDiff = %12.6e\n", TotalAbsDiff);
printf(" MaxRelDiff = %12.6e\n\n", MaxRelDiff);
// Timing information
printf("\nElapsed time = %10.2f (s)\n", elapsed_time);
printf("Grind time (us/z/c) = %10.8g (per dom) (%10.8g overall)\n", grindTime1,
grindTime2);
printf("FOM = %10.8g (z/s)\n\n",
1000.0 / grindTime2); // zones per second
return;
}
+422
Dosyayı Görüntüle
@@ -0,0 +1,422 @@
#include "lulesh.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#ifdef VIZ_MESH
# ifdef __cplusplus
extern "C"
{
# endif
# include "silo.h"
# if USE_MPI
# include "pmpio.h"
# endif
# ifdef __cplusplus
}
# endif
// Function prototypes
static void
DumpDomainToVisit(DBfile* db, Domain& domain, int myRank);
static
# if USE_MPI
// For some reason, earlier versions of g++ (e.g. 4.2) won't let me
// put the 'static' qualifier on this prototype, even if it's done
// consistently in the prototype and definition
void
DumpMultiblockObjects(DBfile* db, PMPIO_baton_t* bat, char basename[], int numRanks);
// Callback prototypes for PMPIO interface (only useful if we're
// running parallel)
static void*
LULESH_PMPIO_Create(const char* fname, const char* dname, void* udata);
static void*
LULESH_PMPIO_Open(const char* fname, const char* dname, PMPIO_iomode_t ioMode,
void* udata);
static void
LULESH_PMPIO_Close(void* file, void* udata);
# else
void
DumpMultiblockObjects(DBfile* db, char basename[], int numRanks);
# endif
/**********************************************************************/
void
DumpToVisit(Domain& domain, int numFiles, int myRank, int numRanks)
{
char subdirName[32];
char basename[32];
DBfile* db;
sprintf(basename, "lulesh_plot_c%d", domain.cycle());
sprintf(subdirName, "data_%d", myRank);
# if USE_MPI
PMPIO_baton_t* bat =
PMPIO_Init(numFiles, PMPIO_WRITE, MPI_COMM_WORLD, 10101, LULESH_PMPIO_Create,
LULESH_PMPIO_Open, LULESH_PMPIO_Close, NULL);
int myiorank = PMPIO_GroupRank(bat, myRank);
char fileName[64];
if(myiorank == 0)
strcpy(fileName, basename);
else
sprintf(fileName, "%s.%03d", basename, myiorank);
db = (DBfile*) PMPIO_WaitForBaton(bat, fileName, subdirName);
DumpDomainToVisit(db, domain, myRank);
// Processor 0 writes out bit of extra data to its file that
// describes how to stitch all the pieces together
if(myRank == 0)
{
DumpMultiblockObjects(db, bat, basename, numRanks);
}
PMPIO_HandOffBaton(bat, db);
PMPIO_Finish(bat);
# else
db = (DBfile*) DBCreate(basename, DB_CLOBBER, DB_LOCAL, NULL, DB_HDF5X);
if(db)
{
DBMkDir(db, subdirName);
DBSetDir(db, subdirName);
DumpDomainToVisit(db, domain, myRank);
DumpMultiblockObjects(db, basename, numRanks);
}
else
{
printf("Error writing out viz file - rank %d\n", myRank);
}
# endif
}
/**********************************************************************/
static void
DumpDomainToVisit(DBfile* db, Domain& domain, int myRank)
{
int ok = 0;
/* Create an option list that will give some hints to VisIt for
* printing out the cycle and time in the annotations */
DBoptlist* optlist;
/* Write out the mesh connectivity in fully unstructured format */
int shapetype[1] = { DB_ZONETYPE_HEX };
int shapesize[1] = { 8 };
int shapecnt[1] = { domain.numElem() };
int* conn = Allocate<int>(domain.numElem() * 8);
int ci = 0;
for(int ei = 0; ei < domain.numElem(); ++ei)
{
Index_t* elemToNode = domain.nodelist(ei);
for(int ni = 0; ni < 8; ++ni)
{
conn[ci++] = elemToNode[ni];
}
}
ok += DBPutZonelist2(db, "connectivity", domain.numElem(), 3, conn,
domain.numElem() * 8, 0, 0, 0, /* Not carrying ghost zones */
shapetype, shapesize, shapecnt, 1, NULL);
Release<int>(&conn);
/* Write out the mesh coordinates associated with the mesh */
const char* coordnames[3] = { "X", "Y", "Z" };
float* coords[3];
coords[0] = Allocate<float>(domain.numNode());
coords[1] = Allocate<float>(domain.numNode());
coords[2] = Allocate<float>(domain.numNode());
for(int ni = 0; ni < domain.numNode(); ++ni)
{
coords[0][ni] = float(domain.x(ni));
coords[1][ni] = float(domain.y(ni));
coords[2][ni] = float(domain.z(ni));
}
optlist = DBMakeOptlist(2);
ok += DBAddOption(optlist, DBOPT_DTIME, &domain.time());
ok += DBAddOption(optlist, DBOPT_CYCLE, &domain.cycle());
ok += DBPutUcdmesh(db, "mesh", 3, (char**) &coordnames[0], (float**) coords,
domain.numNode(), domain.numElem(), "connectivity", 0, DB_FLOAT,
optlist);
ok += DBFreeOptlist(optlist);
Release<float>(&coords[2]);
Release<float>(&coords[1]);
Release<float>(&coords[0]);
/* Write out the materials */
int* matnums = Allocate<int>(domain.numReg());
int dims[1] = { domain.numElem() }; // No mixed elements
for(int i = 0; i < domain.numReg(); ++i)
matnums[i] = i + 1;
ok += DBPutMaterial(db, "regions", "mesh", domain.numReg(), matnums,
domain.regNumList(), dims, 1, NULL, NULL, NULL, NULL, 0, DB_FLOAT,
NULL);
Release<int>(&matnums);
/* Write out pressure, energy, relvol, q */
float* e = Allocate<float>(domain.numElem());
for(int ei = 0; ei < domain.numElem(); ++ei)
{
e[ei] = float(domain.e(ei));
}
ok += DBPutUcdvar1(db, "e", "mesh", e, domain.numElem(), NULL, 0, DB_FLOAT,
DB_ZONECENT, NULL);
Release<float>(&e);
float* p = Allocate<float>(domain.numElem());
for(int ei = 0; ei < domain.numElem(); ++ei)
{
p[ei] = float(domain.p(ei));
}
ok += DBPutUcdvar1(db, "p", "mesh", p, domain.numElem(), NULL, 0, DB_FLOAT,
DB_ZONECENT, NULL);
Release<float>(&p);
float* v = Allocate<float>(domain.numElem());
for(int ei = 0; ei < domain.numElem(); ++ei)
{
v[ei] = float(domain.v(ei));
}
ok += DBPutUcdvar1(db, "v", "mesh", v, domain.numElem(), NULL, 0, DB_FLOAT,
DB_ZONECENT, NULL);
Release<float>(&v);
float* q = Allocate<float>(domain.numElem());
for(int ei = 0; ei < domain.numElem(); ++ei)
{
q[ei] = float(domain.q(ei));
}
ok += DBPutUcdvar1(db, "q", "mesh", q, domain.numElem(), NULL, 0, DB_FLOAT,
DB_ZONECENT, NULL);
Release<float>(&q);
/* Write out nodal speed, velocities */
float* zd = Allocate<float>(domain.numNode());
float* yd = Allocate<float>(domain.numNode());
float* xd = Allocate<float>(domain.numNode());
float* speed = Allocate<float>(domain.numNode());
for(int ni = 0; ni < domain.numNode(); ++ni)
{
xd[ni] = float(domain.xd(ni));
yd[ni] = float(domain.yd(ni));
zd[ni] = float(domain.zd(ni));
speed[ni] =
float(sqrt((xd[ni] * xd[ni]) + (yd[ni] * yd[ni]) + (zd[ni] * zd[ni])));
}
ok += DBPutUcdvar1(db, "speed", "mesh", speed, domain.numNode(), NULL, 0, DB_FLOAT,
DB_NODECENT, NULL);
Release<float>(&speed);
ok += DBPutUcdvar1(db, "xd", "mesh", xd, domain.numNode(), NULL, 0, DB_FLOAT,
DB_NODECENT, NULL);
Release<float>(&xd);
ok += DBPutUcdvar1(db, "yd", "mesh", yd, domain.numNode(), NULL, 0, DB_FLOAT,
DB_NODECENT, NULL);
Release<float>(&yd);
ok += DBPutUcdvar1(db, "zd", "mesh", zd, domain.numNode(), NULL, 0, DB_FLOAT,
DB_NODECENT, NULL);
Release<float>(&zd);
if(ok != 0)
{
printf("Error writing out viz file - rank %d\n", myRank);
}
}
/**********************************************************************/
# if USE_MPI
void
DumpMultiblockObjects(DBfile* db, PMPIO_baton_t* bat, char basename[], int numRanks)
# else
void
DumpMultiblockObjects(DBfile* db, char basename[], int numRanks)
# endif
{
/* MULTIBLOCK objects to tie together multiple files */
char** multimeshObjs;
char** multimatObjs;
char*** multivarObjs;
int* blockTypes;
int* varTypes;
int ok = 0;
// Make sure this list matches what's written out above
char vars[][10] = { "p", "e", "v", "q", "speed", "xd", "yd", "zd" };
int numvars = sizeof(vars) / sizeof(vars[0]);
// Reset to the root directory of the silo file
DBSetDir(db, "/");
// Allocate a bunch of space for building up the string names
multimeshObjs = Allocate<char*>(numRanks);
multimatObjs = Allocate<char*>(numRanks);
multivarObjs = Allocate<char**>(numvars);
blockTypes = Allocate<int>(numRanks);
varTypes = Allocate<int>(numRanks);
for(int v = 0; v < numvars; ++v)
{
multivarObjs[v] = Allocate<char*>(numRanks);
}
for(int i = 0; i < numRanks; ++i)
{
multimeshObjs[i] = Allocate<char>(64);
multimatObjs[i] = Allocate<char>(64);
for(int v = 0; v < numvars; ++v)
{
multivarObjs[v][i] = Allocate<char>(64);
}
blockTypes[i] = DB_UCDMESH;
varTypes[i] = DB_UCDVAR;
}
// Build up the multiobject names
for(int i = 0; i < numRanks; ++i)
{
# if USE_MPI
int iorank = PMPIO_GroupRank(bat, i);
# else
int iorank = 0;
# endif
// delete multivarObjs[i];
if(iorank == 0)
{
snprintf(multimeshObjs[i], 64, "/data_%d/mesh", i);
snprintf(multimatObjs[i], 64, "/data_%d/regions", i);
for(int v = 0; v < numvars; ++v)
{
snprintf(multivarObjs[v][i], 64, "/data_%d/%s", i, vars[v]);
}
}
else
{
snprintf(multimeshObjs[i], 64, "%s.%03d:/data_%d/mesh", basename, iorank, i);
snprintf(multimatObjs[i], 64, "%s.%03d:/data_%d/regions", basename, iorank,
i);
for(int v = 0; v < numvars; ++v)
{
snprintf(multivarObjs[v][i], 64, "%s.%03d:/data_%d/%s", basename, iorank,
i, vars[v]);
}
}
}
// Now write out the objects
ok += DBPutMultimesh(db, "mesh", numRanks, (char**) multimeshObjs, blockTypes, NULL);
ok += DBPutMultimat(db, "regions", numRanks, (char**) multimatObjs, NULL);
for(int v = 0; v < numvars; ++v)
{
ok += DBPutMultivar(db, vars[v], numRanks, (char**) multivarObjs[v], varTypes,
NULL);
}
for(int v = 0; v < numvars; ++v)
{
for(int i = 0; i < numRanks; i++)
{
Release<char>(&multivarObjs[v][i]);
}
Release<char*>(&multivarObjs[v]);
}
// Clean up
for(int i = 0; i < numRanks; i++)
{
Release<char>(&multimeshObjs[i]);
Release<char>(&multimatObjs[i]);
}
Release<char*>(&multimeshObjs);
Release<char*>(&multimatObjs);
Release<char**>(&multivarObjs);
Release<int>(&blockTypes);
Release<int>(&varTypes);
if(ok != 0)
{
printf("Error writing out multiXXX objs to viz file - rank 0\n");
}
}
# if USE_MPI
/**********************************************************************/
static void*
LULESH_PMPIO_Create(const char* fname, const char* dname, void* udata)
{
/* Create the file */
DBfile* db = DBCreate(fname, DB_CLOBBER, DB_LOCAL, NULL, DB_HDF5X);
/* Put the data in a subdirectory, so VisIt only sees the multimesh
* objects we write out in the base file */
if(db)
{
DBMkDir(db, dname);
DBSetDir(db, dname);
}
return (void*) db;
}
/**********************************************************************/
static void*
LULESH_PMPIO_Open(const char* fname, const char* dname, PMPIO_iomode_t ioMode,
void* udata)
{
/* Open the file */
DBfile* db = DBOpen(fname, DB_UNKNOWN, DB_APPEND);
/* Put the data in a subdirectory, so VisIt only sees the multimesh
* objects we write out in the base file */
if(db)
{
DBMkDir(db, dname);
DBSetDir(db, dname);
}
return (void*) db;
}
/**********************************************************************/
static void
LULESH_PMPIO_Close(void* file, void* udata)
{
DBfile* db = (DBfile*) file;
if(db)
DBClose(db);
}
# endif
#else
void
DumpToVisit(Domain& domain, int numFiles, int myRank, int numRanks)
{
if(myRank == 0)
{
printf("Must enable -DVIZ_MESH at compile time to call DumpDomain\n");
}
}
#endif
Dosya farkı çok büyük olduğundan ihmal edildi Fark Yükle
+836
Dosyayı Görüntüle
@@ -0,0 +1,836 @@
#if !defined(USE_MPI)
# error "You should specify USE_MPI=0 or USE_MPI=1 on the compile line"
#endif
// OpenMP will be compiled in if this flag is set to 1 AND the compiler beging
// used supports it (i.e. the _OPENMP symbol is defined)
#define USE_OMP 1
#if USE_MPI
# include <mpi.h>
/*
define one of these three symbols:
SEDOV_SYNC_POS_VEL_NONE
SEDOV_SYNC_POS_VEL_EARLY
SEDOV_SYNC_POS_VEL_LATE
*/
# define SEDOV_SYNC_POS_VEL_EARLY 1
#endif
#include <Kokkos_Core.hpp>
#include <Kokkos_Vector.hpp>
#include <math.h>
#include <vector>
//**************************************************
// Allow flexibility for arithmetic representations
//**************************************************
#define MAX(a, b) (((a) > (b)) ? (a) : (b))
// Precision specification
typedef float real4;
typedef double real8;
typedef long double real10; // 10 bytes on x86
typedef int Index_t; // array subscript and loop index
typedef real8 Real_t; // floating point representation
typedef int Int_t; // integer representation
enum
{
VolumeError = -1,
QStopError = -2
};
inline real4
SQRT(real4 arg)
{
return sqrtf(arg);
}
inline real8
SQRT(real8 arg)
{
return sqrt(arg);
}
inline real10
SQRT(real10 arg)
{
return sqrtl(arg);
}
inline real4
CBRT(real4 arg)
{
return cbrtf(arg);
}
inline real8
CBRT(real8 arg)
{
return cbrt(arg);
}
inline real10
CBRT(real10 arg)
{
return cbrtl(arg);
}
inline real4
FABS(real4 arg)
{
return fabsf(arg);
}
inline real8
FABS(real8 arg)
{
return fabs(arg);
}
inline real10
FABS(real10 arg)
{
return fabsl(arg);
}
// Stuff needed for boundary conditions
// 2 BCs on each of 6 hexahedral faces (12 bits)
#define XI_M 0x00007
#define XI_M_SYMM 0x00001
#define XI_M_FREE 0x00002
#define XI_M_COMM 0x00004
#define XI_P 0x00038
#define XI_P_SYMM 0x00008
#define XI_P_FREE 0x00010
#define XI_P_COMM 0x00020
#define ETA_M 0x001c0
#define ETA_M_SYMM 0x00040
#define ETA_M_FREE 0x00080
#define ETA_M_COMM 0x00100
#define ETA_P 0x00e00
#define ETA_P_SYMM 0x00200
#define ETA_P_FREE 0x00400
#define ETA_P_COMM 0x00800
#define ZETA_M 0x07000
#define ZETA_M_SYMM 0x01000
#define ZETA_M_FREE 0x02000
#define ZETA_M_COMM 0x04000
#define ZETA_P 0x38000
#define ZETA_P_SYMM 0x08000
#define ZETA_P_FREE 0x10000
#define ZETA_P_COMM 0x20000
// MPI Message Tags
#define MSG_COMM_SBN 1024
#define MSG_SYNC_POS_VEL 2048
#define MSG_MONOQ 3072
#define MAX_FIELDS_PER_MPI_COMM 6
// Assume 128 byte coherence
// Assume Real_t is an "integral power of 2" bytes wide
#define CACHE_COHERENCE_PAD_REAL (128 / sizeof(Real_t))
#define CACHE_ALIGN_REAL(n) \
(((n) + (CACHE_COHERENCE_PAD_REAL - 1)) & ~(CACHE_COHERENCE_PAD_REAL - 1))
//////////////////////////////////////////////////////
// Primary data structure
//////////////////////////////////////////////////////
/*
* The implementation of the data abstraction used for lulesh
* resides entirely in the Domain class below. You can change
* grouping and interleaving of fields here to maximize data layout
* efficiency for your underlying architecture or compiler.
*
* For example, fields can be implemented as STL objects or
* raw array pointers. As another example, individual fields
* m_x, m_y, m_z could be budled into
*
* struct { Real_t x, y, z ; } *m_coord ;
*
* allowing accessor functions such as
*
* "Real_t &x(Index_t idx) { return m_coord[idx].x ; }"
* "Real_t &y(Index_t idx) { return m_coord[idx].y ; }"
* "Real_t &z(Index_t idx) { return m_coord[idx].z ; }"
*/
class Domain
{
public:
// Constructor
Domain(Int_t numRanks, Index_t colLoc, Index_t rowLoc, Index_t planeLoc, Index_t nx,
Int_t tp, Int_t nr, Int_t balance, Int_t cost);
// Destructor
~Domain();
//
// ALLOCATION
//
void AllocateNodePersistent(Int_t numNode) // Node-centered
{
m_x.resize(numNode); // coordinates
m_y.resize(numNode);
m_z.resize(numNode);
m_xd.resize(numNode); // velocities
m_yd.resize(numNode);
m_zd.resize(numNode);
m_xdd.resize(numNode); // accelerations
m_ydd.resize(numNode);
m_zdd.resize(numNode);
m_fx.resize(numNode); // forces
m_fy.resize(numNode);
m_fz.resize(numNode);
m_nodalMass.resize(numNode); // mass
m_c_x = m_x.d_view;
m_c_y = m_y.d_view;
m_c_z = m_z.d_view;
m_c_xd = m_xd.d_view;
m_c_yd = m_yd.d_view;
m_c_zd = m_zd.d_view;
}
void AllocateElemPersistent(Int_t numElem) // Elem-centered
{
m_nodelist.resize(8 * numElem);
// elem connectivities through face
m_lxim.resize(numElem);
m_lxip.resize(numElem);
m_letam.resize(numElem);
m_letap.resize(numElem);
m_lzetam.resize(numElem);
m_lzetap.resize(numElem);
m_elemBC.resize(numElem);
m_e.resize(numElem);
m_p.resize(numElem);
m_q.resize(numElem);
m_ql.resize(numElem);
m_qq.resize(numElem);
m_v.resize(numElem);
m_volo.resize(numElem);
m_delv.resize(numElem);
m_vdov.resize(numElem);
m_arealg.resize(numElem);
m_ss.resize(numElem);
m_elemMass.resize(numElem);
m_vnew.resize(numElem);
m_c_e = m_e.d_view;
m_c_p = m_p.d_view;
m_c_q = m_q.d_view;
m_c_ql = m_ql.d_view;
m_c_qq = m_qq.d_view;
m_c_delv = m_delv.d_view;
}
void AllocateGradients(Int_t numElem, Int_t allElem)
{
// Position gradients
m_delx_xi.resize(numElem);
m_delx_eta.resize(numElem);
m_delx_zeta.resize(numElem);
// Velocity gradients
m_delv_xi.resize(allElem);
m_delv_eta.resize(allElem);
m_delv_zeta.resize(allElem);
}
void DeallocateGradients()
{
m_delx_zeta.clear();
m_delx_eta.clear();
m_delx_xi.clear();
m_delv_zeta.clear();
m_delv_eta.clear();
m_delv_xi.clear();
}
void AllocateStrains(Int_t numElem)
{
m_dxx.resize(numElem);
m_dyy.resize(numElem);
m_dzz.resize(numElem);
}
void DeallocateStrains()
{
m_dzz.clear();
m_dyy.clear();
m_dxx.clear();
}
//
// ACCESSORS
//
// Node-centered
// Nodal coordinates
KOKKOS_INLINE_FUNCTION Real_t& x(const Index_t idx) const { return m_x[idx]; }
KOKKOS_INLINE_FUNCTION Real_t& y(const Index_t idx) const { return m_y[idx]; }
KOKKOS_INLINE_FUNCTION Real_t& z(const Index_t idx) const { return m_z[idx]; }
KOKKOS_INLINE_FUNCTION Real_t c_x(const Index_t idx) const { return m_c_x[idx]; }
KOKKOS_INLINE_FUNCTION Real_t c_y(const Index_t idx) const { return m_c_y[idx]; }
KOKKOS_INLINE_FUNCTION Real_t c_z(const Index_t idx) const { return m_c_z[idx]; }
// Nodal velocities
KOKKOS_INLINE_FUNCTION Real_t& xd(const Index_t idx) const { return m_xd[idx]; }
KOKKOS_INLINE_FUNCTION Real_t& yd(const Index_t idx) const { return m_yd[idx]; }
KOKKOS_INLINE_FUNCTION Real_t& zd(const Index_t idx) const { return m_zd[idx]; }
KOKKOS_INLINE_FUNCTION Real_t c_xd(const Index_t idx) const { return m_c_xd[idx]; }
KOKKOS_INLINE_FUNCTION Real_t c_yd(const Index_t idx) const { return m_c_yd[idx]; }
KOKKOS_INLINE_FUNCTION Real_t c_zd(const Index_t idx) const { return m_c_zd[idx]; }
// Nodal accelerations
KOKKOS_INLINE_FUNCTION Real_t& xdd(const Index_t idx) const { return m_xdd[idx]; }
KOKKOS_INLINE_FUNCTION Real_t& ydd(const Index_t idx) const { return m_ydd[idx]; }
KOKKOS_INLINE_FUNCTION Real_t& zdd(const Index_t idx) const { return m_zdd[idx]; }
// Nodal forces
KOKKOS_INLINE_FUNCTION Real_t& fx(const Index_t idx) const { return m_fx[idx]; }
KOKKOS_INLINE_FUNCTION Real_t& fy(const Index_t idx) const { return m_fy[idx]; }
KOKKOS_INLINE_FUNCTION Real_t& fz(const Index_t idx) const { return m_fz[idx]; }
// Nodal mass
KOKKOS_INLINE_FUNCTION Real_t& nodalMass(const Index_t idx) const
{
return m_nodalMass[idx];
}
// Nodes on symmertry planes
Index_t symmX(const Index_t idx) const { return m_symmX[idx]; }
Index_t symmY(const Index_t idx) const { return m_symmY[idx]; }
Index_t symmZ(const Index_t idx) const { return m_symmZ[idx]; }
bool symmXempty() { return m_symmX.empty(); }
bool symmYempty() { return m_symmY.empty(); }
bool symmZempty() { return m_symmZ.empty(); }
//
// Element-centered
//
Index_t& regElemSize(Index_t idx) { return m_regElemSize[idx]; }
Index_t& regNumList(Index_t idx) { return m_regNumList[idx]; }
Index_t* regNumList() { return &m_regNumList[0]; }
Index_t* regElemlist(Int_t r) { return m_regElemlist[r]; }
Index_t& regElemlist(const Int_t r, Index_t idx) const
{
return m_regElemlist[r][idx];
}
Index_t* nodelist(Index_t idx) const { return &m_nodelist[Index_t(8) * idx]; }
// elem connectivities through face
Index_t& lxim(const Index_t idx) const { return m_lxim[idx]; }
Index_t& lxip(const Index_t idx) const { return m_lxip[idx]; }
Index_t& letam(const Index_t idx) const { return m_letam[idx]; }
Index_t& letap(const Index_t idx) const { return m_letap[idx]; }
Index_t& lzetam(const Index_t idx) const { return m_lzetam[idx]; }
Index_t& lzetap(const Index_t idx) const { return m_lzetap[idx]; }
// elem face symm/free-surface flag
Int_t& elemBC(const Index_t idx) const { return m_elemBC[idx]; }
// Principal strains - temporary
KOKKOS_INLINE_FUNCTION Real_t& dxx(const Index_t idx) const { return m_dxx[idx]; }
KOKKOS_INLINE_FUNCTION Real_t& dyy(const Index_t idx) const { return m_dyy[idx]; }
KOKKOS_INLINE_FUNCTION Real_t& dzz(const Index_t idx) const { return m_dzz[idx]; }
// New relative volume - temporary
KOKKOS_INLINE_FUNCTION Real_t& vnew(const Index_t idx) const { return m_vnew[idx]; }
// Velocity gradient - temporary
KOKKOS_INLINE_FUNCTION Real_t& delv_xi(const Index_t idx) const
{
return m_delv_xi[idx];
}
KOKKOS_INLINE_FUNCTION Real_t& delv_eta(const Index_t idx) const
{
return m_delv_eta[idx];
}
KOKKOS_INLINE_FUNCTION Real_t& delv_zeta(const Index_t idx) const
{
return m_delv_zeta[idx];
}
// Position gradient - temporary
KOKKOS_INLINE_FUNCTION Real_t& delx_xi(const Index_t idx) const
{
return m_delx_xi[idx];
}
KOKKOS_INLINE_FUNCTION Real_t& delx_eta(const Index_t idx) const
{
return m_delx_eta[idx];
}
KOKKOS_INLINE_FUNCTION Real_t& delx_zeta(const Index_t idx) const
{
return m_delx_zeta[idx];
}
// Energy
KOKKOS_INLINE_FUNCTION Real_t& e(const Index_t idx) const { return m_e[idx]; }
KOKKOS_INLINE_FUNCTION Real_t c_e(const Index_t idx) const { return m_c_e[idx]; }
// Pressure
KOKKOS_INLINE_FUNCTION Real_t& p(const Index_t idx) const { return m_p[idx]; }
KOKKOS_INLINE_FUNCTION Real_t c_p(const Index_t idx) const { return m_c_p[idx]; }
// Artificial viscosity
KOKKOS_INLINE_FUNCTION Real_t& q(const Index_t idx) const { return m_q[idx]; }
KOKKOS_INLINE_FUNCTION Real_t c_q(const Index_t idx) const { return m_c_q[idx]; }
// Linear term for q
KOKKOS_INLINE_FUNCTION Real_t& ql(const Index_t idx) const { return m_ql[idx]; }
KOKKOS_INLINE_FUNCTION Real_t c_ql(const Index_t idx) const { return m_c_ql[idx]; }
// Quadratic term for q
KOKKOS_INLINE_FUNCTION Real_t& qq(const Index_t idx) const { return m_qq[idx]; }
KOKKOS_INLINE_FUNCTION Real_t c_qq(const Index_t idx) const { return m_c_qq[idx]; }
// Relative volume
KOKKOS_INLINE_FUNCTION Real_t& v(const Index_t idx) const { return m_v[idx]; }
KOKKOS_INLINE_FUNCTION Real_t& delv(const Index_t idx) const { return m_delv[idx]; }
KOKKOS_INLINE_FUNCTION Real_t c_delv(const Index_t idx) const
{
return m_c_delv[idx];
}
// Reference volume
KOKKOS_INLINE_FUNCTION Real_t& volo(Index_t idx) const { return m_volo[idx]; }
// volume derivative over volume
KOKKOS_INLINE_FUNCTION Real_t& vdov(Index_t idx) const { return m_vdov[idx]; }
// Element characteristic length
KOKKOS_INLINE_FUNCTION Real_t& arealg(Index_t idx) const { return m_arealg[idx]; }
// Sound speed
KOKKOS_INLINE_FUNCTION Real_t& ss(const Index_t idx) const { return m_ss[idx]; }
// Element mass
KOKKOS_INLINE_FUNCTION Real_t& elemMass(const Index_t idx) const
{
return m_elemMass[idx];
}
KOKKOS_INLINE_FUNCTION Index_t nodeElemCount(Index_t idx) const
{
return m_nodeElemStart[idx + 1] - m_nodeElemStart[idx];
}
KOKKOS_INLINE_FUNCTION Index_t* nodeElemCornerList(Index_t idx) const
{
return &m_nodeElemCornerList[m_nodeElemStart[idx]];
}
// Parameters
// Cutoffs
KOKKOS_INLINE_FUNCTION Real_t u_cut() const { return m_u_cut; }
KOKKOS_INLINE_FUNCTION Real_t e_cut() const { return m_e_cut; }
KOKKOS_INLINE_FUNCTION Real_t p_cut() const { return m_p_cut; }
KOKKOS_INLINE_FUNCTION Real_t q_cut() const { return m_q_cut; }
KOKKOS_INLINE_FUNCTION Real_t v_cut() const { return m_v_cut; }
// Other constants (usually are settable via input file in real codes)
KOKKOS_INLINE_FUNCTION Real_t hgcoef() const { return m_hgcoef; }
KOKKOS_INLINE_FUNCTION Real_t qstop() const { return m_qstop; }
KOKKOS_INLINE_FUNCTION Real_t monoq_max_slope() const { return m_monoq_max_slope; }
KOKKOS_INLINE_FUNCTION Real_t monoq_limiter_mult() const
{
return m_monoq_limiter_mult;
}
KOKKOS_INLINE_FUNCTION Real_t ss4o3() const { return m_ss4o3; }
KOKKOS_INLINE_FUNCTION Real_t qlc_monoq() const { return m_qlc_monoq; }
KOKKOS_INLINE_FUNCTION Real_t qqc_monoq() const { return m_qqc_monoq; }
KOKKOS_INLINE_FUNCTION Real_t qqc() const { return m_qqc; }
KOKKOS_INLINE_FUNCTION Real_t eosvmax() const { return m_eosvmax; }
KOKKOS_INLINE_FUNCTION Real_t eosvmin() const { return m_eosvmin; }
KOKKOS_INLINE_FUNCTION Real_t pmin() const { return m_pmin; }
KOKKOS_INLINE_FUNCTION Real_t emin() const { return m_emin; }
KOKKOS_INLINE_FUNCTION Real_t dvovmax() const { return m_dvovmax; }
KOKKOS_INLINE_FUNCTION Real_t refdens() const { return m_refdens; }
// Timestep controls, etc...
Real_t& time() { return m_time; }
Real_t& deltatime() { return m_deltatime; }
Real_t& deltatimemultlb() { return m_deltatimemultlb; }
Real_t& deltatimemultub() { return m_deltatimemultub; }
Real_t& stoptime() { return m_stoptime; }
Real_t& dtcourant() { return m_dtcourant; }
Real_t& dthydro() { return m_dthydro; }
Real_t& dtmax() { return m_dtmax; }
Real_t& dtfixed() { return m_dtfixed; }
Int_t& cycle() { return m_cycle; }
Index_t& numRanks() { return m_numRanks; }
Index_t& colLoc() { return m_colLoc; }
Index_t& rowLoc() { return m_rowLoc; }
Index_t& planeLoc() { return m_planeLoc; }
Index_t& tp() { return m_tp; }
Index_t& sizeX() { return m_sizeX; }
Index_t& sizeY() { return m_sizeY; }
Index_t& sizeZ() { return m_sizeZ; }
Index_t& numReg() { return m_numReg; }
Int_t& cost() { return m_cost; }
Index_t& numElem() { return m_numElem; }
Index_t& numNode() { return m_numNode; }
Index_t& maxPlaneSize() { return m_maxPlaneSize; }
Index_t& maxEdgeSize() { return m_maxEdgeSize; }
//
// MPI-Related additional data
//
#if USE_MPI
// Communication Work space
Real_t* commDataSend;
Real_t* commDataRecv;
// Maximum number of block neighbors
MPI_Request recvRequest[26]; // 6 faces + 12 edges + 8 corners
MPI_Request sendRequest[26]; // 6 faces + 12 edges + 8 corners
#endif
private:
void BuildMesh(Int_t nx, Int_t edgeNodes, Int_t edgeElems);
void SetupThreadSupportStructures();
void CreateRegionIndexSets(Int_t nreg, Int_t balance);
void SetupCommBuffers(Int_t edgeNodes);
void SetupSymmetryPlanes(Int_t edgeNodes);
void SetupElementConnectivities(Int_t edgeElems);
void SetupBoundaryConditions(Int_t edgeElems);
//
// IMPLEMENTATION
//
/* Node-centered */
Kokkos::vector<Real_t> m_x; /* coordinates */
Kokkos::vector<Real_t> m_y;
Kokkos::vector<Real_t> m_z;
Kokkos::View<const Real_t*, Kokkos::MemoryTraits<Kokkos::RandomAccess>>
m_c_x; /* coordinates */
Kokkos::View<const Real_t*, Kokkos::MemoryTraits<Kokkos::RandomAccess>>
m_c_y; /* coordinates */
Kokkos::View<const Real_t*, Kokkos::MemoryTraits<Kokkos::RandomAccess>>
m_c_z; /* coordinates */
Kokkos::vector<Real_t> m_xd; /* velocities */
Kokkos::vector<Real_t> m_yd;
Kokkos::vector<Real_t> m_zd;
Kokkos::View<const Real_t*, Kokkos::MemoryTraits<Kokkos::RandomAccess>>
m_c_xd; /* coordinates */
Kokkos::View<const Real_t*, Kokkos::MemoryTraits<Kokkos::RandomAccess>>
m_c_yd; /* coordinates */
Kokkos::View<const Real_t*, Kokkos::MemoryTraits<Kokkos::RandomAccess>>
m_c_zd; /* coordinates */
Kokkos::vector<Real_t> m_xdd; /* accelerations */
Kokkos::vector<Real_t> m_ydd;
Kokkos::vector<Real_t> m_zdd;
Kokkos::vector<Real_t> m_fx; /* forces */
Kokkos::vector<Real_t> m_fy;
Kokkos::vector<Real_t> m_fz;
Kokkos::vector<Real_t> m_nodalMass; /* mass */
Kokkos::vector<Index_t> m_symmX; /* symmetry plane nodesets */
Kokkos::vector<Index_t> m_symmY;
Kokkos::vector<Index_t> m_symmZ;
// Element-centered
// Region information
Int_t m_numReg;
Int_t m_cost; // imbalance cost
Index_t* m_regElemSize; // Size of region sets
Index_t* m_regNumList; // Region number per domain element
Index_t** m_regElemlist; // region indexset
Kokkos::vector<Index_t> m_nodelist; /* elemToNode connectivity */
Kokkos::vector<Index_t> m_lxim; /* element connectivity across each face */
Kokkos::vector<Index_t> m_lxip;
Kokkos::vector<Index_t> m_letam;
Kokkos::vector<Index_t> m_letap;
Kokkos::vector<Index_t> m_lzetam;
Kokkos::vector<Index_t> m_lzetap;
Kokkos::vector<Int_t> m_elemBC; /* symmetry/free-surface flags for each elem face */
Kokkos::vector<Real_t> m_dxx; /* principal strains -- temporary */
Kokkos::vector<Real_t> m_dyy;
Kokkos::vector<Real_t> m_dzz;
Kokkos::vector<Real_t> m_delv_xi; /* velocity gradient -- temporary */
Kokkos::vector<Real_t> m_delv_eta;
Kokkos::vector<Real_t> m_delv_zeta;
Kokkos::vector<Real_t> m_delx_xi; /* coordinate gradient -- temporary */
Kokkos::vector<Real_t> m_delx_eta;
Kokkos::vector<Real_t> m_delx_zeta;
Kokkos::vector<Real_t> m_e; /* energy */
Kokkos::vector<Real_t> m_p; /* pressure */
Kokkos::vector<Real_t> m_q; /* q */
Kokkos::vector<Real_t> m_ql; /* linear term for q */
Kokkos::vector<Real_t> m_qq; /* quadratic term for q */
Kokkos::vector<Real_t> m_v; /* relative volume */
Kokkos::vector<Real_t> m_volo; /* reference volume */
Kokkos::vector<Real_t> m_vnew; /* new relative volume -- temporary */
Kokkos::vector<Real_t> m_delv; /* m_vnew - m_v */
Kokkos::vector<Real_t> m_vdov; /* volume derivative over volume */
Kokkos::View<const Real_t*, Kokkos::MemoryTraits<Kokkos::RandomAccess>>
m_c_e; /* coordinates */
Kokkos::View<const Real_t*, Kokkos::MemoryTraits<Kokkos::RandomAccess>>
m_c_p; /* coordinates */
Kokkos::View<const Real_t*, Kokkos::MemoryTraits<Kokkos::RandomAccess>>
m_c_q; /* coordinates */
Kokkos::View<const Real_t*, Kokkos::MemoryTraits<Kokkos::RandomAccess>>
m_c_ql; /* coordinates */
Kokkos::View<const Real_t*, Kokkos::MemoryTraits<Kokkos::RandomAccess>>
m_c_qq; /* coordinates */
Kokkos::View<const Real_t*, Kokkos::MemoryTraits<Kokkos::RandomAccess>>
m_c_delv; /* coordinates */
Kokkos::vector<Real_t> m_arealg; /* characteristic length of an element */
Kokkos::vector<Real_t> m_ss; /* "sound speed" */
Kokkos::vector<Real_t> m_elemMass; /* mass */
// Cutoffs (treat as constants)
const Real_t m_e_cut; // energy tolerance
const Real_t m_p_cut; // pressure tolerance
const Real_t m_q_cut; // q tolerance
const Real_t m_v_cut; // relative volume tolerance
const Real_t m_u_cut; // velocity tolerance
// Other constants (usually setable, but hardcoded in this proxy app)
const Real_t m_hgcoef; // hourglass control
const Real_t m_ss4o3;
const Real_t m_qstop; // excessive q indicator
const Real_t m_monoq_max_slope;
const Real_t m_monoq_limiter_mult;
const Real_t m_qlc_monoq; // linear term coef for q
const Real_t m_qqc_monoq; // quadratic term coef for q
const Real_t m_qqc;
const Real_t m_eosvmax;
const Real_t m_eosvmin;
const Real_t m_pmin; // pressure floor
const Real_t m_emin; // energy floor
const Real_t m_dvovmax; // maximum allowable volume change
const Real_t m_refdens; // reference density
// Variables to keep track of timestep, simulation time, and cycle
Real_t m_dtcourant; // courant constraint
Real_t m_dthydro; // volume change constraint
Int_t m_cycle; // iteration count for simulation
Real_t m_dtfixed; // fixed time increment
Real_t m_time; // current time
Real_t m_deltatime; // variable time increment
Real_t m_deltatimemultlb;
Real_t m_deltatimemultub;
Real_t m_dtmax; // maximum allowable time increment
Real_t m_stoptime; // end time for simulation
Int_t m_numRanks;
Index_t m_colLoc;
Index_t m_rowLoc;
Index_t m_planeLoc;
Index_t m_tp;
Index_t m_sizeX;
Index_t m_sizeY;
Index_t m_sizeZ;
Index_t m_numElem;
Index_t m_numNode;
Index_t m_maxPlaneSize;
Index_t m_maxEdgeSize;
// OMP hack
Index_t* m_nodeElemStart;
Index_t* m_nodeElemCornerList;
// Used in setup
Index_t m_rowMin, m_rowMax;
Index_t m_colMin, m_colMax;
Index_t m_planeMin, m_planeMax;
};
typedef Real_t& (Domain::*Domain_member)(Index_t) const;
struct cmdLineOpts
{
Int_t its; // -i
Int_t nx; // -s
Int_t numReg; // -r
Int_t numFiles; // -f
Int_t showProg; // -p
Int_t quiet; // -q
Int_t viz; // -v
Int_t cost; // -c
Int_t balance; // -b
Int_t do_atomic; // -a
};
// Function Prototypes
// lulesh-par
/*Real_t CalcElemVolume( const Real_t x[8],
const Real_t y[8],
const Real_t z[8]);*/
// lulesh-util
void
ParseCommandLineOptions(int argc, char* argv[], Int_t myRank, struct cmdLineOpts* opts);
void
VerifyAndWriteFinalOutput(Real_t elapsed_time, Domain& locDom, Int_t nx, Int_t numRanks);
// lulesh-viz
void
DumpToVisit(Domain& domain, int numFiles, int myRank, int numRanks);
// lulesh-comm
void
CommRecv(Domain& domain, Int_t msgType, Index_t xferFields, Index_t dx, Index_t dy,
Index_t dz, bool doRecv, bool planeOnly);
void
CommSend(Domain& domain, Int_t msgType, Index_t xferFields, Domain_member* fieldData,
Index_t dx, Index_t dy, Index_t dz, bool doSend, bool planeOnly);
void
CommSBN(Domain& domain, Int_t xferFields, Domain_member* fieldData);
void
CommSyncPosVel(Domain& domain);
void
CommMonoQ(Domain& domain);
// lulesh-init
void
InitMeshDecomp(Int_t numRanks, Int_t myRank, Int_t* col, Int_t* row, Int_t* plane,
Int_t* side);
/*********************************/
/* Data structure implementation */
/*********************************/
/* might want to add access methods so that memory can be */
/* better managed, as in luleshFT */
template <typename T>
T*
Allocate(size_t size)
{
return static_cast<T*>(Kokkos::kokkos_malloc<>(sizeof(T) * size));
}
template <typename T>
void
Release(T** ptr)
{
if(*ptr != NULL)
{
Kokkos::kokkos_free<>(*ptr);
*ptr = NULL;
}
}
struct MinFinder
{
Real_t val;
int i;
KOKKOS_INLINE_FUNCTION
MinFinder()
: val(100000000000000000000.0000)
, i(-1)
{}
KOKKOS_INLINE_FUNCTION
MinFinder(const double& val_, const int& i_)
: val(val_)
, i(i_)
{}
KOKKOS_INLINE_FUNCTION
MinFinder(const MinFinder& src)
: val(src.val)
, i(src.i)
{}
// overloading += operator to do the max assignment
KOKKOS_INLINE_FUNCTION
void operator+=(MinFinder& src)
{
if(src.val < val)
{
val = src.val;
i = src.i;
}
}
KOKKOS_INLINE_FUNCTION
void operator+=(const volatile MinFinder& src) volatile
{
if(src.val < val)
{
val = src.val;
i = src.i;
}
}
};
struct reduce_double3
{
double x, y, z;
KOKKOS_INLINE_FUNCTION
reduce_double3()
{
x = 0.0;
y = 0.0;
z = 0.0;
}
KOKKOS_INLINE_FUNCTION
void operator+=(const reduce_double3& src)
{
x += src.x;
y += src.y;
z += src.z;
}
};
+651
Dosyayı Görüntüle
@@ -0,0 +1,651 @@
#if !defined(USE_MPI)
# error "You should specify USE_MPI=0 or USE_MPI=1 on the compile line"
#endif
// OpenMP will be compiled in if this flag is set to 1 AND the compiler beging
// used supports it (i.e. the _OPENMP symbol is defined)
#define USE_OMP 1
#if USE_MPI
# include <mpi.h>
#endif
#include <mpi.h>
/*
define one of these three symbols:
SEDOV_SYNC_POS_VEL_NONE
SEDOV_SYNC_POS_VEL_EARLY
SEDOV_SYNC_POS_VEL_LATE
*/
#define SEDOV_SYNC_POS_VEL_EARLY 1
#include <math.h>
#include <vector>
//**************************************************
// Allow flexibility for arithmetic representations
//**************************************************
#define MAX(a, b) (((a) > (b)) ? (a) : (b))
// Precision specification
typedef float real4;
typedef double real8;
typedef long double real10; // 10 bytes on x86
typedef int Index_t; // array subscript and loop index
typedef real8 Real_t; // floating point representation
typedef int Int_t; // integer representation
enum
{
VolumeError = -1,
QStopError = -2
};
inline real4
SQRT(real4 arg)
{
return sqrtf(arg);
}
inline real8
SQRT(real8 arg)
{
return sqrt(arg);
}
inline real10
SQRT(real10 arg)
{
return sqrtl(arg);
}
inline real4
CBRT(real4 arg)
{
return cbrtf(arg);
}
inline real8
CBRT(real8 arg)
{
return cbrt(arg);
}
inline real10
CBRT(real10 arg)
{
return cbrtl(arg);
}
inline real4
FABS(real4 arg)
{
return fabsf(arg);
}
inline real8
FABS(real8 arg)
{
return fabs(arg);
}
inline real10
FABS(real10 arg)
{
return fabsl(arg);
}
// Stuff needed for boundary conditions
// 2 BCs on each of 6 hexahedral faces (12 bits)
#define XI_M 0x00007
#define XI_M_SYMM 0x00001
#define XI_M_FREE 0x00002
#define XI_M_COMM 0x00004
#define XI_P 0x00038
#define XI_P_SYMM 0x00008
#define XI_P_FREE 0x00010
#define XI_P_COMM 0x00020
#define ETA_M 0x001c0
#define ETA_M_SYMM 0x00040
#define ETA_M_FREE 0x00080
#define ETA_M_COMM 0x00100
#define ETA_P 0x00e00
#define ETA_P_SYMM 0x00200
#define ETA_P_FREE 0x00400
#define ETA_P_COMM 0x00800
#define ZETA_M 0x07000
#define ZETA_M_SYMM 0x01000
#define ZETA_M_FREE 0x02000
#define ZETA_M_COMM 0x04000
#define ZETA_P 0x38000
#define ZETA_P_SYMM 0x08000
#define ZETA_P_FREE 0x10000
#define ZETA_P_COMM 0x20000
// MPI Message Tags
#define MSG_COMM_SBN 1024
#define MSG_SYNC_POS_VEL 2048
#define MSG_MONOQ 3072
#define MAX_FIELDS_PER_MPI_COMM 6
// Assume 128 byte coherence
// Assume Real_t is an "integral power of 2" bytes wide
#define CACHE_COHERENCE_PAD_REAL (128 / sizeof(Real_t))
#define CACHE_ALIGN_REAL(n) \
(((n) + (CACHE_COHERENCE_PAD_REAL - 1)) & ~(CACHE_COHERENCE_PAD_REAL - 1))
//////////////////////////////////////////////////////
// Primary data structure
//////////////////////////////////////////////////////
/*
* The implementation of the data abstraction used for lulesh
* resides entirely in the Domain class below. You can change
* grouping and interleaving of fields here to maximize data layout
* efficiency for your underlying architecture or compiler.
*
* For example, fields can be implemented as STL objects or
* raw array pointers. As another example, individual fields
* m_x, m_y, m_z could be budled into
*
* struct { Real_t x, y, z ; } *m_coord ;
*
* allowing accessor functions such as
*
* "Real_t &x(Index_t idx) { return m_coord[idx].x ; }"
* "Real_t &y(Index_t idx) { return m_coord[idx].y ; }"
* "Real_t &z(Index_t idx) { return m_coord[idx].z ; }"
*/
class Domain
{
public:
// Constructor
Domain(Int_t numRanks, Index_t colLoc, Index_t rowLoc, Index_t planeLoc, Index_t nx,
Int_t tp, Int_t nr, Int_t balance, Int_t cost);
//
// ALLOCATION
//
void AllocateNodePersistent(Int_t numNode) // Node-centered
{
m_coord.resize(numNode); // coordinates
m_vel.resize(numNode); // velocities
m_acc.resize(numNode); // accelerations
m_force.resize(numNode); // forces
m_nodalMass.resize(numNode); // mass
}
void AllocateElemPersistent(Int_t numElem) // Elem-centered
{
m_nodelist.resize(8 * numElem);
// elem connectivities through face
m_faceToElem.resize(numElem);
m_elemBC.resize(numElem);
m_e.resize(numElem);
m_pq.resize(numElem);
m_qlqq.resize(numElem);
m_vol.resize(numElem);
m_delv.resize(numElem);
m_vdov.resize(numElem);
m_arealg.resize(numElem);
m_ss.resize(numElem);
m_elemMass.resize(numElem);
}
void AllocateGradients(Int_t numElem, Int_t allElem)
{
// Position gradients
m_delx_xi.resize(numElem);
m_delx_eta.resize(numElem);
m_delx_zeta.resize(numElem);
// Velocity gradients
m_delv_xi.resize(allElem);
m_delv_eta.resize(allElem);
m_delv_zeta.resize(allElem);
}
void DeallocateGradients()
{
m_delx_zeta.clear();
m_delx_eta.clear();
m_delx_xi.clear();
m_delv_zeta.clear();
m_delv_eta.clear();
m_delv_xi.clear();
}
void AllocateStrains(Int_t numElem)
{
m_dxx.resize(numElem);
m_dyy.resize(numElem);
m_dzz.resize(numElem);
}
void DeallocateStrains()
{
m_dzz.clear();
m_dyy.clear();
m_dxx.clear();
}
//
// ACCESSORS
//
// Node-centered
// Nodal coordinates
Real_t& x(Index_t idx) { return m_coord[idx].x; }
Real_t& y(Index_t idx) { return m_coord[idx].y; }
Real_t& z(Index_t idx) { return m_coord[idx].z; }
// Nodal velocities
Real_t& xd(Index_t idx) { return m_vel[idx].x; }
Real_t& yd(Index_t idx) { return m_vel[idx].y; }
Real_t& zd(Index_t idx) { return m_vel[idx].z; }
// Nodal accelerations
Real_t& xdd(Index_t idx) { return m_acc[idx].x; }
Real_t& ydd(Index_t idx) { return m_acc[idx].y; }
Real_t& zdd(Index_t idx) { return m_acc[idx].z; }
// Nodal forces
Real_t& fx(Index_t idx) { return m_force[idx].x; }
Real_t& fy(Index_t idx) { return m_force[idx].y; }
Real_t& fz(Index_t idx) { return m_force[idx].z; }
// Nodal mass
Real_t& nodalMass(Index_t idx) { return m_nodalMass[idx]; }
// Nodes on symmertry planes
Index_t symmX(Index_t idx) { return m_symmX[idx]; }
Index_t symmY(Index_t idx) { return m_symmY[idx]; }
Index_t symmZ(Index_t idx) { return m_symmZ[idx]; }
bool symmXempty() { return m_symmX.empty(); }
bool symmYempty() { return m_symmY.empty(); }
bool symmZempty() { return m_symmZ.empty(); }
//
// Element-centered
//
Index_t& regElemSize(Index_t idx) { return m_regElemSize[idx]; }
Index_t& regNumList(Index_t idx) { return m_regNumList[idx]; }
Index_t* regNumList() { return &m_regNumList[0]; }
Index_t* regElemlist(Int_t r) { return m_regElemlist[r]; }
Index_t& regElemlist(Int_t r, Index_t idx) { return m_regElemlist[r][idx]; }
Index_t* nodelist(Index_t idx) { return &m_nodelist[Index_t(8) * idx]; }
// elem connectivities through face
Index_t& lxim(Index_t idx) { return m_faceToElem[idx].lxim; }
Index_t& lxip(Index_t idx) { return m_faceToElem[idx].lxip; }
Index_t& letam(Index_t idx) { return m_faceToElem[idx].letam; }
Index_t& letap(Index_t idx) { return m_faceToElem[idx].letap; }
Index_t& lzetam(Index_t idx) { return m_faceToElem[idx].lzetam; }
Index_t& lzetap(Index_t idx) { return m_faceToElem[idx].lzetap; }
// elem face symm/free-surface flag
Int_t& elemBC(Index_t idx) { return m_elemBC[idx]; }
// Principal strains - temporary
Real_t& dxx(Index_t idx) { return m_dxx[idx]; }
Real_t& dyy(Index_t idx) { return m_dyy[idx]; }
Real_t& dzz(Index_t idx) { return m_dzz[idx]; }
// Velocity gradient - temporary
Real_t& delv_xi(Index_t idx) { return m_delv_xi[idx]; }
Real_t& delv_eta(Index_t idx) { return m_delv_eta[idx]; }
Real_t& delv_zeta(Index_t idx) { return m_delv_zeta[idx]; }
// Position gradient - temporary
Real_t& delx_xi(Index_t idx) { return m_delx_xi[idx]; }
Real_t& delx_eta(Index_t idx) { return m_delx_eta[idx]; }
Real_t& delx_zeta(Index_t idx) { return m_delx_zeta[idx]; }
// Energy
Real_t& e(Index_t idx) { return m_e[idx]; }
// Pressure
Real_t& p(Index_t idx) { return m_pq[idx].p; }
// Artificial viscosity
Real_t& q(Index_t idx) { return m_pq[idx].q; }
// Linear term for q
Real_t& ql(Index_t idx) { return m_qlqq[idx].ql; }
// Quadratic term for q
Real_t& qq(Index_t idx) { return m_qlqq[idx].qq; }
Real_t& delv(Index_t idx) { return m_delv[idx]; }
// Relative volume
Real_t& v(Index_t idx) { return m_vol[idx].v; }
// Reference volume
Real_t& volo(Index_t idx) { return m_vol[idx].volo; }
// volume derivative over volume
Real_t& vdov(Index_t idx) { return m_vdov[idx]; }
// Element characteristic length
Real_t& arealg(Index_t idx) { return m_arealg[idx]; }
// Sound speed
Real_t& ss(Index_t idx) { return m_ss[idx]; }
// Element mass
Real_t& elemMass(Index_t idx) { return m_elemMass[idx]; }
Index_t nodeElemCount(Index_t idx)
{
return m_nodeElemStart[idx + 1] - m_nodeElemStart[idx];
}
Index_t* nodeElemCornerList(Index_t idx)
{
return &m_nodeElemCornerList[m_nodeElemStart[idx]];
}
// Parameters
// Cutoffs
Real_t u_cut() const { return m_u_cut; }
Real_t e_cut() const { return m_e_cut; }
Real_t p_cut() const { return m_p_cut; }
Real_t q_cut() const { return m_q_cut; }
Real_t v_cut() const { return m_v_cut; }
// Other constants (usually are settable via input file in real codes)
Real_t hgcoef() const { return m_hgcoef; }
Real_t qstop() const { return m_qstop; }
Real_t monoq_max_slope() const { return m_monoq_max_slope; }
Real_t monoq_limiter_mult() const { return m_monoq_limiter_mult; }
Real_t ss4o3() const { return m_ss4o3; }
Real_t qlc_monoq() const { return m_qlc_monoq; }
Real_t qqc_monoq() const { return m_qqc_monoq; }
Real_t qqc() const { return m_qqc; }
Real_t eosvmax() const { return m_eosvmax; }
Real_t eosvmin() const { return m_eosvmin; }
Real_t pmin() const { return m_pmin; }
Real_t emin() const { return m_emin; }
Real_t dvovmax() const { return m_dvovmax; }
Real_t refdens() const { return m_refdens; }
// Timestep controls, etc...
Real_t& time() { return m_time; }
Real_t& deltatime() { return m_deltatime; }
Real_t& deltatimemultlb() { return m_deltatimemultlb; }
Real_t& deltatimemultub() { return m_deltatimemultub; }
Real_t& stoptime() { return m_stoptime; }
Real_t& dtcourant() { return m_dtcourant; }
Real_t& dthydro() { return m_dthydro; }
Real_t& dtmax() { return m_dtmax; }
Real_t& dtfixed() { return m_dtfixed; }
Int_t& cycle() { return m_cycle; }
Index_t& numRanks() { return m_numRanks; }
Index_t& colLoc() { return m_colLoc; }
Index_t& rowLoc() { return m_rowLoc; }
Index_t& planeLoc() { return m_planeLoc; }
Index_t& tp() { return m_tp; }
Index_t& sizeX() { return m_sizeX; }
Index_t& sizeY() { return m_sizeY; }
Index_t& sizeZ() { return m_sizeZ; }
Index_t& numReg() { return m_numReg; }
Int_t& cost() { return m_cost; }
Index_t& numElem() { return m_numElem; }
Index_t& numNode() { return m_numNode; }
Index_t& maxPlaneSize() { return m_maxPlaneSize; }
Index_t& maxEdgeSize() { return m_maxEdgeSize; }
//
// MPI-Related additional data
//
#if USE_MPI
// Communication Work space
Real_t* commDataSend;
Real_t* commDataRecv;
// Maximum number of block neighbors
MPI_Request recvRequest[26]; // 6 faces + 12 edges + 8 corners
MPI_Request sendRequest[26]; // 6 faces + 12 edges + 8 corners
#endif
private:
void BuildMesh(Int_t nx, Int_t edgeNodes, Int_t edgeElems);
void SetupThreadSupportStructures();
void CreateRegionIndexSets(Int_t nreg, Int_t balance);
void SetupCommBuffers(Int_t edgeNodes);
void SetupSymmetryPlanes(Int_t edgeNodes);
void SetupElementConnectivities(Int_t edgeElems);
void SetupBoundaryConditions(Int_t edgeElems);
//
// IMPLEMENTATION
//
/* Node-centered */
struct Tuple3
{
Real_t x, y, z;
};
Kokkos::vector<Tuple3> m_coord; /* coordinates */
Kokkos::vector<Tuple3> m_vel; /* velocities */
Kokkos::vector<Tuple3> m_acc; /* accelerations */
Kokkos::vector<Tuple3> m_force; /* forces */
Kokkos::vector<Real_t> m_nodalMass; /* mass */
Kokkos::vector<Index_t> m_symmX; /* symmetry plane nodesets */
Kokkos::vector<Index_t> m_symmY;
Kokkos::vector<Index_t> m_symmZ;
// Element-centered
// Region information
Int_t m_numReg;
Int_t m_cost; // imbalance cost
Index_t* m_regElemSize; // Size of region sets
Index_t* m_regNumList; // Region number per domain element
Index_t** m_regElemlist; // region indexset
Kokkos::vector<Index_t> m_nodelist; /* elemToNode connectivity */
struct FaceElemConn
{
Index_t lxim, lxip, letam, letap, lzetam, lzetap;
};
Kokkos::vector<FaceElemConn> m_faceToElem; /* element conn across faces */
Kokkos::vector<Int_t> m_elemBC; /* symmetry/free-surface flags for each elem face */
Kokkos::vector<Real_t> m_dxx; /* principal strains -- temporary */
Kokkos::vector<Real_t> m_dyy;
Kokkos::vector<Real_t> m_dzz;
Kokkos::vector<Real_t> m_delv_xi; /* velocity gradient -- temporary */
Kokkos::vector<Real_t> m_delv_eta;
Kokkos::vector<Real_t> m_delv_zeta;
Kokkos::vector<Real_t> m_delx_xi; /* coordinate gradient -- temporary */
Kokkos::vector<Real_t> m_delx_eta;
Kokkos::vector<Real_t> m_delx_zeta;
Kokkos::vector<Real_t> m_e; /* energy */
struct Pcomponents
{
Real_t p, q;
};
Kokkos::vector<Pcomponents> m_pq; /* pressure and artificial viscosity */
struct Qcomponents
{
Real_t ql, qq;
};
Kokkos::vector<Qcomponents> m_qlqq; /* linear and quadratic terms for q */
struct Volume
{
Real_t v, volo;
};
Kokkos::vector<Volume> m_vol; /* relative and reference volume */
Kokkos::vector<Real_t> m_vnew; /* new relative volume -- temporary */
Kokkos::vector<Real_t> m_delv; /* m_vnew - m_v */
Kokkos::vector<Real_t> m_vdov; /* volume derivative over volume */
Kokkos::vector<Real_t> m_arealg; /* characteristic length of an element */
Kokkos::vector<Real_t> m_ss; /* "sound speed" */
Kokkos::vector<Real_t> m_elemMass; /* mass */
// Cutoffs (treat as constants)
const Real_t m_e_cut; // energy tolerance
const Real_t m_p_cut; // pressure tolerance
const Real_t m_q_cut; // q tolerance
const Real_t m_v_cut; // relative volume tolerance
const Real_t m_u_cut; // velocity tolerance
// Other constants (usually setable, but hardcoded in this proxy app)
const Real_t m_hgcoef; // hourglass control
const Real_t m_ss4o3;
const Real_t m_qstop; // excessive q indicator
const Real_t m_monoq_max_slope;
const Real_t m_monoq_limiter_mult;
const Real_t m_qlc_monoq; // linear term coef for q
const Real_t m_qqc_monoq; // quadratic term coef for q
const Real_t m_qqc;
const Real_t m_eosvmax;
const Real_t m_eosvmin;
const Real_t m_pmin; // pressure floor
const Real_t m_emin; // energy floor
const Real_t m_dvovmax; // maximum allowable volume change
const Real_t m_refdens; // reference density
// Variables to keep track of timestep, simulation time, and cycle
Real_t m_dtcourant; // courant constraint
Real_t m_dthydro; // volume change constraint
Int_t m_cycle; // iteration count for simulation
Real_t m_dtfixed; // fixed time increment
Real_t m_time; // current time
Real_t m_deltatime; // variable time increment
Real_t m_deltatimemultlb;
Real_t m_deltatimemultub;
Real_t m_dtmax; // maximum allowable time increment
Real_t m_stoptime; // end time for simulation
Int_t m_numRanks;
Index_t m_colLoc;
Index_t m_rowLoc;
Index_t m_planeLoc;
Index_t m_tp;
Index_t m_sizeX;
Index_t m_sizeY;
Index_t m_sizeZ;
Index_t m_numElem;
Index_t m_numNode;
Index_t m_maxPlaneSize;
Index_t m_maxEdgeSize;
// OMP hack
Index_t* m_nodeElemStart;
Index_t* m_nodeElemCornerList;
// Used in setup
Index_t m_rowMin, m_rowMax;
Index_t m_colMin, m_colMax;
Index_t m_planeMin, m_planeMax;
};
typedef Real_t& (Domain::*Domain_member)(Index_t);
struct cmdLineOpts
{
Int_t its; // -i
Int_t nx; // -s
Int_t numReg; // -r
Int_t numFiles; // -f
Int_t showProg; // -p
Int_t quiet; // -q
Int_t viz; // -v
Int_t cost; // -c
Int_t balance; // -b
};
// Function Prototypes
// lulesh-par
Real_t
CalcElemVolume(const Real_t x[8], const Real_t y[8], const Real_t z[8]);
// lulesh-util
void
ParseCommandLineOptions(int argc, char* argv[], Int_t myRank, struct cmdLineOpts* opts);
void
VerifyAndWriteFinalOutput(Real_t elapsed_time, Domain& locDom, Int_t nx, Int_t numRanks);
// lulesh-viz
void
DumpToVisit(Domain& domain, int numFiles, int myRank, int numRanks);
// lulesh-comm
void
CommRecv(Domain& domain, Int_t msgType, Index_t xferFields, Index_t dx, Index_t dy,
Index_t dz, bool doRecv, bool planeOnly);
void
CommSend(Domain& domain, Int_t msgType, Index_t xferFields, Domain_member* fieldData,
Index_t dx, Index_t dy, Index_t dz, bool doSend, bool planeOnly);
void
CommSBN(Domain& domain, Int_t xferFields, Domain_member* fieldData);
void
CommSyncPosVel(Domain& domain);
void
CommMonoQ(Domain& domain);
// lulesh-init
void
InitMeshDecomp(Int_t numRanks, Int_t myRank, Int_t* col, Int_t* row, Int_t* plane,
Int_t* side);
+1 -1
Dosyayı Görüntüle
@@ -1,4 +1,4 @@
cmake_minimum_required(VERSION 3.13 FATAL_ERROR)
cmake_minimum_required(VERSION 3.15 FATAL_ERROR)
project(omnitrace-parallel-overhead LANGUAGES CXX)
+5 -1
Dosyayı Görüntüle
@@ -36,6 +36,9 @@ main(int argc, char** argv)
if(argc > 2) nthread = atol(argv[2]);
if(argc > 3) nitr = atol(argv[3]);
printf("[%s] Threads: %zu\n[%s] Iterations: %zu\n[%s] fibonacci(%li)...\n", argv[0],
nthread, argv[0], nitr, argv[0], nfib);
std::vector<std::thread> threads{};
for(size_t i = 0; i < nthread; ++i)
{
@@ -43,10 +46,11 @@ main(int argc, char** argv)
threads.emplace_back(&run, _nitr, nfib);
}
run(nitr - 0.25 * nitr, nfib - 0.1 * nfib);
for(auto& itr : threads)
itr.join();
printf("fibonacci(%li) x %lu = %li\n", nfib, nthread, total.load());
printf("[%s] fibonacci(%li) x %lu = %li\n", argv[0], nfib, nthread, total.load());
return 0;
}
+1 -1
Dosyayı Görüntüle
@@ -1,4 +1,4 @@
cmake_minimum_required(VERSION 3.13 FATAL_ERROR)
cmake_minimum_required(VERSION 3.15 FATAL_ERROR)
project(omnitrace-transpose LANGUAGES CXX)
+54 -31
Dosyayı Görüntüle
@@ -21,6 +21,7 @@ THE SOFTWARE.
*/
#include "hip/hip_runtime.h"
#include <cfloat>
#include <chrono>
#include <cmath>
@@ -29,14 +30,19 @@ THE SOFTWARE.
#include <fstream>
#include <iomanip>
#include <iostream>
#include <mutex>
#include <thread>
#include <vector>
static std::mutex print_lock{};
using auto_lock_t = std::unique_lock<std::mutex>;
#define HIP_API_CALL(CALL) \
{ \
hipError_t error_ = (CALL); \
if(error_ != hipSuccess) \
{ \
auto_lock_t _lk{ print_lock }; \
fprintf(stderr, "%s:%d :: HIP error : %s\n", __FILE__, __LINE__, \
hipGetErrorString(error_)); \
exit(EXIT_FAILURE); \
@@ -49,6 +55,7 @@ check_hip_error(void)
hipError_t err = hipGetLastError();
if(err != hipSuccess)
{
auto_lock_t _lk{ print_lock };
std::cerr << "Error: " << hipGetErrorString(err) << std::endl;
exit(err);
}
@@ -63,6 +70,7 @@ verify(int* in, int* out, int M, int N)
int col = rand() % N;
if(in[row * N + col] != out[col * M + row])
{
auto_lock_t _lk{ print_lock };
std::cout << "mismatch: " << row << ", " << col << " : " << in[row * N + col]
<< " | " << out[col * M + row] << "\n";
}
@@ -85,19 +93,23 @@ transpose_a(int* in, int* out, int M, int N)
}
void
run(int rank, int argc, char** argv)
run(int rank, int tid, hipStream_t stream, int argc, char** argv)
{
(void) argc;
(void) argv;
unsigned int M = 4960 * 2;
unsigned int N = 4960 * 2;
size_t nitr = 5000;
unsigned int M = 4960 * 2;
unsigned int N = 4960 * 2;
if(argc > 2) nitr = atoll(argv[2]);
auto_lock_t _lk{ print_lock };
std::cout << "[" << rank << "][" << tid << "] M: " << M << " N: " << N << std::endl;
_lk.unlock();
std::cout << "[" << rank << "] M: " << M << " N: " << N << std::endl;
size_t size = sizeof(int) * M * N;
int* matrix = (int*) malloc(size);
int* matrix = new int[size];
for(size_t i = 0; i < M * N; i++)
matrix[i] = rand() % 1002;
int *in, *out;
int* in = nullptr;
int* out = nullptr;
std::chrono::high_resolution_clock::time_point t1, t2;
@@ -106,37 +118,36 @@ run(int rank, int argc, char** argv)
HIP_API_CALL(hipMemset(in, 0, size));
HIP_API_CALL(hipMemset(out, 0, size));
HIP_API_CALL(hipMemcpy(in, matrix, size, hipMemcpyHostToDevice));
HIP_API_CALL(hipDeviceSynchronize());
hipDeviceProp_t props;
HIP_API_CALL(hipGetDeviceProperties(&props, 0));
dim3 grid(M / 32, N / 32, 1);
dim3 block(32, 32, 1); // transpose_a
t1 = std::chrono::high_resolution_clock::now();
const unsigned times = 10000;
auto _func = [&](hipStream_t stream) {
for(size_t i = 0; i < times / 2; i++)
{
transpose_a<<<grid, block, 0, stream>>>(in, out, M, N);
check_hip_error();
}
HIP_API_CALL(hipStreamSynchronize(stream));
};
hipStream_t _stream{};
HIP_API_CALL(hipStreamCreate(&_stream));
std::thread _t{ _func, _stream };
_t.join();
_func(0);
HIP_API_CALL(hipDeviceSynchronize());
t1 = std::chrono::high_resolution_clock::now();
for(size_t i = 0; i < nitr; i++)
{
transpose_a<<<grid, block, 0, stream>>>(in, out, M, N);
check_hip_error();
}
HIP_API_CALL(hipStreamSynchronize(stream));
t2 = std::chrono::high_resolution_clock::now();
double time =
std::chrono::duration_cast<std::chrono::duration<double>>(t2 - t1).count();
float GB = (float) size * times * 2 / (1 << 30);
std::cout << "[" << rank << "] Runtime of transpose is " << time << " sec\n"
float GB = (float) size * nitr * 2 / (1 << 30);
print_lock.lock();
std::cout << "[" << rank << "][" << tid << "] Runtime of transpose is " << time
<< " sec\n"
<< "The average performance of transpose is " << GB / time << " GBytes/sec"
<< std::endl;
print_lock.unlock();
int* out_matrix = (int*) malloc(size);
HIP_API_CALL(hipDeviceSynchronize());
int* out_matrix = new int[size];
HIP_API_CALL(hipMemcpy(out_matrix, out, size, hipMemcpyDeviceToHost));
// cpu_transpose(matrix, out_matrix, M, N);
@@ -145,8 +156,8 @@ run(int rank, int argc, char** argv)
HIP_API_CALL(hipFree(in));
HIP_API_CALL(hipFree(out));
free(matrix);
free(out_matrix);
delete[] matrix;
delete[] out_matrix;
}
#if defined(USE_MPI)
@@ -174,12 +185,16 @@ main(int argc, char** argv)
int rank = 0;
int size = 1;
int nthreads = 2;
int nitr = 5000;
if(argc > 1) nthreads = atoi(argv[1]);
if(argc > 2) nitr = atoi(argv[2]);
#if defined(USE_MPI)
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
#else
(void) size;
#endif
// this is a temporary workaround in omnitrace when HIP + MPI is enabled
int ndevice = 0;
@@ -193,16 +208,24 @@ main(int argc, char** argv)
if(rank == devid && rank < ndevice)
{
std::vector<std::thread> _threads{};
std::vector<hipStream_t> _streams(nthreads);
for(int i = 0; i < nthreads; ++i)
HIP_API_CALL(hipStreamCreate(&_streams.at(i)));
for(int i = 1; i < nthreads; ++i)
_threads.emplace_back(run, rank, argc, argv);
run(rank, argc, argv);
_threads.emplace_back(run, rank, i, _streams.at(i), argc, argv);
run(rank, 0, _streams.at(0), argc, argv);
for(auto& itr : _threads)
itr.join();
for(int i = 0; i < nthreads; ++i)
HIP_API_CALL(hipStreamDestroy(_streams.at(i)));
}
#if defined(USE_MPI)
MPI_Barrier(MPI_COMM_WORLD);
do_a2a(rank);
MPI_Finalize();
#endif
HIP_API_CALL(hipDeviceSynchronize());
HIP_API_CALL(hipDeviceReset());
return 0;
}