From bd68d98536a85e397ea53e0f197d8b768a62f743 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Mon, 2 Feb 2026 13:59:19 +0100 Subject: [PATCH 01/22] first pass --- .../linearAlgebra/CMakeLists.txt | 1 + .../interfaces/hypre/HypreMGR.cpp | 6 + ...asePoromechanicsConformingFracturesALM.hpp | 134 ++++++++++++++++++ .../utilities/LinearSolverParameters.hpp | 1 + ...asePoromechanicsConformingFracturesALM.hpp | 9 +- 5 files changed, 149 insertions(+), 2 deletions(-) create mode 100644 src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp diff --git a/src/coreComponents/linearAlgebra/CMakeLists.txt b/src/coreComponents/linearAlgebra/CMakeLists.txt index 01452771db9..957748d8f08 100644 --- a/src/coreComponents/linearAlgebra/CMakeLists.txt +++ b/src/coreComponents/linearAlgebra/CMakeLists.txt @@ -168,6 +168,7 @@ if( ENABLE_HYPRE ) interfaces/hypre/mgrStrategies/MultiphasePoromechanics.hpp interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFractures.hpp + interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsReservoirFVM.hpp interfaces/hypre/mgrStrategies/SinglePhaseReservoirFVM.hpp interfaces/hypre/mgrStrategies/SinglePhaseReservoirHybridFVM.hpp diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp index 6c24b749868..e70692272a9 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp @@ -37,6 +37,7 @@ #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanics.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFractures.hpp" +#include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsReservoirFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirHybridFVM.hpp" @@ -186,6 +187,11 @@ void hypre::mgr::createMGR( LinearSolverParameters const & params, setStrategy< SinglePhasePoromechanicsConformingFractures >( params.mgr, numComponentsPerField, precond, mgrData ); break; } + case LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsConformingFracturesALM: + { + setStrategy< SinglePhasePoromechanicsConformingFracturesALM >( params.mgr, numComponentsPerField, precond, mgrData ); + break; + } case LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsReservoirFVM: { setStrategy< SinglePhasePoromechanicsReservoirFVM >( params.mgr, numComponentsPerField, precond, mgrData ); diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp new file mode 100644 index 00000000000..ab86e0c8054 --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -0,0 +1,134 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file SinglePhasePoromechanicsConformingFractures.hpp + */ + +#ifndef GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_ +#define GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_ + +#include "linearAlgebra/interfaces/hypre/HypreMGR.hpp" + +namespace geos +{ + +namespace hypre +{ + +namespace mgr +{ + +/** + * @brief SinglePhasePoromechanicsConformingFractures strategy. + * + * dofLabel: 0 = displacement, x-component + * dofLabel: 1 = displacement, y-component + * dofLabel: 2 = displacement, z-component + * dofLabel: 3 = displacement bubble function, x-component + * dofLabel: 4 = displacement bubble function, y-component + * dofLabel: 5 = displacement bubble function, z-component + * dofLabel: 6 = pressure (cell elem + fracture elems) + + * + * Ingredients: + * 1. Level 1: F-points displacement (3,4,5), C-points pressure (0,1,2,6) + * 2. Level 2: F-points displacement (0,1,2), C-points pressure (6) + * 2. F-points smoother: BoomerAMG, single V-cycle + * 3. C-points coarse-grid/Schur complement solver: BoomerAMG + * 4. Global smoother: none + */ +class SinglePhasePoromechanicsConformingFracturesALM : public MGRStrategyBase< 2 > +{ +public: + + /** + * @brief Constructor. + */ + explicit SinglePhasePoromechanicsConformingFracturesALM( arrayView1d< int const > const & ) + : MGRStrategyBase( 7 ) + { + + // we keep u and p + m_labels[0].push_back( 0 ); + m_labels[0].push_back( 1 ); + m_labels[0].push_back( 2 ); + m_labels[0].push_back( 6 ); + // we keep p + m_labels[1].push_back( 6 ); + + setupLabels(); + + // Level 0 + m_levelFRelaxType[0] = MGRFRelaxationType::l1jacobi; + m_levelFRelaxIters[0] = 1; + m_levelGlobalSmootherType[0] = MGRGlobalSmootherType::none; + m_levelGlobalSmootherIters[0] = 0; + m_levelInterpType[0] = MGRInterpolationType::blockJacobi; + m_levelRestrictType[0] = MGRRestrictionType::injection; + m_levelCoarseGridMethod[0] = MGRCoarseGridMethod::galerkin; + + // Level 1 + m_levelFRelaxType[1] = MGRFRelaxationType::amgVCycle; + m_levelFRelaxIters[1] = 1; + m_levelGlobalSmootherType[1] = MGRGlobalSmootherType::none; + m_levelInterpType[1] = MGRInterpolationType::jacobi; + m_levelRestrictType[1] = MGRRestrictionType::injection; + m_levelCoarseGridMethod[1] = MGRCoarseGridMethod::nonGalerkin; + + } + + /** + * @brief Setup the MGR strategy. + * @param precond preconditioner wrapper + * @param mgrData auxiliary MGR data + */ + void setup( LinearSolverParameters::MGR const & mgrParams, + HyprePrecWrapper & precond, + HypreMGRData & mgrData ) + { + setReduction( precond, mgrData ); + + // Configure the BoomerAMG solver used as F-relaxation for the second level + // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGCreate( &mgrData.mechSolver.ptr ) ); + // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetTol( mgrData.mechSolver.ptr, 0.0 ) ); + // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxIter( mgrData.mechSolver.ptr, 1 ) ); + // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxRowSum( mgrData.mechSolver.ptr, 1.0 ) ); + // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetStrongThreshold( mgrData.mechSolver.ptr, 0.6 ) ); + // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPrintLevel( mgrData.mechSolver.ptr, 0 ) ); + // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumFunctions( mgrData.mechSolver.ptr, 3 ) ); + setDisplacementAMG(mgrData.mechSolver, mgrParams.separateComponents); + +#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCoarsenType( mgrData.mechSolver.ptr, hypre::getAMGCoarseningType( LinearSolverParameters::AMG::CoarseningType::PMIS ) ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxType( mgrData.mechSolver.ptr, hypre::getAMGRelaxationType( LinearSolverParameters::AMG::SmootherType::chebyshev ) ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumSweeps( mgrData.mechSolver.ptr, 1 ) ); +#else + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxOrder( mgrData.mechSolver.ptr, 1 ) ); +#endif + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetFSolverAtLevel( precond.ptr, mgrData.mechSolver.ptr, 1 ) ); + + // Configure the BoomerAMG solver used as mgr coarse solver for the pressure reduced system + setPressureAMG( mgrData.coarseSolver ); + } +}; + +} // namespace mgr + +} // namespace hypre + +} // namespace geos + +#endif /*GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_*/ \ No newline at end of file diff --git a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp index 9b491528ecb..c542bbce097 100644 --- a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp +++ b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp @@ -303,6 +303,7 @@ struct LinearSolverParameters hybridSinglePhasePoromechanics, ///< single phase poromechanics with hybrid finite volume single phase flow singlePhasePoromechanicsEmbeddedFractures, ///< single phase poromechanics with FV embedded fractures singlePhasePoromechanicsConformingFractures, ///< single phase poromechanics with conforming fractures + singlePhasePoromechanicsConformingFracturesALM, ///< single phase poromechanics with conforming fractures for ALM singlePhasePoromechanicsReservoirFVM, ///< single phase poromechanics with finite volume single phase flow with wells compositionalMultiphaseFVM, ///< finite volume compositional multiphase flow compositionalMultiphaseHybridFVM, ///< hybrid finite volume compositional multiphase flow diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp index 15f24d838d9..b7905a9f385 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -103,8 +103,13 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec virtual void setMGRStrategy() override final { - if( this->m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::mgr ) - GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}", this->getName(), this->getCatalogName())); + LinearSolverParameters & linearSolverParameters = this->m_linearSolverParameters.get(); + if( this->m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::mgr ){ + + linearSolverParameters.mgr.strategy = LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsConformingFracturesALM; + linearSolverParameters.mgr.separateComponents = true; + } + // GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}", this->getName(), this->getCatalogName())); } /**@}*/ From 8b2fc7f10411c61baf85f5d7ceca8087840586dc Mon Sep 17 00:00:00 2001 From: npillardou Date: Wed, 4 Feb 2026 11:12:15 +0100 Subject: [PATCH 02/22] Add FIM coupling for single phase fracture ALM --- ...asePoromechanicsConformingFracturesALM.cpp | 1110 ++++++++++++++++- ...asePoromechanicsConformingFracturesALM.hpp | 63 +- ...mechanicsConformingFracturesALMKernels.hpp | 1079 ++++++++++++++++ 3 files changed, 2219 insertions(+), 33 deletions(-) create mode 100644 src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp index ff7d9662a60..65da565f2b4 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp @@ -19,7 +19,19 @@ #include "SinglePhasePoromechanicsConformingFracturesALM.hpp" +#include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" +#include "constitutive/solid/CoupledSolidBase.hpp" +#include "constitutive/solid/PorousSolid.hpp" +#include "constitutive/contact/HydraulicApertureBase.hpp" +#include "constitutive/contact/HydraulicApertureRelationSelector.hpp" #include "finiteVolume/FluxApproximationBase.hpp" +#include "mesh/SurfaceElementRegion.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp" +#include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp" +#include "physicsSolvers/solidMechanics/contact/ContactFields.hpp" namespace geos { @@ -40,10 +52,13 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupCouplin { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain, dofManager ); - - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + // 1. Poromechanical coupling in the bulk (from base class) + Base::setupCoupling( domain, dofManager ); + // 2. Pressure - bubble displacement coupling in the fracture + dofManager.addCoupling( m_pressureKey, + contact::totalBubbleDisplacement::key(), + DofManager::Connector::Elem ); } @@ -55,13 +70,87 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupSystem( ParallelVector & solution, bool const setSparsity ) { - GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain, dofManager, localMatrix, rhs, solution, setSparsity ); + if( this->m_precond ) + { + this->m_precond->clear(); + } + + // Set domain on DofManager + dofManager.setDomain( domain ); + + // Initialize ALM contact solver internal data structures + // These must be called before DOF setup and assembly + this->solidMechanicsSolver()->createFaceTypeList( domain ); + this->solidMechanicsSolver()->updateStickSlipList( domain ); + this->solidMechanicsSolver()->createBubbleCellList( domain ); + + // Setup DOFs for all sub-solvers and coupling (uses PoromechanicsSolver::setupDofs) + // This adds: displacement DOFs, bubble DOFs, pressure DOFs, and all couplings + this->setupDofs( domain, dofManager ); + + // Reorder DOFs to optimize matrix structure + dofManager.reorderByRank(); + + if( setSparsity ) + { + // Get sparsity pattern from flow solver (could include wells) + SparsityPattern< globalIndex > patternOriginal; + this->flowSolver()->setSparsityPattern( domain, dofManager, localMatrix, patternOriginal ); + + // Get the original row lengths + array1d< localIndex > rowLengths( patternOriginal.numRows() ); + for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) + { + rowLengths[localRow] = patternOriginal.numNonZeros( localRow ); + } + + // Add the number of nonzeros induced by coupling + addTransmissibilityCouplingNNZ( domain, dofManager, rowLengths.toView() ); + addPressureForceCouplingNNZ( domain, dofManager, rowLengths.toView() ); + addMatrixPressureBubbleCouplingNNZ( domain, dofManager, rowLengths.toView() ); + + // Create a new pattern with enough capacity for coupled matrix + SparsityPattern< globalIndex > pattern; + pattern.resizeFromRowCapacities< parallelHostPolicy >( patternOriginal.numRows(), + patternOriginal.numColumns(), + rowLengths.data() ); + + // Copy the original nonzeros + for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) + { + globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow ) ); + } + + // Add the nonzeros from coupling + addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView() ); + addPressureForceCouplingPattern( domain, dofManager, pattern.toView() ); + addMatrixPressureBubbleCouplingPattern( domain, dofManager, pattern.toView() ); + + // Assemble the full sparsity pattern using the solid mechanics solver + this->solidMechanicsSolver()->setSparsityPattern( domain, dofManager, localMatrix, pattern ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + // Set up the derivative flux residual matrix + setUpDflux_dApertureMatrix( domain, dofManager, localMatrix ); + // Assimilate the sparsity pattern into the local matrix + localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) ); + } + + localMatrix.setName( this->getName() + "/matrix" ); + + rhs.setName( this->getName() + "/rhs" ); + rhs.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); + + solution.setName( this->getName() + "/solution" ); + solution.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); + + if( !this->m_precond && this->m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) + { + this->m_precond = this->createPreconditioner( domain ); + } } template< typename FLOW_SOLVER > @@ -72,12 +161,25 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleSyst CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { - GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( time_n, dt, domain, dofManager, localMatrix, localRhs ); + // Synchronize fracture state + this->solidMechanicsSolver()->synchronizeFractureState( domain ); + + // Assemble element-based contributions (mechanics + flow accumulation) + assembleElementBasedContributions( time_n, dt, domain, dofManager, localMatrix, localRhs ); + + // Assemble flux terms and get dFluidResidual/dAperture + this->flowSolver()->assembleHydrofracFluxTerms( time_n, + dt, + domain, + dofManager, + localMatrix, + localRhs, + getDerivativeFluxResidual_dNormalJump() ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + // Assemble coupling terms (must be after flux assembly) + assembleCouplingTerms( time_n, dt, domain, dofManager, localMatrix, localRhs ); } template< typename FLOW_SOLVER > @@ -90,10 +192,23 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleElem { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( time_n, dt, domain, dofManager, localMatrix, localRhs ); + // Assemble poromechanics terms (from base class) + Base::assembleElementBasedTerms( time_n, dt, domain, dofManager, localMatrix, localRhs ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + // Flow accumulation for fractures + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const, + FaceElementSubRegion const & subRegion ) + { + this->flowSolver()->accumulationAssemblyLaunch( dofManager, subRegion, localMatrix, localRhs ); + } ); + } ); + // Assemble contact terms (ALM) - note: assembleContact requires time and dt + this->solidMechanicsSolver()->assembleContact( time_n, dt, domain, dofManager, localMatrix, localRhs ); } template< typename FLOW_SOLVER > @@ -105,27 +220,164 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleCoup arrayView1d< real64 > const & localRhs ) { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain, dofManager, localMatrix, localRhs, time_n, dt ); + GEOS_UNUSED_VAR( time_n, dt ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + // These steps must occur after the fluxes are assembled because that's when DerivativeFluxResidual_dAperture is filled. + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshName, + MeshLevel const & mesh, + string_array const & regionNames ) + { + // Assemble Force Residual w.r.t. pressure (Aup) - fracture pressure contribution + assembleForceResidualDerivativeWrtPressure( meshName, mesh, regionNames, dofManager, localMatrix, localRhs ); + + // Assemble Fluid mass residual w.r.t. displacement (Apu) + assembleFluidMassResidualDerivativeWrtDisplacement( meshName, mesh, regionNames, dofManager, localMatrix, localRhs ); + } ); + + // Assemble matrix cell pressure contribution on bubble DOFs (Abp_matrix) + // This must be outside the lambda because it uses regionBasedKernelApplication + assembleMatrixPressureBubbleContribution( dt, const_cast< DomainPartition & >( domain ), dofManager, localMatrix, localRhs ); } template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::updateState( DomainPartition & domain ) { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + + // Call base poromechanics update + Base::updateState( domain ); + + // Need to call solid mechanics update separately to compute face displacement jump + this->solidMechanicsSolver()->updateState( domain ); + + // Remove the contribution of the hydraulic aperture from the stencil weights + this->flowSolver()->prepareStencilWeights( domain ); + + // Update hydraulic aperture and fracture permeability + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + ElementRegionManager & elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, + [&]( localIndex const, + FaceElementSubRegion & subRegion ) + { + arrayView2d< real64 const > const dispJump = subRegion.getField< contact::dispJump >(); + arrayView1d< real64 const > const area = subRegion.getElementArea(); + arrayView1d< real64 const > const volume = subRegion.getElementVolume(); + arrayView2d< real64 const > const fractureTraction = subRegion.getField< contact::traction >(); + arrayView1d< real64 const > const pressure = subRegion.getField< flow::pressure >(); + arrayView1d< real64 const > const oldHydraulicAperture = subRegion.getField< flow::aperture0 >(); + + arrayView1d< real64 > const aperture = subRegion.getElementAperture(); + arrayView1d< real64 > const hydraulicAperture = subRegion.getField< flow::hydraulicAperture >(); + arrayView1d< real64 > const deltaVolume = subRegion.getField< flow::deltaVolume >(); + arrayView1d< integer > const & fractureState = subRegion.getField< contact::fractureState >(); + + string const & hydraulicApertureRelationName = subRegion.template getReference< string >( Base::viewKeyStruct::hydraulicApertureRelationNameString() ); + HydraulicApertureBase const & hydraulicApertureModel = this->template getConstitutiveModel< HydraulicApertureBase >( subRegion, hydraulicApertureRelationName ); + + string const porousSolidName = subRegion.getReference< string >( FlowSolverBase::viewKeyStruct::solidNamesString() ); + CoupledSolidBase & porousSolid = subRegion.template getConstitutiveModel< CoupledSolidBase >( porousSolidName ); + + constitutiveUpdatePassThru( hydraulicApertureModel, [&] ( auto & castedHydraulicAperture ) + { + using HydraulicApertureType = TYPEOFREF( castedHydraulicAperture ); + typename HydraulicApertureType::KernelWrapper hydraulicApertureWrapper = castedHydraulicAperture.createKernelWrapper(); + + ConstitutivePassThru< CompressibleSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid ) + { + typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousMaterialWrapper = castedPorousSolid.createKernelUpdates(); + + poromechanicsFracturesKernels::StateUpdateKernel:: + launch< parallelDevicePolicy<> >( subRegion.size(), + porousMaterialWrapper, + hydraulicApertureWrapper, + dispJump, + pressure, + area, + volume, + deltaVolume, + aperture, + oldHydraulicAperture, + hydraulicAperture, + fractureTraction, + fractureState ); + } ); + } ); + } ); + } ); + + // Update the stencil weights using the updated hydraulic aperture + this->flowSolver()->updateStencilWeights( domain ); } template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: setUpDflux_dApertureMatrix( DomainPartition & domain, - DofManager const & GEOS_UNUSED_PARAM( dofManager ), + DofManager const & dofManager, CRSMatrix< real64, globalIndex > & localMatrix ) { - GEOS_UNUSED_VAR( domain, localMatrix ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + GEOS_UNUSED_VAR( dofManager ); + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & regionNames ) + { + std::unique_ptr< CRSMatrix< real64, localIndex > > & derivativeFluxResidual_dAperture = getRefDerivativeFluxResidual_dAperture(); + + { + // Calculate number of fracture elements + localIndex numRows = 0; + mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, + [&]( localIndex const, FaceElementSubRegion const & subRegion ) + { + numRows += subRegion.size(); + } ); + + // Number of columns (derivatives) = number of fracture elements + localIndex numCol = numRows; + + derivativeFluxResidual_dAperture = std::make_unique< CRSMatrix< real64, localIndex > >( numRows, numCol ); + derivativeFluxResidual_dAperture->setName( this->getName() + "/derivativeFluxResidual_dAperture" ); + + derivativeFluxResidual_dAperture->reserveNonZeros( localMatrix.numNonZeros() ); + localIndex maxRowSize = -1; + for( localIndex row = 0; row < localMatrix.numRows(); ++row ) + { + localIndex const rowSize = localMatrix.numNonZeros( row ); + maxRowSize = maxRowSize > rowSize ? maxRowSize : rowSize; + } + + for( localIndex row = 0; row < numRows; ++row ) + { + derivativeFluxResidual_dAperture->reserveNonZeros( row, maxRowSize ); + } + } + + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); + + fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) + { + localIndex const numFluxElems = stencil.stencilSize( iconn ); + typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); + + for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) + { + for( localIndex k1 = 0; k1 < numFluxElems; ++k1 ) + { + derivativeFluxResidual_dAperture->insertNonZero( sei[iconn][k0], sei[iconn][k1], 0.0 ); + } + } + } + } ); + } ); } template< typename FLOW_SOLVER > @@ -136,9 +388,62 @@ addTransmissibilityCouplingNNZ( DomainPartition const & domain, { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain, dofManager, rowLengths ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & ) + { + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const flowDofKey = dofManager.getKey( m_pressureKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); + + fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) + { + localIndex const numFluxElems = stencil.stencilSize( iconn ); + typename SurfaceElementStencil::IndexContainerViewConstType const & seri = stencil.getElementRegionIndices(); + typename SurfaceElementStencil::IndexContainerViewConstType const & sesri = stencil.getElementSubRegionIndices(); + typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); + + FaceElementSubRegion const & elementSubRegion = + elemManager.getRegion( seri[iconn][0] ).getSubRegion< FaceElementSubRegion >( sesri[iconn][0] ); + + ArrayOfArraysView< localIndex const > const elemsToNodes = elementSubRegion.nodeList().toViewConst(); + + arrayView1d< globalIndex const > const faceElementDofNumber = + elementSubRegion.getReference< array1d< globalIndex > >( flowDofKey ); + + for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) + { + globalIndex const activeFlowDOF = faceElementDofNumber[sei[iconn][k0]]; + globalIndex const rowNumber = activeFlowDOF - rankOffset; + if( rowNumber >= 0 && rowNumber < rowLengths.size() ) + { + for( localIndex k1 = 0; k1 < numFluxElems; ++k1 ) + { + // The coupling with the nodal displacements of the cell itself has already been added by the dofManager + // so we only add the coupling with the nodal displacements of the neighbors. + if( k1 != k0 ) + { + localIndex const numNodesPerElement = elemsToNodes[sei[iconn][k1]].size(); + // Nodal displacement DOFs (3 per node, 2 faces) + rowLengths[rowNumber] += 3 * numNodesPerElement; + // Bubble DOFs (3 per face, 2 faces = 6 total) + rowLengths[rowNumber] += 6; + } + } + } + } + } + } ); + } ); } template< typename FLOW_SOLVER > @@ -149,38 +454,783 @@ addTransmissibilityCouplingPattern( DomainPartition const & domain, { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain, dofManager, pattern ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + using namespace contact; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + string const flowDofKey = dofManager.getKey( m_pressureKey ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const & + bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + + // Get the finite volume method used to compute the fluxes + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fvDiscretization = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); + + SurfaceElementRegion const & fractureRegion = + elemManager.getRegion< SurfaceElementRegion >( this->solidMechanicsSolver()->getUniqueFractureRegionName() ); + FaceElementSubRegion const & fractureSubRegion = + fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); + + GEOS_ERROR_IF( !fractureSubRegion.hasWrapper( flow::pressure::key() ), + this->getDataContext() << ": The fracture subregion must contain pressure field." ); + + arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst(); + + arrayView1d< globalIndex const > const & + flowDofNumber = fractureSubRegion.getReference< globalIndex_array >( flowDofKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + fvDiscretization.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + forAll< serialPolicy >( stencil.size(), [=] ( localIndex const iconn ) + { + localIndex const numFluxElems = stencil.stencilSize( iconn ); + + // A fracture connector has to be an edge shared by two faces + if( numFluxElems == 2 ) + { + typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); + + // First index: face element. Second index: node + for( localIndex kf = 0; kf < 2; ++kf ) + { + // Set row DOF index + // Note that the 1-kf index is intentional, as this is coupling the pressure of one face cell + // to the nodes of the adjacent cell + localIndex const rowIndex = flowDofNumber[sei[iconn][1 - kf]] - rankOffset; + + if( rowIndex >= 0 && rowIndex < pattern.numRows() ) + { + // Get fracture, face and region/subregion/element indices (for elements on both sides) + localIndex const fractureIndex = sei[iconn][kf]; + + // Get the number of nodes + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elem2dToFaces[fractureIndex][0] ); + + // Loop over the two sides of each fracture element + for( localIndex kf1 = 0; kf1 < 2; ++kf1 ) + { + localIndex const faceIndex = elem2dToFaces[fractureIndex][kf1]; + + // Save the list of DOF associated with nodes (displacement) + for( localIndex a = 0; a < numNodesPerFace; ++a ) + { + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const colIndex = dispDofNumber[faceToNodeMap( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i ); + pattern.insertNonZero( rowIndex, colIndex ); + } + } + + // Save the list of DOF associated with bubble displacement + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const colIndex = bubbleDofNumber[faceIndex] + LvArray::integerConversion< globalIndex >( i ); + pattern.insertNonZero( rowIndex, colIndex ); + } + } + } + } + } + } ); + } ); + } ); +} + +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: +addPressureForceCouplingNNZ( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< localIndex > const & rowLengths ) const +{ + GEOS_MARK_FUNCTION; + + using namespace contact; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const & + bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + + globalIndex const rankOffset = dofManager.rankOffset(); + + string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); + SurfaceElementRegion const & fractureRegion = + elemManager.getRegion< SurfaceElementRegion >( fractureRegionName ); + FaceElementSubRegion const & fractureSubRegion = + fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); + + arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst(); + + // For each fracture element, add NNZ for (displacement_row, pressure_col) and (bubble_row, pressure_col) + forAll< serialPolicy >( fractureSubRegion.size(), [=, &rowLengths] ( localIndex const kfe ) + { + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elem2dToFaces[kfe][0] ); + + // For displacement DOFs: add 1 pressure column per displacement DOF row + for( localIndex kf = 0; kf < 2; ++kf ) + { + localIndex const faceIndex = elem2dToFaces[kfe][kf]; + for( localIndex a = 0; a < numNodesPerFace; ++a ) + { + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowNumber = dispDofNumber[faceToNodeMap( faceIndex, a )] + i - rankOffset; + if( rowNumber >= 0 && rowNumber < rowLengths.size() ) + { + rowLengths[rowNumber] += 1; // One pressure column + } + } + } + } + + // For bubble DOFs: add 1 pressure column per bubble DOF row + for( localIndex kf = 0; kf < 2; ++kf ) + { + localIndex const faceIndex = elem2dToFaces[kfe][kf]; + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowNumber = bubbleDofNumber[faceIndex] + i - rankOffset; + if( rowNumber >= 0 && rowNumber < rowLengths.size() ) + { + rowLengths[rowNumber] += 1; // One pressure column + } + } + } + } ); + } ); +} + +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: +addPressureForceCouplingPattern( DomainPartition const & domain, + DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern ) const +{ + GEOS_MARK_FUNCTION; + + using namespace contact; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + string const flowDofKey = dofManager.getKey( m_pressureKey ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const & + bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + + globalIndex const rankOffset = dofManager.rankOffset(); + + string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); + SurfaceElementRegion const & fractureRegion = + elemManager.getRegion< SurfaceElementRegion >( fractureRegionName ); + FaceElementSubRegion const & fractureSubRegion = + fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); + + arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst(); + arrayView1d< globalIndex const > const & + flowDofNumber = fractureSubRegion.getReference< globalIndex_array >( flowDofKey ); + + // For each fracture element, add pattern for (displacement_row, pressure_col) and (bubble_row, pressure_col) + forAll< serialPolicy >( fractureSubRegion.size(), [=] ( localIndex const kfe ) + { + globalIndex const pressureColIndex = flowDofNumber[kfe]; + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elem2dToFaces[kfe][0] ); + + // For displacement DOFs + for( localIndex kf = 0; kf < 2; ++kf ) + { + localIndex const faceIndex = elem2dToFaces[kfe][kf]; + for( localIndex a = 0; a < numNodesPerFace; ++a ) + { + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowIndex = dispDofNumber[faceToNodeMap( faceIndex, a )] + i - rankOffset; + if( rowIndex >= 0 && rowIndex < pattern.numRows() ) + { + pattern.insertNonZero( rowIndex, pressureColIndex ); + } + } + } + } + + // For bubble DOFs + for( localIndex kf = 0; kf < 2; ++kf ) + { + localIndex const faceIndex = elem2dToFaces[kfe][kf]; + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowIndex = bubbleDofNumber[faceIndex] + i - rankOffset; + if( rowIndex >= 0 && rowIndex < pattern.numRows() ) + { + pattern.insertNonZero( rowIndex, pressureColIndex ); + } + } + } + } ); + } ); } template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: assembleForceResidualDerivativeWrtPressure( string const & meshName, MeshLevel const & mesh, - arrayView1d< string const > const & regionNames, + string_array const & regionNames, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( meshName, mesh, regionNames, dofManager, localMatrix, localRhs ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + using namespace contact; + + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + + string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const & bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + string const & flowDofKey = dofManager.getKey( m_pressureKey ); + + arrayView1d< globalIndex const > const dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + + string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); + + // Use the same kernel launch pattern as SolidMechanicsAugmentedLagrangianContact::assembleForceResidualPressureContribution + this->solidMechanicsSolver()->template forFiniteElementOnFractureSubRegions( meshName, + [&] ( string const &, + finiteElement::FiniteElementBase const & subRegionFE, + arrayView1d< localIndex const > const & faceElementList ) + { + // Get pressure DOF number from the fracture subregion + SurfaceElementRegion const & fractureRegion = mesh.getElemManager().getRegion< SurfaceElementRegion >( fractureRegionName ); + FaceElementSubRegion const & fractureSubRegion = fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); + arrayView1d< globalIndex const > const pressureDofNumber = fractureSubRegion.getReference< array1d< globalIndex > >( flowDofKey ); + + poromechanicsALMKernels::AssembleForceResidualDerivativeWrtPressureFactory + kernelFactory( dispDofNumber, + bubbleDofNumber, + dofManager.rankOffset(), + localMatrix, + localRhs, + 0.0, // dt not used + faceElementList, + pressureDofNumber ); + + // Note: const_cast is needed because interfaceBasedKernelApplication takes non-const mesh + // even though it only modifies the matrix/rhs which are passed separately + real64 maxResidual = finiteElement:: + interfaceBasedKernelApplication + < parallelDevicePolicy<>, + constitutive::NullModel >( const_cast< MeshLevel & >( mesh ), + fractureRegionName, + faceElementList, + subRegionFE, + "", + kernelFactory ); + + GEOS_UNUSED_VAR( maxResidual ); + } ); + + GEOS_UNUSED_VAR( regionNames ); } template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: -assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, - arrayView1d< string const > const & regionNames, +assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, + MeshLevel const & mesh, + string_array const & regionNames, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & GEOS_UNUSED_PARAM( localRhs ) ) + arrayView1d< real64 > const & localRhs ) +{ + GEOS_MARK_FUNCTION; + GEOS_UNUSED_VAR( regionNames ); + + using namespace contact; + + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + + CRSMatrixView< real64 const, localIndex const > const & + dFluxResidual_dNormalJump = getDerivativeFluxResidual_dNormalJump().toViewConst(); + + string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const & bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + string const & presDofKey = dofManager.getKey( m_pressureKey ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const & + bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + globalIndex const rankOffset = dofManager.rankOffset(); + + string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); + + // Maximum DOF sizes + constexpr localIndex maxNumNodesPerFace = m_maxFaceNodes; + constexpr localIndex maxNumUdofs = 2 * 3 * maxNumNodesPerFace; // 66 + constexpr localIndex numBdofs = 6; + + // Get the fracture subRegion + SurfaceElementRegion const & fractureRegion = elemManager.getRegion< SurfaceElementRegion >( fractureRegionName ); + FaceElementSubRegion const & subRegion = fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); + + localIndex const numElems = subRegion.size(); + + // Allocate temporary storage for aperture derivatives (computed via kernels) + array2d< real64 > dAperturedU( numElems, maxNumUdofs ); + array2d< real64 > dAperturedB( numElems, numBdofs ); + + // Initialize to zero and move to device for kernel access + dAperturedU.zero(); + dAperturedB.zero(); + dAperturedU.move( parallelDeviceMemorySpace, true ); + dAperturedB.move( parallelDeviceMemorySpace, true ); + + // Launch the ComputeApertureDerivatives kernel to fill dAperturedU and dAperturedB + // This is called for each element type (tri, quad, etc.) + this->solidMechanicsSolver()->template forFiniteElementOnFractureSubRegions( meshName, + [&] ( string const &, + finiteElement::FiniteElementBase const & subRegionFE, + arrayView1d< localIndex const > const & faceElementList ) + { + poromechanicsALMKernels::ComputeApertureDerivativesFactory + kernelFactory( dispDofNumber, + bubbleDofNumber, + dofManager.rankOffset(), + localMatrix, + localRhs, + 0.0, // dt not used + faceElementList, + dAperturedU.toView(), + dAperturedB.toView() ); + + real64 maxResidual = finiteElement:: + interfaceBasedKernelApplication + < parallelDevicePolicy<>, + constitutive::NullModel >( const_cast< MeshLevel & >( mesh ), + fractureRegionName, + faceElementList, + subRegionFE, + "", + kernelFactory ); + + GEOS_UNUSED_VAR( maxResidual ); + } ); + + // Move data to host for serial assembly + dAperturedU.move( hostMemorySpace ); + dAperturedB.move( hostMemorySpace ); + + // Now assemble using the pre-computed derivatives + string const & fluidName = subRegion.getReference< string >( FlowSolverBase::viewKeyStruct::fluidNamesString() ); + + SingleFluidBase const & fluid = this->template getConstitutiveModel< SingleFluidBase >( subRegion, fluidName ); + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const & density = fluid.density(); + + arrayView1d< globalIndex const > const & presDofNumber = subRegion.getReference< array1d< globalIndex > >( presDofKey ); + + arrayView2d< localIndex const > const & elemsToFaces = subRegion.faceList().toViewConst(); + + arrayView1d< integer const > const fractureState = subRegion.getField< contact::fractureState >(); + + // Get element area for proper scaling + // Note: dAperturedU/dB are computed as (1/area) * unitNormal^T * Atu/Atb + // For accumulation: dR_accum/du = density * unitNormal^T * Atu (no 1/area factor) + // For flux: dR_flux/du = dR/dAperture * (1/area) * unitNormal^T * Atu (with 1/area factor) + // So accumulation needs to multiply by area to cancel the 1/area in dAperturedU + arrayView1d< real64 const > const area = subRegion.getElementArea().toViewConst(); + + forAll< serialPolicy >( numElems, [&]( localIndex const kfe ) + { + localIndex const kf0 = elemsToFaces[kfe][0]; + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 ); + localIndex const numUdofs = 2 * 3 * numNodesPerFace; + + globalIndex nodeDOF[maxNumUdofs]; + globalIndex elemDOF[1]; + elemDOF[0] = presDofNumber[kfe]; + + stackArray1d< real64, maxNumUdofs > dRdU( maxNumUdofs ); + + bool const isFractureOpen = ( fractureState[kfe] == FractureState::Open ); + + // ==== Part 1: Apu - Nodal displacement contribution ==== + // Get DOF indices for displacement + for( localIndex kf = 0; kf < 2; ++kf ) + { + for( localIndex a = 0; a < numNodesPerFace; ++a ) + { + for( localIndex i = 0; i < 3; ++i ) + { + nodeDOF[kf * 3 * numNodesPerFace + 3 * a + i] = dispDofNumber[faceToNodeMap( elemsToFaces[kfe][kf], a )] + + LvArray::integerConversion< globalIndex >( i ); + } + } + } + + // Accumulation derivative w.r.t. nodal displacement + // dR_accum/du = density * unitNormal^T * Atu = density * area * dAperturedU + // (dAperturedU already has 1/area factor, so multiply by area to cancel it) + // Note: The ALM Atu matrix has opposite sign convention compared to LM's dAper/dU, + // so we need to negate the result to match the LM sign convention. + if( isFractureOpen ) + { + for( localIndex j = 0; j < numUdofs; ++j ) + { + dRdU( j ) = -density[kfe][0] * dAperturedU( kfe, j ) * area[kfe]; + } + + localIndex const localRow = LvArray::integerConversion< localIndex >( elemDOF[0] - rankOffset ); + + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( localRow, + nodeDOF, + dRdU.data(), + numUdofs ); + } + } + + // Flux derivative w.r.t. nodal displacement + localIndex const numColumns = dFluxResidual_dNormalJump.numNonZeros( kfe ); + arraySlice1d< localIndex const > const & columns = dFluxResidual_dNormalJump.getColumns( kfe ); + arraySlice1d< real64 const > const & values = dFluxResidual_dNormalJump.getEntries( kfe ); + + for( localIndex kfe1 = 0; kfe1 < numColumns; ++kfe1 ) + { + real64 const dR_dAper = values[kfe1]; + localIndex const kfe2 = columns[kfe1]; + + bool const isOpen = ( fractureState[kfe2] == FractureState::Open ); + if( !isOpen && !isFractureOpen ) + continue; + + localIndex const kf0_2 = elemsToFaces[kfe2][0]; + localIndex const numNodesPerFace2 = faceToNodeMap.sizeOfArray( kf0_2 ); + localIndex const numUdofs2 = 2 * 3 * numNodesPerFace2; + + // Get DOF indices for element kfe2 + globalIndex nodeDOF2[maxNumUdofs]; + for( localIndex kf = 0; kf < 2; ++kf ) + { + for( localIndex a = 0; a < numNodesPerFace2; ++a ) + { + for( localIndex i = 0; i < 3; ++i ) + { + nodeDOF2[kf * 3 * numNodesPerFace2 + 3 * a + i] = dispDofNumber[faceToNodeMap( elemsToFaces[kfe2][kf], a )] + + LvArray::integerConversion< globalIndex >( i ); + } + } + } + + // dR_flux/du = dR_flux/dAper * dAper/du (pre-computed for element kfe2) + // Note: Negate to match LM sign convention (ALM Atu has opposite sign) + stackArray1d< real64, maxNumUdofs > dRdU2( maxNumUdofs ); + for( localIndex j = 0; j < numUdofs2; ++j ) + { + dRdU2( j ) = -dR_dAper * dAperturedU( kfe2, j ); + } + + localIndex const localRow = LvArray::integerConversion< localIndex >( elemDOF[0] - rankOffset ); + + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( localRow, + nodeDOF2, + dRdU2.data(), + numUdofs2 ); + } + } + + // ==== Part 2: Apb - Bubble displacement contribution ==== + globalIndex bubbleDOF[numBdofs]; + + // Get DOF indices for bubble + for( localIndex kf = 0; kf < 2; ++kf ) + { + localIndex const faceIndex = elemsToFaces[kfe][kf]; + for( localIndex i = 0; i < 3; ++i ) + { + bubbleDOF[kf * 3 + i] = bubbleDofNumber[faceIndex] + LvArray::integerConversion< globalIndex >( i ); + } + } + + // Accumulation derivative w.r.t. bubble displacement + // dR_accum/db = density * unitNormal^T * Atb = density * area * dAperturedB + // (dAperturedB already has 1/area factor, so multiply by area to cancel it) + // Note: Negate to match LM sign convention (ALM Atb has opposite sign) + if( isFractureOpen ) + { + real64 dRdB[numBdofs]; + for( localIndex j = 0; j < numBdofs; ++j ) + { + dRdB[j] = -density[kfe][0] * dAperturedB( kfe, j ) * area[kfe]; + } + + localIndex const localRow = LvArray::integerConversion< localIndex >( elemDOF[0] - rankOffset ); + + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( localRow, + bubbleDOF, + dRdB, + numBdofs ); + } + } + + // Flux derivative w.r.t. bubble displacement + for( localIndex kfe1 = 0; kfe1 < numColumns; ++kfe1 ) + { + real64 const dR_dAper = values[kfe1]; + localIndex const kfe2 = columns[kfe1]; + + bool const isOpen = ( fractureState[kfe2] == FractureState::Open ); + if( !isOpen && !isFractureOpen ) + continue; + + // Get DOF indices for bubble of element kfe2 + globalIndex bubbleDOF2[numBdofs]; + for( localIndex kf = 0; kf < 2; ++kf ) + { + localIndex const faceIndex = elemsToFaces[kfe2][kf]; + for( localIndex i = 0; i < 3; ++i ) + { + bubbleDOF2[kf * 3 + i] = bubbleDofNumber[faceIndex] + LvArray::integerConversion< globalIndex >( i ); + } + } + + // dR_flux/db = dR_flux/dAper * dAper/db (pre-computed for element kfe2) + // Note: Negate to match LM sign convention (ALM Atb has opposite sign) + real64 dRdB2[numBdofs]; + for( localIndex j = 0; j < numBdofs; ++j ) + { + dRdB2[j] = -dR_dAper * dAperturedB( kfe2, j ); + } + + localIndex const localRow = LvArray::integerConversion< localIndex >( elemDOF[0] - rankOffset ); + + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( localRow, + bubbleDOF2, + dRdB2, + numBdofs ); + } + } + } ); +} + +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: +addMatrixPressureBubbleCouplingNNZ( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< localIndex > const & rowLengths ) const +{ + GEOS_MARK_FUNCTION; + + using namespace contact; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & regionNames ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + // Loop over matrix cell regions that have bubbles + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, + [&]( localIndex const, CellElementSubRegion const & subRegion ) + { + if( !subRegion.hasWrapper( "bubbleElementsList" ) ) + return; + + arrayView1d< localIndex const > const bubbleElems = subRegion.bubbleElementsList(); + arrayView2d< localIndex const > const elemsToFaces = subRegion.faceElementsList(); + + forAll< serialPolicy >( bubbleElems.size(), [=, &rowLengths]( localIndex const kk ) + { + localIndex const faceIndex = elemsToFaces[kk][0]; + + // Add 1 pressure column for each of the 3 bubble DOFs + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowNumber = bubbleDofNumber[faceIndex] + i - rankOffset; + if( rowNumber >= 0 && rowNumber < rowLengths.size() ) + { + rowLengths[rowNumber] += 1; // One pressure DOF from matrix cell + } + } + } ); + } ); + } ); +} + +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: +addMatrixPressureBubbleCouplingPattern( DomainPartition const & domain, + DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern ) const { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( mesh, regionNames, dofManager, localMatrix ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + using namespace contact; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & regionNames ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + string const flowDofKey = dofManager.getKey( m_pressureKey ); + + arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + // Loop over matrix cell regions that have bubbles + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, + [&]( localIndex const, CellElementSubRegion const & subRegion ) + { + if( !subRegion.hasWrapper( "bubbleElementsList" ) ) + return; + + arrayView1d< localIndex const > const bubbleElems = subRegion.bubbleElementsList(); + arrayView2d< localIndex const > const elemsToFaces = subRegion.faceElementsList(); + arrayView1d< globalIndex const > const pressureDofNumber = subRegion.getReference< array1d< globalIndex > >( flowDofKey ); + + forAll< serialPolicy >( bubbleElems.size(), [=]( localIndex const kk ) + { + localIndex const k = bubbleElems[kk]; + localIndex const faceIndex = elemsToFaces[kk][0]; + globalIndex const pressureColIndex = pressureDofNumber[k]; + + // Add pattern for (bubble_row, pressure_col) + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowIndex = bubbleDofNumber[faceIndex] + i - rankOffset; + if( rowIndex >= 0 && rowIndex < pattern.numRows() ) + { + pattern.insertNonZero( rowIndex, pressureColIndex ); + } + } + } ); + } ); + } ); +} + +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: +assembleMatrixPressureBubbleContribution( real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + GEOS_MARK_FUNCTION; + + using namespace contact; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + NodeManager const & nodeManager = mesh.getNodeManager(); + FaceManager const & faceManager = mesh.getFaceManager(); + ElementRegionManager & elemManager = mesh.getElemManager(); + + string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const & bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + string const & flowDofKey = dofManager.getKey( m_pressureKey ); + + arrayView1d< globalIndex const > const dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + + // Identify poromechanics regions (cells with porous solid) + set< string > poromechanicsRegions; + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, + [&]( localIndex const regionIndex, CellElementSubRegion const & subRegion ) + { + if( subRegion.hasWrapper( FlowSolverBase::viewKeyStruct::solidNamesString() ) ) + { + poromechanicsRegions.insert( regionNames[regionIndex] ); + } + } ); + + string_array poromechanicsRegionNames; + poromechanicsRegionNames.reserve( poromechanicsRegions.size() ); + for( auto const & region : poromechanicsRegions ) + { + poromechanicsRegionNames.emplace_back( region ); + } + + // Launch the kernel on matrix cells with bubbles + poromechanicsMatrixBubbleKernels::MatrixPressureBubbleFactory kernelFactory( dispDofNumber, + bubbleDofNumber, + dofManager.rankOffset(), + localMatrix, + localRhs, + dt, + flowDofKey ); + + real64 maxResidual = finiteElement::regionBasedKernelApplication + < parallelDevicePolicy<>, + constitutive::PorousSolid< constitutive::ElasticIsotropic, constitutive::ConstantPermeability >, + CellElementSubRegion >( mesh, + poromechanicsRegionNames, + this->solidMechanicsSolver()->getDiscretizationName(), + FlowSolverBase::viewKeyStruct::solidNamesString(), + kernelFactory ); + GEOS_UNUSED_VAR( maxResidual ); + } ); } diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp index 15f24d838d9..9a427168a21 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -141,13 +141,14 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec void assembleForceResidualDerivativeWrtPressure( string const & meshName, MeshLevel const & mesh, - arrayView1d< string const > const & regionNames, + string_array const & regionNames, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ); - void assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, - arrayView1d< string const > const & regionNames, + void assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, + MeshLevel const & mesh, + string_array const & regionNames, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ); @@ -172,6 +173,62 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec DofManager const & dofManager, SparsityPatternView< globalIndex > const & pattern ) const; + /** + * @Brief add the nnz induced by the pressure-force coupling (Aup, Abp) + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param rowLenghts the nnz in each row + */ + void addPressureForceCouplingNNZ( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< localIndex > const & rowLengths ) const; + + /** + * @Brief add the sparsity pattern induced by the pressure-force coupling (Aup, Abp) + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param pattern the sparsity pattern + */ + void addPressureForceCouplingPattern( DomainPartition const & domain, + DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern ) const; + + /** + * @Brief add the nnz induced by the matrix pressure-bubble coupling (Abp_matrix) + * This handles the contribution of matrix cell pressure on bubble DOFs. + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param rowLengths the nnz in each row + */ + void addMatrixPressureBubbleCouplingNNZ( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< localIndex > const & rowLengths ) const; + + /** + * @Brief add the sparsity pattern induced by the matrix pressure-bubble coupling + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param pattern the sparsity pattern + */ + void addMatrixPressureBubbleCouplingPattern( DomainPartition const & domain, + DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern ) const; + + /** + * @Brief assemble the contribution of matrix cell pressure on bubble DOFs + * with full Jacobian for fully-implicit coupling. + * @param dt the time step size + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param localMatrix the local system matrix + * @param localRhs the local system right-hand side vector + */ + void assembleMatrixPressureBubbleContribution( real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + /** * @brief Set up the Dflux_dApertureMatrix object * diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp new file mode 100644 index 00000000000..f069b7f7b90 --- /dev/null +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp @@ -0,0 +1,1079 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file SinglePhasePoromechanicsConformingFracturesALMKernels.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALMKERNELS_HPP_ +#define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALMKERNELS_HPP_ + +#include "finiteElement/kernelInterface/ImplicitKernelBase.hpp" +#include "physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsConformingContactKernelsBase.hpp" +#include "physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsConformingContactKernelsHelper.hpp" +#include "physicsSolvers/solidMechanics/contact/ContactFields.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" + +namespace geos +{ + +namespace poromechanicsALMKernels +{ + +/** + * @brief Kernel to assemble the pressure force contribution to displacement equations + * AND the Jacobian derivative dR_u/dP for fully implicit coupling. + * + * This kernel is similar to SolidMechanicsConformingPressureContributionKernels::AssemblePressureContribution + * but also adds the Jacobian terms needed for fully implicit poromechanics coupling. + * The sign convention matches the existing sequential coupling kernel. + */ +template< typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +class AssembleForceResidualDerivativeWrtPressure : + public solidMechanicsConformingContactKernels::ConformingContactKernelsBase< CONSTITUTIVE_TYPE, + FE_TYPE > +{ + +public: + /// Alias for the base class + using Base = solidMechanicsConformingContactKernels::ConformingContactKernelsBase< CONSTITUTIVE_TYPE, + FE_TYPE >; + + /// Number of nodes per element + static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem; + + /// The number of displacement dofs per element. + static constexpr int numUdofs = Base::numUdofs; + + /// The number of bubble dofs per element. + static constexpr int numBdofs = Base::numBdofs; + + using Base::m_elemsToFaces; + using Base::m_faceToNodes; + using Base::m_X; + using Base::m_finiteElementSpace; + using Base::m_dofNumber; + using Base::m_bDofNumber; + using Base::m_dofRankOffset; + using Base::m_rotationMatrix; + using Base::m_matrix; + using Base::m_rhs; + + /** + * @brief Constructor + */ + AssembleForceResidualDerivativeWrtPressure( NodeManager const & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + FaceElementSubRegion & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + arrayView1d< globalIndex const > const uDofNumber, + arrayView1d< globalIndex const > const bDofNumber, + globalIndex const rankOffset, + CRSMatrixView< real64, globalIndex const > const inputMatrix, + arrayView1d< real64 > const inputRhs, + real64 const inputDt, + arrayView1d< localIndex const > const & faceElementList, + arrayView1d< globalIndex const > const & pressureDofNumber ): + Base( nodeManager, + edgeManager, + faceManager, + targetRegionIndex, + elementSubRegion, + finiteElementSpace, + inputConstitutiveType, + uDofNumber, + bDofNumber, + rankOffset, + inputMatrix, + inputRhs, + inputDt, + faceElementList ), + m_pressure( elementSubRegion.getField< fields::flow::pressure >().toView() ), + m_pressureDofNumber( pressureDofNumber ) + {} + + /** + * @copydoc finiteElement::InterfaceKernelBase::StackVariables + */ + struct StackVariables : public Base::StackVariables + { +public: + + GEOS_HOST_DEVICE + StackVariables(): + Base::StackVariables(), + dispEqnRowIndices{}, + bEqnRowIndices{}, + pressureColIndex( 0 ), + unitNormal{}, + localRu{}, + localRb{}, + dRudP{}, + dRbdP{} + {} + + /// Row indices for displacement equations + globalIndex dispEqnRowIndices[numUdofs]; + + /// Row indices for bubble equations + globalIndex bEqnRowIndices[numBdofs]; + + /// Column index for pressure + globalIndex pressureColIndex; + + /// Unit normal from rotation matrix + real64 unitNormal[3]; + + /// Local residual for displacement + real64 localRu[numUdofs]; + + /// Local residual for bubble + real64 localRb[numBdofs]; + + /// Jacobian dRu/dP + real64 dRudP[numUdofs]; + + /// Jacobian dRb/dP + real64 dRbdP[numBdofs]; + }; + + template< typename POLICY, + typename KERNEL_TYPE > + static + real64 + kernelLaunch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ) + { + return Base::template kernelLaunch< POLICY, KERNEL_TYPE >( numElems, kernelComponent ); + } + + /** + * @brief Setup for each element + */ + GEOS_HOST_DEVICE + inline + void setup( localIndex const k, + StackVariables & stack ) const + { + constexpr int shift = numNodesPerElem * 3; + + int permutation[numNodesPerElem]; + m_finiteElementSpace.getPermutation( permutation ); + + localIndex const kf0 = m_elemsToFaces[k][0]; + localIndex const kf1 = m_elemsToFaces[k][1]; + + for( localIndex a = 0; a < numNodesPerElem; ++a ) + { + localIndex const kn0 = m_faceToNodes( kf0, a ); + localIndex const kn1 = m_faceToNodes( kf1, a ); + for( int i = 0; i < 3; ++i ) + { + stack.dispEqnRowIndices[a * 3 + i] = m_dofNumber[kn0] + i - m_dofRankOffset; + stack.dispEqnRowIndices[shift + a * 3 + i] = m_dofNumber[kn1] + i - m_dofRankOffset; + stack.X[a][i] = m_X[m_faceToNodes( kf0, permutation[a] )][i]; + } + } + + for( int i = 0; i < 3; ++i ) + { + stack.unitNormal[i] = m_rotationMatrix( k, i, 0 ); + } + + for( int i = 0; i < 3; ++i ) + { + stack.bEqnRowIndices[i] = m_bDofNumber[kf0] + i - m_dofRankOffset; + stack.bEqnRowIndices[3 + i] = m_bDofNumber[kf1] + i - m_dofRankOffset; + } + + stack.pressureColIndex = m_pressureDofNumber[k]; + } + + /** + * @brief Complete - compute residuals and Jacobians + */ + GEOS_HOST_DEVICE + inline + real64 complete( localIndex const k, + StackVariables & stack ) const + { + // Compute residual and Jacobian for displacement DOFs: + // R_u = transpose(Atu) * unitNormal * pressure + // dR_u/dP = transpose(Atu) * unitNormal + // Note: This matches the existing AssemblePressureContribution kernel used for sequential coupling. + LvArray::tensorOps::Ri_eq_AjiBj< numUdofs, 3 >( stack.localRu, stack.localAtu, stack.unitNormal ); + LvArray::tensorOps::Ri_eq_AjiBj< numUdofs, 3 >( stack.dRudP, stack.localAtu, stack.unitNormal ); + + // Scale residual by pressure + LvArray::tensorOps::scale< numUdofs >( stack.localRu, m_pressure[k] ); + + // Compute residual and Jacobian for bubble DOFs: + // R_b = transpose(Atb) * unitNormal * pressure + // dR_b/dP = transpose(Atb) * unitNormal + LvArray::tensorOps::Ri_eq_AjiBj< numBdofs, 3 >( stack.localRb, stack.localAtb, stack.unitNormal ); + LvArray::tensorOps::Ri_eq_AjiBj< numBdofs, 3 >( stack.dRbdP, stack.localAtb, stack.unitNormal ); + + // Scale residual by pressure + LvArray::tensorOps::scale< numBdofs >( stack.localRb, m_pressure[k] ); + + // Assemble displacement contributions + for( localIndex i = 0; i < numUdofs; ++i ) + { + localIndex const dof = LvArray::integerConversion< localIndex >( stack.dispEqnRowIndices[i] ); + + if( dof < 0 || dof >= m_matrix.numRows() ) + continue; + + // Add residual + RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[dof], stack.localRu[i] ); + + // Add Jacobian + m_matrix.template addToRow< parallelDeviceAtomic >( dof, + &stack.pressureColIndex, + &stack.dRudP[i], + 1 ); + } + + // Assemble bubble contributions + for( localIndex i = 0; i < numBdofs; ++i ) + { + localIndex const dof = LvArray::integerConversion< localIndex >( stack.bEqnRowIndices[i] ); + + if( dof < 0 || dof >= m_matrix.numRows() ) + continue; + + // Add residual + RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[dof], stack.localRb[i] ); + + // Add Jacobian + m_matrix.template addToRow< parallelDeviceAtomic >( dof, + &stack.pressureColIndex, + &stack.dRbdP[i], + 1 ); + } + + return 0.0; + } + +protected: + + /// Pressure array + arrayView1d< real64 const > const m_pressure; + + /// Pressure DOF numbers + arrayView1d< globalIndex const > const m_pressureDofNumber; +}; + +/// Factory for AssembleForceResidualDerivativeWrtPressure kernel +using AssembleForceResidualDerivativeWrtPressureFactory = + finiteElement::InterfaceKernelFactory< AssembleForceResidualDerivativeWrtPressure, + arrayView1d< globalIndex const > const, + arrayView1d< globalIndex const > const, + globalIndex const, + CRSMatrixView< real64, globalIndex const > const, + arrayView1d< real64 > const, + real64 const, + arrayView1d< localIndex const > const, + arrayView1d< globalIndex const > const >; + + +/** + * @brief Kernel to assemble the Jacobian derivative dR_p/du and dR_p/db for fully implicit coupling. + * + * This kernel computes the derivative of the flow accumulation residual with respect to + * nodal displacement and bubble displacement, using the properly computed Atu/Atb matrices. + */ +template< typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +class AssembleFluidMassResidualDerivativeWrtDisplacement : + public solidMechanicsConformingContactKernels::ConformingContactKernelsBase< CONSTITUTIVE_TYPE, + FE_TYPE > +{ + +public: + /// Alias for the base class + using Base = solidMechanicsConformingContactKernels::ConformingContactKernelsBase< CONSTITUTIVE_TYPE, + FE_TYPE >; + + /// Number of nodes per element + static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem; + + /// The number of displacement dofs per element. + static constexpr int numUdofs = Base::numUdofs; + + /// The number of bubble dofs per element. + static constexpr int numBdofs = Base::numBdofs; + + using Base::m_elemsToFaces; + using Base::m_faceToNodes; + using Base::m_X; + using Base::m_finiteElementSpace; + using Base::m_dofNumber; + using Base::m_bDofNumber; + using Base::m_dofRankOffset; + using Base::m_rotationMatrix; + using Base::m_matrix; + using Base::m_rhs; + + /** + * @brief Constructor + */ + AssembleFluidMassResidualDerivativeWrtDisplacement( NodeManager const & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + FaceElementSubRegion & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + arrayView1d< globalIndex const > const uDofNumber, + arrayView1d< globalIndex const > const bDofNumber, + globalIndex const rankOffset, + CRSMatrixView< real64, globalIndex const > const inputMatrix, + arrayView1d< real64 > const inputRhs, + real64 const inputDt, + arrayView1d< localIndex const > const & faceElementList, + arrayView1d< globalIndex const > const & pressureDofNumber, + arrayView1d< real64 const > const & density, + arrayView1d< integer const > const & fractureState ): + Base( nodeManager, + edgeManager, + faceManager, + targetRegionIndex, + elementSubRegion, + finiteElementSpace, + inputConstitutiveType, + uDofNumber, + bDofNumber, + rankOffset, + inputMatrix, + inputRhs, + inputDt, + faceElementList ), + m_pressureDofNumber( pressureDofNumber ), + m_density( density ), + m_fractureState( fractureState ) + {} + + /** + * @copydoc finiteElement::InterfaceKernelBase::StackVariables + */ + struct StackVariables : public Base::StackVariables + { +public: + + GEOS_HOST_DEVICE + StackVariables(): + Base::StackVariables(), + dispColIndices{}, + bColIndices{}, + pressureRowIndex( 0 ), + unitNormal{}, + dRpdU{}, + dRpdB{} + {} + + /// Column indices for displacement + globalIndex dispColIndices[numUdofs]; + + /// Column indices for bubble + globalIndex bColIndices[numBdofs]; + + /// Row index for pressure equation + globalIndex pressureRowIndex; + + /// Unit normal from rotation matrix + real64 unitNormal[3]; + + /// Jacobian dRp/dU + real64 dRpdU[numUdofs]; + + /// Jacobian dRp/dB + real64 dRpdB[numBdofs]; + }; + + template< typename POLICY, + typename KERNEL_TYPE > + static + real64 + kernelLaunch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ) + { + return Base::template kernelLaunch< POLICY, KERNEL_TYPE >( numElems, kernelComponent ); + } + + /** + * @brief Setup for each element + */ + GEOS_HOST_DEVICE + inline + void setup( localIndex const k, + StackVariables & stack ) const + { + constexpr int shift = numNodesPerElem * 3; + + int permutation[numNodesPerElem]; + m_finiteElementSpace.getPermutation( permutation ); + + localIndex const kf0 = m_elemsToFaces[k][0]; + localIndex const kf1 = m_elemsToFaces[k][1]; + + for( localIndex a = 0; a < numNodesPerElem; ++a ) + { + localIndex const kn0 = m_faceToNodes( kf0, a ); + localIndex const kn1 = m_faceToNodes( kf1, a ); + for( int i = 0; i < 3; ++i ) + { + stack.dispColIndices[a * 3 + i] = m_dofNumber[kn0] + i; + stack.dispColIndices[shift + a * 3 + i] = m_dofNumber[kn1] + i; + stack.X[a][i] = m_X[m_faceToNodes( kf0, permutation[a] )][i]; + } + } + + for( int i = 0; i < 3; ++i ) + { + stack.unitNormal[i] = m_rotationMatrix( k, i, 0 ); + } + + for( int i = 0; i < 3; ++i ) + { + stack.bColIndices[i] = m_bDofNumber[kf0] + i; + stack.bColIndices[3 + i] = m_bDofNumber[kf1] + i; + } + + stack.pressureRowIndex = m_pressureDofNumber[k] - m_dofRankOffset; + } + + /** + * @brief Complete - compute Jacobians dRp/dU and dRp/dB + */ + GEOS_HOST_DEVICE + inline + real64 complete( localIndex const k, + StackVariables & stack ) const + { + // Only assemble for open fractures + if( m_fractureState[k] != fields::contact::FractureState::Open ) + { + return 0.0; + } + + // Compute Jacobian for pressure equation w.r.t. displacement: + // dR_p/du = density * unitNormal^T * Atu + // This is: dRpdU[j] = sum_i( unitNormal[i] * Atu[i][j] ) * density + // Use Ri_eq_AjiBj for transpose: R[i] = sum_j(A[j][i] * B[j]) + // Note: The ALM Atu matrix has opposite sign convention compared to LM's dAper/dU, + // so we negate the result to match the LM sign convention. + LvArray::tensorOps::Ri_eq_AjiBj< numUdofs, 3 >( stack.dRpdU, stack.localAtu, stack.unitNormal ); + LvArray::tensorOps::scale< numUdofs >( stack.dRpdU, -m_density[k] ); + + // Compute Jacobian for pressure equation w.r.t. bubble displacement: + // dR_p/db = density * unitNormal^T * Atb + // Note: Same sign convention fix as for displacement. + LvArray::tensorOps::Ri_eq_AjiBj< numBdofs, 3 >( stack.dRpdB, stack.localAtb, stack.unitNormal ); + LvArray::tensorOps::scale< numBdofs >( stack.dRpdB, -m_density[k] ); + + // Assemble into matrix + localIndex const localRow = LvArray::integerConversion< localIndex >( stack.pressureRowIndex ); + + if( localRow >= 0 && localRow < m_matrix.numRows() ) + { + // Add dRp/dU + m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( localRow, + stack.dispColIndices, + stack.dRpdU, + numUdofs ); + + // Add dRp/dB + m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( localRow, + stack.bColIndices, + stack.dRpdB, + numBdofs ); + } + + return 0.0; + } + +protected: + + /// Pressure DOF numbers + arrayView1d< globalIndex const > const m_pressureDofNumber; + + /// Fluid density + arrayView1d< real64 const > const m_density; + + /// Fracture state + arrayView1d< integer const > const m_fractureState; +}; + +/// Factory for AssembleFluidMassResidualDerivativeWrtDisplacement kernel +using AssembleFluidMassResidualDerivativeWrtDisplacementFactory = + finiteElement::InterfaceKernelFactory< AssembleFluidMassResidualDerivativeWrtDisplacement, + arrayView1d< globalIndex const > const, + arrayView1d< globalIndex const > const, + globalIndex const, + CRSMatrixView< real64, globalIndex const > const, + arrayView1d< real64 > const, + real64 const, + arrayView1d< localIndex const > const, + arrayView1d< globalIndex const > const, + arrayView1d< real64 const > const, + arrayView1d< integer const > const >; + + +/** + * @brief Kernel to compute d(aperture)/du and d(aperture)/db for each fracture element. + * + * This kernel computes the derivative of aperture with respect to displacement and bubble DOFs + * using the properly computed Atu/Atb matrices. The results are stored for later use in + * flux derivative computation. + * + * d(aperture)/du = (1/area) * unitNormal^T * Atu + * d(aperture)/db = (1/area) * unitNormal^T * Atb + * + * IMPORTANT: The ALM Atu/Atb matrices have opposite sign convention compared to the LM + * dAper/dU formulas. When using these pre-computed values, they must be negated to match + * the LM sign convention. The assembly in assembleFluidMassResidualDerivativeWrtDisplacement + * handles this negation. + */ +template< typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +class ComputeApertureDerivatives : + public solidMechanicsConformingContactKernels::ConformingContactKernelsBase< CONSTITUTIVE_TYPE, + FE_TYPE > +{ + +public: + /// Alias for the base class + using Base = solidMechanicsConformingContactKernels::ConformingContactKernelsBase< CONSTITUTIVE_TYPE, + FE_TYPE >; + + /// Number of nodes per element + static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem; + + /// The number of displacement dofs per element. + static constexpr int numUdofs = Base::numUdofs; + + /// The number of bubble dofs per element. + static constexpr int numBdofs = Base::numBdofs; + + using Base::m_elemsToFaces; + using Base::m_faceToNodes; + using Base::m_X; + using Base::m_finiteElementSpace; + using Base::m_dofNumber; + using Base::m_bDofNumber; + using Base::m_dofRankOffset; + using Base::m_rotationMatrix; + using Base::m_matrix; + using Base::m_rhs; + + /** + * @brief Constructor + */ + ComputeApertureDerivatives( NodeManager const & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + FaceElementSubRegion & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + arrayView1d< globalIndex const > const uDofNumber, + arrayView1d< globalIndex const > const bDofNumber, + globalIndex const rankOffset, + CRSMatrixView< real64, globalIndex const > const inputMatrix, + arrayView1d< real64 > const inputRhs, + real64 const inputDt, + arrayView1d< localIndex const > const & faceElementList, + arrayView2d< real64 > const & dAperturedU, + arrayView2d< real64 > const & dAperturedB ): + Base( nodeManager, + edgeManager, + faceManager, + targetRegionIndex, + elementSubRegion, + finiteElementSpace, + inputConstitutiveType, + uDofNumber, + bDofNumber, + rankOffset, + inputMatrix, + inputRhs, + inputDt, + faceElementList ), + m_area( elementSubRegion.getElementArea().toViewConst() ), + m_dAperturedU( dAperturedU ), + m_dAperturedB( dAperturedB ) + {} + + /** + * @copydoc finiteElement::InterfaceKernelBase::StackVariables + */ + struct StackVariables : public Base::StackVariables + { +public: + + GEOS_HOST_DEVICE + StackVariables(): + Base::StackVariables(), + unitNormal{}, + localDAperturedU{}, + localDAperturedB{} + {} + + /// Unit normal from rotation matrix + real64 unitNormal[3]; + + /// Local d(aperture)/dU + real64 localDAperturedU[numUdofs]; + + /// Local d(aperture)/dB + real64 localDAperturedB[numBdofs]; + }; + + template< typename POLICY, + typename KERNEL_TYPE > + static + real64 + kernelLaunch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ) + { + return Base::template kernelLaunch< POLICY, KERNEL_TYPE >( numElems, kernelComponent ); + } + + /** + * @brief Setup for each element + */ + GEOS_HOST_DEVICE + inline + void setup( localIndex const k, + StackVariables & stack ) const + { + int permutation[numNodesPerElem]; + m_finiteElementSpace.getPermutation( permutation ); + + localIndex const kf0 = m_elemsToFaces[k][0]; + + for( localIndex a = 0; a < numNodesPerElem; ++a ) + { + for( int i = 0; i < 3; ++i ) + { + stack.X[a][i] = m_X[m_faceToNodes( kf0, permutation[a] )][i]; + } + } + + for( int i = 0; i < 3; ++i ) + { + stack.unitNormal[i] = m_rotationMatrix( k, i, 0 ); + } + } + + /** + * @brief Complete - compute and store d(aperture)/dU and d(aperture)/dB + */ + GEOS_HOST_DEVICE + inline + real64 complete( localIndex const k, + StackVariables & stack ) const + { + real64 const invArea = 1.0 / m_area[k]; + + // Compute d(aperture)/dU = (1/area) * unitNormal^T * Atu + // This gives: dAperturedU[j] = (1/area) * sum_i(unitNormal[i] * Atu[i][j]) + // Use Ri_eq_AjiBj for transpose: R[i] = sum_j(A[j][i] * B[j]) + LvArray::tensorOps::Ri_eq_AjiBj< numUdofs, 3 >( stack.localDAperturedU, stack.localAtu, stack.unitNormal ); + LvArray::tensorOps::scale< numUdofs >( stack.localDAperturedU, invArea ); + + // Compute d(aperture)/dB = (1/area) * unitNormal^T * Atb + LvArray::tensorOps::Ri_eq_AjiBj< numBdofs, 3 >( stack.localDAperturedB, stack.localAtb, stack.unitNormal ); + LvArray::tensorOps::scale< numBdofs >( stack.localDAperturedB, invArea ); + + // Store the results + for( localIndex i = 0; i < numUdofs; ++i ) + { + m_dAperturedU( k, i ) = stack.localDAperturedU[i]; + } + for( localIndex i = 0; i < numBdofs; ++i ) + { + m_dAperturedB( k, i ) = stack.localDAperturedB[i]; + } + + return 0.0; + } + +protected: + + /// Element area + arrayView1d< real64 const > const m_area; + + /// Storage for d(aperture)/dU + arrayView2d< real64 > const m_dAperturedU; + + /// Storage for d(aperture)/dB + arrayView2d< real64 > const m_dAperturedB; +}; + +/// Factory for ComputeApertureDerivatives kernel +using ComputeApertureDerivativesFactory = + finiteElement::InterfaceKernelFactory< ComputeApertureDerivatives, + arrayView1d< globalIndex const > const, + arrayView1d< globalIndex const > const, + globalIndex const, + CRSMatrixView< real64, globalIndex const > const, + arrayView1d< real64 > const, + real64 const, + arrayView1d< localIndex const > const, + arrayView2d< real64 > const, + arrayView2d< real64 > const >; + +} // namespace poromechanicsALMKernels + + +namespace poromechanicsMatrixBubbleKernels +{ + +/** + * @brief Kernel to compute the contribution of matrix cell pressure on bubble DOFs + * with full Jacobian for fully-implicit poromechanics coupling. + * + * This kernel is based on SolidMechanicsPressureFaceBubbleKernels but adds the Jacobian + * term dR_b/dP_matrix needed for fully implicit coupling. + * + * Physics: + * R_b = -∫ B_b^T * (biot * p * I) * detJ dΩ + * dR_b/dP = -∫ B_b^T * (biot * I) * detJ dΩ + * + * where B_b is the strain-displacement matrix for bubble functions. + */ +template< typename SUBREGION_TYPE, + typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +class MatrixPressureBubbleKernels : + public finiteElement::ImplicitKernelBase< SUBREGION_TYPE, + CONSTITUTIVE_TYPE, + FE_TYPE, + 3, + 3 > +{ +public: + + /// Alias for the base class + using Base = finiteElement::ImplicitKernelBase< SUBREGION_TYPE, + CONSTITUTIVE_TYPE, + FE_TYPE, + 3, + 3 >; + + /// Number of nodes per element + static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem; + + /// Compile time value for the number of faces per element. + static constexpr int numFacesPerElem = FE_TYPE::numFaces; + + /// Compile time value for the number of quadrature points per element. + static constexpr int numQuadraturePointsPerElem = FE_TYPE::numQuadraturePoints; + + using Base::m_elemsToNodes; + using Base::m_finiteElementSpace; + using Base::m_constitutiveUpdate; + using Base::m_dofNumber; + using Base::m_dofRankOffset; + using Base::m_matrix; + using Base::m_rhs; + using Base::m_dt; + + /** + * @brief Constructor + */ + MatrixPressureBubbleKernels( NodeManager const & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + SUBREGION_TYPE const & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + arrayView1d< globalIndex const > const uDofNumber, + arrayView1d< globalIndex const > const bDofNumber, + globalIndex const rankOffset, + CRSMatrixView< real64, globalIndex const > const inputMatrix, + arrayView1d< real64 > const inputRhs, + real64 const inputDt, + string const & pressureDofKey ): + Base( nodeManager, + edgeManager, + faceManager, + targetRegionIndex, + elementSubRegion, + finiteElementSpace, + inputConstitutiveType, + uDofNumber, + rankOffset, + inputMatrix, + inputRhs, + inputDt ), + m_X( nodeManager.referencePosition()), + m_bDofNumber( bDofNumber ), + m_pDofNumber( elementSubRegion.template getReference< array1d< globalIndex > >( pressureDofKey ) ), + m_bubbleElems( elementSubRegion.bubbleElementsList() ), + m_elemsToFaces( elementSubRegion.faceElementsList() ), + m_pressure( elementSubRegion.template getField< fields::flow::pressure >().toViewConst() ) + {} + + //*************************************************************************** + + /** + * @copydoc finiteElement::ImplicitKernelBase::StackVariables + */ + struct StackVariables + { +public: + + /// The number of displacement dofs per element. + static constexpr int numUdofs = numNodesPerElem * 3; + + /// The number of bubble dofs per element (all faces). + static constexpr int numBubbleUdofs = numFacesPerElem * 3; + + /// The number of pressure dofs per element. + static constexpr int numPdofs = 1; + + /** + * Default constructor + */ + GEOS_HOST_DEVICE + StackVariables(): + bEqnRowIndices{}, + pColIndex( 0 ), + localRb{}, + localdRbdP{}, + X{ {} }, + pLocal{} + {} + + /// C-array storage for the element local row degrees of freedom (bubble). + globalIndex bEqnRowIndices[3]; + + /// Column index for pressure DOF + globalIndex pColIndex; + + /// C-array storage for the element local Rb residual vector. + real64 localRb[numBubbleUdofs]; + + /// C-array storage for the element local dRb/dP Jacobian. + real64 localdRbdP[numBubbleUdofs]; + + /// local nodal coordinates + real64 X[ numNodesPerElem ][ 3 ]; + + /// local pressure + real64 pLocal[numPdofs]; + + }; + + //*************************************************************************** + + /** + * @copydoc ::geos::finiteElement::KernelBase::kernelLaunch + * + * @detail it uses the kernelLaunch interface of KernelBase but it only launches the kernel + * on the set of elements that have bubble dof within the subregion. + */ + template< typename POLICY, + typename KERNEL_TYPE > + static + real64 + kernelLaunch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ) + { + GEOS_MARK_FUNCTION; + GEOS_UNUSED_VAR( numElems ); + + // Define a RAJA reduction variable to get the maximum residual contribution. + RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxResidual( 0 ); + + forAll< POLICY >( kernelComponent.m_bubbleElems.size(), + [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + typename KERNEL_TYPE::StackVariables stack; + + kernelComponent.setup( i, stack ); + for( integer q=0; q( dBubbleNdX, 0 ); + + real64 const detJ = m_finiteElementSpace.calcGradFaceBubbleN( q, stack.X, dBubbleNdX ); + + real64 biotCoefficient; + m_constitutiveUpdate.getBiotCoefficient( k, biotCoefficient ); + + // Assemble the strain-displacement matrix for bubble functions: B_b + real64 strainBubbleMatrix[6][nBubbleUdof]; + solidMechanicsConformingContactKernelsHelper:: + assembleStrainOperator< 6, nBubbleUdof, numFacesPerElem >( strainBubbleMatrix, dBubbleNdX ); + + // Compute biotPressure = -biot * p * I (Voigt notation: [1,1,1,0,0,0]) + real64 biotPressure[6] = {0}; + LvArray::tensorOps::symAddIdentity< 3 >( biotPressure, -biotCoefficient * stack.pLocal[0] ); + + // Compute biotIdentity = -biot * I (for Jacobian) + real64 biotIdentity[6] = {0}; + LvArray::tensorOps::symAddIdentity< 3 >( biotIdentity, -biotCoefficient ); + + // Residual contribution: R_b += B_b^T * biotPressure * detJ + // (note: scaledAdd with -detJ because biotPressure already has negative sign) + real64 Rb_gauss[nBubbleUdof]; + LvArray::tensorOps::Ri_eq_AjiBj< nBubbleUdof, 6 >( Rb_gauss, strainBubbleMatrix, biotPressure ); + LvArray::tensorOps::scaledAdd< nBubbleUdof >( stack.localRb, Rb_gauss, -detJ ); + + // Jacobian contribution: dR_b/dP += B_b^T * biotIdentity * detJ + real64 dRbdP_gauss[nBubbleUdof]; + LvArray::tensorOps::Ri_eq_AjiBj< nBubbleUdof, 6 >( dRbdP_gauss, strainBubbleMatrix, biotIdentity ); + LvArray::tensorOps::scaledAdd< nBubbleUdof >( stack.localdRbdP, dRbdP_gauss, -detJ ); + } + + /** + * @brief Complete - assemble residual and Jacobian into global system + */ + GEOS_HOST_DEVICE + inline + real64 complete( localIndex const kk, + StackVariables & stack ) const + { + localIndex const parentFaceIndex = m_elemsToFaces[kk][1]; + + // Extract only the components corresponding to the parent face + // on which the bubble function was applied. + real64 localRb[3]; + real64 localdRbdP[3]; + for( localIndex i = 0; i < 3; ++i ) + { + localRb[i] = stack.localRb[parentFaceIndex*3+i]; + localdRbdP[i] = stack.localdRbdP[parentFaceIndex*3+i]; + } + + for( localIndex i=0; i < 3; ++i ) + { + localIndex const dof = LvArray::integerConversion< localIndex >( stack.bEqnRowIndices[ i ] ); + + if( dof < 0 || dof >= m_matrix.numRows() ) + continue; + + // Add residual + RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[dof], localRb[i] ); + + // Add Jacobian dR_b/dP + m_matrix.template addToRow< parallelDeviceAtomic >( dof, + &stack.pColIndex, + &localdRbdP[i], + 1 ); + } + + return 0.0; + } + +protected: + + /// The array containing the nodal position array. + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const m_X; + + /// The global degree of freedom number of bubble + arrayView1d< globalIndex const > const m_bDofNumber; + + /// The global degree of freedom number of pressure + arrayView1d< globalIndex const > const m_pDofNumber; + + /// The array containing the list of bubble elements. + arrayView1d< localIndex const > const m_bubbleElems; + + /// The array containing the element to bubble face map. + arrayView2d< localIndex const > const m_elemsToFaces; + + /// The array containing the pressure of each element. + arrayView1d< real64 const > const m_pressure; + +}; + +/// The factory used to construct MatrixPressureBubbleKernels. +using MatrixPressureBubbleFactory = finiteElement::KernelFactory< MatrixPressureBubbleKernels, + arrayView1d< globalIndex const > const, + arrayView1d< globalIndex const > const, + globalIndex const, + CRSMatrixView< real64, globalIndex const > const, + arrayView1d< real64 > const, + real64 const, + string const & >; + +} // namespace poromechanicsMatrixBubbleKernels + +} // namespace geos + +#endif /* GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALMKERNELS_HPP_ */ From dedd57ee45b4298fa41dc26952a0aeaec19d1af6 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 4 Feb 2026 14:52:11 +0100 Subject: [PATCH 03/22] missing guards --- .../SinglePhasePoromechanicsConformingFracturesALM.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp index ab86e0c8054..e9fb7db8bd6 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -17,8 +17,8 @@ * @file SinglePhasePoromechanicsConformingFractures.hpp */ -#ifndef GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_ -#define GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_ +#ifndef GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_ +#define GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_ #include "linearAlgebra/interfaces/hypre/HypreMGR.hpp" @@ -131,4 +131,4 @@ class SinglePhasePoromechanicsConformingFracturesALM : public MGRStrategyBase< 2 } // namespace geos -#endif /*GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_*/ \ No newline at end of file +#endif /*GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_*/ \ No newline at end of file From 9143dd14822f7a9910118867fc79ea1ac3603b95 Mon Sep 17 00:00:00 2001 From: npillardou Date: Thu, 5 Feb 2026 13:16:03 +0100 Subject: [PATCH 04/22] Add fix to fracture state in case of high cohesion --- .../kernels/SolidMechanicsALMKernelsBase.hpp | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp index 65841f9c66c..39e515378d2 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp @@ -144,13 +144,24 @@ struct UpdateStateKernel if( fractureState[k] == fields::contact::FractureState::Open ) { - LvArray::tensorOps::fill< 3 >( localTractionNew, 0.0 ); } else if( LvArray::math::abs( localTractionNew[ 0 ] ) < normalTractionTolerance[k] ) { - LvArray::tensorOps::fill< 3 >( localTractionNew, 0.0 ); - fractureState[k] = fields::contact::FractureState::Slip; + // When normal traction is very small, check if cohesion allows stick state + // before deciding to transition to slip (fix for high cohesion materials) + real64 dLimitTau_dTraction = 0.0; + real64 const limitTau = contactWrapper.computeLimitTangentialTractionNorm( 0.0, dLimitTau_dTraction ); + real64 const currentTau = LvArray::math::sqrt( localTractionNew[1] * localTractionNew[1] + + localTractionNew[2] * localTractionNew[2] ); + + if( currentTau > limitTau ) + { + // Tangential traction exceeds cohesion limit: transition to slip + LvArray::tensorOps::fill< 3 >( localTractionNew, 0.0 ); + fractureState[k] = fields::contact::FractureState::Slip; + } + // else: cohesion allows stick state, keep traction as computed by updateTraction() } LvArray::tensorOps::copy< 3 >( traction[k], localTractionNew ); From 616626be592d1b9c57f5e9a4635c62f466bb147d Mon Sep 17 00:00:00 2001 From: npillardou Date: Thu, 5 Feb 2026 13:35:21 +0100 Subject: [PATCH 05/22] Correct traction in limit --- .../contact/kernels/SolidMechanicsALMKernelsBase.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp index 39e515378d2..12c61d17ccd 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp @@ -151,7 +151,7 @@ struct UpdateStateKernel // When normal traction is very small, check if cohesion allows stick state // before deciding to transition to slip (fix for high cohesion materials) real64 dLimitTau_dTraction = 0.0; - real64 const limitTau = contactWrapper.computeLimitTangentialTractionNorm( 0.0, dLimitTau_dTraction ); + real64 const limitTau = contactWrapper.computeLimitTangentialTractionNorm( localTractionNew[ 0 ], dLimitTau_dTraction ); real64 const currentTau = LvArray::math::sqrt( localTractionNew[1] * localTractionNew[1] + localTractionNew[2] * localTractionNew[2] ); From 24cf0efea26da40ac1d3463c0d84d6e0bd90f285 Mon Sep 17 00:00:00 2001 From: npillardou Date: Thu, 5 Feb 2026 16:22:49 +0100 Subject: [PATCH 06/22] Add fix on synchronization for ghost issue remaining --- ...lidMechanicsAugmentedLagrangianContact.cpp | 24 +++++++++++++------ 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index 2d18a5782a0..b0e00a3e322 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -912,6 +912,21 @@ void SolidMechanicsAugmentedLagrangianContact::applySystemSolution( DofManager c contact::incrementalBubbleDisplacement::key(), scalingFactor ); + // Synchronize bubble displacements before computing displacement jump + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + string_array const & ) + { + FieldIdentifiers fieldsToBeSync; + fieldsToBeSync.addFields( FieldLocation::Face, + { contact::incrementalBubbleDisplacement::key(), + contact::totalBubbleDisplacement::key() } ); + + CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, + mesh, + domain.getNeighbors(), + true ); + } ); // Loop for updating the displacement jump forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshName, @@ -962,19 +977,14 @@ void SolidMechanicsAugmentedLagrangianContact::applySystemSolution( DofManager c } ); } ); + // Synchronize displacement jump after computation forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, string_array const & ) - { FieldIdentifiers fieldsToBeSync; - fieldsToBeSync.addFields( FieldLocation::Face, - { contact::incrementalBubbleDisplacement::key(), - contact::totalBubbleDisplacement::key() } ); - - fieldsToBeSync.addElementFields( { contact::traction::key(), - contact::dispJump::key(), + fieldsToBeSync.addElementFields( { contact::dispJump::key(), contact::deltaDispJump::key() }, { getUniqueFractureRegionName() } ); From 68d789cc665fcb3936010366a3a319c3370d2059 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Thu, 5 Feb 2026 13:28:21 -0600 Subject: [PATCH 07/22] uncrustify --- .../contact/kernels/SolidMechanicsALMKernelsBase.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp index 12c61d17ccd..7c237a8231c 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp @@ -153,7 +153,7 @@ struct UpdateStateKernel real64 dLimitTau_dTraction = 0.0; real64 const limitTau = contactWrapper.computeLimitTangentialTractionNorm( localTractionNew[ 0 ], dLimitTau_dTraction ); real64 const currentTau = LvArray::math::sqrt( localTractionNew[1] * localTractionNew[1] + - localTractionNew[2] * localTractionNew[2] ); + localTractionNew[2] * localTractionNew[2] ); if( currentTau > limitTau ) { From c381048936ba462d37d56c6971662d2eaddea1bd Mon Sep 17 00:00:00 2001 From: npillardou Date: Wed, 4 Feb 2026 11:12:15 +0100 Subject: [PATCH 08/22] Update with other branch --- ...asePoromechanicsConformingFracturesALM.cpp | 1110 ++++++++++++++++- ...asePoromechanicsConformingFracturesALM.hpp | 63 +- ...mechanicsConformingFracturesALMKernels.hpp | 1079 ++++++++++++++++ 3 files changed, 2219 insertions(+), 33 deletions(-) create mode 100644 src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp index ff7d9662a60..65da565f2b4 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp @@ -19,7 +19,19 @@ #include "SinglePhasePoromechanicsConformingFracturesALM.hpp" +#include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" +#include "constitutive/solid/CoupledSolidBase.hpp" +#include "constitutive/solid/PorousSolid.hpp" +#include "constitutive/contact/HydraulicApertureBase.hpp" +#include "constitutive/contact/HydraulicApertureRelationSelector.hpp" #include "finiteVolume/FluxApproximationBase.hpp" +#include "mesh/SurfaceElementRegion.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp" +#include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp" +#include "physicsSolvers/solidMechanics/contact/ContactFields.hpp" namespace geos { @@ -40,10 +52,13 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupCouplin { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain, dofManager ); - - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + // 1. Poromechanical coupling in the bulk (from base class) + Base::setupCoupling( domain, dofManager ); + // 2. Pressure - bubble displacement coupling in the fracture + dofManager.addCoupling( m_pressureKey, + contact::totalBubbleDisplacement::key(), + DofManager::Connector::Elem ); } @@ -55,13 +70,87 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupSystem( ParallelVector & solution, bool const setSparsity ) { - GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain, dofManager, localMatrix, rhs, solution, setSparsity ); + if( this->m_precond ) + { + this->m_precond->clear(); + } + + // Set domain on DofManager + dofManager.setDomain( domain ); + + // Initialize ALM contact solver internal data structures + // These must be called before DOF setup and assembly + this->solidMechanicsSolver()->createFaceTypeList( domain ); + this->solidMechanicsSolver()->updateStickSlipList( domain ); + this->solidMechanicsSolver()->createBubbleCellList( domain ); + + // Setup DOFs for all sub-solvers and coupling (uses PoromechanicsSolver::setupDofs) + // This adds: displacement DOFs, bubble DOFs, pressure DOFs, and all couplings + this->setupDofs( domain, dofManager ); + + // Reorder DOFs to optimize matrix structure + dofManager.reorderByRank(); + + if( setSparsity ) + { + // Get sparsity pattern from flow solver (could include wells) + SparsityPattern< globalIndex > patternOriginal; + this->flowSolver()->setSparsityPattern( domain, dofManager, localMatrix, patternOriginal ); + + // Get the original row lengths + array1d< localIndex > rowLengths( patternOriginal.numRows() ); + for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) + { + rowLengths[localRow] = patternOriginal.numNonZeros( localRow ); + } + + // Add the number of nonzeros induced by coupling + addTransmissibilityCouplingNNZ( domain, dofManager, rowLengths.toView() ); + addPressureForceCouplingNNZ( domain, dofManager, rowLengths.toView() ); + addMatrixPressureBubbleCouplingNNZ( domain, dofManager, rowLengths.toView() ); + + // Create a new pattern with enough capacity for coupled matrix + SparsityPattern< globalIndex > pattern; + pattern.resizeFromRowCapacities< parallelHostPolicy >( patternOriginal.numRows(), + patternOriginal.numColumns(), + rowLengths.data() ); + + // Copy the original nonzeros + for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) + { + globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow ) ); + } + + // Add the nonzeros from coupling + addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView() ); + addPressureForceCouplingPattern( domain, dofManager, pattern.toView() ); + addMatrixPressureBubbleCouplingPattern( domain, dofManager, pattern.toView() ); + + // Assemble the full sparsity pattern using the solid mechanics solver + this->solidMechanicsSolver()->setSparsityPattern( domain, dofManager, localMatrix, pattern ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + // Set up the derivative flux residual matrix + setUpDflux_dApertureMatrix( domain, dofManager, localMatrix ); + // Assimilate the sparsity pattern into the local matrix + localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) ); + } + + localMatrix.setName( this->getName() + "/matrix" ); + + rhs.setName( this->getName() + "/rhs" ); + rhs.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); + + solution.setName( this->getName() + "/solution" ); + solution.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); + + if( !this->m_precond && this->m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) + { + this->m_precond = this->createPreconditioner( domain ); + } } template< typename FLOW_SOLVER > @@ -72,12 +161,25 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleSyst CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { - GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( time_n, dt, domain, dofManager, localMatrix, localRhs ); + // Synchronize fracture state + this->solidMechanicsSolver()->synchronizeFractureState( domain ); + + // Assemble element-based contributions (mechanics + flow accumulation) + assembleElementBasedContributions( time_n, dt, domain, dofManager, localMatrix, localRhs ); + + // Assemble flux terms and get dFluidResidual/dAperture + this->flowSolver()->assembleHydrofracFluxTerms( time_n, + dt, + domain, + dofManager, + localMatrix, + localRhs, + getDerivativeFluxResidual_dNormalJump() ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + // Assemble coupling terms (must be after flux assembly) + assembleCouplingTerms( time_n, dt, domain, dofManager, localMatrix, localRhs ); } template< typename FLOW_SOLVER > @@ -90,10 +192,23 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleElem { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( time_n, dt, domain, dofManager, localMatrix, localRhs ); + // Assemble poromechanics terms (from base class) + Base::assembleElementBasedTerms( time_n, dt, domain, dofManager, localMatrix, localRhs ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + // Flow accumulation for fractures + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const, + FaceElementSubRegion const & subRegion ) + { + this->flowSolver()->accumulationAssemblyLaunch( dofManager, subRegion, localMatrix, localRhs ); + } ); + } ); + // Assemble contact terms (ALM) - note: assembleContact requires time and dt + this->solidMechanicsSolver()->assembleContact( time_n, dt, domain, dofManager, localMatrix, localRhs ); } template< typename FLOW_SOLVER > @@ -105,27 +220,164 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleCoup arrayView1d< real64 > const & localRhs ) { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain, dofManager, localMatrix, localRhs, time_n, dt ); + GEOS_UNUSED_VAR( time_n, dt ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + // These steps must occur after the fluxes are assembled because that's when DerivativeFluxResidual_dAperture is filled. + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshName, + MeshLevel const & mesh, + string_array const & regionNames ) + { + // Assemble Force Residual w.r.t. pressure (Aup) - fracture pressure contribution + assembleForceResidualDerivativeWrtPressure( meshName, mesh, regionNames, dofManager, localMatrix, localRhs ); + + // Assemble Fluid mass residual w.r.t. displacement (Apu) + assembleFluidMassResidualDerivativeWrtDisplacement( meshName, mesh, regionNames, dofManager, localMatrix, localRhs ); + } ); + + // Assemble matrix cell pressure contribution on bubble DOFs (Abp_matrix) + // This must be outside the lambda because it uses regionBasedKernelApplication + assembleMatrixPressureBubbleContribution( dt, const_cast< DomainPartition & >( domain ), dofManager, localMatrix, localRhs ); } template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::updateState( DomainPartition & domain ) { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + + // Call base poromechanics update + Base::updateState( domain ); + + // Need to call solid mechanics update separately to compute face displacement jump + this->solidMechanicsSolver()->updateState( domain ); + + // Remove the contribution of the hydraulic aperture from the stencil weights + this->flowSolver()->prepareStencilWeights( domain ); + + // Update hydraulic aperture and fracture permeability + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + ElementRegionManager & elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, + [&]( localIndex const, + FaceElementSubRegion & subRegion ) + { + arrayView2d< real64 const > const dispJump = subRegion.getField< contact::dispJump >(); + arrayView1d< real64 const > const area = subRegion.getElementArea(); + arrayView1d< real64 const > const volume = subRegion.getElementVolume(); + arrayView2d< real64 const > const fractureTraction = subRegion.getField< contact::traction >(); + arrayView1d< real64 const > const pressure = subRegion.getField< flow::pressure >(); + arrayView1d< real64 const > const oldHydraulicAperture = subRegion.getField< flow::aperture0 >(); + + arrayView1d< real64 > const aperture = subRegion.getElementAperture(); + arrayView1d< real64 > const hydraulicAperture = subRegion.getField< flow::hydraulicAperture >(); + arrayView1d< real64 > const deltaVolume = subRegion.getField< flow::deltaVolume >(); + arrayView1d< integer > const & fractureState = subRegion.getField< contact::fractureState >(); + + string const & hydraulicApertureRelationName = subRegion.template getReference< string >( Base::viewKeyStruct::hydraulicApertureRelationNameString() ); + HydraulicApertureBase const & hydraulicApertureModel = this->template getConstitutiveModel< HydraulicApertureBase >( subRegion, hydraulicApertureRelationName ); + + string const porousSolidName = subRegion.getReference< string >( FlowSolverBase::viewKeyStruct::solidNamesString() ); + CoupledSolidBase & porousSolid = subRegion.template getConstitutiveModel< CoupledSolidBase >( porousSolidName ); + + constitutiveUpdatePassThru( hydraulicApertureModel, [&] ( auto & castedHydraulicAperture ) + { + using HydraulicApertureType = TYPEOFREF( castedHydraulicAperture ); + typename HydraulicApertureType::KernelWrapper hydraulicApertureWrapper = castedHydraulicAperture.createKernelWrapper(); + + ConstitutivePassThru< CompressibleSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid ) + { + typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousMaterialWrapper = castedPorousSolid.createKernelUpdates(); + + poromechanicsFracturesKernels::StateUpdateKernel:: + launch< parallelDevicePolicy<> >( subRegion.size(), + porousMaterialWrapper, + hydraulicApertureWrapper, + dispJump, + pressure, + area, + volume, + deltaVolume, + aperture, + oldHydraulicAperture, + hydraulicAperture, + fractureTraction, + fractureState ); + } ); + } ); + } ); + } ); + + // Update the stencil weights using the updated hydraulic aperture + this->flowSolver()->updateStencilWeights( domain ); } template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: setUpDflux_dApertureMatrix( DomainPartition & domain, - DofManager const & GEOS_UNUSED_PARAM( dofManager ), + DofManager const & dofManager, CRSMatrix< real64, globalIndex > & localMatrix ) { - GEOS_UNUSED_VAR( domain, localMatrix ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + GEOS_UNUSED_VAR( dofManager ); + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & regionNames ) + { + std::unique_ptr< CRSMatrix< real64, localIndex > > & derivativeFluxResidual_dAperture = getRefDerivativeFluxResidual_dAperture(); + + { + // Calculate number of fracture elements + localIndex numRows = 0; + mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, + [&]( localIndex const, FaceElementSubRegion const & subRegion ) + { + numRows += subRegion.size(); + } ); + + // Number of columns (derivatives) = number of fracture elements + localIndex numCol = numRows; + + derivativeFluxResidual_dAperture = std::make_unique< CRSMatrix< real64, localIndex > >( numRows, numCol ); + derivativeFluxResidual_dAperture->setName( this->getName() + "/derivativeFluxResidual_dAperture" ); + + derivativeFluxResidual_dAperture->reserveNonZeros( localMatrix.numNonZeros() ); + localIndex maxRowSize = -1; + for( localIndex row = 0; row < localMatrix.numRows(); ++row ) + { + localIndex const rowSize = localMatrix.numNonZeros( row ); + maxRowSize = maxRowSize > rowSize ? maxRowSize : rowSize; + } + + for( localIndex row = 0; row < numRows; ++row ) + { + derivativeFluxResidual_dAperture->reserveNonZeros( row, maxRowSize ); + } + } + + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); + + fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) + { + localIndex const numFluxElems = stencil.stencilSize( iconn ); + typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); + + for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) + { + for( localIndex k1 = 0; k1 < numFluxElems; ++k1 ) + { + derivativeFluxResidual_dAperture->insertNonZero( sei[iconn][k0], sei[iconn][k1], 0.0 ); + } + } + } + } ); + } ); } template< typename FLOW_SOLVER > @@ -136,9 +388,62 @@ addTransmissibilityCouplingNNZ( DomainPartition const & domain, { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain, dofManager, rowLengths ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & ) + { + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const flowDofKey = dofManager.getKey( m_pressureKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); + + fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) + { + localIndex const numFluxElems = stencil.stencilSize( iconn ); + typename SurfaceElementStencil::IndexContainerViewConstType const & seri = stencil.getElementRegionIndices(); + typename SurfaceElementStencil::IndexContainerViewConstType const & sesri = stencil.getElementSubRegionIndices(); + typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); + + FaceElementSubRegion const & elementSubRegion = + elemManager.getRegion( seri[iconn][0] ).getSubRegion< FaceElementSubRegion >( sesri[iconn][0] ); + + ArrayOfArraysView< localIndex const > const elemsToNodes = elementSubRegion.nodeList().toViewConst(); + + arrayView1d< globalIndex const > const faceElementDofNumber = + elementSubRegion.getReference< array1d< globalIndex > >( flowDofKey ); + + for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) + { + globalIndex const activeFlowDOF = faceElementDofNumber[sei[iconn][k0]]; + globalIndex const rowNumber = activeFlowDOF - rankOffset; + if( rowNumber >= 0 && rowNumber < rowLengths.size() ) + { + for( localIndex k1 = 0; k1 < numFluxElems; ++k1 ) + { + // The coupling with the nodal displacements of the cell itself has already been added by the dofManager + // so we only add the coupling with the nodal displacements of the neighbors. + if( k1 != k0 ) + { + localIndex const numNodesPerElement = elemsToNodes[sei[iconn][k1]].size(); + // Nodal displacement DOFs (3 per node, 2 faces) + rowLengths[rowNumber] += 3 * numNodesPerElement; + // Bubble DOFs (3 per face, 2 faces = 6 total) + rowLengths[rowNumber] += 6; + } + } + } + } + } + } ); + } ); } template< typename FLOW_SOLVER > @@ -149,38 +454,783 @@ addTransmissibilityCouplingPattern( DomainPartition const & domain, { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain, dofManager, pattern ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + using namespace contact; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + string const flowDofKey = dofManager.getKey( m_pressureKey ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const & + bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + + // Get the finite volume method used to compute the fluxes + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fvDiscretization = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); + + SurfaceElementRegion const & fractureRegion = + elemManager.getRegion< SurfaceElementRegion >( this->solidMechanicsSolver()->getUniqueFractureRegionName() ); + FaceElementSubRegion const & fractureSubRegion = + fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); + + GEOS_ERROR_IF( !fractureSubRegion.hasWrapper( flow::pressure::key() ), + this->getDataContext() << ": The fracture subregion must contain pressure field." ); + + arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst(); + + arrayView1d< globalIndex const > const & + flowDofNumber = fractureSubRegion.getReference< globalIndex_array >( flowDofKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + fvDiscretization.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + forAll< serialPolicy >( stencil.size(), [=] ( localIndex const iconn ) + { + localIndex const numFluxElems = stencil.stencilSize( iconn ); + + // A fracture connector has to be an edge shared by two faces + if( numFluxElems == 2 ) + { + typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); + + // First index: face element. Second index: node + for( localIndex kf = 0; kf < 2; ++kf ) + { + // Set row DOF index + // Note that the 1-kf index is intentional, as this is coupling the pressure of one face cell + // to the nodes of the adjacent cell + localIndex const rowIndex = flowDofNumber[sei[iconn][1 - kf]] - rankOffset; + + if( rowIndex >= 0 && rowIndex < pattern.numRows() ) + { + // Get fracture, face and region/subregion/element indices (for elements on both sides) + localIndex const fractureIndex = sei[iconn][kf]; + + // Get the number of nodes + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elem2dToFaces[fractureIndex][0] ); + + // Loop over the two sides of each fracture element + for( localIndex kf1 = 0; kf1 < 2; ++kf1 ) + { + localIndex const faceIndex = elem2dToFaces[fractureIndex][kf1]; + + // Save the list of DOF associated with nodes (displacement) + for( localIndex a = 0; a < numNodesPerFace; ++a ) + { + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const colIndex = dispDofNumber[faceToNodeMap( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i ); + pattern.insertNonZero( rowIndex, colIndex ); + } + } + + // Save the list of DOF associated with bubble displacement + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const colIndex = bubbleDofNumber[faceIndex] + LvArray::integerConversion< globalIndex >( i ); + pattern.insertNonZero( rowIndex, colIndex ); + } + } + } + } + } + } ); + } ); + } ); +} + +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: +addPressureForceCouplingNNZ( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< localIndex > const & rowLengths ) const +{ + GEOS_MARK_FUNCTION; + + using namespace contact; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const & + bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + + globalIndex const rankOffset = dofManager.rankOffset(); + + string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); + SurfaceElementRegion const & fractureRegion = + elemManager.getRegion< SurfaceElementRegion >( fractureRegionName ); + FaceElementSubRegion const & fractureSubRegion = + fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); + + arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst(); + + // For each fracture element, add NNZ for (displacement_row, pressure_col) and (bubble_row, pressure_col) + forAll< serialPolicy >( fractureSubRegion.size(), [=, &rowLengths] ( localIndex const kfe ) + { + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elem2dToFaces[kfe][0] ); + + // For displacement DOFs: add 1 pressure column per displacement DOF row + for( localIndex kf = 0; kf < 2; ++kf ) + { + localIndex const faceIndex = elem2dToFaces[kfe][kf]; + for( localIndex a = 0; a < numNodesPerFace; ++a ) + { + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowNumber = dispDofNumber[faceToNodeMap( faceIndex, a )] + i - rankOffset; + if( rowNumber >= 0 && rowNumber < rowLengths.size() ) + { + rowLengths[rowNumber] += 1; // One pressure column + } + } + } + } + + // For bubble DOFs: add 1 pressure column per bubble DOF row + for( localIndex kf = 0; kf < 2; ++kf ) + { + localIndex const faceIndex = elem2dToFaces[kfe][kf]; + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowNumber = bubbleDofNumber[faceIndex] + i - rankOffset; + if( rowNumber >= 0 && rowNumber < rowLengths.size() ) + { + rowLengths[rowNumber] += 1; // One pressure column + } + } + } + } ); + } ); +} + +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: +addPressureForceCouplingPattern( DomainPartition const & domain, + DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern ) const +{ + GEOS_MARK_FUNCTION; + + using namespace contact; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + string const flowDofKey = dofManager.getKey( m_pressureKey ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const & + bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + + globalIndex const rankOffset = dofManager.rankOffset(); + + string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); + SurfaceElementRegion const & fractureRegion = + elemManager.getRegion< SurfaceElementRegion >( fractureRegionName ); + FaceElementSubRegion const & fractureSubRegion = + fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); + + arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst(); + arrayView1d< globalIndex const > const & + flowDofNumber = fractureSubRegion.getReference< globalIndex_array >( flowDofKey ); + + // For each fracture element, add pattern for (displacement_row, pressure_col) and (bubble_row, pressure_col) + forAll< serialPolicy >( fractureSubRegion.size(), [=] ( localIndex const kfe ) + { + globalIndex const pressureColIndex = flowDofNumber[kfe]; + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elem2dToFaces[kfe][0] ); + + // For displacement DOFs + for( localIndex kf = 0; kf < 2; ++kf ) + { + localIndex const faceIndex = elem2dToFaces[kfe][kf]; + for( localIndex a = 0; a < numNodesPerFace; ++a ) + { + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowIndex = dispDofNumber[faceToNodeMap( faceIndex, a )] + i - rankOffset; + if( rowIndex >= 0 && rowIndex < pattern.numRows() ) + { + pattern.insertNonZero( rowIndex, pressureColIndex ); + } + } + } + } + + // For bubble DOFs + for( localIndex kf = 0; kf < 2; ++kf ) + { + localIndex const faceIndex = elem2dToFaces[kfe][kf]; + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowIndex = bubbleDofNumber[faceIndex] + i - rankOffset; + if( rowIndex >= 0 && rowIndex < pattern.numRows() ) + { + pattern.insertNonZero( rowIndex, pressureColIndex ); + } + } + } + } ); + } ); } template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: assembleForceResidualDerivativeWrtPressure( string const & meshName, MeshLevel const & mesh, - arrayView1d< string const > const & regionNames, + string_array const & regionNames, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( meshName, mesh, regionNames, dofManager, localMatrix, localRhs ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + using namespace contact; + + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + + string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const & bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + string const & flowDofKey = dofManager.getKey( m_pressureKey ); + + arrayView1d< globalIndex const > const dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + + string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); + + // Use the same kernel launch pattern as SolidMechanicsAugmentedLagrangianContact::assembleForceResidualPressureContribution + this->solidMechanicsSolver()->template forFiniteElementOnFractureSubRegions( meshName, + [&] ( string const &, + finiteElement::FiniteElementBase const & subRegionFE, + arrayView1d< localIndex const > const & faceElementList ) + { + // Get pressure DOF number from the fracture subregion + SurfaceElementRegion const & fractureRegion = mesh.getElemManager().getRegion< SurfaceElementRegion >( fractureRegionName ); + FaceElementSubRegion const & fractureSubRegion = fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); + arrayView1d< globalIndex const > const pressureDofNumber = fractureSubRegion.getReference< array1d< globalIndex > >( flowDofKey ); + + poromechanicsALMKernels::AssembleForceResidualDerivativeWrtPressureFactory + kernelFactory( dispDofNumber, + bubbleDofNumber, + dofManager.rankOffset(), + localMatrix, + localRhs, + 0.0, // dt not used + faceElementList, + pressureDofNumber ); + + // Note: const_cast is needed because interfaceBasedKernelApplication takes non-const mesh + // even though it only modifies the matrix/rhs which are passed separately + real64 maxResidual = finiteElement:: + interfaceBasedKernelApplication + < parallelDevicePolicy<>, + constitutive::NullModel >( const_cast< MeshLevel & >( mesh ), + fractureRegionName, + faceElementList, + subRegionFE, + "", + kernelFactory ); + + GEOS_UNUSED_VAR( maxResidual ); + } ); + + GEOS_UNUSED_VAR( regionNames ); } template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: -assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, - arrayView1d< string const > const & regionNames, +assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, + MeshLevel const & mesh, + string_array const & regionNames, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & GEOS_UNUSED_PARAM( localRhs ) ) + arrayView1d< real64 > const & localRhs ) +{ + GEOS_MARK_FUNCTION; + GEOS_UNUSED_VAR( regionNames ); + + using namespace contact; + + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + + CRSMatrixView< real64 const, localIndex const > const & + dFluxResidual_dNormalJump = getDerivativeFluxResidual_dNormalJump().toViewConst(); + + string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const & bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + string const & presDofKey = dofManager.getKey( m_pressureKey ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const & + bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + globalIndex const rankOffset = dofManager.rankOffset(); + + string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); + + // Maximum DOF sizes + constexpr localIndex maxNumNodesPerFace = m_maxFaceNodes; + constexpr localIndex maxNumUdofs = 2 * 3 * maxNumNodesPerFace; // 66 + constexpr localIndex numBdofs = 6; + + // Get the fracture subRegion + SurfaceElementRegion const & fractureRegion = elemManager.getRegion< SurfaceElementRegion >( fractureRegionName ); + FaceElementSubRegion const & subRegion = fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); + + localIndex const numElems = subRegion.size(); + + // Allocate temporary storage for aperture derivatives (computed via kernels) + array2d< real64 > dAperturedU( numElems, maxNumUdofs ); + array2d< real64 > dAperturedB( numElems, numBdofs ); + + // Initialize to zero and move to device for kernel access + dAperturedU.zero(); + dAperturedB.zero(); + dAperturedU.move( parallelDeviceMemorySpace, true ); + dAperturedB.move( parallelDeviceMemorySpace, true ); + + // Launch the ComputeApertureDerivatives kernel to fill dAperturedU and dAperturedB + // This is called for each element type (tri, quad, etc.) + this->solidMechanicsSolver()->template forFiniteElementOnFractureSubRegions( meshName, + [&] ( string const &, + finiteElement::FiniteElementBase const & subRegionFE, + arrayView1d< localIndex const > const & faceElementList ) + { + poromechanicsALMKernels::ComputeApertureDerivativesFactory + kernelFactory( dispDofNumber, + bubbleDofNumber, + dofManager.rankOffset(), + localMatrix, + localRhs, + 0.0, // dt not used + faceElementList, + dAperturedU.toView(), + dAperturedB.toView() ); + + real64 maxResidual = finiteElement:: + interfaceBasedKernelApplication + < parallelDevicePolicy<>, + constitutive::NullModel >( const_cast< MeshLevel & >( mesh ), + fractureRegionName, + faceElementList, + subRegionFE, + "", + kernelFactory ); + + GEOS_UNUSED_VAR( maxResidual ); + } ); + + // Move data to host for serial assembly + dAperturedU.move( hostMemorySpace ); + dAperturedB.move( hostMemorySpace ); + + // Now assemble using the pre-computed derivatives + string const & fluidName = subRegion.getReference< string >( FlowSolverBase::viewKeyStruct::fluidNamesString() ); + + SingleFluidBase const & fluid = this->template getConstitutiveModel< SingleFluidBase >( subRegion, fluidName ); + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const & density = fluid.density(); + + arrayView1d< globalIndex const > const & presDofNumber = subRegion.getReference< array1d< globalIndex > >( presDofKey ); + + arrayView2d< localIndex const > const & elemsToFaces = subRegion.faceList().toViewConst(); + + arrayView1d< integer const > const fractureState = subRegion.getField< contact::fractureState >(); + + // Get element area for proper scaling + // Note: dAperturedU/dB are computed as (1/area) * unitNormal^T * Atu/Atb + // For accumulation: dR_accum/du = density * unitNormal^T * Atu (no 1/area factor) + // For flux: dR_flux/du = dR/dAperture * (1/area) * unitNormal^T * Atu (with 1/area factor) + // So accumulation needs to multiply by area to cancel the 1/area in dAperturedU + arrayView1d< real64 const > const area = subRegion.getElementArea().toViewConst(); + + forAll< serialPolicy >( numElems, [&]( localIndex const kfe ) + { + localIndex const kf0 = elemsToFaces[kfe][0]; + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 ); + localIndex const numUdofs = 2 * 3 * numNodesPerFace; + + globalIndex nodeDOF[maxNumUdofs]; + globalIndex elemDOF[1]; + elemDOF[0] = presDofNumber[kfe]; + + stackArray1d< real64, maxNumUdofs > dRdU( maxNumUdofs ); + + bool const isFractureOpen = ( fractureState[kfe] == FractureState::Open ); + + // ==== Part 1: Apu - Nodal displacement contribution ==== + // Get DOF indices for displacement + for( localIndex kf = 0; kf < 2; ++kf ) + { + for( localIndex a = 0; a < numNodesPerFace; ++a ) + { + for( localIndex i = 0; i < 3; ++i ) + { + nodeDOF[kf * 3 * numNodesPerFace + 3 * a + i] = dispDofNumber[faceToNodeMap( elemsToFaces[kfe][kf], a )] + + LvArray::integerConversion< globalIndex >( i ); + } + } + } + + // Accumulation derivative w.r.t. nodal displacement + // dR_accum/du = density * unitNormal^T * Atu = density * area * dAperturedU + // (dAperturedU already has 1/area factor, so multiply by area to cancel it) + // Note: The ALM Atu matrix has opposite sign convention compared to LM's dAper/dU, + // so we need to negate the result to match the LM sign convention. + if( isFractureOpen ) + { + for( localIndex j = 0; j < numUdofs; ++j ) + { + dRdU( j ) = -density[kfe][0] * dAperturedU( kfe, j ) * area[kfe]; + } + + localIndex const localRow = LvArray::integerConversion< localIndex >( elemDOF[0] - rankOffset ); + + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( localRow, + nodeDOF, + dRdU.data(), + numUdofs ); + } + } + + // Flux derivative w.r.t. nodal displacement + localIndex const numColumns = dFluxResidual_dNormalJump.numNonZeros( kfe ); + arraySlice1d< localIndex const > const & columns = dFluxResidual_dNormalJump.getColumns( kfe ); + arraySlice1d< real64 const > const & values = dFluxResidual_dNormalJump.getEntries( kfe ); + + for( localIndex kfe1 = 0; kfe1 < numColumns; ++kfe1 ) + { + real64 const dR_dAper = values[kfe1]; + localIndex const kfe2 = columns[kfe1]; + + bool const isOpen = ( fractureState[kfe2] == FractureState::Open ); + if( !isOpen && !isFractureOpen ) + continue; + + localIndex const kf0_2 = elemsToFaces[kfe2][0]; + localIndex const numNodesPerFace2 = faceToNodeMap.sizeOfArray( kf0_2 ); + localIndex const numUdofs2 = 2 * 3 * numNodesPerFace2; + + // Get DOF indices for element kfe2 + globalIndex nodeDOF2[maxNumUdofs]; + for( localIndex kf = 0; kf < 2; ++kf ) + { + for( localIndex a = 0; a < numNodesPerFace2; ++a ) + { + for( localIndex i = 0; i < 3; ++i ) + { + nodeDOF2[kf * 3 * numNodesPerFace2 + 3 * a + i] = dispDofNumber[faceToNodeMap( elemsToFaces[kfe2][kf], a )] + + LvArray::integerConversion< globalIndex >( i ); + } + } + } + + // dR_flux/du = dR_flux/dAper * dAper/du (pre-computed for element kfe2) + // Note: Negate to match LM sign convention (ALM Atu has opposite sign) + stackArray1d< real64, maxNumUdofs > dRdU2( maxNumUdofs ); + for( localIndex j = 0; j < numUdofs2; ++j ) + { + dRdU2( j ) = -dR_dAper * dAperturedU( kfe2, j ); + } + + localIndex const localRow = LvArray::integerConversion< localIndex >( elemDOF[0] - rankOffset ); + + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( localRow, + nodeDOF2, + dRdU2.data(), + numUdofs2 ); + } + } + + // ==== Part 2: Apb - Bubble displacement contribution ==== + globalIndex bubbleDOF[numBdofs]; + + // Get DOF indices for bubble + for( localIndex kf = 0; kf < 2; ++kf ) + { + localIndex const faceIndex = elemsToFaces[kfe][kf]; + for( localIndex i = 0; i < 3; ++i ) + { + bubbleDOF[kf * 3 + i] = bubbleDofNumber[faceIndex] + LvArray::integerConversion< globalIndex >( i ); + } + } + + // Accumulation derivative w.r.t. bubble displacement + // dR_accum/db = density * unitNormal^T * Atb = density * area * dAperturedB + // (dAperturedB already has 1/area factor, so multiply by area to cancel it) + // Note: Negate to match LM sign convention (ALM Atb has opposite sign) + if( isFractureOpen ) + { + real64 dRdB[numBdofs]; + for( localIndex j = 0; j < numBdofs; ++j ) + { + dRdB[j] = -density[kfe][0] * dAperturedB( kfe, j ) * area[kfe]; + } + + localIndex const localRow = LvArray::integerConversion< localIndex >( elemDOF[0] - rankOffset ); + + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( localRow, + bubbleDOF, + dRdB, + numBdofs ); + } + } + + // Flux derivative w.r.t. bubble displacement + for( localIndex kfe1 = 0; kfe1 < numColumns; ++kfe1 ) + { + real64 const dR_dAper = values[kfe1]; + localIndex const kfe2 = columns[kfe1]; + + bool const isOpen = ( fractureState[kfe2] == FractureState::Open ); + if( !isOpen && !isFractureOpen ) + continue; + + // Get DOF indices for bubble of element kfe2 + globalIndex bubbleDOF2[numBdofs]; + for( localIndex kf = 0; kf < 2; ++kf ) + { + localIndex const faceIndex = elemsToFaces[kfe2][kf]; + for( localIndex i = 0; i < 3; ++i ) + { + bubbleDOF2[kf * 3 + i] = bubbleDofNumber[faceIndex] + LvArray::integerConversion< globalIndex >( i ); + } + } + + // dR_flux/db = dR_flux/dAper * dAper/db (pre-computed for element kfe2) + // Note: Negate to match LM sign convention (ALM Atb has opposite sign) + real64 dRdB2[numBdofs]; + for( localIndex j = 0; j < numBdofs; ++j ) + { + dRdB2[j] = -dR_dAper * dAperturedB( kfe2, j ); + } + + localIndex const localRow = LvArray::integerConversion< localIndex >( elemDOF[0] - rankOffset ); + + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( localRow, + bubbleDOF2, + dRdB2, + numBdofs ); + } + } + } ); +} + +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: +addMatrixPressureBubbleCouplingNNZ( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< localIndex > const & rowLengths ) const +{ + GEOS_MARK_FUNCTION; + + using namespace contact; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & regionNames ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + // Loop over matrix cell regions that have bubbles + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, + [&]( localIndex const, CellElementSubRegion const & subRegion ) + { + if( !subRegion.hasWrapper( "bubbleElementsList" ) ) + return; + + arrayView1d< localIndex const > const bubbleElems = subRegion.bubbleElementsList(); + arrayView2d< localIndex const > const elemsToFaces = subRegion.faceElementsList(); + + forAll< serialPolicy >( bubbleElems.size(), [=, &rowLengths]( localIndex const kk ) + { + localIndex const faceIndex = elemsToFaces[kk][0]; + + // Add 1 pressure column for each of the 3 bubble DOFs + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowNumber = bubbleDofNumber[faceIndex] + i - rankOffset; + if( rowNumber >= 0 && rowNumber < rowLengths.size() ) + { + rowLengths[rowNumber] += 1; // One pressure DOF from matrix cell + } + } + } ); + } ); + } ); +} + +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: +addMatrixPressureBubbleCouplingPattern( DomainPartition const & domain, + DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern ) const { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( mesh, regionNames, dofManager, localMatrix ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + using namespace contact; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & regionNames ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + string const flowDofKey = dofManager.getKey( m_pressureKey ); + + arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + // Loop over matrix cell regions that have bubbles + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, + [&]( localIndex const, CellElementSubRegion const & subRegion ) + { + if( !subRegion.hasWrapper( "bubbleElementsList" ) ) + return; + + arrayView1d< localIndex const > const bubbleElems = subRegion.bubbleElementsList(); + arrayView2d< localIndex const > const elemsToFaces = subRegion.faceElementsList(); + arrayView1d< globalIndex const > const pressureDofNumber = subRegion.getReference< array1d< globalIndex > >( flowDofKey ); + + forAll< serialPolicy >( bubbleElems.size(), [=]( localIndex const kk ) + { + localIndex const k = bubbleElems[kk]; + localIndex const faceIndex = elemsToFaces[kk][0]; + globalIndex const pressureColIndex = pressureDofNumber[k]; + + // Add pattern for (bubble_row, pressure_col) + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowIndex = bubbleDofNumber[faceIndex] + i - rankOffset; + if( rowIndex >= 0 && rowIndex < pattern.numRows() ) + { + pattern.insertNonZero( rowIndex, pressureColIndex ); + } + } + } ); + } ); + } ); +} + +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: +assembleMatrixPressureBubbleContribution( real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + GEOS_MARK_FUNCTION; + + using namespace contact; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + NodeManager const & nodeManager = mesh.getNodeManager(); + FaceManager const & faceManager = mesh.getFaceManager(); + ElementRegionManager & elemManager = mesh.getElemManager(); + + string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const & bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + string const & flowDofKey = dofManager.getKey( m_pressureKey ); + + arrayView1d< globalIndex const > const dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + + // Identify poromechanics regions (cells with porous solid) + set< string > poromechanicsRegions; + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, + [&]( localIndex const regionIndex, CellElementSubRegion const & subRegion ) + { + if( subRegion.hasWrapper( FlowSolverBase::viewKeyStruct::solidNamesString() ) ) + { + poromechanicsRegions.insert( regionNames[regionIndex] ); + } + } ); + + string_array poromechanicsRegionNames; + poromechanicsRegionNames.reserve( poromechanicsRegions.size() ); + for( auto const & region : poromechanicsRegions ) + { + poromechanicsRegionNames.emplace_back( region ); + } + + // Launch the kernel on matrix cells with bubbles + poromechanicsMatrixBubbleKernels::MatrixPressureBubbleFactory kernelFactory( dispDofNumber, + bubbleDofNumber, + dofManager.rankOffset(), + localMatrix, + localRhs, + dt, + flowDofKey ); + + real64 maxResidual = finiteElement::regionBasedKernelApplication + < parallelDevicePolicy<>, + constitutive::PorousSolid< constitutive::ElasticIsotropic, constitutive::ConstantPermeability >, + CellElementSubRegion >( mesh, + poromechanicsRegionNames, + this->solidMechanicsSolver()->getDiscretizationName(), + FlowSolverBase::viewKeyStruct::solidNamesString(), + kernelFactory ); + GEOS_UNUSED_VAR( maxResidual ); + } ); } diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp index 15f24d838d9..9a427168a21 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -141,13 +141,14 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec void assembleForceResidualDerivativeWrtPressure( string const & meshName, MeshLevel const & mesh, - arrayView1d< string const > const & regionNames, + string_array const & regionNames, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ); - void assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, - arrayView1d< string const > const & regionNames, + void assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, + MeshLevel const & mesh, + string_array const & regionNames, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ); @@ -172,6 +173,62 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec DofManager const & dofManager, SparsityPatternView< globalIndex > const & pattern ) const; + /** + * @Brief add the nnz induced by the pressure-force coupling (Aup, Abp) + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param rowLenghts the nnz in each row + */ + void addPressureForceCouplingNNZ( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< localIndex > const & rowLengths ) const; + + /** + * @Brief add the sparsity pattern induced by the pressure-force coupling (Aup, Abp) + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param pattern the sparsity pattern + */ + void addPressureForceCouplingPattern( DomainPartition const & domain, + DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern ) const; + + /** + * @Brief add the nnz induced by the matrix pressure-bubble coupling (Abp_matrix) + * This handles the contribution of matrix cell pressure on bubble DOFs. + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param rowLengths the nnz in each row + */ + void addMatrixPressureBubbleCouplingNNZ( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< localIndex > const & rowLengths ) const; + + /** + * @Brief add the sparsity pattern induced by the matrix pressure-bubble coupling + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param pattern the sparsity pattern + */ + void addMatrixPressureBubbleCouplingPattern( DomainPartition const & domain, + DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern ) const; + + /** + * @Brief assemble the contribution of matrix cell pressure on bubble DOFs + * with full Jacobian for fully-implicit coupling. + * @param dt the time step size + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param localMatrix the local system matrix + * @param localRhs the local system right-hand side vector + */ + void assembleMatrixPressureBubbleContribution( real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + /** * @brief Set up the Dflux_dApertureMatrix object * diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp new file mode 100644 index 00000000000..f069b7f7b90 --- /dev/null +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp @@ -0,0 +1,1079 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file SinglePhasePoromechanicsConformingFracturesALMKernels.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALMKERNELS_HPP_ +#define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALMKERNELS_HPP_ + +#include "finiteElement/kernelInterface/ImplicitKernelBase.hpp" +#include "physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsConformingContactKernelsBase.hpp" +#include "physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsConformingContactKernelsHelper.hpp" +#include "physicsSolvers/solidMechanics/contact/ContactFields.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" + +namespace geos +{ + +namespace poromechanicsALMKernels +{ + +/** + * @brief Kernel to assemble the pressure force contribution to displacement equations + * AND the Jacobian derivative dR_u/dP for fully implicit coupling. + * + * This kernel is similar to SolidMechanicsConformingPressureContributionKernels::AssemblePressureContribution + * but also adds the Jacobian terms needed for fully implicit poromechanics coupling. + * The sign convention matches the existing sequential coupling kernel. + */ +template< typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +class AssembleForceResidualDerivativeWrtPressure : + public solidMechanicsConformingContactKernels::ConformingContactKernelsBase< CONSTITUTIVE_TYPE, + FE_TYPE > +{ + +public: + /// Alias for the base class + using Base = solidMechanicsConformingContactKernels::ConformingContactKernelsBase< CONSTITUTIVE_TYPE, + FE_TYPE >; + + /// Number of nodes per element + static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem; + + /// The number of displacement dofs per element. + static constexpr int numUdofs = Base::numUdofs; + + /// The number of bubble dofs per element. + static constexpr int numBdofs = Base::numBdofs; + + using Base::m_elemsToFaces; + using Base::m_faceToNodes; + using Base::m_X; + using Base::m_finiteElementSpace; + using Base::m_dofNumber; + using Base::m_bDofNumber; + using Base::m_dofRankOffset; + using Base::m_rotationMatrix; + using Base::m_matrix; + using Base::m_rhs; + + /** + * @brief Constructor + */ + AssembleForceResidualDerivativeWrtPressure( NodeManager const & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + FaceElementSubRegion & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + arrayView1d< globalIndex const > const uDofNumber, + arrayView1d< globalIndex const > const bDofNumber, + globalIndex const rankOffset, + CRSMatrixView< real64, globalIndex const > const inputMatrix, + arrayView1d< real64 > const inputRhs, + real64 const inputDt, + arrayView1d< localIndex const > const & faceElementList, + arrayView1d< globalIndex const > const & pressureDofNumber ): + Base( nodeManager, + edgeManager, + faceManager, + targetRegionIndex, + elementSubRegion, + finiteElementSpace, + inputConstitutiveType, + uDofNumber, + bDofNumber, + rankOffset, + inputMatrix, + inputRhs, + inputDt, + faceElementList ), + m_pressure( elementSubRegion.getField< fields::flow::pressure >().toView() ), + m_pressureDofNumber( pressureDofNumber ) + {} + + /** + * @copydoc finiteElement::InterfaceKernelBase::StackVariables + */ + struct StackVariables : public Base::StackVariables + { +public: + + GEOS_HOST_DEVICE + StackVariables(): + Base::StackVariables(), + dispEqnRowIndices{}, + bEqnRowIndices{}, + pressureColIndex( 0 ), + unitNormal{}, + localRu{}, + localRb{}, + dRudP{}, + dRbdP{} + {} + + /// Row indices for displacement equations + globalIndex dispEqnRowIndices[numUdofs]; + + /// Row indices for bubble equations + globalIndex bEqnRowIndices[numBdofs]; + + /// Column index for pressure + globalIndex pressureColIndex; + + /// Unit normal from rotation matrix + real64 unitNormal[3]; + + /// Local residual for displacement + real64 localRu[numUdofs]; + + /// Local residual for bubble + real64 localRb[numBdofs]; + + /// Jacobian dRu/dP + real64 dRudP[numUdofs]; + + /// Jacobian dRb/dP + real64 dRbdP[numBdofs]; + }; + + template< typename POLICY, + typename KERNEL_TYPE > + static + real64 + kernelLaunch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ) + { + return Base::template kernelLaunch< POLICY, KERNEL_TYPE >( numElems, kernelComponent ); + } + + /** + * @brief Setup for each element + */ + GEOS_HOST_DEVICE + inline + void setup( localIndex const k, + StackVariables & stack ) const + { + constexpr int shift = numNodesPerElem * 3; + + int permutation[numNodesPerElem]; + m_finiteElementSpace.getPermutation( permutation ); + + localIndex const kf0 = m_elemsToFaces[k][0]; + localIndex const kf1 = m_elemsToFaces[k][1]; + + for( localIndex a = 0; a < numNodesPerElem; ++a ) + { + localIndex const kn0 = m_faceToNodes( kf0, a ); + localIndex const kn1 = m_faceToNodes( kf1, a ); + for( int i = 0; i < 3; ++i ) + { + stack.dispEqnRowIndices[a * 3 + i] = m_dofNumber[kn0] + i - m_dofRankOffset; + stack.dispEqnRowIndices[shift + a * 3 + i] = m_dofNumber[kn1] + i - m_dofRankOffset; + stack.X[a][i] = m_X[m_faceToNodes( kf0, permutation[a] )][i]; + } + } + + for( int i = 0; i < 3; ++i ) + { + stack.unitNormal[i] = m_rotationMatrix( k, i, 0 ); + } + + for( int i = 0; i < 3; ++i ) + { + stack.bEqnRowIndices[i] = m_bDofNumber[kf0] + i - m_dofRankOffset; + stack.bEqnRowIndices[3 + i] = m_bDofNumber[kf1] + i - m_dofRankOffset; + } + + stack.pressureColIndex = m_pressureDofNumber[k]; + } + + /** + * @brief Complete - compute residuals and Jacobians + */ + GEOS_HOST_DEVICE + inline + real64 complete( localIndex const k, + StackVariables & stack ) const + { + // Compute residual and Jacobian for displacement DOFs: + // R_u = transpose(Atu) * unitNormal * pressure + // dR_u/dP = transpose(Atu) * unitNormal + // Note: This matches the existing AssemblePressureContribution kernel used for sequential coupling. + LvArray::tensorOps::Ri_eq_AjiBj< numUdofs, 3 >( stack.localRu, stack.localAtu, stack.unitNormal ); + LvArray::tensorOps::Ri_eq_AjiBj< numUdofs, 3 >( stack.dRudP, stack.localAtu, stack.unitNormal ); + + // Scale residual by pressure + LvArray::tensorOps::scale< numUdofs >( stack.localRu, m_pressure[k] ); + + // Compute residual and Jacobian for bubble DOFs: + // R_b = transpose(Atb) * unitNormal * pressure + // dR_b/dP = transpose(Atb) * unitNormal + LvArray::tensorOps::Ri_eq_AjiBj< numBdofs, 3 >( stack.localRb, stack.localAtb, stack.unitNormal ); + LvArray::tensorOps::Ri_eq_AjiBj< numBdofs, 3 >( stack.dRbdP, stack.localAtb, stack.unitNormal ); + + // Scale residual by pressure + LvArray::tensorOps::scale< numBdofs >( stack.localRb, m_pressure[k] ); + + // Assemble displacement contributions + for( localIndex i = 0; i < numUdofs; ++i ) + { + localIndex const dof = LvArray::integerConversion< localIndex >( stack.dispEqnRowIndices[i] ); + + if( dof < 0 || dof >= m_matrix.numRows() ) + continue; + + // Add residual + RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[dof], stack.localRu[i] ); + + // Add Jacobian + m_matrix.template addToRow< parallelDeviceAtomic >( dof, + &stack.pressureColIndex, + &stack.dRudP[i], + 1 ); + } + + // Assemble bubble contributions + for( localIndex i = 0; i < numBdofs; ++i ) + { + localIndex const dof = LvArray::integerConversion< localIndex >( stack.bEqnRowIndices[i] ); + + if( dof < 0 || dof >= m_matrix.numRows() ) + continue; + + // Add residual + RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[dof], stack.localRb[i] ); + + // Add Jacobian + m_matrix.template addToRow< parallelDeviceAtomic >( dof, + &stack.pressureColIndex, + &stack.dRbdP[i], + 1 ); + } + + return 0.0; + } + +protected: + + /// Pressure array + arrayView1d< real64 const > const m_pressure; + + /// Pressure DOF numbers + arrayView1d< globalIndex const > const m_pressureDofNumber; +}; + +/// Factory for AssembleForceResidualDerivativeWrtPressure kernel +using AssembleForceResidualDerivativeWrtPressureFactory = + finiteElement::InterfaceKernelFactory< AssembleForceResidualDerivativeWrtPressure, + arrayView1d< globalIndex const > const, + arrayView1d< globalIndex const > const, + globalIndex const, + CRSMatrixView< real64, globalIndex const > const, + arrayView1d< real64 > const, + real64 const, + arrayView1d< localIndex const > const, + arrayView1d< globalIndex const > const >; + + +/** + * @brief Kernel to assemble the Jacobian derivative dR_p/du and dR_p/db for fully implicit coupling. + * + * This kernel computes the derivative of the flow accumulation residual with respect to + * nodal displacement and bubble displacement, using the properly computed Atu/Atb matrices. + */ +template< typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +class AssembleFluidMassResidualDerivativeWrtDisplacement : + public solidMechanicsConformingContactKernels::ConformingContactKernelsBase< CONSTITUTIVE_TYPE, + FE_TYPE > +{ + +public: + /// Alias for the base class + using Base = solidMechanicsConformingContactKernels::ConformingContactKernelsBase< CONSTITUTIVE_TYPE, + FE_TYPE >; + + /// Number of nodes per element + static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem; + + /// The number of displacement dofs per element. + static constexpr int numUdofs = Base::numUdofs; + + /// The number of bubble dofs per element. + static constexpr int numBdofs = Base::numBdofs; + + using Base::m_elemsToFaces; + using Base::m_faceToNodes; + using Base::m_X; + using Base::m_finiteElementSpace; + using Base::m_dofNumber; + using Base::m_bDofNumber; + using Base::m_dofRankOffset; + using Base::m_rotationMatrix; + using Base::m_matrix; + using Base::m_rhs; + + /** + * @brief Constructor + */ + AssembleFluidMassResidualDerivativeWrtDisplacement( NodeManager const & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + FaceElementSubRegion & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + arrayView1d< globalIndex const > const uDofNumber, + arrayView1d< globalIndex const > const bDofNumber, + globalIndex const rankOffset, + CRSMatrixView< real64, globalIndex const > const inputMatrix, + arrayView1d< real64 > const inputRhs, + real64 const inputDt, + arrayView1d< localIndex const > const & faceElementList, + arrayView1d< globalIndex const > const & pressureDofNumber, + arrayView1d< real64 const > const & density, + arrayView1d< integer const > const & fractureState ): + Base( nodeManager, + edgeManager, + faceManager, + targetRegionIndex, + elementSubRegion, + finiteElementSpace, + inputConstitutiveType, + uDofNumber, + bDofNumber, + rankOffset, + inputMatrix, + inputRhs, + inputDt, + faceElementList ), + m_pressureDofNumber( pressureDofNumber ), + m_density( density ), + m_fractureState( fractureState ) + {} + + /** + * @copydoc finiteElement::InterfaceKernelBase::StackVariables + */ + struct StackVariables : public Base::StackVariables + { +public: + + GEOS_HOST_DEVICE + StackVariables(): + Base::StackVariables(), + dispColIndices{}, + bColIndices{}, + pressureRowIndex( 0 ), + unitNormal{}, + dRpdU{}, + dRpdB{} + {} + + /// Column indices for displacement + globalIndex dispColIndices[numUdofs]; + + /// Column indices for bubble + globalIndex bColIndices[numBdofs]; + + /// Row index for pressure equation + globalIndex pressureRowIndex; + + /// Unit normal from rotation matrix + real64 unitNormal[3]; + + /// Jacobian dRp/dU + real64 dRpdU[numUdofs]; + + /// Jacobian dRp/dB + real64 dRpdB[numBdofs]; + }; + + template< typename POLICY, + typename KERNEL_TYPE > + static + real64 + kernelLaunch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ) + { + return Base::template kernelLaunch< POLICY, KERNEL_TYPE >( numElems, kernelComponent ); + } + + /** + * @brief Setup for each element + */ + GEOS_HOST_DEVICE + inline + void setup( localIndex const k, + StackVariables & stack ) const + { + constexpr int shift = numNodesPerElem * 3; + + int permutation[numNodesPerElem]; + m_finiteElementSpace.getPermutation( permutation ); + + localIndex const kf0 = m_elemsToFaces[k][0]; + localIndex const kf1 = m_elemsToFaces[k][1]; + + for( localIndex a = 0; a < numNodesPerElem; ++a ) + { + localIndex const kn0 = m_faceToNodes( kf0, a ); + localIndex const kn1 = m_faceToNodes( kf1, a ); + for( int i = 0; i < 3; ++i ) + { + stack.dispColIndices[a * 3 + i] = m_dofNumber[kn0] + i; + stack.dispColIndices[shift + a * 3 + i] = m_dofNumber[kn1] + i; + stack.X[a][i] = m_X[m_faceToNodes( kf0, permutation[a] )][i]; + } + } + + for( int i = 0; i < 3; ++i ) + { + stack.unitNormal[i] = m_rotationMatrix( k, i, 0 ); + } + + for( int i = 0; i < 3; ++i ) + { + stack.bColIndices[i] = m_bDofNumber[kf0] + i; + stack.bColIndices[3 + i] = m_bDofNumber[kf1] + i; + } + + stack.pressureRowIndex = m_pressureDofNumber[k] - m_dofRankOffset; + } + + /** + * @brief Complete - compute Jacobians dRp/dU and dRp/dB + */ + GEOS_HOST_DEVICE + inline + real64 complete( localIndex const k, + StackVariables & stack ) const + { + // Only assemble for open fractures + if( m_fractureState[k] != fields::contact::FractureState::Open ) + { + return 0.0; + } + + // Compute Jacobian for pressure equation w.r.t. displacement: + // dR_p/du = density * unitNormal^T * Atu + // This is: dRpdU[j] = sum_i( unitNormal[i] * Atu[i][j] ) * density + // Use Ri_eq_AjiBj for transpose: R[i] = sum_j(A[j][i] * B[j]) + // Note: The ALM Atu matrix has opposite sign convention compared to LM's dAper/dU, + // so we negate the result to match the LM sign convention. + LvArray::tensorOps::Ri_eq_AjiBj< numUdofs, 3 >( stack.dRpdU, stack.localAtu, stack.unitNormal ); + LvArray::tensorOps::scale< numUdofs >( stack.dRpdU, -m_density[k] ); + + // Compute Jacobian for pressure equation w.r.t. bubble displacement: + // dR_p/db = density * unitNormal^T * Atb + // Note: Same sign convention fix as for displacement. + LvArray::tensorOps::Ri_eq_AjiBj< numBdofs, 3 >( stack.dRpdB, stack.localAtb, stack.unitNormal ); + LvArray::tensorOps::scale< numBdofs >( stack.dRpdB, -m_density[k] ); + + // Assemble into matrix + localIndex const localRow = LvArray::integerConversion< localIndex >( stack.pressureRowIndex ); + + if( localRow >= 0 && localRow < m_matrix.numRows() ) + { + // Add dRp/dU + m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( localRow, + stack.dispColIndices, + stack.dRpdU, + numUdofs ); + + // Add dRp/dB + m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( localRow, + stack.bColIndices, + stack.dRpdB, + numBdofs ); + } + + return 0.0; + } + +protected: + + /// Pressure DOF numbers + arrayView1d< globalIndex const > const m_pressureDofNumber; + + /// Fluid density + arrayView1d< real64 const > const m_density; + + /// Fracture state + arrayView1d< integer const > const m_fractureState; +}; + +/// Factory for AssembleFluidMassResidualDerivativeWrtDisplacement kernel +using AssembleFluidMassResidualDerivativeWrtDisplacementFactory = + finiteElement::InterfaceKernelFactory< AssembleFluidMassResidualDerivativeWrtDisplacement, + arrayView1d< globalIndex const > const, + arrayView1d< globalIndex const > const, + globalIndex const, + CRSMatrixView< real64, globalIndex const > const, + arrayView1d< real64 > const, + real64 const, + arrayView1d< localIndex const > const, + arrayView1d< globalIndex const > const, + arrayView1d< real64 const > const, + arrayView1d< integer const > const >; + + +/** + * @brief Kernel to compute d(aperture)/du and d(aperture)/db for each fracture element. + * + * This kernel computes the derivative of aperture with respect to displacement and bubble DOFs + * using the properly computed Atu/Atb matrices. The results are stored for later use in + * flux derivative computation. + * + * d(aperture)/du = (1/area) * unitNormal^T * Atu + * d(aperture)/db = (1/area) * unitNormal^T * Atb + * + * IMPORTANT: The ALM Atu/Atb matrices have opposite sign convention compared to the LM + * dAper/dU formulas. When using these pre-computed values, they must be negated to match + * the LM sign convention. The assembly in assembleFluidMassResidualDerivativeWrtDisplacement + * handles this negation. + */ +template< typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +class ComputeApertureDerivatives : + public solidMechanicsConformingContactKernels::ConformingContactKernelsBase< CONSTITUTIVE_TYPE, + FE_TYPE > +{ + +public: + /// Alias for the base class + using Base = solidMechanicsConformingContactKernels::ConformingContactKernelsBase< CONSTITUTIVE_TYPE, + FE_TYPE >; + + /// Number of nodes per element + static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem; + + /// The number of displacement dofs per element. + static constexpr int numUdofs = Base::numUdofs; + + /// The number of bubble dofs per element. + static constexpr int numBdofs = Base::numBdofs; + + using Base::m_elemsToFaces; + using Base::m_faceToNodes; + using Base::m_X; + using Base::m_finiteElementSpace; + using Base::m_dofNumber; + using Base::m_bDofNumber; + using Base::m_dofRankOffset; + using Base::m_rotationMatrix; + using Base::m_matrix; + using Base::m_rhs; + + /** + * @brief Constructor + */ + ComputeApertureDerivatives( NodeManager const & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + FaceElementSubRegion & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + arrayView1d< globalIndex const > const uDofNumber, + arrayView1d< globalIndex const > const bDofNumber, + globalIndex const rankOffset, + CRSMatrixView< real64, globalIndex const > const inputMatrix, + arrayView1d< real64 > const inputRhs, + real64 const inputDt, + arrayView1d< localIndex const > const & faceElementList, + arrayView2d< real64 > const & dAperturedU, + arrayView2d< real64 > const & dAperturedB ): + Base( nodeManager, + edgeManager, + faceManager, + targetRegionIndex, + elementSubRegion, + finiteElementSpace, + inputConstitutiveType, + uDofNumber, + bDofNumber, + rankOffset, + inputMatrix, + inputRhs, + inputDt, + faceElementList ), + m_area( elementSubRegion.getElementArea().toViewConst() ), + m_dAperturedU( dAperturedU ), + m_dAperturedB( dAperturedB ) + {} + + /** + * @copydoc finiteElement::InterfaceKernelBase::StackVariables + */ + struct StackVariables : public Base::StackVariables + { +public: + + GEOS_HOST_DEVICE + StackVariables(): + Base::StackVariables(), + unitNormal{}, + localDAperturedU{}, + localDAperturedB{} + {} + + /// Unit normal from rotation matrix + real64 unitNormal[3]; + + /// Local d(aperture)/dU + real64 localDAperturedU[numUdofs]; + + /// Local d(aperture)/dB + real64 localDAperturedB[numBdofs]; + }; + + template< typename POLICY, + typename KERNEL_TYPE > + static + real64 + kernelLaunch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ) + { + return Base::template kernelLaunch< POLICY, KERNEL_TYPE >( numElems, kernelComponent ); + } + + /** + * @brief Setup for each element + */ + GEOS_HOST_DEVICE + inline + void setup( localIndex const k, + StackVariables & stack ) const + { + int permutation[numNodesPerElem]; + m_finiteElementSpace.getPermutation( permutation ); + + localIndex const kf0 = m_elemsToFaces[k][0]; + + for( localIndex a = 0; a < numNodesPerElem; ++a ) + { + for( int i = 0; i < 3; ++i ) + { + stack.X[a][i] = m_X[m_faceToNodes( kf0, permutation[a] )][i]; + } + } + + for( int i = 0; i < 3; ++i ) + { + stack.unitNormal[i] = m_rotationMatrix( k, i, 0 ); + } + } + + /** + * @brief Complete - compute and store d(aperture)/dU and d(aperture)/dB + */ + GEOS_HOST_DEVICE + inline + real64 complete( localIndex const k, + StackVariables & stack ) const + { + real64 const invArea = 1.0 / m_area[k]; + + // Compute d(aperture)/dU = (1/area) * unitNormal^T * Atu + // This gives: dAperturedU[j] = (1/area) * sum_i(unitNormal[i] * Atu[i][j]) + // Use Ri_eq_AjiBj for transpose: R[i] = sum_j(A[j][i] * B[j]) + LvArray::tensorOps::Ri_eq_AjiBj< numUdofs, 3 >( stack.localDAperturedU, stack.localAtu, stack.unitNormal ); + LvArray::tensorOps::scale< numUdofs >( stack.localDAperturedU, invArea ); + + // Compute d(aperture)/dB = (1/area) * unitNormal^T * Atb + LvArray::tensorOps::Ri_eq_AjiBj< numBdofs, 3 >( stack.localDAperturedB, stack.localAtb, stack.unitNormal ); + LvArray::tensorOps::scale< numBdofs >( stack.localDAperturedB, invArea ); + + // Store the results + for( localIndex i = 0; i < numUdofs; ++i ) + { + m_dAperturedU( k, i ) = stack.localDAperturedU[i]; + } + for( localIndex i = 0; i < numBdofs; ++i ) + { + m_dAperturedB( k, i ) = stack.localDAperturedB[i]; + } + + return 0.0; + } + +protected: + + /// Element area + arrayView1d< real64 const > const m_area; + + /// Storage for d(aperture)/dU + arrayView2d< real64 > const m_dAperturedU; + + /// Storage for d(aperture)/dB + arrayView2d< real64 > const m_dAperturedB; +}; + +/// Factory for ComputeApertureDerivatives kernel +using ComputeApertureDerivativesFactory = + finiteElement::InterfaceKernelFactory< ComputeApertureDerivatives, + arrayView1d< globalIndex const > const, + arrayView1d< globalIndex const > const, + globalIndex const, + CRSMatrixView< real64, globalIndex const > const, + arrayView1d< real64 > const, + real64 const, + arrayView1d< localIndex const > const, + arrayView2d< real64 > const, + arrayView2d< real64 > const >; + +} // namespace poromechanicsALMKernels + + +namespace poromechanicsMatrixBubbleKernels +{ + +/** + * @brief Kernel to compute the contribution of matrix cell pressure on bubble DOFs + * with full Jacobian for fully-implicit poromechanics coupling. + * + * This kernel is based on SolidMechanicsPressureFaceBubbleKernels but adds the Jacobian + * term dR_b/dP_matrix needed for fully implicit coupling. + * + * Physics: + * R_b = -∫ B_b^T * (biot * p * I) * detJ dΩ + * dR_b/dP = -∫ B_b^T * (biot * I) * detJ dΩ + * + * where B_b is the strain-displacement matrix for bubble functions. + */ +template< typename SUBREGION_TYPE, + typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +class MatrixPressureBubbleKernels : + public finiteElement::ImplicitKernelBase< SUBREGION_TYPE, + CONSTITUTIVE_TYPE, + FE_TYPE, + 3, + 3 > +{ +public: + + /// Alias for the base class + using Base = finiteElement::ImplicitKernelBase< SUBREGION_TYPE, + CONSTITUTIVE_TYPE, + FE_TYPE, + 3, + 3 >; + + /// Number of nodes per element + static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem; + + /// Compile time value for the number of faces per element. + static constexpr int numFacesPerElem = FE_TYPE::numFaces; + + /// Compile time value for the number of quadrature points per element. + static constexpr int numQuadraturePointsPerElem = FE_TYPE::numQuadraturePoints; + + using Base::m_elemsToNodes; + using Base::m_finiteElementSpace; + using Base::m_constitutiveUpdate; + using Base::m_dofNumber; + using Base::m_dofRankOffset; + using Base::m_matrix; + using Base::m_rhs; + using Base::m_dt; + + /** + * @brief Constructor + */ + MatrixPressureBubbleKernels( NodeManager const & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + SUBREGION_TYPE const & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + arrayView1d< globalIndex const > const uDofNumber, + arrayView1d< globalIndex const > const bDofNumber, + globalIndex const rankOffset, + CRSMatrixView< real64, globalIndex const > const inputMatrix, + arrayView1d< real64 > const inputRhs, + real64 const inputDt, + string const & pressureDofKey ): + Base( nodeManager, + edgeManager, + faceManager, + targetRegionIndex, + elementSubRegion, + finiteElementSpace, + inputConstitutiveType, + uDofNumber, + rankOffset, + inputMatrix, + inputRhs, + inputDt ), + m_X( nodeManager.referencePosition()), + m_bDofNumber( bDofNumber ), + m_pDofNumber( elementSubRegion.template getReference< array1d< globalIndex > >( pressureDofKey ) ), + m_bubbleElems( elementSubRegion.bubbleElementsList() ), + m_elemsToFaces( elementSubRegion.faceElementsList() ), + m_pressure( elementSubRegion.template getField< fields::flow::pressure >().toViewConst() ) + {} + + //*************************************************************************** + + /** + * @copydoc finiteElement::ImplicitKernelBase::StackVariables + */ + struct StackVariables + { +public: + + /// The number of displacement dofs per element. + static constexpr int numUdofs = numNodesPerElem * 3; + + /// The number of bubble dofs per element (all faces). + static constexpr int numBubbleUdofs = numFacesPerElem * 3; + + /// The number of pressure dofs per element. + static constexpr int numPdofs = 1; + + /** + * Default constructor + */ + GEOS_HOST_DEVICE + StackVariables(): + bEqnRowIndices{}, + pColIndex( 0 ), + localRb{}, + localdRbdP{}, + X{ {} }, + pLocal{} + {} + + /// C-array storage for the element local row degrees of freedom (bubble). + globalIndex bEqnRowIndices[3]; + + /// Column index for pressure DOF + globalIndex pColIndex; + + /// C-array storage for the element local Rb residual vector. + real64 localRb[numBubbleUdofs]; + + /// C-array storage for the element local dRb/dP Jacobian. + real64 localdRbdP[numBubbleUdofs]; + + /// local nodal coordinates + real64 X[ numNodesPerElem ][ 3 ]; + + /// local pressure + real64 pLocal[numPdofs]; + + }; + + //*************************************************************************** + + /** + * @copydoc ::geos::finiteElement::KernelBase::kernelLaunch + * + * @detail it uses the kernelLaunch interface of KernelBase but it only launches the kernel + * on the set of elements that have bubble dof within the subregion. + */ + template< typename POLICY, + typename KERNEL_TYPE > + static + real64 + kernelLaunch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ) + { + GEOS_MARK_FUNCTION; + GEOS_UNUSED_VAR( numElems ); + + // Define a RAJA reduction variable to get the maximum residual contribution. + RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxResidual( 0 ); + + forAll< POLICY >( kernelComponent.m_bubbleElems.size(), + [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + typename KERNEL_TYPE::StackVariables stack; + + kernelComponent.setup( i, stack ); + for( integer q=0; q( dBubbleNdX, 0 ); + + real64 const detJ = m_finiteElementSpace.calcGradFaceBubbleN( q, stack.X, dBubbleNdX ); + + real64 biotCoefficient; + m_constitutiveUpdate.getBiotCoefficient( k, biotCoefficient ); + + // Assemble the strain-displacement matrix for bubble functions: B_b + real64 strainBubbleMatrix[6][nBubbleUdof]; + solidMechanicsConformingContactKernelsHelper:: + assembleStrainOperator< 6, nBubbleUdof, numFacesPerElem >( strainBubbleMatrix, dBubbleNdX ); + + // Compute biotPressure = -biot * p * I (Voigt notation: [1,1,1,0,0,0]) + real64 biotPressure[6] = {0}; + LvArray::tensorOps::symAddIdentity< 3 >( biotPressure, -biotCoefficient * stack.pLocal[0] ); + + // Compute biotIdentity = -biot * I (for Jacobian) + real64 biotIdentity[6] = {0}; + LvArray::tensorOps::symAddIdentity< 3 >( biotIdentity, -biotCoefficient ); + + // Residual contribution: R_b += B_b^T * biotPressure * detJ + // (note: scaledAdd with -detJ because biotPressure already has negative sign) + real64 Rb_gauss[nBubbleUdof]; + LvArray::tensorOps::Ri_eq_AjiBj< nBubbleUdof, 6 >( Rb_gauss, strainBubbleMatrix, biotPressure ); + LvArray::tensorOps::scaledAdd< nBubbleUdof >( stack.localRb, Rb_gauss, -detJ ); + + // Jacobian contribution: dR_b/dP += B_b^T * biotIdentity * detJ + real64 dRbdP_gauss[nBubbleUdof]; + LvArray::tensorOps::Ri_eq_AjiBj< nBubbleUdof, 6 >( dRbdP_gauss, strainBubbleMatrix, biotIdentity ); + LvArray::tensorOps::scaledAdd< nBubbleUdof >( stack.localdRbdP, dRbdP_gauss, -detJ ); + } + + /** + * @brief Complete - assemble residual and Jacobian into global system + */ + GEOS_HOST_DEVICE + inline + real64 complete( localIndex const kk, + StackVariables & stack ) const + { + localIndex const parentFaceIndex = m_elemsToFaces[kk][1]; + + // Extract only the components corresponding to the parent face + // on which the bubble function was applied. + real64 localRb[3]; + real64 localdRbdP[3]; + for( localIndex i = 0; i < 3; ++i ) + { + localRb[i] = stack.localRb[parentFaceIndex*3+i]; + localdRbdP[i] = stack.localdRbdP[parentFaceIndex*3+i]; + } + + for( localIndex i=0; i < 3; ++i ) + { + localIndex const dof = LvArray::integerConversion< localIndex >( stack.bEqnRowIndices[ i ] ); + + if( dof < 0 || dof >= m_matrix.numRows() ) + continue; + + // Add residual + RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[dof], localRb[i] ); + + // Add Jacobian dR_b/dP + m_matrix.template addToRow< parallelDeviceAtomic >( dof, + &stack.pColIndex, + &localdRbdP[i], + 1 ); + } + + return 0.0; + } + +protected: + + /// The array containing the nodal position array. + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const m_X; + + /// The global degree of freedom number of bubble + arrayView1d< globalIndex const > const m_bDofNumber; + + /// The global degree of freedom number of pressure + arrayView1d< globalIndex const > const m_pDofNumber; + + /// The array containing the list of bubble elements. + arrayView1d< localIndex const > const m_bubbleElems; + + /// The array containing the element to bubble face map. + arrayView2d< localIndex const > const m_elemsToFaces; + + /// The array containing the pressure of each element. + arrayView1d< real64 const > const m_pressure; + +}; + +/// The factory used to construct MatrixPressureBubbleKernels. +using MatrixPressureBubbleFactory = finiteElement::KernelFactory< MatrixPressureBubbleKernels, + arrayView1d< globalIndex const > const, + arrayView1d< globalIndex const > const, + globalIndex const, + CRSMatrixView< real64, globalIndex const > const, + arrayView1d< real64 > const, + real64 const, + string const & >; + +} // namespace poromechanicsMatrixBubbleKernels + +} // namespace geos + +#endif /* GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALMKERNELS_HPP_ */ From edac1f968c67fbfa08bce42566a19be95e84dd61 Mon Sep 17 00:00:00 2001 From: npillardou Date: Wed, 11 Feb 2026 11:46:05 +0100 Subject: [PATCH 09/22] Correct sign convention + hardcoded constitutive --- ...PhasePoromechanicsConformingFracturesALM.cpp | 15 +++++---------- ...romechanicsConformingFracturesALMKernels.hpp | 17 +++++++---------- 2 files changed, 12 insertions(+), 20 deletions(-) diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp index 65da565f2b4..aeb9f8f1b5d 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp @@ -912,13 +912,11 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, // Accumulation derivative w.r.t. nodal displacement // dR_accum/du = density * unitNormal^T * Atu = density * area * dAperturedU // (dAperturedU already has 1/area factor, so multiply by area to cancel it) - // Note: The ALM Atu matrix has opposite sign convention compared to LM's dAper/dU, - // so we need to negate the result to match the LM sign convention. if( isFractureOpen ) { for( localIndex j = 0; j < numUdofs; ++j ) { - dRdU( j ) = -density[kfe][0] * dAperturedU( kfe, j ) * area[kfe]; + dRdU( j ) = density[kfe][0] * dAperturedU( kfe, j ) * area[kfe]; } localIndex const localRow = LvArray::integerConversion< localIndex >( elemDOF[0] - rankOffset ); @@ -965,11 +963,10 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, } // dR_flux/du = dR_flux/dAper * dAper/du (pre-computed for element kfe2) - // Note: Negate to match LM sign convention (ALM Atu has opposite sign) stackArray1d< real64, maxNumUdofs > dRdU2( maxNumUdofs ); for( localIndex j = 0; j < numUdofs2; ++j ) { - dRdU2( j ) = -dR_dAper * dAperturedU( kfe2, j ); + dRdU2( j ) = dR_dAper * dAperturedU( kfe2, j ); } localIndex const localRow = LvArray::integerConversion< localIndex >( elemDOF[0] - rankOffset ); @@ -999,13 +996,12 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, // Accumulation derivative w.r.t. bubble displacement // dR_accum/db = density * unitNormal^T * Atb = density * area * dAperturedB // (dAperturedB already has 1/area factor, so multiply by area to cancel it) - // Note: Negate to match LM sign convention (ALM Atb has opposite sign) if( isFractureOpen ) { real64 dRdB[numBdofs]; for( localIndex j = 0; j < numBdofs; ++j ) { - dRdB[j] = -density[kfe][0] * dAperturedB( kfe, j ) * area[kfe]; + dRdB[j] = density[kfe][0] * dAperturedB( kfe, j ) * area[kfe]; } localIndex const localRow = LvArray::integerConversion< localIndex >( elemDOF[0] - rankOffset ); @@ -1041,11 +1037,10 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, } // dR_flux/db = dR_flux/dAper * dAper/db (pre-computed for element kfe2) - // Note: Negate to match LM sign convention (ALM Atb has opposite sign) real64 dRdB2[numBdofs]; for( localIndex j = 0; j < numBdofs; ++j ) { - dRdB2[j] = -dR_dAper * dAperturedB( kfe2, j ); + dRdB2[j] = dR_dAper * dAperturedB( kfe2, j ); } localIndex const localRow = LvArray::integerConversion< localIndex >( elemDOF[0] - rankOffset ); @@ -1222,7 +1217,7 @@ assembleMatrixPressureBubbleContribution( real64 const dt, real64 maxResidual = finiteElement::regionBasedKernelApplication < parallelDevicePolicy<>, - constitutive::PorousSolid< constitutive::ElasticIsotropic, constitutive::ConstantPermeability >, + constitutive::PorousSolidBase, CellElementSubRegion >( mesh, poromechanicsRegionNames, this->solidMechanicsSolver()->getDiscretizationName(), diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp index f069b7f7b90..14b89821395 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp @@ -474,19 +474,16 @@ class AssembleFluidMassResidualDerivativeWrtDisplacement : } // Compute Jacobian for pressure equation w.r.t. displacement: - // dR_p/du = density * unitNormal^T * Atu + // dR_p/du = density * unitNormal^T * Atu = density * area * dAperture/dU // This is: dRpdU[j] = sum_i( unitNormal[i] * Atu[i][j] ) * density // Use Ri_eq_AjiBj for transpose: R[i] = sum_j(A[j][i] * B[j]) - // Note: The ALM Atu matrix has opposite sign convention compared to LM's dAper/dU, - // so we negate the result to match the LM sign convention. LvArray::tensorOps::Ri_eq_AjiBj< numUdofs, 3 >( stack.dRpdU, stack.localAtu, stack.unitNormal ); - LvArray::tensorOps::scale< numUdofs >( stack.dRpdU, -m_density[k] ); + LvArray::tensorOps::scale< numUdofs >( stack.dRpdU, m_density[k] ); // Compute Jacobian for pressure equation w.r.t. bubble displacement: // dR_p/db = density * unitNormal^T * Atb - // Note: Same sign convention fix as for displacement. LvArray::tensorOps::Ri_eq_AjiBj< numBdofs, 3 >( stack.dRpdB, stack.localAtb, stack.unitNormal ); - LvArray::tensorOps::scale< numBdofs >( stack.dRpdB, -m_density[k] ); + LvArray::tensorOps::scale< numBdofs >( stack.dRpdB, m_density[k] ); // Assemble into matrix localIndex const localRow = LvArray::integerConversion< localIndex >( stack.pressureRowIndex ); @@ -546,10 +543,10 @@ using AssembleFluidMassResidualDerivativeWrtDisplacementFactory = * d(aperture)/du = (1/area) * unitNormal^T * Atu * d(aperture)/db = (1/area) * unitNormal^T * Atb * - * IMPORTANT: The ALM Atu/Atb matrices have opposite sign convention compared to the LM - * dAper/dU formulas. When using these pre-computed values, they must be negated to match - * the LM sign convention. The assembly in assembleFluidMassResidualDerivativeWrtDisplacement - * handles this negation. + * The Atu/Atb matrices encode the sign convention directly: face 0 gets -N*detJ + * and face 1 gets +N*detJ. This means (1/area) * Atu^T * unitNormal already gives + * the correct dAperture/dU matching the LM convention (negative for face 0 nodes, + * positive for face 1 nodes). */ template< typename CONSTITUTIVE_TYPE, typename FE_TYPE > From 889de70c715539b551175f0178a3eef1a23d8a7b Mon Sep 17 00:00:00 2001 From: npillardou Date: Thu, 12 Feb 2026 11:23:53 +0100 Subject: [PATCH 10/22] Correct sparsity pattern --- ...asePoromechanicsConformingFracturesALM.cpp | 53 +++++++++++-------- 1 file changed, 31 insertions(+), 22 deletions(-) diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp index aeb9f8f1b5d..ddde6807a4e 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp @@ -95,43 +95,52 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupSystem( if( setSparsity ) { - // Get sparsity pattern from flow solver (could include wells) - SparsityPattern< globalIndex > patternOriginal; - this->flowSolver()->setSparsityPattern( domain, dofManager, localMatrix, patternOriginal ); - - // Get the original row lengths - array1d< localIndex > rowLengths( patternOriginal.numRows() ); - for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) + // Step 1: Get the flow sparsity pattern (includes DofManager diagonal blocks + flux stencil connections) + SparsityPattern< globalIndex > flowPattern; + this->flowSolver()->setSparsityPattern( domain, dofManager, localMatrix, flowPattern ); + + // Step 2: Get the solid mechanics sparsity pattern (includes DofManager diagonal blocks + mechanics coupling) + // Note: setSparsityPattern replaces its output pattern, so we use a separate variable. + SparsityPattern< globalIndex > mechPattern; + this->solidMechanicsSolver()->setSparsityPattern( domain, dofManager, localMatrix, mechPattern ); + + // Step 3: Compute combined row lengths (overestimate due to shared DofManager diagonal entries) + array1d< localIndex > rowLengths( flowPattern.numRows() ); + for( localIndex localRow = 0; localRow < flowPattern.numRows(); ++localRow ) { - rowLengths[localRow] = patternOriginal.numNonZeros( localRow ); + rowLengths[localRow] = flowPattern.numNonZeros( localRow ) + mechPattern.numNonZeros( localRow ); } - // Add the number of nonzeros induced by coupling + // Step 4: Add the number of nonzeros induced by flow-mechanics coupling addTransmissibilityCouplingNNZ( domain, dofManager, rowLengths.toView() ); addPressureForceCouplingNNZ( domain, dofManager, rowLengths.toView() ); addMatrixPressureBubbleCouplingNNZ( domain, dofManager, rowLengths.toView() ); - // Create a new pattern with enough capacity for coupled matrix + // Step 5: Create a new pattern with enough capacity for the full coupled matrix SparsityPattern< globalIndex > pattern; - pattern.resizeFromRowCapacities< parallelHostPolicy >( patternOriginal.numRows(), - patternOriginal.numColumns(), + pattern.resizeFromRowCapacities< parallelHostPolicy >( flowPattern.numRows(), + flowPattern.numColumns(), rowLengths.data() ); - // Copy the original nonzeros - for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) + // Step 6: Copy flow pattern entries + for( localIndex localRow = 0; localRow < flowPattern.numRows(); ++localRow ) { - globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); - pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow ) ); + globalIndex const * cols = flowPattern.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + flowPattern.numNonZeros( localRow ) ); } - // Add the nonzeros from coupling + // Step 7: Copy mechanics pattern entries (duplicates with flow diagonal blocks handled by insertNonZeros) + for( localIndex localRow = 0; localRow < mechPattern.numRows(); ++localRow ) + { + globalIndex const * cols = mechPattern.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + mechPattern.numNonZeros( localRow ) ); + } + + // Step 8: Add the flow-mechanics coupling patterns addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView() ); addPressureForceCouplingPattern( domain, dofManager, pattern.toView() ); addMatrixPressureBubbleCouplingPattern( domain, dofManager, pattern.toView() ); - // Assemble the full sparsity pattern using the solid mechanics solver - this->solidMechanicsSolver()->setSparsityPattern( domain, dofManager, localMatrix, pattern ); - // Set up the derivative flux residual matrix setUpDflux_dApertureMatrix( domain, dofManager, localMatrix ); @@ -1082,7 +1091,7 @@ addMatrixPressureBubbleCouplingNNZ( DomainPartition const & domain, elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, CellElementSubRegion const & subRegion ) { - if( !subRegion.hasWrapper( "bubbleElementsList" ) ) + if( !subRegion.hasWrapper( CellElementSubRegion::viewKeyStruct::bubbleCellsString() ) ) return; arrayView1d< localIndex const > const bubbleElems = subRegion.bubbleElementsList(); @@ -1134,7 +1143,7 @@ addMatrixPressureBubbleCouplingPattern( DomainPartition const & domain, elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, CellElementSubRegion const & subRegion ) { - if( !subRegion.hasWrapper( "bubbleElementsList" ) ) + if( !subRegion.hasWrapper( CellElementSubRegion::viewKeyStruct::bubbleCellsString() ) ) return; arrayView1d< localIndex const > const bubbleElems = subRegion.bubbleElementsList(); From b87790363d9025ec71c9c1d50213b7c2907f05ff Mon Sep 17 00:00:00 2001 From: npillardou Date: Thu, 12 Feb 2026 13:25:35 +0100 Subject: [PATCH 11/22] Change MGR config --- ...asePoromechanicsConformingFracturesALM.hpp | 65 ++++++++++--------- 1 file changed, 35 insertions(+), 30 deletions(-) diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp index e9fb7db8bd6..5f4f908b140 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -14,7 +14,7 @@ */ /** - * @file SinglePhasePoromechanicsConformingFractures.hpp + * @file SinglePhasePoromechanicsConformingFracturesALM.hpp */ #ifndef GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_ @@ -32,23 +32,25 @@ namespace mgr { /** - * @brief SinglePhasePoromechanicsConformingFractures strategy. + * @brief SinglePhasePoromechanicsConformingFracturesALM strategy. + * + * Labels assigned by DofManager (PoromechanicsSolver::setupDofs calls + * solidMechanicsSolver first, then flowSolver): * * dofLabel: 0 = displacement, x-component * dofLabel: 1 = displacement, y-component * dofLabel: 2 = displacement, z-component - * dofLabel: 3 = displacement bubble function, x-component - * dofLabel: 4 = displacement bubble function, y-component - * dofLabel: 5 = displacement bubble function, z-component + * dofLabel: 3 = bubble displacement, x-component + * dofLabel: 4 = bubble displacement, y-component + * dofLabel: 5 = bubble displacement, z-component * dofLabel: 6 = pressure (cell elem + fracture elems) - * * Ingredients: - * 1. Level 1: F-points displacement (3,4,5), C-points pressure (0,1,2,6) - * 2. Level 2: F-points displacement (0,1,2), C-points pressure (6) - * 2. F-points smoother: BoomerAMG, single V-cycle - * 3. C-points coarse-grid/Schur complement solver: BoomerAMG - * 4. Global smoother: none + * 1. Level 0: F-points bubble (3,4,5), C-points displacement+pressure (0,1,2,6) + * 2. Level 1: F-points displacement (0,1,2), C-points pressure (6) + * 3. F-points smoother at Level 1: BoomerAMG, single V-cycle + * 4. C-points coarse-grid/Schur complement solver: BoomerAMG + * 5. Global smoother: none */ class SinglePhasePoromechanicsConformingFracturesALM : public MGRStrategyBase< 2 > { @@ -61,32 +63,33 @@ class SinglePhasePoromechanicsConformingFracturesALM : public MGRStrategyBase< 2 : MGRStrategyBase( 7 ) { - // we keep u and p + // Level 0: keep displacement (0,1,2) + pressure (6), eliminate bubble (3,4,5) m_labels[0].push_back( 0 ); m_labels[0].push_back( 1 ); m_labels[0].push_back( 2 ); m_labels[0].push_back( 6 ); - // we keep p + // Level 1: keep pressure (6), eliminate displacement (0,1,2) m_labels[1].push_back( 6 ); setupLabels(); - // Level 0 + // Level 0: eliminate bubble DOFs (3,4,5), keep displacement (0,1,2) + pressure (6) m_levelFRelaxType[0] = MGRFRelaxationType::l1jacobi; m_levelFRelaxIters[0] = 1; m_levelGlobalSmootherType[0] = MGRGlobalSmootherType::none; m_levelGlobalSmootherIters[0] = 0; - m_levelInterpType[0] = MGRInterpolationType::blockJacobi; + m_levelInterpType[0] = MGRInterpolationType::injection; m_levelRestrictType[0] = MGRRestrictionType::injection; m_levelCoarseGridMethod[0] = MGRCoarseGridMethod::galerkin; // Level 1 - m_levelFRelaxType[1] = MGRFRelaxationType::amgVCycle; - m_levelFRelaxIters[1] = 1; - m_levelGlobalSmootherType[1] = MGRGlobalSmootherType::none; - m_levelInterpType[1] = MGRInterpolationType::jacobi; - m_levelRestrictType[1] = MGRRestrictionType::injection; - m_levelCoarseGridMethod[1] = MGRCoarseGridMethod::nonGalerkin; + m_levelFRelaxType[1] = MGRFRelaxationType::amgVCycle; + m_levelFRelaxIters[1] = 1; + m_levelGlobalSmootherType[1] = MGRGlobalSmootherType::ilu0; + m_levelGlobalSmootherIters[1] = 1; + m_levelInterpType[1] = MGRInterpolationType::jacobi; + m_levelRestrictType[1] = MGRRestrictionType::injection; + m_levelCoarseGridMethod[1] = MGRCoarseGridMethod::nonGalerkin; } @@ -95,21 +98,23 @@ class SinglePhasePoromechanicsConformingFracturesALM : public MGRStrategyBase< 2 * @param precond preconditioner wrapper * @param mgrData auxiliary MGR data */ - void setup( LinearSolverParameters::MGR const & mgrParams, + void setup( LinearSolverParameters::MGR const &, HyprePrecWrapper & precond, HypreMGRData & mgrData ) { setReduction( precond, mgrData ); // Configure the BoomerAMG solver used as F-relaxation for the second level - // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGCreate( &mgrData.mechSolver.ptr ) ); - // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetTol( mgrData.mechSolver.ptr, 0.0 ) ); - // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxIter( mgrData.mechSolver.ptr, 1 ) ); - // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxRowSum( mgrData.mechSolver.ptr, 1.0 ) ); - // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetStrongThreshold( mgrData.mechSolver.ptr, 0.6 ) ); - // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPrintLevel( mgrData.mechSolver.ptr, 0 ) ); - // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumFunctions( mgrData.mechSolver.ptr, 3 ) ); - setDisplacementAMG(mgrData.mechSolver, mgrParams.separateComponents); + // Note: we do NOT use setDisplacementAMG() here because HYPRE_MGRSetFSolverAtLevel + // transfers ownership of the solver to MGR. MGR will destroy it during HYPRE_MGRDestroy. + // Setting mechSolver.destroy would cause a double-free in HyprePreconditioner::clear(). + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGCreate( &mgrData.mechSolver.ptr ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetTol( mgrData.mechSolver.ptr, 0.0 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxIter( mgrData.mechSolver.ptr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxRowSum( mgrData.mechSolver.ptr, 1.0 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetStrongThreshold( mgrData.mechSolver.ptr, 0.6 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPrintLevel( mgrData.mechSolver.ptr, 0 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumFunctions( mgrData.mechSolver.ptr, 3 ) ); #if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCoarsenType( mgrData.mechSolver.ptr, hypre::getAMGCoarseningType( LinearSolverParameters::AMG::CoarseningType::PMIS ) ) ); From 02cf3e1b48b275c35ba235288ba6dbf2a24f28f8 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Tue, 10 Mar 2026 23:50:57 +0100 Subject: [PATCH 12/22] Revert "Change MGR config" This reverts commit b87790363d9025ec71c9c1d50213b7c2907f05ff. --- ...asePoromechanicsConformingFracturesALM.hpp | 65 +++++++++---------- 1 file changed, 30 insertions(+), 35 deletions(-) diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp index 5f4f908b140..e9fb7db8bd6 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -14,7 +14,7 @@ */ /** - * @file SinglePhasePoromechanicsConformingFracturesALM.hpp + * @file SinglePhasePoromechanicsConformingFractures.hpp */ #ifndef GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_ @@ -32,25 +32,23 @@ namespace mgr { /** - * @brief SinglePhasePoromechanicsConformingFracturesALM strategy. - * - * Labels assigned by DofManager (PoromechanicsSolver::setupDofs calls - * solidMechanicsSolver first, then flowSolver): + * @brief SinglePhasePoromechanicsConformingFractures strategy. * * dofLabel: 0 = displacement, x-component * dofLabel: 1 = displacement, y-component * dofLabel: 2 = displacement, z-component - * dofLabel: 3 = bubble displacement, x-component - * dofLabel: 4 = bubble displacement, y-component - * dofLabel: 5 = bubble displacement, z-component + * dofLabel: 3 = displacement bubble function, x-component + * dofLabel: 4 = displacement bubble function, y-component + * dofLabel: 5 = displacement bubble function, z-component * dofLabel: 6 = pressure (cell elem + fracture elems) + * * Ingredients: - * 1. Level 0: F-points bubble (3,4,5), C-points displacement+pressure (0,1,2,6) - * 2. Level 1: F-points displacement (0,1,2), C-points pressure (6) - * 3. F-points smoother at Level 1: BoomerAMG, single V-cycle - * 4. C-points coarse-grid/Schur complement solver: BoomerAMG - * 5. Global smoother: none + * 1. Level 1: F-points displacement (3,4,5), C-points pressure (0,1,2,6) + * 2. Level 2: F-points displacement (0,1,2), C-points pressure (6) + * 2. F-points smoother: BoomerAMG, single V-cycle + * 3. C-points coarse-grid/Schur complement solver: BoomerAMG + * 4. Global smoother: none */ class SinglePhasePoromechanicsConformingFracturesALM : public MGRStrategyBase< 2 > { @@ -63,33 +61,32 @@ class SinglePhasePoromechanicsConformingFracturesALM : public MGRStrategyBase< 2 : MGRStrategyBase( 7 ) { - // Level 0: keep displacement (0,1,2) + pressure (6), eliminate bubble (3,4,5) + // we keep u and p m_labels[0].push_back( 0 ); m_labels[0].push_back( 1 ); m_labels[0].push_back( 2 ); m_labels[0].push_back( 6 ); - // Level 1: keep pressure (6), eliminate displacement (0,1,2) + // we keep p m_labels[1].push_back( 6 ); setupLabels(); - // Level 0: eliminate bubble DOFs (3,4,5), keep displacement (0,1,2) + pressure (6) + // Level 0 m_levelFRelaxType[0] = MGRFRelaxationType::l1jacobi; m_levelFRelaxIters[0] = 1; m_levelGlobalSmootherType[0] = MGRGlobalSmootherType::none; m_levelGlobalSmootherIters[0] = 0; - m_levelInterpType[0] = MGRInterpolationType::injection; + m_levelInterpType[0] = MGRInterpolationType::blockJacobi; m_levelRestrictType[0] = MGRRestrictionType::injection; m_levelCoarseGridMethod[0] = MGRCoarseGridMethod::galerkin; // Level 1 - m_levelFRelaxType[1] = MGRFRelaxationType::amgVCycle; - m_levelFRelaxIters[1] = 1; - m_levelGlobalSmootherType[1] = MGRGlobalSmootherType::ilu0; - m_levelGlobalSmootherIters[1] = 1; - m_levelInterpType[1] = MGRInterpolationType::jacobi; - m_levelRestrictType[1] = MGRRestrictionType::injection; - m_levelCoarseGridMethod[1] = MGRCoarseGridMethod::nonGalerkin; + m_levelFRelaxType[1] = MGRFRelaxationType::amgVCycle; + m_levelFRelaxIters[1] = 1; + m_levelGlobalSmootherType[1] = MGRGlobalSmootherType::none; + m_levelInterpType[1] = MGRInterpolationType::jacobi; + m_levelRestrictType[1] = MGRRestrictionType::injection; + m_levelCoarseGridMethod[1] = MGRCoarseGridMethod::nonGalerkin; } @@ -98,23 +95,21 @@ class SinglePhasePoromechanicsConformingFracturesALM : public MGRStrategyBase< 2 * @param precond preconditioner wrapper * @param mgrData auxiliary MGR data */ - void setup( LinearSolverParameters::MGR const &, + void setup( LinearSolverParameters::MGR const & mgrParams, HyprePrecWrapper & precond, HypreMGRData & mgrData ) { setReduction( precond, mgrData ); // Configure the BoomerAMG solver used as F-relaxation for the second level - // Note: we do NOT use setDisplacementAMG() here because HYPRE_MGRSetFSolverAtLevel - // transfers ownership of the solver to MGR. MGR will destroy it during HYPRE_MGRDestroy. - // Setting mechSolver.destroy would cause a double-free in HyprePreconditioner::clear(). - GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGCreate( &mgrData.mechSolver.ptr ) ); - GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetTol( mgrData.mechSolver.ptr, 0.0 ) ); - GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxIter( mgrData.mechSolver.ptr, 1 ) ); - GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxRowSum( mgrData.mechSolver.ptr, 1.0 ) ); - GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetStrongThreshold( mgrData.mechSolver.ptr, 0.6 ) ); - GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPrintLevel( mgrData.mechSolver.ptr, 0 ) ); - GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumFunctions( mgrData.mechSolver.ptr, 3 ) ); + // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGCreate( &mgrData.mechSolver.ptr ) ); + // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetTol( mgrData.mechSolver.ptr, 0.0 ) ); + // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxIter( mgrData.mechSolver.ptr, 1 ) ); + // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxRowSum( mgrData.mechSolver.ptr, 1.0 ) ); + // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetStrongThreshold( mgrData.mechSolver.ptr, 0.6 ) ); + // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPrintLevel( mgrData.mechSolver.ptr, 0 ) ); + // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumFunctions( mgrData.mechSolver.ptr, 3 ) ); + setDisplacementAMG(mgrData.mechSolver, mgrParams.separateComponents); #if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCoarsenType( mgrData.mechSolver.ptr, hypre::getAMGCoarseningType( LinearSolverParameters::AMG::CoarseningType::PMIS ) ) ); From 1759535c06b2c568cff7c591fb7d8787fc95fcb7 Mon Sep 17 00:00:00 2001 From: npillardou Date: Fri, 10 Jul 2026 16:26:11 +0200 Subject: [PATCH 13/22] Add missing A^bp coupling term --- ...asePoromechanicsConformingFracturesALM.cpp | 27 +++++++- ...mechanicsConformingFracturesALMKernels.hpp | 67 +++++++++++++++++-- 2 files changed, 87 insertions(+), 7 deletions(-) diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp index ddde6807a4e..0ff22a0961e 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp @@ -1083,6 +1083,7 @@ addMatrixPressureBubbleCouplingNNZ( DomainPartition const & domain, ElementRegionManager const & elemManager = mesh.getElemManager(); string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); + string const flowDofKey = dofManager.getKey( m_pressureKey ); arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); globalIndex const rankOffset = dofManager.rankOffset(); @@ -1096,12 +1097,14 @@ addMatrixPressureBubbleCouplingNNZ( DomainPartition const & domain, arrayView1d< localIndex const > const bubbleElems = subRegion.bubbleElementsList(); arrayView2d< localIndex const > const elemsToFaces = subRegion.faceElementsList(); + arrayView1d< globalIndex const > const pressureDofNumber = subRegion.getReference< array1d< globalIndex > >( flowDofKey ); forAll< serialPolicy >( bubbleElems.size(), [=, &rowLengths]( localIndex const kk ) { + localIndex const k = bubbleElems[kk]; localIndex const faceIndex = elemsToFaces[kk][0]; - // Add 1 pressure column for each of the 3 bubble DOFs + // (bubble_row, pressure_col): 1 pressure column for each of the 3 bubble DOFs for( localIndex i = 0; i < 3; ++i ) { globalIndex const rowNumber = bubbleDofNumber[faceIndex] + i - rankOffset; @@ -1110,6 +1113,13 @@ addMatrixPressureBubbleCouplingNNZ( DomainPartition const & domain, rowLengths[rowNumber] += 1; // One pressure DOF from matrix cell } } + + // (pressure_row, bubble_col): the matrix cell pressure couples to its 3 bubble DOFs (A_pb) + globalIndex const pRow = pressureDofNumber[k] - rankOffset; + if( pRow >= 0 && pRow < rowLengths.size() ) + { + rowLengths[pRow] += 3; // Three bubble DOFs + } } ); } ); } ); @@ -1156,7 +1166,7 @@ addMatrixPressureBubbleCouplingPattern( DomainPartition const & domain, localIndex const faceIndex = elemsToFaces[kk][0]; globalIndex const pressureColIndex = pressureDofNumber[k]; - // Add pattern for (bubble_row, pressure_col) + // (bubble_row, pressure_col) : A_bp for( localIndex i = 0; i < 3; ++i ) { globalIndex const rowIndex = bubbleDofNumber[faceIndex] + i - rankOffset; @@ -1165,6 +1175,16 @@ addMatrixPressureBubbleCouplingPattern( DomainPartition const & domain, pattern.insertNonZero( rowIndex, pressureColIndex ); } } + + // (pressure_row, bubble_col) : A_pb -- transpose location + globalIndex const pRow = pressureDofNumber[k] - rankOffset; + if( pRow >= 0 && pRow < pattern.numRows() ) + { + for( localIndex i = 0; i < 3; ++i ) + { + pattern.insertNonZero( pRow, bubbleDofNumber[faceIndex] + i ); + } + } } ); } ); } ); @@ -1222,7 +1242,8 @@ assembleMatrixPressureBubbleContribution( real64 const dt, localMatrix, localRhs, dt, - flowDofKey ); + flowDofKey, + FlowSolverBase::viewKeyStruct::fluidNamesString() ); real64 maxResidual = finiteElement::regionBasedKernelApplication < parallelDevicePolicy<>, diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp index 14b89821395..b05cbc61d35 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp @@ -25,6 +25,7 @@ #include "physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsConformingContactKernelsHelper.hpp" #include "physicsSolvers/solidMechanics/contact/ContactFields.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" +#include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" namespace geos { @@ -810,7 +811,8 @@ class MatrixPressureBubbleKernels : CRSMatrixView< real64, globalIndex const > const inputMatrix, arrayView1d< real64 > const inputRhs, real64 const inputDt, - string const & pressureDofKey ): + string const pressureDofKey, + string const fluidModelKey ): Base( nodeManager, edgeManager, faceManager, @@ -828,7 +830,10 @@ class MatrixPressureBubbleKernels : m_pDofNumber( elementSubRegion.template getReference< array1d< globalIndex > >( pressureDofKey ) ), m_bubbleElems( elementSubRegion.bubbleElementsList() ), m_elemsToFaces( elementSubRegion.faceElementsList() ), - m_pressure( elementSubRegion.template getField< fields::flow::pressure >().toViewConst() ) + m_pressure( elementSubRegion.template getField< fields::flow::pressure >().toViewConst() ), + m_incrBubbleDisp( faceManager.getField< fields::contact::incrementalBubbleDisplacement >().toViewConst() ), + m_fluidDensity( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( + elementSubRegion.template getReference< string >( fluidModelKey ) ).density() ) {} //*************************************************************************** @@ -855,9 +860,12 @@ class MatrixPressureBubbleKernels : GEOS_HOST_DEVICE StackVariables(): bEqnRowIndices{}, + bColIndices{}, pColIndex( 0 ), + pRowIndex( 0 ), localRb{}, localdRbdP{}, + localdRpdB{}, X{ {} }, pLocal{} {} @@ -865,15 +873,24 @@ class MatrixPressureBubbleKernels : /// C-array storage for the element local row degrees of freedom (bubble). globalIndex bEqnRowIndices[3]; - /// Column index for pressure DOF + /// C-array storage for the (global) bubble column degrees of freedom (for A_pb). + globalIndex bColIndices[3]; + + /// Column index for pressure DOF (for A_bp) globalIndex pColIndex; + /// Row index (local) for pressure DOF (for A_pb) + globalIndex pRowIndex; + /// C-array storage for the element local Rb residual vector. real64 localRb[numBubbleUdofs]; /// C-array storage for the element local dRb/dP Jacobian. real64 localdRbdP[numBubbleUdofs]; + /// C-array storage for the element local dRp/dB Jacobian (A_pb bulk term). + real64 localdRpdB[numBubbleUdofs]; + /// local nodal coordinates real64 X[ numNodesPerElem ][ 3 ]; @@ -944,9 +961,11 @@ class MatrixPressureBubbleKernels : for( int i=0; i<3; ++i ) { stack.bEqnRowIndices[i] = m_bDofNumber[localFaceIndex] + i - m_dofRankOffset; + stack.bColIndices[i] = m_bDofNumber[localFaceIndex] + i; // global column for A_pb } stack.pColIndex = m_pDofNumber[k]; + stack.pRowIndex = m_pDofNumber[k] - m_dofRankOffset; stack.pLocal[0] = m_pressure( k ); } @@ -995,6 +1014,15 @@ class MatrixPressureBubbleKernels : real64 dRbdP_gauss[nBubbleUdof]; LvArray::tensorOps::Ri_eq_AjiBj< nBubbleUdof, 6 >( dRbdP_gauss, strainBubbleMatrix, biotIdentity ); LvArray::tensorOps::scaledAdd< nBubbleUdof >( stack.localdRbdP, dRbdP_gauss, -detJ ); + + // ---- A_pb^Omega : transpose Biot term for the fluid-mass equation ---- + // The bubble mode contributes to the cell volumetric strain, hence to the + // fluid mass storage. With dPorosity_dVolStrain = biot (BiotPorosity), + // dFluidMassIncrement_dVolStrainIncrement = biot * fluidDensity, + // so dR_p/db = fluidDensity * dR_b/dp (same geometric factor biot*div(beta), + // times the fluid density that puts it in mass units). We reuse dRbdP_gauss. + real64 const rhof = m_fluidDensity( k, 0 ); + LvArray::tensorOps::scaledAdd< nBubbleUdof >( stack.localdRpdB, dRbdP_gauss, -detJ * rhof ); } /** @@ -1011,12 +1039,15 @@ class MatrixPressureBubbleKernels : // on which the bubble function was applied. real64 localRb[3]; real64 localdRbdP[3]; + real64 localdRpdB[3]; for( localIndex i = 0; i < 3; ++i ) { localRb[i] = stack.localRb[parentFaceIndex*3+i]; localdRbdP[i] = stack.localdRbdP[parentFaceIndex*3+i]; + localdRpdB[i] = stack.localdRpdB[parentFaceIndex*3+i]; } + // ---- A_bp : bubble-momentum row (dR_b/dP) ---- for( localIndex i=0; i < 3; ++i ) { localIndex const dof = LvArray::integerConversion< localIndex >( stack.bEqnRowIndices[ i ] ); @@ -1034,6 +1065,27 @@ class MatrixPressureBubbleKernels : 1 ); } + // ---- A_pb^Omega : fluid-mass row (dR_p/db) + storage residual ---- + localIndex const pRow = LvArray::integerConversion< localIndex >( stack.pRowIndex ); + if( pRow >= 0 && pRow < m_matrix.numRows() ) + { + // Residual: R_p += fluidDensity * biot * tr(grad beta) * delta_b (bubble storage), + // using the incremental bubble displacement of the bubble face. + localIndex const bFace = m_elemsToFaces[kk][0]; + real64 rp = 0.0; + for( localIndex i = 0; i < 3; ++i ) + { + rp += localdRpdB[i] * m_incrBubbleDisp[bFace][i]; + } + RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[pRow], rp ); + + // Jacobian dR_p/db (three bubble columns of the parent face) + m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( pRow, + stack.bColIndices, + localdRpdB, + 3 ); + } + return 0.0; } @@ -1057,6 +1109,12 @@ class MatrixPressureBubbleKernels : /// The array containing the pressure of each element. arrayView1d< real64 const > const m_pressure; + /// Incremental bubble displacement (face field) -- used for the fluid-mass storage residual. + arrayView2d< real64 const > const m_incrBubbleDisp; + + /// Fluid density [elem][q]. + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_fluidDensity; + }; /// The factory used to construct MatrixPressureBubbleKernels. @@ -1067,7 +1125,8 @@ using MatrixPressureBubbleFactory = finiteElement::KernelFactory< MatrixPressure CRSMatrixView< real64, globalIndex const > const, arrayView1d< real64 > const, real64 const, - string const & >; + string const, + string const >; } // namespace poromechanicsMatrixBubbleKernels From 11d61333940123b013b525acc05d2b35ab80b79a Mon Sep 17 00:00:00 2001 From: "Victor A. P. Magri" Date: Thu, 6 Aug 2026 17:17:25 -0400 Subject: [PATCH 14/22] Fix VTK mesh redistribution --- .../mesh/generators/VTKMeshGeneratorTools.cpp | 7 + .../mesh/generators/VTKUtilities.cpp | 169 ++++++++++++++++++ 2 files changed, 176 insertions(+) diff --git a/src/coreComponents/mesh/generators/VTKMeshGeneratorTools.cpp b/src/coreComponents/mesh/generators/VTKMeshGeneratorTools.cpp index dac362159a9..a2a6e3a3a77 100644 --- a/src/coreComponents/mesh/generators/VTKMeshGeneratorTools.cpp +++ b/src/coreComponents/mesh/generators/VTKMeshGeneratorTools.cpp @@ -19,6 +19,8 @@ #include "VTKMeshGeneratorTools.hpp" +#include "LvArray/src/system.hpp" + #include #include #include @@ -33,6 +35,11 @@ vtkSmartPointer< vtkUnstructuredGrid > redistribute( vtkPartitionedDataSet & localParts, MPI_Comm mpiComm ) { + // VTK's XML writer, used by DIY to serialize the exchanged grids, can + // evaluate empty arrays while reporting progress. Do not let those + // internal operations trip GEOS' floating-point exception handler. + LvArray::system::FloatingPointExceptionGuard guard; + // The code below is modified from vtkDIYKdTreeUtilities::Exchange(): // https://gitlab.kitware.com/vtk/vtk/-/blob/7037a148605bf9628710d8b729c22f27dd0ede93/Filters/ParallelDIY2/vtkDIYKdTreeUtilities.cxx#L289 // We cannot call that function directly because vtkDIYKdTreeUtilities.hpp is a private header in VTK. diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp index 1dca5f69d5f..cd29e657a75 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp @@ -77,6 +77,8 @@ #endif #include +#include +#include #include @@ -995,6 +997,143 @@ static void classifyCellsByDimension( vtkDataSet & mesh, cells3DIndices.size(), cells2DIndices.size(), numCells ) ); } +using VTKPointCoordinate = std::array< real64, 3 >; + +struct VTKPointCoordinateHash +{ + std::size_t operator()( VTKPointCoordinate const & point ) const + { + std::size_t hash = 0; + for( real64 const coordinate: point ) + { + hash ^= std::hash< real64 >{}( coordinate ) + 0x9e3779b9 + ( hash << 6 ) + ( hash >> 2 ); + } + return hash; + } +}; + +static VTKPointCoordinate getPointCoordinate( vtkDataSet & mesh, vtkIdType const pointId ) +{ + double const * const point = mesh.GetPoint( pointId ); + return { point[0], point[1], point[2] }; +} + +/** + * @brief Find 3D cells whose faces match a 2D cell geometrically. + * + * Some mesh generators keep a separate point for each side of a split + * interface. In that case, the point ids of an embedded 2D cell do not + * match the point ids of either adjacent 3D face, even though the + * coordinates do. This routine performs the co-location lookup through + * point coordinates and then verifies that all points belong to one face. + * + * @param[in] mesh Original mesh + * @param[in] pointIds2D Point ids of the 2D cell + * @param[in] meshIdxToGlobalId3D Map from 3D mesh indices to global ids + * @param[in] pointsByCoordinate All mesh point ids grouped by coordinates + * @param[in,out] pointTo3DCells Cache of 3D cells incident to a mesh point + * @return Global ids of geometrically matching 3D cells + */ +static stdVector< int64_t > find2DTo3DNeighborsByCoordinates( + vtkDataSet & mesh, + vtkIdList * pointIds2D, + stdUnorderedMap< vtkIdType, int64_t > const & meshIdxToGlobalId3D, + stdUnorderedMap< VTKPointCoordinate, stdVector< vtkIdType >, VTKPointCoordinateHash > const & pointsByCoordinate, + stdUnorderedMap< vtkIdType, stdVector< vtkIdType > > & pointTo3DCells ) +{ + localIndex const numPoints = pointIds2D->GetNumberOfIds(); + + stdVector< VTKPointCoordinate > targetCoordinates; + targetCoordinates.reserve( numPoints ); + + stdUnorderedMap< vtkIdType, localIndex > candidateCounts; + + vtkNew< vtkIdList > pointCells; + for( vtkIdType i = 0; i < pointIds2D->GetNumberOfIds(); ++i ) + { + VTKPointCoordinate const coordinate = getPointCoordinate( mesh, pointIds2D->GetId( i ) ); + targetCoordinates.emplace_back( coordinate ); + + auto const coordinateIt = pointsByCoordinate.find( coordinate ); + if( coordinateIt == pointsByCoordinate.end() ) + { + continue; + } + + // A 3D cell can contain more than one mesh point with the same + // coordinate in a degenerate input. Count it only once for this 2D + // point when intersecting the incident-cell lists. + stdUnorderedMap< vtkIdType, bool > cellsAtCoordinate; + for( vtkIdType const meshPointId: coordinateIt->second ) + { + auto pointCellsIt = pointTo3DCells.find( meshPointId ); + if( pointCellsIt == pointTo3DCells.end() ) + { + stdVector< vtkIdType > & cachedCells = pointTo3DCells.get_inserted( meshPointId ); + pointCells->Reset(); + mesh.GetPointCells( meshPointId, pointCells ); + cachedCells.reserve( pointCells->GetNumberOfIds() ); + for( vtkIdType j = 0; j < pointCells->GetNumberOfIds(); ++j ) + { + vtkIdType const cellId = pointCells->GetId( j ); + if( meshIdxToGlobalId3D.count( cellId ) > 0 ) + { + cachedCells.emplace_back( cellId ); + } + } + pointCellsIt = pointTo3DCells.find( meshPointId ); + } + + for( vtkIdType const cellId: pointCellsIt->second ) + { + cellsAtCoordinate.emplace( cellId, true ); + } + } + + for( auto const & cell: cellsAtCoordinate ) + { + ++candidateCounts.get_inserted( cell.first ); + } + } + + std::sort( targetCoordinates.begin(), targetCoordinates.end() ); + + stdVector< int64_t > neighbors; + for( auto const & candidate: candidateCounts ) + { + if( candidate.second != numPoints ) + { + continue; + } + + vtkCell * const cell3D = mesh.GetCell( candidate.first ); + for( localIndex faceIndex = 0; faceIndex < cell3D->GetNumberOfFaces(); ++faceIndex ) + { + vtkCell * const face = cell3D->GetFace( faceIndex ); + if( face->GetNumberOfPoints() != numPoints ) + { + continue; + } + + stdVector< VTKPointCoordinate > faceCoordinates; + faceCoordinates.reserve( numPoints ); + for( vtkIdType i = 0; i < face->GetNumberOfPoints(); ++i ) + { + faceCoordinates.emplace_back( getPointCoordinate( mesh, face->GetPointId( i ) ) ); + } + std::sort( faceCoordinates.begin(), faceCoordinates.end() ); + + if( faceCoordinates == targetCoordinates ) + { + neighbors.emplace_back( meshIdxToGlobalId3D.at( candidate.first ) ); + break; + } + } + } + + return neighbors; +} + /** * @brief Build mapping from 2D cells to their neighboring 3D cells using indices * @@ -1028,6 +1167,18 @@ build2DTo3DNeighbors( vtkDataSet & mesh, ArrayOfArrays< localIndex, int64_t > neighbors2Dto3D; neighbors2Dto3D.reserve( cells2DIndices.size() ); + // Build a coordinate lookup for meshes with duplicated/collocated points. + // This is only needed when the input contains 2D cells; the lookup is + // deliberately kept local to this redistribution step. + stdUnorderedMap< VTKPointCoordinate, stdVector< vtkIdType >, VTKPointCoordinateHash > pointsByCoordinate; + pointsByCoordinate.reserve( mesh.GetNumberOfPoints() ); + for( vtkIdType pointId = 0; pointId < mesh.GetNumberOfPoints(); ++pointId ) + { + pointsByCoordinate.get_inserted( getPointCoordinate( mesh, pointId ) ).emplace_back( pointId ); + } + stdUnorderedMap< vtkIdType, stdVector< vtkIdType > > pointTo3DCells; + pointTo3DCells.reserve( mesh.GetNumberOfPoints() ); + // Topology statistics localIndex numStandalone = 0; localIndex numBoundary = 0; // 1 neighbor @@ -1063,6 +1214,24 @@ build2DTo3DNeighbors( vtkDataSet & mesh, // Non-3D neighbors (2D/1D/0D) are silently skipped } + // If point ids were duplicated at a split interface, the exact lookup + // above can miss one or both physical neighbors. Prefer the geometric + // connectivity whenever it is available, and retain the exact result + // as a fallback for meshes whose coordinates are not bitwise identical. + stdVector< int64_t > coordinateNeighborGlobalIds = find2DTo3DNeighborsByCoordinates( + mesh, + pointIds2D, + meshIdxToGlobalId3D, + pointsByCoordinate, + pointTo3DCells ); + if( !coordinateNeighborGlobalIds.empty() ) + { + neighbor3DGlobalIds.resize( coordinateNeighborGlobalIds.size() ); + std::copy( coordinateNeighborGlobalIds.begin(), + coordinateNeighborGlobalIds.end(), + neighbor3DGlobalIds.begin() ); + } + // Update topology statistics localIndex const numNeighbors = neighbor3DGlobalIds.size(); switch( numNeighbors ) From 39f04b1fd2ea781ad98f676ec2ae9221b6b63cd1 Mon Sep 17 00:00:00 2001 From: "Victor A. P. Magri" Date: Thu, 6 Aug 2026 17:18:20 -0400 Subject: [PATCH 15/22] Improve ALM HypreDrive support --- .../interfaces/hypre/hypredrive.cpp | 31 +++++++++++++++++-- ...asePoromechanicsConformingFracturesALM.cpp | 4 +-- 2 files changed, 30 insertions(+), 5 deletions(-) diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp index 1db7a6e581c..67e731d391a 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp @@ -23,6 +23,7 @@ #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseHybridFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanics.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFractures.hpp" +#include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsReservoirFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirFVM.hpp" @@ -40,6 +41,7 @@ #include #include +#include #include #include #include @@ -57,12 +59,26 @@ namespace { void checkHypredriveCall( uint32_t const errorCode, - char const * const call ) + char const * const call, + std::string const & context = {} ) { if( errorCode != 0 ) { HYPREDRV_ErrorCodeDescribe( errorCode ); - GEOS_ERROR( GEOS_FMT( "Error in call to {}", call ) ); + std::fflush( stderr ); + HYPREDRV_ErrorCodeClear(); + + std::string const contextMessage = context.empty() + ? std::string{} + : GEOS_FMT( "\nContext: {}", context ); + GEOS_ERROR( GEOS_FMT( "Error in call to {} (HYPREDRV error code: 0x{:08x}, decimal: {}).{}\n" + "HypreDrive error details are written to stderr; use '2>&1' to merge them into stdout. " + "For additional library tracing, set HYPREDRV_LOG_LEVEL=3 and " + "HYPREDRV_LOG_STREAM=stdout.", + call, + errorCode, + errorCode, + contextMessage ) ); } } @@ -401,6 +417,7 @@ bool strategyUsesCompositionalSemanticLabels( LinearSolverParameters::MGR::Strat case StrategyType::hybridSinglePhasePoromechanics: case StrategyType::singlePhasePoromechanicsEmbeddedFractures: case StrategyType::singlePhasePoromechanicsConformingFractures: + case StrategyType::singlePhasePoromechanicsConformingFracturesALM: case StrategyType::singlePhasePoromechanicsReservoirFVM: case StrategyType::thermalSinglePhasePoromechanicsReservoirFVM: case StrategyType::hydrofracture: @@ -952,6 +969,7 @@ MGRSpecialization getSpecialization( LinearSolverParameters::MGR::StrategyType c } case StrategyType::singlePhasePoromechanicsEmbeddedFractures: case StrategyType::singlePhasePoromechanicsConformingFractures: + case StrategyType::singlePhasePoromechanicsConformingFracturesALM: { MGRSpecialization specialization; specialization.coarseFlavor = AMGFlavor::pressure; @@ -1269,6 +1287,8 @@ bool buildMGRPreconditionerYaml( LinearSolverParameters const & params, return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsEmbeddedFractures >( params, labelNames, numComponentsPerField, preconditionerYaml ); case StrategyType::singlePhasePoromechanicsConformingFractures: return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsConformingFractures >( params, labelNames, numComponentsPerField, preconditionerYaml ); + case StrategyType::singlePhasePoromechanicsConformingFracturesALM: + return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsConformingFracturesALM >( params, labelNames, numComponentsPerField, preconditionerYaml ); case StrategyType::singlePhasePoromechanicsReservoirFVM: return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsReservoirFVM >( params, labelNames, numComponentsPerField, preconditionerYaml ); case StrategyType::thermalSinglePhasePoromechanicsReservoirFVM: @@ -1603,8 +1623,13 @@ void HypredriveSolver::createHypredrive( HypreMatrix const & mat, checkHypredriveCall( HYPREDRV_SetLibraryMode( m_hypredrive ), "HYPREDRV_SetLibraryMode" ); char * argv[] = { const_cast< char * >( parseTarget.argument.c_str() ) }; + std::string const parseContext = + parseTarget.source == hypre::hypredrive::InputSource::authoritativeFile + ? GEOS_FMT( "authoritative YAML file '{}'", parseTarget.argument ) + : "YAML generated by GEOS"; checkHypredriveCall( HYPREDRV_InputArgsParse( 1, argv, m_hypredrive ), - "HYPREDRV_InputArgsParse" ); + "HYPREDRV_InputArgsParse", + parseContext ); if( parseTarget.source == hypre::hypredrive::InputSource::generatedFallback && m_hasExecutionContext && !m_executionContext.solverName.empty() ) diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp index 65da565f2b4..8f55138a5a4 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp @@ -733,7 +733,7 @@ assembleForceResidualDerivativeWrtPressure( string const & meshName, string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); // Use the same kernel launch pattern as SolidMechanicsAugmentedLagrangianContact::assembleForceResidualPressureContribution - this->solidMechanicsSolver()->template forFiniteElementOnFractureSubRegions( meshName, + this->solidMechanicsSolver()->forFiniteElementOnFractureSubRegions( meshName, [&] ( string const &, finiteElement::FiniteElementBase const & subRegionFE, arrayView1d< localIndex const > const & faceElementList ) @@ -829,7 +829,7 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, // Launch the ComputeApertureDerivatives kernel to fill dAperturedU and dAperturedB // This is called for each element type (tri, quad, etc.) - this->solidMechanicsSolver()->template forFiniteElementOnFractureSubRegions( meshName, + this->solidMechanicsSolver()->forFiniteElementOnFractureSubRegions( meshName, [&] ( string const &, finiteElement::FiniteElementBase const & subRegionFE, arrayView1d< localIndex const > const & faceElementList ) From 68fe9d571746c2be5197c71be7d78327e7b1c027 Mon Sep 17 00:00:00 2001 From: "Victor A. P. Magri" Date: Fri, 7 Aug 2026 12:18:21 -0400 Subject: [PATCH 16/22] Improved MGR strategy for SinglePhasePoromechanicsConformingFracturesALM --- .../linearAlgebra/common/LinearSolverBase.hpp | 9 +++ .../interfaces/hypre/hypredrive.cpp | 69 +++++++++++++++++-- .../interfaces/hypre/hypredrive.hpp | 7 ++ ...asePoromechanicsConformingFracturesALM.hpp | 49 +++++++------ .../unitTests/testHypredrive.cpp | 33 +++++++++ .../physicsSolvers/PhysicsSolverBase.cpp | 1 + .../physicsSolvers/PhysicsSolverBase.hpp | 8 +++ .../multiphysics/SinglePhasePoromechanics.cpp | 43 ++++++++++++ .../multiphysics/SinglePhasePoromechanics.hpp | 9 +++ ...asePoromechanicsConformingFracturesALM.cpp | 2 + 10 files changed, 201 insertions(+), 29 deletions(-) diff --git a/src/coreComponents/linearAlgebra/common/LinearSolverBase.hpp b/src/coreComponents/linearAlgebra/common/LinearSolverBase.hpp index 8318c306f3e..d7136035463 100644 --- a/src/coreComponents/linearAlgebra/common/LinearSolverBase.hpp +++ b/src/coreComponents/linearAlgebra/common/LinearSolverBase.hpp @@ -93,6 +93,15 @@ class LinearSolverBase : public PreconditionerBase< LAI > GEOS_UNUSED_VAR( context ); } + /** + * @brief Supply near-null-space vectors used to construct physics-aware preconditioners. + * @param nearNullKernel Distributed vectors spanning the near null space. + */ + virtual void setNearNullKernel( arrayView1d< Vector const > const & nearNullKernel ) + { + GEOS_UNUSED_VAR( nearNullKernel ); + } + /** * @brief @return parameters of the solver. */ diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp index 67e731d391a..4a13b115410 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp @@ -886,7 +886,8 @@ enum class AMGFlavor pressure, pressureTemperature, displacementFiltered, - displacement + displacement, + almDisplacement }; struct LevelAMGBlock @@ -923,6 +924,7 @@ struct MGRSpecialization stdVector< LevelAMGBlock > fRelaxAMGLevels; HYPRE_Int pmax = 0; HYPRE_Int coarseMinCoarseSize = -1; + char const * cycle = nullptr; }; MGRSpecialization getSpecialization( LinearSolverParameters::MGR::StrategyType const strategy ) @@ -969,13 +971,20 @@ MGRSpecialization getSpecialization( LinearSolverParameters::MGR::StrategyType c } case StrategyType::singlePhasePoromechanicsEmbeddedFractures: case StrategyType::singlePhasePoromechanicsConformingFractures: - case StrategyType::singlePhasePoromechanicsConformingFracturesALM: { MGRSpecialization specialization; specialization.coarseFlavor = AMGFlavor::pressure; specialization.fRelaxAMGLevels = { LevelAMGBlock{ 1, AMGFlavor::displacement } }; return specialization; } + case StrategyType::singlePhasePoromechanicsConformingFracturesALM: + { + MGRSpecialization specialization; + specialization.coarseFlavor = AMGFlavor::pressure; + specialization.fRelaxAMGLevels = { LevelAMGBlock{ 1, AMGFlavor::almDisplacement } }; + specialization.cycle = "v(1,0)"; + return specialization; + } case StrategyType::invalid: case StrategyType::singlePhaseReservoirFVM: case StrategyType::singlePhaseHybridFVM: @@ -1020,7 +1029,8 @@ void appendAMGHeader( std::ostringstream & stream, void appendDisplacementAMG( std::ostringstream & stream, integer const indentLevel, integer const separateComponents, - bool const filterFunctions ) + bool const filterFunctions, + bool const useALMSmoother = false ) { appendAMGHeader( stream, indentLevel ); appendLine( stream, indentLevel + 1, "coarsening:" ); @@ -1037,7 +1047,18 @@ void appendDisplacementAMG( std::ostringstream & stream, appendLine( stream, indentLevel + 2, "num_sweeps: 1" ); #else appendLine( stream, indentLevel + 1, "relaxation:" ); - appendLine( stream, indentLevel + 2, "order: 1" ); + if( useALMSmoother ) + { + appendLine( stream, indentLevel + 2, "down_type: l1sym-hgs" ); + appendLine( stream, indentLevel + 2, "up_type: l1sym-hgs" ); + appendLine( stream, indentLevel + 2, "coarse_type: ge" ); + appendLine( stream, indentLevel + 2, "num_sweeps: 1" ); + appendLine( stream, indentLevel + 2, "order: 0" ); + } + else + { + appendLine( stream, indentLevel + 2, "order: 1" ); + } #endif } @@ -1133,6 +1154,11 @@ void appendAMGByFlavor( std::ostringstream & stream, appendDisplacementAMG( stream, indentLevel, 0, false ); break; } + case AMGFlavor::almDisplacement: + { + appendDisplacementAMG( stream, indentLevel, mgrParams.separateComponents, true, true ); + break; + } } } @@ -1182,6 +1208,10 @@ bool buildStrategyYaml( LinearSolverParameters const & params, appendLine( stream, 2, "tolerance: 0.0" ); appendLine( stream, 2, "max_iter: 1" ); appendLine( stream, 2, GEOS_FMT( "print_level: {}", getMGRPrintLevel( params.logLevel ) ) ); + if( specialization.cycle != nullptr ) + { + appendLine( stream, 2, GEOS_FMT( "cycle: {}", specialization.cycle ) ); + } appendLine( stream, 2, GEOS_FMT( "non_c_to_f: {}", 1 ) ); appendLine( stream, 2, GEOS_FMT( "nonglk_max_elmts: {}", 1 ) ); appendLine( stream, 2, GEOS_FMT( "pmax: {}", specialization.pmax ) ); @@ -1600,6 +1630,11 @@ void HypredriveSolver::setExecutionContext( LinearSolverExecutionContext const & m_hasExecutionContext = true; } +void HypredriveSolver::setNearNullKernel( arrayView1d< HypreVector const > const & nearNullKernel ) +{ + m_nearNullKernel = nearNullKernel; +} + void HypredriveSolver::setup( HypreMatrix const & mat ) { Base::setup( mat ); @@ -1675,6 +1710,32 @@ void HypredriveSolver::refreshBoundObjects( HypreMatrix const & mat, pointMarkers.data() ), "HYPREDRV_LinearSystemSetDofmap" ); } + + if( !m_nearNullKernel.empty() ) + { + localIndex const numEntries = mat.numLocalRows(); + localIndex const numModes = m_nearNullKernel.size(); + std::vector< HYPRE_Complex > values( numEntries * numModes ); + + for( localIndex mode = 0; mode < numModes; ++mode ) + { + GEOS_ERROR_IF( m_nearNullKernel[mode].localSize() != numEntries, + "HypreDrive near-null-space vector size does not match the matrix local size" ); + GEOS_LAI_CHECK_ERROR( + HYPRE_IJVectorGetValues( m_nearNullKernel[mode].unwrappedIJ(), + LvArray::integerConversion< HYPRE_Int >( numEntries ), + nullptr, + values.data() + mode * numEntries ) ); + } + + checkHypredriveCall( + HYPREDRV_LinearSystemSetNearNullSpace( + m_hypredrive, + LvArray::integerConversion< int >( numEntries ), + LvArray::integerConversion< int >( numModes ), + values.data() ), + "HYPREDRV_LinearSystemSetNearNullSpace" ); + } } void HypredriveSolver::setupLegacy( HypreMatrix const & mat ) diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.hpp index 5d70d0f124e..ab7adef728b 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.hpp @@ -98,6 +98,12 @@ class HypredriveSolver final : public LinearSolverBase< HypreInterface > */ void setExecutionContext( LinearSolverExecutionContext const & context ) override; + /** + * @brief Set near-null-space modes to pass to HypreDrive during setup. + * @param nearNullKernel Full-system distributed near-null-space vectors. + */ + void setNearNullKernel( arrayView1d< HypreVector const > const & nearNullKernel ) override; + /** * @brief Build or refresh the solver/preconditioner for a matrix. * @param mat Matrix that defines the system structure and coefficients. @@ -164,6 +170,7 @@ class HypredriveSolver final : public LinearSolverBase< HypreInterface > bool m_hasExecutionContext = false; bool m_timestepScopeActive = false; bool m_newtonScopeActive = false; + arrayView1d< HypreVector const > m_nearNullKernel; size_t m_hypredriveGeneration = 0; HYPREDRV_t m_hypredrive{}; mutable HypreVector m_dummyRhs; diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp index e9fb7db8bd6..8799a8277a3 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -14,7 +14,7 @@ */ /** - * @file SinglePhasePoromechanicsConformingFractures.hpp + * @file SinglePhasePoromechanicsConformingFracturesALM.hpp */ #ifndef GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_ @@ -32,7 +32,7 @@ namespace mgr { /** - * @brief SinglePhasePoromechanicsConformingFractures strategy. + * @brief SinglePhasePoromechanicsConformingFracturesALM strategy. * * dofLabel: 0 = displacement, x-component * dofLabel: 1 = displacement, y-component @@ -44,11 +44,12 @@ namespace mgr * * Ingredients: - * 1. Level 1: F-points displacement (3,4,5), C-points pressure (0,1,2,6) - * 2. Level 2: F-points displacement (0,1,2), C-points pressure (6) - * 2. F-points smoother: BoomerAMG, single V-cycle - * 3. C-points coarse-grid/Schur complement solver: BoomerAMG - * 4. Global smoother: none + * 1. Level 0 eliminates bubble displacement (3,4,5) with L1-Jacobi. + * 2. Level 1 eliminates nodal displacement (0,1,2) with one BoomerAMG V-cycle. + * 3. The displacement AMG uses three functions, separate-component filtering, + * and symmetric L1 hybrid Gauss-Seidel relaxation on CPU. + * 4. The pressure Schur complement is solved with BoomerAMG. + * 5. The MGR V-cycle uses pre-relaxation only and no global smoother. */ class SinglePhasePoromechanicsConformingFracturesALM : public MGRStrategyBase< 2 > { @@ -61,13 +62,9 @@ class SinglePhasePoromechanicsConformingFracturesALM : public MGRStrategyBase< 2 : MGRStrategyBase( 7 ) { - // we keep u and p - m_labels[0].push_back( 0 ); - m_labels[0].push_back( 1 ); - m_labels[0].push_back( 2 ); - m_labels[0].push_back( 6 ); - // we keep p - m_labels[1].push_back( 6 ); + // Eliminate bubble displacement first, then nodal displacement. + m_labels[0] = { 0, 1, 2, 6 }; + m_labels[1] = { 6 }; setupLabels(); @@ -100,23 +97,25 @@ class SinglePhasePoromechanicsConformingFracturesALM : public MGRStrategyBase< 2 HypreMGRData & mgrData ) { setReduction( precond, mgrData ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetCycleType( precond.ptr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetFRelaxCycle( precond.ptr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetGlobalSmoothCycle( precond.ptr, 1 ) ); - // Configure the BoomerAMG solver used as F-relaxation for the second level - // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGCreate( &mgrData.mechSolver.ptr ) ); - // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetTol( mgrData.mechSolver.ptr, 0.0 ) ); - // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxIter( mgrData.mechSolver.ptr, 1 ) ); - // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxRowSum( mgrData.mechSolver.ptr, 1.0 ) ); - // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetStrongThreshold( mgrData.mechSolver.ptr, 0.6 ) ); - // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPrintLevel( mgrData.mechSolver.ptr, 0 ) ); - // GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumFunctions( mgrData.mechSolver.ptr, 3 ) ); - setDisplacementAMG(mgrData.mechSolver, mgrParams.separateComponents); + // Configure the BoomerAMG solver used as level-1 F-relaxation. + setDisplacementAMG( mgrData.mechSolver, mgrParams.separateComponents ); #if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCoarsenType( mgrData.mechSolver.ptr, hypre::getAMGCoarseningType( LinearSolverParameters::AMG::CoarseningType::PMIS ) ) ); GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxType( mgrData.mechSolver.ptr, hypre::getAMGRelaxationType( LinearSolverParameters::AMG::SmootherType::chebyshev ) ) ); GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumSweeps( mgrData.mechSolver.ptr, 1 ) ); #else - GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxOrder( mgrData.mechSolver.ptr, 1 ) ); + HYPRE_Int constexpr l1SymmetricHybridGaussSeidel = 89; + HYPRE_Int constexpr gaussianElimination = 9; + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCycleRelaxType( mgrData.mechSolver.ptr, l1SymmetricHybridGaussSeidel, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCycleRelaxType( mgrData.mechSolver.ptr, l1SymmetricHybridGaussSeidel, 2 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCycleRelaxType( mgrData.mechSolver.ptr, gaussianElimination, 3 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumSweeps( mgrData.mechSolver.ptr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxOrder( mgrData.mechSolver.ptr, 0 ) ); #endif GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetFSolverAtLevel( precond.ptr, mgrData.mechSolver.ptr, 1 ) ); @@ -131,4 +130,4 @@ class SinglePhasePoromechanicsConformingFracturesALM : public MGRStrategyBase< 2 } // namespace geos -#endif /*GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_*/ \ No newline at end of file +#endif /*GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_*/ diff --git a/src/coreComponents/linearAlgebra/unitTests/testHypredrive.cpp b/src/coreComponents/linearAlgebra/unitTests/testHypredrive.cpp index c004631a83d..d157a821ffc 100644 --- a/src/coreComponents/linearAlgebra/unitTests/testHypredrive.cpp +++ b/src/coreComponents/linearAlgebra/unitTests/testHypredrive.cpp @@ -261,6 +261,39 @@ TEST( HypredriveYaml, BuildsGeneratedYamlForEveryMGRStrategy ) } } +TEST( HypredriveYaml, BuildsSelectedALMPoromechanicsMGRStrategy ) +{ + stdVector< string > const fieldNames = { "totalDisplacement", "totalBubbleDisplacement", "pressure" }; + array1d< int > numComponentsPerField( 3 ); + numComponentsPerField[0] = 3; + numComponentsPerField[1] = 3; + numComponentsPerField[2] = 1; + + hypre::hypredrive::InputArgsParseTarget target; + ASSERT_TRUE( hypre::hypredrive::buildInputArgsParseTarget( + makeMgrParameters( LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsConformingFracturesALM ), + fieldNames, + numComponentsPerField, + target ) ); + + EXPECT_NE( target.argument.find( "num_levels: 3" ), std::string::npos ); + EXPECT_NE( target.argument.find( "cycle: v(1,0)" ), std::string::npos ); + EXPECT_NE( target.argument.find( "f_dofs: [totalBubbleDisplacement_0, totalBubbleDisplacement_1, totalBubbleDisplacement_2]" ), + std::string::npos ); + EXPECT_NE( target.argument.find( "f_dofs: [totalDisplacement_0, totalDisplacement_1, totalDisplacement_2]" ), + std::string::npos ); + EXPECT_NE( target.argument.find( "filter_functions: 1" ), std::string::npos ); +#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP + EXPECT_NE( target.argument.find( "down_type: 16" ), std::string::npos ); + EXPECT_NE( target.argument.find( "up_type: 16" ), std::string::npos ); +#else + EXPECT_NE( target.argument.find( "down_type: l1sym-hgs" ), std::string::npos ); + EXPECT_NE( target.argument.find( "up_type: l1sym-hgs" ), std::string::npos ); + EXPECT_NE( target.argument.find( "coarse_type: ge" ), std::string::npos ); + EXPECT_NE( target.argument.find( "order: 0" ), std::string::npos ); +#endif +} + TEST( HypredriveYaml, UsesSemanticCompositionalLabelsWhenFieldNamesAreAvailable ) { stdVector< string > const fieldNames = { "totalDisplacement", "compositionalVariables" }; diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp index 8d05dfd6958..f5dc5cdb3b1 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp @@ -1444,6 +1444,7 @@ void PhysicsSolverBase::solveLinearSystem( DofManager const & dofManager, executionContext.nonlinearIteration = nonlinearIteration; executionContext.systemSetupTimestamp = getSystemSetupTimestamp(); m_linearSolver->setExecutionContext( executionContext ); + m_linearSolver->setNearNullKernel( getLinearSolverNearNullKernel() ); if( isSetupNeeded ) { diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp index a174b485e67..107bea65a88 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp @@ -1016,6 +1016,14 @@ class PhysicsSolverBase : public ExecutableGroup virtual void postInputInitialization() override; + /** + * @brief Return optional near-null-space modes for the monolithic linear solver. + */ + virtual arrayView1d< ParallelVector const > getLinearSolverNearNullKernel() const + { + return {}; + } + /** * @brief Eisenstat-Walker adaptive tolerance * diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp index 468117b09dd..06591ccd0d4 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp @@ -25,6 +25,10 @@ #include "linearAlgebra/multiscale/MultiscalePreconditioner.hpp" #include "linearAlgebra/solvers/BlockPreconditioner.hpp" #include "linearAlgebra/solvers/SeparateComponentPreconditioner.hpp" +#include "linearAlgebra/utilities/LAIHelperFunctions.hpp" +#ifdef GEOS_USE_HYPREDRV +#include "linearAlgebra/interfaces/hypre/hypredrive.hpp" +#endif #include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsDamage.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp" @@ -78,6 +82,7 @@ void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::setupSystem( Dom // setup monolithic coupled system PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity ); + setupLinearSolverNearNullKernel( domain, dofManager ); if( !this->m_precond && this->m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) { @@ -85,6 +90,44 @@ void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::setupSystem( Dom } } +template< typename FLOW_SOLVER, typename MECHANICS_SOLVER > +void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::setupLinearSolverNearNullKernel( + DomainPartition & domain, + DofManager const & dofManager ) +{ + m_linearSolverNearNullKernel.resize( 0 ); + +#ifdef GEOS_USE_HYPREDRV + if( !hypre::hypredrive::shouldUse( this->m_linearSolverParameters.get() ) ) + { + return; + } + + this->solidMechanicsSolver()->forDiscretizationOnMeshTargets( + domain.getMeshBodies(), + [&]( string const &, MeshLevel const & mesh, string_array const & ) + { + if( !m_linearSolverNearNullKernel.empty() ) + { + return; + } + + NodeManager const & nodeManager = mesh.getNodeManager(); + arrayView1d< globalIndex const > const dofIndex = + nodeManager.getReference< array1d< globalIndex > >( + dofManager.getKey( solidMechanics::totalDisplacement::key() ) ); + m_linearSolverNearNullKernel = + LAIHelperFunctions::computeRigidBodyModes< ParallelVector >( + nodeManager.referencePosition(), + dofIndex, + dofManager.rankOffset(), + dofManager.numLocalDofs() ); + } ); +#else + GEOS_UNUSED_VAR( domain, dofManager ); +#endif +} + template< typename FLOW_SOLVER, typename MECHANICS_SOLVER > void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::initializePostInitialConditionsPreSubGroups() { diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp index d7e40171f79..518c7cac693 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp @@ -116,6 +116,14 @@ class SinglePhasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHAN virtual void initializePostInitialConditionsPreSubGroups() override; + void setupLinearSolverNearNullKernel( DomainPartition & domain, + DofManager const & dofManager ); + + arrayView1d< ParallelVector const > getLinearSolverNearNullKernel() const override + { + return m_linearSolverNearNullKernel.toViewConst(); + } + virtual void setMGRStrategy() override { if( this->m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::mgr ) @@ -155,6 +163,7 @@ class SinglePhasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHAN virtual string getFlowDofKey() const override { return SinglePhaseBase::viewKeyStruct::elemDofFieldString(); } + array1d< ParallelVector > m_linearSolverNearNullKernel; integer m_damageFlag; }; diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp index 8f55138a5a4..159f906c0de 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp @@ -147,6 +147,8 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupSystem( solution.setName( this->getName() + "/solution" ); solution.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); + this->setupLinearSolverNearNullKernel( domain, dofManager ); + if( !this->m_precond && this->m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) { this->m_precond = this->createPreconditioner( domain ); From bf486974ddf7f8232274d841f32eb7ff7c779aa6 Mon Sep 17 00:00:00 2001 From: "Victor A. P. Magri" Date: Fri, 7 Aug 2026 18:41:36 -0700 Subject: [PATCH 17/22] Fix compilation error --- src/coreComponents/common/format/LogPart.hpp | 2 -- src/coreComponents/common/format/table/TableData.hpp | 2 -- 2 files changed, 4 deletions(-) diff --git a/src/coreComponents/common/format/LogPart.hpp b/src/coreComponents/common/format/LogPart.hpp index f69729b25bc..68a83bc22b8 100644 --- a/src/coreComponents/common/format/LogPart.hpp +++ b/src/coreComponents/common/format/LogPart.hpp @@ -197,8 +197,6 @@ void LogPart::addDescriptionBySection( Description & description, FormattedDescr size_t & formattedDescriptionMaxWidth = formattedDescription.m_maxValueWidth; size_t & formattedDescriptionNameWidth = formattedDescription.m_maxNameWidth; ( [&] { - static_assert( has_formatter_v< decltype(args) >, - "Argument passed cannot be converted to string" ); string const value = GEOS_FMT( "{}", args ); stdVector< string_view > dividedDescriptionValues = diff --git a/src/coreComponents/common/format/table/TableData.hpp b/src/coreComponents/common/format/table/TableData.hpp index e36afb037c5..bfadcacd8d0 100644 --- a/src/coreComponents/common/format/table/TableData.hpp +++ b/src/coreComponents/common/format/table/TableData.hpp @@ -267,7 +267,6 @@ void TableData::addRow( Args const &... args ) { stdVector< CellData > cells; ( [&] { - static_assert( has_formatter_v< decltype(args) > || isCellType< std::decay_t< decltype(args) > >, "Argument passed in addRow cannot be converted to string nor a CellType" ); if constexpr (std::is_same_v< Args, CellType >) { cells.push_back( { args, string() } ); } @@ -297,7 +296,6 @@ void TableData::addRow( Args const &... args ) template< typename T > void TableData2D::addCell( real64 const rowValue, real64 const columnValue, T const & value ) { - static_assert( has_formatter_v< decltype(value) >, "Argument passed in addCell cannot be converted to string" ); m_columnValues.insert( columnValue ); m_data.get_inserted( rowValue ).get_inserted( columnValue ) = GEOS_FMT( "{}", value ); } From 4c29442b66a58b4c76b2d2ffd52fdf96978601a3 Mon Sep 17 00:00:00 2001 From: "Victor A. P. Magri" Date: Fri, 7 Aug 2026 18:41:57 -0700 Subject: [PATCH 18/22] Add missing MGR strategy bridge --- .../linearAlgebra/interfaces/hypre/hypredrive.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp index 1db7a6e581c..04afb7d4aef 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp @@ -401,6 +401,7 @@ bool strategyUsesCompositionalSemanticLabels( LinearSolverParameters::MGR::Strat case StrategyType::hybridSinglePhasePoromechanics: case StrategyType::singlePhasePoromechanicsEmbeddedFractures: case StrategyType::singlePhasePoromechanicsConformingFractures: + case StrategyType::singlePhasePoromechanicsConformingFracturesALM: case StrategyType::singlePhasePoromechanicsReservoirFVM: case StrategyType::thermalSinglePhasePoromechanicsReservoirFVM: case StrategyType::hydrofracture: @@ -952,6 +953,7 @@ MGRSpecialization getSpecialization( LinearSolverParameters::MGR::StrategyType c } case StrategyType::singlePhasePoromechanicsEmbeddedFractures: case StrategyType::singlePhasePoromechanicsConformingFractures: + case StrategyType::singlePhasePoromechanicsConformingFracturesALM: { MGRSpecialization specialization; specialization.coarseFlavor = AMGFlavor::pressure; @@ -1269,6 +1271,8 @@ bool buildMGRPreconditionerYaml( LinearSolverParameters const & params, return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsEmbeddedFractures >( params, labelNames, numComponentsPerField, preconditionerYaml ); case StrategyType::singlePhasePoromechanicsConformingFractures: return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsConformingFractures >( params, labelNames, numComponentsPerField, preconditionerYaml ); + case StrategyType::singlePhasePoromechanicsConformingFracturesALM: + return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsConformingFractures >( params, labelNames, numComponentsPerField, preconditionerYaml ); case StrategyType::singlePhasePoromechanicsReservoirFVM: return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsReservoirFVM >( params, labelNames, numComponentsPerField, preconditionerYaml ); case StrategyType::thermalSinglePhasePoromechanicsReservoirFVM: From 35205fac1584e135c9e1df4da2df85b24a915521 Mon Sep 17 00:00:00 2001 From: "Victor A. P. Magri" Date: Fri, 7 Aug 2026 18:42:15 -0700 Subject: [PATCH 19/22] Add compatibility with vtk 9.6.2 --- .../fileIO/vtk/VTKPolyDataWriterInterface.cpp | 7 +++++++ src/coreComponents/mesh/generators/VTKUtilities.cpp | 6 ++++++ 2 files changed, 13 insertions(+) diff --git a/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp b/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp index bb3c0182b12..cb199d66e11 100644 --- a/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp +++ b/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp @@ -32,6 +32,7 @@ #include #include #include +#include #include #include // System includes @@ -308,7 +309,9 @@ getWell( WellElementSubRegion const & subRegion, localIndex const numPoints = subRegion.size() > 0 ? subRegion.size() + 1 : 0; points->SetNumberOfPoints( numPoints ); auto cellsArray = vtkSmartPointer< vtkCellArray >::New(); +#if VTK_VERSION_NUMBER < VTK_VERSION_CHECK( 9, 6, 0 ) cellsArray->SetNumberOfCells( subRegion.size() ); +#endif localIndex const numberOfNodesPerElement = subRegion.numNodesPerElement(); GEOS_ERROR_IF_NE( numberOfNodesPerElement, 2 ); stdVector< vtkIdType > connectivity( numberOfNodesPerElement ); @@ -358,7 +361,9 @@ getSurface( FaceElementSubRegion const & subRegion, auto & faceToNodes = faceManager.nodeList(); auto cellArray = vtkSmartPointer< vtkCellArray >::New(); +#if VTK_VERSION_NUMBER < VTK_VERSION_CHECK( 9, 6, 0 ) cellArray->SetNumberOfCells( subRegion.size() ); +#endif stdVector< int > cellTypes; cellTypes.reserve( subRegion.size() ); @@ -599,7 +604,9 @@ static ParticleData getVtkCells( ParticleRegion const & region ) { vtkSmartPointer< vtkCellArray > cellsArray = vtkCellArray::New(); +#if VTK_VERSION_NUMBER < VTK_VERSION_CHECK( 9, 6, 0 ) cellsArray->SetNumberOfCells( region.getNumberOfParticles< ParticleRegion >() ); +#endif stdVector< int > cellType; cellType.reserve( region.getNumberOfParticles< ParticleRegion >() ); diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp index 1dca5f69d5f..3c59612ddd0 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp @@ -34,6 +34,7 @@ #include #include #include +#include #include #include #include @@ -68,6 +69,7 @@ #include #include #include +#include #ifdef GEOS_USE_MPI #include @@ -980,7 +982,11 @@ static void classifyCellsByDimension( vtkDataSet & mesh, // Single pass: classify and populate for( vtkIdType i = 0; i < numCells; ++i ) { +#if VTK_VERSION_NUMBER >= VTK_VERSION_CHECK( 9, 6, 0 ) + int const dim = vtkCellTypeUtilities::GetDimension( mesh.GetCellType( i ) ); +#else int const dim = vtkCellTypes::GetDimension( mesh.GetCellType( i ) ); +#endif if( dim == 3 ) { cells3DIndices.emplace_back( i ); From 8fe08f2367c790f59a10ae816b6a75457aa5a025 Mon Sep 17 00:00:00 2001 From: "Victor A. P. Magri" Date: Fri, 7 Aug 2026 23:23:39 -0400 Subject: [PATCH 20/22] Fix bug in the fully coupled ALM + single-phase sparsity construction MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - The mechanics solver replaced the existing flow sparsity pattern, discarding flow and custom FIM couplings. - Pressure–bubble couplings checked for the nonexistent wrapper bubbleElementsList instead of the registered bubbleCells key. Fixed both by unioning the flow and mechanics patterns before adding FIM couplings, and by using the correct wrapper key --- ...asePoromechanicsConformingFracturesALM.cpp | 45 +++++++++++-------- 1 file changed, 26 insertions(+), 19 deletions(-) diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp index 159f906c0de..349f01e7af7 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp @@ -95,15 +95,20 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupSystem( if( setSparsity ) { - // Get sparsity pattern from flow solver (could include wells) - SparsityPattern< globalIndex > patternOriginal; - this->flowSolver()->setSparsityPattern( domain, dofManager, localMatrix, patternOriginal ); - - // Get the original row lengths - array1d< localIndex > rowLengths( patternOriginal.numRows() ); - for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) + // Start from both subsolver patterns. The flow pattern may contain well and + // flux couplings, while the mechanics pattern contains the nodal-bubble couplings. + SparsityPattern< globalIndex > flowPattern; + this->flowSolver()->setSparsityPattern( domain, dofManager, localMatrix, flowPattern ); + + SparsityPattern< globalIndex > mechanicsPattern; + this->solidMechanicsSolver()->setSparsityPattern( domain, dofManager, localMatrix, mechanicsPattern ); + + // Sum the row lengths to provide enough capacity for the union. Entries + // common to both patterns are deduplicated when they are inserted below. + array1d< localIndex > rowLengths( flowPattern.numRows() ); + for( localIndex localRow = 0; localRow < flowPattern.numRows(); ++localRow ) { - rowLengths[localRow] = patternOriginal.numNonZeros( localRow ); + rowLengths[localRow] = flowPattern.numNonZeros( localRow ) + mechanicsPattern.numNonZeros( localRow ); } // Add the number of nonzeros induced by coupling @@ -113,15 +118,20 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupSystem( // Create a new pattern with enough capacity for coupled matrix SparsityPattern< globalIndex > pattern; - pattern.resizeFromRowCapacities< parallelHostPolicy >( patternOriginal.numRows(), - patternOriginal.numColumns(), + pattern.resizeFromRowCapacities< parallelHostPolicy >( flowPattern.numRows(), + flowPattern.numColumns(), rowLengths.data() ); - // Copy the original nonzeros - for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) + // Copy the flow and mechanics nonzeros. + for( localIndex localRow = 0; localRow < flowPattern.numRows(); ++localRow ) + { + globalIndex const * cols = flowPattern.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + flowPattern.numNonZeros( localRow ) ); + } + for( localIndex localRow = 0; localRow < mechanicsPattern.numRows(); ++localRow ) { - globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); - pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow ) ); + globalIndex const * cols = mechanicsPattern.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + mechanicsPattern.numNonZeros( localRow ) ); } // Add the nonzeros from coupling @@ -129,9 +139,6 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupSystem( addPressureForceCouplingPattern( domain, dofManager, pattern.toView() ); addMatrixPressureBubbleCouplingPattern( domain, dofManager, pattern.toView() ); - // Assemble the full sparsity pattern using the solid mechanics solver - this->solidMechanicsSolver()->setSparsityPattern( domain, dofManager, localMatrix, pattern ); - // Set up the derivative flux residual matrix setUpDflux_dApertureMatrix( domain, dofManager, localMatrix ); @@ -1089,7 +1096,7 @@ addMatrixPressureBubbleCouplingNNZ( DomainPartition const & domain, elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, CellElementSubRegion const & subRegion ) { - if( !subRegion.hasWrapper( "bubbleElementsList" ) ) + if( !subRegion.hasWrapper( CellElementSubRegion::viewKeyStruct::bubbleCellsString() ) ) return; arrayView1d< localIndex const > const bubbleElems = subRegion.bubbleElementsList(); @@ -1141,7 +1148,7 @@ addMatrixPressureBubbleCouplingPattern( DomainPartition const & domain, elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, CellElementSubRegion const & subRegion ) { - if( !subRegion.hasWrapper( "bubbleElementsList" ) ) + if( !subRegion.hasWrapper( CellElementSubRegion::viewKeyStruct::bubbleCellsString() ) ) return; arrayView1d< localIndex const > const bubbleElems = subRegion.bubbleElementsList(); From 1cba5630be57ab95647551f5e7a5438252428104 Mon Sep 17 00:00:00 2001 From: "Victor A. P. Magri" Date: Sat, 8 Aug 2026 14:19:07 -0400 Subject: [PATCH 21/22] Fix SIGFPE issue --- src/coreComponents/mesh/generators/VTKUtilities.cpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp index cd29e657a75..236efbb73c6 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp @@ -1971,7 +1971,14 @@ redistributeByKdTree( vtkDataSet & mesh ) vtkNew< vtkRedistributeDataSetFilter > rdsf; rdsf->SetInputDataObject( &mesh ); rdsf->SetNumberOfPartitions( MpiWrapper::commSize() ); - rdsf->Update(); + { + // vtkRedistributeDataSetFilter uses VTK's XML writer internally to + // serialize datasets exchanged by DIY. The writer may raise floating-point + // exceptions while calculating progress for empty arrays. These exceptions + // are harmless to VTK, but GEOS' enabled FPE traps turn them into SIGFPE. + LvArray::system::FloatingPointExceptionGuard guard; + rdsf->Update(); + } vtkSmartPointer< vtkDataSet > result = vtkDataSet::SafeDownCast( rdsf->GetOutputDataObject( 0 ) ); From 43b7c06f302126620ae32ae78dfa6a477cbc9459 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Thu, 13 Aug 2026 14:01:41 +0200 Subject: [PATCH 22/22] missing guards --- .../linearAlgebra/interfaces/hypre/hypredrive.cpp | 7 ------- src/coreComponents/mesh/generators/VTKUtilities.cpp | 2 ++ 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp index 706d70003fc..875f177d088 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp @@ -972,13 +972,6 @@ MGRSpecialization getSpecialization( LinearSolverParameters::MGR::StrategyType c case StrategyType::singlePhasePoromechanicsEmbeddedFractures: case StrategyType::singlePhasePoromechanicsConformingFractures: case StrategyType::singlePhasePoromechanicsConformingFracturesALM: - { - MGRSpecialization specialization; - specialization.coarseFlavor = AMGFlavor::pressure; - specialization.fRelaxAMGLevels = { LevelAMGBlock{ 1, AMGFlavor::displacement } }; - return specialization; - } - case StrategyType::singlePhasePoromechanicsConformingFracturesALM: { MGRSpecialization specialization; specialization.coarseFlavor = AMGFlavor::pressure; diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp index b6705582caa..d8e66284329 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp @@ -34,7 +34,9 @@ #include #include #include +#if VTK_VERSION_NUMBER >= VTK_VERSION_CHECK( 9, 6, 0 ) #include +#endif #include #include #include