fix: host/device memory handling to run integrated tests on GPU

src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.cpp

@@ -1444,7 +1444,15 @@ void HypreMatrix::write( string const & filename,
       {
         // Copy distributed parcsr matrix in a local CSR matrix on every process with at least one row
         // Warning: works for a parcsr matrix that is smaller than 2^31-1
-        hypre_CSRMatrix * const fullMatrix = hypre_ParCSRMatrixToCSRMatrixAll( m_parcsr_mat );
+        // On GPU builds the parcsr data may live on device: gather from a host clone.
+        bool const onHost = hypre_ParCSRMatrixMemoryLocation( m_parcsr_mat ) == HYPRE_MEMORY_HOST;
+        hypre_ParCSRMatrix * const hostMat = onHost ? m_parcsr_mat
+                                                    : hypre_ParCSRMatrixClone_v2( m_parcsr_mat, 1, HYPRE_MEMORY_HOST );
+        hypre_CSRMatrix * const fullMatrix = hypre_ParCSRMatrixToCSRMatrixAll( hostMat );
+        if( !onHost )
+        {
+          GEOS_LAI_CHECK_ERROR( hypre_ParCSRMatrixDestroy( hostMat ) );
+        }
 
         // Identify the smallest process where CSRmatrix exists
         int const printRank = MpiWrapper::min( fullMatrix ? rank : MpiWrapper::commSize( comm() ), comm() );

src/coreComponents/mesh/WellElementSubRegion.cpp

@@ -1126,7 +1126,8 @@ void WellElementSubRegion::setElementStatus( arrayView1d< integer > const & loca
   {
     sortedPerfsPerElem[i]=activePerfsPerElement[indices[i]];
   }
-  // Set segment status
+  // Set segment status (updated on host; make sure the host copy is current and flagged as modified)
+  m_wellElementStatus.move( hostMemorySpace, true );
   std::fill( m_wellElementStatus.begin(), m_wellElementStatus.end(), 0 );
   if( numElements > 0 )
   {

src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp

@@ -1524,6 +1524,10 @@ real64 PhysicsSolverBase::scalingForSystemSolution( DomainPartition & GEOS_UNUSE
                              GEOS_FMT( "        {}: oscillation detected, scaling factor set to {}", getName(), scalingFactor ) );
     }
 
+    // The solution is produced on device and the history may have been captured by the
+    // detectOscillations kernel; stage both on host before the host-side append/erase.
+    localSolution.move( hostMemorySpace, false );
+    m_solutionHistory.move( hostMemorySpace, true );
     m_solutionHistory.appendArray( localSolution.begin(), localSolution.end());
     if( m_solutionHistory.size() > m_nonlinearSolverParameters.m_oscillationCheckDepth )
     {

src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp

@@ -825,6 +825,9 @@ void WellControls::setPerforationStatus( real64 const & time_n, WellElementSubRe
   PerforationData & perforationData = *subRegion.getPerforationData();
   string_array const & perfStatusTableName = perforationData.getPerfStatusTableName();
   arrayView1d< integer > perfStatus = perforationData.getLocalPerfStatus();
+  // The status is updated on host but read by device kernels (perforation fluxes).
+  // Without unified memory, the host copy must be current and flagged as modified.
+  perfStatus.move( hostMemorySpace, true );
   // for now set to open
   for( integer i=0; i<perforationData.size(); i++ )
   {
@@ -871,7 +874,8 @@ void WellControls::setPerforationStatus( real64 const & time_n, WellElementSubRe
   numOpenElements>0 ?  setWellStatus( time_n, WellControls::Status::OPEN ) :  setWellStatus( time_n, WellControls::Status::CLOSED );
 
 
-  // Set local well element status array
+  // Set local well element status array (read by device kernels: stage and touch on host)
+  localElemStatus.move( hostMemorySpace, true );
   for( integer i=0; i<subRegion.size(); i++ )
   {
     integer gi = globalWellElementIndex[i];

src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp

@@ -133,7 +133,9 @@ assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh,
 
     arrayView1d< integer const > const fractureState = subRegion.getField< contact::fractureState >();
 
-    forAll< serialPolicy >( subRegion.size(), [&]( localIndex const kfe )
+    // Capture by value: copying the views inside the host RAJA launch makes CHAI move them
+    // (and mark the matrix entries as touched) on host.
+    forAll< serialPolicy >( subRegion.size(), [=, this]( localIndex const kfe )
     {
       localIndex const kf0 = elemsToFaces[kfe][0], kf1 = elemsToFaces[kfe][1];
       localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 );

src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp

@@ -547,7 +547,9 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER
         constitutive::CoupledSolidBase & solid =
           this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, solidName );
 
-        arrayView1d< const real64 > const & averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
+        arrayView1d< real64 const > const averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
+        // The average is computed on device; without unified memory the view may hold a device pointer.
+        averageMeanTotalStressIncrement_k.move( hostMemorySpace, false );
         for( localIndex k = 0; k < localIndex( averageMeanTotalStressIncrement_k.size()); k++ )
         {
           averageMeanTotalStressIncrement.emplace_back( averageMeanTotalStressIncrement_k[k] );
@@ -571,11 +573,13 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER
         string const & solidName = subRegion.template getReference< string >( viewKeyStruct::porousMaterialNamesString() );
         constitutive::CoupledSolidBase & solid =
           this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, solidName );
-        auto & porosityModel = dynamic_cast< constitutive::BiotPorosity const & >( solid.getBasePorosityModel());
-        arrayView1d< real64 > const & averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
+        arrayView1d< real64 > const averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
+        // Write on host through a view moved (and touched) on host so that CHAI copies the
+        // updated values back to device before the next kernel (no unified memory).
+        averageMeanTotalStressIncrement_k.move( hostMemorySpace, true );
         for( localIndex k = 0; k < localIndex( averageMeanTotalStressIncrement_k.size()); k++ )
         {
-          porosityModel.updateAverageMeanTotalStressIncrement( k, averageMeanTotalStressIncrement[i] );
+          averageMeanTotalStressIncrement_k[k] = averageMeanTotalStressIncrement[i];
           i++;
         }
       } );

src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp

@@ -102,7 +102,9 @@ assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh,
 
     arrayView1d< integer const > const fractureState = subRegion.getField< contact::fractureState >();
 
-    forAll< serialPolicy >( subRegion.size(), [&]( localIndex const kfe )
+    // Capture by value: copying the views inside the host RAJA launch makes CHAI move them
+    // (and mark the matrix entries as touched) on host.
+    forAll< serialPolicy >( subRegion.size(), [=, this]( localIndex const kfe )
     {
       localIndex const kf0 = elemsToFaces[kfe][0], kf1 = elemsToFaces[kfe][1];
       localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 );

src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp

@@ -961,6 +961,19 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName,
   // So accumulation needs to multiply by area to cancel the 1/area in dAperturedU
   arrayView1d< real64 const > const area = subRegion.getElementArea().toViewConst();
 
+  // The loop below runs on host and captures by reference, so CHAI does not move the views.
+  // Without unified memory, stage everything it reads on host and mark the matrix entries as
+  // modified on host so that they are copied back before the next device access.
+  faceToNodeMap.move( hostMemorySpace, false );
+  dispDofNumber.move( hostMemorySpace, false );
+  bubbleDofNumber.move( hostMemorySpace, false );
+  presDofNumber.move( hostMemorySpace, false );
+  density.move( hostMemorySpace, false );
+  elemsToFaces.move( hostMemorySpace, false );
+  fractureState.move( hostMemorySpace, false );
+  area.move( hostMemorySpace, false );
+  localMatrix.move( hostMemorySpace, true );
+
   forAll< serialPolicy >( numElems, [&]( localIndex const kfe )
   {
     localIndex const kf0 = elemsToFaces[kfe][0];

src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp

@@ -1011,6 +1011,10 @@ void SolidMechanicsLagrangeContactBubbleStab::updateStickSlipList( DomainPartiti
       // elements of the same type are adjacent in the vals list.
       // This arrangement allows for efficient copying into the container
       // by leveraging parallelism.
+      // RAJA sorts receive raw pointers: CHAI does not migrate these views automatically.
+      // Make sure the sort operates on the device copies written by the kernel above.
+      keys_v.move( parallelDeviceMemorySpace, true );
+      vals_v.move( parallelDeviceMemorySpace, true );
       RAJA::sort_pairs< parallelDevicePolicy<> >( keys_v, vals_v );
 
       stickList.resize( nStick );
@@ -1030,8 +1034,10 @@ void SolidMechanicsLagrangeContactBubbleStab::updateStickSlipList( DomainPartiti
         slipList_v[kfe] = vals_v[nStick+kfe];
       } );
 
-      this->m_faceTypesToFaceElementsStick.get_inserted( meshName ).get_inserted( finiteElementName ) = stickList;
-      this->m_faceTypesToFaceElementsSlip.get_inserted( meshName ).get_inserted( finiteElementName ) = slipList;
+      // The lists were filled on device. Copy-assignment would deep-copy through host pointers
+      // (stale without unified memory); move-assignment transfers the CHAI buffers instead.
+      this->m_faceTypesToFaceElementsStick.get_inserted( meshName ).get_inserted( finiteElementName ) = std::move( stickList );
+      this->m_faceTypesToFaceElementsSlip.get_inserted( meshName ).get_inserted( finiteElementName ) = std::move( slipList );
 
       GEOS_LOG_LEVEL_RANK_0( logInfo::ConfigurationStatistics, GEOS_FMT( "# stick elements: {}, # slip elements: {}", nStick, nSlip ))
     } );
@@ -1105,6 +1111,9 @@ void SolidMechanicsLagrangeContactBubbleStab::createFaceTypeList( DomainPartitio
     // elements of the same type are adjacent in the vals list.
     // This arrangement allows for efficient copying into the container
     // by leveraging parallelism.
+    // Make sure the sort operates on the device copies written by the kernel above.
+    keys_v.move( parallelDeviceMemorySpace, true );
+    vals_v.move( parallelDeviceMemorySpace, true );
     RAJA::sort_pairs< parallelDevicePolicy<> >( keys_v, vals_v );
 
     quadList.resize( nQuad );
@@ -1122,8 +1131,8 @@ void SolidMechanicsLagrangeContactBubbleStab::createFaceTypeList( DomainPartitio
       quadList_v[kfe] = vals_v[nTri+kfe];
     } );
 
-    this->m_faceTypesToFaceElements.get_inserted( meshName ).get_inserted( "Quadrilateral" ) = quadList;
-    this->m_faceTypesToFaceElements.get_inserted( meshName ).get_inserted( "Triangle" ) = triList;
+    this->m_faceTypesToFaceElements.get_inserted( meshName ).get_inserted( "Quadrilateral" ) = std::move( quadList );
+    this->m_faceTypesToFaceElements.get_inserted( meshName ).get_inserted( "Triangle" ) = std::move( triList );
   } );
 
 }
@@ -1165,6 +1174,7 @@ void SolidMechanicsLagrangeContactBubbleStab::createBubbleCellList( DomainPartit
     }
 
     // Sort indexes to enable efficient searching using binary search.
+    tmpSpace_v.move( parallelDeviceMemorySpace, true );
     RAJA::stable_sort< parallelDevicePolicy<> >( tmpSpace_v );
     // Move data back to host: after the device sort the buffer lives on the GPU,
     // but SortedArray::insert is a host-only operation that dereferences the iterators.
@@ -1223,6 +1233,8 @@ void SolidMechanicsLagrangeContactBubbleStab::createBubbleCellList( DomainPartit
       // This arrangement allows for efficient copying into the container
       // by leveraging parallelism.
       localIndex nBubElems = static_cast< localIndex >(nBubElems_r.get());
+      keys_v.move( parallelDeviceMemorySpace, true );
+      perms_v.move( parallelDeviceMemorySpace, true );
       RAJA::sort_pairs< parallelDevicePolicy<> >( keys_v, perms_v );
 
       array1d< localIndex > bubbleElemsList;
@@ -1239,6 +1251,7 @@ void SolidMechanicsLagrangeContactBubbleStab::createBubbleCellList( DomainPartit
       {
         bubbleElemsList_v[k] = keys_v[k];
       } );
+      bubbleElemsList.move( hostMemorySpace, false );
       cellElementSubRegion.setBubbleElementsList( bubbleElemsList.toViewConst());
 
       forAll< parallelDevicePolicy<> >( n_max, [ = ] GEOS_HOST_DEVICE ( localIndex const k )
@@ -1259,6 +1272,7 @@ void SolidMechanicsLagrangeContactBubbleStab::createBubbleCellList( DomainPartit
         faceElemsList_v[k][0] = elemsToFaces[kfe][keys_v[k]];
         faceElemsList_v[k][1] = keys_v[k];
       } );
+      faceElemsList.move( hostMemorySpace, false );
       cellElementSubRegion.setFaceElementsList( faceElemsList.toViewConst());
 
     } );