From 52d57e2e491b286cf0c4c3cc104d5b90c168cc5c Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Mon, 26 Oct 2026 11:47:21 -0500 Subject: [PATCH] Draft. --- .../interfaces/hypre/HypreMatrix.cpp | 10 ++++++++- .../mesh/WellElementSubRegion.cpp | 3 ++- .../physicsSolvers/PhysicsSolverBase.cpp | 4 ++++ .../fluidFlow/wells/WellControls.cpp | 6 ++++- ...iphasePoromechanicsConformingFractures.cpp | 4 +++- .../multiphysics/PoromechanicsSolver.hpp | 12 ++++++---- ...ePhasePoromechanicsConformingFractures.cpp | 4 +++- ...asePoromechanicsConformingFracturesALM.cpp | 13 +++++++++++ ...olidMechanicsLagrangeContactBubbleStab.cpp | 22 +++++++++++++++---- 9 files changed, 65 insertions(+), 13 deletions(-) diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.cpp index a5f1b78c2d8..b93279abef6 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.cpp +++ b/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() ); diff --git a/src/coreComponents/mesh/WellElementSubRegion.cpp b/src/coreComponents/mesh/WellElementSubRegion.cpp index 0f94e6ebd84..59c2febd3cd 100644 --- a/src/coreComponents/mesh/WellElementSubRegion.cpp +++ b/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 ) { diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp index f5dc5cdb3b1..5aeaf3f5b41 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp +++ b/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 ) { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index c55bcd2c2b7..f2191063ef6 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/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; i0 ? 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 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 ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index 1b34ab94f65..d73c2bf09a3 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/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++; } } ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp index 947a0a6ef4b..ced151cc807 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp +++ b/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 ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp index 19962b85313..4833c0d63f3 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp +++ b/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]; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp index 34417df75db..f120d4646b1 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp +++ b/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()); } );