diff --git a/packages/stk/CHANGELOG.md b/packages/stk/CHANGELOG.md index 57f2d9d358e1..cf98fdc51eaf 100644 --- a/packages/stk/CHANGELOG.md +++ b/packages/stk/CHANGELOG.md @@ -1,5 +1,9 @@ # CHANGELOG +5.21.5 (STK_VERSION 5210500) 9/25/2024 + general: Fixed MI300A unified memory build errors. + stk_search: Turned off sorted results by default. + 5.21.4-1 (STK_VERSION 5210401) 9/04/2024 Fix cmake configuration errors that occurred on AMD MI300A platform diff --git a/packages/stk/stk_balance/stk_balance/balanceUtils.hpp b/packages/stk/stk_balance/stk_balance/balanceUtils.hpp index 99804f05c484..0f88cdf641ea 100644 --- a/packages/stk/stk_balance/stk_balance/balanceUtils.hpp +++ b/packages/stk/stk_balance/stk_balance/balanceUtils.hpp @@ -34,7 +34,6 @@ #define STK_BALANCE_UTILS #include "mpi.h" -#include #include #include #include "stk_mesh/base/Field.hpp" // for field_data @@ -44,6 +43,8 @@ namespace stk { +namespace mesh { class BulkData; } + namespace balance { diff --git a/packages/stk/stk_balance/stk_balance/internal/NodeBalancer.cpp b/packages/stk/stk_balance/stk_balance/internal/NodeBalancer.cpp index 2d0939249f08..ffa379c12733 100644 --- a/packages/stk/stk_balance/stk_balance/internal/NodeBalancer.cpp +++ b/packages/stk/stk_balance/stk_balance/internal/NodeBalancer.cpp @@ -80,10 +80,8 @@ void NodeBalancer::getGlobalLoadImbalance(double &loadFactor, int& numLocallyOwnedNodes) { stk::mesh::Selector localSelector = m_metaData.locally_owned_part(); - stk::mesh::EntityVector ownedNodes; - stk::mesh::get_entities(m_bulkData, stk::topology::NODE_RANK, localSelector, ownedNodes); + numLocallyOwnedNodes = stk::mesh::count_entities(m_bulkData, stk::topology::NODE_RANK, localSelector); - numLocallyOwnedNodes = ownedNodes.size(); int maxLocallyOwned = 0; int minLocallyOwned = 0; stk::all_reduce_max(m_bulkData.parallel(), &numLocallyOwnedNodes, &maxLocallyOwned, 1); diff --git a/packages/stk/stk_balance/stk_balance/internal/StkBalanceUtils.cpp b/packages/stk/stk_balance/stk_balance/internal/StkBalanceUtils.cpp index 52003dc033da..a03b2c69429f 100644 --- a/packages/stk/stk_balance/stk_balance/internal/StkBalanceUtils.cpp +++ b/packages/stk/stk_balance/stk_balance/internal/StkBalanceUtils.cpp @@ -175,7 +175,11 @@ SearchElemPairs getBBIntersectionsForFacesParticles(stk::mesh::BulkData& stkMesh stk::balance::internal::SearchElemPairs searchResults; - stk::search::coarse_search(faceBoxes, faceBoxes, stk::search::KDTREE, stkMeshBulkData.parallel(), searchResults); + bool enforceSearchResultSymmetry = true; + bool autoSwapDomainAndRange = true; + bool sortSearchResults = true; + stk::search::coarse_search(faceBoxes, faceBoxes, stk::search::KDTREE, stkMeshBulkData.parallel(), searchResults, + enforceSearchResultSymmetry, autoSwapDomainAndRange, sortSearchResults); stk::balance::internal::SearchElemPairs::iterator iter = std::unique(searchResults.begin(), searchResults.end()); searchResults.resize(iter - searchResults.begin()); diff --git a/packages/stk/stk_balance/stk_balance/io/BalanceIO.hpp b/packages/stk/stk_balance/stk_balance/io/BalanceIO.hpp index e8f924066471..faa5006c3e8e 100644 --- a/packages/stk/stk_balance/stk_balance/io/BalanceIO.hpp +++ b/packages/stk/stk_balance/stk_balance/io/BalanceIO.hpp @@ -37,7 +37,6 @@ #include "mpi.h" #include "stk_mesh/base/MetaData.hpp" -#include "stk_mesh/base/BulkData.hpp" #include "stk_io/StkMeshIoBroker.hpp" #include "stk_balance/balanceUtils.hpp" @@ -46,6 +45,9 @@ #include namespace stk { + +namespace mesh { class BulkData; } + namespace balance { class BalanceIO diff --git a/packages/stk/stk_balance/stk_balance/mesh/BalanceMesh.cpp b/packages/stk/stk_balance/stk_balance/mesh/BalanceMesh.cpp index 3f2c643e8be8..839097376daf 100644 --- a/packages/stk/stk_balance/stk_balance/mesh/BalanceMesh.cpp +++ b/packages/stk/stk_balance/stk_balance/mesh/BalanceMesh.cpp @@ -32,6 +32,7 @@ // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include "BalanceMesh.hpp" +#include namespace stk { namespace balance { diff --git a/packages/stk/stk_balance/stk_balance/mesh/BalanceMesh.hpp b/packages/stk/stk_balance/stk_balance/mesh/BalanceMesh.hpp index 3d74f60a2487..9195d3a04327 100644 --- a/packages/stk/stk_balance/stk_balance/mesh/BalanceMesh.hpp +++ b/packages/stk/stk_balance/stk_balance/mesh/BalanceMesh.hpp @@ -34,10 +34,11 @@ #ifndef STK_BALANCE_BALANCE_MESH_HPP #define STK_BALANCE_BALANCE_MESH_HPP -#include "stk_mesh/base/MetaData.hpp" -#include "stk_mesh/base/BulkData.hpp" - namespace stk { + +namespace mesh { class MetaData; } +namespace mesh { class BulkData; } + namespace balance { class BalanceMesh diff --git a/packages/stk/stk_doc_tests/stk_io/restartTestUtils.hpp b/packages/stk/stk_doc_tests/stk_io/restartTestUtils.hpp index 7cc7c1506278..fc03f0e4a597 100644 --- a/packages/stk/stk_doc_tests/stk_io/restartTestUtils.hpp +++ b/packages/stk/stk_doc_tests/stk_io/restartTestUtils.hpp @@ -39,11 +39,12 @@ #include #include #include -#include #include #include #include +namespace stk { namespace mesh { class BulkData; } } + inline void checkFileForNodalVarNames(const std::string &exodusFilename, const std::vector& nodalVarNames) { Ioss::DatabaseIO *iossDb = Ioss::IOFactory::create("exodus", exodusFilename, Ioss::READ_MODEL, MPI_COMM_WORLD); @@ -95,25 +96,6 @@ inline void putDataOnTestField(stk::mesh::BulkData &stkMeshBulkData, const doubl } } -//inline void putDataOnTriStateField(stk::mesh::BulkData &bulkData, stk::mesh::FieldBase *triStateField, -// const double stateNp1Value, -// const double stateNValue, -// const double stateNm1Value) -//{ -// stk::mesh::FieldBase *statedFieldNp1 = -// triStateField->field_state(stk::mesh::StateNP1); -// putDataOnTestField(bulkData, stateNp1Value, -// *statedFieldNp1); -// stk::mesh::FieldBase *statedFieldN = -// triStateField->field_state(stk::mesh::StateN); -// putDataOnTestField(bulkData, stateNValue, -// *statedFieldN); -// stk::mesh::FieldBase *statedFieldNm1 = -// triStateField->field_state(stk::mesh::StateNM1); -// putDataOnTestField(bulkData, stateNm1Value, -// *statedFieldNm1); -//} - inline void testDataOnField(stk::mesh::BulkData &stkMeshBulkData, const double goldValue, stk::mesh::FieldBase &field) { std::vector nodes; diff --git a/packages/stk/stk_doc_tests/stk_mesh/stkMeshTestUtils.hpp b/packages/stk/stk_doc_tests/stk_mesh/stkMeshTestUtils.hpp index 78ce388573d2..28ab13b9e52b 100644 --- a/packages/stk/stk_doc_tests/stk_mesh/stkMeshTestUtils.hpp +++ b/packages/stk/stk_doc_tests/stk_mesh/stkMeshTestUtils.hpp @@ -36,11 +36,12 @@ #define stkMeshTestUtilsHpp #include -#include #include #include #include +namespace stk { namespace mesh { class BulkData; } } + namespace testUtils { diff --git a/packages/stk/stk_doc_tests/stk_transfer/howToUseCopyTransfer.cpp b/packages/stk/stk_doc_tests/stk_transfer/howToUseCopyTransfer.cpp index 242a42e95ac5..6ae186c4bb98 100644 --- a/packages/stk/stk_doc_tests/stk_transfer/howToUseCopyTransfer.cpp +++ b/packages/stk/stk_doc_tests/stk_transfer/howToUseCopyTransfer.cpp @@ -72,14 +72,14 @@ void change_mesh_decomposition(stk::mesh::BulkData& mesh) int myProc = mesh.parallel_rank(); int otherProc = 1-myProc; - stk::mesh::for_each_entity_run(mesh, stk::topology::ELEM_RANK, - mesh.mesh_meta_data().locally_owned_part(), - [&entityProcPairs, &otherProc](const stk::mesh::BulkData& bulkData, const stk::mesh::Entity& elem) + stk::mesh::for_each_entity_run_no_threads(mesh, stk::topology::ELEM_RANK, + mesh.mesh_meta_data().locally_owned_part(), + [&entityProcPairs, &otherProc](const stk::mesh::BulkData& bulkData, const stk::mesh::Entity& elem) { entityProcPairs.emplace_back(elem, otherProc); }); - stk::mesh::for_each_entity_run(mesh, stk::topology::NODE_RANK, + stk::mesh::for_each_entity_run_no_threads(mesh, stk::topology::NODE_RANK, mesh.mesh_meta_data().locally_owned_part(), [&entityProcPairs, &otherProc](const stk::mesh::BulkData& bulkData, const stk::mesh::Entity& node) { diff --git a/packages/stk/stk_expreval/stk_expreval/Eval.cpp b/packages/stk/stk_expreval/stk_expreval/Eval.cpp index 6c118a2f7621..ae20873d80d0 100644 --- a/packages/stk/stk_expreval/stk_expreval/Eval.cpp +++ b/packages/stk/stk_expreval/stk_expreval/Eval.cpp @@ -377,12 +377,9 @@ Eval::initialize_function_map() m_functionMap["weibull_pdf"] = FunctionType::WEIBULL_PDF; m_functionMap["gamma_pdf"] = FunctionType::GAMMA_PDF; - m_functionMap["rand"] = FunctionType::RAND; - m_functionMap["srand"] = FunctionType::SRAND; - m_functionMap["random"] = FunctionType::RANDOM; + m_functionMap["ts_random"] = FunctionType::TS_RANDOM; m_functionMap["ts_normal"] = FunctionType::TS_NORMAL; - m_functionMap["time"] = FunctionType::TIME; } Eval & diff --git a/packages/stk/stk_expreval/stk_expreval/Function.hpp b/packages/stk/stk_expreval/stk_expreval/Function.hpp index d04d74632406..da1373e8c568 100644 --- a/packages/stk/stk_expreval/stk_expreval/Function.hpp +++ b/packages/stk/stk_expreval/stk_expreval/Function.hpp @@ -111,24 +111,12 @@ enum class FunctionType { WEIBULL_PDF, GAMMA_PDF, - RAND, - SRAND, - RANDOM, TS_RANDOM, TS_NORMAL, - TIME, UNDEFINED }; -constexpr bool is_function_supported_on_device(FunctionType type) -{ - return type != FunctionType::RAND && - type != FunctionType::SRAND && - type != FunctionType::RANDOM && - type != FunctionType::TIME && - type != FunctionType::UNDEFINED; -} KOKKOS_INLINE_FUNCTION double cycloidal_ramp(double t, double t1, double t2) diff --git a/packages/stk/stk_expreval/stk_expreval/NgpNode.hpp b/packages/stk/stk_expreval/stk_expreval/NgpNode.hpp index 8f378be36940..d57c16746808 100644 --- a/packages/stk/stk_expreval/stk_expreval/NgpNode.hpp +++ b/packages/stk/stk_expreval/stk_expreval/NgpNode.hpp @@ -627,30 +627,6 @@ class NgpNode STK_NGP_ThrowErrorMsg("Incorrect number of arguments for gamma_pdf function"); break; } - case FunctionType::RAND : { - if (argumentCount == 0) { - return real_rand(); - } - STK_NGP_ThrowErrorMsg("Incorrect number of arguments for rand function"); - break; - } - case FunctionType::SRAND : { - if (argumentCount == 1) { - return real_srand(arguments[0]); - } - STK_NGP_ThrowErrorMsg("Incorrect number of arguments for srand function"); - break; - } - case FunctionType::RANDOM : { - if (argumentCount == 0) { - return random0(); - } - else if (argumentCount == 1) { - return random1(arguments[0]); - } - STK_NGP_ThrowErrorMsg("Incorrect number of arguments for random function"); - break; - } case FunctionType::TS_RANDOM : { if (argumentCount == 4) { return time_space_random(arguments[0], arguments[1], arguments[2], arguments[3]); @@ -666,13 +642,6 @@ class NgpNode STK_NGP_ThrowErrorMsg("Incorrect number of arguments for ts_normal function"); break; } - case FunctionType::TIME : { - if (argumentCount == 0) { - return current_time(); - } - STK_NGP_ThrowErrorMsg("Incorrect number of arguments for time function"); - break; - } case FunctionType::UNDEFINED : { STK_NGP_ThrowErrorMsg("Undefined function type"); break; diff --git a/packages/stk/stk_expreval/stk_expreval/ParsedEval.hpp b/packages/stk/stk_expreval/stk_expreval/ParsedEval.hpp index 911901045355..8c6ef916911b 100644 --- a/packages/stk/stk_expreval/stk_expreval/ParsedEval.hpp +++ b/packages/stk/stk_expreval/stk_expreval/ParsedEval.hpp @@ -73,6 +73,8 @@ class ParsedEval : public ParsedEvalBase m_hostNodes(i) = NgpNode(*(eval.get_node(i))); } Kokkos::deep_copy(m_deviceNodes, m_hostNodes); + + check_for_errors(); } KOKKOS_DEFAULTED_FUNCTION ParsedEval(const ParsedEval&) = default; @@ -81,20 +83,13 @@ class ParsedEval : public ParsedEvalBase virtual int get_result_buffer_size() override { return RESULT_BUFFER_SIZE; } - void check_for_errors(bool will_run_on_device) const override + STK_DEPRECATED_MSG("check_for_errors is now called by the constructor. No need to call it yourself") + void check_for_errors(bool /*will_run_on_device*/) const override { - for (size_t i=0; i < m_hostNodes.size(); ++i) - { - const NgpNode& node = m_hostNodes(i); - if (node.m_opcode == OPCODE_FUNCTION) - { - FunctionType funcType = node.m_data.function.functionType; - STK_ThrowRequireMsg(funcType != FunctionType::UNDEFINED, "user defined functions are not supported by ParsedEval"); - STK_ThrowRequireMsg(will_run_on_device && is_function_supported_on_device(funcType), "random number generation and time functions not supported on device"); - } - } + check_for_errors(); } + KOKKOS_INLINE_FUNCTION int get_num_variables() const { return m_numVariables; } @@ -136,6 +131,19 @@ class ParsedEval : public ParsedEvalBase template friend class DeviceVariableMap; + void check_for_errors() const + { + for (size_t i=0; i < m_hostNodes.size(); ++i) + { + const NgpNode& node = m_hostNodes(i); + if (node.m_opcode == OPCODE_FUNCTION) + { + FunctionType funcType = node.m_data.function.functionType; + STK_ThrowRequireMsg(funcType != FunctionType::UNDEFINED, "user defined functions and system functions (rand(), time() etc.) are not supported by ParsedEval"); + } + } + } + int m_numVariables; int m_requiredResultBufferSize; Variable::ArrayOffset m_arrayOffsetType; diff --git a/packages/stk/stk_expreval/stk_expreval/ParsedEvalBase.hpp b/packages/stk/stk_expreval/stk_expreval/ParsedEvalBase.hpp index 5f54b52ed888..ba6fedb6bc28 100644 --- a/packages/stk/stk_expreval/stk_expreval/ParsedEvalBase.hpp +++ b/packages/stk/stk_expreval/stk_expreval/ParsedEvalBase.hpp @@ -36,6 +36,7 @@ #define PARSEDEVALBASE_HPP #include +#include "stk_util/stk_config.h" namespace stk { namespace expreval { @@ -51,6 +52,8 @@ class ParsedEvalBase virtual int get_result_buffer_size() = 0; + // remove on 10/29/2024 + STK_DEPRECATED_MSG("check_for_errors is now called by the constructor. No need to call it yourself") virtual void check_for_errors(bool will_run_on_device) const = 0; }; diff --git a/packages/stk/stk_integration_tests/cmake_install_test/build_stk_no_stk_mesh_using_cmake b/packages/stk/stk_integration_tests/cmake_install_test/build_stk_no_stk_mesh_using_cmake index a4b5e7a52dd3..09024d035cf5 100755 --- a/packages/stk/stk_integration_tests/cmake_install_test/build_stk_no_stk_mesh_using_cmake +++ b/packages/stk/stk_integration_tests/cmake_install_test/build_stk_no_stk_mesh_using_cmake @@ -64,7 +64,7 @@ if [ $? -ne 0 ] ; then fi printf "Now building trilinos/stk using make...\n"; -exe "make VERBOSE=1 -j8 >& ${stk_make_log}"; +exe "make VERBOSE=1 -j16 >& ${stk_make_log}"; if [ $? -ne 0 ] ; then echo "!! error in make, check output in ${stk_make_log} !!"; exit 1; diff --git a/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk b/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk index 4d4b62032503..e7573111ec67 100755 --- a/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk +++ b/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk @@ -77,7 +77,6 @@ cmake \ -DKokkos_ENABLE_CUDA_UVM:BOOL=OFF \ -DKokkos_ENABLE_CUDA_RELOCATABLE_DEVICE_CODE:BOOL=ON \ -DKokkos_ARCH_VOLTA70=${cuda_on_or_off} \ --DTrilinos_ENABLE_KokkosKernels:BOOL=ON \ -DTrilinos_ENABLE_Fortran:BOOL=ON \ -DCMAKE_CXX_STANDARD:STRING=17 \ -DCMAKE_CXX_FLAGS:STRING="-D${fortran_macro} ${cmake_cxx_flags} -Werror=dangling-else" \ diff --git a/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk_no_stk_mesh b/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk_no_stk_mesh index a9fbde405c0b..179f70e26365 100755 --- a/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk_no_stk_mesh +++ b/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk_no_stk_mesh @@ -34,6 +34,7 @@ cmake \ -DTPL_ENABLE_BLAS=ON \ -DTPL_ENABLE_LAPACK=ON \ -DSTK_ENABLE_TESTS:BOOL=ON \ +-DTrilinos_ENABLE_Intrepid2:BOOL=ON \ -DTrilinos_ENABLE_STK:BOOL=ON \ -DTrilinos_ENABLE_STKMesh:BOOL=OFF \ -DTrilinos_ENABLE_STKUtil:BOOL=ON \ diff --git a/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk_standalone b/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk_standalone index ba6d4d045a7c..8288e48c362a 100755 --- a/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk_standalone +++ b/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk_standalone @@ -47,7 +47,6 @@ cmake \ -DShards_DIR=/fgs/william/Trilinos/install_gcc/lib64/cmake/Shards \ -DZoltan2Core_DIR=/fgs/william/Trilinos/install_gcc/lib64/cmake/Zoltan2Core \ -DKokkos_DIR=/fgs/william/Trilinos/install_gcc/lib64/cmake/Kokkos \ --DKokkosKernels_DIR=/fgs/william/kokkos-kernels/install/lib64/cmake/KokkosKernels \ -DGTest_DIR=/fgs/william/googletest/install/lib64/cmake/GTest \ ${stk_src_dir}/ diff --git a/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk_user_facing b/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk_user_facing index 46486de3164f..4bc294d8bd15 100755 --- a/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk_user_facing +++ b/packages/stk/stk_integration_tests/cmake_install_test/run_cmake_stk_user_facing @@ -50,7 +50,6 @@ cmake \ -DKokkos_ENABLE_RELOCATABLE_DEVICE_CODE:BOOL=OFF \ -DKokkos_ARCH_VOLTA70=${cuda_on_or_off} \ -DTpetra_ENABLE_CUDA:BOOL=${cuda_on_or_off} \ --DTrilinos_ENABLE_KokkosKernels:BOOL=ON \ -DTrilinos_ENABLE_Zoltan:BOOL=ON \ -DTrilinos_ENABLE_Fortran:BOOL=OFF \ -DCMAKE_CXX_STANDARD:STRING=17 \ diff --git a/packages/stk/stk_io/stk_io/IOHelpers.hpp b/packages/stk/stk_io/stk_io/IOHelpers.hpp index dd974f0662df..c869ded988f9 100644 --- a/packages/stk/stk_io/stk_io/IOHelpers.hpp +++ b/packages/stk/stk_io/stk_io/IOHelpers.hpp @@ -43,7 +43,6 @@ #include // for DatabasePurpose #include #include // for MeshField, etc -#include // for BulkData #include // for Selector #include // for ParallelMachine #include // for Type diff --git a/packages/stk/stk_io/stk_io/IossBridge.cpp b/packages/stk/stk_io/stk_io/IossBridge.cpp index 451354b23099..57cd09b1ab59 100644 --- a/packages/stk/stk_io/stk_io/IossBridge.cpp +++ b/packages/stk/stk_io/stk_io/IossBridge.cpp @@ -2096,6 +2096,35 @@ const stk::mesh::FieldBase *declare_stk_field_internal(stk::mesh::MetaData &meta return !omitted; } + int64_t get_side_offset(const Ioss::ElementTopology* sideTopo, + const Ioss::ElementTopology* parentTopo) + { + int64_t sideOffset = 0; + if (( sideTopo != nullptr) && ( sideTopo->name() != "unknown") && + (parentTopo != nullptr) && (parentTopo->name() != "unknown")) { + int sideTopoDim = sideTopo->parametric_dimension(); + int elemTopoDim = parentTopo->parametric_dimension(); + int elemSpatDim = parentTopo->spatial_dimension(); + + if (sideTopoDim + 1 < elemSpatDim && sideTopoDim < elemTopoDim) { + sideOffset = parentTopo->number_faces(); + } + } + return sideOffset; + } + + int64_t get_side_offset(const Ioss::SideBlock* sb) + { + if(nullptr != sb) { + const Ioss::ElementTopology *sideTopo = sb->topology(); + const Ioss::ElementTopology *parentTopo = sb->parent_element_topology(); + return get_side_offset(sideTopo, parentTopo); + } + + return 0; + } + + namespace { stk::mesh::EntityRank get_output_rank(stk::io::OutputParams& params) @@ -2208,22 +2237,6 @@ const stk::mesh::FieldBase *declare_stk_field_internal(stk::mesh::MetaData &meta ioss_add_fields(part, stk::topology::NODE_RANK, ns, Ioss::Field::ATTRIBUTE); } - int64_t get_side_offset(const Ioss::ElementTopology* sideTopo, - const Ioss::ElementTopology* parentTopo) - { - int64_t sideOffset = 0; - if ((sideTopo != nullptr) && (parentTopo != nullptr)) { - int sideTopoDim = sideTopo->parametric_dimension(); - int elemTopoDim = parentTopo->parametric_dimension(); - int elemSpatDim = parentTopo->spatial_dimension(); - - if (sideTopoDim + 1 < elemSpatDim && sideTopoDim < elemTopoDim) { - sideOffset = parentTopo->number_faces(); - } - } - return sideOffset; - } - std::tuple get_touching_element_block_topology_from_side_block_by_tokenization(stk::io::OutputParams ¶ms, const stk::mesh::Part& part) { @@ -2261,6 +2274,18 @@ const stk::mesh::FieldBase *declare_stk_field_internal(stk::mesh::MetaData &meta } } + const stk::mesh::MetaData& meta = params.bulk_data().mesh_meta_data(); + std::vector touchingParts = meta.get_blocks_touching_surface(&part); + + if(touchingParts.size() == 1u) { + stk::topology stkTouchingTopology = touchingParts[0]->topology(); + std::string touchingTopoName = map_stk_topology_to_ioss(stkTouchingTopology); + const Ioss::ElementTopology *touchingTopo = Ioss::ElementTopology::factory(touchingTopoName, true); + if(touchingTopo != elementTopo) { + elementTopo = nullptr; + } + } + if (elementTopo != nullptr) { elementTopoName = elementTopo->name(); stkElementTopology = map_ioss_topology_to_stk(elementTopo, bulk.mesh_meta_data().spatial_dimension()); @@ -2320,6 +2345,13 @@ const stk::mesh::FieldBase *declare_stk_field_internal(stk::mesh::MetaData &meta Ioss::SideBlock *sideBlock = sset->get_side_block(name); if(sideBlock == nullptr) { + if("unknown" == ioTopo) { + // Special case of heterogenous sideset encapsulated in one side block + // There is no side topology so we cannot specify an element topology + // due to internal IOSS code that gives wrong SideBlock offset + elementTopoName = "unknown"; + } + sideBlock = new Ioss::SideBlock(sset->get_database(), name, ioTopo, elementTopoName, sideCount); sset->add(sideBlock); } diff --git a/packages/stk/stk_io/stk_io/IossBridge.hpp b/packages/stk/stk_io/stk_io/IossBridge.hpp index 1bfe68d84553..c28a91e58114 100644 --- a/packages/stk/stk_io/stk_io/IossBridge.hpp +++ b/packages/stk/stk_io/stk_io/IossBridge.hpp @@ -52,7 +52,6 @@ #include "Ioss_GroupingEntity.h" // for GroupingEntity #include "SidesetTranslator.hpp" // for fill_element_and_side_ids #include "stk_io/OutputParams.hpp" // for OutputParams -#include "stk_mesh/base/BulkData.hpp" // for BulkData #include "stk_mesh/base/FieldState.hpp" // for FieldState #include "stk_mesh/base/Part.hpp" // for Part #include "stk_topology/topology.hpp" // for topology @@ -642,6 +641,11 @@ const stk::mesh::Part* get_parent_element_block(const stk::mesh::BulkData &bulk, const Ioss::Region &ioRegion, const std::string& name); +int64_t get_side_offset(const Ioss::ElementTopology* sideTopo, + const Ioss::ElementTopology* parentTopo); + +int64_t get_side_offset(const Ioss::SideBlock* sb); + template void fill_data_for_side_block( OutputParams ¶ms, Ioss::GroupingEntity & io , @@ -661,7 +665,12 @@ void fill_data_for_side_block( OutputParams ¶ms, } // An offset required to translate Ioss's interpretation of shell ordinals - INT sideOrdOffset = (io.type() == Ioss::SIDEBLOCK) ? Ioss::Utils::get_side_offset(dynamic_cast(&io)) : 0; + INT sideOrdOffset = 0; + if(io.type() == Ioss::SIDEBLOCK) { + Ioss::SideBlock* sb = dynamic_cast(&io); + sideOrdOffset = get_side_offset(sb); + } + fill_element_and_side_ids(params, part, parentElementBlock, stk_elem_topology, sides, elem_side_ids, sideOrdOffset); } diff --git a/packages/stk/stk_io/stk_io/ProcessSetsOrBlocks.cpp b/packages/stk/stk_io/stk_io/ProcessSetsOrBlocks.cpp index b275bf697931..76939854e99f 100644 --- a/packages/stk/stk_io/stk_io/ProcessSetsOrBlocks.cpp +++ b/packages/stk/stk_io/stk_io/ProcessSetsOrBlocks.cpp @@ -208,6 +208,112 @@ void create_processed_face(stk::mesh::BulkData& bulk, } } +Ioss::ElementTopology* get_side_block_topology_from_entries(Ioss::DatabaseIO* db, Ioss::SideBlock* sb) +{ + Ioss::Region* region = db->get_region(); + if(nullptr == region) return nullptr; + + Ioss::Int64Vector elementSide; + if (db->int_byte_size_api() == 4) { + Ioss::IntVector es32; + sb->get_field_data("element_side", es32); + elementSide.resize(es32.size()); + std::copy(es32.begin(), es32.end(), elementSide.begin()); + } + else { + sb->get_field_data("element_side", elementSide); + } + + int heterogenousTopo = 0; + Ioss::ElementTopology* blockSideTopo = nullptr; + size_t number_sides = elementSide.size() / 2; + Ioss::ElementBlock *block = nullptr; + std::int64_t sideSetOffset = Ioss::Utils::get_side_offset(sb); + for (size_t iel = 0; iel < number_sides; iel++) { + int64_t elemId = elementSide[2 * iel]; + int64_t elemSide = elementSide[2 * iel + 1] + sideSetOffset - 1; + elemId = db->element_global_to_local(elemId); + if (block == nullptr || !block->contains(elemId)) { + block = region->get_element_block(elemId); + assert(block != nullptr); + } + + int64_t oneBasedElemSide = elemSide + 1; + Ioss::ElementTopology* sideTopo = block->topology()->boundary_type(oneBasedElemSide); + + if(nullptr != blockSideTopo && sideTopo != blockSideTopo) { + blockSideTopo = nullptr; + heterogenousTopo = 1; + break; + } + + blockSideTopo = sideTopo; + } + + int topoId = (blockSideTopo != nullptr) ? Ioss::ElementTopology::get_unique_id(blockSideTopo->name()) : 0; + + if (db->is_parallel()) { + std::vector topoVec{topoId, heterogenousTopo}; + db->util().global_array_minmax(topoVec, Ioss::ParallelUtils::DO_MAX); + topoId = topoVec[0]; + heterogenousTopo = topoVec[1]; + } + + blockSideTopo = heterogenousTopo ? nullptr : Ioss::ElementTopology::factory(topoId); + + return blockSideTopo; +} + +bool set_sideset_topology(const Ioss::SideSet *ss, stk::mesh::Part* part, + const Ioss::ElementTopology* sbTopo, stk::topology stkTopology, + bool printWarning = false) +{ + if (stkTopology == stk::topology::INVALID_TOPOLOGY) { + if (printWarning) { + std::ostringstream os; + os<<"stk_io WARNING: failed to obtain sensible topology for sideset: " << ss->name()<<", iossTopology: "<name()<<", stk-part: "<name()<<", rank: "<primary_entity_rank()<<", stk-topology: "<subsets()) { + if(subsetPart->topology() != stkTopology) { + return false; + } + } + + stk::mesh::set_topology(*part, stkTopology); + return true; + } + + return false; +} + +void set_sideset_topology(const Ioss::SideSet *ss, stk::mesh::Part *part, const stk::mesh::MetaData &meta) +{ + if(nullptr == ss) return; + if(nullptr == part) return; + if(ss->side_block_count() != 1) return; + + Ioss::DatabaseIO* db = ss->get_database(); + Ioss::Region* region = db->get_region(); + if(nullptr == region) return; + + Ioss::SideBlock* sb = ss->get_block(0); + + if(sb->name() != ss->name()) { + stk::topology stkTopology = map_ioss_topology_to_stk(sb->topology(), meta.spatial_dimension()); + if(set_sideset_topology(ss, part, sb->topology(), stkTopology, true)) return; + } + + Ioss::ElementTopology* blockSideTopo = get_side_block_topology_from_entries(db, sb); + + if(nullptr != blockSideTopo) { + stk::topology stkTopology = map_ioss_topology_to_stk(blockSideTopo, meta.spatial_dimension()); + set_sideset_topology(ss, part, blockSideTopo, stkTopology); + } +} + template void process_surface_entity(const Ioss::SideSet* sset, stk::mesh::BulkData & bulk, std::vector &sidesToMove, stk::io::StkMeshIoBroker::SideSetFaceCreationBehavior behavior) { diff --git a/packages/stk/stk_io/stk_io/SidesetTranslator.hpp b/packages/stk/stk_io/stk_io/SidesetTranslator.hpp index c76a6a9ae86d..0b4b3d64c007 100644 --- a/packages/stk/stk_io/stk_io/SidesetTranslator.hpp +++ b/packages/stk/stk_io/stk_io/SidesetTranslator.hpp @@ -161,11 +161,11 @@ void fill_element_and_side_ids_from_connectivity(stk::io::OutputParams ¶ms, size_t num_sides = stk::io::get_entities(params, *part, type, allSides, false); elem_side_ids.reserve(num_sides * 2); + std::vector side_elements; + std::vector side_nodes; for(size_t i = 0; i < num_sides; ++i) { stk::mesh::Entity side = allSides[i]; - std::vector side_elements; - std::vector side_nodes; fill_side_elements_and_nodes(bulk_data, side, side_elements, side_nodes); diff --git a/packages/stk/stk_io/stk_io/StkMeshIoBroker.hpp b/packages/stk/stk_io/stk_io/StkMeshIoBroker.hpp index 27678a95c270..918cbe1e796a 100644 --- a/packages/stk/stk_io/stk_io/StkMeshIoBroker.hpp +++ b/packages/stk/stk_io/stk_io/StkMeshIoBroker.hpp @@ -45,7 +45,6 @@ #include // for STKIORequire, FieldNameTo... #include // for MeshField, MeshField::CLO... #include // for OutputFile -#include // for BulkData #include // for Selector #include // for ParallelMachine #include // for Parameter, Type diff --git a/packages/stk/stk_mesh/Jamfile b/packages/stk/stk_mesh/Jamfile index 16b592a4b46d..91216f9cf3c0 100644 --- a/packages/stk/stk_mesh/Jamfile +++ b/packages/stk/stk_mesh/Jamfile @@ -52,11 +52,9 @@ project votd STK_SHOW_DEPRECATED_WARNINGS STK_HIDE_DEPRECATED_CODE SIERRA_MIGRATION - STK_16BIT_UPWARDCONN_INDEX_TYPE $(stk_mesh-root-inc) : usage-requirements SIERRA_MIGRATION - STK_16BIT_UPWARDCONN_INDEX_TYPE $(stk_mesh-root-inc) : build-dir $(stk_mesh-builddir) ; diff --git a/packages/stk/stk_mesh/stk_mesh/base/BulkData.cpp b/packages/stk/stk_mesh/stk_mesh/base/BulkData.cpp index 7fcbc2b88efd..10a3401d3b53 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/BulkData.cpp +++ b/packages/stk/stk_mesh/stk_mesh/base/BulkData.cpp @@ -152,7 +152,7 @@ void BulkData::check_if_entity_from_other_proc_exists_on_this_proc_and_update_in } } -void removeEntitiesNotSelected(stk::mesh::BulkData &mesh, stk::mesh::Selector selected, stk::mesh::EntityVector &entities) +void removeEntitiesNotSelected(stk::mesh::BulkData &mesh, const stk::mesh::Selector& selected, stk::mesh::EntityVector &entities) { if(selected != stk::mesh::Selector(mesh.mesh_meta_data().universal_part())) { @@ -1311,7 +1311,7 @@ uint64_t BulkData::get_max_allowed_id() const { const RelationVector& BulkData::aux_relations(Entity entity) const { - STK_ThrowAssert(m_add_fmwk_data); + STK_ThrowAssert(add_fmwk_data()); STK_ThrowAssert(entity.local_offset() > 0); if (m_fmwk_aux_relations[entity.local_offset()] == NULL) { @@ -1323,7 +1323,7 @@ BulkData::aux_relations(Entity entity) const RelationVector& BulkData::aux_relations(Entity entity) { - STK_ThrowAssert(m_add_fmwk_data); + STK_ThrowAssert(add_fmwk_data()); STK_ThrowAssert(entity.local_offset() > 0); if (m_fmwk_aux_relations[entity.local_offset()] == NULL) { @@ -4060,7 +4060,7 @@ void BulkData::determineEntitiesThatNeedGhosting(stk::mesh::Entity edge, void BulkData::find_upward_connected_entities_to_ghost_onto_other_processors(EntityProcVec& entitiesToGhostOntoOtherProcessors, EntityRank entity_rank, - stk::mesh::Selector selected, + const stk::mesh::Selector& selected, bool connectFacesToPreexistingGhosts) { if(entity_rank == stk::topology::NODE_RANK) { return; } @@ -4147,7 +4147,7 @@ void BulkData::internal_finish_modification_end(ModEndOptimizationFlag opt) notify_finished_mod_end(); } -bool BulkData::internal_modification_end_for_skin_mesh( EntityRank entity_rank, ModEndOptimizationFlag opt, stk::mesh::Selector selectedToSkin, +bool BulkData::internal_modification_end_for_skin_mesh( EntityRank entity_rank, ModEndOptimizationFlag opt, const stk::mesh::Selector& selectedToSkin, const Selector * only_consider_second_element_from_this_selector) { // The two states are MODIFIABLE and SYNCHRONiZED @@ -4188,7 +4188,7 @@ bool BulkData::internal_modification_end_for_skin_mesh( EntityRank entity_rank, return true ; } -void BulkData::resolve_incremental_ghosting_for_entity_creation_or_skin_mesh(EntityRank entity_rank, stk::mesh::Selector selectedToSkin, bool connectFacesToPreexistingGhosts) +void BulkData::resolve_incremental_ghosting_for_entity_creation_or_skin_mesh(EntityRank entity_rank, const stk::mesh::Selector& selectedToSkin, bool connectFacesToPreexistingGhosts) { EntityProcVec sendGhosts; find_upward_connected_entities_to_ghost_onto_other_processors(sendGhosts, entity_rank, selectedToSkin, connectFacesToPreexistingGhosts); @@ -5750,6 +5750,73 @@ void BulkData::mark_entities_as_deleted(stk::mesh::Bucket * bucket) } } +void +BulkData::internal_check_unpopulated_relations(Entity entity, EntityRank rank) const +{ +#if !defined(NDEBUG) && !defined(__HIP_DEVICE_COMPILE__) + if (m_check_invalid_rels) { + const MeshIndex &mesh_idx = mesh_index(entity); + const Bucket &b = *mesh_idx.bucket; + const unsigned bucket_ord = mesh_idx.bucket_ordinal; + STK_ThrowAssertMsg(count_valid_connectivity(entity, rank) == b.num_connectivity(bucket_ord, rank), + count_valid_connectivity(entity,rank) << " = count_valid_connectivity("<entity_rank(), + dst_mesh_idx.bucket->bucket_id(), + dst_mesh_idx.bucket_ordinal, + src_mesh_idx.bucket->bucket_id(), + src_mesh_idx.bucket_ordinal); +} + +bool +BulkData::relation_exist( const Entity entity, EntityRank subcell_rank, RelationIdentifier subcell_id ) +{ + bool found = false; + Entity const * rel_entity_it = bucket(entity).begin(bucket_ordinal(entity),subcell_rank); + const unsigned num_rel = bucket(entity).num_connectivity(bucket_ordinal(entity),subcell_rank); + ConnectivityOrdinal const * rel_ord_it = bucket(entity).begin_ordinals(bucket_ordinal(entity),subcell_rank); + + for (unsigned i=0 ; i < num_rel ; ++i) { + if (rel_ord_it[i] == static_cast(subcell_id) && is_valid(rel_entity_it[i])) { + found = true; + break; + } + } + + return found; +} + void BulkData::create_side_entities(const SideSet &sideSet, const stk::mesh::PartVector& parts) { if(has_face_adjacent_element_graph()) diff --git a/packages/stk/stk_mesh/stk_mesh/base/BulkData.hpp b/packages/stk/stk_mesh/stk_mesh/base/BulkData.hpp index 6d1da2382806..3d2323ecf578 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/BulkData.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/BulkData.hpp @@ -1,5 +1,7 @@ // Copyright 2002 - 2008, 2010, 2011 National Technology Engineering +// // Solutions of Sandia, LLC (NTESS). Under the terms of Contract +// // DE-NA0003525 with NTESS, the U.S. Government retains certain rights // in this software. // @@ -38,7 +40,6 @@ #include // for size_t #include // for uint16_t #include // for max -#include // for operator<<, basic_ostream, etc #include #include // for Entity, etc #include // for EntityCommDatabase @@ -131,7 +132,7 @@ stk::mesh::Entity clone_element_side(stk::mesh::BulkData &bulk, stk::mesh::ConnectivityOrdinal ord, const stk::mesh::PartVector &parts); void remove_ghosts_from_remote_procs(stk::mesh::BulkData& bulk, EntityVector& recvGhostsToRemove); -void destroy_elements(stk::mesh::BulkData &bulk, stk::mesh::EntityVector &elementsToDestroy, stk::mesh::Selector orphansToDelete); +void destroy_elements(stk::mesh::BulkData &bulk, stk::mesh::EntityVector &elementsToDestroy, const stk::mesh::Selector& orphansToDelete); namespace impl { stk::mesh::Entity connect_side_to_element(stk::mesh::BulkData& bulkData, stk::mesh::Entity element, @@ -315,7 +316,7 @@ class BulkData { * - Fields that exist on the dest that don't exist on the src will * be zeroed or initialized with the Field-specified initial-value. */ - inline void copy_entity_fields( Entity src, Entity dst); // CLEANUP: only used in usecases and perf tests - REMOVE from class + void copy_entity_fields( Entity src, Entity dst); //------------------------------------ /** \brief Query all buckets of a given entity rank @@ -540,7 +541,7 @@ class BulkData { void declare_relation( Entity e_from , const EntityVector& to_entities); - + //it's ugly to have 3 scratch-space vectors in the API, but for now //it is a big performance improvement. TODO: improve the internals to remove //the need for these vectors. @@ -930,7 +931,7 @@ class BulkData { bool internal_modification_end_for_skin_mesh( EntityRank entity_rank, ModEndOptimizationFlag opt, - stk::mesh::Selector selectedToSkin, + const stk::mesh::Selector& selectedToSkin, const stk::mesh::Selector * only_consider_second_element_from_this_selector); // Mod Mark bool inputs_ok_and_need_ghosting(Ghosting & ghosts , @@ -1241,7 +1242,7 @@ class BulkData { void internal_change_entity_key(EntityKey old_key, EntityKey new_key, Entity entity); // Mod Mark - void resolve_incremental_ghosting_for_entity_creation_or_skin_mesh(EntityRank entity_rank, stk::mesh::Selector selectedToSkin, bool connectFacesToPreexistingGhosts); + void resolve_incremental_ghosting_for_entity_creation_or_skin_mesh(EntityRank entity_rank, const stk::mesh::Selector& selectedToSkin, bool connectFacesToPreexistingGhosts); void internal_finish_modification_end(ModEndOptimizationFlag opt); // Mod Mark @@ -1428,7 +1429,7 @@ class BulkData { void require_entity_owner( const Entity entity, int owner) const ; - inline bool is_valid_connectivity(Entity entity, EntityRank rank) const; + bool is_valid_connectivity(Entity entity, EntityRank rank) const; friend class MeshBuilder; friend class ::stk::mesh::MetaData; @@ -1465,7 +1466,7 @@ class BulkData { const stk::mesh::PartVector &parts); friend void remove_ghosts_from_remote_procs(stk::mesh::BulkData& bulk, EntityVector& recvGhostsToRemove); - friend void destroy_elements(stk::mesh::BulkData &bulk, stk::mesh::EntityVector &elementsToDestroy, stk::mesh::Selector orphansToDelete); + friend void destroy_elements(stk::mesh::BulkData &bulk, stk::mesh::EntityVector &elementsToDestroy, const stk::mesh::Selector& orphansToDelete); friend stk::mesh::Entity impl::connect_side_to_element(stk::mesh::BulkData& bulkData, stk::mesh::Entity element, stk::mesh::EntityId side_global_id, stk::mesh::ConnectivityOrdinal side_ordinal, @@ -1484,7 +1485,7 @@ class BulkData { void find_upward_connected_entities_to_ghost_onto_other_processors(EntityProcVec& entitiesToGhostOntoOtherProcessors, EntityRank entity_rank, - stk::mesh::Selector selected, + const stk::mesh::Selector& selected, bool connectFacesToPreexistingGhosts); void reset_add_node_sharing() { m_add_node_sharing_called = false; } @@ -1990,56 +1991,6 @@ BulkData::in_send_ghost( Entity entity) const return false; } -inline void -BulkData::internal_check_unpopulated_relations(Entity entity, EntityRank rank) const -{ -#if !defined(NDEBUG) && !defined(__HIP_DEVICE_COMPILE__) - if (m_check_invalid_rels) { - const MeshIndex &mesh_idx = mesh_index(entity); - const Bucket &b = *mesh_idx.bucket; - const unsigned bucket_ord = mesh_idx.bucket_ordinal; - STK_ThrowAssertMsg(count_valid_connectivity(entity, rank) == b.num_connectivity(bucket_ord, rank), - count_valid_connectivity(entity,rank) << " = count_valid_connectivity("<entity_rank(), - dst_mesh_idx.bucket->bucket_id(), - dst_mesh_idx.bucket_ordinal, - src_mesh_idx.bucket->bucket_id(), - src_mesh_idx.bucket_ordinal); -} - -inline bool -BulkData::relation_exist( const Entity entity, EntityRank subcell_rank, RelationIdentifier subcell_id ) -{ - bool found = false; - Entity const * rel_entity_it = bucket(entity).begin(bucket_ordinal(entity),subcell_rank); - const unsigned num_rel = bucket(entity).num_connectivity(bucket_ordinal(entity),subcell_rank); - ConnectivityOrdinal const * rel_ord_it = bucket(entity).begin_ordinals(bucket_ordinal(entity),subcell_rank); - - for (unsigned i=0 ; i < num_rel ; ++i) { - if (rel_ord_it[i] == static_cast(subcell_id) && is_valid(rel_entity_it[i])) { - found = true; - break; - } - } - - return found; -} - #ifndef STK_HIDE_DEPRECATED_CODE // Delete after Sept 2024 STK_DEPRECATED inline VolatileFastSharedCommMapOneRank const& BulkData::volatile_fast_shared_comm_map(EntityRank rank) const @@ -2063,9 +2014,15 @@ BulkData::volatile_fast_shared_comm_map(EntityRank rank, int proc) const internal_update_ngp_fast_comm_maps(); } - const size_t dataBegin = m_ngpMeshHostData->hostVolatileFastSharedCommMapOffset[rank](proc); - const size_t dataEnd = m_ngpMeshHostData->hostVolatileFastSharedCommMapOffset[rank](proc+1); - return Kokkos::subview(m_ngpMeshHostData->hostVolatileFastSharedCommMap[rank], Kokkos::pair(dataBegin, dataEnd)); + if (parallel_size() > 1) + { + const size_t dataBegin = m_ngpMeshHostData->hostVolatileFastSharedCommMapOffset[rank](proc); + const size_t dataEnd = m_ngpMeshHostData->hostVolatileFastSharedCommMapOffset[rank](proc+1); + return Kokkos::subview(m_ngpMeshHostData->hostVolatileFastSharedCommMap[rank], Kokkos::pair(dataBegin, dataEnd)); + } else + { + return HostCommMapIndices("empty comm map indices", 0); + } } inline Part& @@ -2276,24 +2233,6 @@ BulkData::set_local_id(Entity entity, unsigned id) m_local_ids[entity.local_offset()] = id; } -inline void -BulkData::log_created_parallel_copy(Entity entity) -{ - if (state(entity) == Unchanged) { - set_state(entity, Modified); - } -} - -inline bool -BulkData::is_valid_connectivity(Entity entity, EntityRank rank) const -{ - if (!is_valid(entity)) return false; - if (bucket_ptr(entity) == NULL) return false; - internal_check_unpopulated_relations(entity, rank); - return true; -} - - namespace impl { inline NgpMeshBase * get_ngp_mesh(const BulkData & bulk) { return bulk.get_ngp_mesh(); diff --git a/packages/stk/stk_mesh/stk_mesh/base/DataTraits.hpp b/packages/stk/stk_mesh/stk_mesh/base/DataTraits.hpp index 96b505bdd6c8..bb869c133be6 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/DataTraits.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/DataTraits.hpp @@ -6,15 +6,15 @@ // Redistribution and use in source and binary forms, with or without // modification, are permitted provided that the following conditions are // met: -// +// // * Redistributions of source code must retain the above copyright // notice, this list of conditions and the following disclaimer. -// +// // * Redistributions in binary form must reproduce the above // copyright notice, this list of conditions and the following // disclaimer in the documentation and/or other materials provided // with the distribution. -// +// // * Neither the name of NTESS nor the names of its contributors // may be used to endorse or promote products derived from this // software without specific prior written permission. @@ -30,7 +30,7 @@ // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -// +// #ifndef STK_MESH_DATA_TRAITS_HPP #define STK_MESH_DATA_TRAITS_HPP @@ -195,6 +195,16 @@ class DataTraits { DataTraits & operator = ( const DataTraits & ); }; +inline bool operator==(const DataTraits& lhs, const DataTraits& rhs) +{ + return lhs.type_info == rhs.type_info; +} + +inline bool operator!=(const DataTraits& lhs, const DataTraits& rhs) +{ + return lhs.type_info != rhs.type_info; +} + } // namespace mesh } // namespace stk diff --git a/packages/stk/stk_mesh/stk_mesh/base/DestroyElements.cpp b/packages/stk/stk_mesh/stk_mesh/base/DestroyElements.cpp index df4522b5edd4..79547434a93c 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/DestroyElements.cpp +++ b/packages/stk/stk_mesh/stk_mesh/base/DestroyElements.cpp @@ -91,7 +91,7 @@ void destroy_elements(stk::mesh::BulkData &bulk, stk::mesh::EntityVector &elemen destroy_elements(bulk, elementsToDestroy, orphansToDelete); } -void destroy_elements(stk::mesh::BulkData &bulk, stk::mesh::EntityVector &elementsToDestroy, stk::mesh::Selector orphansToDelete) +void destroy_elements(stk::mesh::BulkData &bulk, stk::mesh::EntityVector &elementsToDestroy, const stk::mesh::Selector& orphansToDelete) { bulk.modification_begin(); bulk.m_bucket_repository.set_remove_mode_tracking(); @@ -119,7 +119,7 @@ void get_all_related_entities(BulkData& bulk, EntityVector& elements, const Sele stk::util::sort_and_unique(relatedEntities, stk::mesh::EntityLess(bulk)); } -void destroy_elements_no_mod_cycle(stk::mesh::BulkData &bulk, stk::mesh::EntityVector &elementsToDestroy, stk::mesh::Selector orphansToDelete) +void destroy_elements_no_mod_cycle(stk::mesh::BulkData &bulk, stk::mesh::EntityVector &elementsToDestroy, const stk::mesh::Selector& orphansToDelete) { for(stk::mesh::Entity element : elementsToDestroy) { if(!bulk.is_valid(element)) diff --git a/packages/stk/stk_mesh/stk_mesh/base/DestroyElements.hpp b/packages/stk/stk_mesh/stk_mesh/base/DestroyElements.hpp index c1f515cf4fa1..af53cca3ae0d 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/DestroyElements.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/DestroyElements.hpp @@ -36,15 +36,16 @@ #define STK_MESH_DESTROY_ELEMENTS_HPP #include -#include -#include namespace stk { namespace mesh { +class BulkData; +class Selector; + void destroy_elements(stk::mesh::BulkData &bulk, stk::mesh::EntityVector &elementsToDestroy); -void destroy_elements(stk::mesh::BulkData &bulk, stk::mesh::EntityVector &elementsToDestroy, stk::mesh::Selector orphansToDelete); -void destroy_elements_no_mod_cycle(stk::mesh::BulkData &bulk, stk::mesh::EntityVector &elementsToDestroy, stk::mesh::Selector orphansToDelete); +void destroy_elements(stk::mesh::BulkData &bulk, stk::mesh::EntityVector &elementsToDestroy, const stk::mesh::Selector& orphansToDelete); +void destroy_elements_no_mod_cycle(stk::mesh::BulkData &bulk, stk::mesh::EntityVector &elementsToDestroy, const stk::mesh::Selector& orphansToDelete); void destroy_upward_connected_aura_entities(stk::mesh::BulkData &bulk, stk::mesh::Entity connectedEntity, stk::mesh::EntityRank conRank); }} // namespace stk::mesh diff --git a/packages/stk/stk_mesh/stk_mesh/base/Field.hpp b/packages/stk/stk_mesh/stk_mesh/base/Field.hpp index e3cf3a05eb3f..40cd9d3940a6 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/Field.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/Field.hpp @@ -73,8 +73,8 @@ namespace mesh { * using CoordFieldType stk::mesh::Field CoordFieldType; * CoordFieldType& coord_field = meta.declare_field("", ); * - * - Field restrictions: Defines the set of entities that have a field and the dimensions of the - * field (maximum of 7 dimensions). FieldRestrictions are applied to Parts and entities pick + * - Field restrictions: Defines the set of entities that have a field and the memory layout of + * the data (e.g. a vector-3). FieldRestrictions are applied to Parts and entities pick * up the field by being placed in the associated part. Also, FieldRestrictions allow the same * Field to be applied to different ranks of entities. E.g., you could apply a velocity Field * to both Faces and Edges through two different FieldRestrictions for each of the different ranks. @@ -136,44 +136,30 @@ namespace mesh { * - stk_mesh/base/FieldParallel - Free functions for some parallel operations * - stk_mesh/base/FieldRestriction - Defines FieldRestriction class * - stk_mesh/base/FieldState - Defines enums to refer to field states - * - stk_mesh/base/FieldTraits - Defines API for querying field-types - * - stk_mesh/base/CoordinateSystems - Defines common ArrayDimTags (used for defining Field types) - * - stk_mesh/base/TopologyDimensions - Contains useful Field types / ArrayDimTags for setting up field relations * - Type-related in ArrayDim, DataTraits * * - TODO Describe relationship with Buckets */ -// Implementation Details: -// The template arguments below describe the field type. Scalar is the scalar -// type of data contained by the field. The TagN describe each dimension of the -// Field, these are expected to be ArrayDimTags. Unused dimensions can be ignored. -template< typename Scalar , class Tag1 , class Tag2 , class Tag3 , class Tag4 , - class Tag5 , class Tag6 , class Tag7 > + +template class Field : public FieldBase { public: using value_type = Scalar; - Field( - MetaData * arg_mesh_meta_data , - stk::topology::rank_t arg_entity_rank , - unsigned arg_ordinal , - const std::string & arg_name , - const DataTraits & arg_traits , - unsigned arg_rank, - const shards::ArrayDimTag * const * arg_dim_tags, - unsigned arg_number_of_states , - FieldState arg_this_state - ) + Field(MetaData * arg_mesh_meta_data, + stk::topology::rank_t arg_entity_rank, + unsigned arg_ordinal, + const std::string & arg_name, + const DataTraits & arg_traits, + unsigned arg_number_of_states, + FieldState arg_this_state) : FieldBase(arg_mesh_meta_data, - arg_entity_rank, - arg_ordinal, - arg_name, - arg_traits, - arg_rank, - arg_dim_tags, - arg_number_of_states, - arg_this_state - ) + arg_entity_rank, + arg_ordinal, + arg_name, + arg_traits, + arg_number_of_states, + arg_this_state) {} /** \brief Query this field for a given field state. */ @@ -185,7 +171,7 @@ class Field : public FieldBase { #endif } - virtual ~Field(){} + virtual ~Field() = default; virtual std::ostream& print_data(std::ostream& out, void* data, unsigned size_per_entity) const override { @@ -260,8 +246,6 @@ class Field : public FieldBase { fieldRepo.get_fields().size(), fieldNames[i], data_traits(), - m_field_rank, - m_dim_tags, number_of_states(), static_cast(i)); fieldRepo.add_field(f[i]); diff --git a/packages/stk/stk_mesh/stk_mesh/base/FieldBLAS.hpp b/packages/stk/stk_mesh/stk_mesh/base/FieldBLAS.hpp index c49ee0838715..8c8ec50cf796 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/FieldBLAS.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/FieldBLAS.hpp @@ -48,7 +48,6 @@ #include #include -#include #include #if defined(_OPENMP) && !defined(__INTEL_COMPILER) @@ -618,13 +617,15 @@ void field_copy(const FieldBase& xField, const FieldBase& yField) field_copy(xField,yField,selector); } -template +template inline -Scalar field_dot(const Field & xField, const Field & yField, const Selector& selector, const MPI_Comm comm) +Scalar field_dot(const Field & xField, const Field & yField, + const Selector& selector, const MPI_Comm comm) { STK_ThrowRequire(is_compatible(xField, yField)); - BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), selector & xField.mesh_meta_data().locally_owned_part()); + BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), + selector & xField.mesh_meta_data().locally_owned_part()); Scalar local_result(0.0); @@ -646,13 +647,16 @@ Scalar field_dot(const Field & xField, const Field< return glob_result; } -template +template inline -std::complex field_dot(const Field,T1,T2,T3,T4,T5,T6,T7>& xField, const Field,T1,T2,T3,T4,T5,T6,T7>& yField, const Selector& selector, const MPI_Comm comm) +std::complex field_dot(const Field>& xField, + const Field>& yField, + const Selector& selector, const MPI_Comm comm) { STK_ThrowRequire(is_compatible(xField, yField)); - BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), selector & xField.mesh_meta_data().locally_owned_part()); + BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), + selector & xField.mesh_meta_data().locally_owned_part()); Scalar local_result_r (0.0); Scalar local_result_i (0.0); @@ -680,17 +684,17 @@ std::complex field_dot(const Field,T1,T2,T3,T4,T5,T return std::complex (glob_result_ri[0],glob_result_ri[1]); } -template +template inline -Scalar field_dot(const Field & xField, const Field & yField, const Selector& selector) +Scalar field_dot(const Field & xField, const Field & yField, const Selector& selector) { const MPI_Comm comm = xField.get_mesh().parallel(); return field_dot(xField,yField,selector,comm); } -template +template inline -Scalar field_dot(const Field & xField, const Field & yField) +Scalar field_dot(const Field & xField, const Field & yField) { const Selector selector = selectField(xField) & selectField(yField); return field_dot(xField,yField,selector); @@ -935,11 +939,12 @@ void field_swap(const FieldBase& xField, const FieldBase& yField) field_swap(xField,yField,selector); } -template +template inline -Scalar field_nrm2(const Field & xField, const Selector& selector, const MPI_Comm comm) +Scalar field_nrm2(const Field & xField, const Selector& selector, const MPI_Comm comm) { - BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), selector & xField.mesh_meta_data().locally_owned_part()); + BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), + selector & xField.mesh_meta_data().locally_owned_part()); Scalar local_result(0.0); @@ -958,9 +963,10 @@ Scalar field_nrm2(const Field & xField, const Selec return std::sqrt(glob_result); } -template +template inline -std::complex field_nrm2(const Field< std::complex,T1,T2,T3,T4,T5,T6,T7> & xField, const Selector& selector, const MPI_Comm comm) +std::complex field_nrm2(const Field> & xField, const Selector& selector, + const MPI_Comm comm) { BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), selector & xField.mesh_meta_data().locally_owned_part()); @@ -981,17 +987,17 @@ std::complex field_nrm2(const Field< std::complex,T1,T2,T3,T4,T5 return std::complex(std::sqrt(glob_result),0.0); } -template +template inline -Scalar field_nrm2(const Field & xField, const Selector& selector) +Scalar field_nrm2(const Field & xField, const Selector& selector) { const MPI_Comm comm = xField.get_mesh().parallel(); return field_nrm2(xField,selector,comm); } -template +template inline -Scalar field_nrm2(const Field & xField) +Scalar field_nrm2(const Field & xField) { const Selector selector = selectField(xField); return field_nrm2(xField,selector); @@ -1063,11 +1069,12 @@ void field_nrm2(Scalar& result, const FieldBase& xField) field_nrm2(result,xField,selector); } -template +template inline -Scalar field_asum(const Field & xField, const Selector& selector, const MPI_Comm comm) +Scalar field_asum(const Field & xField, const Selector& selector, const MPI_Comm comm) { - BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), selector & xField.mesh_meta_data().locally_owned_part()); + BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), + selector & xField.mesh_meta_data().locally_owned_part()); Scalar local_result(0.0); @@ -1086,11 +1093,13 @@ Scalar field_asum(const Field & xField, const Selec return glob_result; } -template +template inline -std::complex field_asum(const Field,T1,T2,T3,T4,T5,T6,T7> & xField, const Selector& selector, const MPI_Comm comm) +std::complex field_asum(const Field> & xField, const Selector& selector, + const MPI_Comm comm) { - BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), selector & xField.mesh_meta_data().locally_owned_part()); + BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), + selector & xField.mesh_meta_data().locally_owned_part()); Scalar local_result(0.0); @@ -1109,17 +1118,17 @@ std::complex field_asum(const Field,T1,T2,T3,T4,T5, return std::complex(glob_result,0.0); } -template +template inline -Scalar field_asum(const Field & xField, const Selector& selector) +Scalar field_asum(const Field & xField, const Selector& selector) { const MPI_Comm comm = xField.get_mesh().parallel(); return field_asum(xField,selector,comm); } -template +template inline -Scalar field_asum(const Field & xField) +Scalar field_asum(const Field & xField) { const Selector selector = selectField(xField); return field_asum(xField,selector); @@ -1190,9 +1199,9 @@ void field_asum(Scalar& result, const FieldBase& xField) field_asum(result,xField,selector); } -template +template inline -Scalar field_amax(const Field & xField, const Selector& selector, const MPI_Comm comm) +Scalar field_amax(const Field & xField, const Selector& selector, const MPI_Comm comm) { BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), selector & xField.mesh_meta_data().locally_owned_part()); @@ -1218,11 +1227,13 @@ Scalar field_amax(const Field & xField, const Selec return global_amax; } -template +template inline -std::complex field_amax(const Field,T1,T2,T3,T4,T5,T6,T7> & xField, const Selector& selector, const MPI_Comm comm) +std::complex field_amax(const Field> & xField, const Selector& selector, + const MPI_Comm comm) { - BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), selector & xField.mesh_meta_data().locally_owned_part()); + BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), + selector & xField.mesh_meta_data().locally_owned_part()); Scalar priv_tmp; Scalar local_amax(0.0); @@ -1246,17 +1257,17 @@ std::complex field_amax(const Field,T1,T2,T3,T4,T5, return std::complex(glob_amax,0.0); } -template +template inline -Scalar field_amax(const Field & xField, const Selector& selector) +Scalar field_amax(const Field & xField, const Selector& selector) { const MPI_Comm comm = xField.get_mesh().parallel(); return field_amax(xField,selector,comm); } -template +template inline -Scalar field_amax(const Field & xField) +Scalar field_amax(const Field & xField) { const Selector selector = selectField(xField); return field_amax(xField,selector); @@ -1396,9 +1407,9 @@ Entity field_eamax(const FieldBase& xField) return field_eamax(xField,selector); } -template +template inline -Entity field_eamax(const Field & xField) +Entity field_eamax(const Field & xField) { const Selector selector = selectField(xField); return field_eamax(xField,selector); @@ -1480,11 +1491,12 @@ void field_amax(Scalar& result, const FieldBase& xField) field_amax(result,xField,selector); } -template +template inline -Entity field_eamin(const Field,T1,T2,T3,T4,T5,T6,T7> & xField, const Selector& selector) +Entity field_eamin(const Field> & xField, const Selector& selector) { - BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), selector & xField.mesh_meta_data().locally_owned_part()); + BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), + selector & xField.mesh_meta_data().locally_owned_part()); int priv_iamin; Scalar priv_amin; @@ -1535,11 +1547,12 @@ Entity field_eamin(const Field,T1,T2,T3,T4,T5,T6,T7> & xFie return glob_result; } -template +template inline -Entity field_eamin(const Field & xField, const Selector& selector) +Entity field_eamin(const Field & xField, const Selector& selector) { - BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), selector & xField.mesh_meta_data().locally_owned_part()); + BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), + selector & xField.mesh_meta_data().locally_owned_part()); int priv_iamin; Scalar priv_amin; @@ -1590,19 +1603,21 @@ Entity field_eamin(const Field & xField, const Sele return glob_result; } -template +template inline -Entity field_eamin(const Field & xField) +Entity field_eamin(const Field & xField) { const Selector selector = selectField(xField); return field_eamin(xField,selector); } -template +template inline -std::complex field_amin(const Field,T1,T2,T3,T4,T5,T6,T7> & xField, const Selector& selector, const MPI_Comm comm) +std::complex field_amin(const Field> & xField, const Selector& selector, + const MPI_Comm comm) { - BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), selector & xField.mesh_meta_data().locally_owned_part()); + BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), + selector & xField.mesh_meta_data().locally_owned_part()); Scalar priv_tmp; Scalar local_amin = std::numeric_limits::max(); @@ -1626,9 +1641,9 @@ std::complex field_amin(const Field,T1,T2,T3,T4,T5, return std::sqrt(glob_amin); } -template +template inline -Scalar field_amin(const Field & xField, const Selector& selector, const MPI_Comm comm) +Scalar field_amin(const Field & xField, const Selector& selector, const MPI_Comm comm) { BucketVector const& buckets = xField.get_mesh().get_buckets(xField.entity_rank(), selector & xField.mesh_meta_data().locally_owned_part()); @@ -1654,17 +1669,17 @@ Scalar field_amin(const Field & xField, const Selec return glob_amin; } -template +template inline -Scalar field_amin(const Field & xField, const Selector& selector) +Scalar field_amin(const Field & xField, const Selector& selector) { const MPI_Comm comm = xField.get_mesh().parallel(); return field_amin(xField,selector,comm); } -template +template inline -Scalar field_amin(const Field & xField) +Scalar field_amin(const Field & xField) { const Selector selector = selectField(xField); return field_amin(xField,selector); diff --git a/packages/stk/stk_mesh/stk_mesh/base/FieldBase.cpp b/packages/stk/stk_mesh/stk_mesh/base/FieldBase.cpp index 11c5abc3e604..d96dcada0712 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/FieldBase.cpp +++ b/packages/stk/stk_mesh/stk_mesh/base/FieldBase.cpp @@ -150,19 +150,18 @@ void FieldBase::set_initial_value(const void* new_initial_value, unsigned num_sc data_traits().copy(init_val, new_initial_value, num_scalars); } -void FieldBase::insert_restriction( - const char * arg_method , - const Part & arg_part , - const unsigned arg_num_scalars_per_entity , - const unsigned arg_first_dimension , - const void* arg_init_value ) +void FieldBase::insert_restriction(const char * arg_method, + const Part & arg_part, + const unsigned arg_num_scalars_per_entity, + const unsigned arg_first_dimension, + const void* arg_init_value) { FieldRestriction tmp( arg_part ); tmp.set_num_scalars_per_entity(arg_num_scalars_per_entity); tmp.set_dimension(arg_first_dimension); - if (arg_init_value != NULL) { + if (arg_init_value != nullptr) { //insert_restriction can be called multiple times for the same field, giving //the field different lengths on different mesh-parts. //We will only store one initial-value array, we need to store the one with @@ -182,7 +181,7 @@ void FieldBase::insert_restriction( size_t nbytes = sizeof_scalar * num_scalars; size_t old_nbytes = 0; - if (get_initial_value() != NULL) { + if (get_initial_value() != nullptr) { old_nbytes = get_initial_value_num_bytes(); } if (nbytes > old_nbytes) { @@ -269,19 +268,18 @@ void FieldBase::insert_restriction( } } -void FieldBase::insert_restriction( - const char * arg_method , - const Selector & arg_selector , - const unsigned arg_num_scalars_per_entity , - const unsigned arg_first_dimension , - const void* arg_init_value ) +void FieldBase::insert_restriction(const char * arg_method, + const Selector & arg_selector, + const unsigned arg_num_scalars_per_entity, + const unsigned arg_first_dimension, + const void* arg_init_value) { FieldRestriction tmp( arg_selector ); tmp.set_num_scalars_per_entity(arg_num_scalars_per_entity); tmp.set_dimension(arg_first_dimension); - if (arg_init_value != NULL) { + if (arg_init_value != nullptr) { //insert_restriction can be called multiple times for the same field, giving //the field different lengths on different mesh-parts. //We will only store one initial-value array, we need to store the one with @@ -301,7 +299,7 @@ void FieldBase::insert_restriction( size_t nbytes = sizeof_scalar * num_scalars; size_t old_nbytes = 0; - if (get_initial_value() != NULL) { + if (get_initial_value() != nullptr) { old_nbytes = get_initial_value_num_bytes(); } if (nbytes > old_nbytes) { @@ -440,7 +438,7 @@ void FieldBase::verify_and_clean_restrictions(const Part& superset, const Part& void FieldBase::set_mesh(stk::mesh::BulkData* bulk) { - if (m_mesh == NULL || bulk == NULL) { + if (m_mesh == nullptr || bulk == nullptr) { m_mesh = bulk; } else { diff --git a/packages/stk/stk_mesh/stk_mesh/base/FieldBase.hpp b/packages/stk/stk_mesh/stk_mesh/base/FieldBase.hpp index e4c12e29d702..99c56712f671 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/FieldBase.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/FieldBase.hpp @@ -36,13 +36,13 @@ #define stk_mesh_base_FieldBase_hpp #include -#include // for NULL +#include // for size_t #include // for ostream #include // for Bucket #include // for BulkData #include // for FieldRestriction, etc #include // for FieldState -#include // for FieldTraits, EntityRank, etc +#include // for EntityRank, etc #include #include #include @@ -55,8 +55,6 @@ #include #include -namespace shards { class ArrayDimTag; } - namespace stk::mesh { @@ -321,13 +319,13 @@ class FieldBase const Part & arg_part , const unsigned arg_num_scalars_per_entity , const unsigned arg_first_dimension , - const void* arg_init_value = NULL); + const void* arg_init_value = nullptr); void insert_restriction(const char * arg_method , const Selector & arg_selector , const unsigned arg_num_scalars_per_entity , const unsigned arg_first_dimension , - const void* arg_init_value = NULL); + const void* arg_init_value = nullptr); void verify_and_clean_restrictions( const Part& superset, const Part& subset ); @@ -344,7 +342,7 @@ class FieldBase template class NgpDebugger> friend class HostField; template class NgpDebugger> friend class DeviceField; - template friend class Field; + template friend class Field; protected: FieldBase(MetaData * arg_mesh_meta_data, @@ -352,15 +350,12 @@ class FieldBase unsigned arg_ordinal, const std::string & arg_name, const DataTraits & arg_traits, - unsigned arg_rank, - const shards::ArrayDimTag * const * arg_dim_tags, unsigned arg_number_of_states, FieldState arg_this_state) : m_mesh(nullptr), m_entity_rank(arg_entity_rank), m_name(arg_name), m_num_states(arg_number_of_states), - m_field_rank( arg_rank ), m_restrictions(), m_initial_value(nullptr), m_initial_value_num_bytes(0), @@ -376,13 +371,7 @@ class FieldBase m_execSpace(Kokkos::DefaultExecutionSpace()) { FieldBase * const pzero = nullptr ; - const shards::ArrayDimTag * const dzero = nullptr ; Copy( m_field_states , pzero ); - Copy( m_dim_tags , dzero ); - - for ( unsigned i = 0 ; i < arg_rank ; ++i ) { - m_dim_tags[i] = arg_dim_tags[i]; - } } private: @@ -392,17 +381,15 @@ class FieldBase EntityRank m_entity_rank; const std::string m_name; const unsigned m_num_states; - unsigned m_field_rank; FieldBase * m_field_states[ MaximumFieldStates ]; - const shards::ArrayDimTag * m_dim_tags[ MaximumFieldDimension ]; - FieldRestrictionVector m_restrictions; + FieldRestrictionVector m_restrictions; void* m_initial_value; unsigned m_initial_value_num_bytes; CSet m_attribute; - const DataTraits & m_data_traits ; - MetaData * const m_meta_data ; - const unsigned m_ordinal ; - const FieldState m_this_state ; + const DataTraits & m_data_traits; + MetaData * const m_meta_data; + const unsigned m_ordinal; + const FieldState m_this_state; mutable NgpFieldBase * m_ngpField; mutable size_t m_numSyncsToHost; mutable size_t m_numSyncsToDevice; @@ -670,7 +657,7 @@ field_data(const FieldType & f, const unsigned bucket_id, unsigned bucket_ord, c STK_ThrowAssertMsg(f.get_meta_data_for_field()[bucket_id].m_bytesPerEntity >= knownSize, "field name= " << f.name() << "knownSize= " << knownSize << ", m_bytesPerEntity= " << f.get_meta_data_for_field()[bucket_id].m_bytesPerEntity); - STK_ThrowAssert(f.get_meta_data_for_field()[bucket_id].m_data != NULL); + STK_ThrowAssert(f.get_meta_data_for_field()[bucket_id].m_data != nullptr); debug_stale_access_entity_check(static_cast(f), bucket_id, bucket_ord, fileName, lineNumber); diff --git a/packages/stk/stk_mesh/stk_mesh/base/ForEachEntity.hpp b/packages/stk/stk_mesh/stk_mesh/base/ForEachEntity.hpp index 96a4b492225f..847c33a595b3 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/ForEachEntity.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/ForEachEntity.hpp @@ -37,13 +37,14 @@ #include #include #include -#include #include #include namespace stk { namespace mesh { +class BulkData; + template void for_each_entity_run(const BulkData &mesh, stk::topology::rank_t rank, diff --git a/packages/stk/stk_mesh/stk_mesh/base/Mesh.dox b/packages/stk/stk_mesh/stk_mesh/base/Mesh.dox index ad9d4de9b0b6..cf44f013add4 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/Mesh.dox +++ b/packages/stk/stk_mesh/stk_mesh/base/Mesh.dox @@ -46,14 +46,6 @@ * and derived field relationships. */ -/** \ingroup stk_mesh_module - * \defgroup stk_mesh_field_dimension_tags Meta Data Field Dimension Tags - * \brief Mesh specific - * \ref shards::ArrayDimTag "Array dimension tags" for defining - * \ref shards::Array "multi-dimensional array" - * \ref stk_mesh_field_data "field data". - */ - //---------------------------------------------------------------------- /** \ingroup stk_mesh_module diff --git a/packages/stk/stk_mesh/stk_mesh/base/MetaData.hpp b/packages/stk/stk_mesh/stk_mesh/base/MetaData.hpp index b33b572b66a9..0585b10fbd52 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/MetaData.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/MetaData.hpp @@ -67,7 +67,6 @@ #include "stk_util/diag/StringUtil.hpp" #include -namespace shards { class ArrayDimTag; } namespace shards { class CellTopologyManagedData; } namespace stk { namespace mesh { class BulkData; } } namespace stk { namespace mesh { class MetaData; } } @@ -77,8 +76,8 @@ namespace mesh { template struct is_field : std::false_type {}; -template -struct is_field> : std::true_type {}; +template +struct is_field> : std::true_type {}; template constexpr bool is_field_v = is_field::value; @@ -91,16 +90,6 @@ template constexpr bool is_field_base_v = is_field_base::value; -template -struct is_simple_field : std::false_type {}; - -template -struct is_simple_field> : std::true_type {}; - -template -constexpr bool is_simple_field_v = is_simple_field::value; - - /** \addtogroup stk_mesh_module * \{ */ @@ -784,16 +773,14 @@ Field * MetaData::get_field(stk::mesh::EntityRank arg_entity_rank, const char * fileName, int lineNumber) const { - static_assert(not is_field_v, "You must use a datatype as the template parameter to MetaData::get_field()," - "and not the Field itself"); + static_assert(not is_field_v && not is_field_base_v, + "You must use a datatype as the template parameter to MetaData::get_field(), " + "and not the Field itself"); const DataTraits & dt = data_traits(); const DataTraits & dt_void = data_traits(); - const int fieldRank = 0; - - std::array tags {nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr}; - FieldBase * const field = m_field_repo.get_field(arg_entity_rank, name, dt, fieldRank, tags.data(), 0); + FieldBase * const field = m_field_repo.get_field(arg_entity_rank, name, dt, 0); STK_ThrowRequireMsg(field == nullptr || field->data_traits().type_info == dt.type_info || @@ -812,13 +799,11 @@ MetaData::declare_field(stk::topology::rank_t arg_entity_rank, const char * fileName, int lineNumber) { - static_assert(not is_field_v, "You must use a datatype as the template parameter to MetaData::declare_field()," - "and not the Field itself"); + static_assert(not is_field_v && not is_field_base_v, + "You must use a datatype as the template parameter to MetaData::declare_field(), " + "and not the Field itself"); const DataTraits & traits = data_traits(); - const int fieldRank = 0; - - std::array tags {nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr}; const char** reservedStateSuffix = reserved_state_suffix(); @@ -841,8 +826,7 @@ MetaData::declare_field(stk::topology::rank_t arg_entity_rank, Field * f[MaximumFieldStates] = {nullptr}; - FieldBase* rawField = m_field_repo.get_field(arg_entity_rank, name, - traits, fieldRank, tags.data(), number_of_states); + FieldBase* rawField = m_field_repo.get_field(arg_entity_rank, name, traits, number_of_states); f[0] = dynamic_cast*>(rawField); @@ -881,8 +865,6 @@ MetaData::declare_field(stk::topology::rank_t arg_entity_rank, m_field_repo.get_fields().size(), field_names[i], traits, - fieldRank, - tags.data(), number_of_states, static_cast(i)); @@ -905,10 +887,6 @@ field_type & put_field_on_mesh(field_type & field, const Part & part, const typename field_type::value_type* init_value) { - static_assert(is_simple_field_v || is_field_base_v, - "You must only call put_field_on_mesh() with a simple field argument (i.e. without template parameters" - " beyond the datatype"); - MetaData & meta = MetaData::get(field); unsigned numScalarsPerEntity = 1; @@ -924,10 +902,6 @@ field_type & put_field_on_mesh(field_type & field, const Selector & selector, const typename field_type::value_type* init_value) { - static_assert(is_simple_field_v || is_field_base_v, - "You must only call put_field_on_mesh() with a simple field argument (i.e. without template parameters" - " beyond the datatype"); - MetaData & meta = MetaData::get(field); unsigned numScalarsPerEntity = 1; @@ -1102,8 +1076,9 @@ Field * get_field_by_name(const std::string & name, const char * fileName = HOST_DEBUG_FILE_NAME, int lineNumber = HOST_DEBUG_LINE_NUMBER) { - static_assert(not is_field_v, "You must use a datatype as the template parameter to get_field_by_name()," - "and not the Field itself"); + static_assert(not is_field_v && not is_field_base_v, + "You must use a datatype as the template parameter to get_field_by_name(), " + "and not the Field itself"); Field* field = nullptr; unsigned num_nonnull_fields = 0; diff --git a/packages/stk/stk_mesh/stk_mesh/base/NgpFieldBLAS.hpp b/packages/stk/stk_mesh/stk_mesh/base/NgpFieldBLAS.hpp index 84d8132f08bf..6ff9aaa6fdf9 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/NgpFieldBLAS.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/NgpFieldBLAS.hpp @@ -6,15 +6,15 @@ // Redistribution and use in source and binary forms, with or without // modification, are permitted provided that the following conditions are // met: -// +// // * Redistributions of source code must retain the above copyright // notice, this list of conditions and the following disclaimer. -// +// // * Redistributions in binary form must reproduce the above // copyright notice, this list of conditions and the following // disclaimer in the documentation and/or other materials provided // with the distribution. -// +// // * Neither the name of NTESS nor the names of its contributors // may be used to endorse or promote products derived from this // software without specific prior written permission. @@ -30,7 +30,7 @@ // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -// +// #ifndef STK_MESH_BASE_NGPFIELDBLAS_HPP #define STK_MESH_BASE_NGPFIELDBLAS_HPP @@ -118,60 +118,144 @@ inline void field_axpby(const stk::mesh::BulkData& mesh, const stk::mesh::FieldBase & yField, const stk::mesh::Selector & selector, const EXEC_SPACE& execSpace, - bool IsDeviceExecSpaceUserOverride = (!std::is_same_v)) + bool isDeviceExecSpaceUserOverride = (!std::is_same_v)) { // y = a*x + b*y - - if constexpr (ngp_field_blas::impl::operate_on_ngp_mesh()) { - ngp_field_blas::impl::apply_functor_on_field( - mesh, yField, xField, yField, alpha, beta, selector); - } - else { - xField.sync_to_host(); - yField.sync_to_host(); - stk::mesh::field_axpby(alpha, xField, beta, yField, selector); - } - - yField.clear_sync_state(); - if (ngp_field_blas::impl::mark_modified_on_device(execSpace, IsDeviceExecSpaceUserOverride)) { - yField.modify_on_device(); - } - else { - yField.modify_on_host(); - } + field_axpbyz(mesh, alpha, xField, beta, yField, yField, selector, execSpace, isDeviceExecSpaceUserOverride); } template -inline void field_axpbyz(const stk::mesh::BulkData& mesh, +inline void field_axpby(const stk::mesh::BulkData& mesh, const DataType alpha, const stk::mesh::FieldBase & xField, const DataType beta, const stk::mesh::FieldBase & yField, + const EXEC_SPACE& execSpace, + bool isDeviceExecSpaceUserOverride = (!std::is_same_v)) +{ + // y = a*x + b*y + field_axpbyz(mesh, alpha, xField, beta, yField, yField, execSpace, isDeviceExecSpaceUserOverride); +} + +template +inline void field_axpbyz(const stk::mesh::BulkData& mesh, + const Scalar alpha, + const stk::mesh::FieldBase & xField, + const Scalar beta, + const stk::mesh::FieldBase & yField, const stk::mesh::FieldBase & zField, const stk::mesh::Selector & selector, const EXEC_SPACE& execSpace, - bool IsDeviceExecSpaceUserOverride = (!std::is_same_v)) + bool isDeviceExecSpaceUserOverride = (!std::is_same_v)) { // z = a*x + b*y + ngp_field_blas::impl::field_axpbyz_impl(mesh, alpha, xField, beta, yField, zField, &selector, execSpace, isDeviceExecSpaceUserOverride); +} + +template +inline void field_axpbyz(const stk::mesh::BulkData& mesh, + const Scalar alpha, + const stk::mesh::FieldBase & xField, + const Scalar beta, + const stk::mesh::FieldBase & yField, + const stk::mesh::FieldBase & zField, + const EXEC_SPACE& execSpace, + bool isDeviceExecSpaceUserOverride = (!std::is_same_v)) +{ + // z = a*x + b*y + ngp_field_blas::impl::field_axpbyz_impl(mesh, alpha, xField, beta, yField, zField, nullptr, execSpace, isDeviceExecSpaceUserOverride); +} + + +template +inline void field_axpy(const stk::mesh::BulkData& mesh, + const Scalar alpha, + const stk::mesh::FieldBase & xField, + const stk::mesh::FieldBase & yField, + const stk::mesh::Selector & selector, + const EXEC_SPACE& execSpace, + bool isDeviceExecSpaceUserOverride = (!std::is_same_v)) +{ + Scalar beta = 0; + ngp_field_blas::impl::field_axpbyz_impl(mesh, alpha, xField, beta, yField, yField, &selector, execSpace, isDeviceExecSpaceUserOverride); +} + +template +inline void field_axpy(const stk::mesh::BulkData& mesh, + const Scalar alpha, + const stk::mesh::FieldBase & xField, + const stk::mesh::FieldBase & yField, + const EXEC_SPACE& execSpace, + bool isDeviceExecSpaceUserOverride = (!std::is_same_v)) +{ + Scalar beta = 0; + ngp_field_blas::impl::field_axpbyz_impl(mesh, alpha, xField, beta, yField, yField, nullptr, execSpace, isDeviceExecSpaceUserOverride); +} - if constexpr (ngp_field_blas::impl::operate_on_ngp_mesh()) { - ngp_field_blas::impl::apply_functor_on_field( - mesh, zField, xField, yField, alpha, beta, selector); - } - else { - xField.sync_to_host(); - yField.sync_to_host(); - stk::mesh::field_copy(yField, zField, selector); - stk::mesh::field_axpby(alpha, xField, beta, zField, selector); - } - - zField.clear_sync_state(); - if (ngp_field_blas::impl::mark_modified_on_device(execSpace, IsDeviceExecSpaceUserOverride)) { - zField.modify_on_device(); - } - else { - zField.modify_on_host(); - } +template +inline void field_product(const stk::mesh::BulkData& mesh, + const stk::mesh::FieldBase & xField, + const stk::mesh::FieldBase & yField, + const stk::mesh::FieldBase & zField, + const stk::mesh::Selector & selector, + const EXEC_SPACE& execSpace, + bool isDeviceExecSpaceUserOverride = (!std::is_same_v)) +{ + ngp_field_blas::impl::field_product_impl(mesh, xField, yField, zField, &selector, execSpace, isDeviceExecSpaceUserOverride); +} + +template +inline void field_product(const stk::mesh::BulkData& mesh, + const stk::mesh::FieldBase & xField, + const stk::mesh::FieldBase & yField, + const stk::mesh::FieldBase & zField, + const EXEC_SPACE& execSpace, + bool isDeviceExecSpaceUserOverride = (!std::is_same_v)) +{ + ngp_field_blas::impl::field_product_impl(mesh, xField, yField, zField, nullptr, execSpace, isDeviceExecSpaceUserOverride); +} + + +template +inline void field_scale(const stk::mesh::BulkData& mesh, + const Scalar alpha, + const stk::mesh::FieldBase & xField, + const stk::mesh::Selector & selector, + const EXEC_SPACE& execSpace, + bool isDeviceExecSpaceUserOverride = (!std::is_same_v)) +{ + ngp_field_blas::impl::field_scale_impl(mesh, alpha, xField, &selector, execSpace, isDeviceExecSpaceUserOverride); +} + +template +inline void field_scale(const stk::mesh::BulkData& mesh, + const Scalar alpha, + const stk::mesh::FieldBase & xField, + const EXEC_SPACE& execSpace, + bool isDeviceExecSpaceUserOverride = (!std::is_same_v)) +{ + ngp_field_blas::impl::field_scale_impl(mesh, alpha, xField, nullptr, execSpace, isDeviceExecSpaceUserOverride); +} + +template +inline void field_swap(const stk::mesh::BulkData& mesh, + const stk::mesh::FieldBase & xField, + const stk::mesh::FieldBase & yField, + const stk::mesh::Selector & selector, + const EXEC_SPACE& execSpace, + bool isDeviceExecSpaceUserOverride = (!std::is_same_v)) +{ + ngp_field_blas::impl::field_swap_impl(mesh, xField, yField, &selector, execSpace, isDeviceExecSpaceUserOverride); +} + +template +inline void field_swap(const stk::mesh::BulkData& mesh, + const stk::mesh::FieldBase & xField, + const stk::mesh::FieldBase & yField, + const EXEC_SPACE& execSpace, + bool isDeviceExecSpaceUserOverride = (!std::is_same_v)) +{ + ngp_field_blas::impl::field_swap_impl(mesh, xField, yField, nullptr, execSpace, isDeviceExecSpaceUserOverride); } } // mesh diff --git a/packages/stk/stk_mesh/stk_mesh/base/Selector.cpp b/packages/stk/stk_mesh/stk_mesh/base/Selector.cpp index d92ae7446dc0..99fb4d2960b6 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/Selector.cpp +++ b/packages/stk/stk_mesh/stk_mesh/base/Selector.cpp @@ -583,18 +583,40 @@ stk::mesh::Selector Selector::clone_for_different_mesh(const stk::mesh::MetaData newSelector.m_meta = &differentMeta; for(SelectorNode &selectorNode : newSelector.m_expr) { - if(selectorNode.m_type == SelectorNodeType::PART) + if(selectorNode.m_type == SelectorNodeType::PART || + selectorNode.m_type == SelectorNodeType::PART_UNION || + selectorNode.m_type == SelectorNodeType::PART_INTERSECTION) { + if (selectorNode.m_partOrd != InvalidPartOrdinal) { const std::string& oldPartName = oldMeta.get_part(selectorNode.part()).name(); Part* differentPart = differentMeta.get_part(oldPartName); STK_ThrowRequireMsg(differentPart != nullptr, "Attempting to clone selector into mesh with different parts"); selectorNode.m_partOrd = differentPart->mesh_meta_data_ordinal(); + } + for(unsigned i=0; imesh_meta_data_ordinal(); + } + } + for(unsigned i=0; imesh_meta_data_ordinal(); + } + } } else if(selectorNode.m_type == SelectorNodeType::FIELD) { unsigned ord = selectorNode.m_field_ptr->mesh_meta_data_ordinal(); STK_ThrowRequireMsg(selectorNode.m_field_ptr->name() == differentMeta.get_fields()[ord]->name(), - "Attepting to clone selector into mesh with different parts"); + "Attepting to clone selector into mesh with different fields"); selectorNode.m_field_ptr = differentMeta.get_fields()[ord]; } } diff --git a/packages/stk/stk_mesh/stk_mesh/base/Selector.hpp b/packages/stk/stk_mesh/stk_mesh/base/Selector.hpp index 1c368a1a5adc..0d0c686000f3 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/Selector.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/Selector.hpp @@ -47,9 +47,6 @@ namespace stk { namespace mesh { class BulkData; } } namespace stk { namespace mesh { class Bucket; } } namespace stk { namespace mesh { class FieldBase; } } - - - namespace stk { namespace mesh { // Specify what can be in the expression tree of a selector diff --git a/packages/stk/stk_mesh/stk_mesh/base/SideSetHelper.hpp b/packages/stk/stk_mesh/stk_mesh/base/SideSetHelper.hpp index d005e6825b07..39ef8455df78 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/SideSetHelper.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/SideSetHelper.hpp @@ -38,7 +38,6 @@ #include #include #include -#include #include #include #include diff --git a/packages/stk/stk_mesh/stk_mesh/base/SideSetUtil.cpp b/packages/stk/stk_mesh/stk_mesh/base/SideSetUtil.cpp index b7c9d1f2173f..580d6412a9f8 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/SideSetUtil.cpp +++ b/packages/stk/stk_mesh/stk_mesh/base/SideSetUtil.cpp @@ -11,6 +11,7 @@ #include "stk_mesh/baseImpl/elementGraph/GraphEdgeData.hpp" #include "stk_mesh/baseImpl/elementGraph/GraphTypes.hpp" #include "stk_mesh/baseImpl/SideSetUtilImpl.hpp" +#include "stk_mesh/baseImpl/MeshImplUtils.hpp" #include "stk_util/util/SortAndUnique.hpp" #include "stk_util/diag/StringUtil.hpp" // for Type, etc #include "stk_util/util/string_case_compare.hpp" @@ -182,9 +183,9 @@ bool is_face_represented_in_sideset(const stk::mesh::BulkData& bulk, const stk:: STK_ThrowRequire(bulk.entity_rank(face) == sideRank); std::vector side_elements; - std::vector side_nodes(bulk.begin_nodes(face), bulk.end_nodes(face)); - stk::mesh::get_entities_through_relations(bulk, side_nodes, stk::topology::ELEMENT_RANK, side_elements); + stk::mesh::ConnectedEntities sideNodes = bulk.get_connected_entities(face, stk::topology::NODE_RANK); + stk::mesh::impl::find_entities_these_nodes_have_in_common(bulk, stk::topology::ELEM_RANK, sideNodes.size(), sideNodes.data(), side_elements); bool found = false; diff --git a/packages/stk/stk_mesh/stk_mesh/base/SkinMeshUtil.hpp b/packages/stk/stk_mesh/stk_mesh/base/SkinMeshUtil.hpp index 5968d61e7c9d..9fb50661a380 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/SkinMeshUtil.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/SkinMeshUtil.hpp @@ -35,7 +35,6 @@ #ifndef SKINMESHUTIL_HPP_ #define SKINMESHUTIL_HPP_ -#include #include #include #include @@ -47,6 +46,8 @@ namespace stk { namespace mesh { +class BulkData; + class SkinMeshUtil { public: SkinMeshUtil(ElemElemGraph& elemElemGraph, diff --git a/packages/stk/stk_mesh/stk_mesh/base/Types.hpp b/packages/stk/stk_mesh/stk_mesh/base/Types.hpp index 510964a5ff2c..4457df3f1796 100644 --- a/packages/stk/stk_mesh/stk_mesh/base/Types.hpp +++ b/packages/stk/stk_mesh/stk_mesh/base/Types.hpp @@ -83,18 +83,8 @@ typedef std::vector< unsigned > PermutationIndexVector; typedef std::vector EntityVector; typedef std::vector EntityKeyVector; -template< typename Scalar = void , - class Tag1 = void , class Tag2 = void , - class Tag3 = void , class Tag4 = void , - class Tag5 = void , class Tag6 = void , - class Tag7 = void > - class Field ; - -/** \brief Maximum - * \ref "multi-dimensional array" dimension of a - * \ref stk::mesh::Field "field" - */ -enum { MaximumFieldDimension = 7 }; +template +class Field; enum class Operation { @@ -118,28 +108,9 @@ enum EntityState : char { Unchanged = 0 , Created = 1 , Modified = 2 , Deleted = 3 }; -inline -std::ostream& operator<<(std::ostream& os, EntityState state) -{ - switch(state) { - case Unchanged: os<<"Unchanged"; break; - case Created: os<<"Created"; break; - case Modified: os<<"Modified"; break; - case Deleted: os<<"Deleted"; break; - default: break; - }; - return os; -} - -template< class FieldType > struct STK_DEPRECATED FieldTraits ; - -namespace legacy { -template< class FieldType > struct FieldTraits; -} - +// //MeshIndex describes an Entity's location in the mesh, specifying which bucket, //and the offset (ordinal) into that bucket. -//Ultimately we want this struct to contain two ints rather than a pointer and an int... struct MeshIndex { Bucket* bucket; diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/DeleteEdgesOnFaces.cpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/DeleteEdgesOnFaces.cpp index 9386303e9cf9..525e5a70ff78 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/DeleteEdgesOnFaces.cpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/DeleteEdgesOnFaces.cpp @@ -1,8 +1,11 @@ #include "stk_mesh/baseImpl/DeleteEdgesOnFaces.hpp" #include -#include "stk_mesh/base/GetEntities.hpp" #include "stk_mesh/base/Types.hpp" +#include "stk_mesh/base/GetEntities.hpp" +#include "stk_mesh/base/Entity.hpp" +#include "stk_mesh/base/BulkData.hpp" +#include "stk_mesh/base/Selector.hpp" #include "stk_topology/topology.hpp" #include "stk_util/util/ReportHandler.hpp" @@ -79,4 +82,4 @@ int EdgesOnFacesDeleter::get_edge_idx(Entity face, Entity edge) } } -} \ No newline at end of file +} diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/DeleteEdgesOnFaces.hpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/DeleteEdgesOnFaces.hpp index 2552af61e9ef..9d1bf2346533 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/DeleteEdgesOnFaces.hpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/DeleteEdgesOnFaces.hpp @@ -35,10 +35,15 @@ #ifndef STK_MESH_BASE_DELETE_EDGES_ON_FACES_H #define STK_MESH_BASE_DELETE_EDGES_ON_FACES_H -#include "stk_mesh/base/BulkData.hpp" +#include namespace stk { namespace mesh { + +class BulkData; +class Part; +struct Entity; + namespace impl { class EdgesOnFacesDeleter @@ -64,4 +69,4 @@ class EdgesOnFacesDeleter } } -#endif \ No newline at end of file +#endif diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/ElementTopologyDeletions.cpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/ElementTopologyDeletions.cpp index d6f5e76232eb..3b0c0afa893b 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/ElementTopologyDeletions.cpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/ElementTopologyDeletions.cpp @@ -1,4 +1,6 @@ #include +#include +#include namespace stk { diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/ElementTopologyDeletions.hpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/ElementTopologyDeletions.hpp index d65d863fb38f..fc5283f49aa8 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/ElementTopologyDeletions.hpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/ElementTopologyDeletions.hpp @@ -1,16 +1,50 @@ +// Copyright 2002 - 2008, 2010, 2011 National Technology Engineering +// Solutions of Sandia, LLC (NTESS). Under the terms of Contract +// DE-NA0003525 with NTESS, the U.S. Government retains certain rights +// in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following +// disclaimer in the documentation and/or other materials provided +// with the distribution. +// +// * Neither the name of NTESS nor the names of its contributors +// may be used to endorse or promote products derived from this +// software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + #ifndef STK_ELEMENT_TOPOLOGY_DELETIONS_H #define STK_ELEMENT_TOPOLOGY_DELETIONS_H -#include "stk_mesh/base/BulkData.hpp" +#include "stk_mesh/base/Types.hpp" #include "stk_mesh/base/Entity.hpp" -#include "stk_mesh/base/MetaData.hpp" #include "stk_mesh/base/Selector.hpp" -#include "stk_mesh/base/Types.hpp" namespace stk { namespace mesh { +class MetaData; +class BulkData; + namespace impl { @@ -44,37 +78,4 @@ class ElementTopologyDeletions } } -// Copyright 2002 - 2008, 2010, 2011 National Technology Engineering -// Solutions of Sandia, LLC (NTESS). Under the terms of Contract -// DE-NA0003525 with NTESS, the U.S. Government retains certain rights -// in this software. -// - // Redistribution and use in source and binary forms, with or without - // modification, are permitted provided that the following conditions are - // met: - // - // * Redistributions of source code must retain the above copyright - // notice, this list of conditions and the following disclaimer. - // - // * Redistributions in binary form must reproduce the above - // copyright notice, this list of conditions and the following - // disclaimer in the documentation and/or other materials provided - // with the distribution. - // -// * Neither the name of NTESS nor the names of its contributors -// may be used to endorse or promote products derived from this -// software without specific prior written permission. -// - // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS - // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT - // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR - // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT - // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, - // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT - // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, - // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY - // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE - // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - #endif diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityGhostData.hpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityGhostData.hpp index 056bbf2804e9..b038bcd24cbd 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityGhostData.hpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityGhostData.hpp @@ -37,10 +37,10 @@ #include #include -#include namespace stk { namespace mesh { +class BulkData; namespace impl { diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityKeyMapping.cpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityKeyMapping.cpp index 33d36f4aac81..f510717d7ea9 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityKeyMapping.cpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/EntityKeyMapping.cpp @@ -244,7 +244,15 @@ EntityKeyMapping::internal_create_entity( const EntityKey & key) Entity EntityKeyMapping::get_entity(const EntityKey &key) const { + if (!key.is_valid()) { + return Entity(); + } + EntityRank rank = key.rank(); + if (rank >= entity_rank_count()) { + return Entity(); + } + if (!m_destroy_cache[rank].empty()) { const std::vector& destroyed = m_destroy_cache[rank]; if (destroyed.size() < 64) { @@ -260,13 +268,6 @@ Entity EntityKeyMapping::get_entity(const EntityKey &key) const clear_updated_entity_cache(rank); - STK_ThrowErrorMsgIf( ! key.is_valid(), - "Invalid key: " << key.rank() << " " << key.id()); - - if (rank >= entity_rank_count()) { - return Entity(); - } - entity_iterator ent = get_from_cache(key); if (ent != m_create_cache[rank].end()) { return ent->second; diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/FieldRepository.cpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/FieldRepository.cpp index 78bcf4bbb8dd..1bdc38ee1b46 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/FieldRepository.cpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/FieldRepository.cpp @@ -37,7 +37,6 @@ #include // for ostringstream #include // for operator<<, basic_ostream, etc #include // for equal_case -#include "Shards_Array.hpp" // for ArrayDimTag #include "stk_mesh/base/DataTraits.hpp" // for DataTraits #include "stk_mesh/base/MetaData.hpp" #include "stk_mesh/base/FieldState.hpp" // for ::MaximumFieldStates, etc @@ -51,19 +50,10 @@ namespace impl { namespace { -std::string print_field_type(const DataTraits & arg_traits , - unsigned arg_rank , - const shards::ArrayDimTag * const * arg_tags ) +std::string print_field_type(const DataTraits & arg_traits) { - STK_ThrowRequireMsg(arg_rank < 8, "Invalid field rank: " << arg_rank); - std::ostringstream oss; - oss << "FieldBase<" ; - oss << arg_traits.name ; - for ( unsigned i = 0 ; i < arg_rank ; ++i ) { - oss << "," << arg_tags[i]->name(); - } - oss << ">" ; + oss << "Field<" << arg_traits.name << ">"; return oss.str(); } @@ -75,11 +65,9 @@ std::string print_field_type(const DataTraits & arg_traits , // 3) Dimension must be different by at most one rank, // where the tags match for the smaller rank. void -FieldRepository::verify_field_type(const FieldBase & arg_field, - const DataTraits & arg_traits, - unsigned arg_rank, - const shards::ArrayDimTag * const * arg_dim_tags, - unsigned arg_num_states) const +FieldRepository::verify_field_type(const FieldBase & arg_field, + const DataTraits & arg_traits, + unsigned arg_num_states) const { const bool ok_traits = arg_traits.is_void || &arg_traits == &arg_field.data_traits(); @@ -88,25 +76,21 @@ FieldRepository::verify_field_type(const FieldBase & arg_field STK_ThrowErrorMsgIf(not ok_traits || not ok_number_states, " verify_field_type FAILED: Existing field = " << - print_field_type(arg_field.data_traits(), arg_field.m_field_rank, arg_field.m_dim_tags) << - "[ name = \"" << arg_field.name() << + print_field_type(arg_field.data_traits()) << "[ name = \"" << arg_field.name() << "\" , #states = " << arg_field.number_of_states() << " ]" << - " Expected field info = " << - print_field_type(arg_traits, arg_rank, arg_dim_tags) << + " Expected field info = " << print_field_type(arg_traits) << "[ #states = " << arg_num_states << " ]"); } //---------------------------------------------------------------------- -FieldBase * FieldRepository::get_field(stk::topology::rank_t arg_entity_rank, - const std::string & arg_name, - const DataTraits & arg_traits, - unsigned arg_array_rank, - const shards::ArrayDimTag * const * arg_dim_tags, - unsigned arg_num_states) const +FieldBase * FieldRepository::get_field(stk::topology::rank_t arg_entity_rank, + const std::string & arg_name, + const DataTraits & arg_traits, + unsigned arg_num_states) const { for (FieldBase * field : m_rankedFields[arg_entity_rank]) { if (equal_case(field->name(), arg_name)) { - verify_field_type(*field, arg_traits, arg_array_rank, arg_dim_tags, arg_num_states); + verify_field_type(*field, arg_traits, arg_num_states); return field; } } diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/FieldRepository.hpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/FieldRepository.hpp index f2c1d375f30b..4da062f9270d 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/FieldRepository.hpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/FieldRepository.hpp @@ -43,7 +43,6 @@ #include "stk_mesh/base/Selector.hpp" // for Selector #include "stk_topology/topology.hpp" // for topology, topology::rank_t, etc #include "stk_util/util/ReportHandler.hpp" // for ThrowAssert -namespace shards { class ArrayDimTag; } namespace stk { namespace mesh { class DataTraits; } } namespace stk { namespace mesh { class MetaData; } } namespace stk { namespace mesh { class Part; } } @@ -67,14 +66,10 @@ class FieldRepository { m_rankedFields[new_field->entity_rank()].push_back(new_field); } - FieldBase * get_field( - stk::topology::rank_t arg_entity_rank , - const std::string & arg_name , - const DataTraits & arg_traits , - unsigned arg_array_rank , - const shards::ArrayDimTag * const * arg_dim_tags , - unsigned arg_num_states - ) const; + FieldBase * get_field(stk::topology::rank_t arg_entity_rank, + const std::string & arg_name, + const DataTraits & arg_traits, + unsigned arg_num_states) const; void verify_and_clean_restrictions(const Part& superset, const Part& subset); @@ -106,7 +101,7 @@ class FieldRepository { template bool remove_attribute( FieldBase & f , const T * a ) { - return f.remove_attribute( a ); + return f.remove_attribute( a ); } void declare_field_restriction( @@ -136,11 +131,9 @@ class FieldRepository { MetaData & mesh_meta_data() { return m_meta; } private: - void verify_field_type(const FieldBase & arg_field, - const DataTraits & arg_traits, - unsigned arg_rank, - const shards::ArrayDimTag * const * arg_dim_tags, - unsigned arg_num_states) const; + void verify_field_type(const FieldBase & arg_field, + const DataTraits & arg_traits, + unsigned arg_num_states) const; MetaData & m_meta; FieldVector m_fields; diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/ForEachEntityLoopAbstractions.hpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/ForEachEntityLoopAbstractions.hpp index ba7a470f30ff..8eb28dc32016 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/ForEachEntityLoopAbstractions.hpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/ForEachEntityLoopAbstractions.hpp @@ -39,14 +39,15 @@ #include // for Bucket #include #include +#include namespace stk { namespace mesh { namespace impl { -inline +template void for_each_selected_entity_run(const BulkData &mesh, stk::topology::rank_t rank, const Selector &selector, - const std::function& functor) + const Functor& functor) { const stk::mesh::BucketVector & buckets = mesh.get_buckets(rank, selector); const size_t numBuckets = buckets.size(); @@ -58,26 +59,12 @@ void for_each_selected_entity_run(const BulkData &mesh, stk::topology::rank_t ra stk::mesh::Bucket *bucket = buckets[j]; for(size_t i=0; isize(); i++) { + if constexpr (std::is_invocable_v) { functor(mesh, stk::mesh::MeshIndex({bucket,i})); - } - } -} - -inline -void for_each_selected_entity_run(const BulkData &mesh, stk::topology::rank_t rank, const Selector &selector, - const std::function& functor) -{ - const stk::mesh::BucketVector & buckets = mesh.get_buckets(rank, selector); - const size_t numBuckets = buckets.size(); -#ifdef _OPENMP -#pragma omp parallel for -#endif - for(size_t j=0; jsize(); i++) - { + } + else { functor(mesh, (*bucket)[i]); + } } } } @@ -106,8 +93,8 @@ void for_each_selected_entity_run_with_nodes(const BulkData &mesh, stk::topology } } -inline -void for_each_selected_entity_run_no_threads(const BulkData &mesh, stk::topology::rank_t rank, const stk::mesh::Selector &selector, const std::function &functor) +template +void for_each_selected_entity_run_no_threads(const BulkData &mesh, stk::topology::rank_t rank, const stk::mesh::Selector &selector, const Functor &functor) { const stk::mesh::BucketVector & buckets = mesh.get_buckets(rank, selector); for(size_t j=0; jsize(); i++) { + if constexpr (std::is_invocable_v) { functor(mesh, stk::mesh::MeshIndex({bucket,i})); - } - } -} - -inline -void for_each_selected_entity_run_no_threads(const BulkData &mesh, stk::topology::rank_t rank, const Selector &selector, - const std::function& functor) -{ - const stk::mesh::BucketVector & buckets = mesh.get_buckets(rank, selector); - const size_t numBuckets = buckets.size(); - - for(size_t j=0; jsize(); i++) - { + } + else { functor(mesh, (*bucket)[i]); + } } } } diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/MeshCommVerify.hpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/MeshCommVerify.hpp index 43a6b797de5c..d07852bfb0b7 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/MeshCommVerify.hpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/MeshCommVerify.hpp @@ -39,7 +39,6 @@ #include #include -#include #include #include @@ -50,6 +49,8 @@ namespace stk { namespace mesh { +class BulkData; + namespace impl { //---------------------------------------------------------------------- diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/MeshModLogObserver.cpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/MeshModLogObserver.cpp index 6bf79c36b1cb..492baed7f814 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/MeshModLogObserver.cpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/MeshModLogObserver.cpp @@ -32,11 +32,11 @@ // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include +#include #include #include #include #include -#include namespace stk { namespace mesh { diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/NgpFieldBLASImpl.hpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/NgpFieldBLASImpl.hpp index bbcc1183b770..60b2f0b264f1 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/NgpFieldBLASImpl.hpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/NgpFieldBLASImpl.hpp @@ -6,15 +6,15 @@ // Redistribution and use in source and binary forms, with or without // modification, are permitted provided that the following conditions are // met: -// +// // * Redistributions of source code must retain the above copyright // notice, this list of conditions and the following disclaimer. -// +// // * Redistributions in binary form must reproduce the above // copyright notice, this list of conditions and the following // disclaimer in the documentation and/or other materials provided // with the distribution. -// +// // * Neither the name of NTESS nor the names of its contributors // may be used to endorse or promote products derived from this // software without specific prior written permission. @@ -30,7 +30,7 @@ // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -// +// #ifndef STK_MESH_BASEIMPL_NGPFIELDBLASIMPL_HPP #define STK_MESH_BASEIMPL_NGPFIELDBLASIMPL_HPP @@ -87,6 +87,18 @@ bool mark_modified_on_device( return operate_on_ngp_mesh() || isDeviceExecSpaceUserOverride; } +template +void mark_field_modified(const mesh::FieldBase& field, EXEC_SPACE execSpace, bool isDeviceExecSpaceUserOverride) +{ + field.clear_sync_state(); + if (ngp_field_blas::impl::mark_modified_on_device(execSpace, isDeviceExecSpaceUserOverride)) { + field.modify_on_device(); + } + else { + field.modify_on_host(); + } +} + template class FieldFill { public: @@ -156,9 +168,9 @@ void field_fill_impl(const Scalar alpha, field.clear_sync_state(); std::unique_ptr fieldSelector; - if (selectorPtr == nullptr) { + if (selectorPtr == nullptr) { fieldSelector = std::make_unique(field); - } + } const stk::mesh::Selector& selector = selectorPtr != nullptr ? *selectorPtr : *(fieldSelector.get()); if constexpr (operate_on_ngp_mesh()) { @@ -172,12 +184,7 @@ void field_fill_impl(const Scalar alpha, field_fill_for_each_entity(hostMesh, hostField, alpha, component, selector, execSpace); } - if (mark_modified_on_device(execSpace, isDeviceExecSpaceUserOverride)) { - field.modify_on_device(); - } - else { - field.modify_on_host(); - } + mark_field_modified(field, execSpace, isDeviceExecSpaceUserOverride); } template @@ -258,31 +265,26 @@ void field_copy_impl(const stk::mesh::FieldBase& xField, bool isDeviceExecSpaceUserOverride) { std::unique_ptr fieldSelector; - if (selectorPtr == nullptr) { + if (selectorPtr == nullptr) { fieldSelector = std::make_unique(stk::mesh::Selector(xField) & stk::mesh::Selector(yField)); } const stk::mesh::Selector& selector = selectorPtr != nullptr ? *selectorPtr : *(fieldSelector.get()); field_copy_no_mark_mod(xField, yField, selector, execSpace); - yField.clear_sync_state(); - if (mark_modified_on_device(execSpace, isDeviceExecSpaceUserOverride)) { - yField.modify_on_device(); - } - else { - yField.modify_on_host(); - } + mark_field_modified(yField, execSpace, isDeviceExecSpaceUserOverride); } -template +template void apply_functor_on_field(const stk::mesh::BulkData& mesh, const stk::mesh::FieldBase & zField, const stk::mesh::FieldBase & xField, const stk::mesh::FieldBase & yField, - const DataType alpha, - const DataType beta, + typename Functor::value_type alpha, + typename Functor::value_type beta, const stk::mesh::Selector & select) { + using DataType = typename Functor::value_type; const stk::mesh::Selector selector = select & stk::mesh::selectField(zField) & stk::mesh::selectField(xField) & stk::mesh::selectField(yField); stk::mesh::EntityRank entityRank = zField.entity_rank(); @@ -298,24 +300,26 @@ void apply_functor_on_field(const stk::mesh::BulkData& mesh, zField.modify_on_device(); } - -struct FieldAXPBYFunctor +template +struct FieldAXPBYZFunctor { - FieldAXPBYFunctor(stk::mesh::NgpMesh & my_ngp_mesh, - stk::mesh::NgpField & output_field, - stk::mesh::NgpField & input_x, - stk::mesh::NgpField & input_y, - const double & alpha, - const double & beta) + using value_type = Scalar; + + FieldAXPBYZFunctor(stk::mesh::NgpMesh & my_ngp_mesh, + stk::mesh::NgpField & output_field, + stk::mesh::NgpField & input_x, + stk::mesh::NgpField & input_y, + const Scalar & alpha, + const Scalar & beta) : my_mesh(my_ngp_mesh), output(output_field), x(input_x), y(input_y), a(alpha), b(beta) { } stk::mesh::NgpMesh my_mesh; - stk::mesh::NgpField output; - stk::mesh::NgpField x; - stk::mesh::NgpField y; - const double a; - const double b; + stk::mesh::NgpField output; + stk::mesh::NgpField x; + stk::mesh::NgpField y; + const Scalar a; + const Scalar b; KOKKOS_FUNCTION void operator()(stk::mesh::FastMeshIndex f) const { @@ -331,6 +335,317 @@ struct FieldAXPBYFunctor } }; +template +struct FieldProductFunctor +{ + using value_type = Scalar; + + FieldProductFunctor(stk::mesh::NgpMesh & my_ngp_mesh, + stk::mesh::NgpField & output_field, + stk::mesh::NgpField & input_x, + stk::mesh::NgpField & input_y, + const Scalar a, + const Scalar b) + : my_mesh(my_ngp_mesh), + output(output_field), + x(input_x), y(input_y) + { + } + stk::mesh::NgpMesh my_mesh; + stk::mesh::NgpField output; + stk::mesh::NgpField x; + stk::mesh::NgpField y; + + KOKKOS_FUNCTION + void operator()(stk::mesh::FastMeshIndex f) const + { + unsigned num_components = output.get_num_components_per_entity(f); + unsigned other = x.get_num_components_per_entity(f); + num_components = (other < num_components) ? other : num_components; + other = y.get_num_components_per_entity(f); + num_components = (other < num_components) ? other : num_components; + for (unsigned i = 0; i < num_components; ++i) + { + output.get(f, i) = x.get(f, i) * y.get(f, i); + } + } +}; + +template +struct FieldScaleFunctor +{ + using value_type = Scalar; + + FieldScaleFunctor(stk::mesh::NgpMesh & my_ngp_mesh, + stk::mesh::NgpField & output_field, + stk::mesh::NgpField & unused_field1, + stk::mesh::NgpField & unused_field2, + const Scalar a, + const Scalar b) + : my_mesh(my_ngp_mesh), + output(output_field), + alpha(a), + m_unused_field1(unused_field1), m_unused_field2(unused_field2), + m_beta(b) + + { + } + stk::mesh::NgpMesh my_mesh; + stk::mesh::NgpField output; + Scalar alpha; + + KOKKOS_FUNCTION + void operator()(stk::mesh::FastMeshIndex f) const + { + unsigned num_components = output.get_num_components_per_entity(f); + for (unsigned i = 0; i < num_components; ++i) + { + output.get(f, i) = alpha * output.get(f, i); + } + } + +private: + stk::mesh::NgpField m_unused_field1; + stk::mesh::NgpField m_unused_field2; + Scalar m_beta; +}; + +template +struct FieldSwapFunctor +{ + using value_type = Scalar; + + FieldSwapFunctor(stk::mesh::NgpMesh & my_ngp_mesh, + stk::mesh::NgpField & xFieldInput, + stk::mesh::NgpField & yFieldInput, + stk::mesh::NgpField & unused_field2, + const Scalar a, + const Scalar b) + : my_mesh(my_ngp_mesh), + xField(xFieldInput), + yField(yFieldInput), + m_unused_field2(unused_field2), + m_alpha(a), + m_beta(b) + + { + } + + stk::mesh::NgpMesh my_mesh; + stk::mesh::NgpField xField; + stk::mesh::NgpField yField; + + KOKKOS_FUNCTION + void operator()(stk::mesh::FastMeshIndex f) const + { + unsigned num_components = xField.get_num_components_per_entity(f); + unsigned other = yField.get_num_components_per_entity(f); + num_components = (other < num_components) ? other : num_components; + for (unsigned i = 0; i < num_components; ++i) + { + Scalar tmp = xField.get(f, i); + xField.get(f, i) = yField.get(f, i); + yField.get(f, i) = tmp; + } + } + +private: + stk::mesh::NgpField m_unused_field2; + Scalar m_alpha; + Scalar m_beta; +}; + +template +inline void field_axpbyz_impl(const stk::mesh::BulkData& mesh, + const DataType alpha, + const stk::mesh::FieldBase & xField, + const DataType beta, + const stk::mesh::FieldBase & yField, + const stk::mesh::FieldBase & zField, + const stk::mesh::Selector* selectorPtr, + const EXEC_SPACE& execSpace, + bool isDeviceExecSpaceUserOverride) +{ + const stk::mesh::DataTraits& dataTraits = xField.data_traits(); + STK_ThrowRequireMsg(dataTraits == yField.data_traits(), "xField and yField must have same datatype"); + STK_ThrowRequireMsg(dataTraits == zField.data_traits(), "xField and zField must have same datatype"); + + stk::mesh::Selector fieldSelector; + if (selectorPtr == nullptr) { + fieldSelector = stk::mesh::Selector(xField) & stk::mesh::Selector(yField); + } else + { + fieldSelector = *selectorPtr; + } + + if constexpr (ngp_field_blas::impl::operate_on_ngp_mesh()) { + if (dataTraits.type_info == typeid(double)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, zField, xField, yField, alpha, beta, fieldSelector); + } else if (dataTraits.type_info == typeid(float)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, zField, xField, yField, alpha, beta, fieldSelector); + } else if (dataTraits.type_info == typeid(int)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, zField, xField, yField, alpha, beta, fieldSelector); + } else if (dataTraits.type_info == typeid(unsigned)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, zField, xField, yField, alpha, beta, fieldSelector); + } else + { + STK_ThrowErrorMsg("axpbyz doesn't yet support fields of type "< +inline void field_product_impl(const stk::mesh::BulkData& mesh, + const stk::mesh::FieldBase & xField, + const stk::mesh::FieldBase & yField, + const stk::mesh::FieldBase & zField, + const stk::mesh::Selector* selectorPtr, + const EXEC_SPACE& execSpace, + bool isDeviceExecSpaceUserOverride) +{ + const stk::mesh::DataTraits& dataTraits = xField.data_traits(); + STK_ThrowRequireMsg(dataTraits == yField.data_traits(), "xField and yField must have same datatype"); + STK_ThrowRequireMsg(dataTraits == zField.data_traits(), "xField and zField must have same datatype"); + + stk::mesh::Selector fieldSelector; + if (selectorPtr == nullptr) { + fieldSelector = stk::mesh::Selector(xField) & stk::mesh::Selector(yField); + } else + { + fieldSelector = *selectorPtr; + } + + if constexpr (ngp_field_blas::impl::operate_on_ngp_mesh()) { + if (dataTraits.type_info == typeid(double)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, zField, xField, yField, 1, 1, fieldSelector); + } else if (dataTraits.type_info == typeid(float)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, zField, xField, yField, 1, 1, fieldSelector); + } else if (dataTraits.type_info == typeid(int)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, zField, xField, yField, 1, 1, fieldSelector); + } else if (dataTraits.type_info == typeid(unsigned)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, zField, xField, yField, 1, 1, fieldSelector); + } else + { + STK_ThrowErrorMsg("field_product doesn't yet support fields of type "< +inline void field_scale_impl(const stk::mesh::BulkData& mesh, + const Scalar alpha, + const stk::mesh::FieldBase & xField, + const stk::mesh::Selector* selectorPtr, + const EXEC_SPACE& execSpace, + bool isDeviceExecSpaceUserOverride) +{ + const stk::mesh::DataTraits& dataTraits = xField.data_traits(); + + stk::mesh::Selector fieldSelector; + if (selectorPtr == nullptr) { + fieldSelector = stk::mesh::Selector(xField); + } else + { + fieldSelector = *selectorPtr; + } + + if constexpr (ngp_field_blas::impl::operate_on_ngp_mesh()) { + if (dataTraits.type_info == typeid(double)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, xField, xField, xField, alpha, alpha, fieldSelector); + } else if (dataTraits.type_info == typeid(float)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, xField, xField, xField, alpha, alpha, fieldSelector); + } else if (dataTraits.type_info == typeid(int)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, xField, xField, xField, alpha, alpha, fieldSelector); + } else if (dataTraits.type_info == typeid(unsigned)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, xField, xField, xField, alpha, alpha, fieldSelector); + } else + { + STK_ThrowErrorMsg("field_product doesn't yet support fields of type "< +inline void field_swap_impl(const stk::mesh::BulkData& mesh, + const stk::mesh::FieldBase & xField, + const stk::mesh::FieldBase & yField, + const stk::mesh::Selector* selectorPtr, + const EXEC_SPACE& execSpace, + bool isDeviceExecSpaceUserOverride) +{ + const stk::mesh::DataTraits& dataTraits = xField.data_traits(); + STK_ThrowRequireMsg(dataTraits == yField.data_traits(), "xField and yField must have same datatype"); + + stk::mesh::Selector fieldSelector; + if (selectorPtr == nullptr) { + fieldSelector = stk::mesh::Selector(xField); + } else + { + fieldSelector = *selectorPtr; + } + + if constexpr (ngp_field_blas::impl::operate_on_ngp_mesh()) { + if (dataTraits.type_info == typeid(double)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, xField, yField, xField, 0, 0, fieldSelector); + } else if (dataTraits.type_info == typeid(float)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, xField, xField, xField, 0, 0, fieldSelector); + } else if (dataTraits.type_info == typeid(int)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, xField, xField, xField, 0, 0, fieldSelector); + } else if (dataTraits.type_info == typeid(unsigned)) { + ngp_field_blas::impl::apply_functor_on_field>( + mesh, xField, yField, xField, 0, 0, fieldSelector); + } else + { + STK_ThrowErrorMsg("field_product doesn't yet support fields of type "< +#include + +namespace stk { +namespace mesh { + +inline +std::ostream& operator<<(std::ostream& os, EntityState state) +{ + switch(state) { + case Unchanged: os<<"Unchanged"; break; + case Created: os<<"Created"; break; + case Modified: os<<"Modified"; break; + case Deleted: os<<"Deleted"; break; + default: break; + }; + return os; +} + +} // namespace mesh +} // namespace stk + +#endif diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/SideSetImpl.hpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/SideSetImpl.hpp index 89819c6cb9e3..a33a98edd186 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/SideSetImpl.hpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/SideSetImpl.hpp @@ -42,7 +42,6 @@ #include // for uint16_t #include // for max #include // for less, equal_to -#include // for operator<<, basic_ostream, etc #include // for list #include // for map, map<>::value_compare #include diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/elementGraph/ElemFilter.hpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/elementGraph/ElemFilter.hpp index a4e0b5874f87..3783e03c4340 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/elementGraph/ElemFilter.hpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/elementGraph/ElemFilter.hpp @@ -37,13 +37,14 @@ #include #include #include -#include #include #include #include "BulkDataIdMapper.hpp" #include namespace stk { namespace mesh { +class BulkData; + namespace impl { inline bool is_solid_shell_connection(stk::topology t1, stk::topology t2) diff --git a/packages/stk/stk_mesh/stk_mesh/baseImpl/elementGraph/RemoteSelectorUpdater.hpp b/packages/stk/stk_mesh/stk_mesh/baseImpl/elementGraph/RemoteSelectorUpdater.hpp index fff017f04b06..504bf8c54a55 100644 --- a/packages/stk/stk_mesh/stk_mesh/baseImpl/elementGraph/RemoteSelectorUpdater.hpp +++ b/packages/stk/stk_mesh/stk_mesh/baseImpl/elementGraph/RemoteSelectorUpdater.hpp @@ -35,7 +35,6 @@ #define STK_REMOTE_SELECTOR_UPDATER_HPP #include // for size_t -#include // for BulkData #include #include // for MetaData #include // for topology, etc @@ -51,6 +50,7 @@ #include "ElemElemGraphImpl.hpp" namespace stk { namespace mesh { +class BulkData; class RemoteSelectorUpdater : public stk::mesh::ModificationObserver { diff --git a/packages/stk/stk_middle_mesh/stk_middle_mesh/mesh_gather_to_root.cpp b/packages/stk/stk_middle_mesh/stk_middle_mesh/mesh_gather_to_root.cpp index 11f959543088..d7c503cc0710 100644 --- a/packages/stk/stk_middle_mesh/stk_middle_mesh/mesh_gather_to_root.cpp +++ b/packages/stk/stk_middle_mesh/stk_middle_mesh/mesh_gather_to_root.cpp @@ -271,7 +271,7 @@ void MeshGatherToRoot::unpack_info(utils::impl::ParallelExchange& e for (auto& elInfo : recvBuf) { - std::array edges; + std::array edges = {0}; for (int i = 0; i < 4; ++i) if (elInfo.edgeOwnerRanks[i] != -1) { diff --git a/packages/stk/stk_performance_tests/stk_mesh/multi_block.hpp b/packages/stk/stk_performance_tests/stk_mesh/multi_block.hpp index 6e6a07bcd7f6..4d6fc21210d9 100644 --- a/packages/stk/stk_performance_tests/stk_mesh/multi_block.hpp +++ b/packages/stk/stk_performance_tests/stk_mesh/multi_block.hpp @@ -35,10 +35,11 @@ #ifndef STK_PERFORMANCE_MULTI_BLOCK_HPP #define STK_PERFORMANCE_MULTI_BLOCK_HPP -#include "stk_mesh/base/MetaData.hpp" -#include "stk_mesh/base/BulkData.hpp" - namespace stk { + +namespace mesh { class MetaData; } +namespace mesh { class BulkData; } + namespace performance_tests { void setup_multiple_blocks(stk::mesh::MetaData& meta, unsigned numBlocks); diff --git a/packages/stk/stk_performance_tests/stk_search/SurfaceToSurface.cpp b/packages/stk/stk_performance_tests/stk_search/SurfaceToSurface.cpp index 84b4ef702f81..8598268e934e 100644 --- a/packages/stk/stk_performance_tests/stk_search/SurfaceToSurface.cpp +++ b/packages/stk/stk_performance_tests/stk_search/SurfaceToSurface.cpp @@ -653,12 +653,14 @@ void run_imported_surface_to_surface_test_pll_local_with_views(const std::string Kokkos::deep_copy(diceBoxes, diceBoxesHost); Kokkos::deep_copy(toolBoxes, toolBoxesHost); + std::shared_ptr searchData; for (unsigned run = 0; run < NUM_RUNS; ++run) { batchTimer.start_batch_timer(); for (int i = 0; i < numIterations; ++i) { Kokkos::View searchResults; - stk::search::local_coarse_search(diceBoxes, toolBoxes, searchMethod, searchResults, ExecSpace{}); + constexpr bool sortResults = false; + searchData = stk::search::local_coarse_search(diceBoxes, toolBoxes, searchMethod, searchResults, ExecSpace{}, sortResults, searchData); } batchTimer.stop_batch_timer(); } diff --git a/packages/stk/stk_search/cmake/Dependencies.cmake b/packages/stk/stk_search/cmake/Dependencies.cmake index a8e3e549be78..249a303c0afc 100644 --- a/packages/stk/stk_search/cmake/Dependencies.cmake +++ b/packages/stk/stk_search/cmake/Dependencies.cmake @@ -1,4 +1,4 @@ -SET(LIB_REQUIRED_DEP_PACKAGES STKUtil Kokkos KokkosKernels STKMath) +SET(LIB_REQUIRED_DEP_PACKAGES STKUtil Kokkos STKMath) SET(LIB_OPTIONAL_DEP_PACKAGES) SET(TEST_REQUIRED_DEP_PACKAGES) SET(TEST_OPTIONAL_DEP_PACKAGES) diff --git a/packages/stk/stk_search/stk_search/BoxIdent.hpp b/packages/stk/stk_search/stk_search/BoxIdent.hpp index d3252f74f028..55cf5bb5488c 100644 --- a/packages/stk/stk_search/stk_search/BoxIdent.hpp +++ b/packages/stk/stk_search/stk_search/BoxIdent.hpp @@ -103,6 +103,16 @@ struct IdentProcIntersection }; +template +struct Comparator { + +KOKKOS_FUNCTION +bool operator()(T a, T b) const { + return a < b; +} + +}; + } #endif // BOXIDENT_HPP diff --git a/packages/stk/stk_search/stk_search/CMakeLists.txt b/packages/stk/stk_search/stk_search/CMakeLists.txt index 4a64bd432acc..3dc52cbc2578 100644 --- a/packages/stk/stk_search/stk_search/CMakeLists.txt +++ b/packages/stk/stk_search/stk_search/CMakeLists.txt @@ -66,8 +66,6 @@ else() ${${PROJECT_NAME}_INSTALL_INCLUDE_DIR}/stk_search/arborx) endif() -# find_package(KokkosKernels REQUIRED) -# target_link_libraries(stk_search PUBLIC KokkosKernels::all_libs) endif() target_include_directories(stk_search PUBLIC diff --git a/packages/stk/stk_search/stk_search/CoarseSearch.hpp b/packages/stk/stk_search/stk_search/CoarseSearch.hpp index 4bae8287d1af..d81a24eeed01 100644 --- a/packages/stk/stk_search/stk_search/CoarseSearch.hpp +++ b/packages/stk/stk_search/stk_search/CoarseSearch.hpp @@ -83,12 +83,13 @@ void coarse_search(std::vector> co stk::ParallelMachine comm, std::vector>& intersections, bool enforceSearchResultSymmetry = true, - bool autoSwapDomainAndRange = true) + bool autoSwapDomainAndRange = true, + bool sortSearchResults = false) { switch (method) { case ARBORX: { #ifdef STK_HAS_ARBORX - coarse_search_arborx(domain, range, comm, intersections, enforceSearchResultSymmetry); + coarse_search_arborx(domain, range, comm, intersections, enforceSearchResultSymmetry, sortSearchResults); #else STK_ThrowErrorMsg("STK(stk_search) was not configured with ARBORX enabled. Please use KDTREE or MORTON_LBVH."); #endif @@ -96,15 +97,15 @@ void coarse_search(std::vector> co } case KDTREE: { if (autoSwapDomainAndRange) { - coarse_search_kdtree_driver(domain, range, comm, intersections, enforceSearchResultSymmetry); + coarse_search_kdtree_driver(domain, range, comm, intersections, enforceSearchResultSymmetry, sortSearchResults); } else { - coarse_search_kdtree(domain, range, comm, intersections, enforceSearchResultSymmetry); + coarse_search_kdtree(domain, range, comm, intersections, enforceSearchResultSymmetry, sortSearchResults); } break; } case MORTON_LBVH: { - coarse_search_morton_lbvh(domain, range, comm, intersections, enforceSearchResultSymmetry); + coarse_search_morton_lbvh(domain, range, comm, intersections, enforceSearchResultSymmetry, sortSearchResults); break; } default: { @@ -121,7 +122,8 @@ void coarse_search(DomainView const & domain, ResultView& intersections, ExecutionSpace const& execSpace = ExecutionSpace{}, bool enforceSearchResultSymmetry = true, - bool autoSwapDomainAndRange = true) + bool autoSwapDomainAndRange = true, + bool sortSearchResults = false) { check_coarse_search_types_parallel(); Kokkos::Profiling::pushRegion("STK coarse search with Views"); @@ -129,7 +131,7 @@ void coarse_search(DomainView const & domain, switch (method) { case ARBORX: { #ifdef STK_HAS_ARBORX - coarse_search_arborx(domain, range, comm, intersections, execSpace, enforceSearchResultSymmetry); + coarse_search_arborx(domain, range, comm, intersections, execSpace, enforceSearchResultSymmetry, sortSearchResults); #else STK_ThrowErrorMsg("STK(stk_search) was not configured with ARBORX enabled. Please use KDTREE or MORTON_LBVH."); #endif @@ -140,7 +142,7 @@ void coarse_search(DomainView const & domain, break; } case MORTON_LBVH: { - coarse_search_morton_lbvh(domain, range, comm, intersections, execSpace, enforceSearchResultSymmetry); + coarse_search_morton_lbvh(domain, range, comm, intersections, execSpace, enforceSearchResultSymmetry, sortSearchResults); break; } default: { diff --git a/packages/stk/stk_search/stk_search/CommonSearchUtil.hpp b/packages/stk/stk_search/stk_search/CommonSearchUtil.hpp index 836a088b0853..bb807ebb7559 100644 --- a/packages/stk/stk_search/stk_search/CommonSearchUtil.hpp +++ b/packages/stk/stk_search/stk_search/CommonSearchUtil.hpp @@ -47,6 +47,8 @@ #include "stk_search/kdtree/KDTree_BoundingBox.hpp" #include "stk_search/kdtree/KDTree.hpp" #include "DeviceMPIUtils.hpp" +#include +#include namespace stk::search { @@ -431,6 +433,11 @@ class SearchResultCommunication DeviceBuffers m_recvBuffers; }; +struct SearchData +{ + std::any data; +}; + } // end namespace stk::search #endif diff --git a/packages/stk/stk_search/stk_search/LocalCoarseSearch.hpp b/packages/stk/stk_search/stk_search/LocalCoarseSearch.hpp index cbc0ba445782..f355788ea527 100644 --- a/packages/stk/stk_search/stk_search/LocalCoarseSearch.hpp +++ b/packages/stk/stk_search/stk_search/LocalCoarseSearch.hpp @@ -53,23 +53,24 @@ void local_coarse_search( std::vector> const & domain, std::vector> const & range, SearchMethod method, - std::vector> & intersections) + std::vector> & intersections, + bool sortSearchResults = false) { switch (method) { case ARBORX: { #ifdef STK_HAS_ARBORX - local_coarse_search_arborx(domain, range, intersections); + local_coarse_search_arborx(domain, range, intersections, sortSearchResults); #else STK_ThrowErrorMsg("STK(stk_search) was not configured with ARBORX. Please use KDTREE or MORTON_LBVH."); #endif break; } case KDTREE: { - local_coarse_search_kdtree_driver(domain, range, intersections); + local_coarse_search_kdtree_driver(domain, range, intersections, sortSearchResults); break; } case MORTON_LBVH: { - local_coarse_search_morton_lbvh(domain, range, intersections); + local_coarse_search_morton_lbvh(domain, range, intersections, sortSearchResults); break; } default: { @@ -81,20 +82,22 @@ void local_coarse_search( template -void local_coarse_search( +std::shared_ptr + local_coarse_search( DomainView const & domain, RangeView const & range, SearchMethod method, ResultView & intersections, - ExecutionSpace const& execSpace = ExecutionSpace{}) - + ExecutionSpace const& execSpace = ExecutionSpace{}, + bool sortSearchResults = false, + std::shared_ptr searchData = nullptr) { check_coarse_search_types_local(); switch (method) { case ARBORX: { #ifdef STK_HAS_ARBORX - local_coarse_search_arborx(domain, range, intersections, execSpace); + local_coarse_search_arborx(domain, range, intersections, execSpace, sortSearchResults); #else STK_ThrowErrorMsg("STK(stk_search) was not configured with ARBORX. Please use KDTREE or MORTON_LBVH."); #endif @@ -105,7 +108,7 @@ void local_coarse_search( break; } case MORTON_LBVH: { - local_coarse_search_morton_lbvh(domain, range, intersections, execSpace); + searchData = local_coarse_search_morton_lbvh(domain, range, intersections, execSpace, sortSearchResults, searchData); break; } default: { @@ -113,6 +116,7 @@ void local_coarse_search( } } + return searchData; } } diff --git a/packages/stk/stk_search/stk_search/Plane.hpp b/packages/stk/stk_search/stk_search/Plane.hpp index f9fa2473462d..1a865a55b84f 100644 --- a/packages/stk/stk_search/stk_search/Plane.hpp +++ b/packages/stk/stk_search/stk_search/Plane.hpp @@ -36,7 +36,8 @@ #define STK_SEARCH_PLANE_HPP #include -#include +#include +#include #include namespace stk { namespace search { @@ -49,7 +50,7 @@ class Plane typedef Point point_type; static const int Dim = 3; - static KOKKOS_FUNCTION constexpr value_type max() { return Kokkos::Details::ArithTraits::max() ;} + static KOKKOS_FUNCTION constexpr value_type max() { return Kokkos::Experimental::finite_max_v ;} KOKKOS_FUNCTION point_type cross(const point_type& a, const point_type& b) const{ return point_type(a[1]*b[2] - a[2]*b[1], @@ -58,7 +59,7 @@ class Plane } KOKKOS_FUNCTION void normalize(point_type& a) const{ - const value_type vec_len = Kokkos::Details::ArithTraits::sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]); + const value_type vec_len = Kokkos::sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]); const value_type denom = vec_len == 0.0 ? 1.0 : vec_len; const value_type vec_len_inv = 1.0 / denom; a[0] *= vec_len_inv; diff --git a/packages/stk/stk_search/stk_search/arborx/CoarseSearchArborX.hpp b/packages/stk/stk_search/stk_search/arborx/CoarseSearchArborX.hpp index e7ad3f271847..7f8ffb299fe5 100644 --- a/packages/stk/stk_search/stk_search/arborx/CoarseSearchArborX.hpp +++ b/packages/stk/stk_search/stk_search/arborx/CoarseSearchArborX.hpp @@ -194,7 +194,8 @@ inline void coarse_search_arborx(std::vector> const& localRange, MPI_Comm comm, std::vector>& searchResults, - bool enforceSearchResultSymmetry = true) + bool enforceSearchResultSymmetry = true, + bool sortSearchResults = false) { using HostSpace = Kokkos::DefaultHostExecutionSpace; using MemSpace = typename HostSpace::memory_space; @@ -237,6 +238,11 @@ inline void coarse_search_arborx(std::vector(); using HostSpace = Kokkos::DefaultHostExecutionSpace; @@ -318,10 +325,12 @@ inline void coarse_search_arborx( Kokkos::resize(Kokkos::WithoutInitializing, searchResults, searchResultsHost.extent(0)); Kokkos::deep_copy(searchResults, searchResultsHost); Kokkos::Profiling::popRegion(); - - Kokkos::Profiling::pushRegion("sorting searchResults"); - Kokkos::sort(searchResults); - Kokkos::Profiling::popRegion(); + + if (sortSearchResults) { + Kokkos::Profiling::pushRegion("Sort searchResults"); + Kokkos::sort(searchResults, Comparator()); + Kokkos::Profiling::popRegion(); + } } diff --git a/packages/stk/stk_search/stk_search/arborx/LocalCoarseSearchArborX.hpp b/packages/stk/stk_search/stk_search/arborx/LocalCoarseSearchArborX.hpp index e95c6f4223b5..c5db6c6b3faf 100644 --- a/packages/stk/stk_search/stk_search/arborx/LocalCoarseSearchArborX.hpp +++ b/packages/stk/stk_search/stk_search/arborx/LocalCoarseSearchArborX.hpp @@ -156,7 +156,8 @@ class SearchResultsInserter template inline void local_coarse_search_arborx(const std::vector>& localDomain, const std::vector>& localRange, - std::vector>& searchResults) + std::vector>& searchResults, + bool sortSearchResults = false) { using ExecSpace = Kokkos::DefaultHostExecutionSpace; using MemSpace = typename ExecSpace::memory_space; @@ -196,6 +197,12 @@ inline void local_coarse_search_arborx(const std::vector bvh(execSpace, localRangeViewWrapped); bvh.query(execSpace, localDomainViewWrapped, callback); execSpace.fence(); + + if (sortSearchResults) { + Kokkos::Profiling::pushRegion("Sort searchResults"); + std::sort(searchResults.begin(), searchResults.end()); + Kokkos::Profiling::popRegion(); + } } template (); using MemSpace = typename ExecutionSpace::memory_space; @@ -258,6 +266,13 @@ inline void local_coarse_search_arborx( Kokkos::fence(); Kokkos::Profiling::popRegion(); + + if (sortSearchResults) { + Kokkos::Profiling::pushRegion("Sort searchResults"); + Kokkos::sort(searchResults, Comparator()); + Kokkos::Profiling::popRegion(); + } + } } // namespace stk::search diff --git a/packages/stk/stk_search/stk_search/kdtree/CoarseSearchKdTree.hpp b/packages/stk/stk_search/stk_search/kdtree/CoarseSearchKdTree.hpp index 95a13b5b70cb..936478f1bcfb 100644 --- a/packages/stk/stk_search/stk_search/kdtree/CoarseSearchKdTree.hpp +++ b/packages/stk/stk_search/stk_search/kdtree/CoarseSearchKdTree.hpp @@ -58,7 +58,8 @@ inline void coarse_search_kdtree(std::vector< std::pair > const & local_range, MPI_Comm comm, std::vector >& searchResults, - bool enforceSearchResultSymmetry=true) + bool enforceSearchResultSymmetry = true, + bool sortSearchResults = false) { int num_procs = -1; @@ -146,6 +147,9 @@ inline void coarse_search_kdtree(std::vector< std::pair, RangeIdentifier > > const & local_range, MPI_Comm comm, std::vector >& searchResults, - bool enforceSearchResultSymmetry=true) + bool enforceSearchResultSymmetry=true, + bool sortSearchResults = false) { int num_procs = -1; int proc_id = -1; @@ -232,6 +237,9 @@ inline void coarse_search_kdtree(std::vector< std::pair > const & local_range, MPI_Comm comm, std::vector >& searchResults, - bool enforceSearchResultSymmetry=true) + bool enforceSearchResultSymmetry = true, + bool sortSearchResults = false) { const size_t local_sizes[2] = {local_domain.size(), local_range.size()}; size_t global_sizes[2]; @@ -249,12 +258,12 @@ inline void coarse_search_kdtree_driver(std::vector< std::pair= global_sizes[1]); if(domain_has_more_boxes) { - coarse_search_kdtree(local_domain, local_range, comm, searchResults, enforceSearchResultSymmetry); + coarse_search_kdtree(local_domain, local_range, comm, searchResults, enforceSearchResultSymmetry, sortSearchResults); } else { std::vector > tempSearchResults; - coarse_search_kdtree(local_range, local_domain, comm, tempSearchResults, enforceSearchResultSymmetry); + coarse_search_kdtree(local_range, local_domain, comm, tempSearchResults, enforceSearchResultSymmetry, sortSearchResults); const int p_rank = stk::parallel_machine_rank(comm); searchResults.reserve(tempSearchResults.size()); for(size_t i=0; i()(pair.second) == p_rank) searchResults.emplace_back(pair.second, pair.first); } + + if (sortSearchResults) { + std::sort(searchResults.begin(), searchResults.end()); + } } } diff --git a/packages/stk/stk_search/stk_search/kdtree/LocalCoarseSearchKdTree.hpp b/packages/stk/stk_search/stk_search/kdtree/LocalCoarseSearchKdTree.hpp index 3541f8f05d82..dc3f37e5ef6a 100644 --- a/packages/stk/stk_search/stk_search/kdtree/LocalCoarseSearchKdTree.hpp +++ b/packages/stk/stk_search/stk_search/kdtree/LocalCoarseSearchKdTree.hpp @@ -50,7 +50,8 @@ template > const & local_domain, std::vector< std::pair > const & local_range, std::vector >& searchResults, - bool enforceSearchResultSymmetry=true) + bool enforceSearchResultSymmetry = true, + bool sortSearchResults = false) { #ifdef _OPENMP @@ -116,6 +117,10 @@ inline void local_coarse_search_kdtree(std::vector< std::pair > const & local_domain, std::vector< std::pair, RangeIdentifier > > const & local_range, std::vector >& searchResults, - bool enforceSearchResultSymmetry=true) + bool enforceSearchResultSymmetry = true, + bool sortSearchResults = false) { searchResults.clear(); @@ -186,28 +192,37 @@ inline void local_coarse_search_kdtree(std::vector< std::pair inline void local_coarse_search_kdtree_driver(std::vector< std::pair > const & local_domain, std::vector< std::pair > const & local_range, - std::vector >& searchResults) + std::vector >& searchResults, + bool sortSearchResults = false) { const bool domain_has_more_boxes = local_domain.size() > local_range.size(); if(domain_has_more_boxes) { - local_coarse_search_kdtree(local_domain, local_range, searchResults); + local_coarse_search_kdtree(local_domain, local_range, searchResults, sortSearchResults); } else { std::vector > tempSearchResults; - local_coarse_search_kdtree(local_range, local_domain, tempSearchResults); + local_coarse_search_kdtree(local_range, local_domain, tempSearchResults, sortSearchResults); searchResults.reserve(tempSearchResults.size()); for (auto& [range_id, domain_id] : tempSearchResults) { searchResults.emplace_back(domain_id, range_id); } + + if (sortSearchResults) { + std::sort(searchResults.begin(), searchResults.end()); + } } } diff --git a/packages/stk/stk_search/stk_search/morton_lbvh/CoarseSearchMortonLBVH.hpp b/packages/stk/stk_search/stk_search/morton_lbvh/CoarseSearchMortonLBVH.hpp index 4cc6a3b6a6a0..36094b37baa0 100644 --- a/packages/stk/stk_search/stk_search/morton_lbvh/CoarseSearchMortonLBVH.hpp +++ b/packages/stk/stk_search/stk_search/morton_lbvh/CoarseSearchMortonLBVH.hpp @@ -133,7 +133,8 @@ inline void coarse_search_morton_lbvh(std::vector> const & localRange, MPI_Comm comm, std::vector> & searchResults, - bool enforceSearchResultSymmetry = true) + bool enforceSearchResultSymmetry = true, + bool sortSearchResults = false) { using HostSpace = Kokkos::DefaultHostExecutionSpace; using Callback = impl::MortonCoarseSearchVectorCallback; @@ -141,31 +142,44 @@ inline void coarse_search_morton_lbvh(std::vector), "The domain and range boxes must have the same floating-point precision"); - using ValueType = typename DomainBoxType::value_type; + using DomainValueType = typename DomainBoxType::value_type; + using RangeValueType = typename RangeBoxType::value_type; Kokkos::Profiling::pushRegion("Parallel consistency: extend range box list"); const auto [extendedRangeBoxes, remoteRangeIdentProcs] = morton_extend_local_range_with_remote_boxes_that_might_intersect(localDomain, localRange, comm, HostSpace{}); Kokkos::Profiling::popRegion(); - Kokkos::Profiling::pushRegion("Fill domain and range trees"); - stk::search::MortonAabbTree domainTree("Domain Tree", localDomain.size()); - stk::search::MortonAabbTree rangeTree("Range Tree", extendedRangeBoxes.size()); + using StkDomainBoxType = stk::search::Box; + using StkRangeBoxType = stk::search::Box; - stk::search::export_from_box_ident_proc_vec_to_morton_tree(localDomain, domainTree); - stk::search::export_from_box_vec_to_morton_tree(extendedRangeBoxes, rangeTree); + Kokkos::Profiling::pushRegion("Fill domain and range trees"); + using DomainViewType = Kokkos::View*,HostSpace>; + using RangeViewType = Kokkos::View*,HostSpace>; + using DomainTreeType = stk::search::MortonAabbTree; + using RangeTreeType = stk::search::MortonAabbTree; + DomainTreeType domainTree("Domain Tree", localDomain.size()); + RangeTreeType rangeTree("Range Tree", extendedRangeBoxes.size()); + + stk::search::export_from_box_ident_proc_vec_to_morton_tree(localDomain, domainTree); + stk::search::export_from_box_vec_to_morton_tree(extendedRangeBoxes, rangeTree); domainTree.sync_to_device(); rangeTree.sync_to_device(); Kokkos::Profiling::popRegion(); Kokkos::Profiling::pushRegion("Perform Morton query"); Callback callback(localDomain, localRange, extendedRangeBoxes, remoteRangeIdentProcs, searchResults); - stk::search::morton_lbvh_search(domainTree, rangeTree, callback, HostSpace{}); + stk::search::morton_lbvh_search(domainTree, rangeTree, callback, HostSpace{}); Kokkos::Profiling::popRegion(); if (enforceSearchResultSymmetry) { Kokkos::Profiling::pushRegion("Enforce results symmetry"); stk::search::communicate_vector(comm, searchResults, enforceSearchResultSymmetry); + Kokkos::Profiling::popRegion(); + } + + if (sortSearchResults) { + Kokkos::Profiling::pushRegion("Sort searchResults"); std::sort(searchResults.begin(), searchResults.end()); Kokkos::Profiling::popRegion(); } @@ -263,7 +277,6 @@ class BoundingShapeIntersectionCheckFunctor } - template (); @@ -291,6 +305,9 @@ inline void coarse_search_morton_lbvh( static_assert(std::is_same_v, "The domain and range boxes must have the same floating-point precision"); + using MDomainViewType = typename BoxIdentProcViewTrait::ViewType; + using MRangeViewType = typename BoxIdentProcViewTrait::ViewType; + Kokkos::Profiling::pushRegion("Parallel consistency: extend range box list"); ExtendedRangeBoxView extendedRangeBoxes; ExtendedRangeIdentProcView remoteRangeIdentProcs; @@ -318,13 +335,21 @@ inline void coarse_search_morton_lbvh( } Kokkos::Profiling::popRegion(); - Kokkos::Profiling::pushRegion("Fill domain and range trees"); + Kokkos::Profiling::pushRegion("STK Fill domain and range trees"); - bool setBoxesOnHost = false; - stk::search::MortonAabbTree domainTree("Domain Tree", localDomain.extent(0), setBoxesOnHost); - stk::search::MortonAabbTree rangeTree("Range Tree", extendedRangeBoxes.extent(0), setBoxesOnHost); + const bool setBoxesOnHost = false; + using DomainTreeType = stk::search::MortonAabbTree; + using RangeTreeType = stk::search::MortonAabbTree; + DomainTreeType domainTree("Domain Tree", localDomain.extent(0), setBoxesOnHost); + RangeTreeType rangeTree("Range Tree", extendedRangeBoxes.extent(0), setBoxesOnHost); + + if constexpr (std::is_same_v) { + domainTree.m_minMaxs = localDomain; + } + else { + stk::search::export_box_ident_view_to_morton_tree(localDomain, domainTree, execSpace); + } - stk::search::export_box_ident_view_to_morton_tree(localDomain, domainTree, execSpace); stk::search::export_box_view_to_morton_tree(extendedRangeBoxes, rangeTree, execSpace); execSpace.fence(); @@ -341,21 +366,24 @@ inline void coarse_search_morton_lbvh( Kokkos::Profiling::pushRegion("Inner morton search"); BoundingShapeIntersectionChecker intersectionChecker(localDomain, localRange, extendedRangeBoxes, remoteRangeIdentProcs, searchResults); - stk::search::morton_lbvh_search(domainTree, rangeTree, intersectionChecker); + stk::search::morton_lbvh_search(domainTree, rangeTree, intersectionChecker, execSpace); searchResults = intersectionChecker.get_search_results(); Kokkos::Profiling::popRegion(); - Kokkos::Profiling::pushRegion("Kokkos sort"); - Kokkos::sort(searchResults); - Kokkos::Profiling::popRegion(); - if (enforceSearchResultSymmetry) { Kokkos::Profiling::pushRegion("Enforce results symmetry"); SearchResultCommunication resultComm(comm, searchResults, execSpace, enforceSearchResultSymmetry); searchResults = resultComm.run(); Kokkos::Profiling::popRegion(); } + + if (sortSearchResults) { + Kokkos::Profiling::pushRegion("Sort searchResults"); + Kokkos::sort(searchResults, Comparator()); + Kokkos::Profiling::popRegion(); + } + } } diff --git a/packages/stk/stk_search/stk_search/morton_lbvh/LocalCoarseSearchMortonLBVH.hpp b/packages/stk/stk_search/stk_search/morton_lbvh/LocalCoarseSearchMortonLBVH.hpp index ef1f21f0ca56..1d6b39beab82 100644 --- a/packages/stk/stk_search/stk_search/morton_lbvh/LocalCoarseSearchMortonLBVH.hpp +++ b/packages/stk/stk_search/stk_search/morton_lbvh/LocalCoarseSearchMortonLBVH.hpp @@ -40,11 +40,13 @@ #include "stk_search/HelperTraits.hpp" #include "stk_search/morton_lbvh/MortonLBVH_Search.hpp" #include "stk_search/morton_lbvh/MortonLBVH_Tree.hpp" +#include "stk_util/ngp/NgpSpaces.hpp" #include #include #include "Kokkos_Core.hpp" +#include "Kokkos_Sort.hpp" namespace stk::search { @@ -88,6 +90,8 @@ class LocalMortonCoarseSearchVectorCallback #endif } + ResultVec& get_search_results() const { return m_searchResults; } + bool resize_for_second_pass() { return false; @@ -105,7 +109,8 @@ template > & domain, const std::vector> & range, - std::vector> & searchResults) + std::vector> & searchResults, + bool sortSearchResults = false) { STK_ThrowRequireMsg((std::is_same_v), "The domain and range boxes must have the same floating-point precision"); @@ -116,23 +121,36 @@ void local_coarse_search_morton_lbvh( Kokkos::Profiling::pushRegion("Fill domain and range trees"); const bool supportHostBoxes = false; - stk::search::MortonAabbTree domainTree("Domain Tree", domain.size(), supportHostBoxes); - stk::search::MortonAabbTree rangeTree("Range Tree", range.size(), supportHostBoxes); - - stk::search::export_from_box_ident_vector_to_morton_tree(domain, domainTree); - stk::search::export_from_box_ident_vector_to_morton_tree(range, rangeTree); + using MDomainBoxType = stk::search::Box; + using MRangeBoxType = stk::search::Box; + using DomainViewType = Kokkos::View*, HostSpace>; + using RangeViewType = Kokkos::View*, HostSpace>; + using DomainTreeType = stk::search::MortonAabbTree; + using RangeTreeType = stk::search::MortonAabbTree; + + DomainTreeType domainTree("Domain Tree", domain.size(), supportHostBoxes); + RangeTreeType rangeTree("Range Tree", range.size(), supportHostBoxes); + + stk::search::export_from_box_ident_vector_to_morton_tree(domain, domainTree); + stk::search::export_from_box_ident_vector_to_morton_tree(range, rangeTree); Kokkos::Profiling::popRegion(); stk::search::CollisionList collisionList("Collision List"); Callback callback(domain, range, searchResults); - stk::search::morton_lbvh_search(domainTree, rangeTree, callback, HostSpace{}); + stk::search::morton_lbvh_search(domainTree, rangeTree, callback, HostSpace{}); + searchResults = callback.get_search_results(); + + if (sortSearchResults) { + Kokkos::Profiling::pushRegion("Sort searchResults"); + std::sort(searchResults.begin(), searchResults.end()); + Kokkos::Profiling::popRegion(); + } } namespace impl { template class LocalMortonCoarseSearchViewCallback { - using DomainBoxIdent = typename DomainView::value_type; using RangeBoxIdent = typename RangeView::value_type; using DomainBox = typename DomainBoxIdent::box_type; @@ -140,6 +158,8 @@ class LocalMortonCoarseSearchViewCallback static bool constexpr isSearchExact = !(impl::is_stk_sphere || impl::is_stk_sphere); public: + using ResultViewType = ResultView; + LocalMortonCoarseSearchViewCallback(const DomainView& domain, const RangeView& range, ResultView& searchResults) : m_domain(domain), m_range(range), @@ -150,6 +170,14 @@ class LocalMortonCoarseSearchViewCallback Kokkos::deep_copy(m_idx, 0); } + void reset(const DomainView& domain, const RangeView& range, ResultView& resultsView) + { + m_domain = domain; + m_range = range; + m_searchResults = resultsView; + Kokkos::deep_copy(ExecutionSpace{}, m_idx, 0); + } + KOKKOS_INLINE_FUNCTION void operator()(int domainIdx, int rangeIdx) const { @@ -178,10 +206,12 @@ class LocalMortonCoarseSearchViewCallback bool resize_for_second_pass() { unsigned int numResults = 0; - Kokkos::deep_copy(numResults, m_idx); - bool needSecondPass = numResults > m_searchResults.size(); + Kokkos::deep_copy(ExecutionSpace{}, numResults, m_idx); + const bool needSecondPass = numResults > m_searchResults.size(); Kokkos::resize(Kokkos::WithoutInitializing, m_searchResults, numResults); - Kokkos::deep_copy(m_idx, 0); + if (needSecondPass) { + Kokkos::deep_copy(m_idx, 0); + } return needSecondPass; } @@ -197,16 +227,37 @@ class LocalMortonCoarseSearchViewCallback } +template +struct MortonData +{ + using ResultViewType = typename CallbackType::ResultViewType; + template + MortonData(const DomainViewType& domain, const RangeViewType& range, bool supportHostBoxes, ResultViewType& results) + : domainTree("Domain Tree", domain.extent(0), supportHostBoxes), + rangeTree("Range Tree", range.extent(0), supportHostBoxes), + callback(domain, range, results) + {} + + DomainTreeType domainTree; + RangeTreeType rangeTree; + CallbackType callback; +}; + template -void local_coarse_search_morton_lbvh( +std::shared_ptr +local_coarse_search_morton_lbvh( const DomainView & domain, const RangeView & range, ResultView & searchResults, - ExecutionSpace const& execSpace = ExecutionSpace{}) + ExecutionSpace const& execSpace = ExecutionSpace{}, + bool sortSearchResults = false, + std::shared_ptr searchData = nullptr) { Kokkos::Profiling::pushRegion("local_coarse_search_morton_lbvh"); check_coarse_search_types_local(); + using DomainIdentType = typename DomainView::value_type::ident_type; + using RangeIdentType = typename RangeView::value_type::ident_type; using DomainBoxType = typename DomainView::value_type::box_type; using RangeBoxType = typename RangeView::value_type::box_type; STK_ThrowRequireMsg((std::is_same_v), @@ -216,28 +267,87 @@ void local_coarse_search_morton_lbvh( Kokkos::Profiling::pushRegion("STK Fill domain and range trees"); const bool supportHostBoxes = false; - stk::search::MortonAabbTree domainTree("Domain Tree", domain.extent(0), supportHostBoxes); - stk::search::MortonAabbTree rangeTree("Range Tree", range.extent(0), supportHostBoxes); - Kokkos::Profiling::pushRegion("STK Export box ident views to trees"); - stk::search::export_box_ident_view_to_morton_tree(domain, domainTree, execSpace); - stk::search::export_box_ident_view_to_morton_tree(range, rangeTree, execSpace); + using MDomainViewType = typename BoxIdentViewTrait::ViewType; + using MRangeViewType = typename BoxIdentViewTrait::ViewType; + + using DomainTreeType = stk::search::MortonAabbTree; + using RangeTreeType = stk::search::MortonAabbTree; + using MortonDataType = MortonData; + + bool newlyCreatedData = false; + if (searchData == nullptr || !searchData->data.has_value()) { + searchData = std::make_shared(); + searchData->data = std::make_shared(domain, range, supportHostBoxes, searchResults); + newlyCreatedData = true; + } + + std::shared_ptr mortonData; + try { + mortonData = std::any_cast>(searchData->data); + if (!newlyCreatedData) { + mortonData->domainTree.reset(domain.extent(0)); + mortonData->rangeTree.reset(range.extent(0)); + Kokkos::Profiling::pushRegion("Callback.reset"); + mortonData->callback.reset(domain, range, searchResults); + Kokkos::Profiling::popRegion(); + } + } + catch (const std::bad_any_cast& e) { + searchData->data = std::make_shared(domain, range, supportHostBoxes, searchResults); + } + + DomainTreeType& domainTree = mortonData->domainTree; + RangeTreeType& rangeTree = mortonData->rangeTree; + + + if constexpr (std::is_same_v) { + domainTree.m_minMaxs = domain; + } + else { + Kokkos::Profiling::pushRegion("STK Export box ident views to trees"); + stk::search::export_box_ident_view_to_morton_tree(domain, domainTree, execSpace); + execSpace.fence(); + Kokkos::Profiling::popRegion(); + } + if constexpr (std::is_same_v) { + rangeTree.m_minMaxs = range; + } + else { + Kokkos::Profiling::pushRegion("STK Export box ident views to trees"); + stk::search::export_box_ident_view_to_morton_tree(range, rangeTree, execSpace); + execSpace.fence(); + Kokkos::Profiling::popRegion(); + } execSpace.fence(); Kokkos::Profiling::popRegion(); - Kokkos::Profiling::popRegion(); Kokkos::Profiling::pushRegion("STK Morton Search"); + Kokkos::Profiling::pushRegion("Pre-resize results view"); if (searchResults.size() == 0) { size_t sizeGuess = std::max(domain.size(), range.size()) * COLLISION_SCALE_FACTOR; Kokkos::resize(Kokkos::WithoutInitializing, searchResults, sizeGuess); } + Kokkos::Profiling::popRegion(); - Callback callback(domain, range, searchResults); - stk::search::morton_lbvh_search(domainTree, rangeTree, callback, execSpace); + Kokkos::Profiling::pushRegion("Assign callback reference"); + Callback& callback = mortonData->callback; + Kokkos::Profiling::popRegion(); + + stk::search::morton_lbvh_search(domainTree, rangeTree, callback, execSpace); searchResults = callback.get_search_results(); Kokkos::Profiling::popRegion(); + + if (sortSearchResults) { + Kokkos::Profiling::pushRegion("Sort searchResults"); + Kokkos::sort(searchResults, Comparator()); + Kokkos::Profiling::popRegion(); + } + Kokkos::Profiling::popRegion(); + + return searchData; } } diff --git a/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_ParallelConsistencyUtils.hpp b/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_ParallelConsistencyUtils.hpp index d8c8e297c484..a228a4de1080 100644 --- a/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_ParallelConsistencyUtils.hpp +++ b/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_ParallelConsistencyUtils.hpp @@ -132,24 +132,30 @@ morton_extend_local_range_with_remote_boxes_that_might_intersect( ExecutionSpace const& execSpace) { using DomainValueType = typename DomainBoxType::value_type; + using RangeValueType = typename RangeBoxType::value_type; const int numProcs = stk::parallel_machine_size(comm); const int procId = stk::parallel_machine_rank(comm); - const auto globalSupersetBoxes = gather_all_processor_superset_domain_boxes(localDomain, comm); + using StkDomainBoxType = stk::search::Box; + using StkRangeBoxType = stk::search::Box; + const std::vector globalSupersetBoxes = gather_all_processor_superset_domain_boxes(localDomain, comm); - stk::search::MortonAabbTree domainTree("Proc Domain Tree", - localRange.size()); - stk::search::MortonAabbTree rangeTree("Proc Range Tree", - globalSupersetBoxes.size()); + using DomainViewType = Kokkos::View*,ExecutionSpace>; + using RangeViewType = Kokkos::View*,ExecutionSpace>; + using DomainTreeType = stk::search::MortonAabbTree; + using RangeTreeType = stk::search::MortonAabbTree; - export_from_box_ident_proc_vec_to_morton_tree(localRange, domainTree); - export_from_box_vec_to_morton_tree(globalSupersetBoxes, rangeTree); + DomainTreeType domainTree("Proc Domain Tree", localRange.size()); + RangeTreeType rangeTree("Proc Range Tree", globalSupersetBoxes.size()); + + export_from_box_ident_proc_vec_to_morton_tree(localRange, domainTree); + export_from_box_vec_to_morton_tree(globalSupersetBoxes, rangeTree); domainTree.sync_to_device(); rangeTree.sync_to_device(); stk::search::CollisionList collisionList("Proc Collision List"); - stk::search::morton_lbvh_search(domainTree, rangeTree, collisionList, execSpace); + stk::search::morton_lbvh_search(domainTree, rangeTree, collisionList, execSpace); collisionList.sync_from_device(); using GlobalIdType = typename RangeIdentProcType::ident_type; @@ -260,11 +266,17 @@ morton_extend_local_range_with_remote_boxes_that_might_intersect( { check_domain_or_range_view_parallel(); check_domain_or_range_view_parallel(); - using Real = typename DomainView::value_type::box_type::value_type; + using DomainRealType = typename DomainView::value_type::box_type::value_type; + using RangeRealType = typename RangeView::value_type::box_type::value_type; using RangeBoxType = typename RangeView::value_type::box_type; using RangeIdentProcType = typename RangeView::value_type::ident_proc_type; - using ViewType = Kokkos::View*, ExecutionSpace>; + using ViewType = Kokkos::View*, ExecutionSpace>; + using BoxIdentViewType = Kokkos::View,int>*,ExecutionSpace>; + using MRangeBoxType = stk::search::Box; + using MRangeViewType = Kokkos::View*,ExecutionSpace>; + using DomainTreeType = stk::search::MortonAabbTree; + using RangeTreeType = stk::search::MortonAabbTree; using FillMPIBuffersType = impl::FillGhostBoxBuffers; using DeviceBuffers = typename FillMPIBuffersType::DeviceBuffers; @@ -272,19 +284,17 @@ morton_extend_local_range_with_remote_boxes_that_might_intersect( ViewType globalSupersetBoxes = gather_all_processor_superset_domain_boxes(localDomain, execSpace, comm); const bool setBoxesOnHost = false; - stk::search::MortonAabbTree domainTree("Proc Domain Tree", - localRange.size(), setBoxesOnHost); - stk::search::MortonAabbTree rangeTree("Proc Range Tree", - globalSupersetBoxes.size(), setBoxesOnHost); + DomainTreeType domainTree("Proc Domain Tree", localRange.size(), setBoxesOnHost); + RangeTreeType rangeTree("Proc Range Tree", globalSupersetBoxes.size(), setBoxesOnHost); - export_box_ident_view_to_morton_tree(localRange, domainTree, execSpace); - export_box_view_to_morton_tree(globalSupersetBoxes, rangeTree, execSpace); + export_box_ident_view_to_morton_tree(localRange, domainTree, execSpace); + export_box_view_to_morton_tree(globalSupersetBoxes, rangeTree, execSpace); execSpace.fence(); domainTree.sync_to_device(); rangeTree.sync_to_device(); stk::search::CollisionList collisionList("Proc Collision List"); - stk::search::morton_lbvh_search(domainTree, rangeTree, collisionList); + stk::search::morton_lbvh_search(domainTree, rangeTree, collisionList, execSpace); FillMPIBuffersType fill_buffers(collisionList, localRange, execSpace, comm); DeviceBuffers deviceSendBuffers = fill_buffers.run(); diff --git a/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_Search.hpp b/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_Search.hpp index 75d98a2ffb72..07dd1bfd57be 100644 --- a/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_Search.hpp +++ b/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_Search.hpp @@ -49,9 +49,9 @@ namespace stk::search { constexpr size_t COLLISION_SCALE_FACTOR = 16; -template -inline void determine_mas_calling(const MortonAabbTree &partialTree1, - const MortonAabbTree &partialTree2, +template +inline void determine_mas_calling(const Tree1Type &partialTree1, + const Tree2Type &partialTree2, bool &sortTree1, bool &sortTree2, bool &flipTrees) { if (ALWAYS_SORT_MORTON_TREES) { @@ -99,10 +99,10 @@ inline void determine_mas_calling(const MortonAabbTree } } -template +template inline void export_from_box_ident_proc_vec_to_morton_tree( const std::vector> &boxIdentProcVec, - MortonAabbTree &tree) + TreeType &tree) { int numBoxes = static_cast(boxIdentProcVec.size()); tree.reset(numBoxes); @@ -114,33 +114,32 @@ inline void export_from_box_ident_proc_vec_to_morton_tree( } } -template +template inline void export_box_ident_view_to_morton_tree( - const DomainView& boxIdentProcs, - MortonAabbTree& tree, + const ViewType& boxIdentProcs, + TreeType& tree, ExecutionSpace execSpace) { - check_view_is_usable_from(); - static_assert(is_box_ident_proc_container_v || is_box_ident_container_v, + check_view_is_usable_from(); + static_assert(is_box_ident_proc_container_v || is_box_ident_container_v, "view must be a View or View"); - using BoxType = typename DomainView::value_type::box_type; + using BoxType = typename ViewType::value_type::box_type; tree.reset(boxIdentProcs.extent(0)); Kokkos::RangePolicy policy(execSpace, 0, boxIdentProcs.extent(0)); auto func = KOKKOS_LAMBDA(int index) { const BoxType box = boxIdentProcs(index).box; - tree.device_set_box(index, box.get_x_min(), box.get_x_max(), - box.get_y_min(), box.get_y_max(), - box.get_z_min(), box.get_z_max()); + tree.device_set_box(index, box.get_x_min(), box.get_y_min(), box.get_z_min(), + box.get_x_max(), box.get_y_max(), box.get_z_max()); }; Kokkos::parallel_for("export box-ident view to tree", policy, func); } -template +template inline void export_from_box_ident_vector_to_morton_tree( - const std::vector> &boxIdentList, MortonAabbTree &tree) + const std::vector> &boxIdentList, TreeType &tree) { static_assert(Kokkos::SpaceAccessibility::assignable); int numBoxes = static_cast(boxIdentList.size()); @@ -149,15 +148,15 @@ inline void export_from_box_ident_vector_to_morton_tree( Kokkos::parallel_for( Kokkos::RangePolicy(0, numBoxes), KOKKOS_LAMBDA(int index) { const auto &box = boxIdentList[index].first; - tree.device_set_box(index, box.get_x_min(), box.get_x_max(), box.get_y_min(), box.get_y_max(), box.get_z_min(), - box.get_z_max()); + tree.device_set_box(index, box.get_x_min(), box.get_y_min(), box.get_z_min(), + box.get_x_max(), box.get_y_max(), box.get_z_max()); }); } -template +template inline void export_from_box_vec_to_morton_tree(const std::vector &boxVec, - MortonAabbTree &tree) + TreeType &tree) { int numBoxes = static_cast(boxVec.size()); tree.reset(numBoxes); @@ -168,9 +167,9 @@ inline void export_from_box_vec_to_morton_tree(const std::vector &boxVe } } -template +template inline void export_box_view_to_morton_tree(const BoxView boxes, - MortonAabbTree& tree, + TreeType& tree, ExecutionSpace execSpace) { check_view_is_usable_from(); @@ -181,62 +180,66 @@ inline void export_box_view_to_morton_tree(const BoxView boxes, auto func = KOKKOS_LAMBDA (int index) { const BoxType box = boxes(index); - tree.device_set_box(index, box.get_x_min(), box.get_x_max(), - box.get_y_min(), box.get_y_max(), - box.get_z_min(), box.get_z_max()); + tree.device_set_box(index, box.get_x_min(), box.get_y_min(), box.get_z_min(), + box.get_x_max(), box.get_y_max(), box.get_z_max()); }; Kokkos::parallel_for(policy, func); } -template -inline void morton_lbvh_search(MortonAabbTree &tree1, - MortonAabbTree &tree2, +template +inline void morton_lbvh_search(MortonAabbTree &tree1, + MortonAabbTree &tree2, Callback& resultCallback, ExecutionSpace const& execSpace = ExecutionSpace{}) { Kokkos::Profiling::pushRegion("Initialize the trees"); - Kokkos::Profiling::pushRegion("Get global bounds"); - // Get total bounds - TotalBoundsFunctor::apply(tree1, execSpace); - TotalBoundsFunctor::apply(tree2, execSpace); - execSpace.fence(); - Kokkos::Profiling::popRegion(); - Kokkos::Profiling::pushRegion("Determine need to sort/flip trees"); bool sortTree1, sortTree2, flipOrder; - determine_mas_calling(tree1, tree2, sortTree1, sortTree2, flipOrder); + determine_mas_calling, + MortonAabbTree, + ExecutionSpace>(tree1, tree2, sortTree1, sortTree2, flipOrder); + Kokkos::Profiling::popRegion(); + + Kokkos::Profiling::pushRegion("Get total bounds"); + TotalBoundsFunctor::apply(tree1, execSpace); + TotalBoundsFunctor::apply(tree2, execSpace); Kokkos::Profiling::popRegion(); - // Morton encode the centroids of the leaves Kokkos::Profiling::pushRegion("Morton encoding of leaves"); - MortonEncoder::apply(tree1, execSpace, sortTree1); - MortonEncoder::apply(tree2, execSpace, sortTree2); - execSpace.fence(); + MortonEncoder::apply(tree1, execSpace, sortTree1); + MortonEncoder::apply(tree2, execSpace, sortTree2); Kokkos::Profiling::popRegion(); - // Sort the leaves if appropriate - Kokkos::Profiling::pushRegion("Sort the trees"); + execSpace.fence(); + if (sortTree1) { - // printf("Sorting tree with %d leaves\n", tree1.hm_numLeaves()); - SortByCode::apply(tree1, execSpace); + Kokkos::Profiling::pushRegion("Sort the trees"); + SortByCode, ExecutionSpace>::apply(tree1, execSpace); + Kokkos::Profiling::popRegion(); } + if (sortTree2) { - // printf("Sorting tree with %d leaves\n", tree2.hm_numLeaves()); - SortByCode::apply(tree2, execSpace); + Kokkos::Profiling::pushRegion("Sort the trees"); + SortByCode, ExecutionSpace>::apply(tree2, execSpace); + Kokkos::Profiling::popRegion(); } + execSpace.fence(); - Kokkos::Profiling::popRegion(); // Build the tree structures, if appropriate, following Karras's algorithm Kokkos::Profiling::pushRegion("Build the tree structures"); bool buildTree1 = (sortTree1 && flipOrder); bool buildTree2 = (sortTree2 && !flipOrder); + using Box1Type = typename View1Type::value_type::box_type; + using Real1Type = typename Box1Type::value_type; + using Box2Type = typename View2Type::value_type::box_type; + using Real2Type = typename Box2Type::value_type; if (buildTree1) { - BuildRadixTree::apply(tree1, execSpace); + BuildRadixTree, ExecutionSpace>::apply(tree1, execSpace); } if (buildTree2) { - BuildRadixTree::apply(tree2, execSpace); + BuildRadixTree, ExecutionSpace>::apply(tree2, execSpace); } execSpace.fence(); Kokkos::Profiling::popRegion(); @@ -244,10 +247,10 @@ inline void morton_lbvh_search(MortonAabbTree &tree1, // Augment the trees to be bounding volume (box) hierarchies Kokkos::Profiling::pushRegion("Augment the trees to be bounding volume hierarchies"); if (buildTree1) { - UpdateInteriorNodeBVs::apply(tree1, execSpace); + UpdateInteriorNodeBVs::apply(tree1, execSpace); } if (buildTree2) { - UpdateInteriorNodeBVs::apply(tree2, execSpace); + UpdateInteriorNodeBVs::apply(tree2, execSpace); } execSpace.fence(); Kokkos::Profiling::popRegion(); @@ -265,9 +268,9 @@ inline void morton_lbvh_search(MortonAabbTree &tree1, Kokkos::Profiling::popRegion(); } -template -inline void morton_lbvh_search(MortonAabbTree &tree1, - MortonAabbTree &tree2, +template +inline void morton_lbvh_search(Tree1Type &tree1, + Tree2Type &tree2, CollisionList &searchResults, ExecutionSpace const& execSpace = ExecutionSpace{}) { @@ -291,15 +294,17 @@ inline void morton_lbvh_search(const std::vector &boxA, ExecutionSpace const& execSpace = ExecutionSpace{}) { Kokkos::Profiling::pushRegion("morton_lbvh_search: export boxes to trees"); - MortonAabbTree mlbvhA("a"), mlbvhB("b"); - export_from_box_vec_to_morton_tree(boxA, mlbvhA); - export_from_box_vec_to_morton_tree(boxB, mlbvhB); + using ViewType = Kokkos::View*,ExecutionSpace>; + using TreeType = MortonAabbTree; + TreeType mlbvhA("a"), mlbvhB("b"); + export_from_box_vec_to_morton_tree(boxA, mlbvhA); + export_from_box_vec_to_morton_tree(boxB, mlbvhB); mlbvhA.sync_to_device(); mlbvhB.sync_to_device(); Kokkos::Profiling::popRegion(); Kokkos::Profiling::pushRegion("morton_lbvh_search: execute search"); - morton_lbvh_search(mlbvhA, mlbvhB, searchResults, execSpace); + morton_lbvh_search(mlbvhA, mlbvhB, searchResults, execSpace); Kokkos::Profiling::popRegion(); } diff --git a/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_Tree.hpp b/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_Tree.hpp index 1dcf86cca14e..b2f5f0a90c65 100644 --- a/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_Tree.hpp +++ b/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_Tree.hpp @@ -36,6 +36,8 @@ #define STK_SEARCH_MORTON_LBVH_MORTONLBVH_TREE_HPP #include +#include +#include #include #include #include @@ -46,11 +48,51 @@ namespace stk::search { -template +template +KOKKOS_INLINE_FUNCTION +const BoxType& get_box(BoxType& box) { return box; } + +template +KOKKOS_INLINE_FUNCTION +const BoxType& get_box(BoxIdent& boxIdent) { return boxIdent.box; } + +template +KOKKOS_INLINE_FUNCTION +const BoxType& get_box(BoxIdentProc& boxIdentProc) { return boxIdentProc.box; } + + +template struct BoxIdentViewTrait {}; + +template +struct BoxIdentViewTrait> +{ + using InputBoxType = typename ViewValueType::box_type; + using ValueType = typename InputBoxType::value_type; + using IdentType = typename ViewValueType::ident_type; + using ViewType = Kokkos::View,IdentType>*, otherTemplateArgs...>; +}; + +template struct BoxIdentProcViewTrait {}; + +template +struct BoxIdentProcViewTrait> +{ + using InputBoxType = typename ViewValueType::box_type; + using ValueType = typename InputBoxType::value_type; + using IdentProcType = typename ViewValueType::ident_proc_type; + using ViewType = Kokkos::View,IdentProcType>*, otherTemplateArgs...>; +}; + + +template struct MortonAabbTree { using execution_space = ExecutionSpace; - using real_type = RealType; + using view_type = ViewType; + using BoxViewType = ViewType; + using BoxViewType_hmt = typename ViewType::HostMirror; + using BoxType = typename BoxViewType::value_type::box_type; + using real_type = typename BoxType::value_type; using LBVH_types = MortonLbvhTypes; using kokkos_aabb_types = MortonAabbTypes; @@ -80,8 +122,8 @@ struct MortonAabbTree real_type min_z, real_type max_z); KOKKOS_FORCEINLINE_FUNCTION - void device_set_box(LocalOrdinal box_idx, real_type min_x, real_type max_x, real_type min_y, real_type max_y, - real_type min_z, real_type max_z) const; + void device_set_box(LocalOrdinal box_idx, real_type min_x, real_type min_y, real_type min_z, + real_type max_x, real_type max_y, real_type max_z) const; void sync_from_device() const; void sync_to_device() const; @@ -106,8 +148,10 @@ struct MortonAabbTree local_ordinal_scl_hmt hm_numInternalNodes; local_ordinal_scl_tmt tm_numInternalNodes; - bboxes_3d_view_t m_minMaxs; - bboxes_3d_view_hmt hm_minMaxs; + BoxViewType m_minMaxs; + BoxViewType_hmt hm_minMaxs; + bboxes_3d_view_t m_nodeMinMaxs; + bboxes_3d_view_hmt hm_nodeMinMaxs; local_ordinal_pairs_t m_nodeChildren; local_ordinals_t m_nodeParents; @@ -123,8 +167,8 @@ struct MortonAabbTree #endif }; -template -MortonAabbTree::MortonAabbTree(const std::string &baseName, +template +MortonAabbTree::MortonAabbTree(const std::string &baseName, LocalOrdinal numLeaves, bool supportHostBoxes) : m_baseName(baseName), @@ -134,8 +178,8 @@ MortonAabbTree::MortonAabbTree(const std::string &base init(numLeaves); } -template -void MortonAabbTree::init(LocalOrdinal numLeaves) +template +void MortonAabbTree::init(LocalOrdinal numLeaves) { LocalOrdinal numInternalNodes = std::max(numLeaves - 1, 0); LocalOrdinal numNodes = numLeaves + numInternalNodes; @@ -144,7 +188,8 @@ void MortonAabbTree::init(LocalOrdinal numLeaves) Kokkos::Timer timer0; m_numLeaves = no_init(compound_name(m_baseName, "numLeaves")); m_numInternalNodes = no_init(compound_name(m_baseName, "numInternalNodes")); - m_minMaxs = no_init(compound_name(m_baseName, "minMaxs"), numNodes); + m_minMaxs = no_init(compound_name(m_baseName, "minMaxs"), numLeaves); + m_nodeMinMaxs = no_init(compound_name(m_baseName, "nodeMinMaxs"), numInternalNodes); m_nodeChildren = no_init(compound_name(m_baseName, "children"), numNodes); m_nodeParents = no_init(compound_name(m_baseName, "parents"), numNodes); @@ -176,6 +221,7 @@ void MortonAabbTree::init(LocalOrdinal numLeaves) if (m_supportHostBoxes) { hm_minMaxs = Kokkos::create_mirror_view(Kokkos::WithoutInitializing, m_minMaxs); + hm_nodeMinMaxs = Kokkos::create_mirror_view(Kokkos::WithoutInitializing, m_nodeMinMaxs); } #ifdef DEBUG_MORTON_ACCELERATED_SEARCH @@ -192,16 +238,18 @@ void MortonAabbTree::init(LocalOrdinal numLeaves) } } -template -void MortonAabbTree::resize(LocalOrdinal numLeaves) +template +void MortonAabbTree::resize(LocalOrdinal numLeaves) { if ((numLeaves < 0) || (numLeaves == hm_numLeaves())) return; + Kokkos::Profiling::pushRegion("MortonAabbTree::resize"); LocalOrdinal numInternalNodes = std::max(numLeaves - 1, 0); LocalOrdinal numNodes = numLeaves + numInternalNodes; // Allocate views. - Kokkos::resize(m_minMaxs, numNodes); + Kokkos::resize(m_minMaxs, numLeaves); + Kokkos::resize(m_nodeMinMaxs, numInternalNodes); Kokkos::resize(m_nodeChildren, numNodes); Kokkos::resize(m_nodeParents, numNodes); Kokkos::resize(m_atomicFlags, numInternalNodes); @@ -219,6 +267,7 @@ void MortonAabbTree::resize(LocalOrdinal numLeaves) if (m_supportHostBoxes) { // Host mirror views really for debugging hm_minMaxs = Kokkos::create_mirror_view(Kokkos::WithoutInitializing, m_minMaxs); + hm_nodeMinMaxs = Kokkos::create_mirror_view(Kokkos::WithoutInitializing, m_nodeMinMaxs); } #ifdef DEBUG_MORTON_ACCELERATED_SEARCH @@ -227,46 +276,60 @@ void MortonAabbTree::resize(LocalOrdinal numLeaves) hm_leafIds = Kokkos::create_mirror_view(Kokkos::WithoutInitializing, m_leafIds); hm_leafCodes = Kokkos::create_mirror_view(Kokkos::WithoutInitializing, m_leafCodes); #endif + Kokkos::Profiling::popRegion(); } -template +template KOKKOS_INLINE_FUNCTION -void MortonAabbTree::host_set_box(LocalOrdinal boxIdx, +void set_stk_box(Box& box, RealType minX, RealType minY, RealType minZ, + RealType maxX, RealType maxY, RealType maxZ) +{ + box.min_corner()[0] = minX; + box.min_corner()[1] = minY; + box.min_corner()[2] = minZ; + box.max_corner()[0] = maxX; + box.max_corner()[1] = maxY; + box.max_corner()[2] = maxZ; +} + +template +KOKKOS_INLINE_FUNCTION +void MortonAabbTree::host_set_box(LocalOrdinal boxIdx, real_type minX, real_type maxX, real_type minY, real_type maxY, real_type minZ, real_type maxZ) { - hm_minMaxs(boxIdx, 0) = minX; - hm_minMaxs(boxIdx, 1) = minY; - hm_minMaxs(boxIdx, 2) = minZ; - hm_minMaxs(boxIdx, 3) = maxX; - hm_minMaxs(boxIdx, 4) = maxY; - hm_minMaxs(boxIdx, 5) = maxZ; + hm_minMaxs(boxIdx).box.min_corner()[0] = minX; + hm_minMaxs(boxIdx).box.min_corner()[1] = minY; + hm_minMaxs(boxIdx).box.min_corner()[2] = minZ; + hm_minMaxs(boxIdx).box.max_corner()[0] = maxX; + hm_minMaxs(boxIdx).box.max_corner()[1] = maxY; + hm_minMaxs(boxIdx).box.max_corner()[2] = maxZ; } -template +template KOKKOS_FORCEINLINE_FUNCTION -void MortonAabbTree::device_set_box(LocalOrdinal boxIdx, - real_type minX, real_type maxX, - real_type minY, real_type maxY, - real_type minZ, real_type maxZ) const +void MortonAabbTree::device_set_box(LocalOrdinal boxIdx, + real_type minX, real_type minY, real_type minZ, + real_type maxX, real_type maxY, real_type maxZ) const { - m_minMaxs(boxIdx, 0) = minX; - m_minMaxs(boxIdx, 1) = minY; - m_minMaxs(boxIdx, 2) = minZ; - m_minMaxs(boxIdx, 3) = maxX; - m_minMaxs(boxIdx, 4) = maxY; - m_minMaxs(boxIdx, 5) = maxZ; + m_minMaxs(boxIdx).box.min_corner()[0] = minX; + m_minMaxs(boxIdx).box.min_corner()[1] = minY; + m_minMaxs(boxIdx).box.min_corner()[2] = minZ; + m_minMaxs(boxIdx).box.max_corner()[0] = maxX; + m_minMaxs(boxIdx).box.max_corner()[1] = maxY; + m_minMaxs(boxIdx).box.max_corner()[2] = maxZ; } -template -void MortonAabbTree::sync_from_device() const +template +void MortonAabbTree::sync_from_device() const { Kokkos::deep_copy(hm_numLeaves, m_numLeaves); Kokkos::deep_copy(hm_numInternalNodes, m_numInternalNodes); if (m_supportHostBoxes) { Kokkos::deep_copy(hm_minMaxs, m_minMaxs); + Kokkos::deep_copy(hm_nodeMinMaxs, m_nodeMinMaxs); } #ifdef DEBUG_MORTON_ACCELERATED_SEARCH @@ -278,14 +341,15 @@ void MortonAabbTree::sync_from_device() const #endif } -template -void MortonAabbTree::sync_to_device() const +template +void MortonAabbTree::sync_to_device() const { Kokkos::deep_copy(m_numLeaves, hm_numLeaves); Kokkos::deep_copy(m_numInternalNodes, hm_numInternalNodes); if (m_supportHostBoxes) { Kokkos::deep_copy(m_minMaxs, hm_minMaxs); + Kokkos::deep_copy(m_nodeMinMaxs, hm_nodeMinMaxs); } #ifdef DEBUG_MORTON_ACCELERATED_SEARCH @@ -297,8 +361,8 @@ void MortonAabbTree::sync_to_device() const #endif } -template -std::ostream &MortonAabbTree::streamit(std::ostream &os) const +template +std::ostream &MortonAabbTree::streamit(std::ostream &os) const { os << "{MortonAabbTree " << m_baseName << " " << sizeof(morton_code_t) << " (" << m_globalMinPt[0] << " " << m_globalMinPt[1] << " " << m_globalMinPt[2] @@ -317,8 +381,9 @@ std::ostream &MortonAabbTree::streamit(std::ostream &o << "}" << std::endl; } for (LocalOrdinal idx = hm_numLeaves(); idx < hm_numLeaves()+hm_numInternalNodes(); ++idx) { - os << " {(" << hm_minMaxs(idx, 0) << " " << hm_minMaxs(idx, 1) << " " << hm_minMaxs(idx, 2) - << ") (" << hm_minMaxs(idx, 3) << " " << hm_minMaxs(idx, 4) << " " << hm_minMaxs(idx, 5) << ")" + LocalOrdinal nodeIdx = idx - hm_numLeaves(); + os << " {(" << hm_nodeMinMaxs(nodeIdx, 0) << " " << hm_nodeMinMaxs(nodeIdx, 1) << " " << hm_nodeMinMaxs(nodeIdx, 2) + << ") (" << hm_nodeMinMaxs(nodeIdx, 3) << " " << hm_nodeMinMaxs(nodeIdx, 4) << " " << hm_nodeMinMaxs(nodeIdx, 5) << ")" #ifdef DEBUG_MORTON_ACCELERATED_SEARCH << " (" << hm_nodeChildren(idx, 0) << ", " << hm_nodeChildren(idx, 1) << ") " << hm_nodeParents(idx) #endif @@ -329,12 +394,19 @@ std::ostream &MortonAabbTree::streamit(std::ostream &o return os; } -template -std::ostream &MortonAabbTree::streamit(std::ostream &os, size_t idx) const +template +std::ostream &MortonAabbTree::streamit(std::ostream &os, size_t idx) const { if (m_supportHostBoxes) { - os << " {(" << hm_minMaxs(idx, 0) << " " << hm_minMaxs(idx, 1) << " " << hm_minMaxs(idx, 2) - << ") (" << hm_minMaxs(idx, 3) << " " << hm_minMaxs(idx, 4) << " " << hm_minMaxs(idx, 5) << ")}"; + if (idx < hm_numLeaves()) { + os << " {(" << hm_minMaxs(idx, 0) << " " << hm_minMaxs(idx, 1) << " " << hm_minMaxs(idx, 2) + << ") (" << hm_minMaxs(idx, 3) << " " << hm_minMaxs(idx, 4) << " " << hm_minMaxs(idx, 5) << ")}"; + } + else { + size_t nodeIdx = idx - hm_numLeaves(); + os << " {(" << hm_nodeMinMaxs(nodeIdx, 0) << " " << hm_nodeMinMaxs(nodeIdx, 1) << " " << hm_nodeMinMaxs(nodeIdx, 2) + << ") (" << hm_nodeMinMaxs(nodeIdx, 3) << " " << hm_nodeMinMaxs(nodeIdx, 4) << " " << hm_nodeMinMaxs(nodeIdx, 5) << ")}"; + } } else { os << "MortonAabbTree::streamit(.) m_supportHostBoxes = 0. "; @@ -343,8 +415,8 @@ std::ostream &MortonAabbTree::streamit(std::ostream &o } #ifdef DEBUG_MORTON_ACCELERATED_SEARCH -template -void MortonAabbTree::dump(const std::string& prefix) const +template +void MortonAabbTree::dump(const std::string& prefix) const { std::cout << prefix << ": dump, numLeaves=" << hm_numLeaves() << ", numInternalNodes=" << hm_numInternalNodes() << std::endl; @@ -383,6 +455,16 @@ void MortonAabbTree::dump(const std::string& prefix) c of << std::endl; } } + { + std::ofstream of(prefix+"_nodeMinMaxs.txt"); + for (size_t i = 0; i < hm_nodeMinMaxs.extent(0); ++i) { + of << i; + for (size_t j = 0; j < 6; ++j) { + of << " " << std::setprecision(7) << hm_nodeMinMaxs(i,j); + } + of << std::endl; + } + } } #endif diff --git a/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_TreeManipulationUtils.hpp b/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_TreeManipulationUtils.hpp index f9ef1562721e..b080ccf30b20 100644 --- a/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_TreeManipulationUtils.hpp +++ b/packages/stk/stk_search/stk_search/morton_lbvh/MortonLBVH_TreeManipulationUtils.hpp @@ -113,21 +113,21 @@ namespace stk::search { -template +template struct TotalBoundsFunctor { using size_type = typename ExecutionSpace::size_type; + using BoxType = typename ViewType::value_type::box_type; + using RealType = typename BoxType::value_type; using value_type = MortonAABox; - using kokkos_aabb_types = MortonAabbTypes; - using bboxes_const_3d_view_t = typename kokkos_aabb_types::bboxes_const_3d_view_t; - TotalBoundsFunctor(const MortonAabbTree &tree); + TotalBoundsFunctor(const MortonAabbTree &tree); KOKKOS_INLINE_FUNCTION void init(value_type &update) const; - static void apply(MortonAabbTree &tree, ExecutionSpace const& execSpace); + static void apply(MortonAabbTree &tree, ExecutionSpace const& execSpace); KOKKOS_INLINE_FUNCTION void operator()(size_type idx, value_type &update) const; @@ -135,17 +135,17 @@ struct TotalBoundsFunctor KOKKOS_INLINE_FUNCTION void join(value_type &update, const value_type &input) const; - bboxes_const_3d_view_t m_minMaxs; + ViewType m_minMaxs; }; -template -TotalBoundsFunctor::TotalBoundsFunctor(const MortonAabbTree &tree) +template +TotalBoundsFunctor::TotalBoundsFunctor(const MortonAabbTree &tree) : m_minMaxs(tree.m_minMaxs) {} -template +template KOKKOS_INLINE_FUNCTION -void TotalBoundsFunctor::init(value_type &update) const +void TotalBoundsFunctor::init(value_type &update) const { update.m_min[0] = FLT_MAX; update.m_min[1] = FLT_MAX; @@ -156,8 +156,8 @@ void TotalBoundsFunctor::init(value_type &update) cons update.m_max[2] = -FLT_MAX; } -template -void TotalBoundsFunctor::apply(MortonAabbTree &tree, ExecutionSpace const& execSpace) +template +void TotalBoundsFunctor::apply(MortonAabbTree &tree, ExecutionSpace const& execSpace) { value_type retBox; retBox.m_min[0] = FLT_MAX; @@ -184,22 +184,23 @@ void TotalBoundsFunctor::apply(MortonAabbTree +template KOKKOS_INLINE_FUNCTION -void TotalBoundsFunctor::operator()(size_type idx, value_type &update) const +void TotalBoundsFunctor::operator()(size_type idx, value_type &update) const { - update.m_min[0] = fmin(m_minMaxs(idx, 0), update.m_min[0]); - update.m_min[1] = fmin(m_minMaxs(idx, 1), update.m_min[1]); - update.m_min[2] = fmin(m_minMaxs(idx, 2), update.m_min[2]); - - update.m_max[0] = fmax(m_minMaxs(idx, 3), update.m_max[0]); - update.m_max[1] = fmax(m_minMaxs(idx, 4), update.m_max[1]); - update.m_max[2] = fmax(m_minMaxs(idx, 5), update.m_max[2]); + const BoxType& box = get_box(m_minMaxs(idx)); + update.m_min[0] = fmin(box.get_x_min(), update.m_min[0]); + update.m_min[1] = fmin(box.get_y_min(), update.m_min[1]); + update.m_min[2] = fmin(box.get_z_min(), update.m_min[2]); + + update.m_max[0] = fmax(box.get_x_max(), update.m_max[0]); + update.m_max[1] = fmax(box.get_y_max(), update.m_max[1]); + update.m_max[2] = fmax(box.get_z_max(), update.m_max[2]); } -template +template KOKKOS_INLINE_FUNCTION -void TotalBoundsFunctor::join(value_type &update, const value_type &input) const +void TotalBoundsFunctor::join(value_type &update, const value_type &input) const { update.m_min[0] = fmin(update.m_min[0], input.m_min[0]); update.m_min[1] = fmin(update.m_min[1], input.m_min[1]); @@ -211,24 +212,23 @@ void TotalBoundsFunctor::join(value_type &update, cons } -template +template struct MortonEncoder { using value_type = int; + using BoxType = typename ViewType::value_type::box_type; + using RealType = typename BoxType::value_type; using LBVH_types = MortonLbvhTypes; - using kokkos_aabb_types = MortonAabbTypes; - using bboxes_const_3d_view_t = typename kokkos_aabb_types::bboxes_const_3d_view_t; - using bboxes_3d_view_amt = typename kokkos_aabb_types::bboxes_3d_view_amt; - MortonEncoder(const MortonAabbTree &tree, bool reallyEncode); + MortonEncoder(const MortonAabbTree &tree, bool reallyEncode); - static void apply(const MortonAabbTree &tree, ExecutionSpace const& execSpace, bool reallyEncode = true); + static void apply(const MortonAabbTree &tree, ExecutionSpace const& execSpace, bool reallyEncode = true); KOKKOS_INLINE_FUNCTION void operator()(unsigned leafIdx) const; - bboxes_const_3d_view_t m_minMaxs; + ViewType m_minMaxs; typename LBVH_types::local_ordinals_t m_idsOut; typename LBVH_types::aabb_morton_codes_t m_codesOut; const LocalOrdinal m_numPts; @@ -241,8 +241,8 @@ struct MortonEncoder const bool m_reallyDo; }; -template -MortonEncoder::MortonEncoder(const MortonAabbTree &tree, +template +MortonEncoder::MortonEncoder(const MortonAabbTree &tree, bool reallyEncode) : m_minMaxs(tree.m_minMaxs), m_idsOut(tree.m_leafIds), @@ -258,8 +258,8 @@ MortonEncoder::MortonEncoder(const MortonAabbTree -void MortonEncoder::apply(const MortonAabbTree &tree, +template +void MortonEncoder::apply(const MortonAabbTree &tree, ExecutionSpace const& execSpace, bool reallyEncode) { const MortonEncoder op(tree, reallyEncode); @@ -270,13 +270,14 @@ void MortonEncoder::apply(const MortonAabbTree +template KOKKOS_INLINE_FUNCTION void MortonEncoder::operator()(unsigned leafIdx) const { - RealType ctdX = 0.5 * (m_minMaxs(leafIdx, 0) + m_minMaxs(leafIdx, 3)); - RealType ctdY = 0.5 * (m_minMaxs(leafIdx, 1) + m_minMaxs(leafIdx, 4)); - RealType ctdZ = 0.5 * (m_minMaxs(leafIdx, 2) + m_minMaxs(leafIdx, 5)); + const BoxType& box = get_box(m_minMaxs(leafIdx)); + const RealType ctdX = 0.5 * (box.get_x_min() + box.get_x_max()); + const RealType ctdY = 0.5 * (box.get_y_min() + box.get_y_max()); + const RealType ctdZ = 0.5 * (box.get_z_min() + box.get_z_max()); // std::cout << "box(" << leafIdx << ") = (" << m_minMax(leafIdx, 0) << " " // << m_minMax(leafIdx, 1) << " " << m_minMax(leafIdx, 2) @@ -318,16 +319,17 @@ void MortonEncoder::operator()(unsigned leafIdx) const #else // 64 bit Morton codes -template +template KOKKOS_INLINE_FUNCTION -void MortonEncoder::operator()(unsigned leafIdx) const +void MortonEncoder::operator()(unsigned leafIdx) const { m_idsOut(leafIdx) = leafIdx; if (m_reallyDo) { - RealType ctdX = 0.5 * (m_minMaxs(leafIdx, 0) + m_minMaxs(leafIdx, 3)); - RealType ctdY = 0.5 * (m_minMaxs(leafIdx, 1) + m_minMaxs(leafIdx, 4)); - RealType ctdZ = 0.5 * (m_minMaxs(leafIdx, 2) + m_minMaxs(leafIdx, 5)); + const BoxType& box = get_box(m_minMaxs(leafIdx)); + const RealType ctdX = 0.5 * (box.get_x_min() + box.get_x_max()); + const RealType ctdY = 0.5 * (box.get_y_min() + box.get_y_max()); + const RealType ctdZ = 0.5 * (box.get_z_min() + box.get_z_max()); // std::cout << "box(" << leafIdx << ") = (" << m_minMaxs(leafIdx, 0) << " " // << m_minMaxs(leafIdx, 1) << " " << m_minMaxs(leafIdx, 2) @@ -346,7 +348,7 @@ void MortonEncoder::operator()(unsigned leafIdx) const ux = static_cast((ctdX - m_globalXMin) * 2097151.0f / m_xWidth); uy = static_cast((ctdY - m_globalYMin) * 2097151.0f / m_yWidth); - uz = static_cast((ctdZ = - m_globalZMin) * 2097151.0f / m_zWidth); + uz = static_cast((ctdZ - m_globalZMin) * 2097151.0f / m_zWidth); ux = (ux | ux << 32) & 0x001f00000000ffff; ux = (ux | ux << 16) & 0x001f0000ff0000ff; @@ -374,22 +376,22 @@ void MortonEncoder::operator()(unsigned leafIdx) const // Serial sort is the default. -template +template struct SortByCodeIdPair { using LBVH_types = MortonLbvhTypes; - SortByCodeIdPair(const MortonAabbTree &tree); + SortByCodeIdPair(const TreeType &tree); - static void apply(const MortonAabbTree &tree, bool reallyEncode = true); + static void apply(const TreeType &tree, bool reallyEncode = true); std::vector m_buffer; typename LBVH_types::local_ordinals_hmt hm_leafIds; typename LBVH_types::aabb_morton_codes_hmt hm_leafCodes; }; -template -SortByCodeIdPair::SortByCodeIdPair(const MortonAabbTree &tree) +template +SortByCodeIdPair::SortByCodeIdPair(const TreeType &tree) { hm_leafIds = Kokkos::create_mirror_view(Kokkos::WithoutInitializing, tree.m_leafIds); hm_leafCodes = Kokkos::create_mirror_view(Kokkos::WithoutInitializing, tree.m_leafCodes); @@ -405,8 +407,8 @@ SortByCodeIdPair::SortByCodeIdPair(const MortonAabbTre } } -template -void SortByCodeIdPair::apply(const MortonAabbTree &tree, +template +void SortByCodeIdPair::apply(const TreeType &tree, bool reallyEncode) { SortByCodeIdPair tmp(tree); @@ -423,13 +425,13 @@ void SortByCodeIdPair::apply(const MortonAabbTree +template struct SortByCode { - static void apply(const MortonAabbTree &tree, ExecutionSpace const& execSpace) + static void apply(const TreeType &tree, ExecutionSpace const& execSpace) { if constexpr (Kokkos::SpaceAccessibility::accessible) { - SortByCodeIdPair::apply(tree); + SortByCodeIdPair::apply(tree); } else { //#if KOKKOS_VERSION >= 40300 @@ -445,20 +447,20 @@ struct SortByCode //thrust::stable_sort_by_key(rawLeafCodesThr, rawLeafCodesThr + n, rawLeafIdsThr); thrust::sort_by_key(rawLeafCodesThr, rawLeafCodesThr + n, rawLeafIdsThr); #else - STK_ThrowErrorMsg("shouldn't be able to get here"); // SortByCodeIdPair::apply(tree); + STK_ThrowErrorMsg("shouldn't be able to get here"); // SortByCodeIdPair::apply(tree); #endif } } }; -template +template struct BuildRadixTree { using LBVH_types = MortonLbvhTypes; - BuildRadixTree(const MortonAabbTree &tree); + BuildRadixTree(const TreeType &tree); - static void apply(const MortonAabbTree &tree, ExecutionSpace const& execSpace); + static void apply(const TreeType &tree, ExecutionSpace const& execSpace); KOKKOS_INLINE_FUNCTION void operator()(unsigned argIdx) const; @@ -475,8 +477,8 @@ struct BuildRadixTree typename LBVH_types::local_ordinals_t m_atomicFlags; }; -template -BuildRadixTree::BuildRadixTree(const MortonAabbTree &tree) +template +BuildRadixTree::BuildRadixTree(const TreeType &tree) : m_numLeaves(tree.hm_numLeaves()), m_numInternalNodes(tree.hm_numInternalNodes()), tm_leafCodes(tree.m_leafCodes), @@ -486,20 +488,20 @@ BuildRadixTree::BuildRadixTree(const MortonAabbTree -void BuildRadixTree::apply(const MortonAabbTree &tree, ExecutionSpace const& execSpace) +template +void BuildRadixTree::apply(const TreeType &tree, ExecutionSpace const& execSpace) { if (tree.hm_numLeaves() <= 0) { return; } - BuildRadixTree op(tree); + BuildRadixTree op(tree); auto policy = Kokkos::RangePolicy(execSpace, 0, static_cast(tree.hm_numInternalNodes())); Kokkos::parallel_for(policy, op); } -template +template KOKKOS_INLINE_FUNCTION -void BuildRadixTree::operator()(unsigned argIdx) const +void BuildRadixTree::operator()(unsigned argIdx) const { LocalOrdinal idx = static_cast(argIdx); @@ -572,9 +574,9 @@ void BuildRadixTree::operator()(unsigned argIdx) const #ifdef SMALL_MORTON // 32 bit Morton -template +template KOKKOS_INLINE_FUNCTION -int leaves_cpr(LocalOrdinal baseIdx, LocalOrdinal testIdx) const +int BuildRadixTree::leaves_cpr(LocalOrdinal baseIdx, LocalOrdinal testIdx) const { if (testIdx < 0 || testIdx >= m_numLeaves) { return -1; @@ -591,9 +593,9 @@ int leaves_cpr(LocalOrdinal baseIdx, LocalOrdinal testIdx) const #else // 64 bit Morton -template +template KOKKOS_INLINE_FUNCTION -int BuildRadixTree::leaves_cpr(LocalOrdinal baseIdx, LocalOrdinal testIdx) const +int BuildRadixTree::leaves_cpr(LocalOrdinal baseIdx, LocalOrdinal testIdx) const { if (testIdx < 0 || testIdx >= m_numLeaves) { return -1; @@ -610,65 +612,69 @@ int BuildRadixTree::leaves_cpr(LocalOrdinal baseIdx, L #endif // 64 bit Morton -template +template struct UpdateInteriorNodeBVs { using LBVH_types = MortonLbvhTypes; + using BoxType = typename ViewType::value_type::box_type; + using RealType = typename BoxType::value_type; using kokkos_aabb_types = MortonAabbTypes; - using bboxes_const_3d_view_t = typename kokkos_aabb_types::bboxes_const_3d_view_t; + using bboxes_3d_view_amt = typename kokkos_aabb_types::bboxes_3d_view_amt; - UpdateInteriorNodeBVs(const MortonAabbTree &tree); + UpdateInteriorNodeBVs(const MortonAabbTree &tree); - static void apply(const MortonAabbTree &tree, ExecutionSpace const& execSpace); + static void apply(const MortonAabbTree &tree, ExecutionSpace const& execSpace); KOKKOS_INLINE_FUNCTION void operator()(unsigned argIdx) const; - template KOKKOS_FORCEINLINE_FUNCTION - void get_box(RealType bvMinMax[6], LocalOrdinal idx, const BBox3dViewType &boxesMinMax) const; + void get_box(RealType bvMinMax[6], LocalOrdinal idx, const bboxes_3d_view_amt &boxesMinMax) const; + + KOKKOS_FORCEINLINE_FUNCTION + void get_stk_box(RealType bvMinMax[6], LocalOrdinal idx, const ViewType &boxesMinMax) const; const LocalOrdinal m_numLeaves; const LocalOrdinal m_numInternalNodes; typename LBVH_types::local_ordinal_pairs_tmt tm_nodeChildren; typename LBVH_types::local_ordinals_tmt tm_nodeParents; - bboxes_const_3d_view_t m_leafMinMaxs; + ViewType m_leafMinMaxs; // Will write to internal nodes bounding boxes. - typename kokkos_aabb_types::bboxes_3d_view_amt m_nodeMinMaxs; + bboxes_3d_view_amt m_nodeMinMaxs; typename LBVH_types::local_ordinals_t m_atomicFlags; }; -template -UpdateInteriorNodeBVs::UpdateInteriorNodeBVs(const MortonAabbTree &tree) +template +UpdateInteriorNodeBVs::UpdateInteriorNodeBVs(const MortonAabbTree &tree) : m_numLeaves(tree.hm_numLeaves()), m_numInternalNodes(tree.hm_numInternalNodes()), tm_nodeChildren(tree.m_nodeChildren), tm_nodeParents(tree.m_nodeParents), m_leafMinMaxs(tree.m_minMaxs), - m_nodeMinMaxs(tree.m_minMaxs), + m_nodeMinMaxs(tree.m_nodeMinMaxs), m_atomicFlags(tree.m_atomicFlags) {} -template -void UpdateInteriorNodeBVs::apply(const MortonAabbTree &tree, ExecutionSpace const& execSpace) +template +void UpdateInteriorNodeBVs::apply(const MortonAabbTree &tree, ExecutionSpace const& execSpace) { - const UpdateInteriorNodeBVs op(tree); + const UpdateInteriorNodeBVs op(tree); const size_t numLeaves = tree.hm_numLeaves(); auto policy = Kokkos::RangePolicy(execSpace, 0, numLeaves); Kokkos::parallel_for(policy, op); } -template +template KOKKOS_INLINE_FUNCTION -void UpdateInteriorNodeBVs::operator()(unsigned argIdx) const +void UpdateInteriorNodeBVs::operator()(unsigned argIdx) const { if (m_numLeaves > 1) { LocalOrdinal idx = static_cast(argIdx); RealType bvMinMax[6]; - get_box(bvMinMax, idx, m_leafMinMaxs); + get_stk_box(bvMinMax, idx, m_leafMinMaxs); LocalOrdinal parent = tm_nodeParents(idx); RealType sibMinMax[6]; @@ -679,17 +685,28 @@ void UpdateInteriorNodeBVs::operator()(unsigned argIdx sib = tm_nodeChildren(parent, 1); } - get_box(sibMinMax, sib, m_nodeMinMaxs); - - for (LocalOrdinal j = 0; j < 3; ++j) { - bvMinMax[j] = AABB_MIN(bvMinMax[j], sibMinMax[j]); - m_nodeMinMaxs(parent, j) = bvMinMax[j]; + if (sib < m_numLeaves) { + get_stk_box(sibMinMax, sib, m_leafMinMaxs); } - for (LocalOrdinal j = 3; j < 6; ++j) { - bvMinMax[j] = AABB_MAX(bvMinMax[j], sibMinMax[j]); - m_nodeMinMaxs(parent, j) = bvMinMax[j]; + else { + get_box(sibMinMax, sib-m_numLeaves, m_nodeMinMaxs); } + const LocalOrdinal thisIdx = parent - m_numLeaves; + + bvMinMax[0] = AABB_MIN(bvMinMax[0], sibMinMax[0]); + bvMinMax[1] = AABB_MIN(bvMinMax[1], sibMinMax[1]); + bvMinMax[2] = AABB_MIN(bvMinMax[2], sibMinMax[2]); + bvMinMax[3] = AABB_MAX(bvMinMax[3], sibMinMax[3]); + bvMinMax[4] = AABB_MAX(bvMinMax[4], sibMinMax[4]); + bvMinMax[5] = AABB_MAX(bvMinMax[5], sibMinMax[5]); + m_nodeMinMaxs(thisIdx, 0) = bvMinMax[0]; + m_nodeMinMaxs(thisIdx, 1) = bvMinMax[1]; + m_nodeMinMaxs(thisIdx, 2) = bvMinMax[2]; + m_nodeMinMaxs(thisIdx, 3) = bvMinMax[3]; + m_nodeMinMaxs(thisIdx, 4) = bvMinMax[4]; + m_nodeMinMaxs(thisIdx, 5) = bvMinMax[5]; + idx = parent; parent = tm_nodeParents(parent); if (idx == parent) { @@ -699,15 +716,31 @@ void UpdateInteriorNodeBVs::operator()(unsigned argIdx } } -template -template +template KOKKOS_FORCEINLINE_FUNCTION -void UpdateInteriorNodeBVs::get_box(RealType bvMinMax[6], LocalOrdinal idx, - const BBox3dViewType &boxMinMaxs) const +void UpdateInteriorNodeBVs::get_box(RealType bvMinMax[6], LocalOrdinal idx, + const bboxes_3d_view_amt &boxMinMaxs) const { - for (LocalOrdinal j = 0; j < 6; ++j) { - bvMinMax[j] = boxMinMaxs(idx, j); - } + bvMinMax[0] = boxMinMaxs(idx, 0); + bvMinMax[1] = boxMinMaxs(idx, 1); + bvMinMax[2] = boxMinMaxs(idx, 2); + bvMinMax[3] = boxMinMaxs(idx, 3); + bvMinMax[4] = boxMinMaxs(idx, 4); + bvMinMax[5] = boxMinMaxs(idx, 5); +} + +template +KOKKOS_FORCEINLINE_FUNCTION +void UpdateInteriorNodeBVs::get_stk_box(RealType bvMinMax[6], LocalOrdinal idx, + const ViewType &boxMinMaxs) const +{ + const BoxType& box = stk::search::get_box(boxMinMaxs(idx)); + bvMinMax[0] = box.get_x_min(); + bvMinMax[1] = box.get_y_min(); + bvMinMax[2] = box.get_z_min(); + bvMinMax[3] = box.get_x_max(); + bvMinMax[4] = box.get_y_max(); + bvMinMax[5] = box.get_z_max(); } @@ -744,9 +777,12 @@ class CollisionListCallback }; -template +template struct Traverse_MASTB_BVH_Functor { + using DomainBoxType = typename DomainViewType::value_type::box_type; + using RangeBoxType = typename RangeViewType::value_type::box_type; + using RealType = typename RangeBoxType::value_type; using LBVH_types = MortonLbvhTypes; using kokkos_aabb_types = MortonAabbTypes; using local_ordinals_tmt = typename LBVH_types::local_ordinals_tmt; @@ -754,9 +790,9 @@ struct Traverse_MASTB_BVH_Functor using bboxes_const_3d_view_t = typename kokkos_aabb_types::bboxes_const_3d_view_t; using collision_list_type = CollisionList; - Traverse_MASTB_BVH_Functor(bboxes_3d_view_t domainMinMaxs, + Traverse_MASTB_BVH_Functor(DomainViewType domainMinMaxs, local_ordinals_tmt domainIds, - const MortonAabbTree &rangeTree, + const MortonAabbTree &rangeTree, Callback& callback, bool flippedResults = false); @@ -770,7 +806,7 @@ struct Traverse_MASTB_BVH_Functor bool is_range_leaf(LocalOrdinal rangeIdx) const{ return (rangeIdx < m_rangeRoot); } KOKKOS_FORCEINLINE_FUNCTION - void get_box(RealType bvMinMax[6], LocalOrdinal idx, const bboxes_const_3d_view_t &boxMinMaxs) const; + void get_box(RealType bvMinMax[6], LocalOrdinal idx, const DomainViewType &boxMinMaxs) const; std::ostream &stream_pair(LocalOrdinal domainIdx, bool overlap, LocalOrdinal rangeIdx, std::ostream &os) const; @@ -782,11 +818,12 @@ struct Traverse_MASTB_BVH_Functor m_callback(domainIdxFlipped, rangeIdxFlipped); } - bboxes_const_3d_view_t m_domainMinMaxs; + DomainViewType m_domainMinMaxs; typename LBVH_types::local_ordinals_tmt tm_domainIds; const LocalOrdinal m_rangeRoot; - bboxes_const_3d_view_t tm_rangeMinMaxs; + RangeViewType tm_rangeMinMaxs; + bboxes_const_3d_view_t tm_rangeNodeMinMaxs; typename LBVH_types::local_ordinal_pairs_tmt tm_rangeNodeChildren; typename LBVH_types::local_ordinals_tmt tm_rangeLeafIds; @@ -794,17 +831,18 @@ struct Traverse_MASTB_BVH_Functor Callback m_callback; }; -template -Traverse_MASTB_BVH_Functor::Traverse_MASTB_BVH_Functor( - bboxes_3d_view_t domainMinMaxs, +template +Traverse_MASTB_BVH_Functor::Traverse_MASTB_BVH_Functor( + DomainViewType domainMinMaxs, local_ordinals_tmt domainIds, - const MortonAabbTree &rangeTree, + const MortonAabbTree &rangeTree, Callback& callback, bool flippedResults) : m_domainMinMaxs(domainMinMaxs), tm_domainIds(domainIds), m_rangeRoot(rangeTree.hm_numLeaves()), tm_rangeMinMaxs(rangeTree.m_minMaxs), + tm_rangeNodeMinMaxs(rangeTree.m_nodeMinMaxs), tm_rangeNodeChildren(rangeTree.m_nodeChildren), tm_rangeLeafIds(rangeTree.m_leafIds), m_flippedResults(flippedResults), @@ -812,10 +850,10 @@ Traverse_MASTB_BVH_Functor::Traverse_MASTB_B {} -template +template void search_tree( - const MortonAabbTree &domainTree, - const MortonAabbTree &rangeTree, + const Tree1Type &domainTree, + const Tree2Type &rangeTree, Callback& callback, ExecutionSpace const& execSpace, bool flipOutputPairs = false) @@ -826,24 +864,34 @@ void search_tree( return; } - const Traverse_MASTB_BVH_Functor op(domainTree.m_minMaxs, domainTree.m_leafIds, rangeTree, + using DomainViewType = typename Tree1Type::view_type; + using RangeViewType = typename Tree2Type::view_type; + Kokkos::Profiling::pushRegion("construct Traverse_MASTB_BVH_Functor"); + const Traverse_MASTB_BVH_Functor op(domainTree.m_minMaxs, domainTree.m_leafIds, rangeTree, callback, flipOutputPairs); + Kokkos::Profiling::popRegion(); auto policy = Kokkos::RangePolicy(execSpace, 0, domainTree.hm_numLeaves()); - Kokkos::parallel_for("Traverse_MASTB_BVH_Functor", policy, op); + Kokkos::parallel_for("pll-for Traverse_MASTB_BVH_Functor", policy, op); + Kokkos::Profiling::pushRegion("search_tree - fence after pll-for"); execSpace.fence(); + Kokkos::Profiling::popRegion(); + Kokkos::Profiling::pushRegion("search_tree - resize for second pass"); if (callback.resize_for_second_pass()) { - const Traverse_MASTB_BVH_Functor op2(domainTree.m_minMaxs, domainTree.m_leafIds, rangeTree, + const Traverse_MASTB_BVH_Functor op2(domainTree.m_minMaxs, domainTree.m_leafIds, rangeTree, callback, flipOutputPairs); Kokkos::parallel_for("Traverse_MASTB_BVH_Functor - pass2", policy, op2); } + Kokkos::Profiling::popRegion(); + Kokkos::Profiling::pushRegion("search_tree - last fence"); execSpace.fence(); Kokkos::Profiling::popRegion(); + Kokkos::Profiling::popRegion(); } -template -KOKKOS_INLINE_FUNCTION void Traverse_MASTB_BVH_Functor::operator()(unsigned argDomainIdx) const +template +KOKKOS_INLINE_FUNCTION void Traverse_MASTB_BVH_Functor::operator()(unsigned argDomainIdx) const { LocalOrdinal domainIdx = tm_domainIds(argDomainIdx); @@ -907,42 +955,57 @@ KOKKOS_INLINE_FUNCTION void Traverse_MASTB_BVH_Functor +template KOKKOS_FORCEINLINE_FUNCTION -bool Traverse_MASTB_BVH_Functor::overlaps_range(RealType bvMinMax[6], - LocalOrdinal rangeIdx) const +bool Traverse_MASTB_BVH_Functor::overlaps_range(RealType bvMinMax[6], + LocalOrdinal rangeIdx) const { - return (bvMinMax[3] < tm_rangeMinMaxs(rangeIdx, 0) || - bvMinMax[4] < tm_rangeMinMaxs(rangeIdx, 1) || - bvMinMax[5] < tm_rangeMinMaxs(rangeIdx, 2) || - bvMinMax[0] > tm_rangeMinMaxs(rangeIdx, 3) || - bvMinMax[1] > tm_rangeMinMaxs(rangeIdx, 4) || - bvMinMax[2] > tm_rangeMinMaxs(rangeIdx, 5)) ? false : true; + if(rangeIdx < m_rangeRoot) { + const RangeBoxType& box = stk::search::get_box(tm_rangeMinMaxs(rangeIdx)); + + return (bvMinMax[3] < box.get_x_min() || + bvMinMax[4] < box.get_y_min() || + bvMinMax[5] < box.get_z_min() || + bvMinMax[0] > box.get_x_max() || + bvMinMax[1] > box.get_y_max() || + bvMinMax[2] > box.get_z_max()) ? false : true; + } + + const LocalOrdinal rangeNodeIdx = rangeIdx - m_rangeRoot; + return (bvMinMax[3] < tm_rangeNodeMinMaxs(rangeNodeIdx,0) || + bvMinMax[4] < tm_rangeNodeMinMaxs(rangeNodeIdx,1) || + bvMinMax[5] < tm_rangeNodeMinMaxs(rangeNodeIdx,2) || + bvMinMax[0] > tm_rangeNodeMinMaxs(rangeNodeIdx,3) || + bvMinMax[1] > tm_rangeNodeMinMaxs(rangeNodeIdx,4) || + bvMinMax[2] > tm_rangeNodeMinMaxs(rangeNodeIdx,5)) ? false : true; } -template +template KOKKOS_FORCEINLINE_FUNCTION -void Traverse_MASTB_BVH_Functor::get_box(RealType bvMinMax[6], LocalOrdinal idx, - const bboxes_const_3d_view_t &boxMinMaxs) const +void Traverse_MASTB_BVH_Functor::get_box(RealType bvMinMax[6], LocalOrdinal idx, + const DomainViewType &boxMinMaxs) const { - bvMinMax[0] = boxMinMaxs(idx, 0); - bvMinMax[1] = boxMinMaxs(idx, 1); - bvMinMax[2] = boxMinMaxs(idx, 2); - bvMinMax[3] = boxMinMaxs(idx, 3); - bvMinMax[4] = boxMinMaxs(idx, 4); - bvMinMax[5] = boxMinMaxs(idx, 5); + const DomainBoxType& box = stk::search::get_box(boxMinMaxs(idx)); + bvMinMax[0] = box.get_x_min(); + bvMinMax[1] = box.get_y_min(); + bvMinMax[2] = box.get_z_min(); + bvMinMax[3] = box.get_x_max(); + bvMinMax[4] = box.get_y_max(); + bvMinMax[5] = box.get_z_max(); } -template -std::ostream &Traverse_MASTB_BVH_Functor::stream_pair(LocalOrdinal domainIdx, bool overlap, +template +std::ostream &Traverse_MASTB_BVH_Functor::stream_pair(LocalOrdinal domainIdx, bool overlap, LocalOrdinal rangeIdx, std::ostream &os) const { - os << " {(" << m_domainMinMaxs(domainIdx, 0) << "," << m_domainMinMaxs(domainIdx, 1) << "," << m_domainMinMaxs(domainIdx, 2) - << ") (" << m_domainMinMaxs(domainIdx, 3) << "," << m_domainMinMaxs(domainIdx, 4) << "," << m_domainMinMaxs(domainIdx, 5) + auto dbox = get_box(m_domainMinMaxs(domainIdx)); + os << " {(" << dbox.get_x_min() << "," << dbox.get_y_min() << "," <::execution_space; +#ifdef KOKKOS_ENABLE_CUDA return std::is_same::value; +#elif defined(KOKKOS_ENABLE_HIP) + return std::is_same_v; +#endif #else return false; #endif diff --git a/packages/stk/stk_tools/stk_tools/mesh_tools/DetectHingesImpl.cpp b/packages/stk/stk_tools/stk_tools/mesh_tools/DetectHingesImpl.cpp index 7fc1773d64d5..53321ab47256 100644 --- a/packages/stk/stk_tools/stk_tools/mesh_tools/DetectHingesImpl.cpp +++ b/packages/stk/stk_tools/stk_tools/mesh_tools/DetectHingesImpl.cpp @@ -51,11 +51,10 @@ namespace impl { void print_node_count(stk::mesh::BulkData& bulk, const std::string str) { - stk::mesh::EntityVector nodes; - stk::mesh::get_entities(bulk, stk::topology::NODE_RANK, nodes); + const unsigned nodeCount = stk::mesh::count_entities(bulk, stk::topology::NODE_RANK, bulk.mesh_meta_data().universal_part()); std::cout << str << std::endl; - std::cout << "p:" << bulk.parallel_rank() << " node vec size: " << nodes.size() << std::endl; + std::cout << "p:" << bulk.parallel_rank() << " num nodes: " << nodeCount << std::endl; } bool are_equal(const stk::mesh::Entity* nodesPtr, const stk::mesh::EntityVector& nodesVec) @@ -317,8 +316,8 @@ stk::mesh::EntityVector get_mesh_nodes(const stk::mesh::BulkData& bulk, const st parts.push_back(part); } - selector = stk::mesh::selectUnion(parts); } + selector = stk::mesh::selectUnion(parts); if(parts.empty()) { selector = meta.universal_part(); diff --git a/packages/stk/stk_tools/stk_tools/mesh_tools/DetectHingesImpl.hpp b/packages/stk/stk_tools/stk_tools/mesh_tools/DetectHingesImpl.hpp index 15e97382bbd1..e5537b25a5fa 100644 --- a/packages/stk/stk_tools/stk_tools/mesh_tools/DetectHingesImpl.hpp +++ b/packages/stk/stk_tools/stk_tools/mesh_tools/DetectHingesImpl.hpp @@ -35,12 +35,12 @@ #ifndef _DetectHingesImpl_hpp_ #define _DetectHingesImpl_hpp_ -#include #include #include // for MetaData #include namespace stk { +namespace mesh { class BulkData; } namespace tools { namespace impl { diff --git a/packages/stk/stk_tools/stk_tools/mesh_tools/DisconnectBlocksImpl.hpp b/packages/stk/stk_tools/stk_tools/mesh_tools/DisconnectBlocksImpl.hpp index 2a1970a019ac..1fb491adadc6 100644 --- a/packages/stk/stk_tools/stk_tools/mesh_tools/DisconnectBlocksImpl.hpp +++ b/packages/stk/stk_tools/stk_tools/mesh_tools/DisconnectBlocksImpl.hpp @@ -49,8 +49,6 @@ #include #include -namespace stk { namespace mesh { class BulkData; } } - namespace stk { namespace tools { namespace impl { diff --git a/packages/stk/stk_tools/stk_tools/mesh_tools/DisconnectGroup.hpp b/packages/stk/stk_tools/stk_tools/mesh_tools/DisconnectGroup.hpp index 50d035712936..ba1c6926519a 100644 --- a/packages/stk/stk_tools/stk_tools/mesh_tools/DisconnectGroup.hpp +++ b/packages/stk/stk_tools/stk_tools/mesh_tools/DisconnectGroup.hpp @@ -34,7 +34,6 @@ #ifndef _DisconnectGroup_hpp_ #define _DisconnectGroup_hpp_ -#include "stk_mesh/base/BulkData.hpp" #include "stk_mesh/base/Entity.hpp" #include "stk_mesh/base/MetaData.hpp" #include "stk_mesh/base/Part.hpp" @@ -47,6 +46,7 @@ #include namespace stk { +namespace mesh { class BulkData; } namespace tools { namespace impl { diff --git a/packages/stk/stk_tools/stk_tools/transfer_utils/TransientFieldTransferById.hpp b/packages/stk/stk_tools/stk_tools/transfer_utils/TransientFieldTransferById.hpp index 0e4a30995e2f..46e8ced2a53e 100644 --- a/packages/stk/stk_tools/stk_tools/transfer_utils/TransientFieldTransferById.hpp +++ b/packages/stk/stk_tools/stk_tools/transfer_utils/TransientFieldTransferById.hpp @@ -37,7 +37,6 @@ #include #include -#include #include "stk_transfer/copy_by_id/TransferCopyById.hpp" #include "stk_transfer/copy_by_id/TransferCopyByIdStkMeshAdapter.hpp" @@ -48,6 +47,7 @@ namespace Ioss { class Region; } namespace stk { +namespace mesh { class BulkData; } namespace transfer_utils { class RepeatedTransferCopyByIdStkMeshAdapter : public stk::transfer::TransferCopyByIdStkMeshAdapter diff --git a/packages/stk/stk_unit_test_utils/stk_unit_test_utils/TextMesh.cpp b/packages/stk/stk_unit_test_utils/stk_unit_test_utils/TextMesh.cpp index 437be5a1e2be..0369bbfe69a8 100644 --- a/packages/stk/stk_unit_test_utils/stk_unit_test_utils/TextMesh.cpp +++ b/packages/stk/stk_unit_test_utils/stk_unit_test_utils/TextMesh.cpp @@ -301,8 +301,9 @@ class MetaDataInitializer class BulkDataInitializer { public: - BulkDataInitializer(const TextMeshData& d, stk::mesh::BulkData& b) - : m_data(d), + BulkDataInitializer(const StkTopologyMapping& t, const TextMeshData& d, stk::mesh::BulkData& b) + : m_topologyMapping(t), + m_data(d), m_bulk(b), m_meta(m_bulk.mesh_meta_data()) { } @@ -352,6 +353,92 @@ class BulkDataInitializer } } + stk::topology get_side_block_topology_from_entries(const SidesetData& ss, const SideBlockInfo& sideBlock) + { + std::vector sideIndices = ss.get_sideblock_indices_local_to_proc(sideBlock, m_bulk.parallel_rank()); + stk::topology invalidTopo = stk::topology::INVALID_TOPOLOGY; + stk::topology blockSideTopo = stk::topology::INVALID_TOPOLOGY; + int heterogenousTopo = 0; + + for (size_t sideIndex : sideIndices) { + stk::mesh::EntityId elemId = ss.data[sideIndex].first; + int elemSide = ss.data[sideIndex].second - 1; + + auto elemIter = text_mesh::bound_search(m_data.elementDataVec.begin(), m_data.elementDataVec.end(), elemId, ElementDataLess()); + STK_ThrowRequire(elemIter != m_data.elementDataVec.end()); + + Topology blockTopo = elemIter->topology; + Topology sideTopo = blockTopo.side_topology(elemSide); + + if(stk::topology::INVALID_TOPOLOGY != blockSideTopo && sideTopo.topology != blockSideTopo) { + blockSideTopo = stk::topology::INVALID_TOPOLOGY; + heterogenousTopo = 1; + break; + } + + blockSideTopo = sideTopo.topology; + } + + int topoId = (stk::topology::topology_t) blockSideTopo; + + if (m_bulk.parallel_size()) { + std::vector l_topoData{topoId, heterogenousTopo}; + std::vector g_topoData(2); + stk::all_reduce_max(m_bulk.parallel(), l_topoData.data(), g_topoData.data(), 2); + + topoId = g_topoData[0]; + heterogenousTopo = g_topoData[1]; + } + + blockSideTopo = heterogenousTopo ? invalidTopo : stk::topology( (stk::topology::topology_t)topoId ); + + return blockSideTopo; + } + + bool set_sideset_topology(const SidesetData& ss, stk::mesh::Part* part, + const std::string& sideBlockTopo, stk::topology stkTopology, + bool printWarning = false) + { + if (stkTopology == stk::topology::INVALID_TOPOLOGY) { + if (printWarning) { + std::ostringstream os; + os<<"TextMesh WARNING: failed to obtain sensible topology for sideset: " << ss.name<<", iossTopology: "<name()<<", rank: "<primary_entity_rank()<<", stk-topology: "<subsets()) { + if(subsetPart->topology() != stkTopology) { + return false; + } + } + + stk::mesh::set_topology(*part, stkTopology); + return true; + } + + return false; + } + + + void setup_sideset_topology(const SidesetData& ss, stk::mesh::Part* part) + { + if(nullptr == part) return; + + std::vector sideBlocks = ss.get_side_block_info(); + if(sideBlocks.size() != 1u) return; + + const SideBlockInfo& sideBlock = sideBlocks[0]; + + if(sideBlock.name != ss.name) { + stk::topology stkTopology = m_topologyMapping.topology(sideBlock.sideTopology).topology; + if(set_sideset_topology(ss, part, sideBlock.sideTopology, stkTopology, true)) return; + } + + stk::topology blockSideTopo = get_side_block_topology_from_entries(ss, sideBlock); + set_sideset_topology(ss, part, sideBlock.sideTopology, blockSideTopo); + } + void setup_sidesets() { const bool fromInput = true; @@ -359,6 +446,8 @@ class BulkDataInitializer stk::mesh::Part* part = m_meta.get_part(sidesetData.name); stk::mesh::SideSet& stkSideSet = m_bulk.create_sideset(*part, fromInput); + setup_sideset_topology(sidesetData, part); + std::vector sideBlocks = sidesetData.get_side_block_info(); for (const auto& sideBlock : sideBlocks) { @@ -441,6 +530,7 @@ class BulkDataInitializer } } + const StkTopologyMapping &m_topologyMapping; const TextMeshData& m_data; stk::mesh::BulkData& m_bulk; @@ -523,7 +613,7 @@ class TextMesh MetaDataInitializer metaInit(m_topologyMapping, m_data, m_meta); metaInit.setup(); - BulkDataInitializer bulkInit(m_data, m_bulk); + BulkDataInitializer bulkInit(m_topologyMapping, m_data, m_bulk); bulkInit.setup(); CoordinateInitializer coordInit(m_data, m_bulk); diff --git a/packages/stk/stk_unit_test_utils/stk_unit_test_utils/stk_transfer_fixtures/ConservativeTransferFixture.cpp b/packages/stk/stk_unit_test_utils/stk_unit_test_utils/stk_transfer_fixtures/ConservativeTransferFixture.cpp index 9c6c2f63a72d..0610f1db3af8 100644 --- a/packages/stk/stk_unit_test_utils/stk_unit_test_utils/stk_transfer_fixtures/ConservativeTransferFixture.cpp +++ b/packages/stk/stk_unit_test_utils/stk_unit_test_utils/stk_transfer_fixtures/ConservativeTransferFixture.cpp @@ -229,7 +229,7 @@ void ConservativeTransferTests::test_conservation_bidirectional(std::function sendFieldPtr, recvFieldPtr; - double sendIntegral, recvIntegral; + double sendIntegral = 0.0, recvIntegral = 0.0; if (inputMesh1) { sendFieldPtr = mesh::create_field(inputMesh1, mesh::FieldShape(1, 0, 0), 1); @@ -317,4 +317,4 @@ void ConservativeTransferTests::set_field(mesh::FieldPtr fieldPtr, } } -} \ No newline at end of file +} diff --git a/packages/stk/stk_unit_tests/stk_expreval/UnitTestEvaluator.cpp b/packages/stk/stk_unit_tests/stk_expreval/UnitTestEvaluator.cpp index fe69dbaca5c5..16b8daf7da74 100644 --- a/packages/stk/stk_unit_tests/stk_expreval/UnitTestEvaluator.cpp +++ b/packages/stk/stk_unit_tests/stk_expreval/UnitTestEvaluator.cpp @@ -1200,15 +1200,13 @@ TEST(UnitTestEvaluator, testParsedEvalNoUserDefinedFunctions) { stk::expreval::addFunction("my_function", new OneArgFunction()); stk::expreval::Eval eval("my_function(1)"); - EXPECT_ANY_THROW(eval.get_parsed_eval().check_for_errors(false)); - EXPECT_ANY_THROW(eval.get_parsed_eval().check_for_errors(true)); + EXPECT_ANY_THROW(eval.get_parsed_eval()); } TEST(UnitTestEvaluator, testParsedEvalNoRandom) { stk::expreval::Eval eval("rand()"); - EXPECT_ANY_THROW(eval.get_parsed_eval().check_for_errors(false)); - EXPECT_ANY_THROW(eval.get_parsed_eval().check_for_errors(true)); + EXPECT_ANY_THROW(eval.get_parsed_eval()); } @@ -3258,22 +3256,6 @@ TEST(UnitTestEvaluator, testFunction_rand) testRandom("rand()"); } -void Ngp_testRandom(const char * expression) -{ - const int NUM_SAMPLES = 10000; - std::vector results(NUM_SAMPLES); - for (int i = 0; i < NUM_SAMPLES; ++i) { - results[i] = device_evaluate(expression); - } - checkUniformDist(results); -} - -#if !defined(STK_ENABLE_GPU) && !defined(KOKKOS_ENABLE_SYCL) -TEST(UnitTestEvaluator, Ngp_testFunction_rand) -{ - Ngp_testRandom("rand()"); -} -#endif TEST(UnitTestEvaluator, testFunction_srand_repeatability) { @@ -3289,34 +3271,11 @@ TEST(UnitTestEvaluator, testFunction_srand_repeatability) } } -#if !defined(STK_ENABLE_GPU) && !defined(KOKKOS_ENABLE_SYCL) -TEST(UnitTestEvaluator, Ngp_testFunction_srand_repeatability) -{ - std::vector result(10); - device_evaluate("srand(123.)"); - for (unsigned i = 0; i < result.size(); ++i) { - result[i] = device_evaluate("rand()"); - } - - device_evaluate("srand(123.)"); - for (unsigned i = 0; i < result.size(); ++i) { - EXPECT_DOUBLE_EQ(device_evaluate("rand()"), result[i]); - } -} -#endif - TEST(UnitTestEvaluator, testFunction_random) { testRandom("random()"); } -#if !defined(STK_ENABLE_GPU) && !defined(KOKKOS_ENABLE_SYCL) -TEST(UnitTestEvaluator, Ngp_testFunction_random) -{ - Ngp_testRandom("random()"); -} -#endif - TEST(UnitTestEvaluator, testFunction_random1_repeatability) { std::vector result(10); @@ -3331,22 +3290,6 @@ TEST(UnitTestEvaluator, testFunction_random1_repeatability) } } -#if !defined(STK_ENABLE_GPU) && !defined(KOKKOS_ENABLE_SYCL) -TEST(UnitTestEvaluator, Ngp_testFunction_random1_repeatability) -{ - std::vector result(10); - device_evaluate("random(123.)"); - for (unsigned i = 0; i < result.size(); ++i) { - result[i] = device_evaluate("random()"); - } - - device_evaluate("random(123.)"); - for (unsigned i = 0; i < result.size(); ++i) { - EXPECT_DOUBLE_EQ(device_evaluate("random()"), result[i]); - } -} -#endif - TEST(UnitTestEvaluator, testFunction_ts_random_distribution) { const int NX = 50; @@ -3496,11 +3439,5 @@ TEST(UnitTestEvaluator, testFunction_time) EXPECT_NEAR(evaluate("time()"), std::time(nullptr), 1.1); } -#if !defined(STK_ENABLE_GPU) && !defined(KOKKOS_ENABLE_SYCL) -TEST(UnitTestEvaluator, Ngp_testFunction_time) -{ - EXPECT_NEAR(device_evaluate("time()"), std::time(nullptr), 1.1); -} -#endif } // namespace diff --git a/packages/stk/stk_unit_tests/stk_io/UnitTestReadWriteSideSets.cpp b/packages/stk/stk_unit_tests/stk_io/UnitTestReadWriteSideSets.cpp index 21cf866f92ff..5ed5e7632770 100644 --- a/packages/stk/stk_unit_tests/stk_io/UnitTestReadWriteSideSets.cpp +++ b/packages/stk/stk_unit_tests/stk_io/UnitTestReadWriteSideSets.cpp @@ -573,6 +573,32 @@ TEST(StkIo, modify_sideset) } } + +TEST(StkIo, skinned_sideset_from_badly_named_element_block) +{ + stk::ParallelMachine pm = MPI_COMM_WORLD; + + if(stk::parallel_machine_size(pm) > 1) GTEST_SKIP(); + + std::string meshDesc = "textmesh: 0,1,HEX_8,1,2,3,4,5,6,7,8, shell_a" + "|coordinates: 0,0,0, 1,0,0, 1,1,0, 0,1,0, 0,0,1, 1,0,1, 1,1,1, 0,1,1" + "|sideset:name=surface_skinned_shell_a; data=1,1"; + + const unsigned spatialDim = 3; + std::shared_ptr bulk = build_mesh_no_simple_fields(spatialDim, pm, stk::mesh::BulkData::NO_AUTO_AURA); + EXPECT_NO_THROW(stk::io::fill_mesh(meshDesc, *bulk)); + + std::string fileName("missing_skinned_shell.g"); + stk::io::write_mesh(fileName, *bulk); + + std::shared_ptr bulk2 = build_mesh_no_simple_fields(spatialDim, pm, stk::mesh::BulkData::NO_AUTO_AURA); + stk::unit_test_util::sideset::read_exo_file(*bulk2, fileName, stk::unit_test_util::sideset::READ_SERIAL_AND_DECOMPOSE); + stk::unit_test_util::sideset::SideSetIdAndElemIdSidesVector expectedInputSideset{ {1, {{1, 0}}} }; + stk::unit_test_util::sideset::compare_sidesets(fileName, *bulk2, expectedInputSideset); + + unlink(fileName.c_str()); +} + TEST(StkIo, parallel_transform_AA_to_ADA_to_ARA) { stk::ParallelMachine pm = MPI_COMM_WORLD; diff --git a/packages/stk/stk_unit_tests/stk_mesh/UnitTestEntityKeyMapping.cpp b/packages/stk/stk_unit_tests/stk_mesh/UnitTestEntityKeyMapping.cpp new file mode 100644 index 000000000000..b0e7817974b3 --- /dev/null +++ b/packages/stk/stk_unit_tests/stk_mesh/UnitTestEntityKeyMapping.cpp @@ -0,0 +1,53 @@ +// Copyright 2002 - 2008, 2010, 2011 National Technology Engineering +// Solutions of Sandia, LLC (NTESS). Under the terms of Contract +// DE-NA0003525 with NTESS, the U.S. Government retains certain rights +// in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following +// disclaimer in the documentation and/or other materials provided +// with the distribution. +// +// * Neither the name of NTESS nor the names of its contributors +// may be used to endorse or promote products derived from this +// software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// + +#include +#include +#include +#include +#include + +namespace { + +TEST(UnitTestEntityKeyMapping, invalidEntityKey) +{ + stk::mesh::impl::EntityKeyMapping entityKeyMapping; + entityKeyMapping.update_num_ranks(1); + stk::mesh::EntityKey invalidKey; + stk::mesh::Entity entity = entityKeyMapping.get_entity(invalidKey); + EXPECT_EQ(stk::mesh::Entity::InvalidEntity, entity.local_offset()); +} + +}//namespace + diff --git a/packages/stk/stk_unit_tests/stk_mesh/UnitTestFieldBLAS.cpp b/packages/stk/stk_unit_tests/stk_mesh/UnitTestFieldBLAS.cpp index 756105b9976b..b7b9695921b9 100644 --- a/packages/stk/stk_unit_tests/stk_mesh/UnitTestFieldBLAS.cpp +++ b/packages/stk/stk_unit_tests/stk_mesh/UnitTestFieldBLAS.cpp @@ -1939,8 +1939,8 @@ BLASFixture3d::~BLASFixture3d() { delete stkMeshIoBroker; } -template -bool test3dfield(const stk::mesh::Field & field,const A* expected_value,const double tol=1.5e-3) +template +bool test3dfield(const stk::mesh::Field & field,const A* expected_value,const double tol=1.5e-3) { bool result=true; const stk::mesh::BucketVector& buckets_init = field.get_mesh().get_buckets(field.entity_rank(),field.mesh_meta_data().universal_part() & stk::mesh::selectField(field)); @@ -1959,8 +1959,9 @@ bool test3dfield(const stk::mesh::Field & field,const A* return result; } -template -bool test3dfield(const stk::mesh::Field,T1,T2,T3,T4,T5,T6,T7> & field,const std::complex* expected_value,const double tol=1.5e-3) +template +bool test3dfield(const stk::mesh::Field> & field, + const std::complex* expected_value, const double tol=1.5e-3) { bool result=true; const stk::mesh::BucketVector& buckets_init = field.get_mesh().get_buckets(field.entity_rank(),field.mesh_meta_data().universal_part() & stk::mesh::selectField(field)); diff --git a/packages/stk/stk_unit_tests/stk_mesh/UnitTestSelector.cpp b/packages/stk/stk_unit_tests/stk_mesh/UnitTestSelector.cpp index 71aed117df7c..cfe71cca2202 100644 --- a/packages/stk/stk_unit_tests/stk_mesh/UnitTestSelector.cpp +++ b/packages/stk/stk_unit_tests/stk_mesh/UnitTestSelector.cpp @@ -907,7 +907,31 @@ TEST(Selector, get_parts_intersection_ranked) EXPECT_EQ(rankedPart.mesh_meta_data_ordinal(), selector2Parts[0]->mesh_meta_data_ordinal()); } +TEST(Selector, selectUnion_clone_for_different_mesh) +{ + if (stk::parallel_machine_size(MPI_COMM_WORLD) != 1) { GTEST_SKIP(); } + + constexpr unsigned spatialDim = 3; + stk::mesh::MetaData meta1(spatialDim); + std::string part1Name("myPart1"), part2Name("myPart2"); + stk::mesh::PartVector parts1 = { &meta1.declare_part(part1Name), &meta1.declare_part(part2Name)}; + stk::mesh::Selector select1 = stk::mesh::selectUnion(parts1); + + stk::mesh::MetaData meta2; + meta2.declare_part(part1Name); + meta2.declare_part(part2Name); + meta2.initialize(spatialDim); + stk::mesh::Selector select2 = select1.clone_for_different_mesh(meta2); + stk::mesh::PartVector parts2; + select2.get_parts(parts2); + + EXPECT_EQ(2u, parts2.size()); + EXPECT_NE(parts2[0]->mesh_meta_data_ordinal(), parts1[0]->mesh_meta_data_ordinal()); + EXPECT_NE(parts2[1]->mesh_meta_data_ordinal(), parts1[1]->mesh_meta_data_ordinal()); + EXPECT_EQ(parts2[0]->name(), part1Name); + EXPECT_EQ(parts2[1]->name(), part2Name); +} TEST( UnitTestRootTopology, bucketAlsoHasAutoCreatedRootParts ) { diff --git a/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpFieldTestUtils.hpp b/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpFieldTestUtils.hpp index 08329bcac769..17768ef7edfc 100644 --- a/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpFieldTestUtils.hpp +++ b/packages/stk/stk_unit_tests/stk_mesh/ngp/NgpFieldTestUtils.hpp @@ -82,6 +82,50 @@ void check_field_data_on_host(const stk::mesh::BulkData& stkMesh, ); } +inline void set_field_data_on_host(const stk::mesh::BulkData& stkMesh, + const stk::mesh::FieldBase& stkField, + const stk::mesh::Selector& selector, + std::function(const double*)> func) +{ + const stk::mesh::FieldBase* coordField = stkMesh.mesh_meta_data().coordinate_field(); + stk::mesh::for_each_entity_run(stkMesh, stkField.entity_rank(), selector, + [&](const stk::mesh::BulkData& bulk, const stk::mesh::Entity entity) { + double* entityCoords = static_cast(stk::mesh::field_data(*coordField, entity)); + auto expectedValues = func(entityCoords); + const int numComponents = stk::mesh::field_scalars_per_entity(stkField, entity); + double* fieldData = static_cast(stk::mesh::field_data(stkField, entity)); + for(int i=0; i& otherFields, + std::function(const double*)> func) +{ + const stk::mesh::FieldBase* coordField = stkMesh.mesh_meta_data().coordinate_field(); + stk::mesh::for_each_entity_run(stkMesh, stkField.entity_rank(), selector, + [&](const stk::mesh::BulkData& bulk, const stk::mesh::Entity entity) { + double* entityCoords = static_cast(stk::mesh::field_data(*coordField, entity)); + auto expectedValues = func(entityCoords); + + unsigned int numComponents = stk::mesh::field_scalars_per_entity(stkField, entity); + for (const stk::mesh::FieldBase* otherField : otherFields) + { + numComponents = std::min(numComponents, stk::mesh::field_scalars_per_entity(*otherField, entity)); + } + const double* fieldData = reinterpret_cast(stk::mesh::field_data(stkField, entity)); + for(unsigned int i=0; i #include @@ -114,8 +114,93 @@ class NgpFieldBLAS : public stk::unit_test_util::MeshFixture stk::mesh::Field* stkField1 = nullptr; stk::mesh::Field* stkField2 = nullptr; stk::mesh::Field* stkField3 = nullptr; + + stk::mesh::Field* stkNodeField1 = nullptr; + stk::mesh::Field* stkNodeField2 = nullptr; + stk::mesh::Field* stkNodeField3 = nullptr; }; +class NgpFieldBLASNode : public stk::unit_test_util::MeshFixture +{ +public: + NgpFieldBLASNode() + { + setup_three_fields_three_hex_three_block_mesh(); + } + + void setup_three_fields_three_hex_three_block_mesh() + { + const unsigned numStates = 1; + const unsigned bucketCapacity = 2; + setup_empty_mesh(stk::mesh::BulkData::NO_AUTO_AURA, bucketCapacity, bucketCapacity); + + stkField1 = &get_meta().declare_field(stk::topology::NODE_RANK, "variableLengthNodeField1", numStates); + stkField2 = &get_meta().declare_field(stk::topology::NODE_RANK, "variableLengthNodeField2", numStates); + stkField3 = &get_meta().declare_field(stk::topology::NODE_RANK, "variableLengthNodeField3", numStates); + + + stk::mesh::Part& block1 = get_meta().declare_part_with_topology("block_1", stk::topology::HEX_8); + stk::mesh::Part& block2 = get_meta().declare_part_with_topology("block_2", stk::topology::HEX_8); + get_meta().declare_part_with_topology("block_3", stk::topology::HEX_8); + + const std::vector init1(numComponent1, -1); + stk::mesh::put_field_on_mesh(*stkField1, block1, numComponent1, init1.data()); + + const std::vector init2(numComponent2, -2); + stk::mesh::put_field_on_mesh(*stkField1, block2, numComponent2, init2.data()); + + const std::vector init3(numComponent1, -1); + stk::mesh::put_field_on_mesh(*stkField2, block1, numComponent1, init3.data()); + + const std::vector init4(numComponent2, -2); + stk::mesh::put_field_on_mesh(*stkField2, block2, numComponent2, init4.data()); + + const std::vector init5(numComponent1, -1); + stk::mesh::put_field_on_mesh(*stkField3, block1, numComponent1, init5.data()); + + const std::vector init6(numComponent2, -2); + stk::mesh::put_field_on_mesh(*stkField3, block2, numComponent2, init6.data()); + + const std::string meshDesc = "0,1,HEX_8,1,2,3,4,5,6,7,8,block_1\n" + "0,2,HEX_8,9,10,11,12,13,14,15,16,block_2\n" + "0,3,HEX_8,17,18,19,20,21,22,23,24,block_3\n" + "|coordinates:" + // elem 1 + "0,0,0, 1,0,0, 1,1,0, 0,1,0, " + "0,0,1, 1,0,1, 1,1,1, 0,1,1, " + // elem 2 + "2,0,0, 3,0,0, 3,1,0, 2,1,0, " + "2,0,1, 3,0,1, 3,1,1, 2,1,1, " + // elem 3 + "4,0,0, 5,0,0, 5,1,0, 4,1,0, " + "4,0,1, 5,0,1, 5,1,1, 4,1,1"; + stk::unit_test_util::setup_text_mesh(get_bulk(), meshDesc); + + EXPECT_FALSE(stkField1->need_sync_to_host()); + } + + const int numComponent1 = 8; + const int numComponent2 = 3; + stk::mesh::Field* stkField1 = nullptr; + stk::mesh::Field* stkField2 = nullptr; + stk::mesh::Field* stkField3 = nullptr; +}; + +namespace { +std::vector func1(const double* coords) +{ + double v = coords[0] + 2*coords[1] + 3*coords[2]; + return {v, v + 1, v + 2, v + 3, v + 4, v + 5, v + 6, v + 7}; +} + +std::vector func2(const double* coords) +{ + double v = coords[0] + 2*coords[1] + 3*coords[2]; + return {v + 8, v + 9, v + 10, v + 11, v + 12, v + 13, v + 14, v + 15}; +} + +} + #ifdef STK_USE_DEVICE_MESH TEST_F(NgpFieldBLAS, field_fill_device) @@ -406,23 +491,517 @@ TEST_F(NgpFieldBLAS, field_copy_device_with_host_build) #endif -TEST_F(NgpFieldBLAS, field_axpbyz) +#ifdef STK_USE_DEVICE_MESH + +TEST_F(NgpFieldBLASNode, field_axpbyz_device) { - if (get_parallel_size() != 1) GTEST_SKIP(); + if (get_parallel_size() > 2) GTEST_SKIP(); - stk::mesh::field_fill(3.0, *stkField1, stk::ngp::ExecSpace()); - stk::mesh::field_fill(10.0, *stkField2, stk::ngp::ExecSpace()); + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); double alpha = 2.0; double beta = 5.0; + auto f_expected = [&](const double* coords) + { + std::vector vals1 = func1(coords); + std::vector vals2 = func2(coords); + std::vector result(vals1.size()); + for (size_t i=0; i < vals1.size(); ++i) + { + result[i] = alpha * vals1[i] + beta * vals2[i]; + } + + return result; + }; + + + stk::mesh::field_axpbyz(get_bulk(), alpha, *stkField1, beta, *stkField2, *stkField3, selectRule, stk::ngp::ExecSpace()); + EXPECT_TRUE(stkField3->need_sync_to_host()); + + stkField3->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField3, selectRule, {stkField1, stkField2}, f_expected); +} + +TEST_F(NgpFieldBLASNode, field_axpbyz_host) +{ + if (get_parallel_size() > 2) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); + + double alpha = 2.0; + double beta = 5.0; + auto f_expected = [&](const double* coords) + { + std::vector vals1 = func1(coords); + std::vector vals2 = func2(coords); + std::vector result(vals1.size()); + for (size_t i=0; i < vals1.size(); ++i) + { + result[i] = alpha * vals1[i] + beta * vals2[i]; + } + + return result; + }; + + + stk::mesh::field_axpbyz(get_bulk(), alpha, *stkField1, beta, *stkField2, *stkField3, selectRule, stk::ngp::HostExecSpace()); + +#ifdef STK_ENABLE_GPU + EXPECT_TRUE(stkField3->need_sync_to_device()); +#else + EXPECT_TRUE(stkField3->need_sync_to_host()); +#endif + + stkField3->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField3, selectRule, {stkField1, stkField2}, f_expected); +} + +#else + +TEST_F(NgpFieldBLASNode, field_axpbyz_exec_space) +{ + if (get_parallel_size() != 1) GTEST_SKIP(); + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); + + double alpha = 2.0; + double beta = 5.0; + auto f_expected = [&](const double* coords) + { + std::vector vals1 = func1(coords); + std::vector vals2 = func2(coords); + std::vector result(vals1.size()); + for (size_t i=0; i < vals1.size(); ++i) + { + result[i] = alpha * vals1[i] + beta * vals2[i]; + } + + return result; + }; + stk::mesh::field_axpbyz(get_bulk(), alpha, *stkField1, beta, *stkField2, *stkField3, selectRule, stk::ngp::ExecSpace()); + EXPECT_FALSE(stkField3->need_sync_to_host()); stkField3->sync_to_host(); - stk::mesh::Selector selector(*stkField3); - constexpr double expectedValue = 56.0; - ngp_field_test_utils::check_field_data_on_host(get_bulk(), *stkField3, selector, expectedValue); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField3, selectRule, {stkField1, stkField2}, f_expected); +} + +TEST_F(NgpFieldBLASNode, field_axpbyz_no_selector) +{ + if (get_parallel_size() != 1) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); + + double alpha = 2.0; + double beta = 5.0; + auto f_expected = [&](const double* coords) + { + std::vector vals1 = func1(coords); + std::vector vals2 = func2(coords); + std::vector result(vals1.size()); + for (size_t i=0; i < vals1.size(); ++i) + { + result[i] = alpha * vals1[i] + beta * vals2[i]; + } + + return result; + }; + + + stk::mesh::field_axpbyz(get_bulk(), alpha, *stkField1, beta, *stkField2, *stkField3, stk::ngp::ExecSpace()); + EXPECT_FALSE(stkField3->need_sync_to_host()); + + stkField3->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField3, selectRule, {stkField1, stkField2}, f_expected); +} + +#endif + +#ifdef STK_USE_DEVICE_MESH + +TEST_F(NgpFieldBLASNode, field_axpby_device) +{ + if (get_parallel_size() > 2) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); + + double alpha = 2.0; + double beta = 5.0; + auto f_expected = [&](const double* coords) + { + std::vector vals1 = func1(coords); + std::vector vals2 = func2(coords); + std::vector result(vals1.size()); + for (size_t i=0; i < vals1.size(); ++i) + { + result[i] = alpha * vals1[i] + beta * vals2[i]; + } + + return result; + }; + + + stk::mesh::field_axpby(get_bulk(), alpha, *stkField1, beta, *stkField2, selectRule, stk::ngp::ExecSpace()); + EXPECT_TRUE(stkField2->need_sync_to_host()); + + stkField2->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField2, selectRule, {stkField1}, f_expected); +} + +TEST_F(NgpFieldBLASNode, field_axpby_host) +{ + if (get_parallel_size() > 2) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); + + double alpha = 2.0; + double beta = 5.0; + auto f_expected = [&](const double* coords) + { + std::vector vals1 = func1(coords); + std::vector vals2 = func2(coords); + std::vector result(vals1.size()); + for (size_t i=0; i < vals1.size(); ++i) + { + result[i] = alpha * vals1[i] + beta * vals2[i]; + } + + return result; + }; + + + stk::mesh::field_axpby(get_bulk(), alpha, *stkField1, beta, *stkField2, selectRule, stk::ngp::HostExecSpace()); + +#ifdef STK_ENABLE_GPU + EXPECT_TRUE(stkField2->need_sync_to_device()); +#else + EXPECT_TRUE(stkField2->need_sync_to_host()); +#endif + + stkField2->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField2, selectRule, {stkField1}, f_expected); +} + +#else + +TEST_F(NgpFieldBLASNode, field_axpby_exec_space) +{ + if (get_parallel_size() != 1) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); + + double alpha = 2.0; + double beta = 5.0; + auto f_expected = [&](const double* coords) + { + std::vector vals1 = func1(coords); + std::vector vals2 = func2(coords); + std::vector result(vals1.size()); + for (size_t i=0; i < vals1.size(); ++i) + { + result[i] = alpha * vals1[i] + beta * vals2[i]; + } + + return result; + }; + + + stk::mesh::field_axpby(get_bulk(), alpha, *stkField1, beta, *stkField2, selectRule, stk::ngp::ExecSpace()); + EXPECT_FALSE(stkField2->need_sync_to_host()); + + stkField2->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField2, selectRule, {stkField1}, f_expected); +} + +TEST_F(NgpFieldBLASNode, field_axpby_no_selector) +{ + if (get_parallel_size() != 1) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); + + double alpha = 2.0; + double beta = 5.0; + auto f_expected = [&](const double* coords) + { + std::vector vals1 = func1(coords); + std::vector vals2 = func2(coords); + std::vector result(vals1.size()); + for (size_t i=0; i < vals1.size(); ++i) + { + result[i] = alpha * vals1[i] + beta * vals2[i]; + } + + return result; + }; + + + stk::mesh::field_axpby(get_bulk(), alpha, *stkField1, beta, *stkField2, stk::ngp::ExecSpace()); + EXPECT_FALSE(stkField2->need_sync_to_host()); + + stkField2->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField2, selectRule, {stkField1}, f_expected); +} + + +#endif + +#ifdef STK_USE_DEVICE_MESH + +TEST_F(NgpFieldBLASNode, field_axpy_device) +{ + if (get_parallel_size() > 2) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + + double alpha = 2.0; + auto f_expected = [&](const double* coords) + { + std::vector result = func1(coords); + for (size_t i=0; i < result.size(); ++i) + { + result[i] = alpha * result[i]; + } + + return result; + }; + + + stk::mesh::field_axpy(get_bulk(), alpha, *stkField1, *stkField2, selectRule, stk::ngp::ExecSpace()); + EXPECT_TRUE(stkField2->need_sync_to_host()); + + stkField2->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField2, selectRule, {stkField1}, f_expected); +} + +TEST_F(NgpFieldBLASNode, field_axpy_host) +{ + if (get_parallel_size() > 2) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + + double alpha = 2.0; + auto f_expected = [&](const double* coords) + { + std::vector result = func1(coords); + for (size_t i=0; i < result.size(); ++i) + { + result[i] = alpha * result[i]; + } + + return result; + }; + + + stk::mesh::field_axpy(get_bulk(), alpha, *stkField1, *stkField2, selectRule, stk::ngp::HostExecSpace()); + +#ifdef STK_ENABLE_GPU + EXPECT_TRUE(stkField2->need_sync_to_device()); +#else + EXPECT_TRUE(stkField2->need_sync_to_host()); +#endif + + stkField2->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField2, selectRule, {stkField1}, f_expected); +} + +#else + +TEST_F(NgpFieldBLASNode, field_axpy_exec_space) +{ + if (get_parallel_size() != 1) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + + double alpha = 2.0; + auto f_expected = [&](const double* coords) + { + std::vector result = func1(coords); + for (size_t i=0; i < result.size(); ++i) + { + result[i] = alpha * result[i]; + } + + return result; + }; + + stk::mesh::field_axpy(get_bulk(), alpha, *stkField1, *stkField2, selectRule, stk::ngp::ExecSpace()); + EXPECT_FALSE(stkField2->need_sync_to_host()); + + stkField2->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField2, selectRule, {stkField1}, f_expected); +} + +TEST_F(NgpFieldBLASNode, field_axpy_no_selector) +{ + if (get_parallel_size() != 1) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + + double alpha = 2.0; + auto f_expected = [&](const double* coords) + { + std::vector result = func1(coords); + for (size_t i=0; i < result.size(); ++i) + { + result[i] = alpha * result[i]; + } + + return result; + }; + + stk::mesh::field_axpy(get_bulk(), alpha, *stkField1, *stkField2, stk::ngp::ExecSpace()); + EXPECT_FALSE(stkField2->need_sync_to_host()); + + stkField2->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField2, selectRule, {stkField1}, f_expected); +} + +#endif + +#ifdef STK_USE_DEVICE_MESH + +TEST_F(NgpFieldBLASNode, field_product_device) +{ + if (get_parallel_size() > 2) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); + + + auto f_expected = [&](const double* coords) + { + std::vector vals1 = func1(coords); + std::vector vals2 = func2(coords); + std::vector result(vals1.size()); + for (size_t i=0; i < vals1.size(); ++i) + { + result[i] = vals1[i] * vals2[i]; + } + + return result; + }; + + + stk::mesh::field_product(get_bulk(), *stkField1, *stkField2, *stkField3, selectRule, stk::ngp::ExecSpace()); + EXPECT_TRUE(stkField3->need_sync_to_host()); + + stkField3->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField3, selectRule, {stkField1, stkField2}, f_expected); +} + +TEST_F(NgpFieldBLASNode, field_product_host) +{ + if (get_parallel_size() > 2) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); + + auto f_expected = [&](const double* coords) + { + std::vector vals1 = func1(coords); + std::vector vals2 = func2(coords); + std::vector result(vals1.size()); + for (size_t i=0; i < vals1.size(); ++i) + { + result[i] = vals1[i] * vals2[i]; + } + + return result; + }; + + + stk::mesh::field_product(get_bulk(), *stkField1, *stkField2, *stkField3, selectRule, stk::ngp::HostExecSpace()); + +#ifdef STK_ENABLE_GPU + EXPECT_TRUE(stkField3->need_sync_to_device()); +#else + EXPECT_TRUE(stkField3->need_sync_to_host()); +#endif + + stkField3->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField3, selectRule, {stkField1, stkField2}, f_expected); +} + +#else + +TEST_F(NgpFieldBLASNode, field_product_exec_space) +{ + if (get_parallel_size() != 1) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); + + auto f_expected = [&](const double* coords) + { + std::vector vals1 = func1(coords); + std::vector vals2 = func2(coords); + std::vector result(vals1.size()); + for (size_t i=0; i < vals1.size(); ++i) + { + result[i] = vals1[i] * vals2[i]; + } + + return result; + }; + + stk::mesh::field_product(get_bulk(), *stkField1, *stkField2, *stkField3, selectRule, stk::ngp::ExecSpace()); + EXPECT_FALSE(stkField3->need_sync_to_host()); + + stkField3->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField3, selectRule, {stkField1, stkField2}, f_expected); +} + +TEST_F(NgpFieldBLASNode, field_product_no_selector) +{ + if (get_parallel_size() != 1) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); + + + auto f_expected = [&](const double* coords) + { + std::vector vals1 = func1(coords); + std::vector vals2 = func2(coords); + std::vector result(vals1.size()); + for (size_t i=0; i < vals1.size(); ++i) + { + result[i] = vals1[i] * vals2[i]; + } + + return result; + }; + + stk::mesh::field_product(get_bulk(), *stkField1, *stkField2, *stkField3, stk::ngp::ExecSpace()); + EXPECT_FALSE(stkField3->need_sync_to_host()); + + stkField3->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField3, selectRule, {stkField1, stkField2}, f_expected); } TEST_F(NgpFieldBLAS, field_axpby) @@ -443,6 +1022,251 @@ TEST_F(NgpFieldBLAS, field_axpby) constexpr double expectedValue = 56.0; ngp_field_test_utils::check_field_data_on_host(get_bulk(), *stkField2, selector, expectedValue); } +#endif + +#ifdef STK_USE_DEVICE_MESH + +TEST_F(NgpFieldBLASNode, field_scale_device) +{ + if (get_parallel_size() > 2) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + + double alpha = 2.0; + auto f_expected = [&](const double* coords) + { + std::vector result = func1(coords); + for (size_t i=0; i < result.size(); ++i) + { + result[i] = alpha * result[i]; + } + + return result; + }; + + + stk::mesh::field_scale(get_bulk(), alpha, *stkField1, selectRule, stk::ngp::ExecSpace()); + EXPECT_TRUE(stkField1->need_sync_to_host()); + + stkField1->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField1, selectRule, {}, f_expected); +} + +TEST_F(NgpFieldBLASNode, field_scale_host) +{ + if (get_parallel_size() > 2) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + + double alpha = 2.0; + auto f_expected = [&](const double* coords) + { + std::vector result = func1(coords); + for (size_t i=0; i < result.size(); ++i) + { + result[i] = alpha * result[i]; + } + + return result; + }; + + + stk::mesh::field_scale(get_bulk(), alpha, *stkField1, selectRule, stk::ngp::HostExecSpace()); + +#ifdef STK_ENABLE_GPU + EXPECT_TRUE(stkField1->need_sync_to_device()); +#else + EXPECT_TRUE(stkField1->need_sync_to_host()); +#endif + + stkField1->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField1, selectRule, {}, f_expected); +} + +#else + +TEST_F(NgpFieldBLASNode, field_scale_exec_space) +{ + if (get_parallel_size() != 1) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + + double alpha = 2.0; + auto f_expected = [&](const double* coords) + { + std::vector result = func1(coords); + for (size_t i=0; i < result.size(); ++i) + { + result[i] = alpha * result[i]; + } + + return result; + }; + + stk::mesh::field_scale(get_bulk(), alpha, *stkField1, selectRule, stk::ngp::ExecSpace()); + EXPECT_FALSE(stkField1->need_sync_to_host()); + + stkField1->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField1, selectRule, {}, f_expected); +} + +TEST_F(NgpFieldBLASNode, field_scale_no_selector) +{ + if (get_parallel_size() != 1) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + + double alpha = 2.0; + auto f_expected = [&](const double* coords) + { + std::vector result = func1(coords); + for (size_t i=0; i < result.size(); ++i) + { + result[i] = alpha * result[i]; + } + + return result; + }; + + stk::mesh::field_scale(get_bulk(), alpha, *stkField1, stk::ngp::ExecSpace()); + EXPECT_FALSE(stkField1->need_sync_to_host()); + + stkField1->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField1, selectRule, {}, f_expected); +} + +#endif + +#ifdef STK_USE_DEVICE_MESH + +TEST_F(NgpFieldBLASNode, field_swap_device) +{ + if (get_parallel_size() > 2) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); + + auto f1_expected = [&](const double* coords) + { + return func2(coords); + }; + + auto f2_expected = [&](const double* coords) + { + return func1(coords); + }; + + stk::mesh::field_swap(get_bulk(), *stkField1, *stkField2, selectRule, stk::ngp::ExecSpace()); + EXPECT_TRUE(stkField1->need_sync_to_host()); + EXPECT_TRUE(stkField2->need_sync_to_host()); + + stkField1->sync_to_host(); + stkField2->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField1, selectRule, {stkField2}, f1_expected); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField2, selectRule, {stkField1}, f2_expected); +} + +TEST_F(NgpFieldBLASNode, field_swap_host) +{ + if (get_parallel_size() > 2) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); + + auto f1_expected = [&](const double* coords) + { + return func2(coords); + }; + + auto f2_expected = [&](const double* coords) + { + return func1(coords); + }; + + + stk::mesh::field_swap(get_bulk(), *stkField1, *stkField2, selectRule, stk::ngp::HostExecSpace()); + +#ifdef STK_ENABLE_GPU + EXPECT_TRUE(stkField1->need_sync_to_device()); + EXPECT_TRUE(stkField2->need_sync_to_device()); +#else + EXPECT_TRUE(stkField1->need_sync_to_host()); + EXPECT_TRUE(stkField2->need_sync_to_host()); +#endif + + stkField1->sync_to_host(); + stkField2->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField1, selectRule, {stkField2}, f1_expected); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField2, selectRule, {stkField1}, f2_expected); +} + +#else + +TEST_F(NgpFieldBLASNode, field_swap_exec_space) +{ + if (get_parallel_size() != 1) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); + + auto f1_expected = [&](const double* coords) + { + return func2(coords); + }; + + auto f2_expected = [&](const double* coords) + { + return func1(coords); + }; + + stk::mesh::field_swap(get_bulk(), *stkField1, *stkField2, selectRule, stk::ngp::ExecSpace()); + EXPECT_FALSE(stkField1->need_sync_to_host()); + EXPECT_FALSE(stkField2->need_sync_to_host()); + + stkField1->sync_to_host(); + stkField2->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField1, selectRule, {stkField2}, f1_expected); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField2, selectRule, {stkField1}, f2_expected); +} + +TEST_F(NgpFieldBLASNode, field_swap_no_selector) +{ + if (get_parallel_size() != 1) GTEST_SKIP(); + + stk::mesh::Selector selectRule(*stkField1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField1, selectRule, &func1); + ngp_field_test_utils::set_field_data_on_host(get_bulk(), *stkField2, selectRule, &func2); + + auto f1_expected = [&](const double* coords) + { + return func2(coords); + }; + + auto f2_expected = [&](const double* coords) + { + return func1(coords); + }; + + stk::mesh::field_swap(get_bulk(), *stkField1, *stkField2, stk::ngp::ExecSpace()); + EXPECT_FALSE(stkField1->need_sync_to_host()); + EXPECT_FALSE(stkField2->need_sync_to_host()); + + stkField1->sync_to_host(); + stkField2->sync_to_host(); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField1, selectRule, {stkField2}, f1_expected); + ngp_field_test_utils::check_field_data_on_host_func(get_bulk(), *stkField2, selectRule, {stkField1}, f2_expected); + +} + +#endif + } diff --git a/packages/stk/stk_unit_tests/stk_middle_mesh/test_active_vert_container.cpp b/packages/stk/stk_unit_tests/stk_middle_mesh/test_active_vert_container.cpp index 2c6d8f38c6b5..17e390e3710c 100644 --- a/packages/stk/stk_unit_tests/stk_middle_mesh/test_active_vert_container.cpp +++ b/packages/stk/stk_unit_tests/stk_middle_mesh/test_active_vert_container.cpp @@ -1,4 +1,5 @@ #include "gtest/gtest.h" +#include "stk_util/util/ReportHandler.hpp" #include "stk_middle_mesh/active_vert_container.hpp" #include "stk_middle_mesh/boundary_fixture.hpp" #include "stk_middle_mesh/create_mesh.hpp" @@ -21,7 +22,7 @@ class AlwaysTrue ActiveVertData& get_closest_patch(stk::middle_mesh::mesh::impl::ActiveVertContainer& container, const Point& pt) { double minDist = std::numeric_limits::max(); - ActiveVertData* minPatch; + ActiveVertData* minPatch = nullptr; for (auto& patch : container.get_active_verts()) { auto disp = patch.get_current_vert()->get_point_orig(0) - pt; @@ -33,6 +34,7 @@ ActiveVertData& get_closest_patch(stk::middle_mesh::mesh::impl::ActiveVertContai } } + STK_ThrowRequireMsg(minPatch != nullptr, "get_closest_patch failed to find patch"); return *minPatch; } @@ -293,4 +295,4 @@ TEST(ActiveVertContainer, VertNotConnectedByLocalElement) EXPECT_EQ(activeVert.get_num_local_verts(), 7); EXPECT_EQ(activeVert.get_num_verts(), 9); } -} \ No newline at end of file +} diff --git a/packages/stk/stk_unit_tests/stk_search/UnitTestCoarseSearch.cpp b/packages/stk/stk_unit_tests/stk_search/UnitTestCoarseSearch.cpp index 15f1ba9c0bf2..290868f7b1f1 100644 --- a/packages/stk/stk_unit_tests/stk_search/UnitTestCoarseSearch.cpp +++ b/packages/stk/stk_unit_tests/stk_search/UnitTestCoarseSearch.cpp @@ -169,7 +169,13 @@ void test_coarse_search_for_algorithm_with_views(stk::search::SearchMethod algor SearchResultsViewType searchResults; - stk::search::coarse_search(domain, range, algorithm, comm, searchResults); + auto execSpace = ExecSpace{}; + bool enforceSearchResultSymmetry = true; + bool autoSwapDomainAndRange = true; + bool sortSearchResults = true; + + stk::search::coarse_search(domain, range, algorithm, comm, searchResults, execSpace, + enforceSearchResultSymmetry, autoSwapDomainAndRange, sortSearchResults); auto searchResultsHost = Kokkos::create_mirror_view_and_copy(HostSpace{}, searchResults); @@ -208,7 +214,11 @@ void test_coarse_search_for_algorithm(stk::search::SearchMethod algorithm, MPI_C SearchResults searchResults; - stk::search::coarse_search(domain, range, algorithm, comm, searchResults); + bool enforceSearchResultSymmetry = true; + bool autoSwapDomainAndRange = true; + bool sortSearchResults = true; + + stk::search::coarse_search(domain, range, algorithm, comm, searchResults, enforceSearchResultSymmetry, autoSwapDomainAndRange, sortSearchResults); expect_search_results(num_procs, proc_id, searchResults); } @@ -360,9 +370,8 @@ void host_local_test_coarse_search_for_algorithm(stk::search::SearchMethod algor LocalSearchResults intersections; - stk::search::local_coarse_search(domain, range, algorithm, intersections); - std::sort(intersections.begin(), intersections.end()); - + bool sortSearchResults = true; + stk::search::local_coarse_search(domain, range, algorithm, intersections, sortSearchResults); local_expect_search_results(intersections); } @@ -387,25 +396,32 @@ void device_local_test_coarse_search_for_algorithm(stk::search::SearchMethod alg auto intersections = Kokkos::View("intersections", 0); - stk::search::local_coarse_search(domain, range, algorithm, intersections); - Kokkos::sort(intersections); + auto execSpace = stk::ngp::ExecSpace{}; + bool sortSearchResults = true; + stk::search::local_coarse_search(domain, range, algorithm, intersections, execSpace, sortSearchResults); Kokkos::View::HostMirror hostIntersections = Kokkos::create_mirror_view(intersections); Kokkos::deep_copy(hostIntersections, intersections); + Kokkos::sort(hostIntersections); + local_expect_search_results(hostIntersections); } TEST(CoarseSearchCorrectness, Ngp_Local_CoarseSearchDoubleBoxes_MORTON_LBVH) { host_local_test_coarse_search_for_algorithm(stk::search::MORTON_LBVH); +std::cout<<"finished local test"<(stk::search::MORTON_LBVH); +std::cout<<"finished device test"<(stk::search::MORTON_LBVH); +std::cout<<"finished local test"<(stk::search::MORTON_LBVH); +std::cout<<"finished device test"<("intersections", 0); - stk::search::local_coarse_search(domain, range, algorithm, intersections); - Kokkos::sort(intersections); + auto execSpace = ExecSpace{}; + bool sortSearchResults = true; + stk::search::local_coarse_search(domain, range, algorithm, intersections, execSpace, sortSearchResults); auto hostIntersections = Kokkos::create_mirror_view(HostSpace{}, intersections); Kokkos::deep_copy(hostIntersections, intersections); @@ -580,11 +597,8 @@ void test_coarse_search_determine_domain_and_range_communicate_on(stk::search::S SearchResults searchResultsDetermineOn; SearchResults searchResultsDetermineOff; - stk::search::coarse_search(local_domain, local_range, algorithm, comm, searchResultsDetermineOn, true, true); - stk::search::coarse_search(local_domain, local_range, algorithm, comm, searchResultsDetermineOff, true, false); - - std::sort(searchResultsDetermineOn.begin(), searchResultsDetermineOn.end()); - std::sort(searchResultsDetermineOff.begin(), searchResultsDetermineOff.end()); + stk::search::coarse_search(local_domain, local_range, algorithm, comm, searchResultsDetermineOn, true, true, true); + stk::search::coarse_search(local_domain, local_range, algorithm, comm, searchResultsDetermineOff, true, false, true); EXPECT_EQ(searchResultsDetermineOn, searchResultsDetermineOff); } @@ -687,7 +701,12 @@ void test_ident_proc_with_search_with_views(stk::search::SearchMethod searchMeth SearchResultsViewType searchResults("", 3); - coarse_search(boxes, boxes, searchMethod, comm, searchResults); + auto execSpace = ExecSpace{}; + bool enforceSearchResultSymmetry = true; + bool autoSwapDomainAndRange = true; + bool sortSearchResults = true; + coarse_search(boxes, boxes, searchMethod, comm, searchResults, execSpace, + enforceSearchResultSymmetry, autoSwapDomainAndRange, sortSearchResults); SearchResultsViewType goldResults("", 3); @@ -706,7 +725,6 @@ void test_ident_proc_with_search_with_views(stk::search::SearchMethod searchMeth ASSERT_EQ(3u, searchResults.extent(0)); } - Kokkos::sort(searchResults); Kokkos::sort(goldResults); for (size_t i = 0; i < goldResults.extent(0); i++) { @@ -739,7 +757,11 @@ void test_ident_proc_with_search(stk::search::SearchMethod searchMethod) SearchResults searchResults; - coarse_search(boxes, boxes, searchMethod, comm, searchResults); + bool enforceSearchResultSymmetry = true; + bool autoSwapDomainAndRange = true; + bool sortSearchResults = true; + coarse_search(boxes, boxes, searchMethod, comm, searchResults, + enforceSearchResultSymmetry, autoSwapDomainAndRange, sortSearchResults); SearchResults goldResults; diff --git a/packages/stk/stk_unit_tests/stk_search/UnitTestCoarseSearchTwoBox.cpp b/packages/stk/stk_unit_tests/stk_search/UnitTestCoarseSearchTwoBox.cpp index 9f785c2d470c..cc8a6635672e 100644 --- a/packages/stk/stk_unit_tests/stk_search/UnitTestCoarseSearchTwoBox.cpp +++ b/packages/stk/stk_unit_tests/stk_search/UnitTestCoarseSearchTwoBox.cpp @@ -251,9 +251,9 @@ void host_local_runTwoBoxTest(stk::search::SearchMethod searchMethod, distanceBetweenBoxCenters, 0, 0, boxSize / 2, 2))); LocalSearchResults intersections; - - stk::search::local_coarse_search(domain, range, searchMethod, intersections); - std::sort(intersections.begin(), intersections.end()); + + bool sortSearchResults = true; + stk::search::local_coarse_search(domain, range, searchMethod, intersections, sortSearchResults); ASSERT_EQ(intersections.size(), expectedNumOverlap); @@ -281,10 +281,11 @@ void device_local_runTwoBoxTest(stk::search::SearchMethod searchMethod, const do auto intersections = Kokkos::View("intersections", 0); - stk::search::local_coarse_search(domain, range, searchMethod, intersections); + auto execSpace = stk::ngp::ExecSpace{}; + bool sortSearchResults = true; + stk::search::local_coarse_search(domain, range, searchMethod, intersections, execSpace, sortSearchResults); Kokkos::View::HostMirror hostIntersections = Kokkos::create_mirror_view(intersections); - Kokkos::sort(intersections); Kokkos::deep_copy(hostIntersections, intersections); ASSERT_EQ(intersections.extent(0), expectedNumOverlap); diff --git a/packages/stk/stk_unit_tests/stk_search/UnitTestCoarseSearchTwoSpheres.cpp b/packages/stk/stk_unit_tests/stk_search/UnitTestCoarseSearchTwoSpheres.cpp index f26d822a29e6..cd3431a63ccd 100644 --- a/packages/stk/stk_unit_tests/stk_search/UnitTestCoarseSearchTwoSpheres.cpp +++ b/packages/stk/stk_unit_tests/stk_search/UnitTestCoarseSearchTwoSpheres.cpp @@ -254,8 +254,8 @@ void host_local_runTwoSpheresTest(stk::search::SearchMethod searchMethod, LocalSearchResults intersections; - stk::search::local_coarse_search(domain, range, searchMethod, intersections); - std::sort(intersections.begin(), intersections.end()); + bool sortSearchResults = true; + stk::search::local_coarse_search(domain, range, searchMethod, intersections, sortSearchResults); ASSERT_EQ(intersections.size(), expectedNumOverlap); @@ -284,10 +284,11 @@ void device_local_runTwoSpheresTest(stk::search::SearchMethod searchMethod, cons auto intersections = Kokkos::View("intersections", 0); - stk::search::local_coarse_search(domain, range, searchMethod, intersections); + auto execSpace = stk::ngp::ExecSpace{}; + bool sortSearchResults = true; + stk::search::local_coarse_search(domain, range, searchMethod, intersections, execSpace, sortSearchResults); Kokkos::View::HostMirror hostIntersections = Kokkos::create_mirror_view(intersections); - Kokkos::sort(intersections); Kokkos::deep_copy(hostIntersections, intersections); ASSERT_EQ(intersections.extent(0), expectedNumOverlap); diff --git a/packages/stk/stk_unit_tests/stk_search_util/CMakeLists.txt b/packages/stk/stk_unit_tests/stk_search_util/CMakeLists.txt index c3830fe00aa5..e2c6e42acf4f 100644 --- a/packages/stk/stk_unit_tests/stk_search_util/CMakeLists.txt +++ b/packages/stk/stk_unit_tests/stk_search_util/CMakeLists.txt @@ -5,7 +5,7 @@ if(HAVE_STK_Trilinos) TRIBITS_ADD_EXECUTABLE_AND_TEST(stk_search_util_utest SOURCES ${SOURCES} TESTONLYLIBS stk_unit_main - ARGS "--gtest_filter=-SearchUtil.Intrepid2*" + ARGS "" COMM serial mpi NUM_MPI_PROCS 1 ) diff --git a/packages/stk/stk_util/stk_util/Version.hpp b/packages/stk/stk_util/stk_util/Version.hpp index fd0d003e2e0a..297266908285 100644 --- a/packages/stk/stk_util/stk_util/Version.hpp +++ b/packages/stk/stk_util/stk_util/Version.hpp @@ -44,7 +44,7 @@ //See the file CHANGELOG.md for a listing that shows the //correspondence between version numbers and API changes. -#define STK_VERSION 5210401 +#define STK_VERSION 5210500 namespace stk diff --git a/packages/stk/stk_util/stk_util/registry/ProductRegistry.cpp b/packages/stk/stk_util/stk_util/registry/ProductRegistry.cpp index 33d28485af59..32c5615621f3 100644 --- a/packages/stk/stk_util/stk_util/registry/ProductRegistry.cpp +++ b/packages/stk/stk_util/stk_util/registry/ProductRegistry.cpp @@ -42,7 +42,7 @@ //In Sierra, STK_VERSION_STRING is provided on the compile line by bake. //For Trilinos stk snapshots, the following macro definition gets populated with //the real version string by the trilinos_snapshot.sh script. -#define STK_VERSION_STRING "5.21.4-85-gd1ee7818" +#define STK_VERSION_STRING "5.21.5-241-g354b4bbc" #endif namespace stk {