diff --git a/src/coreComponents/common/format/Format.hpp b/src/coreComponents/common/format/Format.hpp index 910b7aa99fa..c8d485a4462 100644 --- a/src/coreComponents/common/format/Format.hpp +++ b/src/coreComponents/common/format/Format.hpp @@ -18,6 +18,8 @@ #include #include +#include +#include #define GEOS_USE_FMT diff --git a/src/coreComponents/linearAlgebra/CMakeLists.txt b/src/coreComponents/linearAlgebra/CMakeLists.txt index 00f692bede0..a7464c6b81b 100644 --- a/src/coreComponents/linearAlgebra/CMakeLists.txt +++ b/src/coreComponents/linearAlgebra/CMakeLists.txt @@ -73,6 +73,7 @@ set( linearAlgebra_headers utilities/ComponentMask.hpp utilities/InverseNormalOperator.hpp utilities/LAIHelperFunctions.hpp + utilities/SparsityPatternUtilities.hpp utilities/LinearSolverParameters.hpp utilities/LinearSolverResult.hpp utilities/NormalOperator.hpp @@ -168,6 +169,8 @@ 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/SinglePhasePoromechanicsConformingFracturesALMReservoirFVM.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..9ad2e8e4cb0 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp @@ -37,6 +37,8 @@ #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/SinglePhasePoromechanicsConformingFracturesALMReservoirFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsReservoirFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirHybridFVM.hpp" @@ -188,6 +190,16 @@ 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::singlePhasePoromechanicsConformingFracturesALMReservoirFVM: + { + setStrategy< SinglePhasePoromechanicsConformingFracturesALMReservoirFVM >( 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/HypreMGR.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.hpp index 8457dd8df09..c58d1d929fd 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.hpp @@ -39,6 +39,7 @@ struct HypreMGRData array1d< HYPRE_Int > pointMarkers; ///< array1d of unique tags for local degrees of freedom HyprePrecWrapper coarseSolver; ///< MGR coarse solver pointer and functions HyprePrecWrapper mechSolver; ///< MGR mechanics fine solver pointer and functions + HyprePrecWrapper nestedSolver; ///< Optional nested MGR F-relaxation wrapper }; namespace hypre @@ -98,6 +99,21 @@ class MGRStrategyBase // HYPRE_Int m_numRestrictSweeps{ -1 }; ///< Number of restrict sweeps // HYPRE_Int m_numInterpSweeps{ -1 }; ///< Number of interpolation sweeps + /** + * @brief Total number of dof labels, i.e. the sum of all fields' components. + * @param numComponentsPerField number of components of each field + * @return the total number of blocks + */ + static HYPRE_Int totalNumBlocks( arrayView1d< int const > const & numComponentsPerField ) + { + HYPRE_Int result = 0; + for( localIndex i = 0; i < numComponentsPerField.size(); ++i ) + { + result += LvArray::integerConversion< HYPRE_Int >( numComponentsPerField[i] ); + } + return result; + } + /** * @brief Constructor. * @param numBlocks number of blocks @@ -196,6 +212,58 @@ class MGRStrategyBase solver.destroy = HYPRE_BoomerAMGDestroy; } + /** + * @brief Set up one of the two BoomerAMG instances used by the fully + * coupled single-phase ALM hierarchy. + * @param solver solver wrapper to initialize + * @param separateComponents whether displacement components are filtered + * @param bubbleCoarse true for the inner (bubble-displacement) coarse solve + * + * The CPU values mirror the reference nested-MGR YAML. The device branch + * retains the device-safe smoother choices used by the existing strategies. + */ + void setALMDisplacementAMG( HyprePrecWrapper & solver, + integer const separateComponents, + bool const bubbleCoarse ) + { + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGCreate( &solver.ptr ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetTol( solver.ptr, 0.0 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxIter( solver.ptr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPrintLevel( solver.ptr, 0 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxRowSum( solver.ptr, 1.0 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetStrongThreshold( solver.ptr, bubbleCoarse ? 0.75 : 0.8 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumFunctions( solver.ptr, 3 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetFilterFunctions( solver.ptr, bubbleCoarse ? 0 : separateComponents ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPMaxElmts( solver.ptr, bubbleCoarse ? 10 : 20 ) ); + + if( !bubbleCoarse ) + { + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetAggNumLevels( solver.ptr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCoarsenType( solver.ptr, + hypre::getAMGCoarseningType( LinearSolverParameters::AMG::CoarseningType::Falgout ) ) ); + } + +#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCoarsenType( solver.ptr, + hypre::getAMGCoarseningType( LinearSolverParameters::AMG::CoarseningType::PMIS ) ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxType( solver.ptr, + hypre::getAMGRelaxationType( LinearSolverParameters::AMG::SmootherType::chebyshev ) ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumSweeps( solver.ptr, bubbleCoarse ? 1 : 2 ) ); +#else + HYPRE_Int constexpr l1SymmetricHybridGaussSeidel = 89; + HYPRE_Int constexpr gaussianElimination = 9; + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCycleRelaxType( solver.ptr, l1SymmetricHybridGaussSeidel, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCycleRelaxType( solver.ptr, l1SymmetricHybridGaussSeidel, 2 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCycleRelaxType( solver.ptr, gaussianElimination, 3 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumSweeps( solver.ptr, bubbleCoarse ? 1 : 2 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxOrder( solver.ptr, 0 ) ); +#endif + + solver.setup = HYPRE_BoomerAMGSetup; + solver.solve = HYPRE_BoomerAMGSolve; + solver.destroy = HYPRE_BoomerAMGDestroy; + } + /** * @brief Set up BoomerAMG to perform the solve for the pressure system * @param solver the solver wrapper diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp index ff374dcc982..1a6e2691b34 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp @@ -453,6 +453,10 @@ void HyprePreconditioner::clear() { GEOS_LAI_CHECK_ERROR( m_mgrData->mechSolver.destroy( m_mgrData->mechSolver.ptr ) ); } + if( m_mgrData && m_mgrData->nestedSolver.ptr && m_mgrData->nestedSolver.destroy ) + { + GEOS_LAI_CHECK_ERROR( m_mgrData->nestedSolver.destroy( m_mgrData->nestedSolver.ptr ) ); + } m_precond.reset(); m_mgrData.reset(); } diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp index 9a943892ebc..881514865db 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/hypredrive.cpp @@ -24,6 +24,8 @@ #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/SinglePhasePoromechanicsConformingFracturesALMReservoirFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsReservoirFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirFVM.hpp" @@ -41,6 +43,7 @@ #include #include +#include #include #include #include @@ -59,12 +62,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 ) ); } } @@ -413,6 +430,8 @@ bool strategyUsesCompositionalSemanticLabels( LinearSolverParameters::MGR::Strat case StrategyType::hybridSinglePhasePoromechanics: case StrategyType::singlePhasePoromechanicsEmbeddedFractures: case StrategyType::singlePhasePoromechanicsConformingFractures: + case StrategyType::singlePhasePoromechanicsConformingFracturesALM: + case StrategyType::singlePhasePoromechanicsConformingFracturesALMReservoirFVM: case StrategyType::singlePhasePoromechanicsReservoirFVM: case StrategyType::thermalSinglePhasePoromechanicsReservoirFVM: case StrategyType::hydrofracture: @@ -881,7 +900,8 @@ enum class AMGFlavor pressure, pressureTemperature, displacementFiltered, - displacement + displacement, + almDisplacement }; struct LevelAMGBlock @@ -918,6 +938,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 ) @@ -968,6 +989,22 @@ MGRSpecialization getSpecialization( LinearSolverParameters::MGR::StrategyType c MGRSpecialization specialization; specialization.coarseFlavor = AMGFlavor::pressure; specialization.fRelaxAMGLevels = { LevelAMGBlock{ 1, AMGFlavor::displacement } }; + specialization.cycle = "v(1,0)"; + return specialization; + } + case StrategyType::singlePhasePoromechanicsConformingFracturesALM: + { + MGRSpecialization specialization; + specialization.coarseFlavor = AMGFlavor::pressure; + specialization.cycle = "v(1,0)"; + return specialization; + } + case StrategyType::singlePhasePoromechanicsConformingFracturesALMReservoirFVM: + { + MGRSpecialization specialization; + specialization.coarseFlavor = AMGFlavor::pressure; + specialization.fRelaxAMGLevels = { LevelAMGBlock{ 1, AMGFlavor::almDisplacement } }; + specialization.cycle = "v(1,0)"; return specialization; } case StrategyType::invalid: @@ -1034,7 +1071,8 @@ void appendHypreBoomerAMGCreateMaxCoarseSize( 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:" ); @@ -1052,7 +1090,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 } @@ -1089,6 +1138,93 @@ void appendPressureAMG( std::ostringstream & stream, #endif } +void appendALMDisplacementFineAMG( std::ostringstream & stream, + integer const indentLevel ) +{ + appendAMGHeader( stream, indentLevel ); + appendLine( stream, indentLevel + 1, "aggressive:" ); + appendLine( stream, indentLevel + 2, "num_levels: 1" ); + appendLine( stream, indentLevel + 1, "interpolation:" ); + appendLine( stream, indentLevel + 2, "max_nnz_row: 20" ); + appendLine( stream, indentLevel + 1, "coarsening:" ); +#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP + appendLine( stream, indentLevel + 2, "type: 8" ); +#else + appendLine( stream, indentLevel + 2, "type: falgout" ); +#endif + appendLine( stream, indentLevel + 2, "max_row_sum: 1.0" ); + appendLine( stream, indentLevel + 2, "strong_th: 0.8" ); + appendLine( stream, indentLevel + 2, "num_functions: 3" ); + appendLine( stream, indentLevel + 2, "filter_functions: 1" ); + appendLine( stream, indentLevel + 1, "relaxation:" ); +#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP + appendLine( stream, indentLevel + 2, "down_type: 16" ); + appendLine( stream, indentLevel + 2, "up_type: 16" ); + appendLine( stream, indentLevel + 2, "coarse_type: 16" ); + appendLine( stream, indentLevel + 2, "num_sweeps: 2" ); +#else + 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: 2" ); + appendLine( stream, indentLevel + 2, "order: 0" ); +#endif +} + +void appendALMDisplacementBubbleAMG( std::ostringstream & stream, + integer const indentLevel ) +{ + appendAMGHeader( stream, indentLevel ); + appendLine( stream, indentLevel + 1, "interpolation:" ); + appendLine( stream, indentLevel + 2, "max_nnz_row: 10" ); + appendLine( stream, indentLevel + 1, "coarsening:" ); +#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP + appendLine( stream, indentLevel + 2, "type: 8" ); +#endif + appendLine( stream, indentLevel + 2, "max_row_sum: 1.0" ); + appendLine( stream, indentLevel + 2, "strong_th: 0.75" ); + appendLine( stream, indentLevel + 2, "num_functions: 3" ); + appendLine( stream, indentLevel + 2, "filter_functions: 0" ); + appendLine( stream, indentLevel + 1, "relaxation:" ); +#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP + appendLine( stream, indentLevel + 2, "down_type: 16" ); + appendLine( stream, indentLevel + 2, "up_type: 16" ); + appendLine( stream, indentLevel + 2, "coarse_type: 16" ); + appendLine( stream, indentLevel + 2, "num_sweeps: 1" ); +#else + 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" ); +#endif +} + +void appendALMNestedMGR( std::ostringstream & stream, + integer const indentLevel, + stdVector< string > const & labelNames ) +{ + stdVector< HYPRE_Int > const nodalLabels = { 0, 1, 2 }; + appendLine( stream, indentLevel, "mgr:" ); + appendLine( stream, indentLevel + 1, "tolerance: 0.0" ); + appendLine( stream, indentLevel + 1, "max_iter: 1" ); + appendLine( stream, indentLevel + 1, "print_level: 0" ); + appendLine( stream, indentLevel + 1, "cycle: v(1,0)" ); + appendLine( stream, indentLevel + 1, "num_levels: 2" ); + appendLine( stream, indentLevel + 1, "level:" ); + appendLine( stream, indentLevel + 2, "0:" ); + appendLine( stream, indentLevel + 3, + GEOS_FMT( "f_dofs: [{}]", joinLabelNames( nodalLabels, labelNames ) ) ); + appendLine( stream, indentLevel + 3, "f_relaxation:" ); + appendALMDisplacementFineAMG( stream, indentLevel + 4 ); + appendLine( stream, indentLevel + 3, "g_relaxation: none" ); + appendLine( stream, indentLevel + 3, "prolongation_type: injection" ); + appendLine( stream, indentLevel + 3, "restriction_type: injection" ); + appendLine( stream, indentLevel + 3, "coarse_level_type: rap" ); + appendLine( stream, indentLevel + 1, "coarsest_level:" ); + appendALMDisplacementBubbleAMG( stream, indentLevel + 2 ); +} + void appendPressureTemperatureAMG( std::ostringstream & stream, integer const indentLevel ) { @@ -1200,6 +1336,11 @@ void appendAMGByFlavor( std::ostringstream & stream, appendDisplacementAMG( stream, indentLevel, 0, false ); break; } + case AMGFlavor::almDisplacement: + { + appendDisplacementAMG( stream, indentLevel, mgrParams.separateComponents, true, true ); + break; + } } } @@ -1244,6 +1385,45 @@ bool buildStrategyYaml( LinearSolverParameters const & params, strategy.setup( params.mgr, precond, mgrData ); std::ostringstream stream; + if( params.mgr.strategy == LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsConformingFracturesALM ) + { + // Keep this representation in lockstep with the fully coupled ALM + // strategy: the outer F block is displacement plus bubble displacement, + // and its F-relaxation is a two-level nested MGR. + appendLine( stream, 0, "preconditioner:" ); + appendLine( stream, 1, "mgr:" ); + appendLine( stream, 2, "tolerance: 0.0" ); + appendLine( stream, 2, "max_iter: 1" ); + appendLine( stream, 2, GEOS_FMT( "print_level: {}", getMGRPrintLevel( params.logLevel ) ) ); + appendLine( stream, 2, "cycle: v(1,0)" ); + appendLine( stream, 2, "non_c_to_f: 1" ); + appendLine( stream, 2, "nonglk_max_elmts: 1" ); + appendLine( stream, 2, "pmax: 0" ); + appendLine( stream, 2, GEOS_FMT( "coarse_th: {}", strategy.m_coarseGridThreshold ) ); + appendLine( stream, 2, "num_levels: 2" ); + appendLine( stream, 2, "level:" ); + appendLine( stream, 3, "0:" ); + stdVector< HYPRE_Int > const displacementLabels = { 0, 1, 2, 3, 4, 5 }; + appendLine( stream, 4, + GEOS_FMT( "f_dofs: [{}]", joinLabelNames( displacementLabels, labelNames ) ) ); + appendLine( stream, 4, "f_relaxation:" ); + appendALMNestedMGR( stream, 5, labelNames ); + appendLine( stream, 4, "g_relaxation: none" ); + appendLine( stream, 4, "prolongation_type: 2" ); + appendLine( stream, 4, "restriction_type: injection" ); + appendLine( stream, 4, "coarse_level_type: rap" ); + appendLine( stream, 2, "coarsest_level:" ); + appendPressureAMG( stream, 3, -1 ); + + preconditionerYaml = stream.str(); + + destroyWrapper( mgrData.coarseSolver ); + destroyWrapper( mgrData.mechSolver ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRDestroy( precond.ptr ) ); + destroyWrapper( mgrData.nestedSolver ); + return true; + } + appendLine( stream, 0, "preconditioner:" ); appendLine( stream, 1, "mgr:" ); appendLine( stream, 2, "tolerance: 0.0" ); @@ -1253,6 +1433,10 @@ bool buildStrategyYaml( LinearSolverParameters const & params, // hypre's MGR setup also rewrites F-relax type 0 to 7 unless interp is 12, // so omitting the key matches the effective legacy behavior. 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 ) ); @@ -1319,6 +1503,7 @@ bool buildStrategyYaml( LinearSolverParameters const & params, destroyWrapper( mgrData.coarseSolver ); destroyWrapper( mgrData.mechSolver ); GEOS_LAI_CHECK_ERROR( HYPRE_MGRDestroy( precond.ptr ) ); + destroyWrapper( mgrData.nestedSolver ); return true; } @@ -1359,6 +1544,10 @@ 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::singlePhasePoromechanicsConformingFracturesALMReservoirFVM: + return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsConformingFracturesALMReservoirFVM >( params, labelNames, numComponentsPerField, preconditionerYaml ); case StrategyType::singlePhasePoromechanicsReservoirFVM: return buildStrategyYaml< hypre::mgr::SinglePhasePoromechanicsReservoirFVM >( params, labelNames, numComponentsPerField, preconditionerYaml ); case StrategyType::thermalSinglePhasePoromechanicsReservoirFVM: @@ -1670,6 +1859,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 ); @@ -1702,8 +1896,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() ) @@ -1762,6 +1961,36 @@ void HypredriveSolver::refreshBoundObjects( HypreMatrix const & mat, LvArray::integerConversion< int >( pointMarkers.size() ), pointMarkers.data() ), "HYPREDRV_LinearSystemSetDofmap" ); + + if( !m_nearNullKernel.empty() ) + { + localIndex const numEntries = mat.numLocalRows(); + localIndex const numModes = m_nearNullKernel.size(); + array1d< HYPRE_Complex > values; + values.resizeWithoutInitializationOrDestruction( hypre::memorySpace, 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 ) ); + } + + values.registerTouch( hypre::memorySpace ); + values.move( hostMemorySpace ); + + 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 54cdc909943..c8af254e778 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; bool m_reportedGeneratedYamlFailure = false; size_t m_hypredriveGeneration = 0; HYPREDRV_t m_hypredrive{}; 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..2dddebeeb78 --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -0,0 +1,293 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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" + +#include + +namespace geos +{ + +namespace hypre +{ + +namespace mgr +{ + +/** + * @brief Adapter that lets an MGR instance be used as an MGR F-relaxation + * solver. + * + * HYPRE's public MGR API accepts callbacks for the first F-relaxation, while + * subsequent solves dispatch through the base hypre_Solver layout. Keep the + * callback members first so that the adapter is valid in both paths. The + * nested MGR and its two user-owned AMG solvers are released by destroy(). + */ +struct SinglePhaseALMNestedMGR +{ + HyprePrecWrapper::SetupFunc setup{}; + HyprePrecWrapper::SolveFunc solve{}; + HyprePrecWrapper::DestroyFunc destroy{}; + HYPRE_Int isSetup{ 0 }; + HYPRE_Solver mgr{}; + HyprePrecWrapper displacementSolver{}; + HyprePrecWrapper bubbleSolver{}; + stdVector< HYPRE_Int > pointMarkers; +}; + +inline SinglePhaseALMNestedMGR * singlePhaseALMNestedMGR( HYPRE_Solver solver ) +{ + return reinterpret_cast< SinglePhaseALMNestedMGR * >( solver ); +} + +inline HYPRE_Int singlePhaseALMNestedMGRSetup( HYPRE_Solver solver, + HYPRE_ParCSRMatrix A, + HYPRE_ParVector b, + HYPRE_ParVector x ) +{ + SinglePhaseALMNestedMGR * const nested = singlePhaseALMNestedMGR( solver ); + if( nested == nullptr || nested->mgr == nullptr ) + { + return 1; + } + + HYPRE_Int const ierr = HYPRE_MGRSetup( nested->mgr, A, b, x ); + if( ierr == 0 ) + { + nested->isSetup = 1; + } + return ierr; +} + +inline HYPRE_Int singlePhaseALMNestedMGRSolve( HYPRE_Solver solver, + HYPRE_ParCSRMatrix A, + HYPRE_ParVector b, + HYPRE_ParVector x ) +{ + SinglePhaseALMNestedMGR * const nested = singlePhaseALMNestedMGR( solver ); + if( nested == nullptr || nested->mgr == nullptr ) + { + return 1; + } + return HYPRE_MGRSolve( nested->mgr, A, b, x ); +} + +inline HYPRE_Int singlePhaseALMNestedMGRDestroy( HYPRE_Solver solver ) +{ + SinglePhaseALMNestedMGR * const nested = singlePhaseALMNestedMGR( solver ); + if( nested == nullptr ) + { + return 0; + } + + HYPRE_Int ierr = 0; + if( nested->mgr != nullptr ) + { + ierr |= HYPRE_MGRDestroy( nested->mgr ); + nested->mgr = nullptr; + } + if( nested->displacementSolver.ptr != nullptr && nested->displacementSolver.destroy != nullptr ) + { + ierr |= nested->displacementSolver.destroy( nested->displacementSolver.ptr ); + nested->displacementSolver.ptr = nullptr; + } + if( nested->bubbleSolver.ptr != nullptr && nested->bubbleSolver.destroy != nullptr ) + { + ierr |= nested->bubbleSolver.destroy( nested->bubbleSolver.ptr ); + nested->bubbleSolver.ptr = nullptr; + } + delete nested; + return ierr; +} + +static_assert( offsetof( SinglePhaseALMNestedMGR, isSetup ) == + 3 * sizeof( HYPRE_PtrToSolverFcn ), + "The nested MGR adapter must start with the hypre_Solver callback layout" ); + +/** + * @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) + * + * Well unknowns are not handled here: see + * SinglePhasePoromechanicsConformingFracturesALMReservoirFVM, which adds a + * third level to eliminate the well block before the coarse solve. + * + * Ingredients: + * 1. The outer MGR F-block is the fully coupled displacement/bubble block. + * 2. Its F-relaxation is a nested MGR that eliminates nodal displacement + * (0,1,2) and leaves bubble displacement (3,4,5) for the inner coarse AMG. + * 3. Both displacement AMG solves use three functions and the ALM smoother + * settings used by the reference strategy. + * 4. The pressure Schur complement is solved with BoomerAMG. + * 5. Both MGR V-cycles use pre-relaxation only and no global smoother. + */ +class SinglePhasePoromechanicsConformingFracturesALM : public MGRStrategyBase< 1 > +{ +public: + + /** + * @brief Constructor. + */ + explicit SinglePhasePoromechanicsConformingFracturesALM( arrayView1d< int const > const & numComponentsPerField ) + : MGRStrategyBase( totalNumBlocks( numComponentsPerField ) ) + { + GEOS_ERROR_IF_NE_MSG( numComponentsPerField.size(), 3, + "singlePhasePoromechanicsConformingFracturesALM requires exactly the displacement, " + "bubble-displacement, and pressure fields. Use " + "singlePhasePoromechanicsConformingFracturesALMReservoirFVM when wells are present." ); + GEOS_ERROR_IF_NE_MSG( numComponentsPerField[0], 3, + "singlePhasePoromechanicsConformingFracturesALM requires three displacement components" ); + GEOS_ERROR_IF_NE_MSG( numComponentsPerField[1], 3, + "singlePhasePoromechanicsConformingFracturesALM requires three bubble-displacement components" ); + + // The outer MGR eliminates the fully coupled displacement/bubble block. + // Its F-relaxation is a nested MGR that eliminates nodal displacement and + // leaves the bubble displacement block for the inner coarse AMG solve. + HYPRE_Int const numDisplacementLabels = LvArray::integerConversion< HYPRE_Int >( numComponentsPerField[0] ); + HYPRE_Int const pressureLabel = numDisplacementLabels + + LvArray::integerConversion< HYPRE_Int >( numComponentsPerField[1] ); + for( HYPRE_Int label = pressureLabel; label < m_numBlocks; ++label ) + { + m_labels[0].push_back( label ); + } + + setupLabels(); + + // Level 0: nested MGR F-relaxation over all six displacement unknowns. + // HYPRE requires the AMG F-relaxation enum for the callback-based solver. + m_levelFRelaxType[0] = MGRFRelaxationType::amgVCycle; + m_levelFRelaxIters[0] = 1; + m_levelGlobalSmootherType[0] = MGRGlobalSmootherType::none; + m_levelGlobalSmootherIters[0] = 0; + m_levelInterpType[0] = MGRInterpolationType::jacobi; + m_levelRestrictType[0] = MGRRestrictionType::injection; + m_levelCoarseGridMethod[0] = MGRCoarseGridMethod::galerkin; + + } + + /** + * @brief Setup the MGR strategy. + * @param precond preconditioner wrapper + * @param mgrData auxiliary MGR data + */ + void setup( LinearSolverParameters::MGR const & mgrParams, + HyprePrecWrapper & precond, + HypreMGRData & mgrData ) + { + GEOS_UNUSED_VAR( mgrParams ); + 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 ) ); + + auto * const nested = new SinglePhaseALMNestedMGR; + nested->setup = singlePhaseALMNestedMGRSetup; + nested->solve = singlePhaseALMNestedMGRSolve; + nested->destroy = singlePhaseALMNestedMGRDestroy; + mgrData.nestedSolver.ptr = reinterpret_cast< HYPRE_Solver >( nested ); + mgrData.nestedSolver.setup = nested->setup; + mgrData.nestedSolver.solve = nested->solve; + mgrData.nestedSolver.destroy = nested->destroy; + + // The outer F block contains labels 0..5. The marker order must match + // HYPRE's projected A_FF row order, hence the filtering of the original + // local marker array rather than constructing a synthetic marker list. + nested->pointMarkers.reserve( mgrData.pointMarkers.size() ); + for( HYPRE_Int const marker : mgrData.pointMarkers ) + { + if( marker < 6 ) + { + nested->pointMarkers.push_back( marker ); + } + } + + GEOS_LAI_CHECK_ERROR( HYPRE_MGRCreate( &nested->mgr ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetTol( nested->mgr, 0.0 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetMaxIter( nested->mgr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetPrintLevel( nested->mgr, 0 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetCycleType( nested->mgr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetFRelaxCycle( nested->mgr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetGlobalSmoothCycle( nested->mgr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetNonCpointsToFpoints( nested->mgr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetNonGalerkinMaxElmts( nested->mgr, 1 ) ); + + HYPRE_Int innerCoarseLabels[3] = { 3, 4, 5 }; + HYPRE_Int * innerCoarseLabelsPtr[1] = { innerCoarseLabels }; + HYPRE_Int numInnerCoarseLabels[1] = { 3 }; + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetCpointsByPointMarkerArray( nested->mgr, + 6, 1, + numInnerCoarseLabels, + innerCoarseLabelsPtr, + nested->pointMarkers.data() ) ); + + HYPRE_Int fRelaxType[1] = { static_cast< HYPRE_Int >( MGRFRelaxationType::amgVCycle ) }; + HYPRE_Int fRelaxIters[1] = { 1 }; + HYPRE_Int interpolationType[1] = { static_cast< HYPRE_Int >( MGRInterpolationType::injection ) }; + HYPRE_Int restrictionType[1] = { static_cast< HYPRE_Int >( MGRRestrictionType::injection ) }; + HYPRE_Int coarseGridMethod[1] = { static_cast< HYPRE_Int >( MGRCoarseGridMethod::galerkin ) }; + HYPRE_Int globalSmoothType[1] = { static_cast< HYPRE_Int >( MGRGlobalSmootherType::none ) }; + HYPRE_Int globalSmoothIters[1] = { 0 }; + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetLevelFRelaxType( nested->mgr, fRelaxType ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetLevelNumRelaxSweeps( nested->mgr, fRelaxIters ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetLevelInterpType( nested->mgr, interpolationType ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetLevelRestrictType( nested->mgr, restrictionType ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetCoarseGridMethod( nested->mgr, coarseGridMethod ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetLevelSmoothType( nested->mgr, globalSmoothType ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetLevelSmoothIters( nested->mgr, globalSmoothIters ) ); + + // Inner level 0: AMG for nodal displacement; inner coarsest level: + // AMG for the bubble-displacement block. + setALMDisplacementAMG( nested->displacementSolver, 1, false ); + setALMDisplacementAMG( nested->bubbleSolver, 0, true ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetFSolverAtLevel( nested->mgr, nested->displacementSolver.ptr, 0 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetCoarseSolver( nested->mgr, + nested->bubbleSolver.solve, + nested->bubbleSolver.setup, + nested->bubbleSolver.ptr ) ); + + // HYPRE invokes the callback during setup and subsequently dispatches via + // the adapter's hypre_Solver-compatible first four fields. + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetFSolver( precond.ptr, + nested->solve, + nested->setup, + mgrData.nestedSolver.ptr ) ); + + // Configure the BoomerAMG solver used as the outer 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/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALMReservoirFVM.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALMReservoirFVM.hpp new file mode 100644 index 00000000000..8ec35a44784 --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALMReservoirFVM.hpp @@ -0,0 +1,168 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 SinglePhasePoromechanicsConformingFracturesALMReservoirFVM.hpp + */ + +#ifndef GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALMRESERVOIRFVM_HPP_ +#define GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALMRESERVOIRFVM_HPP_ + +#include "linearAlgebra/interfaces/hypre/HypreMGR.hpp" + +namespace geos +{ + +namespace hypre +{ + +namespace mgr +{ + +/** + * @brief SinglePhasePoromechanicsConformingFracturesALMReservoirFVM 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) + * dofLabel: 7 = well pressure + * dofLabel: 8 = well rate + * + * 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. Level 2 eliminates the well block (7,8), as in + * SinglePhasePoromechanicsReservoirFVM. Leaving the well unknowns in the + * coarse grid would hand BoomerAMG the well rate and BHP constraint rows, + * which have no elliptic structure. + * 4. The displacement AMG uses three functions, separate-component filtering, + * and symmetric L1 hybrid Gauss-Seidel relaxation on CPU. + * 5. The reservoir pressure Schur complement is solved with BoomerAMG. + * 6. The MGR V-cycle uses pre-relaxation only and no global smoother. + */ +class SinglePhasePoromechanicsConformingFracturesALMReservoirFVM : public MGRStrategyBase< 3 > +{ +public: + + /** + * @brief Constructor. + * @param numComponentsPerField array with number of components for each field + */ + explicit SinglePhasePoromechanicsConformingFracturesALMReservoirFVM( arrayView1d< int const > const & numComponentsPerField ) + : MGRStrategyBase( totalNumBlocks( numComponentsPerField ) ) + { + GEOS_ERROR_IF_NE_MSG( numComponentsPerField.size(), 4, + "singlePhasePoromechanicsConformingFracturesALMReservoirFVM requires exactly the " + "displacement, bubble-displacement, pressure, and well fields. Any further field would " + "be swept into the well block eliminated on level 2." ); + + // Eliminate bubble displacement first, then nodal displacement, then the + // well block. Only the reservoir pressure survives on the coarsest grid. + HYPRE_Int const numDisplacementLabels = LvArray::integerConversion< HYPRE_Int >( numComponentsPerField[0] ); + HYPRE_Int const pressureLabel = numDisplacementLabels + + LvArray::integerConversion< HYPRE_Int >( numComponentsPerField[1] ); + HYPRE_Int const wellLabel = pressureLabel + + LvArray::integerConversion< HYPRE_Int >( numComponentsPerField[2] ); + for( HYPRE_Int label = 0; label < numDisplacementLabels; ++label ) + { + m_labels[0].push_back( label ); + } + for( HYPRE_Int label = pressureLabel; label < m_numBlocks; ++label ) + { + m_labels[0].push_back( label ); + m_labels[1].push_back( label ); + } + for( HYPRE_Int label = pressureLabel; label < wellLabel; ++label ) + { + m_labels[2].push_back( label ); + } + + 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; + + // Level 2 + m_levelFRelaxType[2] = MGRFRelaxationType::gsElimWInverse; + m_levelFRelaxIters[2] = 1; + m_levelGlobalSmootherType[2] = MGRGlobalSmootherType::none; + m_levelInterpType[2] = MGRInterpolationType::blockJacobi; + m_levelRestrictType[2] = MGRRestrictionType::injection; + m_levelCoarseGridMethod[2] = MGRCoarseGridMethod::galerkin; + } + + /** + * @brief Setup the MGR strategy. + * @param mgrParams MGR configuration parameters + * @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 reservoir pressure reduced system + setPressureAMG( mgrData.coarseSolver ); + } +}; + +} // namespace mgr + +} // namespace hypre + +} // namespace geos + +#endif /*GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALMRESERVOIRFVM_HPP_*/ diff --git a/src/coreComponents/linearAlgebra/unitTests/testHypredrive.cpp b/src/coreComponents/linearAlgebra/unitTests/testHypredrive.cpp index 30e981e0224..6ca3daa065b 100644 --- a/src/coreComponents/linearAlgebra/unitTests/testHypredrive.cpp +++ b/src/coreComponents/linearAlgebra/unitTests/testHypredrive.cpp @@ -13,6 +13,8 @@ #ifdef GEOS_USE_HYPRE #include "linearAlgebra/interfaces/hypre/HypreInterface.hpp" +#include "linearAlgebra/interfaces/hypre/HypreMGR.hpp" +#include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFracturesALM.hpp" #include "linearAlgebra/unitTests/testLinearAlgebraUtils.hpp" #ifdef GEOS_USE_HYPREDRV @@ -250,13 +252,25 @@ TEST( HypredriveYaml, BuildsGeneratedYamlForEveryMGRStrategy ) ++value ) { StrategyType const strategy = static_cast< StrategyType >( value ); + stdVector< string > strategyFieldNames = fieldNames; + array1d< int > strategyNumComponentsPerField = numComponentsPerField; + if( strategy == StrategyType::singlePhasePoromechanicsConformingFracturesALM ) + { + strategyFieldNames = { "field0", "field1", "field2" }; + strategyNumComponentsPerField.resize( 3 ); + } + else if( strategy == StrategyType::singlePhasePoromechanicsConformingFracturesALMReservoirFVM ) + { + strategyFieldNames = { "field0", "field1", "field2", "field3" }; + strategyNumComponentsPerField.resize( 4 ); + } hypre::hypredrive::InputArgsParseTarget target; try { EXPECT_TRUE( hypre::hypredrive::buildInputArgsParseTarget( makeMgrParameters( strategy ), - fieldNames, - numComponentsPerField, + strategyFieldNames, + strategyNumComponentsPerField, target ) ) << static_cast< int >( value ); EXPECT_EQ( target.source, hypre::hypredrive::InputSource::generatedFallback ) << static_cast< int >( value ); @@ -275,6 +289,79 @@ 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: [totalDisplacement_0, totalDisplacement_1, totalDisplacement_2, " + "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 ); + EXPECT_NE( target.argument.find( "aggressive:" ), std::string::npos ); + EXPECT_NE( target.argument.find( "max_nnz_row: 10" ), 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( HypreMGR, SetsUpFullyCoupledSinglePhaseALM ) +{ + array1d< int > numComponentsPerField( 3 ); + numComponentsPerField[0] = 3; + numComponentsPerField[1] = 3; + numComponentsPerField[2] = 1; + + HypreMatrix matrix; + testing::computeIdentity( MPI_COMM_GEOS, 7, matrix ); + + LinearSolverParameters params = makeMgrParameters( LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsConformingFracturesALM ); + HyprePrecWrapper precond; + ASSERT_EQ( HYPRE_MGRCreate( &precond.ptr ), 0 ); + HypreMGRData mgrData; + mgrData.pointMarkers.resize( 7 ); + for( HYPRE_Int label = 0; label < 7; ++label ) + { + mgrData.pointMarkers[label] = label; + } + + hypre::mgr::SinglePhasePoromechanicsConformingFracturesALM strategy( numComponentsPerField.toView() ); + strategy.setup( params.mgr, precond, mgrData ); + + EXPECT_EQ( HYPRE_MGRSetup( precond.ptr, matrix.unwrapped(), nullptr, nullptr ), 0 ); + + HypreVector rhs; + HypreVector solution; + rhs.create( matrix.numLocalRows(), MPI_COMM_GEOS ); + rhs.set( 1.0 ); + solution.create( matrix.numLocalCols(), MPI_COMM_GEOS ); + solution.zero(); + EXPECT_EQ( HYPRE_MGRSolve( precond.ptr, matrix.unwrapped(), rhs.unwrapped(), solution.unwrapped() ), 0 ); + + EXPECT_EQ( HYPRE_MGRDestroy( precond.ptr ), 0 ); + EXPECT_EQ( mgrData.coarseSolver.destroy( mgrData.coarseSolver.ptr ), 0 ); + EXPECT_EQ( mgrData.nestedSolver.destroy( mgrData.nestedSolver.ptr ), 0 ); +} + size_t countSubstrings( std::string const & text, std::string const & token ) { size_t count = 0; diff --git a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp index d9c6a954a54..3bad3f9d633 100644 --- a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp +++ b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp @@ -295,42 +295,47 @@ 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 + singlePhasePoromechanicsConformingFracturesALMReservoirFVM, ///< single phase poromechanics with conforming fractures for ALM, with + ///< wells + 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 +587,8 @@ ENUM_STRINGS( LinearSolverParameters::MGR::StrategyType, "hybridSinglePhasePoromechanics", "singlePhasePoromechanicsEmbeddedFractures", "singlePhasePoromechanicsConformingFractures", + "singlePhasePoromechanicsConformingFracturesALM", + "singlePhasePoromechanicsConformingFracturesALMReservoirFVM", "singlePhasePoromechanicsReservoirFVM", "thermalSinglePhasePoromechanicsReservoirFVM", "compositionalMultiphaseFVM", diff --git a/src/coreComponents/linearAlgebra/utilities/SparsityPatternUtilities.hpp b/src/coreComponents/linearAlgebra/utilities/SparsityPatternUtilities.hpp new file mode 100644 index 00000000000..e20c26d9b5e --- /dev/null +++ b/src/coreComponents/linearAlgebra/utilities/SparsityPatternUtilities.hpp @@ -0,0 +1,56 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 SparsityPatternUtilities.hpp + */ + +#ifndef GEOS_LINEARALGEBRA_UTILITIES_SPARSITYPATTERNUTILITIES_HPP_ +#define GEOS_LINEARALGEBRA_UTILITIES_SPARSITYPATTERNUTILITIES_HPP_ + +#include "common/DataTypes.hpp" + +namespace geos +{ + +/** + * @brief Append all entries from one sparsity pattern to another. + * + * Coupled solvers build a diagonal pattern from the DofManager and then widen + * it with their coupling entries. Keeping this small operation in one place + * avoids subtly different row-copy loops in every such solver. + * + * Entries already present in @p target are discarded, so the result is the union + * of the two patterns. @p target must have been sized for at least that union: + * a SparsityPattern grows on demand, but growing it row by row is quadratic. + * + * @param target the pattern to append to + * @param source the pattern whose entries are appended + */ +inline void appendSparsityPattern( SparsityPattern< globalIndex > & target, + SparsityPattern< globalIndex > const & source ) +{ + GEOS_ERROR_IF_NE( target.numRows(), source.numRows() ); + GEOS_ERROR_IF_NE( target.numColumns(), source.numColumns() ); + for( localIndex row = 0; row < source.numRows(); ++row ) + { + globalIndex const * const columns = source.getColumns( row ).dataIfContiguous(); + target.insertNonZeros( row, columns, columns + source.numNonZeros( row ) ); + } +} + +} // namespace geos + +#endif // GEOS_LINEARALGEBRA_UTILITIES_SPARSITYPATTERNUTILITIES_HPP_ diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp index 274998b5734..c1f83e559ac 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp @@ -87,6 +87,8 @@ #endif #include +#include +#include #include @@ -1030,6 +1032,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 * @@ -1063,6 +1202,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 @@ -1098,6 +1249,22 @@ 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.clear(); + neighbor3DGlobalIds.insert( 0, coordinateNeighborGlobalIds.begin(), coordinateNeighborGlobalIds.end() ); + } + // 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/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index 5fd4f484ac6..235a129d64c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -1495,10 +1495,14 @@ void CompositionalMultiphaseFVM::assembleHydrofracFluxTerms( real64 const GEOS_U DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs, - CRSMatrixView< real64, localIndex const > const & dR_dAper ) + CRSMatrixView< real64, localIndex const > const & dR_dAper, + stdMap< string, localIndex > const * const dR_dAperOffsets ) { GEOS_MARK_FUNCTION; + GEOS_ERROR_IF( dR_dAperOffsets != nullptr, + "CompositionalMultiphaseFVM does not support mesh-specific dR/dAperture row offsets." ); + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( m_discretizationName ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp index b08b523a353..264b6fefec5 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp @@ -155,7 +155,8 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs, - CRSMatrixView< real64, localIndex const > const & dR_dAper ) override final; + CRSMatrixView< real64, localIndex const > const & dR_dAper, + stdMap< string, localIndex > const * const dR_dAperOffsets ) override final; virtual void updatePhaseMobility( ObjectManagerBase & dataGroup ) const override; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index 75afb4fcd65..e2b99c482ea 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -155,7 +155,13 @@ class FlowSolverBase : public PhysicsSolverBase * @param dofManager degree-of-freedom manager associated with the linear system * @param localMatrix the system matrix * @param localRhs the system right-hand side vector - * @param dR_dAper + * @param dR_dAper derivative of the flux residual with respect to the aperture + * @param dR_dAperOffsets first row of @p dR_dAper belonging to each mesh body, + * keyed by mesh body name. The caller and this solver walk their own + * mesh targets, so the mesh body name is the only key they are + * guaranteed to agree on; the caller must therefore target each body at + * a single discretization level. Pass nullptr when @p dR_dAper uses a + * single index space, i.e. when the caller has a single mesh target. */ virtual void assembleHydrofracFluxTerms( real64 const time_n, real64 const dt, @@ -163,9 +169,10 @@ class FlowSolverBase : public PhysicsSolverBase DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs, - CRSMatrixView< real64, localIndex const > const & dR_dAper ) + CRSMatrixView< real64, localIndex const > const & dR_dAper, + stdMap< string, localIndex > const * const dR_dAperOffsets ) { - GEOS_UNUSED_VAR ( time_n, dt, domain, dofManager, localMatrix, localRhs, dR_dAper ); + GEOS_UNUSED_VAR ( time_n, dt, domain, dofManager, localMatrix, localRhs, dR_dAper, dR_dAperOffsets ); GEOS_ERROR( "Poroelastic fluxes with conforming fractures not yet implemented." ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp index c8c54388630..ad2a994a364 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp @@ -570,7 +570,8 @@ void SinglePhaseFVM< BASE >::assembleHydrofracFluxTerms( real64 const GEOS_UNUSE DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs, - CRSMatrixView< real64, localIndex const > const & dR_dAper ) + CRSMatrixView< real64, localIndex const > const & dR_dAper, + stdMap< string, localIndex > const * const dR_dAperOffsets ) { GEOS_MARK_FUNCTION; @@ -581,7 +582,7 @@ void SinglePhaseFVM< BASE >::assembleHydrofracFluxTerms( real64 const GEOS_UNUSE string const & dofKey = dofManager.getKey( SinglePhaseBase::viewKeyStruct::elemDofFieldString() ); - this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshName, MeshLevel const & mesh, string_array const & ) { @@ -617,6 +618,21 @@ void SinglePhaseFVM< BASE >::assembleHydrofracFluxTerms( real64 const GEOS_UNUSE fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( auto & stencil ) { + if( stencil.size() == 0 ) + { + return; + } + localIndex const dR_dAperOffset = [&]() + { + if( dR_dAperOffsets == nullptr ) + { + return localIndex( 0 ); + } + auto const offsetIt = dR_dAperOffsets->find( meshName ); + GEOS_ERROR_IF( offsetIt == dR_dAperOffsets->end(), + GEOS_FMT( "No dR/dAperture row offset is available for mesh body '{}'", meshName ) ); + return offsetIt->second; + }(); typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); if( m_isThermal ) @@ -630,7 +646,8 @@ void SinglePhaseFVM< BASE >::assembleHydrofracFluxTerms( real64 const GEOS_UNUSE dt, localMatrix.toViewConstSizes(), localRhs.toView(), - dR_dAper ); + dR_dAper, + dR_dAperOffset ); } else { @@ -643,7 +660,8 @@ void SinglePhaseFVM< BASE >::assembleHydrofracFluxTerms( real64 const GEOS_UNUSE dt, localMatrix.toViewConstSizes(), localRhs.toView(), - dR_dAper ); + dR_dAper, + dR_dAperOffset ); } } ); } ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp index a85545c0c43..a0adca9dec5 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp @@ -178,7 +178,8 @@ class SinglePhaseFVM : public BASE DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs, - CRSMatrixView< real64, localIndex const > const & dR_dAper ) override final; + CRSMatrixView< real64, localIndex const > const & dR_dAper, + stdMap< string, localIndex > const * const dR_dAperOffsets ) override final; /**@}*/ diff --git a/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp b/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp index 11a8ed584e4..01c82ea079c 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.cpp @@ -443,9 +443,10 @@ assembleHydrofracFluxTerms( real64 const time_n, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs, - CRSMatrixView< real64, localIndex const > const & dR_dAper ) + CRSMatrixView< real64, localIndex const > const & dR_dAper, + stdMap< string, localIndex > const * const dR_dAperOffsets ) { - flowSolver()->assembleHydrofracFluxTerms( time_n, dt, domain, dofManager, localMatrix, localRhs, dR_dAper ); + flowSolver()->assembleHydrofracFluxTerms( time_n, dt, domain, dofManager, localMatrix, localRhs, dR_dAper, dR_dAperOffsets ); } template< typename RESERVOIR_SOLVER > diff --git a/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.hpp b/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.hpp index 4a460846dba..c7d557b4db4 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.hpp @@ -103,7 +103,8 @@ class CompositionalMultiphaseReservoirAndWells : public CoupledReservoirAndWells DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs, - CRSMatrixView< real64, localIndex const > const & dR_dAper ); + CRSMatrixView< real64, localIndex const > const & dR_dAper, + stdMap< string, localIndex > const * const dR_dAperOffsets ); template< typename SUBREGION_TYPE > void accumulationAssemblyLaunch( DofManager const & dofManager, diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp index 63b6eee82b0..f4619d013e1 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp @@ -22,6 +22,7 @@ #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_COUPLEDRESERVOIRANDWELLSBASE_HPP_ #include "physicsSolvers/multiphysics/CoupledSolver.hpp" +#include "linearAlgebra/utilities/SparsityPatternUtilities.hpp" #include "common/TimingMacros.hpp" #include "constitutive/permeability/PermeabilityFields.hpp" @@ -136,11 +137,7 @@ class CoupledReservoirAndWellsBase : public CoupledSolver< RESERVOIR_SOLVER, WEL pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data()); // Copy the original nonzeros - for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) - { - globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous(); - pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow )); - } + appendSparsityPattern( pattern, patternDiag ); // Add the nonzeros from coupling addCouplingSparsityPattern( domain, dofManager, pattern.toView()); diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp index e26c8abaf34..f6a580505a7 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp @@ -23,6 +23,7 @@ #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" #include "constitutive/fluid/singlefluid/SingleFluidFields.hpp" #include "physicsSolvers/multiphysics/HydrofractureSolverKernels.hpp" +#include "linearAlgebra/utilities/SparsityPatternUtilities.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp" #include "physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" @@ -491,7 +492,7 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::setupCoupling( DomainPartition template< typename POROMECHANICS_SOLVER > void HydrofractureSolver< POROMECHANICS_SOLVER >::setSparsityPattern( DomainPartition & domain, DofManager & dofManager, - CRSMatrix< real64, globalIndex > & localMatrix, + CRSMatrix< real64, globalIndex > & GEOS_UNUSED_PARAM( localMatrix ), SparsityPattern< globalIndex > & pattern ) { SparsityPattern< globalIndex > patternOriginal; @@ -513,16 +514,12 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::setSparsityPattern( DomainPart rowLengths.data()); // Copy the original nonzeros - for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) - { - globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); - pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow )); - } + appendSparsityPattern( pattern, patternOriginal ); // Add the nonzeros from coupling addFluxApertureCouplingSparsityPattern( domain, dofManager, pattern.toView()); - setUpDflux_dApertureMatrix( domain, dofManager, localMatrix ); + setUpDflux_dApertureMatrix( domain ); } template< typename POROMECHANICS_SOLVER > @@ -706,13 +703,26 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::assembleSystem( real64 const t localRhs ); } + if( !getRefDerivativeFluxResidual_dAperture() ) + { + setUpDflux_dApertureMatrix( domain ); + } + // Move without touching: zero() memsets the entries in this space, and + // touching here would mark the immutable sparsity structure dirty on device. + getRefDerivativeFluxResidual_dAperture()->move( parallelDeviceMemorySpace, false ); + getRefDerivativeFluxResidual_dAperture()->zero(); flowSolver()->assembleHydrofracFluxTerms( time, dt, domain, dofManager, localMatrix, localRhs, - getDerivativeFluxResidual_dNormalJump() ); + getDerivativeFluxResidual_dNormalJump(), + nullptr ); + + // Read-only on the host from here on: do not touch, or the next iteration + // has to re-upload the whole matrix. + getRefDerivativeFluxResidual_dAperture()->move( hostMemorySpace, false ); assembleForceResidualDerivativeWrtPressure( domain, localMatrix, localRhs ); @@ -720,7 +730,6 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::assembleSystem( real64 const t assembleFluidLeakSource( time, dt, domain, dofManager, localMatrix, localRhs ); - this->getRefDerivativeFluxResidual_dAperture()->zero(); } template< typename POROMECHANICS_SOLVER > @@ -1056,48 +1065,68 @@ real64 HydrofractureSolver< POROMECHANICS_SOLVER >::setNextDt( real64 const & cu return nextDt; } template< typename POROMECHANICS_SOLVER > -void HydrofractureSolver< POROMECHANICS_SOLVER >::setUpDflux_dApertureMatrix( DomainPartition & domain, - DofManager const & dofManager, - CRSMatrix< real64, globalIndex > & localMatrix ) +void HydrofractureSolver< POROMECHANICS_SOLVER >::setUpDflux_dApertureMatrix( DomainPartition & domain ) { std::unique_ptr< CRSMatrix< real64, localIndex > > & derivativeFluxResidual_dAperture = this->getRefDerivativeFluxResidual_dAperture(); + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( flowSolver()->getDiscretizationName() ); + + // Count the fracture rows and accumulate their capacities in one traversal. + stdVector< localIndex > rowCapacities; + localIndex numMeshTargets = 0; + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshName, + MeshLevel const & mesh, + string_array const & regionNames ) { - localIndex numRows = 0; - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel & mesh, - string_array const & regionNames ) + // The stencil sweeps here, the insertion sweep below and the flux kernel + // (called with a null offset map) all index dR/dAperture by the raw + // per-target surface element index. That is only unambiguous for a single + // target: a second one would alias into the first target's rows. + ++numMeshTargets; + GEOS_ERROR_IF_GT_MSG( numMeshTargets, 1, + GEOS_FMT( "{}: the hydrofracture solver supports a single mesh target; '{}' is the second.", + this->getName(), meshName ) ); + + localIndex numMeshRows = 0; + mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const, + FaceElementSubRegion const & elementSubRegion ) { - mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const, - FaceElementSubRegion const & elementSubRegion ) + numMeshRows += elementSubRegion.size(); + } ); + rowCapacities.resize( rowCapacities.size() + numMeshRows, 0 ); + + fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) { - numRows += elementSubRegion.size(); - } ); + localIndex const numFluxElems = stencil.stencilSize( iconn ); + typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); + + for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) + { + GEOS_ERROR_IF_GE_MSG( sei[iconn][k0], + LvArray::integerConversion< localIndex >( rowCapacities.size() ), + "Surface stencil index exceeds the fracture derivative matrix size." ); + rowCapacities[sei[iconn][k0]] += numFluxElems; + } + } } ); + } ); - derivativeFluxResidual_dAperture = std::make_unique< CRSMatrix< real64, localIndex > >( numRows, numRows ); - derivativeFluxResidual_dAperture->setName( this->getName() + "/derivativeFluxResidual_dAperture" ); + localIndex const numRows = LvArray::integerConversion< localIndex >( rowCapacities.size() ); + derivativeFluxResidual_dAperture = std::make_unique< CRSMatrix< real64, localIndex > >( numRows, numRows ); + 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; - } - // TODO This is way too much. The With the full system rowSize is not a good estimate for this. - for( localIndex row = 0; row < numRows; ++row ) - { - derivativeFluxResidual_dAperture->reserveNonZeros( row, maxRowSize ); - } + if( numRows > 0 ) + { + derivativeFluxResidual_dAperture->resizeFromRowCapacities< parallelHostPolicy >( numRows, + numRows, + rowCapacities.data() ); } - string const presDofKey = dofManager.getKey( SinglePhaseBase::viewKeyStruct::elemDofFieldString() ); - - NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); - FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); - FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( flowSolver()->getDiscretizationName() ); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel const & mesh, string_array const & ) @@ -1108,9 +1137,10 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::setUpDflux_dApertureMatrix( Do { localIndex const numFluxElems = stencil.stencilSize( iconn ); typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); - for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) { + GEOS_ERROR_IF_GE_MSG( sei[iconn][k0], numRows, + "Surface stencil index exceeds the fracture derivative matrix size." ); for( localIndex k1 = 0; k1 < numFluxElems; ++k1 ) { derivativeFluxResidual_dAperture->insertNonZero( sei[iconn][k0], sei[iconn][k1], 0.0 ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp index 02e2e70d830..cf16fbe3665 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp @@ -212,9 +212,7 @@ class HydrofractureSolver : public POROMECHANICS_SOLVER SparsityPatternView< globalIndex > const & pattern ) const; - void setUpDflux_dApertureMatrix( DomainPartition & domain, - DofManager const & dofManager, - CRSMatrix< real64, globalIndex > & localMatrix ); + void setUpDflux_dApertureMatrix( DomainPartition & domain ); virtual void setMGRStrategy() override; diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp index 7218768bc96..95cf15b0861 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp @@ -27,12 +27,14 @@ #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/solidMechanics/contact/ContactFields.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp" +#include "physicsSolvers/multiphysics/PoromechanicsSolver.hpp" #include "constitutive/solid/CoupledSolidBase.hpp" #include "constitutive/contact/HydraulicApertureBase.hpp" #include "constitutive/contact/HydraulicApertureRelationSelector.hpp" #include "finiteVolume/FluxApproximationBase.hpp" #include "common/DataTypes.hpp" #include "mesh/DomainPartition.hpp" +#include "linearAlgebra/utilities/SparsityPatternUtilities.hpp" namespace geos { @@ -86,16 +88,12 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, rowLengths.data()); // Copy the original nonzeros - for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) - { - globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); - pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow )); - } + appendSparsityPattern( pattern, patternOriginal ); // Add the nonzeros from coupling addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView()); - setUpDflux_dApertureMatrix( domain, dofManager, localMatrix ); + setUpDflux_dApertureMatrix( domain ); } virtual void assembleSystem( real64 const time_n, @@ -110,6 +108,16 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, this->solidMechanicsSolver()->synchronizeFractureState( domain ); + // The flux assembly accumulates into this matrix. Clear it before every + // Newton assembly and make the host copy explicit before the host-side + // coupling kernels consume it. + if( !m_derivativeFluxResidual_dAperture ) + { + setUpDflux_dApertureMatrix( domain ); + } + m_derivativeFluxResidual_dAperture->move( parallelDeviceMemorySpace, false ); + m_derivativeFluxResidual_dAperture->zero(); + assembleElementBasedContributions( time_n, dt, domain, @@ -124,7 +132,10 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, dofManager, localMatrix, localRhs, - getDerivativeFluxResidual_dNormalJump() ); + getDerivativeFluxResidual_dNormalJump(), + nullptr ); + + m_derivativeFluxResidual_dAperture->move( hostMemorySpace, false ); // This step must occur after the fluxes are assembled because that's when DerivativeFluxResidual_dAperture is filled. assembleCouplingTerms( time_n, @@ -335,54 +346,76 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, * @brief Set up the Dflux_dApertureMatrix object * * @param domain - * @param dofManager - * @param localMatrix */ - void setUpDflux_dApertureMatrix( DomainPartition & domain, - DofManager const & GEOS_UNUSED_PARAM( dofManager ), - CRSMatrix< real64, globalIndex > & localMatrix ) + void setUpDflux_dApertureMatrix( DomainPartition & domain ) { integer const numComp = numFluidComponents(); + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); - this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + localIndex numMeshTargets = 0; + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshName, MeshLevel const & mesh, string_array const & regionNames ) { std::unique_ptr< CRSMatrix< real64, localIndex > > & derivativeFluxResidual_dAperture = getRefDerivativeFluxResidual_dAperture(); + // The matrix is re-created per target and the flux kernel indexes it by + // the raw per-target surface element index, so only the last target would + // survive and the others would write into its rows. + ++numMeshTargets; + GEOS_ERROR_IF_GT_MSG( numMeshTargets, 1, + GEOS_FMT( "{}: this solver supports a single mesh target; '{}' is the second.", + this->getName(), meshName ) ); + + localIndex numRows = 0; + localIndex numCol = 0; { // 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; + numCol = numRows; // number of rows (equations) = number of fracture elements * number of components numRows *= numComp; 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; - } - // TODO This is way too much. The With the full system rowSize is not a good estimate for this. - for( localIndex row = 0; row < numRows; ++row ) + // array1d's sized constructor value-initializes, so no explicit zero(). + array1d< localIndex > rowCapacities( numRows ); + fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) { - derivativeFluxResidual_dAperture->reserveNonZeros( row, maxRowSize ); + localIndex const numFluxElems = stencil.stencilSize( iconn ); + typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); + + for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) + { + // The stencil sweep covers every SurfaceElementStencil on the mesh, + // while numRows only counts the subregions found in regionNames. + GEOS_ERROR_IF_GE_MSG( sei[iconn][k0] * numComp + numComp - 1, numRows, + "Surface stencil index exceeds the fracture derivative matrix size." ); + for( integer ic = 0; ic < numComp; ic++ ) + { + rowCapacities[sei[iconn][k0] * numComp + ic] += numFluxElems; + } + } } - } + } ); - NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); - FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); - FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); + if( numRows > 0 ) + { + derivativeFluxResidual_dAperture->resizeFromRowCapacities< parallelHostPolicy >( numRows, + numCol, + rowCapacities.data() ); + } fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) { @@ -393,11 +426,15 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) { + GEOS_ERROR_IF_GE_MSG( sei[iconn][k0] * numComp + numComp - 1, numRows, + "Surface stencil index exceeds the fracture derivative matrix size." ); for( localIndex k1 = 0; k1 < numFluxElems; ++k1 ) { - for( integer ic = 0; ic < numComp; ic++ ) + for( integer ic = 0; ic < numComp; ++ic ) { - derivativeFluxResidual_dAperture->insertNonZero( sei[iconn][k0] * numComp + ic, sei[iconn][k1], 0.0 ); + derivativeFluxResidual_dAperture->insertNonZero( sei[iconn][k0] * numComp + ic, + sei[iconn][k1], + 0.0 ); } } } diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp index 468117b09dd..cc73dba16f0 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,48 @@ 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 + LinearSolverParameters const & linearSolverParameters = this->m_linearSolverParameters.get(); + if( !hypre::hypredrive::shouldUse( linearSolverParameters ) ) + { + return; + } + + this->solidMechanicsSolver()->forDiscretizationOnMeshTargets( + domain.getMeshBodies(), + [&]( string const &, MeshLevel const & mesh, string_array const & ) + { + NodeManager const & nodeManager = mesh.getNodeManager(); + arrayView1d< globalIndex const > const dofIndex = + nodeManager.getReference< array1d< globalIndex > >( + dofManager.getKey( solidMechanics::totalDisplacement::key() ) ); + array1d< ParallelVector > modes = + LAIHelperFunctions::computeRigidBodyModes< ParallelVector >( nodeManager.referencePosition(), + dofIndex, + dofManager.rankOffset(), + dofManager.numLocalDofs() ); + + // Each target's modes are supported on that target's displacement dofs and + // vanish elsewhere, so they are appended rather than summed: the near-null + // space of N mechanically independent bodies is 6N-dimensional, and summing + // into six vectors cannot span it. computeRigidBodyModes already returns + // them normalized, so no further scaling is needed. + for( ParallelVector & mode : modes ) + { + m_linearSolverNearNullKernel.emplace_back( std::move( mode ) ); + } + } ); +#else + GEOS_UNUSED_VAR( domain, dofManager ); +#endif +} + template< typename FLOW_SOLVER, typename MECHANICS_SOLVER > void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::initializePostInitialConditionsPreSubGroups() { @@ -202,7 +249,6 @@ void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::assembleElementB // step 1: apply the full poromechanics coupling on the target regions on the poromechanics solver set< string > poromechanicsRegionNames; - this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, string_array const & regionNames ) 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 7407f593662..19962b85313 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp @@ -25,6 +25,7 @@ #include "constitutive/contact/HydraulicApertureBase.hpp" #include "constitutive/contact/HydraulicApertureRelationSelector.hpp" #include "finiteVolume/FluxApproximationBase.hpp" +#include "linearAlgebra/utilities/SparsityPatternUtilities.hpp" #include "mesh/SurfaceElementRegion.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp" @@ -63,103 +64,84 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupCouplin template< typename FLOW_SOLVER > -void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupSystem( DomainPartition & domain, - DofManager & dofManager, - CRSMatrix< real64, globalIndex > & localMatrix, - ParallelVector & rhs, - ParallelVector & solution, - bool const setSparsity ) +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + CRSMatrix< real64, globalIndex > & localMatrix, + SparsityPattern< globalIndex > & pattern ) { GEOS_MARK_FUNCTION; - 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 + // These must be called before assembling the contact-dependent pattern. 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 ); + GEOS_ERROR_IF_NE( flowPattern.numRows(), mechanicsPattern.numRows() ); + GEOS_ERROR_IF_NE( flowPattern.numColumns(), mechanicsPattern.numColumns() ); + + // Count the union of the two sorted row sets. This avoids double-counting + // shared diagonal entries without assuming either pattern contains the + // other. + array1d< localIndex > rowLengths( flowPattern.numRows() ); + rowLengths.zero(); + for( localIndex localRow = 0; localRow < flowPattern.numRows(); ++localRow ) { - // Step 1: Get the flow sparsity pattern (includes DofManager diagonal blocks + flux stencil connections) - SparsityPattern< globalIndex > flowPattern; - this->flowSolver()->setSparsityPattern( domain, dofManager, localMatrix, flowPattern ); - - // Step 2: Get the solid mechanics sparsity pattern (includes DofManager diagonal blocks + mechanics coupling) - // Note: setSparsityPattern replaces its output pattern, so we use a separate variable. - SparsityPattern< globalIndex > mechPattern; - this->solidMechanicsSolver()->setSparsityPattern( domain, dofManager, localMatrix, mechPattern ); - - // Step 3: Compute combined row lengths (overestimate due to shared DofManager diagonal entries) - array1d< localIndex > rowLengths( flowPattern.numRows() ); - for( localIndex localRow = 0; localRow < flowPattern.numRows(); ++localRow ) + arraySlice1d< globalIndex const > const flowColumns = flowPattern.getColumns( localRow ); + arraySlice1d< globalIndex const > const mechanicsColumns = mechanicsPattern.getColumns( localRow ); + localIndex flowColumn = 0; + localIndex mechanicsColumn = 0; + while( flowColumn < flowColumns.size() || mechanicsColumn < mechanicsColumns.size() ) { - rowLengths[localRow] = flowPattern.numNonZeros( localRow ) + mechPattern.numNonZeros( localRow ); - } - - // Step 4: Add the number of nonzeros induced by flow-mechanics coupling - addTransmissibilityCouplingNNZ( domain, dofManager, rowLengths.toView() ); - addPressureForceCouplingNNZ( domain, dofManager, rowLengths.toView() ); - addMatrixPressureBubbleCouplingNNZ( domain, dofManager, rowLengths.toView() ); - - // Step 5: Create a new pattern with enough capacity for the full coupled matrix - SparsityPattern< globalIndex > pattern; - pattern.resizeFromRowCapacities< parallelHostPolicy >( flowPattern.numRows(), - flowPattern.numColumns(), - rowLengths.data() ); - - // Step 6: Copy flow pattern entries - for( localIndex localRow = 0; localRow < flowPattern.numRows(); ++localRow ) - { - globalIndex const * cols = flowPattern.getColumns( localRow ).dataIfContiguous(); - pattern.insertNonZeros( localRow, cols, cols + flowPattern.numNonZeros( localRow ) ); - } - - // Step 7: Copy mechanics pattern entries (duplicates with flow diagonal blocks handled by insertNonZeros) - for( localIndex localRow = 0; localRow < mechPattern.numRows(); ++localRow ) - { - globalIndex const * cols = mechPattern.getColumns( localRow ).dataIfContiguous(); - pattern.insertNonZeros( localRow, cols, cols + mechPattern.numNonZeros( localRow ) ); + if( mechanicsColumn == mechanicsColumns.size() || + ( flowColumn < flowColumns.size() && flowColumns[flowColumn] < mechanicsColumns[mechanicsColumn] ) ) + { + ++flowColumn; + } + else if( flowColumn == flowColumns.size() || mechanicsColumns[mechanicsColumn] < flowColumns[flowColumn] ) + { + ++mechanicsColumn; + } + else + { + ++flowColumn; + ++mechanicsColumn; + } + ++rowLengths[localRow]; } - - // Step 8: Add the flow-mechanics coupling patterns - addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView() ); - addPressureForceCouplingPattern( domain, dofManager, pattern.toView() ); - addMatrixPressureBubbleCouplingPattern( domain, dofManager, pattern.toView() ); - - // Set up the derivative flux residual matrix - setUpDflux_dApertureMatrix( domain, dofManager, localMatrix ); - - // Assimilate the sparsity pattern into the local matrix - localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) ); } - localMatrix.setName( this->getName() + "/matrix" ); - - rhs.setName( this->getName() + "/rhs" ); - rhs.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); - - solution.setName( this->getName() + "/solution" ); - solution.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); - - if( !this->m_precond && this->m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) - { - this->m_precond = this->createPreconditioner( domain ); - } + // Add the number of nonzeros induced by coupling + addTransmissibilityCouplingNNZ( domain, dofManager, rowLengths.toView() ); + addPressureForceCouplingNNZ( domain, dofManager, rowLengths.toView() ); + addMatrixPressureBubbleCouplingNNZ( domain, dofManager, rowLengths.toView() ); + + // Allocate the coupled pattern in one pass. Growing an already populated + // pattern row by row would shift every subsequent row of the contiguous + // column buffer on each call, which is quadratic in the number of rows. + pattern.resizeFromRowCapacities< parallelHostPolicy >( flowPattern.numRows(), + flowPattern.numColumns(), + rowLengths.data() ); + + // Copy both subsolver patterns in. insertNonZeros discards entries that are + // already present, so what remains is the union counted above. + appendSparsityPattern( pattern, flowPattern ); + appendSparsityPattern( pattern, mechanicsPattern ); + + // Add the nonzeros from coupling + addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView() ); + addPressureForceCouplingPattern( domain, dofManager, pattern.toView() ); + addMatrixPressureBubbleCouplingPattern( domain, dofManager, pattern.toView() ); + + // Set up the derivative flux residual matrix + setUpDflux_dApertureMatrix( domain ); } template< typename FLOW_SOLVER > @@ -175,6 +157,18 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleSyst // Synchronize fracture state this->solidMechanicsSolver()->synchronizeFractureState( domain ); + // setSparsityPattern owns the contact lookup tables and this matrix. Rather + // than lazily rebuilding a second copy of that prologue here, require that it + // has run: setupSystem( ..., setSparsity = false ) is not supported. + GEOS_ERROR_IF( !m_derivativeFluxResidual_dAperture, + GEOS_FMT( "{}: setupSystem must be called with sparsity construction enabled before assembling.", + this->getName() ) ); + + // Move without touching: zero() below memsets the entries in this space, and + // touching here would mark the immutable sparsity structure dirty on device. + m_derivativeFluxResidual_dAperture->move( parallelDeviceMemorySpace, false ); + m_derivativeFluxResidual_dAperture->zero(); + // Assemble element-based contributions (mechanics + flow accumulation) assembleElementBasedContributions( time_n, dt, domain, dofManager, localMatrix, localRhs ); @@ -185,10 +179,23 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleSyst dofManager, localMatrix, localRhs, - getDerivativeFluxResidual_dNormalJump() ); + getDerivativeFluxResidual_dNormalJump(), + &m_derivativeFluxResidual_dApertureOffsets ); + + // The flux kernel populates the derivative matrix in device memory. The + // coupling assembly below reads it on the host, so bring it over without + // touching it: touching would force a re-upload on the next assembly. + m_derivativeFluxResidual_dAperture->move( hostMemorySpace, false ); // Assemble coupling terms (must be after flux assembly) assembleCouplingTerms( time_n, dt, domain, dofManager, localMatrix, localRhs ); + + if constexpr ( hasWells ) + { + this->flowSolver()->wellSolver()->assembleSystem( time_n, dt, domain, dofManager, localMatrix, localRhs ); + this->flowSolver()->assembleCouplingTerms( time_n, dt, domain, dofManager, localMatrix, localRhs ); + } + } template< typename FLOW_SOLVER > @@ -325,50 +332,102 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::updateState( template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: -setUpDflux_dApertureMatrix( DomainPartition & domain, - DofManager const & dofManager, - CRSMatrix< real64, globalIndex > & localMatrix ) +setUpDflux_dApertureMatrix( DomainPartition & domain ) { - GEOS_UNUSED_VAR( dofManager ); + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); + string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); - this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + // Build the global row offsets and the row capacities together, so that each + // target is visited only once before the matrix is allocated. + m_derivativeFluxResidual_dApertureOffsets.clear(); + stdVector< localIndex > rowCapacities; + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshName, MeshLevel const & mesh, string_array const & regionNames ) { - std::unique_ptr< CRSMatrix< real64, localIndex > > & derivativeFluxResidual_dAperture = getRefDerivativeFluxResidual_dAperture(); + GEOS_UNUSED_VAR( regionNames ); + ElementRegionManager const & elemManager = mesh.getElemManager(); + // These offsets are consumed by the flow sub-solver, which walks its own + // mesh targets and therefore resolves the discretization level with its own + // discretization name. The mesh body name is the only part of a target both + // solvers are guaranteed to agree on, so it is the key; that in turn + // requires each body to appear exactly once here. + GEOS_ERROR_IF( m_derivativeFluxResidual_dApertureOffsets.find( meshName ) != + m_derivativeFluxResidual_dApertureOffsets.end(), + GEOS_FMT( "{}: mesh body '{}' is targeted at more than one discretization level. The augmented " + "Lagrangian contact formulation supports a single level per mesh body.", + this->getName(), meshName ) ); + + localIndex const rowOffset = rowCapacities.size(); + m_derivativeFluxResidual_dApertureOffsets.get_inserted( meshName ) = rowOffset; + + // The stencil sweeps below index rows by the raw surface-element index, so + // the contact fracture must be the only face-element region on this target: + // a second one would alias into its rows. Embedded-surface regions hold a + // different subregion type and contribute no SurfaceElementStencil here, so + // they are left alone. The region is required rather than optional because + // every consumer of this matrix (assembleCouplingTerms, + // assembleFluidMassResidualDerivativeWrtDisplacement) looks it up + // unconditionally on every target; skipping a target here would also leave + // its offset pointing at the next target's rows. + localIndex numFractureRegions = 0; + elemManager.forElementRegions< SurfaceElementRegion >( [&]( SurfaceElementRegion const & region ) { - // 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 ) + if( region.subRegionType() == SurfaceElementRegion::SurfaceSubRegionType::faceElement ) { - localIndex const rowSize = localMatrix.numNonZeros( row ); - maxRowSize = maxRowSize > rowSize ? maxRowSize : rowSize; + ++numFractureRegions; } + } ); + GEOS_ERROR_IF_NE_MSG( numFractureRegions, 1, + GEOS_FMT( "{}: mesh target '{}' holds {} face-element regions. The augmented Lagrangian " + "contact formulation requires exactly one, named '{}'.", + this->getName(), meshName, numFractureRegions, fractureRegionName ) ); + GEOS_ERROR_IF( !elemManager.hasRegion( fractureRegionName ), + GEOS_FMT( "{}: mesh target '{}' does not hold the fracture region '{}' of the contact solver.", + this->getName(), meshName, fractureRegionName ) ); + + SurfaceElementRegion const & fractureRegion = elemManager.getRegion< SurfaceElementRegion >( fractureRegionName ); + FaceElementSubRegion const & fractureSubRegion = fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); + rowCapacities.resize( rowOffset + fractureSubRegion.size(), 0 ); - for( localIndex row = 0; row < numRows; ++row ) + fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) { - derivativeFluxResidual_dAperture->reserveNonZeros( row, maxRowSize ); + localIndex const numFluxElems = stencil.stencilSize( iconn ); + typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); + for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) + { + localIndex const row = rowOffset + sei[iconn][k0]; + GEOS_ERROR_IF_GE_MSG( row, + LvArray::integerConversion< localIndex >( rowCapacities.size() ), + "Surface stencil index exceeds the fracture derivative matrix size." ); + rowCapacities[row] += numFluxElems; + } } - } + } ); + } ); - NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); - FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); - FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); + std::unique_ptr< CRSMatrix< real64, localIndex > > & derivativeFluxResidual_dAperture = getRefDerivativeFluxResidual_dAperture(); + localIndex const numRows = rowCapacities.size(); + derivativeFluxResidual_dAperture = std::make_unique< CRSMatrix< real64, localIndex > >( numRows, numRows ); + derivativeFluxResidual_dAperture->setName( this->getName() + "/derivativeFluxResidual_dAperture" ); + if( numRows > 0 ) + { + derivativeFluxResidual_dAperture->resizeFromRowCapacities< parallelHostPolicy >( numRows, + numRows, + rowCapacities.data() ); + } + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshName, + MeshLevel const & mesh, + string_array const & regionNames ) + { + GEOS_UNUSED_VAR( regionNames ); + localIndex const rowOffset = m_derivativeFluxResidual_dApertureOffsets.at( meshName ); fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) { @@ -379,9 +438,13 @@ setUpDflux_dApertureMatrix( DomainPartition & domain, for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) { + localIndex const row = rowOffset + sei[iconn][k0]; + GEOS_ERROR_IF_GE_MSG( row, numRows, "Surface stencil index exceeds the fracture derivative matrix size." ); for( localIndex k1 = 0; k1 < numFluxElems; ++k1 ) { - derivativeFluxResidual_dAperture->insertNonZero( sei[iconn][k0], sei[iconn][k1], 0.0 ); + derivativeFluxResidual_dAperture->insertNonZero( row, + rowOffset + sei[iconn][k1], + 0.0 ); } } } @@ -800,8 +863,14 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + // assembleSystem has already brought this matrix to the host after the flux + // assembly; the traversal below only reads it. CRSMatrixView< real64 const, localIndex const > const & dFluxResidual_dNormalJump = getDerivativeFluxResidual_dNormalJump().toViewConst(); + auto const derivativeOffsetIt = m_derivativeFluxResidual_dApertureOffsets.find( meshName ); + GEOS_ERROR_IF( derivativeOffsetIt == m_derivativeFluxResidual_dApertureOffsets.end(), + GEOS_FMT( "No dR/dAperture row offset is available for mesh body '{}'", meshName ) ); + localIndex const derivativeOffset = derivativeOffsetIt->second; string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); string const & bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); @@ -826,22 +895,24 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, localIndex const numElems = subRegion.size(); - // Allocate temporary storage for aperture derivatives (computed via kernels) + // These arrays are temporary workspaces for one assembly call. Keeping them + // local avoids retaining potentially large aperture buffers for the solver's + // lifetime while still giving the device kernels explicit views to capture. 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 ); + arrayView2d< real64 > const dAperturedUView = dAperturedU.toView(); + arrayView2d< real64 > const dAperturedBView = dAperturedB.toView(); // 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 ) + [&, dAperturedUView, dAperturedBView] ( string const &, + finiteElement::FiniteElementBase const & subRegionFE, + arrayView1d< localIndex const > const & faceElementList ) { poromechanicsALMKernels::ComputeApertureDerivativesFactory kernelFactory( dispDofNumber, @@ -851,8 +922,8 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, localRhs, 0.0, // dt not used faceElementList, - dAperturedU.toView(), - dAperturedB.toView() ); + dAperturedUView, + dAperturedBView ); real64 maxResidual = finiteElement:: interfaceBasedKernelApplication @@ -940,14 +1011,14 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, } // 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 ); + localIndex const numColumns = dFluxResidual_dNormalJump.numNonZeros( derivativeOffset + kfe ); + arraySlice1d< localIndex const > const & columns = dFluxResidual_dNormalJump.getColumns( derivativeOffset + kfe ); + arraySlice1d< real64 const > const & values = dFluxResidual_dNormalJump.getEntries( derivativeOffset + kfe ); for( localIndex kfe1 = 0; kfe1 < numColumns; ++kfe1 ) { real64 const dR_dAper = values[kfe1]; - localIndex const kfe2 = columns[kfe1]; + localIndex const kfe2 = columns[kfe1] - derivativeOffset; bool const isOpen = ( fractureState[kfe2] == FractureState::Open ); if( !isOpen && !isFractureOpen ) @@ -1028,7 +1099,7 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, for( localIndex kfe1 = 0; kfe1 < numColumns; ++kfe1 ) { real64 const dR_dAper = values[kfe1]; - localIndex const kfe2 = columns[kfe1]; + localIndex const kfe2 = columns[kfe1] - derivativeOffset; bool const isOpen = ( fractureState[kfe2] == FractureState::Open ); if( !isOpen && !isFractureOpen ) @@ -1092,9 +1163,6 @@ addMatrixPressureBubbleCouplingNNZ( DomainPartition const & domain, 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 ); @@ -1153,9 +1221,6 @@ addMatrixPressureBubbleCouplingPattern( DomainPartition const & domain, 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 ); @@ -1202,6 +1267,9 @@ assembleMatrixPressureBubbleContribution( real64 const dt, using namespace contact; + string const flowDofKey = dofManager.getKey( m_pressureKey ); + string const mechanicsDiscretizationName = this->solidMechanicsSolver()->getDiscretizationName(); + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, string_array const & regionNames ) @@ -1212,8 +1280,6 @@ assembleMatrixPressureBubbleContribution( real64 const dt, 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 ); @@ -1250,7 +1316,7 @@ assembleMatrixPressureBubbleContribution( real64 const dt, constitutive::PorousSolidBase, CellElementSubRegion >( mesh, poromechanicsRegionNames, - this->solidMechanicsSolver()->getDiscretizationName(), + mechanicsDiscretizationName, FlowSolverBase::viewKeyStruct::solidNamesString(), kernelFactory ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp index 865bfd66f1a..010ceb64ec4 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -38,6 +38,14 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec using Base::m_rhs; using Base::m_solution; + /// True when the flow solver carries well degrees of freedom. + static constexpr bool hasWells = std::is_same_v< FLOW_SOLVER, SinglePhaseReservoirAndWells<> >; + + static_assert( hasWells || std::is_same_v< FLOW_SOLVER, SinglePhaseBase >, + "SinglePhasePoromechanicsConformingFracturesALM supports only the SinglePhaseBase and " + "SinglePhaseReservoirAndWells<> flow solvers. Both setMGRStrategy and assembleSystem branch " + "on hasWells, so a new instantiation must be handled in both places." ); + /// String used to form the solverName used to register solvers in CoupledSolver static string coupledSolverAttributePrefix() { return "poromechanicsConformingFracturesALM"; } @@ -85,12 +93,10 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec virtual void setupCoupling( DomainPartition const & domain, DofManager & dofManager ) const override final; - virtual void setupSystem( DomainPartition & domain, - DofManager & dofManager, - CRSMatrix< real64, globalIndex > & localMatrix, - ParallelVector & rhs, - ParallelVector & solution, - bool const setSparsity = true ) override final; + virtual void setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + CRSMatrix< real64, globalIndex > & localMatrix, + SparsityPattern< globalIndex > & pattern ) override final; virtual void assembleSystem( real64 const time, real64 const dt, @@ -103,12 +109,54 @@ 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( linearSolverParameters.preconditionerType != LinearSolverParameters::PreconditionerType::mgr ) + { + return; + } + + // Wells contribute their own dof labels and need an extra reduction level + // to keep the well block out of the coarse grid, so they get a separate + // strategy. + if constexpr ( hasWells ) + { + linearSolverParameters.mgr.strategy = + LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsConformingFracturesALMReservoirFVM; + } + else + { + linearSolverParameters.mgr.strategy = + LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsConformingFracturesALM; + } + linearSolverParameters.mgr.separateComponents = true; + + GEOS_LOG_LEVEL_RANK_0( logInfo::LinearSolver, + GEOS_FMT( "{}: MGR strategy set to {}", this->getName(), + EnumStrings< LinearSolverParameters::MGR::StrategyType >::toString( linearSolverParameters.mgr.strategy ) ) ); } /**@}*/ +protected: + + virtual void initializePreSubGroups() override + { + Base::initializePreSubGroups(); + + // The ALM fracture assembly carries a single flow dof per fracture element: + // the dR/dAperture matrix is sized numElements x numElements and the contact + // kernels have no temperature block. Reject the thermal input rather than + // silently assembling an incomplete Jacobian. + // Checking the flow sub-solver too: PoromechanicsSolver only rejects the + // opposite direction (thermal coupled solver over a non-thermal flow + // solver), so a thermal SinglePhaseFVM under a non-thermal ALM solver would + // otherwise reach the two-equation thermal connector kernel. + GEOS_THROW_IF( this->m_isThermal || this->flowSolver()->isThermal(), + GEOS_FMT( "{}: thermal coupling is not supported by {}", + this->getName(), this->getCatalogName() ), + InputError, this->getDataContext() ); + } + private: struct viewKeyStruct : public Base::viewKeyStruct @@ -233,12 +281,8 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec * @brief Set up the Dflux_dApertureMatrix object * * @param domain - * @param dofManager - * @param localMatrix */ - void setUpDflux_dApertureMatrix( DomainPartition & domain, - DofManager const & dofManager, - CRSMatrix< real64, globalIndex > & localMatrix ); + void setUpDflux_dApertureMatrix( DomainPartition & domain ); std::unique_ptr< CRSMatrix< real64, localIndex > > & getRefDerivativeFluxResidual_dAperture() { @@ -257,6 +301,8 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec std::unique_ptr< CRSMatrix< real64, localIndex > > m_derivativeFluxResidual_dAperture; + stdMap< string, localIndex > m_derivativeFluxResidual_dApertureOffsets; + string const m_pressureKey = SinglePhaseBase::viewKeyStruct::elemDofFieldString(); }; diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp index 98da0ab2737..50d804eaf7f 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp @@ -92,8 +92,9 @@ class SinglePhaseReservoirAndWells : public CoupledReservoirAndWellsBase< RESERV DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs, - CRSMatrixView< real64, localIndex const > const & dR_dAper ) - { flowSolver()->assembleHydrofracFluxTerms( time_n, dt, domain, dofManager, localMatrix, localRhs, dR_dAper ); } + CRSMatrixView< real64, localIndex const > const & dR_dAper, + stdMap< string, localIndex > const * const dR_dAperOffsets ) + { flowSolver()->assembleHydrofracFluxTerms( time_n, dt, domain, dofManager, localMatrix, localRhs, dR_dAper, dR_dAperOffsets ); } template< typename SUBREGION_TYPE > void accumulationAssemblyLaunch( DofManager const & dofManager, diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFractures.hpp index e8ca8d6c3b7..58cbe419ff2 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFractures.hpp @@ -87,7 +87,8 @@ class ConnectorBasedAssemblyKernel : public singlePhaseFVMKernels::FluxComputeKe real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs, - CRSMatrixView< real64, localIndex const > const & dR_dAper ) + CRSMatrixView< real64, localIndex const > const & dR_dAper, + localIndex const dR_dAperOffset ) : Base( rankOffset, stencilWrapper, flowDofNumberAccessor, @@ -98,6 +99,7 @@ class ConnectorBasedAssemblyKernel : public singlePhaseFVMKernels::FluxComputeKe localMatrix, localRhs ), m_dR_dAper( dR_dAper ), + m_dR_dAperOffset( dR_dAperOffset ), m_dPerm_dDispJump( fracturePermeabilityAccessors.get( fields::permeability::dPerm_dDispJump {} ) ) {} @@ -147,7 +149,7 @@ class ConnectorBasedAssemblyKernel : public singlePhaseFVMKernels::FluxComputeKe { stack.dofColIndices[i * numDof + jdof] = offset + jdof; } - stack.localColIndices[ i ] = m_sei( iconn, i ); + stack.localColIndices[ i ] = m_dR_dAperOffset + m_sei( iconn, i ); } } @@ -250,7 +252,7 @@ class ConnectorBasedAssemblyKernel : public singlePhaseFVMKernels::FluxComputeKe localIndex const localRow ) { - localIndex const row = LvArray::integerConversion< localIndex >( m_sei( iconn, i ) ); + localIndex const row = m_dR_dAperOffset + LvArray::integerConversion< localIndex >( m_sei( iconn, i ) ); m_dR_dAper.addToRowBinarySearch< parallelDeviceAtomic >( row, stack.localColIndices.data(), @@ -264,6 +266,7 @@ class ConnectorBasedAssemblyKernel : public singlePhaseFVMKernels::FluxComputeKe private: CRSMatrixView< real64, localIndex const > m_dR_dAper; + localIndex const m_dR_dAperOffset; ElementViewConst< arrayView4d< real64 const > > const m_dPerm_dDispJump; }; @@ -300,7 +303,8 @@ class ConnectorBasedAssemblyKernelFactory real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs, - CRSMatrixView< real64, localIndex const > const & dR_dAper ) + CRSMatrixView< real64, localIndex const > const & dR_dAper, + localIndex const dR_dAperOffset = 0 ) { integer constexpr NUM_DOF = 1; // pressure integer constexpr NUM_EQN = 1; @@ -317,7 +321,7 @@ class ConnectorBasedAssemblyKernelFactory kernelType kernel( rankOffset, stencilWrapper, flowDofNumberAccessor, flowAccessors, fluidAccessors, permAccessors, fracPermAccessors, - dt, localMatrix, localRhs, dR_dAper ); + dt, localMatrix, localRhs, dR_dAper, dR_dAperOffset ); kernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); } diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFractures.hpp index 07f4d0775de..261da1eb2d3 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFractures.hpp @@ -100,7 +100,8 @@ class ConnectorBasedAssemblyKernel : public singlePhasePoromechanicsConformingFr real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs, - CRSMatrixView< real64, localIndex const > const & dR_dAper ) + CRSMatrixView< real64, localIndex const > const & dR_dAper, + localIndex const dR_dAperOffset ) : Base( rankOffset, stencilWrapper, flowDofNumberAccessor, @@ -111,7 +112,8 @@ class ConnectorBasedAssemblyKernel : public singlePhasePoromechanicsConformingFr dt, localMatrix, localRhs, - dR_dAper ), + dR_dAper, + dR_dAperOffset ), m_temp( thermalSinglePhaseFlowAccessors.get( fields::flow::temperature {} ) ), m_enthalpy( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::enthalpy {} ) ), m_dEnthalpy( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy {} ) ), @@ -351,7 +353,8 @@ class ConnectorBasedAssemblyKernelFactory real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs, - CRSMatrixView< real64, localIndex const > const & dR_dAper ) + CRSMatrixView< real64, localIndex const > const & dR_dAper, + localIndex const dR_dAperOffset = 0 ) { integer constexpr NUM_DOF = 2; // pressure + temperature integer constexpr NUM_EQN = 2; // mass balance + energy balance @@ -376,7 +379,7 @@ class ConnectorBasedAssemblyKernelFactory flowDofNumberAccessor, flowAccessors, thermalFlowAccessors, fluidAccessors, thermalFluidAccessors, permAccessors, edfmPermAccessors, thermalConductivityAccessors, - dt, localMatrix, localRhs, dR_dAper ); + dt, localMatrix, localRhs, dR_dAper, dR_dAperOffset ); kernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); } diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index f0b90771230..dc0c204301a 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -20,6 +20,7 @@ #include "SolidMechanicsAugmentedLagrangianContact.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBase.hpp" +#include "linearAlgebra/utilities/SparsityPatternUtilities.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsConformingContactKernelsBase.hpp" @@ -370,11 +371,7 @@ void SolidMechanicsAugmentedLagrangianContact::setSparsityPattern( DomainPartiti pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data()); // Copy the original nonzeros - for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) - { - globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous(); - pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow )); - } + appendSparsityPattern( pattern, patternDiag ); // Add the nonzeros from coupling addCouplingSparsityPattern( domain, dofManager, pattern.toView());