diff --git a/CMakeLists.txt b/CMakeLists.txt index 2ebb5f77c..38eaa5b1e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,7 +1,7 @@ # CMakeLists to build the Spheral library. cmake_minimum_required(VERSION 3.10) -project(spheral LANGUAGES CXX Fortran) +project(spheral LANGUAGES C CXX Fortran) set(SPHERAL_ROOT_DIR ${CMAKE_CURRENT_SOURCE_DIR}) set(SPHERAL_TEST_INSTALL_PREFIX ${CMAKE_INSTALL_PREFIX}) diff --git a/scripts/spack/packages/spheral/package.py b/scripts/spack/packages/spheral/package.py index 11f3881b0..5a75f5386 100644 --- a/scripts/spack/packages/spheral/package.py +++ b/scripts/spack/packages/spheral/package.py @@ -52,8 +52,8 @@ class Spheral(CachedCMakePackage, CudaPackage): depends_on('conduit@0.8.2 +shared +mpi +hdf5 -test ~parmetis', type='build', when='+mpi') depends_on('conduit@0.8.2 +shared ~mpi +hdf5 -test ~parmetis', type='build', when='~mpi') - depends_on('axom@0.5.0 ~shared +mpi +hdf5 -lua -examples -python -fortran -umpire -raja', type='build', when='+mpi') - depends_on('axom@0.5.0 ~shared ~mpi +hdf5 -lua -examples -python -fortran -umpire -raja', type='build', when='~mpi') + depends_on('axom@0.7.0 ~shared +mpi +hdf5 -lua -examples -python -fortran -umpire -raja', type='build', when='+mpi') + depends_on('axom@0.7.0 ~shared ~mpi +hdf5 -lua -examples -python -fortran -umpire -raja', type='build', when='~mpi') depends_on('caliper@2.8.0 ~shared ~adiak ~libdw ~papi ~libunwind +pic', type='build') diff --git a/src/KernelIntegrator/BilinearIndex.cc b/src/KernelIntegrator/BilinearIndex.cc index e28f32b41..c4373d6e6 100644 --- a/src/KernelIntegrator/BilinearIndex.cc +++ b/src/KernelIntegrator/BilinearIndex.cc @@ -36,7 +36,7 @@ computeConnectivity(const DataBase& dataBase) { mConnectivityInitialized = true; } - const auto numNodeLists = dataBase.numNodeLists(); + const auto numNodeLists = dataBase.numFluidNodeLists(); const auto& connectivityMap = dataBase.connectivityMap(); // Clear the index arrays and make a guess about how many indices there will be @@ -128,7 +128,7 @@ computeSurfaceIndexing(const DataBase& dataBase, } // Get DataBase values - const auto numNodeLists = dataBase.numNodeLists(); + const auto numNodeLists = dataBase.numFluidNodeLists(); const auto& connectivityMap = dataBase.connectivityMap(); // Get State values diff --git a/src/KernelIntegrator/FlatConnectivity.cc b/src/KernelIntegrator/FlatConnectivity.cc index 28d22ec67..4645a1476 100644 --- a/src/KernelIntegrator/FlatConnectivity.cc +++ b/src/KernelIntegrator/FlatConnectivity.cc @@ -22,6 +22,23 @@ namespace Spheral { +namespace { // anonymous +template +int +numFluidNeighbors(const std::vector>& connectivity, + const DataBase& dataBase) +{ + auto numNeighbors = 0; + auto nodeListj = 0; + for (auto nodeListItrj = dataBase.fluidNodeListBegin(); + nodeListItrj != dataBase.fluidNodeListEnd(); + ++nodeListItrj, ++nodeListj) { + numNeighbors += connectivity[nodeListj].size(); + } + return numNeighbors; +} +} // end namespace anonymous + //------------------------------------------------------------------------------ // Constructor //------------------------------------------------------------------------------ @@ -46,10 +63,12 @@ template void FlatConnectivity:: computeIndices(const DataBase& dataBase) { + VERIFY(fluidNodeListsFirst(dataBase)); + // Get information from DataBase - const auto numNodeListsDB = dataBase.numNodeLists(); - const auto numNodesDB = dataBase.numNodes(); - const auto numInternalNodesDB = dataBase.numInternalNodes(); + const auto numNodeListsDB = dataBase.numFluidNodeLists(); + const auto numNodesDB = dataBase.numFluidNodes(); + const auto numInternalNodesDB = dataBase.numFluidInternalNodes(); const auto& connectivity = dataBase.connectivityMap(); const auto requireGhostConnectivity = connectivity.buildGhostConnectivity(); @@ -109,10 +128,10 @@ computeIndices(const DataBase& dataBase) { const auto numNodesi = requireGhostConnectivity ? (*nodeListItri)->numNodes() : (*nodeListItri)->numInternalNodes(); for (auto nodei = 0u; nodei < numNodesi; ++nodei) { // Get data from the connectivity map - const auto numNeighborsi = connectivity.numNeighborsForNode(nodeListi, nodei); const auto connectivityi = connectivity.connectivityForNode(nodeListi, nodei); const auto locali = mNodeToLocalIndex[nodeListi][nodei]; - + const auto numNeighborsi = numFluidNeighbors(connectivityi, dataBase); + // Resize the arrays CHECK(locali < mNumConnectivityNodes); mNumNeighbors[locali] = numNeighborsi + 1; @@ -188,9 +207,9 @@ computeOverlapIndices(const DataBase& dataBase) { VERIFY(mIndexingInitialized); // Get information from DataBase - const auto numNodeListsDB = dataBase.numNodeLists(); - const auto numNodesDB = dataBase.numNodes(); - const auto numInternalNodesDB = dataBase.numInternalNodes(); + const auto numNodeListsDB = dataBase.numFluidNodeLists(); + const auto numNodesDB = dataBase.numFluidNodes(); + const auto numInternalNodesDB = dataBase.numFluidInternalNodes(); const auto& connectivity = dataBase.connectivityMap(); const auto requireGhostConnectivity = connectivity.buildGhostConnectivity(); VERIFY(connectivity.buildOverlapConnectivity()); @@ -213,9 +232,9 @@ computeOverlapIndices(const DataBase& dataBase) { const auto numNodesi = requireGhostConnectivity ? (*nodeListItri)->numNodes() : (*nodeListItri)->numInternalNodes(); for (auto nodei = 0u; nodei < numNodesi; ++nodei) { // Get data from the connectivity map - const auto numNeighborsi = connectivity.numOverlapNeighborsForNode(nodeListi, nodei); const auto connectivityi = connectivity.overlapConnectivityForNode(nodeListi, nodei); const auto locali = mNodeToLocalIndex[nodeListi][nodei]; + const auto numNeighborsi = numFluidNeighbors(connectivityi, dataBase); // Resize the arrays mNumOverlapNeighbors[locali] = numNeighborsi + 1; @@ -299,10 +318,10 @@ computeGlobalIndices(const DataBase& dataBase, VERIFY(mIndexingInitialized); // Get information from DataBase - const auto numNodeListsDB = dataBase.numNodeLists(); - const auto numNodesDB = dataBase.numNodes(); - const auto numInternalNodesDB = dataBase.numInternalNodes(); - const auto numGlobalNodesDB = dataBase.globalNumInternalNodes(); + const auto numNodeListsDB = dataBase.numFluidNodeLists(); + const auto numNodesDB = dataBase.numFluidNodes(); + const auto numInternalNodesDB = dataBase.numFluidInternalNodes(); + const auto numGlobalNodesDB = dataBase.globalNumFluidInternalNodes(); // Make sure number of nodes has not changed since computing indices VERIFY(numNodesDB == mNumLocalNodes); @@ -387,9 +406,9 @@ computeSurfaceIndices(const DataBase& dataBase, VERIFY(mGhostIndexingInitialized); // Could consider editing to not require this // Get information from the DataBase and State - const auto numNodeListsDB = dataBase.numNodeLists(); - const auto numNodesDB = dataBase.numNodes(); - const auto numInternalNodesDB = dataBase.numInternalNodes(); + const auto numNodeListsDB = dataBase.numFluidNodeLists(); + const auto numNodesDB = dataBase.numFluidNodes(); + const auto numInternalNodesDB = dataBase.numFluidInternalNodes(); const auto& connectivity = dataBase.connectivityMap(); const auto cells = state.fields(HydroFieldNames::cells, FacetedVolume()); const auto cellFaceFlags = state.fields(HydroFieldNames::cellFaceFlags, std::vector()); @@ -535,9 +554,9 @@ computeBoundaryInformation(const DataBase& dataBase, VERIFY(mIndexingInitialized); // Get information from the dataBase - const auto numNodeListsDB = dataBase.numNodeLists(); - const auto numNodesDB = dataBase.numNodes(); - const auto numInternalNodesDB = dataBase.numInternalNodes(); + const auto numNodeListsDB = dataBase.numFluidNodeLists(); + const auto numNodesDB = dataBase.numFluidNodes(); + const auto numInternalNodesDB = dataBase.numFluidInternalNodes(); // Make sure the sizes haven't changed since the indexing was initialized VERIFY(numNodesDB == mNumLocalNodes); @@ -787,6 +806,27 @@ globalOverlapNeighborIndices(const int locali, CHECK(index == numNonConstNeighbors); } +//------------------------------------------------------------------------------ +// Check whether NodeList ordering is appropriate for this function +//------------------------------------------------------------------------------ +template +bool +FlatConnectivity:: +fluidNodeListsFirst(const DataBase& dataBase) const +{ + auto nodeListi = 0; + auto nodeListItri = dataBase.nodeListBegin(); + auto nodeListItrj = dataBase.fluidNodeListBegin(); + for (; nodeListItrj != dataBase.fluidNodeListEnd(); + ++nodeListItri, ++nodeListItrj, ++nodeListi) { + if (*nodeListItri != *nodeListItrj) + { + return false; + } + } + return true; +} + // //------------------------------------------------------------------------------ // // Check whether two points overlap, given the H values // //------------------------------------------------------------------------------ diff --git a/src/KernelIntegrator/FlatConnectivity.hh b/src/KernelIntegrator/FlatConnectivity.hh index 90d912244..b503daa44 100644 --- a/src/KernelIntegrator/FlatConnectivity.hh +++ b/src/KernelIntegrator/FlatConnectivity.hh @@ -46,10 +46,17 @@ public: int firstGlobalIndex() const; int lastGlobalIndex() const; - // Local and global number of elements + // Number of nodes owned by this processor, including constant boundary nodes int numNodes() const; + + // Number of variable nodes owned by this processor int numInternalNodes() const; + + // Total number of variable nodes among all processors int numGlobalNodes() const; + + // Number of constant boundary nodes on this processor + // This doesn't include ghost nodes that are variable and owned by this or another processor int numBoundaryNodes() const; // Get local index from NodeList and Node indices @@ -64,10 +71,18 @@ public: // Number of neighbors, including self, for this point int numNeighbors(const int locali) const; int numOverlapNeighbors(const int locali) const; + + // Number of neighbors that are constant due to boundaries int numConstNeighbors(const int locali) const; int numConstOverlapNeighbors(const int locali) const; + + // Number of neighbors that are variables, on this proc or another int numNonConstNeighbors(const int locali) const; int numNonConstOverlapNeighbors(const int locali) const; + + // The sum of non-const neighbors over all internal nodes + // Corresponds to the size of the flattened connectivity in a matrix solve + int totalNumNonConstNeighbors() const; // For the point i, for its neighbor j, get the flattened index for point j int localToFlat(const int locali, const int localj) const; // return flatj @@ -157,7 +172,7 @@ public: // const Scalar extent) const; private: - + // Keep track of what has been initialized bool mIndexingInitialized; bool mGhostIndexingInitialized; @@ -218,6 +233,18 @@ private: // Scratch mutable ArrayDim mScratchArray; + + // The functions in this class assume that the FluidNodeLists come first in the + // connectivity order; that is, that the standard NodeLists are at the end of + // the alphabetical ordering in the DataBase. This is not generally true, but + // this class isn't used at the moment when this is not true. This function + // checks whether the standard NodeLists are at the end of the ordering. + // The root of the issue is that the connectivity for a point takes an integer index, + // while we iterate over NodeList iterators. We could either have the point connectivity + // accept a NodeList iterator or we could iterate over FluidNodeList indices. + // This function checks that the NodeLists are ordered as expected. + bool fluidNodeListsFirst(const DataBase& dataBase) const; + }; } // end namespace Spheral diff --git a/src/KernelIntegrator/FlatConnectivityInline.hh b/src/KernelIntegrator/FlatConnectivityInline.hh index 712527e21..b1fe618e9 100644 --- a/src/KernelIntegrator/FlatConnectivityInline.hh +++ b/src/KernelIntegrator/FlatConnectivityInline.hh @@ -243,6 +243,21 @@ numNonConstOverlapNeighbors(const int locali) const { return mNumOverlapNeighbors[locali] - mNumConstantBoundaryOverlapNeighbors[locali]; } +template +inline +int +FlatConnectivity:: +totalNumNonConstNeighbors() const { + CHECK(mIndexingInitialized); + CHECK(mBoundaryInformationInitialized); + auto total = 0; + for (auto locali = 0; locali < mNumInternalLocalNodes; ++locali) + { + total += mNumNeighbors[locali] - mNumConstantBoundaryNeighbors[locali]; + } + return total; +} + //------------------------------------------------------------------------------ // For the point i, for its neighbor j, get the flattened index for point j //------------------------------------------------------------------------------ diff --git a/src/KernelIntegrator/KernelIntegral.cc b/src/KernelIntegrator/KernelIntegral.cc index f5ecd8a7c..93d55f7e1 100644 --- a/src/KernelIntegrator/KernelIntegral.cc +++ b/src/KernelIntegrator/KernelIntegral.cc @@ -60,7 +60,7 @@ LinearKernelVector:: addToIntegral(const KernelIntegrationData& kid) { const auto coeff = this->mCoefficient->evaluateCoefficient(kid); const auto numIndices = kid.indices.size(); - CHECK(kid.dvalues.size() == numIndices); + CHECK(kid.values.size() == numIndices); CHECK(kid.indices.size() == numIndices); for (auto i = 0u; i < numIndices; ++i) { const size_t locali = kid.indices[i]; diff --git a/src/KernelIntegrator/KernelIntegrator.cc b/src/KernelIntegrator/KernelIntegrator.cc index 2334d7b30..a3a838b2f 100644 --- a/src/KernelIntegrator/KernelIntegrator.cc +++ b/src/KernelIntegrator/KernelIntegrator.cc @@ -82,15 +82,15 @@ performIntegration() { VERIFY2(omp_get_num_threads() == 1, "integration fails for > 1 OpenMP thread"); // Get some data out of database and state - const auto numNodeLists = mDataBase.numNodeLists(); + const auto numNodeLists = mDataBase.numFluidNodeLists(); const auto position = mState.fields(HydroFieldNames::position, Vector::zero); const auto H = mState.fields(HydroFieldNames::H, SymTensor::zero); const auto volume = mState.fields(HydroFieldNames::volume, 0.0); const auto cells = mState.fields(HydroFieldNames::cells, FacetedVolume()); // Since we may be receiving state from user, we need to make sure FieldLists aren't empty - VERIFY(position.size() == numNodeLists && - H.size() == numNodeLists && + VERIFY(position.size() >= numNodeLists && + H.size() >= numNodeLists && cells.size() == numNodeLists); // Initialize the integration quadrature for a single triangle