diff --git a/src/coreComponents/common/format/LogPart.hpp b/src/coreComponents/common/format/LogPart.hpp index f69729b25bc..68a83bc22b8 100644 --- a/src/coreComponents/common/format/LogPart.hpp +++ b/src/coreComponents/common/format/LogPart.hpp @@ -197,8 +197,6 @@ void LogPart::addDescriptionBySection( Description & description, FormattedDescr size_t & formattedDescriptionMaxWidth = formattedDescription.m_maxValueWidth; size_t & formattedDescriptionNameWidth = formattedDescription.m_maxNameWidth; ( [&] { - static_assert( has_formatter_v< decltype(args) >, - "Argument passed cannot be converted to string" ); string const value = GEOS_FMT( "{}", args ); stdVector< string_view > dividedDescriptionValues = diff --git a/src/coreComponents/common/format/table/TableData.hpp b/src/coreComponents/common/format/table/TableData.hpp index e36afb037c5..bfadcacd8d0 100644 --- a/src/coreComponents/common/format/table/TableData.hpp +++ b/src/coreComponents/common/format/table/TableData.hpp @@ -267,7 +267,6 @@ void TableData::addRow( Args const &... args ) { stdVector< CellData > cells; ( [&] { - static_assert( has_formatter_v< decltype(args) > || isCellType< std::decay_t< decltype(args) > >, "Argument passed in addRow cannot be converted to string nor a CellType" ); if constexpr (std::is_same_v< Args, CellType >) { cells.push_back( { args, string() } ); } @@ -297,7 +296,6 @@ void TableData::addRow( Args const &... args ) template< typename T > void TableData2D::addCell( real64 const rowValue, real64 const columnValue, T const & value ) { - static_assert( has_formatter_v< decltype(value) >, "Argument passed in addCell cannot be converted to string" ); m_columnValues.insert( columnValue ); m_data.get_inserted( rowValue ).get_inserted( columnValue ) = GEOS_FMT( "{}", value ); } diff --git a/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp b/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp index bb3c0182b12..cb199d66e11 100644 --- a/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp +++ b/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp @@ -32,6 +32,7 @@ #include #include #include +#include #include #include // System includes @@ -308,7 +309,9 @@ getWell( WellElementSubRegion const & subRegion, localIndex const numPoints = subRegion.size() > 0 ? subRegion.size() + 1 : 0; points->SetNumberOfPoints( numPoints ); auto cellsArray = vtkSmartPointer< vtkCellArray >::New(); +#if VTK_VERSION_NUMBER < VTK_VERSION_CHECK( 9, 6, 0 ) cellsArray->SetNumberOfCells( subRegion.size() ); +#endif localIndex const numberOfNodesPerElement = subRegion.numNodesPerElement(); GEOS_ERROR_IF_NE( numberOfNodesPerElement, 2 ); stdVector< vtkIdType > connectivity( numberOfNodesPerElement ); @@ -358,7 +361,9 @@ getSurface( FaceElementSubRegion const & subRegion, auto & faceToNodes = faceManager.nodeList(); auto cellArray = vtkSmartPointer< vtkCellArray >::New(); +#if VTK_VERSION_NUMBER < VTK_VERSION_CHECK( 9, 6, 0 ) cellArray->SetNumberOfCells( subRegion.size() ); +#endif stdVector< int > cellTypes; cellTypes.reserve( subRegion.size() ); @@ -599,7 +604,9 @@ static ParticleData getVtkCells( ParticleRegion const & region ) { vtkSmartPointer< vtkCellArray > cellsArray = vtkCellArray::New(); +#if VTK_VERSION_NUMBER < VTK_VERSION_CHECK( 9, 6, 0 ) cellsArray->SetNumberOfCells( region.getNumberOfParticles< ParticleRegion >() ); +#endif stdVector< int > cellType; cellType.reserve( region.getNumberOfParticles< ParticleRegion >() ); diff --git a/src/coreComponents/linearAlgebra/CMakeLists.txt b/src/coreComponents/linearAlgebra/CMakeLists.txt index af5310e79ea..e4919321958 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/common/LinearSolverBase.hpp b/src/coreComponents/linearAlgebra/common/LinearSolverBase.hpp index 8318c306f3e..d7136035463 100644 --- a/src/coreComponents/linearAlgebra/common/LinearSolverBase.hpp +++ b/src/coreComponents/linearAlgebra/common/LinearSolverBase.hpp @@ -93,6 +93,15 @@ class LinearSolverBase : public PreconditionerBase< LAI > GEOS_UNUSED_VAR( context ); } + /** + * @brief Supply near-null-space vectors used to construct physics-aware preconditioners. + * @param nearNullKernel Distributed vectors spanning the near null space. + */ + virtual void setNearNullKernel( arrayView1d< Vector const > const & nearNullKernel ) + { + GEOS_UNUSED_VAR( nearNullKernel ); + } + /** * @brief @return parameters of the solver. */ diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp index 84912abe16c..defd7e07818 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" @@ -188,6 +189,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/hypredrive.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp index 1db7a6e581c..875f177d088 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp @@ -23,6 +23,7 @@ #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseHybridFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanics.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFractures.hpp" +#include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsReservoirFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirFVM.hpp" @@ -40,6 +41,7 @@ #include #include +#include #include #include #include @@ -57,12 +59,26 @@ namespace { void checkHypredriveCall( uint32_t const errorCode, - char const * const call ) + char const * const call, + std::string const & context = {} ) { if( errorCode != 0 ) { HYPREDRV_ErrorCodeDescribe( errorCode ); - GEOS_ERROR( GEOS_FMT( "Error in call to {}", call ) ); + std::fflush( stderr ); + HYPREDRV_ErrorCodeClear(); + + std::string const contextMessage = context.empty() + ? std::string{} + : GEOS_FMT( "\nContext: {}", context ); + GEOS_ERROR( GEOS_FMT( "Error in call to {} (HYPREDRV error code: 0x{:08x}, decimal: {}).{}\n" + "HypreDrive error details are written to stderr; use '2>&1' to merge them into stdout. " + "For additional library tracing, set HYPREDRV_LOG_LEVEL=3 and " + "HYPREDRV_LOG_STREAM=stdout.", + call, + errorCode, + errorCode, + contextMessage ) ); } } @@ -401,6 +417,7 @@ bool strategyUsesCompositionalSemanticLabels( LinearSolverParameters::MGR::Strat case StrategyType::hybridSinglePhasePoromechanics: case StrategyType::singlePhasePoromechanicsEmbeddedFractures: case StrategyType::singlePhasePoromechanicsConformingFractures: + case StrategyType::singlePhasePoromechanicsConformingFracturesALM: case StrategyType::singlePhasePoromechanicsReservoirFVM: case StrategyType::thermalSinglePhasePoromechanicsReservoirFVM: case StrategyType::hydrofracture: @@ -869,7 +886,8 @@ enum class AMGFlavor pressure, pressureTemperature, displacementFiltered, - displacement + displacement, + almDisplacement }; struct LevelAMGBlock @@ -906,6 +924,7 @@ struct MGRSpecialization stdVector< LevelAMGBlock > fRelaxAMGLevels; HYPRE_Int pmax = 0; HYPRE_Int coarseMinCoarseSize = -1; + char const * cycle = nullptr; }; MGRSpecialization getSpecialization( LinearSolverParameters::MGR::StrategyType const strategy ) @@ -952,10 +971,12 @@ MGRSpecialization getSpecialization( LinearSolverParameters::MGR::StrategyType c } case StrategyType::singlePhasePoromechanicsEmbeddedFractures: case StrategyType::singlePhasePoromechanicsConformingFractures: + case StrategyType::singlePhasePoromechanicsConformingFracturesALM: { MGRSpecialization specialization; specialization.coarseFlavor = AMGFlavor::pressure; - specialization.fRelaxAMGLevels = { LevelAMGBlock{ 1, AMGFlavor::displacement } }; + specialization.fRelaxAMGLevels = { LevelAMGBlock{ 1, AMGFlavor::almDisplacement } }; + specialization.cycle = "v(1,0)"; return specialization; } case StrategyType::invalid: @@ -1002,7 +1023,8 @@ void appendAMGHeader( std::ostringstream & stream, void appendDisplacementAMG( std::ostringstream & stream, integer const indentLevel, integer const separateComponents, - bool const filterFunctions ) + bool const filterFunctions, + bool const useALMSmoother = false ) { appendAMGHeader( stream, indentLevel ); appendLine( stream, indentLevel + 1, "coarsening:" ); @@ -1019,7 +1041,18 @@ void appendDisplacementAMG( std::ostringstream & stream, appendLine( stream, indentLevel + 2, "num_sweeps: 1" ); #else appendLine( stream, indentLevel + 1, "relaxation:" ); - appendLine( stream, indentLevel + 2, "order: 1" ); + if( useALMSmoother ) + { + appendLine( stream, indentLevel + 2, "down_type: l1sym-hgs" ); + appendLine( stream, indentLevel + 2, "up_type: l1sym-hgs" ); + appendLine( stream, indentLevel + 2, "coarse_type: ge" ); + appendLine( stream, indentLevel + 2, "num_sweeps: 1" ); + appendLine( stream, indentLevel + 2, "order: 0" ); + } + else + { + appendLine( stream, indentLevel + 2, "order: 1" ); + } #endif } @@ -1115,6 +1148,11 @@ void appendAMGByFlavor( std::ostringstream & stream, appendDisplacementAMG( stream, indentLevel, 0, false ); break; } + case AMGFlavor::almDisplacement: + { + appendDisplacementAMG( stream, indentLevel, mgrParams.separateComponents, true, true ); + break; + } } } @@ -1164,6 +1202,10 @@ bool buildStrategyYaml( LinearSolverParameters const & params, appendLine( stream, 2, "tolerance: 0.0" ); appendLine( stream, 2, "max_iter: 1" ); appendLine( stream, 2, GEOS_FMT( "print_level: {}", getMGRPrintLevel( params.logLevel ) ) ); + if( specialization.cycle != nullptr ) + { + appendLine( stream, 2, GEOS_FMT( "cycle: {}", specialization.cycle ) ); + } appendLine( stream, 2, GEOS_FMT( "non_c_to_f: {}", 1 ) ); appendLine( stream, 2, GEOS_FMT( "nonglk_max_elmts: {}", 1 ) ); appendLine( stream, 2, GEOS_FMT( "pmax: {}", specialization.pmax ) ); @@ -1269,6 +1311,8 @@ bool buildMGRPreconditionerYaml( LinearSolverParameters const & params, return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsEmbeddedFractures >( params, labelNames, numComponentsPerField, preconditionerYaml ); case StrategyType::singlePhasePoromechanicsConformingFractures: return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsConformingFractures >( params, labelNames, numComponentsPerField, preconditionerYaml ); + case StrategyType::singlePhasePoromechanicsConformingFracturesALM: + return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsConformingFracturesALM >( params, labelNames, numComponentsPerField, preconditionerYaml ); case StrategyType::singlePhasePoromechanicsReservoirFVM: return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsReservoirFVM >( params, labelNames, numComponentsPerField, preconditionerYaml ); case StrategyType::thermalSinglePhasePoromechanicsReservoirFVM: @@ -1580,6 +1624,11 @@ void HypredriveSolver::setExecutionContext( LinearSolverExecutionContext const & m_hasExecutionContext = true; } +void HypredriveSolver::setNearNullKernel( arrayView1d< HypreVector const > const & nearNullKernel ) +{ + m_nearNullKernel = nearNullKernel; +} + void HypredriveSolver::setup( HypreMatrix const & mat ) { Base::setup( mat ); @@ -1603,8 +1652,13 @@ void HypredriveSolver::createHypredrive( HypreMatrix const & mat, checkHypredriveCall( HYPREDRV_SetLibraryMode( m_hypredrive ), "HYPREDRV_SetLibraryMode" ); char * argv[] = { const_cast< char * >( parseTarget.argument.c_str() ) }; + std::string const parseContext = + parseTarget.source == hypre::hypredrive::InputSource::authoritativeFile + ? GEOS_FMT( "authoritative YAML file '{}'", parseTarget.argument ) + : "YAML generated by GEOS"; checkHypredriveCall( HYPREDRV_InputArgsParse( 1, argv, m_hypredrive ), - "HYPREDRV_InputArgsParse" ); + "HYPREDRV_InputArgsParse", + parseContext ); if( parseTarget.source == hypre::hypredrive::InputSource::generatedFallback && m_hasExecutionContext && !m_executionContext.solverName.empty() ) @@ -1650,6 +1704,32 @@ void HypredriveSolver::refreshBoundObjects( HypreMatrix const & mat, pointMarkers.data() ), "HYPREDRV_LinearSystemSetDofmap" ); } + + if( !m_nearNullKernel.empty() ) + { + localIndex const numEntries = mat.numLocalRows(); + localIndex const numModes = m_nearNullKernel.size(); + std::vector< HYPRE_Complex > values( numEntries * numModes ); + + for( localIndex mode = 0; mode < numModes; ++mode ) + { + GEOS_ERROR_IF( m_nearNullKernel[mode].localSize() != numEntries, + "HypreDrive near-null-space vector size does not match the matrix local size" ); + GEOS_LAI_CHECK_ERROR( + HYPRE_IJVectorGetValues( m_nearNullKernel[mode].unwrappedIJ(), + LvArray::integerConversion< HYPRE_Int >( numEntries ), + nullptr, + values.data() + mode * numEntries ) ); + } + + checkHypredriveCall( + HYPREDRV_LinearSystemSetNearNullSpace( + m_hypredrive, + LvArray::integerConversion< int >( numEntries ), + LvArray::integerConversion< int >( numModes ), + values.data() ), + "HYPREDRV_LinearSystemSetNearNullSpace" ); + } } void HypredriveSolver::setupLegacy( HypreMatrix const & mat ) diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.hpp index 5d70d0f124e..ab7adef728b 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.hpp @@ -98,6 +98,12 @@ class HypredriveSolver final : public LinearSolverBase< HypreInterface > */ void setExecutionContext( LinearSolverExecutionContext const & context ) override; + /** + * @brief Set near-null-space modes to pass to HypreDrive during setup. + * @param nearNullKernel Full-system distributed near-null-space vectors. + */ + void setNearNullKernel( arrayView1d< HypreVector const > const & nearNullKernel ) override; + /** * @brief Build or refresh the solver/preconditioner for a matrix. * @param mat Matrix that defines the system structure and coefficients. @@ -164,6 +170,7 @@ class HypredriveSolver final : public LinearSolverBase< HypreInterface > bool m_hasExecutionContext = false; bool m_timestepScopeActive = false; bool m_newtonScopeActive = false; + arrayView1d< HypreVector const > m_nearNullKernel; size_t m_hypredriveGeneration = 0; HYPREDRV_t m_hypredrive{}; mutable HypreVector m_dummyRhs; diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp new file mode 100644 index 00000000000..8799a8277a3 --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -0,0 +1,133 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 SinglePhasePoromechanicsConformingFracturesALM.hpp + */ + +#ifndef GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_ +#define GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_ + +#include "linearAlgebra/interfaces/hypre/HypreMGR.hpp" + +namespace geos +{ + +namespace hypre +{ + +namespace mgr +{ + +/** + * @brief SinglePhasePoromechanicsConformingFracturesALM 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 0 eliminates bubble displacement (3,4,5) with L1-Jacobi. + * 2. Level 1 eliminates nodal displacement (0,1,2) with one BoomerAMG V-cycle. + * 3. The displacement AMG uses three functions, separate-component filtering, + * and symmetric L1 hybrid Gauss-Seidel relaxation on CPU. + * 4. The pressure Schur complement is solved with BoomerAMG. + * 5. The MGR V-cycle uses pre-relaxation only and no global smoother. + */ +class SinglePhasePoromechanicsConformingFracturesALM : public MGRStrategyBase< 2 > +{ +public: + + /** + * @brief Constructor. + */ + explicit SinglePhasePoromechanicsConformingFracturesALM( arrayView1d< int const > const & ) + : MGRStrategyBase( 7 ) + { + + // Eliminate bubble displacement first, then nodal displacement. + m_labels[0] = { 0, 1, 2, 6 }; + m_labels[1] = { 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 ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetCycleType( precond.ptr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetFRelaxCycle( precond.ptr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetGlobalSmoothCycle( precond.ptr, 1 ) ); + + // Configure the BoomerAMG solver used as level-1 F-relaxation. + setDisplacementAMG( mgrData.mechSolver, mgrParams.separateComponents ); + +#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCoarsenType( mgrData.mechSolver.ptr, hypre::getAMGCoarseningType( LinearSolverParameters::AMG::CoarseningType::PMIS ) ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxType( mgrData.mechSolver.ptr, hypre::getAMGRelaxationType( LinearSolverParameters::AMG::SmootherType::chebyshev ) ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumSweeps( mgrData.mechSolver.ptr, 1 ) ); +#else + HYPRE_Int constexpr l1SymmetricHybridGaussSeidel = 89; + HYPRE_Int constexpr gaussianElimination = 9; + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCycleRelaxType( mgrData.mechSolver.ptr, l1SymmetricHybridGaussSeidel, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCycleRelaxType( mgrData.mechSolver.ptr, l1SymmetricHybridGaussSeidel, 2 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCycleRelaxType( mgrData.mechSolver.ptr, gaussianElimination, 3 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumSweeps( mgrData.mechSolver.ptr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxOrder( mgrData.mechSolver.ptr, 0 ) ); +#endif + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetFSolverAtLevel( precond.ptr, mgrData.mechSolver.ptr, 1 ) ); + + // 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_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_*/ diff --git a/src/coreComponents/linearAlgebra/unitTests/testHypredrive.cpp b/src/coreComponents/linearAlgebra/unitTests/testHypredrive.cpp index c004631a83d..d157a821ffc 100644 --- a/src/coreComponents/linearAlgebra/unitTests/testHypredrive.cpp +++ b/src/coreComponents/linearAlgebra/unitTests/testHypredrive.cpp @@ -261,6 +261,39 @@ TEST( HypredriveYaml, BuildsGeneratedYamlForEveryMGRStrategy ) } } +TEST( HypredriveYaml, BuildsSelectedALMPoromechanicsMGRStrategy ) +{ + stdVector< string > const fieldNames = { "totalDisplacement", "totalBubbleDisplacement", "pressure" }; + array1d< int > numComponentsPerField( 3 ); + numComponentsPerField[0] = 3; + numComponentsPerField[1] = 3; + numComponentsPerField[2] = 1; + + hypre::hypredrive::InputArgsParseTarget target; + ASSERT_TRUE( hypre::hypredrive::buildInputArgsParseTarget( + makeMgrParameters( LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsConformingFracturesALM ), + fieldNames, + numComponentsPerField, + target ) ); + + EXPECT_NE( target.argument.find( "num_levels: 3" ), std::string::npos ); + EXPECT_NE( target.argument.find( "cycle: v(1,0)" ), std::string::npos ); + EXPECT_NE( target.argument.find( "f_dofs: [totalBubbleDisplacement_0, totalBubbleDisplacement_1, totalBubbleDisplacement_2]" ), + std::string::npos ); + EXPECT_NE( target.argument.find( "f_dofs: [totalDisplacement_0, totalDisplacement_1, totalDisplacement_2]" ), + std::string::npos ); + EXPECT_NE( target.argument.find( "filter_functions: 1" ), std::string::npos ); +#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP + EXPECT_NE( target.argument.find( "down_type: 16" ), std::string::npos ); + EXPECT_NE( target.argument.find( "up_type: 16" ), std::string::npos ); +#else + EXPECT_NE( target.argument.find( "down_type: l1sym-hgs" ), std::string::npos ); + EXPECT_NE( target.argument.find( "up_type: l1sym-hgs" ), std::string::npos ); + EXPECT_NE( target.argument.find( "coarse_type: ge" ), std::string::npos ); + EXPECT_NE( target.argument.find( "order: 0" ), std::string::npos ); +#endif +} + TEST( HypredriveYaml, UsesSemanticCompositionalLabelsWhenFieldNamesAreAvailable ) { stdVector< string > const fieldNames = { "totalDisplacement", "compositionalVariables" }; diff --git a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp index d9c6a954a54..86479a5f857 100644 --- a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp +++ b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp @@ -295,42 +295,43 @@ struct LinearSolverParameters */ enum class StrategyType : integer { - invalid, ///< default value, to ensure solver sets something - singlePhaseReservoirFVM, ///< finite volume single-phase flow with wells - thermalSinglePhaseReservoirFVM, ///< finite volume thermal single-phase flow with wells - singlePhaseHybridFVM, ///< hybrid finite volume single-phase flow - singlePhaseReservoirHybridFVM, ///< hybrid finite volume single-phase flow with wells - singlePhasePoromechanics, ///< single phase poromechanics with finite volume single phase flow - thermalSinglePhasePoromechanics, ///< thermal single phase poromechanics with finite volume single phase flow - 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 - singlePhasePoromechanicsReservoirFVM, ///< single phase poromechanics with finite volume single phase flow with wells - thermalSinglePhasePoromechanicsReservoirFVM, ///< thermal single phase poromechanics with finite volume single phase flow with wells - compositionalMultiphaseFVM, ///< finite volume compositional multiphase flow - compositionalMultiphaseHybridFVM, ///< hybrid finite volume compositional multiphase flow - compositionalMultiphaseReservoirFVM, ///< finite volume compositional multiphase flow with wells - compositionalMultiphaseReservoirHybridFVM, ///< hybrid finite volume compositional multiphase flow with wells - immiscibleMultiphaseFVM, ///< finite volume immiscible multiphase flow - reactiveCompositionalMultiphaseOBL, ///< finite volume reactive compositional flow with OBL - thermalCompositionalMultiphaseFVM, ///< finite volume thermal compositional multiphase flow - thermalCompositionalMultiphaseReservoirFVM, ///< finite volume thermal compositional multiphase flow - multiphasePoromechanics, ///< multiphase poromechanics with finite volume compositional multiphase flow - multiphasePoromechanicsReservoirFVM, ///< multiphase poromechanics with finite volume compositional multiphase flow with wells - thermalMultiphasePoromechanics, ///< thermal multiphase poromechanics with finite volume compositional multiphase flow - hydrofracture, ///< hydrofracture - lagrangianContactMechanics, ///< Lagrangian contact mechanics - augmentedLagrangianContactMechanics, ///< Augmented Lagrangian contact mechanics - lagrangianContactMechanicsBubbleStab, ///< Lagrangian contact mechanics with bubble stabilization - solidMechanicsEmbeddedFractures ///< Embedded fractures mechanics + invalid, ///< default value, to ensure solver sets something + singlePhaseReservoirFVM, ///< finite volume single-phase flow with wells + thermalSinglePhaseReservoirFVM, ///< finite volume thermal single-phase flow with wells + singlePhaseHybridFVM, ///< hybrid finite volume single-phase flow + singlePhaseReservoirHybridFVM, ///< hybrid finite volume single-phase flow with wells + singlePhasePoromechanics, ///< single phase poromechanics with finite volume single phase flow + thermalSinglePhasePoromechanics, ///< thermal single phase poromechanics with finite volume single phase flow + 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 + thermalSinglePhasePoromechanicsReservoirFVM, ///< thermal single phase poromechanics with finite volume single phase flow with wells + compositionalMultiphaseFVM, ///< finite volume compositional multiphase flow + compositionalMultiphaseHybridFVM, ///< hybrid finite volume compositional multiphase flow + compositionalMultiphaseReservoirFVM, ///< finite volume compositional multiphase flow with wells + compositionalMultiphaseReservoirHybridFVM, ///< hybrid finite volume compositional multiphase flow with wells + immiscibleMultiphaseFVM, ///< finite volume immiscible multiphase flow + reactiveCompositionalMultiphaseOBL, ///< finite volume reactive compositional flow with OBL + thermalCompositionalMultiphaseFVM, ///< finite volume thermal compositional multiphase flow + thermalCompositionalMultiphaseReservoirFVM, ///< finite volume thermal compositional multiphase flow + multiphasePoromechanics, ///< multiphase poromechanics with finite volume compositional multiphase flow + multiphasePoromechanicsReservoirFVM, ///< multiphase poromechanics with finite volume compositional multiphase flow with wells + thermalMultiphasePoromechanics, ///< thermal multiphase poromechanics with finite volume compositional multiphase flow + hydrofracture, ///< hydrofracture + lagrangianContactMechanics, ///< Lagrangian contact mechanics + augmentedLagrangianContactMechanics, ///< Augmented Lagrangian contact mechanics + lagrangianContactMechanicsBubbleStab, ///< Lagrangian contact mechanics with bubble stabilization + solidMechanicsEmbeddedFractures ///< Embedded fractures mechanics }; - StrategyType strategy = StrategyType::invalid; ///< Predefined MGR solution strategy (solver specific) - integer separateComponents = false; ///< Apply a separate displacement component (SDC) filter before AMG construction - integer areWellsShut = false; ///< Flag to let MGR know that wells are shut, and that jacobi can be applied to the well - ///< block + StrategyType strategy = StrategyType::invalid; ///< Predefined MGR solution strategy (solver specific) + integer separateComponents = false; ///< Apply a separate displacement component (SDC) filter before AMG construction + integer areWellsShut = false; ///< Flag to let MGR know that wells are shut, and that jacobi can be applied to the well + ///< block } - mgr; ///< Multigrid reduction (MGR) parameters + mgr; ///< Multigrid reduction (MGR) parameters /// Incomplete factorization parameters struct IFact @@ -582,6 +583,7 @@ ENUM_STRINGS( LinearSolverParameters::MGR::StrategyType, "hybridSinglePhasePoromechanics", "singlePhasePoromechanicsEmbeddedFractures", "singlePhasePoromechanicsConformingFractures", + "singlePhasePoromechanicsConformingFracturesALM", "singlePhasePoromechanicsReservoirFVM", "thermalSinglePhasePoromechanicsReservoirFVM", "compositionalMultiphaseFVM", diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp index 6f82154c03f..d8e66284329 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp @@ -34,6 +34,9 @@ #include #include #include +#if VTK_VERSION_NUMBER >= VTK_VERSION_CHECK( 9, 6, 0 ) +#include +#endif #include #include #include @@ -68,6 +71,7 @@ #include #include #include +#include #ifdef GEOS_USE_MPI #include @@ -77,6 +81,8 @@ #endif #include +#include +#include #include @@ -980,7 +986,11 @@ static void classifyCellsByDimension( vtkDataSet & mesh, // Single pass: classify and populate for( vtkIdType i = 0; i < numCells; ++i ) { +#if VTK_VERSION_NUMBER >= VTK_VERSION_CHECK( 9, 6, 0 ) + int const dim = vtkCellTypeUtilities::GetDimension( mesh.GetCellType( i ) ); +#else int const dim = vtkCellTypes::GetDimension( mesh.GetCellType( i ) ); +#endif if( dim == 3 ) { cells3DIndices.emplace_back( i ); @@ -995,6 +1005,143 @@ static void classifyCellsByDimension( vtkDataSet & mesh, cells3DIndices.size(), cells2DIndices.size(), numCells ) ); } +using VTKPointCoordinate = std::array< real64, 3 >; + +struct VTKPointCoordinateHash +{ + std::size_t operator()( VTKPointCoordinate const & point ) const + { + std::size_t hash = 0; + for( real64 const coordinate: point ) + { + hash ^= std::hash< real64 >{}( coordinate ) + 0x9e3779b9 + ( hash << 6 ) + ( hash >> 2 ); + } + return hash; + } +}; + +static VTKPointCoordinate getPointCoordinate( vtkDataSet & mesh, vtkIdType const pointId ) +{ + double const * const point = mesh.GetPoint( pointId ); + return { point[0], point[1], point[2] }; +} + +/** + * @brief Find 3D cells whose faces match a 2D cell geometrically. + * + * Some mesh generators keep a separate point for each side of a split + * interface. In that case, the point ids of an embedded 2D cell do not + * match the point ids of either adjacent 3D face, even though the + * coordinates do. This routine performs the co-location lookup through + * point coordinates and then verifies that all points belong to one face. + * + * @param[in] mesh Original mesh + * @param[in] pointIds2D Point ids of the 2D cell + * @param[in] meshIdxToGlobalId3D Map from 3D mesh indices to global ids + * @param[in] pointsByCoordinate All mesh point ids grouped by coordinates + * @param[in,out] pointTo3DCells Cache of 3D cells incident to a mesh point + * @return Global ids of geometrically matching 3D cells + */ +static stdVector< int64_t > find2DTo3DNeighborsByCoordinates( + vtkDataSet & mesh, + vtkIdList * pointIds2D, + stdUnorderedMap< vtkIdType, int64_t > const & meshIdxToGlobalId3D, + stdUnorderedMap< VTKPointCoordinate, stdVector< vtkIdType >, VTKPointCoordinateHash > const & pointsByCoordinate, + stdUnorderedMap< vtkIdType, stdVector< vtkIdType > > & pointTo3DCells ) +{ + localIndex const numPoints = pointIds2D->GetNumberOfIds(); + + stdVector< VTKPointCoordinate > targetCoordinates; + targetCoordinates.reserve( numPoints ); + + stdUnorderedMap< vtkIdType, localIndex > candidateCounts; + + vtkNew< vtkIdList > pointCells; + for( vtkIdType i = 0; i < pointIds2D->GetNumberOfIds(); ++i ) + { + VTKPointCoordinate const coordinate = getPointCoordinate( mesh, pointIds2D->GetId( i ) ); + targetCoordinates.emplace_back( coordinate ); + + auto const coordinateIt = pointsByCoordinate.find( coordinate ); + if( coordinateIt == pointsByCoordinate.end() ) + { + continue; + } + + // A 3D cell can contain more than one mesh point with the same + // coordinate in a degenerate input. Count it only once for this 2D + // point when intersecting the incident-cell lists. + stdUnorderedMap< vtkIdType, bool > cellsAtCoordinate; + for( vtkIdType const meshPointId: coordinateIt->second ) + { + auto pointCellsIt = pointTo3DCells.find( meshPointId ); + if( pointCellsIt == pointTo3DCells.end() ) + { + stdVector< vtkIdType > & cachedCells = pointTo3DCells.get_inserted( meshPointId ); + pointCells->Reset(); + mesh.GetPointCells( meshPointId, pointCells ); + cachedCells.reserve( pointCells->GetNumberOfIds() ); + for( vtkIdType j = 0; j < pointCells->GetNumberOfIds(); ++j ) + { + vtkIdType const cellId = pointCells->GetId( j ); + if( meshIdxToGlobalId3D.count( cellId ) > 0 ) + { + cachedCells.emplace_back( cellId ); + } + } + pointCellsIt = pointTo3DCells.find( meshPointId ); + } + + for( vtkIdType const cellId: pointCellsIt->second ) + { + cellsAtCoordinate.emplace( cellId, true ); + } + } + + for( auto const & cell: cellsAtCoordinate ) + { + ++candidateCounts.get_inserted( cell.first ); + } + } + + std::sort( targetCoordinates.begin(), targetCoordinates.end() ); + + stdVector< int64_t > neighbors; + for( auto const & candidate: candidateCounts ) + { + if( candidate.second != numPoints ) + { + continue; + } + + vtkCell * const cell3D = mesh.GetCell( candidate.first ); + for( localIndex faceIndex = 0; faceIndex < cell3D->GetNumberOfFaces(); ++faceIndex ) + { + vtkCell * const face = cell3D->GetFace( faceIndex ); + if( face->GetNumberOfPoints() != numPoints ) + { + continue; + } + + stdVector< VTKPointCoordinate > faceCoordinates; + faceCoordinates.reserve( numPoints ); + for( vtkIdType i = 0; i < face->GetNumberOfPoints(); ++i ) + { + faceCoordinates.emplace_back( getPointCoordinate( mesh, face->GetPointId( i ) ) ); + } + std::sort( faceCoordinates.begin(), faceCoordinates.end() ); + + if( faceCoordinates == targetCoordinates ) + { + neighbors.emplace_back( meshIdxToGlobalId3D.at( candidate.first ) ); + break; + } + } + } + + return neighbors; +} + /** * @brief Build mapping from 2D cells to their neighboring 3D cells using indices * @@ -1028,6 +1175,18 @@ build2DTo3DNeighbors( vtkDataSet & mesh, ArrayOfArrays< localIndex, int64_t > neighbors2Dto3D; neighbors2Dto3D.reserve( cells2DIndices.size() ); + // Build a coordinate lookup for meshes with duplicated/collocated points. + // This is only needed when the input contains 2D cells; the lookup is + // deliberately kept local to this redistribution step. + stdUnorderedMap< VTKPointCoordinate, stdVector< vtkIdType >, VTKPointCoordinateHash > pointsByCoordinate; + pointsByCoordinate.reserve( mesh.GetNumberOfPoints() ); + for( vtkIdType pointId = 0; pointId < mesh.GetNumberOfPoints(); ++pointId ) + { + pointsByCoordinate.get_inserted( getPointCoordinate( mesh, pointId ) ).emplace_back( pointId ); + } + stdUnorderedMap< vtkIdType, stdVector< vtkIdType > > pointTo3DCells; + pointTo3DCells.reserve( mesh.GetNumberOfPoints() ); + // Topology statistics localIndex numStandalone = 0; localIndex numBoundary = 0; // 1 neighbor @@ -1063,6 +1222,24 @@ build2DTo3DNeighbors( vtkDataSet & mesh, // Non-3D neighbors (2D/1D/0D) are silently skipped } + // If point ids were duplicated at a split interface, the exact lookup + // above can miss one or both physical neighbors. Prefer the geometric + // connectivity whenever it is available, and retain the exact result + // as a fallback for meshes whose coordinates are not bitwise identical. + stdVector< int64_t > coordinateNeighborGlobalIds = find2DTo3DNeighborsByCoordinates( + mesh, + pointIds2D, + meshIdxToGlobalId3D, + pointsByCoordinate, + pointTo3DCells ); + if( !coordinateNeighborGlobalIds.empty() ) + { + neighbor3DGlobalIds.resize( coordinateNeighborGlobalIds.size() ); + std::copy( coordinateNeighborGlobalIds.begin(), + coordinateNeighborGlobalIds.end(), + neighbor3DGlobalIds.begin() ); + } + // Update topology statistics localIndex const numNeighbors = neighbor3DGlobalIds.size(); switch( numNeighbors ) diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp index 8d05dfd6958..f5dc5cdb3b1 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp @@ -1444,6 +1444,7 @@ void PhysicsSolverBase::solveLinearSystem( DofManager const & dofManager, executionContext.nonlinearIteration = nonlinearIteration; executionContext.systemSetupTimestamp = getSystemSetupTimestamp(); m_linearSolver->setExecutionContext( executionContext ); + m_linearSolver->setNearNullKernel( getLinearSolverNearNullKernel() ); if( isSetupNeeded ) { diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp index a174b485e67..107bea65a88 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp @@ -1016,6 +1016,14 @@ class PhysicsSolverBase : public ExecutableGroup virtual void postInputInitialization() override; + /** + * @brief Return optional near-null-space modes for the monolithic linear solver. + */ + virtual arrayView1d< ParallelVector const > getLinearSolverNearNullKernel() const + { + return {}; + } + /** * @brief Eisenstat-Walker adaptive tolerance * diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp index 468117b09dd..06591ccd0d4 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp @@ -25,6 +25,10 @@ #include "linearAlgebra/multiscale/MultiscalePreconditioner.hpp" #include "linearAlgebra/solvers/BlockPreconditioner.hpp" #include "linearAlgebra/solvers/SeparateComponentPreconditioner.hpp" +#include "linearAlgebra/utilities/LAIHelperFunctions.hpp" +#ifdef GEOS_USE_HYPREDRV +#include "linearAlgebra/interfaces/hypre/hypredrive.hpp" +#endif #include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsDamage.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp" @@ -78,6 +82,7 @@ void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::setupSystem( Dom // setup monolithic coupled system PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity ); + setupLinearSolverNearNullKernel( domain, dofManager ); if( !this->m_precond && this->m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) { @@ -85,6 +90,44 @@ void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::setupSystem( Dom } } +template< typename FLOW_SOLVER, typename MECHANICS_SOLVER > +void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::setupLinearSolverNearNullKernel( + DomainPartition & domain, + DofManager const & dofManager ) +{ + m_linearSolverNearNullKernel.resize( 0 ); + +#ifdef GEOS_USE_HYPREDRV + if( !hypre::hypredrive::shouldUse( this->m_linearSolverParameters.get() ) ) + { + return; + } + + this->solidMechanicsSolver()->forDiscretizationOnMeshTargets( + domain.getMeshBodies(), + [&]( string const &, MeshLevel const & mesh, string_array const & ) + { + if( !m_linearSolverNearNullKernel.empty() ) + { + return; + } + + NodeManager const & nodeManager = mesh.getNodeManager(); + arrayView1d< globalIndex const > const dofIndex = + nodeManager.getReference< array1d< globalIndex > >( + dofManager.getKey( solidMechanics::totalDisplacement::key() ) ); + m_linearSolverNearNullKernel = + LAIHelperFunctions::computeRigidBodyModes< ParallelVector >( + nodeManager.referencePosition(), + dofIndex, + dofManager.rankOffset(), + dofManager.numLocalDofs() ); + } ); +#else + GEOS_UNUSED_VAR( domain, dofManager ); +#endif +} + template< typename FLOW_SOLVER, typename MECHANICS_SOLVER > void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::initializePostInitialConditionsPreSubGroups() { diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp index d7e40171f79..518c7cac693 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp @@ -116,6 +116,14 @@ class SinglePhasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHAN virtual void initializePostInitialConditionsPreSubGroups() override; + void setupLinearSolverNearNullKernel( DomainPartition & domain, + DofManager const & dofManager ); + + arrayView1d< ParallelVector const > getLinearSolverNearNullKernel() const override + { + return m_linearSolverNearNullKernel.toViewConst(); + } + virtual void setMGRStrategy() override { if( this->m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::mgr ) @@ -155,6 +163,7 @@ class SinglePhasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHAN virtual string getFlowDofKey() const override { return SinglePhaseBase::viewKeyStruct::elemDofFieldString(); } + array1d< ParallelVector > m_linearSolverNearNullKernel; integer m_damageFlag; }; diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp index ff7d9662a60..0aca31b49a9 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,96 @@ 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 ) + { + // Start from both subsolver patterns. The flow pattern may contain well and + // flux couplings, while the mechanics pattern contains the nodal-bubble couplings. + SparsityPattern< globalIndex > flowPattern; + this->flowSolver()->setSparsityPattern( domain, dofManager, localMatrix, flowPattern ); + + SparsityPattern< globalIndex > mechanicsPattern; + this->solidMechanicsSolver()->setSparsityPattern( domain, dofManager, localMatrix, mechanicsPattern ); + + // Sum the row lengths to provide enough capacity for the union. Entries + // common to both patterns are deduplicated when they are inserted below. + array1d< localIndex > rowLengths( flowPattern.numRows() ); + for( localIndex localRow = 0; localRow < flowPattern.numRows(); ++localRow ) + { + rowLengths[localRow] = flowPattern.numNonZeros( localRow ) + mechanicsPattern.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 >( flowPattern.numRows(), + flowPattern.numColumns(), + rowLengths.data() ); + + // Copy the flow and mechanics nonzeros. + for( localIndex localRow = 0; localRow < flowPattern.numRows(); ++localRow ) + { + globalIndex const * cols = flowPattern.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + flowPattern.numNonZeros( localRow ) ); + } + for( localIndex localRow = 0; localRow < mechanicsPattern.numRows(); ++localRow ) + { + globalIndex const * cols = mechanicsPattern.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + mechanicsPattern.numNonZeros( localRow ) ); + } + + // Add the nonzeros from coupling + addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView() ); + addPressureForceCouplingPattern( domain, dofManager, pattern.toView() ); + addMatrixPressureBubbleCouplingPattern( domain, dofManager, pattern.toView() ); - 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 ); + + this->setupLinearSolverNearNullKernel( domain, dofManager ); + + if( !this->m_precond && this->m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) + { + this->m_precond = this->createPreconditioner( domain ); + } } template< typename FLOW_SOLVER > @@ -72,12 +170,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 ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + // Assemble flux terms and get dFluidResidual/dAperture + this->flowSolver()->assembleHydrofracFluxTerms( time_n, + dt, + domain, + dofManager, + localMatrix, + localRhs, + getDerivativeFluxResidual_dNormalJump() ); + + // Assemble coupling terms (must be after flux assembly) + assembleCouplingTerms( time_n, dt, domain, dofManager, localMatrix, localRhs ); } template< typename FLOW_SOLVER > @@ -90,10 +201,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 +229,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 ); + + // 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 ); + } ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + // 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 +397,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 +463,799 @@ 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()->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()->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) + 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) + 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) + 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) + 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; - 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( CellElementSubRegion::viewKeyStruct::bubbleCellsString() ) ) + 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(), [=, &rowLengths]( localIndex const kk ) + { + localIndex const k = bubbleElems[kk]; + localIndex const faceIndex = elemsToFaces[kk][0]; + + // (bubble_row, pressure_col): 1 pressure column for each of the 3 bubble DOFs + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowNumber = bubbleDofNumber[faceIndex] + i - rankOffset; + if( rowNumber >= 0 && rowNumber < rowLengths.size() ) + { + rowLengths[rowNumber] += 1; // One pressure DOF from matrix cell + } + } + + // (pressure_row, bubble_col): the matrix cell pressure couples to its 3 bubble DOFs (A_pb) + globalIndex const pRow = pressureDofNumber[k] - rankOffset; + if( pRow >= 0 && pRow < rowLengths.size() ) + { + rowLengths[pRow] += 3; // Three bubble DOFs + } + } ); + } ); + } ); +} + +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: +addMatrixPressureBubbleCouplingPattern( 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 & 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( CellElementSubRegion::viewKeyStruct::bubbleCellsString() ) ) + 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]; + + // (bubble_row, pressure_col) : A_bp + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowIndex = bubbleDofNumber[faceIndex] + i - rankOffset; + if( rowIndex >= 0 && rowIndex < pattern.numRows() ) + { + pattern.insertNonZero( rowIndex, pressureColIndex ); + } + } + + // (pressure_row, bubble_col) : A_pb -- transpose location + globalIndex const pRow = pressureDofNumber[k] - rankOffset; + if( pRow >= 0 && pRow < pattern.numRows() ) + { + for( localIndex i = 0; i < 3; ++i ) + { + pattern.insertNonZero( pRow, bubbleDofNumber[faceIndex] + i ); + } + } + } ); + } ); + } ); +} + +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, + FlowSolverBase::viewKeyStruct::fluidNamesString() ); + + real64 maxResidual = finiteElement::regionBasedKernelApplication + < parallelDevicePolicy<>, + constitutive::PorousSolidBase, + 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..58285027d01 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())); } /**@}*/ @@ -141,13 +146,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 +178,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..b05cbc61d35 --- /dev/null +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALMKernels.hpp @@ -0,0 +1,1135 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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" +#include "constitutive/fluid/singlefluid/SingleFluidBase.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 = density * area * dAperture/dU + // This is: dRpdU[j] = sum_i( unitNormal[i] * Atu[i][j] ) * density + // Use Ri_eq_AjiBj for transpose: R[i] = sum_j(A[j][i] * B[j]) + 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 + 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 + * + * The Atu/Atb matrices encode the sign convention directly: face 0 gets -N*detJ + * and face 1 gets +N*detJ. This means (1/area) * Atu^T * unitNormal already gives + * the correct dAperture/dU matching the LM convention (negative for face 0 nodes, + * positive for face 1 nodes). + */ +template< typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +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, + string const fluidModelKey ): + 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() ), + m_incrBubbleDisp( faceManager.getField< fields::contact::incrementalBubbleDisplacement >().toViewConst() ), + m_fluidDensity( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( + elementSubRegion.template getReference< string >( fluidModelKey ) ).density() ) + {} + + //*************************************************************************** + + /** + * @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{}, + bColIndices{}, + pColIndex( 0 ), + pRowIndex( 0 ), + localRb{}, + localdRbdP{}, + localdRpdB{}, + X{ {} }, + pLocal{} + {} + + /// C-array storage for the element local row degrees of freedom (bubble). + globalIndex bEqnRowIndices[3]; + + /// C-array storage for the (global) bubble column degrees of freedom (for A_pb). + globalIndex bColIndices[3]; + + /// Column index for pressure DOF (for A_bp) + globalIndex pColIndex; + + /// Row index (local) for pressure DOF (for A_pb) + globalIndex pRowIndex; + + /// C-array storage for the element local Rb residual vector. + real64 localRb[numBubbleUdofs]; + + /// C-array storage for the element local dRb/dP Jacobian. + real64 localdRbdP[numBubbleUdofs]; + + /// C-array storage for the element local dRp/dB Jacobian (A_pb bulk term). + real64 localdRpdB[numBubbleUdofs]; + + /// local nodal coordinates + real64 X[ numNodesPerElem ][ 3 ]; + + /// 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 ); + + // ---- A_pb^Omega : transpose Biot term for the fluid-mass equation ---- + // The bubble mode contributes to the cell volumetric strain, hence to the + // fluid mass storage. With dPorosity_dVolStrain = biot (BiotPorosity), + // dFluidMassIncrement_dVolStrainIncrement = biot * fluidDensity, + // so dR_p/db = fluidDensity * dR_b/dp (same geometric factor biot*div(beta), + // times the fluid density that puts it in mass units). We reuse dRbdP_gauss. + real64 const rhof = m_fluidDensity( k, 0 ); + LvArray::tensorOps::scaledAdd< nBubbleUdof >( stack.localdRpdB, dRbdP_gauss, -detJ * rhof ); + } + + /** + * @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]; + real64 localdRpdB[3]; + for( localIndex i = 0; i < 3; ++i ) + { + localRb[i] = stack.localRb[parentFaceIndex*3+i]; + localdRbdP[i] = stack.localdRbdP[parentFaceIndex*3+i]; + localdRpdB[i] = stack.localdRpdB[parentFaceIndex*3+i]; + } + + // ---- A_bp : bubble-momentum row (dR_b/dP) ---- + for( localIndex i=0; i < 3; ++i ) + { + localIndex const dof = LvArray::integerConversion< localIndex >( stack.bEqnRowIndices[ i ] ); + + 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 ); + } + + // ---- A_pb^Omega : fluid-mass row (dR_p/db) + storage residual ---- + localIndex const pRow = LvArray::integerConversion< localIndex >( stack.pRowIndex ); + if( pRow >= 0 && pRow < m_matrix.numRows() ) + { + // Residual: R_p += fluidDensity * biot * tr(grad beta) * delta_b (bubble storage), + // using the incremental bubble displacement of the bubble face. + localIndex const bFace = m_elemsToFaces[kk][0]; + real64 rp = 0.0; + for( localIndex i = 0; i < 3; ++i ) + { + rp += localdRpdB[i] * m_incrBubbleDisp[bFace][i]; + } + RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[pRow], rp ); + + // Jacobian dR_p/db (three bubble columns of the parent face) + m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( pRow, + stack.bColIndices, + localdRpdB, + 3 ); + } + + return 0.0; + } + +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; + + /// Incremental bubble displacement (face field) -- used for the fluid-mass storage residual. + arrayView2d< real64 const > const m_incrBubbleDisp; + + /// Fluid density [elem][q]. + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_fluidDensity; + +}; + +/// The factory used to construct MatrixPressureBubbleKernels. +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, + string const >; + +} // namespace poromechanicsMatrixBubbleKernels + +} // namespace geos + +#endif /* GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALMKERNELS_HPP_ */