Flat connectivity updates and axom 0.7.0

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})

scripts/spack/packages/spheral/package.py

@@ -52,8 +52,8 @@ class Spheral(CachedCMakePackage, CudaPackage):
     depends_on('[email protected] +shared +mpi +hdf5 -test ~parmetis', type='build', when='+mpi')
     depends_on('[email protected] +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('[email protected] ~shared ~adiak ~libdw ~papi ~libunwind +pic', type='build')
 

src/KernelIntegrator/BilinearIndex.cc

@@ -36,7 +36,7 @@ computeConnectivity(const DataBase<Dimension>& 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<Dimension>& dataBase,
   }
 
   // Get DataBase values
-  const auto numNodeLists = dataBase.numNodeLists();
+  const auto numNodeLists = dataBase.numFluidNodeLists();
   const auto& connectivityMap = dataBase.connectivityMap();
 
   // Get State values

src/KernelIntegrator/FlatConnectivity.cc

@@ -22,6 +22,23 @@
 
 namespace Spheral {
 
+namespace { // anonymous
+template<typename Dimension>
+int
+numFluidNeighbors(const std::vector<std::vector<int>>& connectivity,
+                  const DataBase<Dimension>& 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<typename Dimension>
 void 
 FlatConnectivity<Dimension>::
 computeIndices(const DataBase<Dimension>& 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<Dimension>& 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<Dimension>& 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<Dimension>& 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<Dimension>& 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<Dimension>& 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<CellFaceFlag>());
@@ -535,9 +554,9 @@ computeBoundaryInformation(const DataBase<Dimension>& 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<typename Dimension>
+bool
+FlatConnectivity<Dimension>::
+fluidNodeListsFirst(const DataBase<Dimension>& 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
 // //------------------------------------------------------------------------------

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<Dimension>& dataBase) const;
+
 };
 
 } // end namespace Spheral

src/KernelIntegrator/FlatConnectivityInline.hh

@@ -243,6 +243,21 @@ numNonConstOverlapNeighbors(const int locali) const {
   return mNumOverlapNeighbors[locali] - mNumConstantBoundaryOverlapNeighbors[locali];
 }
 
+template<typename Dimension>
+inline
+int
+FlatConnectivity<Dimension>::
+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
 //------------------------------------------------------------------------------

src/KernelIntegrator/KernelIntegral.cc

@@ -60,7 +60,7 @@ LinearKernelVector<Dimension>::
 addToIntegral(const KernelIntegrationData<Dimension>& 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];

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