From bd68d98536a85e397ea53e0f197d8b768a62f743 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Mon, 2 Feb 2026 13:59:19 +0100 Subject: [PATCH 1/8] 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 2/8] 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 3/8] 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 4/8] 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 5/8] 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 6/8] 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 7/8] 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 b87790363d9025ec71c9c1d50213b7c2907f05ff Mon Sep 17 00:00:00 2001 From: npillardou Date: Thu, 12 Feb 2026 13:25:35 +0100 Subject: [PATCH 8/8] 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 ) ) );