diff --git a/.github/workflows/continuousIntegration.yml b/.github/workflows/continuousIntegration.yml index c6e13bf..326f282 100644 --- a/.github/workflows/continuousIntegration.yml +++ b/.github/workflows/continuousIntegration.yml @@ -80,12 +80,14 @@ jobs: lfs: false - name: code_checks env: - DOCKER_REPOSITORY: geosx/ubuntu:22.04-${{ needs.get_docker_image_tag_hash.outputs.DOCKER_IMAGE_TAG_HASH }} + # Cppcheck 2.13 and LLVM 20 support the project's required C++20 syntax. + DOCKER_REPOSITORY: geosx/ubuntu:24.04-clang20-1c3c049b3f629d9d44838656fd306b2a0c04c9e8 HOST_CONFIG: hostconfigs/environment.cmake - CMAKE_CXX_COMPILER: /usr/bin/clang++ - CMAKE_C_COMPILER: /usr/bin/clang + CMAKE_CXX_COMPILER: /usr/lib/llvm-20/bin/clang++ + CMAKE_C_COMPILER: /usr/lib/llvm-20/bin/clang CMAKE_BUILD_TYPE: Release BUILD_AND_TEST_ARGS: --code-checks + DOCKER_RUN_ARGS: --env PATH=/usr/lib/llvm-20/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin run: ./scripts/ci_build_and_test.sh linux_builds: @@ -154,51 +156,7 @@ jobs: CMAKE_C_COMPILER: /usr/bin/gcc CMAKE_BUILD_TYPE: Debug - - name: ubuntu22-gcc11-cuda11-rel (ubuntu-22, gcc 11.4.0, cuda-11.8.0, Release) - DOCKER_REPOSITORY: geosx/ubuntu:22.04-cuda11.8 - RUNS_ON: streak2 - CMAKE_CXX_COMPILER: /usr/bin/g++ - CMAKE_C_COMPILER: /usr/bin/gcc - CMAKE_BUILD_TYPE: Release - ENABLE_CUDA: ON - CMAKE_CUDA_ARCHITECTURES: "86" - NPROC: 1 - DOCKER_RUN_ARGS: "--cpus=1 --memory=64g --runtime=nvidia --gpus all" - - - name: ubuntu22-gcc11-cuda11-dbg (ubuntu-22, gcc 11.4.0, cuda-11.8.0, Debug) - DOCKER_REPOSITORY: geosx/ubuntu:22.04-cuda11.8 - RUNS_ON: streak2 - CMAKE_CXX_COMPILER: /usr/bin/g++ - CMAKE_C_COMPILER: /usr/bin/gcc - CMAKE_BUILD_TYPE: Debug - ENABLE_CUDA: ON - CMAKE_CUDA_ARCHITECTURES: "86" - NPROC: 1 - DOCKER_RUN_ARGS: "--cpus=1 --memory=64g --runtime=nvidia --gpus all" - - - name: ubuntu22-clang14-cuda11-rel (ubuntu-22, clang-14.0.0, cuda-11.8.0, Release) - DOCKER_REPOSITORY: geosx/ubuntu:22.04-cuda11.8 - RUNS_ON: streak2 - CMAKE_CXX_COMPILER: /usr/bin/clang++ - CMAKE_C_COMPILER: /usr/bin/clang - CMAKE_BUILD_TYPE: Release - ENABLE_CUDA: ON - CMAKE_CUDA_ARCHITECTURES: "86" - NPROC: 1 - DOCKER_RUN_ARGS: "--cpus=1 --memory=64g --runtime=nvidia --gpus all" - - - name: ubuntu22-clang14-cuda11-dbg (ubuntu-22, clang-14.0.0, cuda-11.8.0, Debug) - DOCKER_REPOSITORY: geosx/ubuntu:22.04-cuda11.8 - RUNS_ON: streak2 - CMAKE_CXX_COMPILER: /usr/bin/clang++ - CMAKE_C_COMPILER: /usr/bin/clang - CMAKE_BUILD_TYPE: Debug - ENABLE_CUDA: ON - CMAKE_CUDA_ARCHITECTURES: "86" - NPROC: 1 - DOCKER_RUN_ARGS: "--cpus=1 --memory=64g --runtime=nvidia --gpus all" - - - name: RHEL8.10-clang17-cuda12-rel (ubi8.10, clang-17.0.6, cuda-12.4.1, Release) + - name: RHEL8.10-clang17-cuda12-rel (ubi8.10, clang-17.0.6, cuda-12.8, Release) DOCKER_REPOSITORY: geosx/ubi:8.10-cuda12.8 RUNS_ON: streak2 CMAKE_CXX_COMPILER: clang++ @@ -209,7 +167,7 @@ jobs: NPROC: 1 DOCKER_RUN_ARGS: "--cpus=1 --memory=64g --runtime=nvidia --gpus all" - - name: RHEL8.10-gcc13-cuda12-rel (ubi8.10, gcc 13.2.1, cuda-12.4.1, Release) + - name: RHEL8.10-gcc13-cuda12-rel (ubi8.10, gcc 13.2.1, cuda-12.8, Release) DOCKER_REPOSITORY: geosx/ubi:8.10-cuda12.8 RUNS_ON: streak2 CMAKE_CXX_COMPILER: /opt/rh/gcc-toolset-13/root/bin/g++ @@ -293,4 +251,4 @@ jobs: run: "true" - name: PR Failure if: ${{ contains(env.RETURN_VAL, 'false') }} - run: "false" \ No newline at end of file + run: "false" diff --git a/CMakeLists.txt b/CMakeLists.txt index ef4ae26..d68e971 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -46,11 +46,11 @@ target_include_directories( Shiva INTERFACE "$" # build/include/shiva/ShivaConfig.hpp "$" ) -target_compile_features(Shiva INTERFACE cxx_std_17) +target_compile_features(Shiva INTERFACE cxx_std_20 cuda_std_20) # -------- BLT (optional; usable as top-level or submodule) --------------------- if( SHIVA_ENABLE_BLT ) - set( BLT_CXX_STD "c++17" ) + set( BLT_CXX_STD "c++20" ) # Load BLT if Shiva is standalone if( SHIVA_IS_TOPLEVEL ) set(BLT_SOURCE_DIR "${PROJECT_SOURCE_DIR}/cmake/blt" CACHE PATH "Path to BLT") diff --git a/README.md b/README.md index 6066656..b03f30c 100644 --- a/README.md +++ b/README.md @@ -4,10 +4,20 @@ # Shiva Shiva Discretization Library +Requirements: a C++20 compiler. CUDA builds require CUDA 12 or newer. + +HIP consumers must explicitly require C++20 on each target that compiles HIP sources using Shiva. The Shiva interface propagates the C++ and CUDA standards; HIP standard propagation is unavailable at the project's minimum CMake version. Replace `my_target` with your consuming target: + +```cmake +set_target_properties(my_target PROPERTIES + HIP_STANDARD 20 + HIP_STANDARD_REQUIRED ON) +``` + # Project Goals The goals are: - Provide portable device callable functionality to enable the application of discrete numerical methods such as the Finite Element Method, and Finite Difference Method, etc. - Provide compile time interface to functionality through use of templates and `constexpr` for the special case where compile time knowledge is available. - Specialize for common use cases to achieve maximum performance. For example a Hexahedral element with 1st order Lagrange shape functions and Gauss-Legendre quadrature. -- Support wide range of base geometries (tet, prism, pyramid, hex, polyhedral, etc), and discretization methods ( FEM, FVM, MPM, etc. ) \ No newline at end of file +- Support wide range of base geometries (tet, prism, pyramid, hex, polyhedral, etc), and discretization methods ( FEM, FVM, MPM, etc. ) diff --git a/cmake/Macros.cmake b/cmake/Macros.cmake index 7fc55ec..7ed6605 100644 --- a/cmake/Macros.cmake +++ b/cmake/Macros.cmake @@ -32,19 +32,23 @@ macro(shiva_add_code_checks) endforeach() endif() - set( CPPCHECK_FLAGS --std=c++17 - --enable=all + set( CPPCHECK_FLAGS --std=c++20 + --enable=all --quiet + # Use the existing SHIVA_HAS_BUILTIN fallback during analysis. + -U__has_builtin + # Omit the informational summary: testCppCheck rejects any stderr. + --suppress=checkersReport --suppress=missingIncludeSystem --suppress=unmatchedSuppression - --suppress=missingInclude - --suppress=noConstructor + --suppress=missingInclude + --suppress=noConstructor --suppress=noExplicitConstructor - --suppress=unusedFunction - --suppress=constStatement + --suppress=unusedFunction + --suppress=constStatement --suppress=unusedStructMember --suppress=unknownMacro ) - + blt_add_code_checks( PREFIX ${arg_PREFIX} SOURCES ${_sources} UNCRUSTIFY_CFG_FILE ${PROJECT_SOURCE_DIR}/src/uncrustify.cfg @@ -71,4 +75,4 @@ macro(shiva_add_code_checks) RUNNER ctest -E 'blt_gtest_smoke|testCppCheck|testClangTidy|testUncrustifyCheck|testDoxygenCheck|testCppCheck|testClangTidy' SOURCE_DIRECTORIES ${PROJECT_SOURCE_DIR}/include ) endif() -endmacro(shiva_add_code_checks) \ No newline at end of file +endmacro(shiva_add_code_checks) diff --git a/hostconfigs/LLNL/lassen-clang-13-cuda-11.cmake b/hostconfigs/LLNL/lassen-clang-13-cuda-11.cmake deleted file mode 100644 index d5221c2..0000000 --- a/hostconfigs/LLNL/lassen-clang-13-cuda-11.cmake +++ /dev/null @@ -1,12 +0,0 @@ -set(CONFIG_NAME "lassen-clang-13-cuda-11" CACHE PATH "") - -set(COMPILER_DIR /usr/tce/packages/clang/clang-13.0.1-gcc-8.3.1 ) -set(CMAKE_C_COMPILER ${COMPILER_DIR}/bin/clang CACHE PATH "") -set(CMAKE_CXX_COMPILER ${COMPILER_DIR}/bin/clang++ CACHE PATH "") - -# C++ options -set(CMAKE_CXX_FLAGS_RELEASE "-O3 -DNDEBUG -mcpu=powerpc64le -mtune=powerpc64le" CACHE STRING "") -set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-g ${CMAKE_CXX_FLAGS_RELEASE}" CACHE STRING "") -set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g" CACHE STRING "") - -include(${CMAKE_CURRENT_LIST_DIR}/lassen-cuda-11-base.cmake) \ No newline at end of file diff --git a/hostconfigs/LLNL/lassen-cuda-11-base.cmake b/hostconfigs/LLNL/lassen-cuda-11-base.cmake deleted file mode 100644 index 5448cb4..0000000 --- a/hostconfigs/LLNL/lassen-cuda-11-base.cmake +++ /dev/null @@ -1,13 +0,0 @@ - -# Cuda options -set(ENABLE_CUDA ON CACHE BOOL "") -set(CUDA_TOOLKIT_ROOT_DIR /usr/tce/packages/cuda/cuda-11.8.0 CACHE STRING "") -set(CMAKE_CUDA_HOST_COMPILER ${CMAKE_CXX_COMPILER} CACHE STRING "") -set(CMAKE_CUDA_COMPILER ${CUDA_TOOLKIT_ROOT_DIR}/bin/nvcc CACHE STRING "") -set(CMAKE_CUDA_ARCHITECTURES 70 CACHE STRING "") -set(CMAKE_CUDA_STANDARD 17 CACHE STRING "") -set(CMAKE_CUDA_FLAGS "-restrict --expt-extended-lambda --expt-relaxed-constexpr -Werror cross-execution-space-call,reorder,deprecated-declarations" CACHE STRING "") -#set(CMAKE_CUDA_FLAGS_RELEASE "-O3 -DNDEBUG -Xcompiler -DNDEBUG -Xcompiler -O3 -Xcompiler -mcpu=powerpc64le -Xcompiler -mtune=powerpc64le" CACHE STRING "") -#set(CMAKE_CUDA_FLAGS_RELWITHDEBINFO "-g -lineinfo ${CMAKE_CUDA_FLAGS_RELEASE}" CACHE STRING "") -#set(CMAKE_CUDA_FLAGS_DEBUG "-g -G -O0 -Xcompiler -O0" CACHE STRING "") - diff --git a/hostconfigs/LLNL/lassen-gcc-8-cuda-11.cmake b/hostconfigs/LLNL/lassen-gcc-8-cuda-11.cmake deleted file mode 100644 index 9a46d8b..0000000 --- a/hostconfigs/LLNL/lassen-gcc-8-cuda-11.cmake +++ /dev/null @@ -1,12 +0,0 @@ -set(CONFIG_NAME "lassen-gcc-8-cuda-11" CACHE PATH "") - -set(COMPILER_DIR /usr/tce/packages/gcc/gcc-8.3.1) -set(CMAKE_C_COMPILER ${COMPILER_DIR}/bin/gcc CACHE PATH "") -set(CMAKE_CXX_COMPILER ${COMPILER_DIR}/bin/g++ CACHE PATH "") - -# C++ options -set(CMAKE_CXX_FLAGS_RELEASE "-O3 -DNDEBUG -mcpu=power9 -mtune=power9" CACHE STRING "") -set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-g ${CMAKE_CXX_FLAGS_RELEASE}" CACHE STRING "") -set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g" CACHE STRING "") - -include(${CMAKE_CURRENT_LIST_DIR}/lassen-cuda-11-base.cmake) diff --git a/hostconfigs/TTE/maple_rocky9.cmake b/hostconfigs/TTE/maple_rocky9.cmake index e3f4a83..849cf8d 100644 --- a/hostconfigs/TTE/maple_rocky9.cmake +++ b/hostconfigs/TTE/maple_rocky9.cmake @@ -15,7 +15,7 @@ set(CUDA_TOOLKIT_ROOT_DIR /hrtc/apps/cuda/12.6.20/aarch64/rocky9 CACHE STRING "" set(CMAKE_CUDA_HOST_COMPILER ${CMAKE_CXX_COMPILER} CACHE STRING "") set(CMAKE_CUDA_COMPILER ${CUDA_TOOLKIT_ROOT_DIR}/bin/nvcc CACHE STRING "") set(CMAKE_CUDA_ARCHITECTURES 90 CACHE STRING "") -set(CMAKE_CUDA_STANDARD 17 CACHE STRING "") +set(CMAKE_CUDA_STANDARD 20 CACHE STRING "") set(CMAKE_CUDA_FLAGS "-restrict --expt-extended-lambda --expt-relaxed-constexpr -Werror cross-execution-space-call,reorder,deprecated-declarations" CACHE STRING "") #set(CMAKE_CUDA_FLAGS_RELEASE "-O3 -DNDEBUG -Xcompiler -DNDEBUG -Xcompiler -O3 -Xcompiler -mcpu=powerpc64le -Xcompiler -mtune=powerpc64le" CACHE STRING "") #set(CMAKE_CUDA_FLAGS_RELWITHDEBINFO "-g -lineinfo ${CMAKE_CUDA_FLAGS_RELEASE}" CACHE STRING "") diff --git a/include/shiva/common/ShivaMacros.hpp b/include/shiva/common/ShivaMacros.hpp index add0299..e43753a 100644 --- a/include/shiva/common/ShivaMacros.hpp +++ b/include/shiva/common/ShivaMacros.hpp @@ -120,7 +120,7 @@ void i_g_n_o_r_e( ARGS const & ... ) {} * @brief Define SHIVA_IS_CONST_EVAL() depending on compiler/toolchain */ #if defined(__CUDA_ARCH__) -// Device code (nvcc, hipcc): no support in C++17 +// Keep CUDA device checks on the runtime path. #define SHIVA_IS_CONST_EVAL() (false) #elif SHIVA_HAS_BUILTIN( __builtin_is_constant_evaluated ) diff --git a/include/shiva/functions/bases/BasisProduct.hpp b/include/shiva/functions/bases/BasisProduct.hpp index 2e19417..d3dad7c 100644 --- a/include/shiva/functions/bases/BasisProduct.hpp +++ b/include/shiva/functions/bases/BasisProduct.hpp @@ -70,6 +70,7 @@ struct BasisProduct * coordinate. * @tparam BASIS_FUNCTION_INDICES Pack of indices of the basis function to * evaluate in each dimension. + * @tparam COORD_TYPE Coordinate type supporting indexing in each dimension. * @param parentCoord The parent coordinate at which to evaluate the basis * function. * @return The value of the basis function at the specified parent coordinate. @@ -79,7 +80,7 @@ struct BasisProduct * \Phi_{i_0 i_1 ... i_{(numDims-1)}}(\boldsymbol{\xi}) = * \prod_{k=0}^{(numDims-1)} \phi_{i_k}(\xi_k)\text{, where } \\ * i_j \text{is index of the basis function in the jth dimension, and - * ranges from [0...(order+1)]} + * ranges from 0 to one less than the number of support points in that dimension.} * \f] * */ @@ -89,23 +90,11 @@ struct BasisProduct { static_assert( sizeof...(BASIS_FUNCTION_INDICES) == numDims, "Wrong number of basis function indicies specified" ); - -#if __cplusplus >= 202002L - // expand pack over number of dimensions - return executeSequence< numDims >( [&]< int ... PRODUCT_TERM_INDEX > () constexpr - { - return ( BASIS_TYPE::template value< BASIS_FUNCTION_INDICES >( parentCoord[PRODUCT_TERM_INDEX] ) * ... ); - } ); -#else + // Integral-constant arguments preserve distinct coordinate indices in NVCC device code. return executeSequence< numDims >( [&] ( auto ... PRODUCT_TERM_INDEX ) constexpr { - // fold expression to multiply the value of each BASIS_TYPE in each - // dimension. In other words the fold expands on BASIS_TYPE..., - // BASIS_FUNCTION_INDICES..., and PRODUCT_TERM_INDEX... together. return ( BASIS_TYPES::template value< BASIS_FUNCTION_INDICES >( parentCoord[decltype(PRODUCT_TERM_INDEX)::value] ) * ... ); } ); - -#endif } /** @@ -113,6 +102,7 @@ struct BasisProduct * parent coordinate. * @tparam BASIS_FUNCTION_INDICES Pack of indices of the basis function to * evaluate in each dimension. + * @tparam COORD_TYPE Coordinate type supporting indexing in each dimension. * @param parentCoord The parent coordinate at which to evaluate the basis * function gradient. * @return The gradient of the basis function at the specified parent @@ -131,14 +121,12 @@ struct BasisProduct gradient( COORD_TYPE const & parentCoord ) { static_assert( sizeof...(BASIS_FUNCTION_INDICES) == numDims, "Wrong number of basis function indicies specified" ); - -#if __cplusplus >= 202002L return executeSequence< numDims >( [&]< int ... i > () constexpr->CArrayNd< RealType, numDims > { auto gradientComponent = [&] ( auto const iGrad, auto const ... PRODUCT_TERM_INDICES ) constexpr { - // Ca + // Multiply the directional gradient by the values in every other dimension. return ( gradientComponentHelper< BASIS_TYPES, decltype(iGrad)::value, BASIS_FUNCTION_INDICES, @@ -147,46 +135,26 @@ struct BasisProduct return { (executeSequence< numDims >( gradientComponent, std::integral_constant< int, i >{} ) )... }; } ); -#else - // Expand over the dimensions. - return executeSequence< numDims >( [&] ( auto ... a ) constexpr->CArrayNd< RealType, numDims > - { - // define a lambda that calculates the gradient of the basis function in - // a single dimension/direction. - auto gradientComponent = [&] ( auto GRADIENT_COMPONENT, auto ... PRODUCT_TERM_INDICES ) constexpr - { - // fold expression calling gradientComponentHelper using expanding on - // BASIS_TYPE, BASIS_FUNCTION_INDICES, and PRODUCT_TERM_INDICES. - return ( gradientComponentHelper< BASIS_TYPES, - decltype(GRADIENT_COMPONENT)::value, - BASIS_FUNCTION_INDICES, - decltype(PRODUCT_TERM_INDICES)::value >( parentCoord ) * ... ); - }; - - // execute the gradientComponent lambda on each direction, expand the - // pack on "i" corresponding to each direction of the gradient. - return { (executeSequence< numDims >( gradientComponent, a ) )... }; - } ); -#endif } private: /** * @brief Helper function to return the gradient of a basis function, or the - * value of the basis function depending on whether or not the index of the - * basis function (@p BASIS_FUNCTION) matches the index of the direction - * component index (@p GRADIENT_COMPONENT). This filter is illustrated in the + * value of the basis function depending on whether the coordinate dimension + * (@p COORD_INDEX) matches the gradient direction (@p GRADIENT_COMPONENT). + * This filter is illustrated in the * documentation for the gradient function. * @tparam BASIS_FUNCTION The basis function type * @tparam GRADIENT_COMPONENT The dimension component of the gradient. * @tparam BASIS_FUNCTION_INDEX The index of the basis function that is being * evaluated. * @tparam COORD_INDEX The dimension component of the coordinate. + * @tparam COORD_TYPE Coordinate type supporting indexing in each dimension. * @param parentCoord The parent coordinate at which to evaluate the basis * function gradient. - * @return The gradient component of the basis function at the specified - * parent coordinate. + * @return The one-dimensional basis derivative when the dimensions match, + * otherwise the one-dimensional basis value. */ template< typename BASIS_FUNCTION, int GRADIENT_COMPONENT, int BASIS_FUNCTION_INDEX, int COORD_INDEX, typename COORD_TYPE > SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE RealType diff --git a/include/shiva/functions/bases/LagrangeBasis.hpp b/include/shiva/functions/bases/LagrangeBasis.hpp index 9330885..b32ae99 100644 --- a/include/shiva/functions/bases/LagrangeBasis.hpp +++ b/include/shiva/functions/bases/LagrangeBasis.hpp @@ -38,9 +38,6 @@ namespace functions * @tparam ORDER The order of the basis function * @tparam SPACING_TYPE The spacing type to define the interpolation points for * the Lagrange polynomials. - * @tparam USE_FOR_SEQUENCE If true, the staticFor will be used to calculate - * values and gradient functions. If false, the executeSequence function will be - * used. * * The equation for a Lagrange interpolating basis function is: * @@ -81,19 +78,12 @@ class LagrangeBasis : public SPACING_TYPE< REAL_TYPE, ORDER + 1 > SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE REAL_TYPE value( REAL_TYPE const & coord ) { -#if __cplusplus >= 202002L return executeSequence< numSupportPoints >( [&]< int ... a > () constexpr { // return fold expression that is the product of all the polynomial // factor terms. return ( valueProductTerm< BF_INDEX, a >( coord ) * ... ); } ); -#else - return executeSequence< numSupportPoints >( [&] ( auto const ... a ) constexpr - { - return ( valueProductTerm< BF_INDEX, decltype(a)::value >( coord ) * ... ); - } ); -#endif } /** @@ -116,30 +106,17 @@ class LagrangeBasis : public SPACING_TYPE< REAL_TYPE, ORDER + 1 > gradient( REAL_TYPE const & coord ) { -#if __cplusplus >= 202002L return executeSequence< numSupportPoints >( [&coord]< int ... a > () constexpr { - auto func = [&coord]< int ... b > ( auto aa ) constexpr + REAL_TYPE const values[ numSupportPoints ] = { valueProductTerm< BF_INDEX, a >( coord )... }; + auto func = [&values]< int ... b > ( auto aa ) constexpr { constexpr int aVal = decltype(aa)::value; - return gradientOfValueTerm< BF_INDEX, aVal >() * ( valueProductFactor< BF_INDEX, b, aVal >( coord ) * ... ); + return gradientOfValueTerm< BF_INDEX, aVal >() * ( valueProductFactor< b, aVal >( values ) * ... ); }; return ( executeSequence< numSupportPoints >( func, std::integral_constant< int, a >{} ) + ... ); } ); -#else - return executeSequence< numSupportPoints >( [&coord] ( auto const ... a ) constexpr - { - REAL_TYPE const values[ numSupportPoints ] = { valueProductTerm< BF_INDEX, decltype(a)::value >( coord )... }; - auto func = [&values] ( auto aa, auto ... b ) constexpr - { - constexpr int aVal = decltype(aa)::value; - return gradientOfValueTerm< BF_INDEX, aVal >() * ( valueProductFactor< decltype(b)::value, aVal >( values ) * ... ); - }; - - return ( executeSequence< numSupportPoints >( func, a ) + ... ); - } ); -#endif } @@ -151,10 +128,11 @@ class LagrangeBasis : public SPACING_TYPE< REAL_TYPE, ORDER + 1 > * @tparam BF_INDEX The index of the basis function to calculate the value of. * @tparam TERM_INDEX The index of the term in the sequence product. * @param coord The coordinate to calculate the value at. - * @return The value of the Lagrange polynomial basis function. + * @return One when @p TERM_INDEX equals @p BF_INDEX, otherwise the value + * of the product factor. * * In the equation for the "Lagrange basis", the "value term" is: - * \f[ \frac{ x - x_k }{ x_j - x_k } \f]. The + * \f[ \frac{ x - x_k }{ x_j - x_k } \f] for \f$ k \neq j \f$. * */ template< int BF_INDEX, int TERM_INDEX > @@ -175,18 +153,17 @@ class LagrangeBasis : public SPACING_TYPE< REAL_TYPE, ORDER + 1 > } /** - * @brief Applies an index filter to an array that contains the values of the - * Lagrange polynomial basis function product terms. + * @brief Omits the differentiated factor from the product of Lagrange basis + * terms in a product-rule summand. * @tparam TERM_INDEX The index of the factor to calculate the value of. * @tparam DERIVATIVE_INDEX The index of the derivative term that is calling * this function. * @param values The array that contains the Lagrange polynomial basis terms. - * @return The value of the Lagrange polynomial basis function. + * @return One for the differentiated factor, otherwise the stored value term. * - * In the equation for the "Lagrange basis": - * \f[ P_j(x) = \prod_{ { {k=0} \atop {k\neq j} } }^{n-1} - *\frac{x-x_k}{x_j-x_k} \f], the filter addresses the case where \f$ k = j - *\f$. + * The factor with index @p DERIVATIVE_INDEX is replaced by its derivative + * outside this product. Returning one here omits it without dividing by its + * value, which may be zero at a support point. * */ template< int TERM_INDEX, int DERIVATIVE_INDEX = -1 > @@ -208,7 +185,8 @@ class LagrangeBasis : public SPACING_TYPE< REAL_TYPE, ORDER + 1 > * function for a @p BF_INDEX. * @tparam BF_INDEX The index of the basis function. * @tparam TERM_INDEX The index of the term in the sequence product. - * @return The gradient of the Lagrange polynomial basis function. + * @return Zero when @p TERM_INDEX equals @p BF_INDEX, otherwise the + * derivative of the product factor. * */ template< int BF_INDEX, int TERM_INDEX > diff --git a/tests/common/testSequenceUtilities.cpp b/tests/common/testSequenceUtilities.cpp index c0a6ae1..3ca18f8 100644 --- a/tests/common/testSequenceUtilities.cpp +++ b/tests/common/testSequenceUtilities.cpp @@ -112,8 +112,6 @@ TEST( testSequenceUtilities, testForSequenceLambda ) -#if __cplusplus >= 202002L - void testSequenceExpansionTemplateLambdaHelper() { kernelLaunch( [] SHIVA_HOST_DEVICE () @@ -135,18 +133,18 @@ TEST( testSequenceUtilities, testSequenceExpansionTemplateLambda ) -void testSequenceExpansionTemplateLambdaHelper() +void testNestedSequenceExpansionTemplateLambdaHelper() { kernelLaunch( [] SHIVA_HOST_DEVICE () { - constexpr int staticSum0 = executeSequence< 10 >( [&]< int ... a > () constexpr + constexpr int staticSum0 = + executeSequence< 10 >( [&]< int ... a > () constexpr { return ( executeSequence< 10 > - ( - [ h = Data::h, aa = std::integral_constant< int, a >{} ]< int ... b > () constexpr + ( [ h = Data::h, aa = std::integral_constant< int, a >{} ]< int ... b > () constexpr { return ( (h[aa] * h[b]) + ...); } - ) + ... + ) + ... ); } ); static_assert( staticSum0 == Data::nested_sum_of_h ); @@ -155,7 +153,7 @@ void testSequenceExpansionTemplateLambdaHelper() TEST( testSequenceUtilities, testNestedSequenceExpansionTemplateLambda ) { - testSequenceExpansionTemplateLambdaHelper(); + testNestedSequenceExpansionTemplateLambdaHelper(); } @@ -183,7 +181,6 @@ TEST( testSequenceUtilities, testForSequenceTemplateLambda ) { testForSequenceTemplateLambdaHelper(); } -#endif int main( int argc, char * * argv ) diff --git a/tests/functions/testLagrangeBasis.cpp b/tests/functions/testLagrangeBasis.cpp index 07fcb4e..9e7bbb9 100644 --- a/tests/functions/testLagrangeBasis.cpp +++ b/tests/functions/testLagrangeBasis.cpp @@ -13,6 +13,7 @@ #include "shiva/functions/bases/LagrangeBasis.hpp" +#include "shiva/functions/bases/BasisProduct.hpp" #include "shiva/functions/spacing/Spacing.hpp" #include "shiva/common/ShivaMacros.hpp" #include "shiva/common/pmpl.hpp" @@ -20,6 +21,7 @@ #include #include +#include using namespace shiva; using namespace shiva::functions; @@ -142,6 +144,138 @@ TEST( testSpacing, testLagrangeBasisGaussLobattoSpacing ) testBasisAtRunTime< BasisHelperType >(); } +// Closed-form values and derivatives for the constant, linear, and quadratic +// bases on [-1, 1]. Unused entries pad the lower-order rows with zeros. +template< typename REAL_TYPE, int ORDER > +SHIVA_CONSTEXPR_HOSTDEVICE_FORCEINLINE CArrayNd< REAL_TYPE, 2, 3 > +lowOrderReference( REAL_TYPE const x ) +{ + constexpr REAL_TYPE zero = 0; + constexpr REAL_TYPE one = 1; + constexpr REAL_TYPE half = 0.5; + if constexpr ( ORDER == 0 ) + { + return { one, zero, zero, zero, zero, zero }; + } + else if constexpr ( ORDER == 1 ) + { + return { (1 - x) / 2, (1 + x) / 2, zero, -half, half, zero }; + } + else + { + return { x * (x - 1) / 2, 1 - x * x, x * (x + 1) / 2, + x - half, -2 * x, x + half }; + } +} + +template< typename REAL_TYPE, int ORDER, template< typename, int > typename SPACING_TYPE > +void testLowOrderBasis() +{ + using BasisType = LagrangeBasis< REAL_TYPE, ORDER, SPACING_TYPE >; + constexpr int numPoints = 4; + constexpr int stride = 1 + 2 * BasisType::numSupportPoints; + constexpr REAL_TYPE coords[numPoints] = { -1, 0, REAL_TYPE( 0.3 ), 1 }; + constexpr REAL_TYPE tolerance = 16 * std::numeric_limits< REAL_TYPE >::epsilon(); + REAL_TYPE data[numPoints * stride]{}; + for ( int point = 0; point < numPoints; ++point ) + { + data[point * stride] = coords[point]; + } + + pmpl::genericKernelWrapper( numPoints * stride, data, [] SHIVA_HOST_DEVICE ( REAL_TYPE * const kernelData ) + { + forSequence< numPoints >( [] ( auto const POINT ) constexpr + { + constexpr REAL_TYPE compileTimeCoords[numPoints] = { -1, 0, REAL_TYPE( 0.3 ), 1 }; + constexpr REAL_TYPE coord = compileTimeCoords[POINT]; + constexpr auto expected = lowOrderReference< REAL_TYPE, ORDER >( coord ); + forSequence< BasisType::numSupportPoints >( [&] ( auto const BF_INDEX ) constexpr + { + constexpr REAL_TYPE value = BasisType::template value< BF_INDEX >( coord ); + constexpr REAL_TYPE gradient = BasisType::template gradient< BF_INDEX >( coord ); + static_assert( pmpl::check( value, expected( 0, BF_INDEX ), tolerance ) ); + static_assert( pmpl::check( gradient, expected( 1, BF_INDEX ), tolerance ) ); + } ); + } ); + + for ( int point = 0; point < numPoints; ++point ) + { + int const offset = point * stride; + REAL_TYPE const coord = kernelData[offset]; + forSequence< BasisType::numSupportPoints >( [&] ( auto const BF_INDEX ) + { + kernelData[offset + 1 + BF_INDEX] = BasisType::template value< BF_INDEX >( coord ); + kernelData[offset + 1 + BasisType::numSupportPoints + BF_INDEX] = BasisType::template gradient< BF_INDEX >( coord ); + } ); + } + } ); + + for ( int point = 0; point < numPoints; ++point ) + { + auto const expected = lowOrderReference< REAL_TYPE, ORDER >( coords[point] ); + for ( int basis = 0; basis < BasisType::numSupportPoints; ++basis ) + { + EXPECT_NEAR( data[point * stride + 1 + basis], expected( 0, basis ), tolerance ); + EXPECT_NEAR( data[point * stride + 1 + BasisType::numSupportPoints + basis], expected( 1, basis ), tolerance ); + } + } +} + +TEST( testLagrangeBasis, lowOrderEqualSpacing ) +{ + testLowOrderBasis< float, 0, EqualSpacing >(); + testLowOrderBasis< double, 0, EqualSpacing >(); + testLowOrderBasis< float, 1, EqualSpacing >(); + testLowOrderBasis< double, 1, EqualSpacing >(); + testLowOrderBasis< float, 2, EqualSpacing >(); + testLowOrderBasis< double, 2, EqualSpacing >(); +} + +TEST( testLagrangeBasis, lowOrderGaussLobattoSpacing ) +{ + testLowOrderBasis< float, 1, GaussLobattoSpacing >(); + testLowOrderBasis< double, 1, GaussLobattoSpacing >(); + testLowOrderBasis< float, 2, GaussLobattoSpacing >(); + testLowOrderBasis< double, 2, GaussLobattoSpacing >(); +} + +template< typename REAL_TYPE > +void testAnisotropicBasisProduct() +{ + using Product = BasisProduct< REAL_TYPE, + LagrangeBasis< REAL_TYPE, 1, EqualSpacing >, + LagrangeBasis< REAL_TYPE, 2, GaussLobattoSpacing > >; + constexpr REAL_TYPE tolerance = std::numeric_limits< REAL_TYPE >::epsilon(); + REAL_TYPE data[3] = { REAL_TYPE( 0.5 ), REAL_TYPE( -0.25 ), 0 }; + pmpl::genericKernelWrapper( 3, data, [] SHIVA_HOST_DEVICE ( REAL_TYPE * const kernelData ) + { + // Phi_01(x,y) = (1-x)(1-y*y)/2. + // Distinct magnitudes ensure reusing a coordinate changes the basis value. + constexpr REAL_TYPE coord[2] = { REAL_TYPE( 0.5 ), REAL_TYPE( -0.25 ) }; + constexpr auto value = Product::template value< 0, 1 >( coord ); + constexpr auto gradient = Product::template gradient< 0, 1 >( coord ); + static_assert( pmpl::check( value, REAL_TYPE( 0.234375 ), tolerance ) ); + static_assert( pmpl::check( gradient( 0 ), REAL_TYPE( -0.46875 ), tolerance ) ); + static_assert( pmpl::check( gradient( 1 ), REAL_TYPE( 0.125 ), tolerance ) ); + + REAL_TYPE const runtimeCoord[2] = { kernelData[0], kernelData[1] }; + auto const runtimeGradient = Product::template gradient< 0, 1 >( runtimeCoord ); + kernelData[0] = Product::template value< 0, 1 >( runtimeCoord ); + kernelData[1] = runtimeGradient( 0 ); + kernelData[2] = runtimeGradient( 1 ); + } ); + + EXPECT_NEAR( data[0], REAL_TYPE( 0.234375 ), tolerance ); + EXPECT_NEAR( data[1], REAL_TYPE( -0.46875 ), tolerance ); + EXPECT_NEAR( data[2], REAL_TYPE( 0.125 ), tolerance ); +} + +TEST( testLagrangeBasis, anisotropicBasisProduct ) +{ + testAnisotropicBasisProduct< float >(); + testAnisotropicBasisProduct< double >(); +} + int main( int argc, char * * argv ) { ::testing::InitGoogleTest( &argc, argv ); diff --git a/tests/geometry/mapping/testScaling.cpp b/tests/geometry/mapping/testScaling.cpp index 16fd29b..942e987 100644 --- a/tests/geometry/mapping/testScaling.cpp +++ b/tests/geometry/mapping/testScaling.cpp @@ -17,7 +17,9 @@ #include "shiva/common/pmpl.hpp" #include "shiva/functions/quadrature/Quadrature.hpp" +#include #include +#include #include using namespace shiva; @@ -152,14 +154,7 @@ constexpr bool checkScalingInterfaceConstexpr() { bool checks[10]{}; checkScalingInterface< REAL_TYPE >( checks ); - for( bool const check : checks ) - { - if( !check ) - { - return false; - } - } - return true; + return std::all_of( std::begin( checks ), std::end( checks ), [] ( bool const check ) constexpr { return check; } ); } static_assert( checkScalingInterfaceConstexpr< float >() );