diff --git a/inputFiles/poromechanicsFractures/Contact/AugmentedLagrangianMultipliers/ThermoPoroElastic_ALM_conforming_base.xml b/inputFiles/poromechanicsFractures/Contact/AugmentedLagrangianMultipliers/ThermoPoroElastic_ALM_conforming_base.xml new file mode 100644 index 00000000000..b4345570747 --- /dev/null +++ b/inputFiles/poromechanicsFractures/Contact/AugmentedLagrangianMultipliers/ThermoPoroElastic_ALM_conforming_base.xml @@ -0,0 +1,200 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/poromechanicsFractures/Contact/AugmentedLagrangianMultipliers/ThermoPoroElastic_ALM_conforming_smoke.xml b/inputFiles/poromechanicsFractures/Contact/AugmentedLagrangianMultipliers/ThermoPoroElastic_ALM_conforming_smoke.xml new file mode 100644 index 00000000000..0feaf82ea81 --- /dev/null +++ b/inputFiles/poromechanicsFractures/Contact/AugmentedLagrangianMultipliers/ThermoPoroElastic_ALM_conforming_smoke.xml @@ -0,0 +1,54 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index e2b99c482ea..1a1fb371bf7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -176,6 +176,18 @@ class FlowSolverBase : public PhysicsSolverBase GEOS_ERROR( "Poroelastic fluxes with conforming fractures not yet implemented." ); } + virtual void assembleHydrofracFluxTermsALM( real64 const time_n, + real64 const dt, + DomainPartition const & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) + { + GEOS_UNUSED_VAR ( time_n, dt, domain, dofManager, localMatrix, localRhs, dR_dAper ); + GEOS_ERROR( "Poroelastic fluxes with conforming fractures ALM not yet implemented." ); + } + void initializeState( DomainPartition & domain ); virtual void initializeFluidState( MeshLevel & mesh, string_array const & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index 4d01e0a5370..ac0d4d51852 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -268,6 +268,12 @@ class SinglePhaseBase : public FlowSolverBase virtual void updateState ( DomainPartition & domain ) override; + /** + * @brief Getter for the number of fluid components (species) + * @return the number of components + */ + integer numFluidComponents() const { return 1; } + /** * @brief Function to update all constitutive state and dependent variables * @param subRegion subregion that contains the fields diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp index ad2a994a364..b2da00aa2cc 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp @@ -45,6 +45,8 @@ #include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEmbeddedFractures.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFractures.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFractures.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALM.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFracturesALM.hpp" /** * @namespace the geos namespace that encapsulates the majority of the code @@ -669,6 +671,92 @@ void SinglePhaseFVM< BASE >::assembleHydrofracFluxTerms( real64 const GEOS_UNUSE } +template< typename BASE > +void SinglePhaseFVM< BASE >::assembleHydrofracFluxTermsALM( real64 const GEOS_UNUSED_PARAM ( time_n ), + real64 const dt, + DomainPartition const & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) +{ + GEOS_MARK_FUNCTION; + + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( m_discretizationName ); + + string const & dofKey = dofManager.getKey( SinglePhaseBase::viewKeyStruct::elemDofFieldString() ); + + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & ) + { + fluxApprox.forStencils< CellElementStencilTPFA, FaceElementToCellStencil >( mesh, [&]( auto & stencil ) + { + typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); + + if( m_isThermal ) + { + thermalSinglePhaseFVMKernels:: + FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(), + dofKey, + this->getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView() ); + } + else + { + singlePhaseFVMKernels:: + FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(), + dofKey, + this->getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView() ); + } + } ); + + fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( auto & stencil ) + { + typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); + + if( m_isThermal ) + { + thermalSinglePhasePoromechanicsConformingFracturesALMKernels:: + ConnectorBasedAssemblyKernelFactory::createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(), + dofKey, + this->getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView(), + dR_dAper ); + } + else + { + singlePhasePoromechanicsConformingFracturesALMKernels:: + ConnectorBasedAssemblyKernelFactory::createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(), + dofKey, + this->getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView(), + dR_dAper ); + } + } ); + } ); +} + template< typename BASE > void SinglePhaseFVM< BASE >::applyBoundaryConditions( real64 const time_n, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp index a0adca9dec5..d42ed1da040 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp @@ -181,6 +181,15 @@ class SinglePhaseFVM : public BASE CRSMatrixView< real64, localIndex const > const & dR_dAper, stdMap< string, localIndex > const * const dR_dAperOffsets ) override final; + virtual void + assembleHydrofracFluxTermsALM( real64 const time_n, + real64 const dt, + DomainPartition const & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) override final; + /**@}*/ virtual void diff --git a/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt b/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt index 36e8a23c9c0..684ab99233a 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt @@ -45,6 +45,7 @@ set( multiPhysicsSolvers_headers poromechanicsKernels/SinglePhasePoromechanics.hpp poromechanicsKernels/SinglePhasePoromechanics_impl.hpp poromechanicsKernels/SinglePhasePoromechanicsConformingFractures.hpp + poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALM.hpp poromechanicsKernels/SinglePhasePoromechanicsDamage.hpp poromechanicsKernels/SinglePhasePoromechanicsDamage_impl.hpp poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp @@ -58,6 +59,7 @@ set( multiPhysicsSolvers_headers poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFractures.hpp + poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFracturesALM.hpp poromechanicsKernels/ThermalSinglePhasePoromechanicsEmbeddedFractures.hpp SinglePhasePoromechanics.hpp SinglePhasePoromechanicsEmbeddedFractures.hpp diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp index de5de39883a..bf06083e779 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp @@ -101,11 +101,7 @@ class MultiphasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHANI protected: - virtual void setMGRStrategy() override - { - if( this->m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::mgr ) - GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}", this->getName(), this->getCatalogName())); - } + GEOS_MGR_STRATEGY_NOT_SUPPORTED() /** * @brief Helper function to recompute the bulk density diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp index 0ef1cb73f36..bb69c154637 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp @@ -83,7 +83,116 @@ assembleSystem( real64 const time_n, template< typename FLOW_SOLVER > void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >:: -assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, +assembleForceResidualDerivativeWrtPressure( string const & GEOS_UNUSED_PARAM( meshName ), + MeshLevel const & mesh, + string_array const & regionNames, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + GEOS_MARK_FUNCTION; + + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + EdgeManager const & edgeManager = mesh.getEdgeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + ArrayOfArraysView< localIndex const > const faceToEdgeMap = faceManager.edgeList().toViewConst(); + arrayView2d< localIndex const > const & edgeToNodeMap = edgeManager.nodeList().toViewConst(); + arrayView2d< real64 const > faceCenters = faceManager.faceCenter(); + arrayView2d< real64 const > const & faceNormal = faceManager.faceNormal(); + arrayView1d< real64 const > faceAreas = faceManager.faceArea(); + + string const & dispDofKey = dofManager.getKey( fields::solidMechanics::totalDisplacement::key() ); + string const & flowDofKey = dofManager.getKey( this->getFlowDofKey() ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + globalIndex const rankOffset = dofManager.rankOffset(); + + // Get the coordinates for all nodes + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition = nodeManager.referencePosition(); + + elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, + [&]( localIndex const, + FaceElementSubRegion const & subRegion ) + { + arrayView1d< globalIndex const > const & + flowDofNumber = subRegion.getReference< globalIndex_array >( flowDofKey ); + arrayView1d< real64 const > const & pressure = subRegion.getReference< array1d< real64 > >( fields::flow::pressure::key() ); + arrayView2d< localIndex const > const & elemsToFaces = subRegion.faceList().toViewConst(); + + forAll< serialPolicy >( subRegion.size(), [=, this]( localIndex const kfe ) + { + localIndex const kf0 = elemsToFaces[kfe][0]; + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 ); + + real64 Nbar[3]; + Nbar[ 0 ] = faceNormal[elemsToFaces[kfe][0]][0] - faceNormal[elemsToFaces[kfe][1]][0]; + Nbar[ 1 ] = faceNormal[elemsToFaces[kfe][0]][1] - faceNormal[elemsToFaces[kfe][1]][1]; + Nbar[ 2 ] = faceNormal[elemsToFaces[kfe][0]][2] - faceNormal[elemsToFaces[kfe][1]][2]; + LvArray::tensorOps::normalize< 3 >( Nbar ); + globalIndex rowDOF[3 * m_maxFaceNodes]; // this needs to be changed when dealing with arbitrary element types + real64 nodeRHS[3 * m_maxFaceNodes]; + stackArray1d< real64, 3 * m_maxFaceNodes > dRdP( 3*m_maxFaceNodes ); + globalIndex colDOF[1]; + colDOF[0] = flowDofNumber[kfe]; // pressure is always first + + for( localIndex kf=0; kf<2; ++kf ) + { + localIndex const faceIndex = elemsToFaces[kfe][kf]; + + // Compute local area contribution for each node + stackArray1d< real64, FaceManager::maxFaceNodes() > nodalArea; + this->solidMechanicsSolver()->computeFaceNodalArea( elemsToFaces[kfe][kf], + nodePosition, + faceToNodeMap, + faceToEdgeMap, + edgeToNodeMap, + faceCenters, + faceNormal, + faceAreas, + nodalArea ); + for( localIndex a=0; a( globalNodalForce, Nbar, nodalForceMag ); + + for( localIndex i=0; i<3; ++i ) + { + rowDOF[3*a+i] = dispDofNumber[faceToNodeMap( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i ); + // Opposite sign w.r.t. theory because of minus sign in stiffness matrix definition (K < 0) + nodeRHS[3*a+i] = +globalNodalForce[i] * pow( -1, kf ); + + // Opposite sign w.r.t. theory because of minus sign in stiffness matrix definition (K < 0) + dRdP( 3*a+i ) = -nodalArea[a] * Nbar[i] * pow( -1, kf ); + } + } + + for( localIndex idof = 0; idof < numNodesPerFace * 3; ++idof ) + { + localIndex const localRow = LvArray::integerConversion< localIndex >( rowDOF[idof] - rankOffset ); + + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + localMatrix.addToRow< parallelHostAtomic >( localRow, + colDOF, + &dRdP[idof], + 1 ); + RAJA::atomicAdd( parallelHostAtomic{}, &localRhs[localRow], nodeRHS[idof] ); + } + } + } + } ); + } ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >:: +assembleFluidMassResidualDerivativeWrtDisplacement( string const & GEOS_UNUSED_PARAM( meshName ), + MeshLevel const & mesh, string_array const & regionNames, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.hpp index 43052453f85..bfb2d2d37ad 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.hpp @@ -74,6 +74,8 @@ class MultiphasePoromechanicsConformingFractures : public PoromechanicsConformin */ string getCatalogName() const override { return catalogName(); } + GEOS_MGR_STRATEGY_NOT_SUPPORTED()//TODO should we keep ? + void assembleSystem( real64 const time_n, real64 const dt, DomainPartition & domain, @@ -85,14 +87,20 @@ class MultiphasePoromechanicsConformingFractures : public PoromechanicsConformin virtual void initializePreSubGroups() override; - virtual void assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, + virtual void assembleForceResidualDerivativeWrtPressure( string const & meshName, + MeshLevel const & mesh, + string_array const & regionNames, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) override final; + + virtual void assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, + MeshLevel const & mesh, string_array const & regionNames, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) override; - virtual integer numFluidComponents() const override { return this->flowSolver()->numFluidComponents(); } - virtual string getFlowDofKey() const override { return CompositionalMultiphaseBase::viewKeyStruct::elemDofFieldString(); } }; diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFracturesALM.cpp index ef17e087d1d..8c86d8b2ed8 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFracturesALM.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFracturesALM.cpp @@ -19,7 +19,6 @@ #include "MultiphasePoromechanicsConformingFracturesALM.hpp" - namespace geos { @@ -29,151 +28,56 @@ using namespace fields; template< typename FLOW_SOLVER > MultiphasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::MultiphasePoromechanicsConformingFracturesALM( const string & name, - Group * const parent ) + Group * const parent ) : Base( name, parent ) {} template< typename FLOW_SOLVER > -void MultiphasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupCoupling( DomainPartition const & domain, - DofManager & dofManager ) const -{ - GEOS_MARK_FUNCTION; - - GEOS_UNUSED_VAR( domain, dofManager ); - - GEOS_ERROR( "MultiphasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); - -} - - -template< typename FLOW_SOLVER > -void MultiphasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupSystem( DomainPartition & domain, - DofManager & dofManager, - CRSMatrix< real64, globalIndex > & localMatrix, - ParallelVector & rhs, - ParallelVector & solution, - bool const setSparsity ) +void MultiphasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + CRSMatrix< real64, globalIndex > & localMatrix, + SparsityPattern< globalIndex > & pattern ) { - GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain, dofManager, localMatrix, rhs, solution, setSparsity ); - - GEOS_ERROR( "MultiphasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); - -} - -template< typename FLOW_SOLVER > -void MultiphasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleSystem( real64 const time_n, - real64 const dt, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) -{ - - GEOS_MARK_FUNCTION; - - GEOS_UNUSED_VAR( time_n, dt, domain, dofManager, localMatrix, localRhs ); - - GEOS_ERROR( "MultiphasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); -} - -template< typename FLOW_SOLVER > -void MultiphasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleElementBasedContributions( real64 const time_n, - real64 const dt, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) -{ - GEOS_MARK_FUNCTION; - - GEOS_UNUSED_VAR( time_n, dt, domain, dofManager, localMatrix, localRhs ); - - GEOS_ERROR( "MultiphasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); - -} - -template< typename FLOW_SOLVER > -void MultiphasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleCouplingTerms( real64 const time_n, - real64 const dt, - DomainPartition const & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) -{ - GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain, dofManager, localMatrix, localRhs, time_n, dt ); - - GEOS_ERROR( "MultiphasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); -} - -template< typename FLOW_SOLVER > -void MultiphasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::updateState( DomainPartition & domain ) -{ - GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain ); - GEOS_ERROR( "MultiphasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); -} - - -template< typename FLOW_SOLVER > -void MultiphasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: -addTransmissibilityCouplingNNZ( DomainPartition const & domain, - DofManager const & dofManager, - arrayView1d< localIndex > const & rowLengths ) const -{ - GEOS_MARK_FUNCTION; - - GEOS_UNUSED_VAR( domain, dofManager, rowLengths ); - GEOS_ERROR( "MultiphasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + // Recompute fracture face/element geometry and rebuild the ALM contact solver's internal lists. + // These must happen before assembling the contact-dependent pattern; setSparsityPattern() is the + // only place this coupled solver hooks into for that, since it does not route through the contact + // sub-solver's own setupSystem(). + this->solidMechanicsSolver()->updateFractureGeometry( domain ); + Base::setSparsityPattern( domain, dofManager, localMatrix, pattern ); } template< typename FLOW_SOLVER > void MultiphasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: -addTransmissibilityCouplingPattern( DomainPartition const & domain, - DofManager const & dofManager, - SparsityPatternView< globalIndex > const & pattern ) const +assembleForceResidualDerivativeWrtPressure( string const & GEOS_UNUSED_PARAM( meshName ), + MeshLevel const & GEOS_UNUSED_PARAM( mesh ), + string_array const & GEOS_UNUSED_PARAM( regionNames ), + DofManager const & GEOS_UNUSED_PARAM( dofManager ), + CRSMatrixView< real64, globalIndex const > const & GEOS_UNUSED_PARAM( localMatrix ), + arrayView1d< real64 > const & GEOS_UNUSED_PARAM( localRhs ) ) { - GEOS_MARK_FUNCTION; - - GEOS_UNUSED_VAR( domain, dofManager, pattern ); - GEOS_ERROR( "MultiphasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + // Only reached through the monolithic assembleSystem(), which is exercised by FullyImplicit + // coupling. The compositional-flow ALM force/pressure coupling kernels are not implemented yet; + // this solver is sequential-only for now (see postInputInitialization/couplingType checks upstream). + GEOS_ERROR( GEOS_FMT( "{}: FullyImplicit coupling is not supported by {}", this->getName(), this->getCatalogName() ) ); } template< typename FLOW_SOLVER > void MultiphasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: -assembleForceResidualDerivativeWrtPressure( string const & meshName, - MeshLevel const & mesh, - arrayView1d< string const > const & regionNames, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) -{ - GEOS_MARK_FUNCTION; - - GEOS_UNUSED_VAR( meshName, mesh, regionNames, dofManager, localMatrix, localRhs ); - GEOS_ERROR( "MultiphasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); -} - -template< typename FLOW_SOLVER > -void MultiphasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: -assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, - arrayView1d< string const > const & regionNames, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, +assembleFluidMassResidualDerivativeWrtDisplacement( string const & GEOS_UNUSED_PARAM( meshName ), + MeshLevel const & GEOS_UNUSED_PARAM( mesh ), + string_array const & GEOS_UNUSED_PARAM( regionNames ), + DofManager const & GEOS_UNUSED_PARAM( dofManager ), + CRSMatrixView< real64, globalIndex const > const & GEOS_UNUSED_PARAM( localMatrix ), arrayView1d< real64 > const & GEOS_UNUSED_PARAM( localRhs ) ) { - GEOS_MARK_FUNCTION; - - GEOS_UNUSED_VAR( mesh, regionNames, dofManager, localMatrix ); - GEOS_ERROR( "MultiphasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); - + // See assembleForceResidualDerivativeWrtPressure: FullyImplicit-only path, not implemented for + // compositional flow yet. + GEOS_ERROR( GEOS_FMT( "{}: FullyImplicit coupling is not supported by {}", this->getName(), this->getCatalogName() ) ); } - template class MultiphasePoromechanicsConformingFracturesALM<>; template class MultiphasePoromechanicsConformingFracturesALM< CompositionalMultiphaseReservoirAndWells<> >; diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFracturesALM.hpp index 121b2916758..ab93d4d8238 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFracturesALM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFracturesALM.hpp @@ -20,6 +20,7 @@ #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_MULTIPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_ #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_MULTIPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_ +#include "physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp" #include "physicsSolvers/multiphysics/MultiphasePoromechanics.hpp" #include "physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp" @@ -27,16 +28,17 @@ namespace geos { template< typename FLOW_SOLVER = CompositionalMultiphaseBase > -class MultiphasePoromechanicsConformingFracturesALM : public MultiphasePoromechanics< FLOW_SOLVER, SolidMechanicsAugmentedLagrangianContact > +class MultiphasePoromechanicsConformingFracturesALM : public PoromechanicsConformingFractures< MultiphasePoromechanics, FLOW_SOLVER, SolidMechanicsAugmentedLagrangianContact > { public: - using Base = MultiphasePoromechanics< FLOW_SOLVER, SolidMechanicsAugmentedLagrangianContact >; + using Base = PoromechanicsConformingFractures< MultiphasePoromechanics, FLOW_SOLVER, SolidMechanicsAugmentedLagrangianContact >; using Base::m_solvers; using Base::m_dofManager; using Base::m_localMatrix; using Base::m_rhs; using Base::m_solution; + using Base::m_maxFaceNodes; /// String used to form the solverName used to register solvers in CoupledSolver static string coupledSolverAttributePrefix() { return "poromechanicsConformingFracturesALM"; } @@ -82,99 +84,48 @@ class MultiphasePoromechanicsConformingFracturesALM : public MultiphasePoromecha */ /**@{*/ - virtual void setupCoupling( DomainPartition const & domain, - DofManager & dofManager ) const override final; + GEOS_MGR_STRATEGY_NOT_SUPPORTED()//TODO: no MGR strategy exists yet for multiphase ALM - 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, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) override final; + /**@}*/ - virtual void updateState( DomainPartition & domain ) override final; +protected: - virtual void setMGRStrategy() override final + virtual void initializePreSubGroups() override { - if( this->m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::mgr ) - GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}", this->getName(), this->getCatalogName())); + Base::initializePreSubGroups(); + + GEOS_THROW_IF( this->m_isThermal || this->flowSolver()->isThermal(), + GEOS_FMT( "{}: thermal coupling is not supported by {}", + this->getName(), this->getCatalogName() ), + InputError, this->getDataContext() ); } - /**@}*/ + virtual void assembleForceResidualDerivativeWrtPressure( string const & meshName, + MeshLevel const & mesh, + string_array const & regionNames, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) override final; + + virtual void assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, + MeshLevel const & mesh, + string_array const & regionNames, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) override final; + + virtual string getFlowDofKey() const override { return CompositionalMultiphaseBase::viewKeyStruct::elemDofFieldString(); } private: struct viewKeyStruct : public Base::viewKeyStruct {}; - static const localIndex m_maxFaceNodes=11; // Maximum number of nodes on a contact face - - /** - * @Brief assemble the element-based contributions - * @param time_n the current time - * @param dt the time step - * @param domain the physical domain object - * @param dofManager degree-of-freedom manager associated with the linear system - * @param localMatrix the local system matrix - * @param localRhs the local system right-hand side vector - */ - void assembleElementBasedContributions( real64 const time_n, - real64 const dt, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ); - - virtual void assembleCouplingTerms( real64 const time_n, - real64 const dt, - DomainPartition const & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) override final; - - void assembleForceResidualDerivativeWrtPressure( string const & meshName, - MeshLevel const & mesh, - arrayView1d< string const > const & regionNames, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ); - - void assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, - arrayView1d< string const > const & regionNames, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ); - - /** - * @Brief add the nnz induced by the flux-aperture coupling - * @param domain the physical domain object - * @param dofManager degree-of-freedom manager associated with the linear system - * @param rowLenghts the nnz in each row - */ - void addTransmissibilityCouplingNNZ( DomainPartition const & domain, - DofManager const & dofManager, - arrayView1d< localIndex > const & rowLengths ) const; - - /** - * @Brief add the sparsity pattern induced by the flux-aperture coupling - * @param domain the physical domain object - * @param dofManager degree-of-freedom manager associated with the linear system - * @param pattern the sparsity pattern - */ - void addTransmissibilityCouplingPattern( DomainPartition const & domain, - DofManager const & dofManager, - SparsityPatternView< globalIndex > const & pattern ) const; - - - string const m_pressureKey = CompositionalMultiphaseBase::viewKeyStruct::elemDofFieldString(); - }; } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp index 95cf15b0861..9522fb58b9a 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp @@ -21,7 +21,11 @@ #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSCONFORMINGFRACTURES_HPP_ #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSCONFORMINGFRACTURES_HPP_ +#include "common/logger/Logger.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" +#include "physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp" #include "physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp" +#include "physicsSolvers/solidMechanics/contact/ContactFields.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBase.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" @@ -39,11 +43,11 @@ namespace geos { -template< template< typename, typename > class POROMECHANICS_BASE, typename FLOW_SOLVER > -class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, SolidMechanicsLagrangeContact > +template< template< typename, typename > class POROMECHANICS_BASE, typename FLOW_SOLVER = SinglePhaseBase , typename CONTACT_SOLVER = SolidMechanicsLagrangeContact > +class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, CONTACT_SOLVER > { public: - using Base = POROMECHANICS_BASE< FLOW_SOLVER, SolidMechanicsLagrangeContact >; + using Base = POROMECHANICS_BASE< FLOW_SOLVER, CONTACT_SOLVER >; PoromechanicsConformingFractures( const string & name, dataRepository::Group * const parent ) @@ -54,13 +58,28 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, DofManager & dofManager ) const override { /// We need to add 2 coupling terms: - // 1. Poromechanical coupling in the bulk + // 1. Poromechanical coupling in the bulk (p<->disp) + GEOS_LOG_RANK_0("[debug] \t \t Coupling flow to displacement (call super from poromeca)\n"); Base::setupCoupling( domain, dofManager ); - // 2. Traction - pressure coupling in the fracture + + + if constexpr (CONTACT_SOLVER::hasContactStabilization) { + + GEOS_LOG_RANK_0("[debug] \t \t Coupling flow to bubble\n"); + // 2. Pressure - bubble displacement coupling in the fracture + dofManager.addCoupling( this->getFlowDofKey(), + fields::contact::totalBubbleDisplacement::key(), + DofManager::Connector::Elem ); + + } + else { + // 2. Traction - pressure coupling in the fracture + GEOS_LOG_RANK_0("[debug] \t \t Coupling flow to traction\n"); dofManager.addCoupling( this->getFlowDofKey(), fields::contact::traction::key(), DofManager::Connector::Elem ); + } } virtual void setSparsityPattern( DomainPartition & domain, @@ -72,15 +91,32 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, SparsityPattern< globalIndex > patternOriginal; this->flowSolver()->setSparsityPattern( domain, dofManager, localMatrix, patternOriginal ); + SparsityPattern< globalIndex > mechanicsPattern; + if constexpr (CONTACT_SOLVER::hasContactStabilization) { + this->solidMechanicsSolver()->setSparsityPattern( domain, dofManager, localMatrix, mechanicsPattern ); + } + // Get the original row lengths (diagonal blocks only) array1d< localIndex > rowLengths( patternOriginal.numRows()); for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) { rowLengths[localRow] = patternOriginal.numNonZeros( localRow ); + if constexpr ( CONTACT_SOLVER::hasContactStabilization ) + { + rowLengths[localRow] += mechanicsPattern.numNonZeros( localRow ); // simple sum, see note below + } } // Add the number of nonzeros induced by coupling + //displacement (and opt. bubble) to flow coupling addTransmissibilityCouplingNNZ( domain, dofManager, rowLengths.toView()); + if constexpr (CONTACT_SOLVER::hasContactStabilization) { + //bubble to displacement coupling + GEOS_LOG_RANK_0("[debug] \t \t ALM nnz sizing"); + addPressureForceCouplingNNZ( domain, dofManager, rowLengths.toView() ); + addMatrixPressureBubbleCouplingNNZ( domain, dofManager, rowLengths.toView() );//TODO should be brought by CONTACT::STABILIZATION + } + // Create a new pattern with enough capacity for coupled matrix pattern.resizeFromRowCapacities< parallelHostPolicy >( patternOriginal.numRows(), @@ -89,13 +125,258 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, // Copy the original nonzeros appendSparsityPattern( pattern, patternOriginal ); + //ALM appendSparsityPattern( pattern, flowPattern ); + if constexpr (CONTACT_SOLVER::hasContactStabilization) { + appendSparsityPattern( pattern, mechanicsPattern ); + } // Add the nonzeros from coupling + //displacement (and opt. bubble) to flow coupling addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView()); + if constexpr (CONTACT_SOLVER::hasContactStabilization) { + GEOS_LOG_RANK_0("[debug] \t \t ALM coupling patterns"); + addPressureForceCouplingPattern( domain, dofManager, pattern.toView() ); + addMatrixPressureBubbleCouplingPattern( domain, dofManager, pattern.toView() ); + } + GEOS_LOG_RANK_0("[debug] \t \t Setting Apperture derivatives"); setUpDflux_dApertureMatrix( domain ); + + +// for( localIndex row = 0; row < pattern.numRows(); ++row ) +// { +// GEOS_LOG_RANK( GEOS_FMT( "row {} nnz {} cols {}", row, pattern.numNonZeros( row ), +// stringutilities::join( pattern.getColumns( row ).begin(), pattern.getColumns( row ).end(), "," ) ) ); +// } + } + //Stabilization specific + //TODO see refacto with below +void addPressureForceCouplingNNZ( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< localIndex > const & rowLengths ) const +{ + GEOS_MARK_FUNCTION; + + integer const numComp = this->flowSolver()->numFluidComponents(); + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const bubbleDofKey = dofManager.getKey( fields::contact::totalBubbleDisplacement::key() ); + + arrayView1d< globalIndex const > const & + bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); + SurfaceElementRegion const & fractureRegion = + elemManager.getRegion< SurfaceElementRegion >( fractureRegionName ); + FaceElementSubRegion const & fractureSubRegion = + fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); + + arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst(); + + // For each fracture element, add NNZ for (displacement_row, pressure_col) and (bubble_row, pressure_col) + forAll< serialPolicy >( fractureSubRegion.size(), [=, &rowLengths] ( localIndex const kfe ) + { + // For bubble DOFs: add 1 pressure column per bubble DOF row + for( localIndex kf = 0; kf < 2; ++kf ) + { + localIndex const faceIndex = elem2dToFaces[kfe][kf]; + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowNumber = bubbleDofNumber[faceIndex] + i - rankOffset; + if( rowNumber >= 0 && rowNumber < rowLengths.size() ) + { + rowLengths[rowNumber] += numComp; // One pressure column + } + } + } + } ); + } ); +} + +void addPressureForceCouplingPattern( DomainPartition const & domain, + DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern ) const +{ + GEOS_MARK_FUNCTION; + + integer const numComp = this->flowSolver()->numFluidComponents(); + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const bubbleDofKey = dofManager.getKey( fields::contact::totalBubbleDisplacement::key() ); + string const flowDofKey = dofManager.getKey( this->getFlowDofKey() ); + + arrayView1d< globalIndex const > const & + bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); + SurfaceElementRegion const & fractureRegion = + elemManager.getRegion< SurfaceElementRegion >( fractureRegionName ); + FaceElementSubRegion const & fractureSubRegion = + fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); + + arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst(); + arrayView1d< globalIndex const > const & + flowDofNumber = fractureSubRegion.getReference< globalIndex_array >( flowDofKey ); + + // For each fracture element, add pattern for (displacement_row, pressure_col) and (bubble_row, pressure_col) + forAll< serialPolicy >( fractureSubRegion.size(), [=] ( localIndex const kfe ) + { + globalIndex const pressureColIndex = flowDofNumber[kfe]; + + // For bubble DOFs + for( localIndex kf = 0; kf < 2; ++kf ) + { + localIndex const faceIndex = elem2dToFaces[kfe][kf]; + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowIndex = bubbleDofNumber[faceIndex] + i - rankOffset; + if( rowIndex >= 0 && rowIndex < pattern.numRows() ) + { + for(integer ic = 0; ic < numComp; ++ic) + pattern.insertNonZero( rowIndex, pressureColIndex + ic ); + } + } + } + } ); + } ); +} + +void addMatrixPressureBubbleCouplingNNZ( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< localIndex > const & rowLengths ) const +{ + GEOS_MARK_FUNCTION; + + integer const numComp = this->flowSolver()->numFluidComponents(); + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & regionNames ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const bubbleDofKey = dofManager.getKey( fields::contact::totalBubbleDisplacement::key() ); + string const flowDofKey = dofManager.getKey( this->getFlowDofKey() ); + arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + // Loop over matrix cell regions that have bubbles + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, + [&]( localIndex const, CellElementSubRegion const & subRegion ) + { + arrayView1d< localIndex const > const bubbleElems = subRegion.bubbleElementsList(); + arrayView2d< localIndex const > const elemsToFaces = subRegion.faceElementsList(); + arrayView1d< globalIndex const > const pressureDofNumber = subRegion.getReference< array1d< globalIndex > >( flowDofKey ); + + forAll< serialPolicy >( bubbleElems.size(), [=, &rowLengths]( localIndex const kk ) + { + localIndex const k = bubbleElems[kk]; + localIndex const faceIndex = elemsToFaces[kk][0]; + + // (bubble_row, pressure_col): 1 pressure column for each of the 3 bubble DOFs + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowNumber = bubbleDofNumber[faceIndex] + i - rankOffset; + if( rowNumber >= 0 && rowNumber < rowLengths.size() ) + { + rowLengths[rowNumber] += numComp; // One pressure DOF from matrix cell + } + } + + // (pressure_row, bubble_col): the matrix cell pressure couples to its 3 bubble DOFs (A_pb) + globalIndex const pRow = pressureDofNumber[k] - rankOffset; + if( pRow >= 0 && pRow < rowLengths.size() ) + { + for( integer ic = 0; ic const & pattern ) const +{ + GEOS_MARK_FUNCTION; + + integer const numComp = this->flowSolver()->numFluidComponents(); + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & regionNames ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const bubbleDofKey = dofManager.getKey( fields::contact::totalBubbleDisplacement::key() ); + string const flowDofKey = dofManager.getKey( this->getFlowDofKey() ); + + arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + // Loop over matrix cell regions that have bubbles + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, + [&]( localIndex const, CellElementSubRegion const & subRegion ) + { + arrayView1d< localIndex const > const bubbleElems = subRegion.bubbleElementsList(); + arrayView2d< localIndex const > const elemsToFaces = subRegion.faceElementsList(); + arrayView1d< globalIndex const > const pressureDofNumber = subRegion.getReference< array1d< globalIndex > >( flowDofKey ); + + forAll< serialPolicy >( bubbleElems.size(), [=]( localIndex const kk ) + { + localIndex const k = bubbleElems[kk]; + localIndex const faceIndex = elemsToFaces[kk][0]; + globalIndex const pressureColIndex = pressureDofNumber[k]; + + // (bubble_row, pressure_col) : A_bp + for( localIndex i = 0; i < 3; ++i ) + { + globalIndex const rowIndex = bubbleDofNumber[faceIndex] + i - rankOffset; + if( rowIndex >= 0 && rowIndex < pattern.numRows() ) + { + for( integer ic = 0; ic < numComp; ++ic ) + pattern.insertNonZero( rowIndex, pressureColIndex + ic); + } + } + + // (pressure_row, bubble_col) : A_pb -- transpose location + globalIndex const pRow = pressureDofNumber[k] - rankOffset; + if( pRow >= 0 && pRow < pattern.numRows() ) + { + for( localIndex i = 0; i < 3; ++i ) + { + for( integer ic = 0; ic < numComp; ++ic ) + pattern.insertNonZero( pRow + ic, bubbleDofNumber[faceIndex] + i ); + } + } + } ); + } ); + } ); +} + + virtual void assembleSystem( real64 const time_n, real64 const dt, DomainPartition & domain, @@ -111,13 +392,16 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, // 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 ) { + GEOS_LOG_RANK_0("[debug] \t \t Assembling dFluxdAperture if not already"); setUpDflux_dApertureMatrix( domain ); } m_derivativeFluxResidual_dAperture->move( parallelDeviceMemorySpace, false ); m_derivativeFluxResidual_dAperture->zero(); + GEOS_LOG_RANK_0("[debug] \t \t Assembling element based"); assembleElementBasedContributions( time_n, dt, domain, @@ -126,6 +410,7 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, localRhs ); // Assemble fluxes 3D/2D and get dFluidResidualDAperture + GEOS_LOG_RANK_0("[debug] \t \t Assembling hydrofrac flux terms"); this->flowSolver()->assembleHydrofracFluxTerms( time_n, dt, domain, @@ -138,12 +423,21 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, m_derivativeFluxResidual_dAperture->move( hostMemorySpace, false ); // This step must occur after the fluxes are assembled because that's when DerivativeFluxResidual_dAperture is filled. + GEOS_LOG_RANK_0("[debug] \t \t Assembling coupling terms"); assembleCouplingTerms( time_n, dt, domain, dofManager, localMatrix, localRhs ); + if constexpr ( std::is_same_v< FLOW_SOLVER, SinglePhaseReservoirAndWells<> > ) + { + this->flowSolver()->wellSolver()->assembleSystem( time_n, dt, domain, dofManager, localMatrix, localRhs ); + this->flowSolver()->assembleCouplingTerms( time_n, dt, domain, dofManager, localMatrix, localRhs ); + } + + // this->getSystemMatrix().write( "alm_branch.mtx", LAIOutputFormat::MATRIX_MARKET ); + } virtual void updateState( DomainPartition & domain ) override @@ -178,7 +472,8 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, { GEOS_MARK_FUNCTION; - integer const numComp = numFluidComponents(); + + integer const numComp = this->flowSolver()->numFluidComponents(); this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, // meshBodyName, MeshLevel const & mesh, @@ -192,6 +487,12 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + + //TODO (jafranc) - remove once ALM-bubble is frame as a stab method - tmp runtime is fine as it is tmp + if(this->solidMechanicsSolver()->hasStabilization())//why is this not done in SolidMech ? + { + + GEOS_LOG_RANK_0("[debug] \t \t Fetching stabilization from meca\n"); FluxApproximationBase const & stabilizationMethod = fvManager.getFluxApproximation( this->solidMechanicsSolver()->getStabilizationName() ); stabilizationMethod.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) @@ -235,7 +536,59 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, } } } ); - } ); + } + //to decide -- reduce duplication -- can we have both ? + if constexpr (CONTACT_SOLVER::hasContactStabilization) { + + GEOS_LOG_RANK_0("[debug] \t \t Fetching fv for bubbles\n"); + + FluxApproximationBase const & fvMethod = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); + + fvMethod.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + for( localIndex iconn=0; iconn( sesri[iconn][0] ); + + ArrayOfArraysView< localIndex const > const elemsToNodes = elementSubRegion.nodeList().toViewConst(); + + arrayView1d< globalIndex const > const faceElementDofNumber = + elementSubRegion.getReference< array1d< globalIndex > >( flowDofKey ); + + for( localIndex k0=0; k0= 0 && rowNumber < rowLengths.size() ) + { + for( localIndex k1=0; k1forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel const & mesh, string_array const & ) { - FaceManager const & faceManager = mesh.getFaceManager(); - NodeManager const & nodeManager = mesh.getNodeManager(); - ElementRegionManager const & elemManager = mesh.getElemManager(); - string const dispDofKey = dofManager.getKey( fields::solidMechanics::totalDisplacement::key() ); - string const flowDofKey = dofManager.getKey( this->getFlowDofKey() ); + NodeManager const & nodeManager = mesh.getNodeManager(); + FaceManager const & faceManager = mesh.getFaceManager(); + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); - arrayView1d< globalIndex const > const & - dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); - ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + if constexpr (CONTACT_SOLVER::hasContactStabilization) + { + GEOS_LOG_RANK_0(GEOS_FMT( "[debug] \t \t Offset tables: {}({}) flow \n {}({}) ut \n {}({}) ub", + dofManager.rankOffset(this->getFlowDofKey()), dofManager.numLocalDofs(this->getFlowDofKey()), + dofManager.rankOffset(fields::solidMechanics::totalDisplacement::key()), dofManager.numLocalDofs(fields::solidMechanics::totalDisplacement::key()), + dofManager.rankOffset(fields::contact::totalBubbleDisplacement::key()), dofManager.numLocalDofs(fields::contact::totalBubbleDisplacement::key()) + )); + } + else + { + GEOS_LOG_RANK_0(GEOS_FMT( "[debug] \t \t Offset tables: {}({}) flow \n {}({}) ut", + dofManager.rankOffset(this->getFlowDofKey()), dofManager.numLocalDofs(this->getFlowDofKey()), + dofManager.rankOffset(fields::solidMechanics::totalDisplacement::key()), dofManager.numLocalDofs(fields::solidMechanics::totalDisplacement::key()) + )); + } + + + GEOS_LOG_RANK_0("[debug] \t \t Adding T-pattern [flow -> utot]\n"); + addTransmissibilityCouplingPattern( domain, mesh, dofManager, pattern, + nodeManager.getReference< globalIndex_array >( dofManager.getKey( fields::solidMechanics::totalDisplacement::key() ) ), + [&faceToNodeMap](localIndex const& faceIndex, localIndex const& a){ return faceToNodeMap(faceIndex,a); }, + [&faceToNodeMap](localIndex const& faceIndex){ return faceToNodeMap.sizeOfArray(faceIndex); } ); + + if constexpr (CONTACT_SOLVER::hasContactStabilization) + { + GEOS_LOG_RANK_0("[debug] \t \t Adding T-pattern [flow -> btot]\n"); + addTransmissibilityCouplingPattern( domain, mesh, dofManager, pattern, + faceManager.getReference< globalIndex_array >( dofManager.getKey( fields::contact::totalBubbleDisplacement::key() ) ), + [](localIndex const & faceIndex, localIndex const& GEOS_UNUSED_PARAM(a)){ return faceIndex; }, + [](localIndex const & GEOS_UNUSED_PARAM(faceIndex)){ return 1; } ); + } + + } ); + } + + template< typename NODE_INDEX_MAP, typename NNODE_PER_FACE > + void addTransmissibilityCouplingPattern( DomainPartition const & domain, + MeshLevel const & mesh, + DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern, + arrayView1d< globalIndex const > const & coupledDisplacementDofNumber, + NODE_INDEX_MAP && dofIndirectionCb, + NNODE_PER_FACE && numNodesPerFace + ) const + { // Get the finite volume method used to compute the stabilization NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); @@ -274,7 +665,7 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, FluxApproximationBase const & fvDiscretization = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); SurfaceElementRegion const & fractureRegion = - elemManager.getRegion< SurfaceElementRegion >( this->solidMechanicsSolver()->getUniqueFractureRegionName() ); + mesh.getElemManager().getRegion< SurfaceElementRegion >( this->solidMechanicsSolver()->getUniqueFractureRegionName() ); FaceElementSubRegion const & fractureSubRegion = fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); @@ -284,16 +675,21 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst(); arrayView1d< globalIndex const > const & - flowDofNumber = fractureSubRegion.getReference< globalIndex_array >( flowDofKey ); - + flowDofNumber = fractureSubRegion.getReference< globalIndex_array >( dofManager.getKey(this->getFlowDofKey()) ); + + //debug + // ArrayOfArraysView< localIndex const > const elemsToNodes = fractureSubRegion.nodeList().toViewConst(); + globalIndex const rankOffset = dofManager.rankOffset(); fvDiscretization.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) { - forAll< serialPolicy >( stencil.size(), [=] ( localIndex const iconn ) + forAll< serialPolicy >( stencil.size(), [=,this] ( localIndex const iconn ) { localIndex const numFluxElems = stencil.stencilSize( iconn ); + + // A fracture connector has to be an edge shared by two faces if( numFluxElems == 2 ) { @@ -312,9 +708,14 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, // Get fracture, face and region/subregion/element indices (for elements on both sides) localIndex const fractureIndex = sei[iconn][kf]; + // integer n1 = numNodesPerFace(elem2dToFaces[fractureIndex][0]), + // n2 = numNodesPerFace(elem2dToFaces[fractureIndex][1]); + // GEOS_LOG_RANK_0( GEOS_FMT( "[debug] elemsToNodes.size={} vs faceToNodeMap.sizeOfArray(face0)={} vs sizeOfArray(face1)={}", + // elemsToNodes[sei[iconn][1]].size(), n1,n2 + // ) ); // Get the number of nodes - localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elem2dToFaces[fractureIndex][0] ); + // localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elem2dToFaces[fractureIndex][0] ); // Loop over the two sides of each fracture element for( localIndex kf1 = 0; kf1 < 2; ++kf1 ) @@ -322,12 +723,12 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, localIndex const faceIndex = elem2dToFaces[fractureIndex][kf1]; // Save the list of DOF associated with nodes - for( localIndex a=0; a( i ); - for( integer ic = 0; ic < numComp; ic++ ) + globalIndex const colIndex = coupledDisplacementDofNumber[dofIndirectionCb( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i ); + for( integer ic = 0; ic < this->flowSolver()->numFluidComponents(); ic++ ) { pattern.insertNonZero( rowIndex + ic, colIndex ); } @@ -339,9 +740,10 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, } } ); } ); - } ); + } + /** * @brief Set up the Dflux_dApertureMatrix object * @@ -349,98 +751,132 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, */ void setUpDflux_dApertureMatrix( DomainPartition & domain ) { - integer const numComp = numFluidComponents(); + integer const numComp = this->flowSolver()->numFluidComponents(); + localIndex numCols = 0.;//number of outerloop pass (not considering innermost component loop) + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); - - localIndex numMeshTargets = 0; + + string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); + // 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 ) + MeshLevel const & mesh, + string_array const & regionNames ) + { + 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 ) { - 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; + if( region.subRegionType() == SurfaceElementRegion::SurfaceSubRegionType::faceElement ) { - // calculate number of fracture elements - mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, - [&]( localIndex const, FaceElementSubRegion const & subRegion ) - { - numRows += subRegion.size(); - } ); - // number of columns (derivatives) = number of fracture elements - 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" ); + ++numFractureRegions; } - - // array1d's sized constructor value-initializes, so no explicit zero(). - array1d< localIndex > rowCapacities( numRows ); - fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + } ); + 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() * numComp, 0 ); + numCols += fractureSubRegion.size(); + + fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) { - for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) + localIndex const numFluxElems = stencil.stencilSize( iconn ); + typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); + for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) { - 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, + for( integer ic = 0; ic < numComp; ++ic){ + localIndex const row = rowOffset + sei[iconn][k0] * numComp; + GEOS_ERROR_IF_GE_MSG( row, + LvArray::integerConversion< localIndex >( rowCapacities.size() ), "Surface stencil index exceeds the fracture derivative matrix size." ); - for( integer ic = 0; ic < numComp; ic++ ) - { - rowCapacities[sei[iconn][k0] * numComp + ic] += numFluxElems; - } + rowCapacities[ row + ic ] += numFluxElems; } } - } ); - - if( numRows > 0 ) - { - derivativeFluxResidual_dAperture->resizeFromRowCapacities< parallelHostPolicy >( numRows, - numCol, - rowCapacities.data() ); } + } ); + } ); + + //write real data in structure + std::unique_ptr< CRSMatrix< real64, localIndex > > & derivativeFluxResidual_dAperture = getRefDerivativeFluxResidual_dAperture(); + localIndex const numRows = rowCapacities.size(); + derivativeFluxResidual_dAperture = std::make_unique< CRSMatrix< real64, localIndex > >( numRows, numCols ); + derivativeFluxResidual_dAperture->setName( this->getName() + "/derivativeFluxResidual_dAperture" ); + if( numRows > 0 ) + { + derivativeFluxResidual_dAperture->resizeFromRowCapacities< parallelHostPolicy >( numRows, + numCols, + 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 ) + fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) { - for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) - { - localIndex const numFluxElems = stencil.stencilSize( iconn ); - typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); + localIndex const numFluxElems = stencil.stencilSize( iconn ); + typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); - for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) + for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) + { + for(integer ic = 0 ; icinsertNonZero( sei[iconn][k0] * numComp + ic, - sei[iconn][k1], - 0.0 ); - } - } + derivativeFluxResidual_dAperture->insertNonZero( row, + rowOffset/numComp + sei[iconn][k1], //as component-independent indexing + 0.0 ); } } - } ); + } + } } ); + } ); } void assembleElementBasedContributions( real64 const time_n, @@ -454,6 +890,7 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, /// 3. assemble Force Residual w.r.t. pressure and Flow mass residual w.r.t. displacement + GEOS_LOG_RANK_0("[debug] \t \t Assembling Poromecha"); Base::assembleElementBasedTerms( time_n, dt, domain, dofManager, localMatrix, localRhs ); // Flow accumulation for fractures @@ -464,11 +901,13 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const, FaceElementSubRegion const & subRegion ) { + GEOS_LOG_RANK_0("[debug] \t \t Assembling Flow accumulation on Fractures"); this->flowSolver()->accumulationAssemblyLaunch( dofManager, subRegion, localMatrix, localRhs ); } ); } ); - this->solidMechanicsSolver()->assembleContact( domain, dofManager, localMatrix, localRhs ); + GEOS_LOG_RANK_0("[debug] \t \t Assembling contact"); + this->solidMechanicsSolver()->assembleContact( time_n, dt, domain, dofManager, localMatrix, localRhs ); } virtual void assembleCouplingTerms( real64 const time_n, @@ -480,127 +919,150 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, { GEOS_UNUSED_VAR( time_n, dt ); // These 2 steps need to occur after the fluxes are assembled because that's when DerivativeFluxResidual_dAperture is filled. - this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshName, MeshLevel const & mesh, string_array const & regionNames ) { /// 3. assemble Force Residual w.r.t. pressure and Flow mass residual w.r.t. displacement - assembleForceResidualDerivativeWrtPressure( mesh, regionNames, dofManager, localMatrix, localRhs ); - assembleFluidMassResidualDerivativeWrtDisplacement( mesh, regionNames, dofManager, localMatrix, localRhs ); + assembleForceResidualDerivativeWrtPressure( meshName, mesh, regionNames, dofManager, localMatrix, localRhs ); + assembleFluidMassResidualDerivativeWrtDisplacement( meshName, mesh, regionNames, dofManager, localMatrix, localRhs ); } ); + + //if hasStabilization via bubble - Apb + if constexpr (CONTACT_SOLVER::hasContactStabilization) + assembleMatrixPressureBubbleContribution( dt, const_cast< DomainPartition & >( domain ), dofManager, localMatrix, localRhs ); + } - void assembleForceResidualDerivativeWrtPressure( MeshLevel const & mesh, + virtual void assembleForceResidualDerivativeWrtPressure( string const & GEOS_UNUSED_PARAM(meshName), + MeshLevel const & mesh, string_array const & regionNames, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) - { - GEOS_MARK_FUNCTION; - - FaceManager const & faceManager = mesh.getFaceManager(); - NodeManager const & nodeManager = mesh.getNodeManager(); - EdgeManager const & edgeManager = mesh.getEdgeManager(); - ElementRegionManager const & elemManager = mesh.getElemManager(); - - ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); - ArrayOfArraysView< localIndex const > const faceToEdgeMap = faceManager.edgeList().toViewConst(); - arrayView2d< localIndex const > const & edgeToNodeMap = edgeManager.nodeList().toViewConst(); - arrayView2d< real64 const > faceCenters = faceManager.faceCenter(); - arrayView2d< real64 const > const & faceNormal = faceManager.faceNormal(); - arrayView1d< real64 const > faceAreas = faceManager.faceArea(); - - string const & dispDofKey = dofManager.getKey( fields::solidMechanics::totalDisplacement::key() ); - string const & flowDofKey = dofManager.getKey( this->getFlowDofKey() ); - - arrayView1d< globalIndex const > const & - dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); - globalIndex const rankOffset = dofManager.rankOffset(); - - // Get the coordinates for all nodes - arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition = nodeManager.referencePosition(); - - elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, - [&]( localIndex const, - FaceElementSubRegion const & subRegion ) - { - arrayView1d< globalIndex const > const & - flowDofNumber = subRegion.getReference< globalIndex_array >( flowDofKey ); - arrayView1d< real64 const > const & pressure = subRegion.getReference< array1d< real64 > >( fields::flow::pressure::key() ); - arrayView2d< localIndex const > const & elemsToFaces = subRegion.faceList().toViewConst(); - - forAll< serialPolicy >( subRegion.size(), [=, this]( localIndex const kfe ) - { - localIndex const kf0 = elemsToFaces[kfe][0]; - localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 ); - - real64 Nbar[3]; - Nbar[ 0 ] = faceNormal[elemsToFaces[kfe][0]][0] - faceNormal[elemsToFaces[kfe][1]][0]; - Nbar[ 1 ] = faceNormal[elemsToFaces[kfe][0]][1] - faceNormal[elemsToFaces[kfe][1]][1]; - Nbar[ 2 ] = faceNormal[elemsToFaces[kfe][0]][2] - faceNormal[elemsToFaces[kfe][1]][2]; - LvArray::tensorOps::normalize< 3 >( Nbar ); - globalIndex rowDOF[3 * m_maxFaceNodes]; // this needs to be changed when dealing with arbitrary element types - real64 nodeRHS[3 * m_maxFaceNodes]; - stackArray1d< real64, 3 * m_maxFaceNodes > dRdP( 3*m_maxFaceNodes ); - globalIndex colDOF[1]; - colDOF[0] = flowDofNumber[kfe]; // pressure is always first - - for( localIndex kf=0; kf<2; ++kf ) - { - localIndex const faceIndex = elemsToFaces[kfe][kf]; - - // Compute local area contribution for each node - stackArray1d< real64, FaceManager::maxFaceNodes() > nodalArea; - this->solidMechanicsSolver()->computeFaceNodalArea( elemsToFaces[kfe][kf], - nodePosition, - faceToNodeMap, - faceToEdgeMap, - edgeToNodeMap, - faceCenters, - faceNormal, - faceAreas, - nodalArea ); - for( localIndex a=0; a( globalNodalForce, Nbar, nodalForceMag ); - - for( localIndex i=0; i<3; ++i ) - { - rowDOF[3*a+i] = dispDofNumber[faceToNodeMap( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i ); - // Opposite sign w.r.t. theory because of minus sign in stiffness matrix definition (K < 0) - nodeRHS[3*a+i] = +globalNodalForce[i] * pow( -1, kf ); - - // Opposite sign w.r.t. theory because of minus sign in stiffness matrix definition (K < 0) - dRdP( 3*a+i ) = -nodalArea[a] * Nbar[i] * pow( -1, kf ); - } - } - - for( localIndex idof = 0; idof < numNodesPerFace * 3; ++idof ) - { - localIndex const localRow = LvArray::integerConversion< localIndex >( rowDOF[idof] - rankOffset ); - - if( localRow >= 0 && localRow < localMatrix.numRows() ) - { - localMatrix.addToRow< parallelHostAtomic >( localRow, - colDOF, - &dRdP[idof], - 1 ); - RAJA::atomicAdd( parallelHostAtomic{}, &localRhs[localRow], nodeRHS[idof] ); - } - } - } - } ); - } ); - } - - virtual void assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, + arrayView1d< real64 > const & localRhs ) = 0; + // { + // GEOS_MARK_FUNCTION; + + // FaceManager const & faceManager = mesh.getFaceManager(); + // NodeManager const & nodeManager = mesh.getNodeManager(); + // EdgeManager const & edgeManager = mesh.getEdgeManager(); + // ElementRegionManager const & elemManager = mesh.getElemManager(); + + // ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + // ArrayOfArraysView< localIndex const > const faceToEdgeMap = faceManager.edgeList().toViewConst(); + // arrayView2d< localIndex const > const & edgeToNodeMap = edgeManager.nodeList().toViewConst(); + // arrayView2d< real64 const > faceCenters = faceManager.faceCenter(); + // arrayView2d< real64 const > const & faceNormal = faceManager.faceNormal(); + // arrayView1d< real64 const > faceAreas = faceManager.faceArea(); + + // string const & dispDofKey = dofManager.getKey( fields::solidMechanics::totalDisplacement::key() ); + // string const & flowDofKey = dofManager.getKey( this->getFlowDofKey() ); + + // arrayView1d< globalIndex const > const & + // dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + // globalIndex const rankOffset = dofManager.rankOffset(); + + // // Get the coordinates for all nodes + // arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition = nodeManager.referencePosition(); + + // elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, + // [&]( localIndex const, + // FaceElementSubRegion const & subRegion ) + // { + // arrayView1d< globalIndex const > const & + // flowDofNumber = subRegion.getReference< globalIndex_array >( flowDofKey ); + // arrayView1d< real64 const > const & pressure = subRegion.getReference< array1d< real64 > >( fields::flow::pressure::key() ); + // arrayView2d< localIndex const > const & elemsToFaces = subRegion.faceList().toViewConst(); + + // forAll< serialPolicy >( subRegion.size(), [=, this]( localIndex const kfe ) + // { + // localIndex const kf0 = elemsToFaces[kfe][0]; + // localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 ); + + // real64 Nbar[3]; + // Nbar[ 0 ] = faceNormal[elemsToFaces[kfe][0]][0] - faceNormal[elemsToFaces[kfe][1]][0]; + // Nbar[ 1 ] = faceNormal[elemsToFaces[kfe][0]][1] - faceNormal[elemsToFaces[kfe][1]][1]; + // Nbar[ 2 ] = faceNormal[elemsToFaces[kfe][0]][2] - faceNormal[elemsToFaces[kfe][1]][2]; + // LvArray::tensorOps::normalize< 3 >( Nbar ); + // globalIndex rowDOF[3 * m_maxFaceNodes]; // this needs to be changed when dealing with arbitrary element types + // real64 nodeRHS[3 * m_maxFaceNodes]; + // stackArray1d< real64, 3 * m_maxFaceNodes > dRdP( 3*m_maxFaceNodes ); + // globalIndex colDOF[1]; + // colDOF[0] = flowDofNumber[kfe]; // pressure is always first + + // for( localIndex kf=0; kf<2; ++kf ) + // { + // localIndex const faceIndex = elemsToFaces[kfe][kf]; + + // // Compute local area contribution for each node + // stackArray1d< real64, FaceManager::maxFaceNodes() > nodalArea; + // this->solidMechanicsSolver()->computeFaceNodalArea( elemsToFaces[kfe][kf], + // nodePosition, + // faceToNodeMap, + // faceToEdgeMap, + // edgeToNodeMap, + // faceCenters, + // faceNormal, + // faceAreas, + // nodalArea ); + // for( localIndex a=0; a( globalNodalForce, Nbar, nodalForceMag ); + + // for( localIndex i=0; i<3; ++i ) + // { + // rowDOF[3*a+i] = dispDofNumber[faceToNodeMap( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i ); + // // Opposite sign w.r.t. theory because of minus sign in stiffness matrix definition (K < 0) + // nodeRHS[3*a+i] = +globalNodalForce[i] * pow( -1, kf ); + + // // Opposite sign w.r.t. theory because of minus sign in stiffness matrix definition (K < 0) + // dRdP( 3*a+i ) = -nodalArea[a] * Nbar[i] * pow( -1, kf ); + // } + // } + + // for( localIndex idof = 0; idof < numNodesPerFace * 3; ++idof ) + // { + // localIndex const localRow = LvArray::integerConversion< localIndex >( rowDOF[idof] - rankOffset ); + + // if( localRow >= 0 && localRow < localMatrix.numRows() ) + // { + // localMatrix.addToRow< parallelHostAtomic >( localRow, + // colDOF, + // &dRdP[idof], + // 1 ); + // RAJA::atomicAdd( parallelHostAtomic{}, &localRhs[localRow], nodeRHS[idof] ); + // } + // } + // } + // } ); + // } ); + // } + + virtual void assembleFluidMassResidualDerivativeWrtDisplacement( string const& meshName, + MeshLevel const & mesh, string_array const & regionNames, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) = 0; + /** + * @Brief assemble the contribution of matrix cell pressure on bubble DOFs + * with full Jacobian for fully-implicit coupling. + * @param dt the time step size + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param localMatrix the local system matrix + * @param localRhs the local system right-hand side vector + */ + virtual void assembleMatrixPressureBubbleContribution( real64 const GEOS_UNUSED_PARAM(dt), + DomainPartition & GEOS_UNUSED_PARAM(domain), + DofManager const & GEOS_UNUSED_PARAM(dofManager), + CRSMatrixView< real64, globalIndex const > const & GEOS_UNUSED_PARAM(localMatrix), + arrayView1d< real64 > const & GEOS_UNUSED_PARAM(localRhs) ) + { GEOS_WARNING("Should override"); }; + virtual void mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) override { GEOS_MARK_FUNCTION; @@ -698,14 +1160,13 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, return m_derivativeFluxResidual_dAperture->toViewConst(); } - virtual integer numFluidComponents() const = 0; - struct viewKeyStruct : public Base::viewKeyStruct {}; static const localIndex m_maxFaceNodes = 11; // Maximum number of nodes on a contact face std::unique_ptr< CRSMatrix< real64, localIndex > > m_derivativeFluxResidual_dAperture; + stdMap< string, localIndex > m_derivativeFluxResidual_dApertureOffsets; }; diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index 1b34ab94f65..1b14e979478 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -50,6 +50,12 @@ ENUM_STRINGS( StabilizationType, "Local" ); } +#define GEOS_MGR_STRATEGY_NOT_SUPPORTED() \ + virtual void setMGRStrategy() override \ + { \ + if( this->m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::mgr ) \ + GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}", this->getName(), this->getCatalogName() ) ); \ + } template< typename FLOW_SOLVER, typename MECHANICS_SOLVER = SolidMechanicsLagrangianFEM > class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER > diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp index 947a0a6ef4b..8a00ec6ba55 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp @@ -39,7 +39,11 @@ void SinglePhasePoromechanicsConformingFractures<>::setMGRStrategy() if( linearSolverParameters.preconditionerType != LinearSolverParameters::PreconditionerType::mgr ) return; - + + if( this->m_isThermal ) + { + GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for thermal {}", getName(), getCatalogName() ) ); + } linearSolverParameters.mgr.separateComponents = true; linearSolverParameters.dofsPerNode = 3; @@ -51,14 +55,21 @@ void SinglePhasePoromechanicsConformingFractures<>::setMGRStrategy() template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFractures< FLOW_SOLVER >:: -assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, +assembleFluidMassResidualDerivativeWrtDisplacement( string const & GEOS_UNUSED_PARAM(meshName), + MeshLevel const & mesh, string_array const & regionNames, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & GEOS_UNUSED_PARAM( localRhs ) ) { + +//TODO refactor with contact browser and kernel + GEOS_MARK_FUNCTION; + // TODO(thermal): getDerivativeFluxResidual_dNormalJump() below is sized/indexed one row per + // fracture element (mass-only) regardless of m_isThermal, should be 2 rows per fracture element for energy balance if m_isThermal + using namespace contact; FaceManager const & faceManager = mesh.getFaceManager(); @@ -225,23 +236,114 @@ assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, } ); } -template<> -void SinglePhasePoromechanicsConformingFractures< SinglePhaseReservoirAndWells<> >:: -assembleSystem( real64 const time_n, - real64 const dt, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) -{ - GEOS_MARK_FUNCTION; +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFractures::assembleForceResidualDerivativeWrtPressure( string const & GEOS_UNUSED_PARAM(meshName), + MeshLevel const & mesh, + string_array const & regionNames, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) + { + GEOS_MARK_FUNCTION; + + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + EdgeManager const & edgeManager = mesh.getEdgeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + ArrayOfArraysView< localIndex const > const faceToEdgeMap = faceManager.edgeList().toViewConst(); + arrayView2d< localIndex const > const & edgeToNodeMap = edgeManager.nodeList().toViewConst(); + arrayView2d< real64 const > faceCenters = faceManager.faceCenter(); + arrayView2d< real64 const > const & faceNormal = faceManager.faceNormal(); + arrayView1d< real64 const > faceAreas = faceManager.faceArea(); + + string const & dispDofKey = dofManager.getKey( fields::solidMechanics::totalDisplacement::key() ); + string const & flowDofKey = dofManager.getKey( this->getFlowDofKey() ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + globalIndex const rankOffset = dofManager.rankOffset(); + + // Get the coordinates for all nodes + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition = nodeManager.referencePosition(); + + elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, + [&]( localIndex const, + FaceElementSubRegion const & subRegion ) + { + arrayView1d< globalIndex const > const & + flowDofNumber = subRegion.getReference< globalIndex_array >( flowDofKey ); + arrayView1d< real64 const > const & pressure = subRegion.getReference< array1d< real64 > >( fields::flow::pressure::key() ); + arrayView2d< localIndex const > const & elemsToFaces = subRegion.faceList().toViewConst(); + + forAll< serialPolicy >( subRegion.size(), [=, this]( localIndex const kfe ) + { + localIndex const kf0 = elemsToFaces[kfe][0]; + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 ); + + real64 Nbar[3]; + Nbar[ 0 ] = faceNormal[elemsToFaces[kfe][0]][0] - faceNormal[elemsToFaces[kfe][1]][0]; + Nbar[ 1 ] = faceNormal[elemsToFaces[kfe][0]][1] - faceNormal[elemsToFaces[kfe][1]][1]; + Nbar[ 2 ] = faceNormal[elemsToFaces[kfe][0]][2] - faceNormal[elemsToFaces[kfe][1]][2]; + LvArray::tensorOps::normalize< 3 >( Nbar ); + globalIndex rowDOF[3 * m_maxFaceNodes]; // this needs to be changed when dealing with arbitrary element types + real64 nodeRHS[3 * m_maxFaceNodes]; + stackArray1d< real64, 3 * m_maxFaceNodes > dRdP( 3*m_maxFaceNodes ); + globalIndex colDOF[1]; + colDOF[0] = flowDofNumber[kfe]; // pressure is always first + + for( localIndex kf=0; kf<2; ++kf ) + { + localIndex const faceIndex = elemsToFaces[kfe][kf]; + + // Compute local area contribution for each node + stackArray1d< real64, FaceManager::maxFaceNodes() > nodalArea; + this->solidMechanicsSolver()->computeFaceNodalArea( elemsToFaces[kfe][kf], + nodePosition, + faceToNodeMap, + faceToEdgeMap, + edgeToNodeMap, + faceCenters, + faceNormal, + faceAreas, + nodalArea ); + for( localIndex a=0; a( globalNodalForce, Nbar, nodalForceMag ); + + for( localIndex i=0; i<3; ++i ) + { + rowDOF[3*a+i] = dispDofNumber[faceToNodeMap( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i ); + // Opposite sign w.r.t. theory because of minus sign in stiffness matrix definition (K < 0) + nodeRHS[3*a+i] = +globalNodalForce[i] * pow( -1, kf ); + + // Opposite sign w.r.t. theory because of minus sign in stiffness matrix definition (K < 0) + dRdP( 3*a+i ) = -nodalArea[a] * Nbar[i] * pow( -1, kf ); + } + } + + for( localIndex idof = 0; idof < numNodesPerFace * 3; ++idof ) + { + localIndex const localRow = LvArray::integerConversion< localIndex >( rowDOF[idof] - rankOffset ); + + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + localMatrix.addToRow< parallelHostAtomic >( localRow, + colDOF, + &dRdP[idof], + 1 ); + RAJA::atomicAdd( parallelHostAtomic{}, &localRhs[localRow], nodeRHS[idof] ); + } + } + } + } ); + } ); + } - Base::assembleSystem( time_n, dt, domain, dofManager, localMatrix, localRhs ); - // assemble well contributions - flowSolver()->wellSolver()->assembleSystem( time_n, dt, domain, dofManager, localMatrix, localRhs ); - flowSolver()->assembleCouplingTerms( time_n, dt, domain, dofManager, localMatrix, localRhs ); -} template class SinglePhasePoromechanicsConformingFractures<>; template class SinglePhasePoromechanicsConformingFractures< SinglePhaseReservoirAndWells<> >; diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp index f0ddd9bc8f1..8d2cca64bac 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp @@ -26,6 +26,7 @@ namespace geos { + template< typename FLOW_SOLVER = SinglePhaseBase > class SinglePhasePoromechanicsConformingFractures : public PoromechanicsConformingFractures< SinglePhasePoromechanics, FLOW_SOLVER > { @@ -81,32 +82,26 @@ class SinglePhasePoromechanicsConformingFractures : public PoromechanicsConformi * These functions provide the primary interface that is required for derived classes */ /**@{*/ - - virtual void setMGRStrategy() override - { - if( this->m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::mgr ) - GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}", this->getName(), this->getCatalogName())); - } - - virtual void assembleSystem( real64 const time_n, - real64 const dt, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) override - { Base::assembleSystem( time_n, dt, domain, dofManager, localMatrix, localRhs ); } + GEOS_MGR_STRATEGY_NOT_SUPPORTED() /**@}*/ protected: - virtual void assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, + virtual void assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, + MeshLevel const & mesh, string_array const & regionNames, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) override; + arrayView1d< real64 > const & localRhs ) override final; + + virtual void assembleForceResidualDerivativeWrtPressure( string const & GEOS_UNUSED_PARAM(meshName), + MeshLevel const & mesh, + string_array const & regionNames, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) override final; - virtual integer numFluidComponents() const override { return 1; } virtual string getFlowDofKey() const override { return SinglePhaseBase::viewKeyStruct::elemDofFieldString(); } diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp index 19962b85313..07f7f6dac1a 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp @@ -47,22 +47,6 @@ SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::SinglePhasePorome : Base( name, parent ) {} -template< typename FLOW_SOLVER > -void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupCoupling( DomainPartition const & domain, - DofManager & dofManager ) const -{ - GEOS_MARK_FUNCTION; - - // 1. Poromechanical coupling in the bulk (from base class) - Base::setupCoupling( domain, dofManager ); - - // 2. Pressure - bubble displacement coupling in the fracture - dofManager.addCoupling( m_pressureKey, - contact::totalBubbleDisplacement::key(), - DofManager::Connector::Elem ); -} - - template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setSparsityPattern( DomainPartition & domain, DofManager & dofManager, @@ -71,188 +55,11 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setSparsityP { GEOS_MARK_FUNCTION; - // Initialize ALM contact solver internal data structures - // These must be called before assembling the contact-dependent pattern. - this->solidMechanicsSolver()->createFaceTypeList( domain ); - this->solidMechanicsSolver()->updateStickSlipList( domain ); - this->solidMechanicsSolver()->createBubbleCellList( domain ); - - // 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 ) - { - 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() ) - { - 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]; - } - } - - // 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 > -void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleSystem( real64 const time_n, - real64 const dt, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) -{ - GEOS_MARK_FUNCTION; - - // 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 ); - - // Assemble flux terms and get dFluidResidual/dAperture - this->flowSolver()->assembleHydrofracFluxTerms( time_n, - dt, - domain, - dofManager, - localMatrix, - localRhs, - 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 > -void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleElementBasedContributions( real64 const time_n, - real64 const dt, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) -{ - GEOS_MARK_FUNCTION; - - // Assemble poromechanics terms (from base class) - Base::assembleElementBasedTerms( time_n, dt, domain, dofManager, localMatrix, localRhs ); - - // Flow accumulation for fractures - this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel & mesh, - string_array const & regionNames ) - { - mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const, - FaceElementSubRegion const & subRegion ) - { - this->flowSolver()->accumulationAssemblyLaunch( dofManager, subRegion, localMatrix, localRhs ); - } ); - } ); - - // Assemble contact terms (ALM) - note: assembleContact requires time and dt - this->solidMechanicsSolver()->assembleContact( time_n, dt, domain, dofManager, localMatrix, localRhs ); -} - -template< typename FLOW_SOLVER > -void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleCouplingTerms( real64 const time_n, - real64 const dt, - DomainPartition const & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) -{ - GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( time_n, dt ); - - // These steps must occur after the fluxes are assembled because that's when DerivativeFluxResidual_dAperture is filled. - this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshName, - MeshLevel const & mesh, - string_array const & regionNames ) - { - // Assemble Force Residual w.r.t. pressure (Aup) - fracture pressure contribution - assembleForceResidualDerivativeWrtPressure( meshName, mesh, regionNames, dofManager, localMatrix, localRhs ); - - // Assemble Fluid mass residual w.r.t. displacement (Apu) - assembleFluidMassResidualDerivativeWrtDisplacement( meshName, mesh, regionNames, dofManager, localMatrix, localRhs ); - } ); + // Recompute fracture face/element geometry and rebuild the ALM contact solver's internal lists. + // These must happen before assembling the contact-dependent pattern; setSparsityPattern() + this->solidMechanicsSolver()->updateFractureGeometry( domain ); - // Assemble matrix cell pressure contribution on bubble DOFs (Abp_matrix) - // This must be outside the lambda because it uses regionBasedKernelApplication - assembleMatrixPressureBubbleContribution( dt, const_cast< DomainPartition & >( domain ), dofManager, localMatrix, localRhs ); + Base::setSparsityPattern(domain,dofManager,localMatrix, pattern); } template< typename FLOW_SOLVER > @@ -330,454 +137,6 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::updateState( this->flowSolver()->updateStencilWeights( domain ); } -template< typename FLOW_SOLVER > -void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: -setUpDflux_dApertureMatrix( DomainPartition & domain ) -{ - NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); - FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); - FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); - string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); - - // 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 ) - { - 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 ) - { - if( region.subRegionType() == SurfaceElementRegion::SurfaceSubRegionType::faceElement ) - { - ++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 ); - - fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) - { - for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) - { - localIndex const numFluxElems = stencil.stencilSize( iconn ); - typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); - for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) - { - 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; - } - } - } ); - } ); - - 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 ) - { - for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) - { - localIndex const numFluxElems = stencil.stencilSize( iconn ); - typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); - - for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) - { - 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( row, - rowOffset + sei[iconn][k1], - 0.0 ); - } - } - } - } ); - } ); -} - -template< typename FLOW_SOLVER > -void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: -addTransmissibilityCouplingNNZ( DomainPartition const & domain, - DofManager const & dofManager, - arrayView1d< localIndex > const & rowLengths ) const -{ - GEOS_MARK_FUNCTION; - - this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel const & mesh, - string_array const & ) - { - ElementRegionManager const & elemManager = mesh.getElemManager(); - - string const flowDofKey = dofManager.getKey( m_pressureKey ); - - globalIndex const rankOffset = dofManager.rankOffset(); - - NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); - FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); - FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); - - fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) - { - for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) - { - localIndex const numFluxElems = stencil.stencilSize( iconn ); - typename SurfaceElementStencil::IndexContainerViewConstType const & seri = stencil.getElementRegionIndices(); - typename SurfaceElementStencil::IndexContainerViewConstType const & sesri = stencil.getElementSubRegionIndices(); - typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); - - FaceElementSubRegion const & elementSubRegion = - elemManager.getRegion( seri[iconn][0] ).getSubRegion< FaceElementSubRegion >( sesri[iconn][0] ); - - ArrayOfArraysView< localIndex const > const elemsToNodes = elementSubRegion.nodeList().toViewConst(); - - arrayView1d< globalIndex const > const faceElementDofNumber = - elementSubRegion.getReference< array1d< globalIndex > >( flowDofKey ); - - for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) - { - globalIndex const activeFlowDOF = faceElementDofNumber[sei[iconn][k0]]; - globalIndex const rowNumber = activeFlowDOF - rankOffset; - - if( rowNumber >= 0 && rowNumber < rowLengths.size() ) - { - for( localIndex k1 = 0; k1 < numFluxElems; ++k1 ) - { - // The coupling with the nodal displacements of the cell itself has already been added by the dofManager - // so we only add the coupling with the nodal displacements of the neighbors. - if( k1 != k0 ) - { - localIndex const numNodesPerElement = elemsToNodes[sei[iconn][k1]].size(); - // Nodal displacement DOFs (3 per node, 2 faces) - rowLengths[rowNumber] += 3 * numNodesPerElement; - // Bubble DOFs (3 per face, 2 faces = 6 total) - rowLengths[rowNumber] += 6; - } - } - } - } - } - } ); - } ); -} - -template< typename FLOW_SOLVER > -void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: -addTransmissibilityCouplingPattern( DomainPartition const & domain, - DofManager const & dofManager, - SparsityPatternView< globalIndex > const & pattern ) const -{ - GEOS_MARK_FUNCTION; - - using namespace contact; - - this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel const & mesh, - string_array const & ) - { - FaceManager const & faceManager = mesh.getFaceManager(); - NodeManager const & nodeManager = mesh.getNodeManager(); - ElementRegionManager const & elemManager = mesh.getElemManager(); - - string const dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); - string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); - string const flowDofKey = dofManager.getKey( m_pressureKey ); - - arrayView1d< globalIndex const > const & - dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); - arrayView1d< globalIndex const > const & - bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); - ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); - - // Get the finite volume method used to compute the fluxes - NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); - FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); - FluxApproximationBase const & fvDiscretization = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); - - SurfaceElementRegion const & fractureRegion = - elemManager.getRegion< SurfaceElementRegion >( this->solidMechanicsSolver()->getUniqueFractureRegionName() ); - FaceElementSubRegion const & fractureSubRegion = - fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); - - GEOS_ERROR_IF( !fractureSubRegion.hasWrapper( flow::pressure::key() ), - this->getDataContext() << ": The fracture subregion must contain pressure field." ); - - arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst(); - - arrayView1d< globalIndex const > const & - flowDofNumber = fractureSubRegion.getReference< globalIndex_array >( flowDofKey ); - - globalIndex const rankOffset = dofManager.rankOffset(); - - fvDiscretization.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) - { - forAll< serialPolicy >( stencil.size(), [=] ( localIndex const iconn ) - { - localIndex const numFluxElems = stencil.stencilSize( iconn ); - - // A fracture connector has to be an edge shared by two faces - if( numFluxElems == 2 ) - { - typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); - - // First index: face element. Second index: node - for( localIndex kf = 0; kf < 2; ++kf ) - { - // Set row DOF index - // Note that the 1-kf index is intentional, as this is coupling the pressure of one face cell - // to the nodes of the adjacent cell - localIndex const rowIndex = flowDofNumber[sei[iconn][1 - kf]] - rankOffset; - - if( rowIndex >= 0 && rowIndex < pattern.numRows() ) - { - // Get fracture, face and region/subregion/element indices (for elements on both sides) - localIndex const fractureIndex = sei[iconn][kf]; - - // Get the number of nodes - localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elem2dToFaces[fractureIndex][0] ); - - // Loop over the two sides of each fracture element - for( localIndex kf1 = 0; kf1 < 2; ++kf1 ) - { - localIndex const faceIndex = elem2dToFaces[fractureIndex][kf1]; - - // Save the list of DOF associated with nodes (displacement) - for( localIndex a = 0; a < numNodesPerFace; ++a ) - { - for( localIndex i = 0; i < 3; ++i ) - { - globalIndex const colIndex = dispDofNumber[faceToNodeMap( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i ); - pattern.insertNonZero( rowIndex, colIndex ); - } - } - - // Save the list of DOF associated with bubble displacement - for( localIndex i = 0; i < 3; ++i ) - { - globalIndex const colIndex = bubbleDofNumber[faceIndex] + LvArray::integerConversion< globalIndex >( i ); - pattern.insertNonZero( rowIndex, colIndex ); - } - } - } - } - } - } ); - } ); - } ); -} - -template< typename FLOW_SOLVER > -void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: -addPressureForceCouplingNNZ( DomainPartition const & domain, - DofManager const & dofManager, - arrayView1d< localIndex > const & rowLengths ) const -{ - GEOS_MARK_FUNCTION; - - using namespace contact; - - this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel const & mesh, - string_array const & ) - { - FaceManager const & faceManager = mesh.getFaceManager(); - NodeManager const & nodeManager = mesh.getNodeManager(); - ElementRegionManager const & elemManager = mesh.getElemManager(); - - string const dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); - string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); - - arrayView1d< globalIndex const > const & - dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); - arrayView1d< globalIndex const > const & - bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); - ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); - - globalIndex const rankOffset = dofManager.rankOffset(); - - string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); - SurfaceElementRegion const & fractureRegion = - elemManager.getRegion< SurfaceElementRegion >( fractureRegionName ); - FaceElementSubRegion const & fractureSubRegion = - fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); - - arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst(); - - // For each fracture element, add NNZ for (displacement_row, pressure_col) and (bubble_row, pressure_col) - forAll< serialPolicy >( fractureSubRegion.size(), [=, &rowLengths] ( localIndex const kfe ) - { - localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elem2dToFaces[kfe][0] ); - - // For displacement DOFs: add 1 pressure column per displacement DOF row - for( localIndex kf = 0; kf < 2; ++kf ) - { - localIndex const faceIndex = elem2dToFaces[kfe][kf]; - for( localIndex a = 0; a < numNodesPerFace; ++a ) - { - for( localIndex i = 0; i < 3; ++i ) - { - globalIndex const rowNumber = dispDofNumber[faceToNodeMap( faceIndex, a )] + i - rankOffset; - if( rowNumber >= 0 && rowNumber < rowLengths.size() ) - { - rowLengths[rowNumber] += 1; // One pressure column - } - } - } - } - - // For bubble DOFs: add 1 pressure column per bubble DOF row - for( localIndex kf = 0; kf < 2; ++kf ) - { - localIndex const faceIndex = elem2dToFaces[kfe][kf]; - for( localIndex i = 0; i < 3; ++i ) - { - globalIndex const rowNumber = bubbleDofNumber[faceIndex] + i - rankOffset; - if( rowNumber >= 0 && rowNumber < rowLengths.size() ) - { - rowLengths[rowNumber] += 1; // One pressure column - } - } - } - } ); - } ); -} - -template< typename FLOW_SOLVER > -void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: -addPressureForceCouplingPattern( DomainPartition const & domain, - DofManager const & dofManager, - SparsityPatternView< globalIndex > const & pattern ) const -{ - GEOS_MARK_FUNCTION; - - using namespace contact; - - this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel const & mesh, - string_array const & ) - { - FaceManager const & faceManager = mesh.getFaceManager(); - NodeManager const & nodeManager = mesh.getNodeManager(); - ElementRegionManager const & elemManager = mesh.getElemManager(); - - string const dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); - string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); - string const flowDofKey = dofManager.getKey( m_pressureKey ); - - arrayView1d< globalIndex const > const & - dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); - arrayView1d< globalIndex const > const & - bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); - ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); - - globalIndex const rankOffset = dofManager.rankOffset(); - - string const & fractureRegionName = this->solidMechanicsSolver()->getUniqueFractureRegionName(); - SurfaceElementRegion const & fractureRegion = - elemManager.getRegion< SurfaceElementRegion >( fractureRegionName ); - FaceElementSubRegion const & fractureSubRegion = - fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); - - arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst(); - arrayView1d< globalIndex const > const & - flowDofNumber = fractureSubRegion.getReference< globalIndex_array >( flowDofKey ); - - // For each fracture element, add pattern for (displacement_row, pressure_col) and (bubble_row, pressure_col) - forAll< serialPolicy >( fractureSubRegion.size(), [=] ( localIndex const kfe ) - { - globalIndex const pressureColIndex = flowDofNumber[kfe]; - localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elem2dToFaces[kfe][0] ); - - // For displacement DOFs - for( localIndex kf = 0; kf < 2; ++kf ) - { - localIndex const faceIndex = elem2dToFaces[kfe][kf]; - for( localIndex a = 0; a < numNodesPerFace; ++a ) - { - for( localIndex i = 0; i < 3; ++i ) - { - globalIndex const rowIndex = dispDofNumber[faceToNodeMap( faceIndex, a )] + i - rankOffset; - if( rowIndex >= 0 && rowIndex < pattern.numRows() ) - { - pattern.insertNonZero( rowIndex, pressureColIndex ); - } - } - } - } - - // For bubble DOFs - for( localIndex kf = 0; kf < 2; ++kf ) - { - localIndex const faceIndex = elem2dToFaces[kfe][kf]; - for( localIndex i = 0; i < 3; ++i ) - { - globalIndex const rowIndex = bubbleDofNumber[faceIndex] + i - rankOffset; - if( rowIndex >= 0 && rowIndex < pattern.numRows() ) - { - pattern.insertNonZero( rowIndex, pressureColIndex ); - } - } - } - } ); - } ); -} template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: @@ -797,7 +156,7 @@ assembleForceResidualDerivativeWrtPressure( string const & meshName, string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); string const & bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); - string const & flowDofKey = dofManager.getKey( m_pressureKey ); + string const & flowDofKey = dofManager.getKey( this->getFlowDofKey() ); arrayView1d< globalIndex const > const dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); @@ -855,6 +214,8 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, GEOS_MARK_FUNCTION; GEOS_UNUSED_VAR( regionNames ); + // TODO(thermal): getDerivativeFluxResidual_dNormalJump() below is sized/indexed one row per + // fracture element (mass-only) regardless of m_isThermal -- extend to 2 rows per fracture element for energy balance if m_isThermal using namespace contact; FaceManager const & faceManager = mesh.getFaceManager(); @@ -866,7 +227,7 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, // 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(); + dFluxResidual_dNormalJump = this->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 ) ); @@ -874,7 +235,7 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); string const & bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); - string const & presDofKey = dofManager.getKey( m_pressureKey ); + string const & presDofKey = dofManager.getKey( this->getFlowDofKey() ); arrayView1d< globalIndex const > const & dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); @@ -895,9 +256,6 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, localIndex const numElems = subRegion.size(); - // 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 ); dAperturedU.zero(); @@ -908,7 +266,6 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, 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, [&, dAperturedUView, dAperturedBView] ( string const &, finiteElement::FiniteElementBase const & subRegionFE, @@ -958,7 +315,6 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, // Note: dAperturedU/dB are computed as (1/area) * unitNormal^T * Atu/Atb // For accumulation: dR_accum/du = density * unitNormal^T * Atu (no 1/area factor) // For flux: dR_flux/du = dR/dAperture * (1/area) * unitNormal^T * Atu (with 1/area factor) - // So accumulation needs to multiply by area to cancel the 1/area in dAperturedU arrayView1d< real64 const > const area = subRegion.getElementArea().toViewConst(); forAll< serialPolicy >( numElems, [&]( localIndex const kfe ) @@ -1136,124 +492,7 @@ assembleFluidMassResidualDerivativeWrtDisplacement( string const & meshName, } ); } -template< typename FLOW_SOLVER > -void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: -addMatrixPressureBubbleCouplingNNZ( DomainPartition const & domain, - DofManager const & dofManager, - arrayView1d< localIndex > const & rowLengths ) const -{ - GEOS_MARK_FUNCTION; - - using namespace contact; - - this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel const & mesh, - string_array const & regionNames ) - { - FaceManager const & faceManager = mesh.getFaceManager(); - ElementRegionManager const & elemManager = mesh.getElemManager(); - - string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); - string const flowDofKey = dofManager.getKey( m_pressureKey ); - arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); - - globalIndex const rankOffset = dofManager.rankOffset(); - - // Loop over matrix cell regions that have bubbles - elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, - [&]( localIndex const, CellElementSubRegion const & subRegion ) - { - arrayView1d< localIndex const > const bubbleElems = subRegion.bubbleElementsList(); - arrayView2d< localIndex const > const elemsToFaces = subRegion.faceElementsList(); - arrayView1d< globalIndex const > const pressureDofNumber = subRegion.getReference< array1d< globalIndex > >( flowDofKey ); - - forAll< serialPolicy >( bubbleElems.size(), [=, &rowLengths]( localIndex const kk ) - { - localIndex const k = bubbleElems[kk]; - localIndex const faceIndex = elemsToFaces[kk][0]; - - // (bubble_row, pressure_col): 1 pressure column for each of the 3 bubble DOFs - for( localIndex i = 0; i < 3; ++i ) - { - globalIndex const rowNumber = bubbleDofNumber[faceIndex] + i - rankOffset; - if( rowNumber >= 0 && rowNumber < rowLengths.size() ) - { - rowLengths[rowNumber] += 1; // One pressure DOF from matrix cell - } - } - - // (pressure_row, bubble_col): the matrix cell pressure couples to its 3 bubble DOFs (A_pb) - globalIndex const pRow = pressureDofNumber[k] - rankOffset; - if( pRow >= 0 && pRow < rowLengths.size() ) - { - rowLengths[pRow] += 3; // Three bubble DOFs - } - } ); - } ); - } ); -} -template< typename FLOW_SOLVER > -void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: -addMatrixPressureBubbleCouplingPattern( DomainPartition const & domain, - DofManager const & dofManager, - SparsityPatternView< globalIndex > const & pattern ) const -{ - GEOS_MARK_FUNCTION; - - using namespace contact; - - this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel const & mesh, - string_array const & regionNames ) - { - FaceManager const & faceManager = mesh.getFaceManager(); - ElementRegionManager const & elemManager = mesh.getElemManager(); - - string const bubbleDofKey = dofManager.getKey( totalBubbleDisplacement::key() ); - string const flowDofKey = dofManager.getKey( m_pressureKey ); - - arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); - - globalIndex const rankOffset = dofManager.rankOffset(); - - // Loop over matrix cell regions that have bubbles - elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, - [&]( localIndex const, CellElementSubRegion const & subRegion ) - { - arrayView1d< localIndex const > const bubbleElems = subRegion.bubbleElementsList(); - arrayView2d< localIndex const > const elemsToFaces = subRegion.faceElementsList(); - arrayView1d< globalIndex const > const pressureDofNumber = subRegion.getReference< array1d< globalIndex > >( flowDofKey ); - - forAll< serialPolicy >( bubbleElems.size(), [=]( localIndex const kk ) - { - localIndex const k = bubbleElems[kk]; - localIndex const faceIndex = elemsToFaces[kk][0]; - globalIndex const pressureColIndex = pressureDofNumber[k]; - - // (bubble_row, pressure_col) : A_bp - for( localIndex i = 0; i < 3; ++i ) - { - globalIndex const rowIndex = bubbleDofNumber[faceIndex] + i - rankOffset; - if( rowIndex >= 0 && rowIndex < pattern.numRows() ) - { - pattern.insertNonZero( rowIndex, pressureColIndex ); - } - } - - // (pressure_row, bubble_col) : A_pb -- transpose location - globalIndex const pRow = pressureDofNumber[k] - rankOffset; - if( pRow >= 0 && pRow < pattern.numRows() ) - { - for( localIndex i = 0; i < 3; ++i ) - { - pattern.insertNonZero( pRow, bubbleDofNumber[faceIndex] + i ); - } - } - } ); - } ); - } ); -} template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: @@ -1267,7 +506,7 @@ assembleMatrixPressureBubbleContribution( real64 const dt, using namespace contact; - string const flowDofKey = dofManager.getKey( m_pressureKey ); + string const flowDofKey = dofManager.getKey( this->getFlowDofKey() ); string const mechanicsDiscretizationName = this->solidMechanicsSolver()->getDiscretizationName(); this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp index 010ceb64ec4..ccacbfa5b97 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -20,23 +20,30 @@ #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_ #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_ -#include "physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp" +#include "common/logger/Logger.hpp" +// #include "physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFracturesALM.hpp" #include "physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp" +#include "physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp" +#include "physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp" namespace geos { template< typename FLOW_SOLVER = SinglePhaseBase > -class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromechanics< FLOW_SOLVER, SolidMechanicsAugmentedLagrangianContact > +class SinglePhasePoromechanicsConformingFracturesALM : public PoromechanicsConformingFractures< SinglePhasePoromechanics, FLOW_SOLVER, SolidMechanicsAugmentedLagrangianContact > { public: - using Base = SinglePhasePoromechanics< FLOW_SOLVER, SolidMechanicsAugmentedLagrangianContact >; + using Base = PoromechanicsConformingFractures< SinglePhasePoromechanics, FLOW_SOLVER , SolidMechanicsAugmentedLagrangianContact >; using Base::m_solvers; using Base::m_dofManager; using Base::m_localMatrix; using Base::m_rhs; using Base::m_solution; + using Base::m_maxFaceNodes; + + using Base::m_derivativeFluxResidual_dAperture; + using Base::m_derivativeFluxResidual_dApertureOffsets; /// True when the flow solver carries well degrees of freedom. static constexpr bool hasWells = std::is_same_v< FLOW_SOLVER, SinglePhaseReservoirAndWells<> >; @@ -90,20 +97,12 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec */ /**@{*/ - virtual void setupCoupling( DomainPartition const & domain, - DofManager & dofManager ) const 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, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) override final; virtual void updateState( DomainPartition & domain ) override final; @@ -118,6 +117,12 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec // 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 (this->m_isThermal) + { + if( this->m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::mgr ) + GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}", this->getName(), this->getCatalogName() ) ); + } + if constexpr ( hasWells ) { linearSolverParameters.mgr.strategy = @@ -128,6 +133,7 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec linearSolverParameters.mgr.strategy = LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsConformingFracturesALM; } + linearSolverParameters.mgr.separateComponents = true; GEOS_LOG_LEVEL_RANK_0( logInfo::LinearSolver, @@ -142,19 +148,6 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec 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: @@ -162,30 +155,8 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec struct viewKeyStruct : public Base::viewKeyStruct {}; - static const localIndex m_maxFaceNodes=11; // Maximum number of nodes on a contact face - /** - * @Brief assemble the element-based contributions - * @param time_n the current time - * @param dt the time step - * @param domain the physical domain object - * @param dofManager degree-of-freedom manager associated with the linear system - * @param localMatrix the local system matrix - * @param localRhs the local system right-hand side vector - */ - void assembleElementBasedContributions( real64 const time_n, - real64 const dt, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ); - - virtual void assembleCouplingTerms( real64 const time_n, - real64 const dt, - DomainPartition const & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) override final; + void assembleForceResidualDerivativeWrtPressure( string const & meshName, MeshLevel const & mesh, @@ -199,111 +170,19 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec string_array const & regionNames, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ); - - /** - * @Brief add the nnz induced by the flux-aperture coupling - * @param domain the physical domain object - * @param dofManager degree-of-freedom manager associated with the linear system - * @param rowLenghts the nnz in each row - */ - void addTransmissibilityCouplingNNZ( DomainPartition const & domain, - DofManager const & dofManager, - arrayView1d< localIndex > const & rowLengths ) const; - - /** - * @Brief add the sparsity pattern induced by the flux-aperture coupling - * @param domain the physical domain object - * @param dofManager degree-of-freedom manager associated with the linear system - * @param pattern the sparsity pattern - */ - void addTransmissibilityCouplingPattern( DomainPartition const & domain, - DofManager const & dofManager, - SparsityPatternView< globalIndex > const & pattern ) const; + arrayView1d< real64 > const & localRhs ) override final; - /** - * @Brief add the nnz induced by the pressure-force coupling (Aup, Abp) - * @param domain the physical domain object - * @param dofManager degree-of-freedom manager associated with the linear system - * @param rowLenghts the nnz in each row - */ - void addPressureForceCouplingNNZ( DomainPartition const & domain, - DofManager const & dofManager, - arrayView1d< localIndex > const & rowLengths ) const; - - /** - * @Brief add the sparsity pattern induced by the pressure-force coupling (Aup, Abp) - * @param domain the physical domain object - * @param dofManager degree-of-freedom manager associated with the linear system - * @param pattern the sparsity pattern - */ - void addPressureForceCouplingPattern( DomainPartition const & domain, - DofManager const & dofManager, - SparsityPatternView< globalIndex > const & pattern ) const; - /** - * @Brief add the nnz induced by the matrix pressure-bubble coupling (Abp_matrix) - * This handles the contribution of matrix cell pressure on bubble DOFs. - * @param domain the physical domain object - * @param dofManager degree-of-freedom manager associated with the linear system - * @param rowLengths the nnz in each row - */ - void addMatrixPressureBubbleCouplingNNZ( DomainPartition const & domain, - DofManager const & dofManager, - arrayView1d< localIndex > const & rowLengths ) const; - /** - * @Brief add the sparsity pattern induced by the matrix pressure-bubble coupling - * @param domain the physical domain object - * @param dofManager degree-of-freedom manager associated with the linear system - * @param pattern the sparsity pattern - */ - void addMatrixPressureBubbleCouplingPattern( DomainPartition const & domain, - DofManager const & dofManager, - SparsityPatternView< globalIndex > const & pattern ) const; - /** - * @Brief assemble the contribution of matrix cell pressure on bubble DOFs - * with full Jacobian for fully-implicit coupling. - * @param dt the time step size - * @param domain the physical domain object - * @param dofManager degree-of-freedom manager associated with the linear system - * @param localMatrix the local system matrix - * @param localRhs the local system right-hand side vector - */ void assembleMatrixPressureBubbleContribution( real64 const dt, DomainPartition & domain, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ); - - /** - * @brief Set up the Dflux_dApertureMatrix object - * - * @param domain - */ - void setUpDflux_dApertureMatrix( DomainPartition & domain ); - - std::unique_ptr< CRSMatrix< real64, localIndex > > & getRefDerivativeFluxResidual_dAperture() - { - return m_derivativeFluxResidual_dAperture; - } - - CRSMatrixView< real64, localIndex const > getDerivativeFluxResidual_dNormalJump() - { - return m_derivativeFluxResidual_dAperture->toViewConstSizes(); - } - - CRSMatrixView< real64 const, localIndex const > getDerivativeFluxResidual_dNormalJump() const - { - return m_derivativeFluxResidual_dAperture->toViewConst(); - } - - std::unique_ptr< CRSMatrix< real64, localIndex > > m_derivativeFluxResidual_dAperture; + arrayView1d< real64 > const & localRhs ) override; - stdMap< string, localIndex > m_derivativeFluxResidual_dApertureOffsets; - string const m_pressureKey = SinglePhaseBase::viewKeyStruct::elemDofFieldString(); + void updateHydraulicApertureAndFracturePermeability( DomainPartition & domain ); }; diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp index 9ac8d570c7f..076eb46bc23 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp @@ -145,6 +145,7 @@ void SinglePhasePoromechanicsEmbeddedFractures::setSparsityPattern( DomainPartit pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data()); // Copy the original nonzeros + //TODO appendSparsityPattern(pattern,patternDiag); for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) { globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous(); diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp index 50d804eaf7f..b2483c31613 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp @@ -111,6 +111,8 @@ class SinglePhaseReservoirAndWells : public CoupledReservoirAndWellsBase< RESERV void updateMass( CellElementSubRegion & subRegion ) const { flowSolver()->updateMass( subRegion ); } + integer numFluidComponents() const { return flowSolver()->numFluidComponents(); }; + protected: virtual void initializePreSubGroups() override; diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALM.hpp new file mode 100644 index 00000000000..d4ea015b384 --- /dev/null +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -0,0 +1,330 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP +#define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP + +#include "physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernel.hpp" +#include "physicsSolvers/fluidFlow/kernels/singlePhase/FluxKernelsHelper.hpp" +#include "codingUtilities/Utilities.hpp" + +namespace geos +{ + +namespace singlePhasePoromechanicsConformingFracturesALMKernels +{ + +template< integer NUM_EQN, integer NUM_DOF > +class ConnectorBasedAssemblyKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_EQN, NUM_DOF, SurfaceElementStencilWrapper > +{ +public: + + /** + * @brief The type for element-based data. Consists entirely of ArrayView's. + * + * Can be converted from ElementRegionManager::ElementViewConstAccessor + * by calling .toView() or .toViewConst() on an accessor instance + */ + template< typename VIEWTYPE > + using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >; + + using AbstractBase = singlePhaseFVMKernels::FluxComputeKernelBase; + using DofNumberAccessor = AbstractBase::DofNumberAccessor; + using SinglePhaseFlowAccessors = AbstractBase::SinglePhaseFlowAccessors; + using SinglePhaseFluidAccessors = AbstractBase::SinglePhaseFluidAccessors; + using PermeabilityAccessors = AbstractBase::PermeabilityAccessors; + using FracturePermeabilityAccessors = StencilMaterialAccessors< constitutive::PermeabilityBase, + fields::permeability::dPerm_dDispJump >; + + using AbstractBase::m_dt; + using AbstractBase::m_rankOffset; + using AbstractBase::m_dofNumber; + using AbstractBase::m_permeability; + using AbstractBase::m_dPerm_dPres; + using AbstractBase::m_gravCoef; + using AbstractBase::m_pres; + using AbstractBase::m_mob; + using AbstractBase::m_dMob; + using AbstractBase::m_dMob_dPres; + using AbstractBase::m_dens; + using AbstractBase::m_dDens; + + + using Base = singlePhaseFVMKernels::FluxComputeKernel< NUM_EQN, NUM_DOF, SurfaceElementStencilWrapper >; + using Base::numDof; + using Base::numEqn; + using Base::maxNumElems; + using Base::maxNumConns; + using Base::maxStencilSize; + using Base::m_stencilWrapper; + using Base::m_seri; + using Base::m_sesri; + using Base::m_sei; + + ConnectorBasedAssemblyKernel( globalIndex const rankOffset, + SurfaceElementStencilWrapper const & stencilWrapper, + DofNumberAccessor const & flowDofNumberAccessor, + SinglePhaseFlowAccessors const & singlePhaseFlowAccessors, + SinglePhaseFluidAccessors const & singlePhaseFluidAccessors, + PermeabilityAccessors const & permeabilityAccessors, + FracturePermeabilityAccessors const & fracturePermeabilityAccessors, + real64 const & dt, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) + : Base( rankOffset, + stencilWrapper, + flowDofNumberAccessor, + singlePhaseFlowAccessors, + singlePhaseFluidAccessors, + permeabilityAccessors, + dt, + localMatrix, + localRhs ), + m_dR_dAper( dR_dAper ), + m_dPerm_dDispJump( fracturePermeabilityAccessors.get( fields::permeability::dPerm_dDispJump {} ) ) + {} + + + /** + * @struct StackVariables + * @brief Kernel variables (dof numbers, jacobian and residual) located on the stack + */ + struct StackVariables : public Base::StackVariables + { +public: + + /** + * @brief Constructor for the stack variables + * @param[in] size size of the stencil for this connection + * @param[in] numElems number of elements for this connection + */ + GEOS_HOST_DEVICE + StackVariables( localIndex const size, localIndex numElems ) + : Base::StackVariables( size, numElems ), + localColIndices( numElems ), + dFlux_dAperture( numElems, size ) + {} + + stackArray1d< localIndex, maxNumElems > localColIndices; + + stackArray2d< real64, maxNumElems * maxStencilSize > dFlux_dAperture; + + /// Derivatives of transmissibility with respect to the dispJump + real64 dTrans_dDispJump[maxNumConns][2][3]{}; + }; + + /** + * @brief Performs the setup phase for the kernel. + * @param[in] iconn the connection index + * @param[in] stack the stack variables + */ + GEOS_HOST_DEVICE + void setup( localIndex const iconn, + StackVariables & stack ) const + { + // set degrees of freedom indices for this face + for( integer i = 0; i < stack.stencilSize; ++i ) + { + globalIndex const offset = m_dofNumber[m_seri( iconn, i )][m_sesri( iconn, i )][m_sei( iconn, i )]; + for( integer jdof = 0; jdof < numDof; ++jdof ) + { + stack.dofColIndices[i * numDof + jdof] = offset + jdof; + } + stack.localColIndices[ i ] = m_sei( iconn, i ); + } + } + + /** + * @brief Compute the local flux contributions to the residual and Jacobian + * @tparam FUNC the type of the function that can be used to customize the computation of the flux + * @param[in] iconn the connection index + * @param[inout] stack the stack variables + * @param[in] NoOpFunc the function used to customize the computation of the flux + */ + template< typename FUNC = NoOpFunc > + GEOS_HOST_DEVICE + void computeFlux( localIndex const iconn, + StackVariables & stack, + FUNC && kernelOp = NoOpFunc{} ) const + + { + + m_stencilWrapper.computeWeights( iconn, + m_permeability, + m_dPerm_dPres, + m_dPerm_dDispJump, + stack.transmissibility, + stack.dTrans_dPres, + stack.dTrans_dDispJump ); + + + localIndex k[2]; + localIndex connectionIndex = 0; + for( k[0]=0; k[0] + GEOS_HOST_DEVICE + void complete( localIndex const iconn, + StackVariables & stack, + FUNC && kernelOp = NoOpFunc{} ) const + { + // Call Base::complete to assemble the mass balance equations + // In the lambda, fill the dR_dAper matrix + Base::complete( iconn, stack, [&] ( integer const i, + localIndex const localRow ) + { + + localIndex const row = LvArray::integerConversion< localIndex >( m_sei( iconn, i ) ); + + m_dR_dAper.addToRowBinarySearch< parallelDeviceAtomic >( row, + stack.localColIndices.data(), + stack.dFlux_dAperture[i].dataIfContiguous(), + stack.stencilSize ); + // call the lambda to assemble additional terms, such as thermal terms + kernelOp( i, localRow ); + } ); + } + +private: + + CRSMatrixView< real64, localIndex const > m_dR_dAper; + + ElementViewConst< arrayView4d< real64 const > > const m_dPerm_dDispJump; +}; + + + +/** + * @class ConnectorBasedAssemblyKernelFactory + */ +class ConnectorBasedAssemblyKernelFactory +{ +public: + + /** + * @brief Create a new kernel and launch + * @tparam POLICY the policy used in the RAJA kernel + * @tparam STENCILWRAPPER the type of the stencil wrapper + * @param[in] rankOffset the offset of my MPI rank + * @param[in] dofKey string to get the element degrees of freedom numbers + * @param[in] solverName name of the solver (to name accessors) + * @param[in] elemManager reference to the element region manager + * @param[in] stencilWrapper reference to the stencil wrapper + * @param[in] dt time step size + * @param[inout] localMatrix the local CRS matrix + * @param[inout] localRhs the local right-hand side vector + */ + template< typename POLICY > + static void + createAndLaunch( globalIndex const rankOffset, + string const & dofKey, + string const & solverName, + ElementRegionManager const & elemManager, + SurfaceElementStencilWrapper const & stencilWrapper, + real64 const & dt, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) + { + integer constexpr NUM_DOF = 1; // pressure + integer constexpr NUM_EQN = 1; + + ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > > flowDofNumberAccessor = + elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); + flowDofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); + + using kernelType = ConnectorBasedAssemblyKernel< NUM_EQN, NUM_DOF >; + typename kernelType::SinglePhaseFlowAccessors flowAccessors( elemManager, solverName ); + typename kernelType::SinglePhaseFluidAccessors fluidAccessors( elemManager, solverName ); + typename kernelType::PermeabilityAccessors permAccessors( elemManager, solverName ); + typename kernelType::FracturePermeabilityAccessors fracPermAccessors( elemManager, solverName ); + + kernelType kernel( rankOffset, stencilWrapper, flowDofNumberAccessor, + flowAccessors, fluidAccessors, permAccessors, fracPermAccessors, + dt, localMatrix, localRhs, dR_dAper ); + + kernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); + } +}; + +} // namespace SinglePhasePoromechanicsConformingFracturesALMKernels + +} // namespace geos + +#endif //GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFracturesALM.hpp new file mode 100644 index 00000000000..66f90840ec1 --- /dev/null +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFracturesALM.hpp @@ -0,0 +1,391 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_THERMALSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP +#define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_THERMALSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP + +#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALM.hpp" +#include "physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernelBase.hpp" + +namespace geos +{ + +namespace thermalSinglePhasePoromechanicsConformingFracturesALMKernels +{ + +template< integer NUM_EQN, integer NUM_DOF > +class ConnectorBasedAssemblyKernel : public singlePhasePoromechanicsConformingFracturesALMKernels::ConnectorBasedAssemblyKernel< NUM_EQN, NUM_DOF > +{ +public: + + /** + * @brief The type for element-based data. Consists entirely of ArrayView's. + * + * Can be converted from ElementRegionManager::ElementViewConstAccessor + * by calling .toView() or .toViewConst() on an accessor instance + */ + template< typename VIEWTYPE > + using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >; + + using SinglePhaseFVMAbstractBase = singlePhaseFVMKernels::FluxComputeKernelBase; + using DofNumberAccessor = SinglePhaseFVMAbstractBase::DofNumberAccessor; + using SinglePhaseFlowAccessors = SinglePhaseFVMAbstractBase::SinglePhaseFlowAccessors; + using SinglePhaseFluidAccessors = SinglePhaseFVMAbstractBase::SinglePhaseFluidAccessors; + using PermeabilityAccessors = SinglePhaseFVMAbstractBase::PermeabilityAccessors; + using FracturePermeabilityAccessors = StencilMaterialAccessors< constitutive::PermeabilityBase, + fields::permeability::dPerm_dDispJump >; + using SinglePhaseFVMAbstractBase::m_dt; + using SinglePhaseFVMAbstractBase::m_rankOffset; + using SinglePhaseFVMAbstractBase::m_dofNumber; + using SinglePhaseFVMAbstractBase::m_gravCoef; + using SinglePhaseFVMAbstractBase::m_mob; + using SinglePhaseFVMAbstractBase::m_dens; + using SinglePhaseFVMAbstractBase::m_dDens; + using SinglePhaseFVMAbstractBase::m_dMob; + + using SinglePhaseFVMBase = singlePhaseFVMKernels::FluxComputeKernel< NUM_EQN, NUM_DOF, SurfaceElementStencilWrapper >; + using SinglePhaseFVMBase::numDof; + using SinglePhaseFVMBase::numEqn; + using SinglePhaseFVMBase::maxNumElems; + using SinglePhaseFVMBase::maxNumConns; + using SinglePhaseFVMBase::maxStencilSize; + using SinglePhaseFVMBase::m_stencilWrapper; + using SinglePhaseFVMBase::m_seri; + using SinglePhaseFVMBase::m_sesri; + using SinglePhaseFVMBase::m_sei; + using SinglePhaseFVMBase::m_ghostRank; + + using Base = singlePhasePoromechanicsConformingFracturesALMKernels::ConnectorBasedAssemblyKernel< NUM_EQN, NUM_DOF >; + + using ThermalSinglePhaseFlowAccessors = + StencilAccessors< fields::flow::temperature >; + + using ThermalSinglePhaseFluidAccessors = + StencilMaterialAccessors< constitutive::SingleFluidBase, + fields::singlefluid::enthalpy, + fields::singlefluid::dEnthalpy >; + + using ThermalConductivityAccessors = + StencilMaterialAccessors< constitutive::SinglePhaseThermalConductivityBase, + fields::thermalconductivity::effectiveConductivity >; + + + + ConnectorBasedAssemblyKernel( globalIndex const rankOffset, + SurfaceElementStencilWrapper const & stencilWrapper, + DofNumberAccessor const & flowDofNumberAccessor, + SinglePhaseFlowAccessors const & singlePhaseFlowAccessors, + ThermalSinglePhaseFlowAccessors const & thermalSinglePhaseFlowAccessors, + SinglePhaseFluidAccessors const & singlePhaseFluidAccessors, + ThermalSinglePhaseFluidAccessors const & thermalSinglePhaseFluidAccessors, + PermeabilityAccessors const & permeabilityAccessors, + FracturePermeabilityAccessors const & edfmPermeabilityAccessors, + ThermalConductivityAccessors const & thermalConductivityAccessors, + real64 const & dt, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) + : Base( rankOffset, + stencilWrapper, + flowDofNumberAccessor, + singlePhaseFlowAccessors, + singlePhaseFluidAccessors, + permeabilityAccessors, + edfmPermeabilityAccessors, + dt, + localMatrix, + localRhs, + dR_dAper ), + m_temp( thermalSinglePhaseFlowAccessors.get( fields::flow::temperature {} ) ), + m_enthalpy( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::enthalpy {} ) ), + m_dEnthalpy( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy {} ) ), + m_thermalConductivity( thermalConductivityAccessors.get( fields::thermalconductivity::effectiveConductivity {} ) ) + {} + + + /** + * @struct StackVariables + * @brief Kernel variables (dof numbers, jacobian and residual) located on the stack + */ + struct StackVariables : public Base::StackVariables + { +public: + + /** + * @brief Constructor for the stack variables + * @param[in] size size of the stencil for this connection + * @param[in] numElems number of elements for this connection + */ + GEOS_HOST_DEVICE + StackVariables( localIndex const size, localIndex numElems ) + : Base::StackVariables( size, numElems ), + energyFlux( 0.0 ), + dEnergyFlux_dTrans( 0.0 ), + dEnergyFlux_dP( size ), + dEnergyFlux_dT( size ), + dEnergyFlux_dDispJump( size, 3 ) + {} + using SinglePhaseFVMBase::StackVariables::stencilSize; + using SinglePhaseFVMBase::StackVariables::numFluxElems; + using SinglePhaseFVMBase::StackVariables::transmissibility; + using SinglePhaseFVMBase::StackVariables::dTrans_dPres; + using SinglePhaseFVMBase::StackVariables::dofColIndices; + using SinglePhaseFVMBase::StackVariables::localFlux; + using SinglePhaseFVMBase::StackVariables::localFluxJacobian; + + // Thermal transmissibility (for now, no derivatives) + + real64 thermalTransmissibility[maxNumConns][2]{}; + + // Energy fluxes and derivatives + + /// Energy fluxes + real64 energyFlux; + /// Derivative of the Energy fluxes wrt transmissibility + real64 dEnergyFlux_dTrans; + /// Derivatives of energy fluxes wrt pressure + stackArray1d< real64, maxStencilSize > dEnergyFlux_dP; + /// Derivatives of energy fluxes wrt temperature + stackArray1d< real64, maxStencilSize > dEnergyFlux_dT; + /// Derivatives of energy fluxes wrt dispJump + stackArray2d< real64, maxStencilSize *3 > dEnergyFlux_dDispJump{}; + + }; + + /** + * @brief Compute the local flux contributions to the residual and Jacobian + * @tparam FUNC the type of the function that can be used to customize the computation of the flux + * @param[in] iconn the connection index + * @param[inout] stack the stack variables + */ + GEOS_HOST_DEVICE + void computeFlux( localIndex const iconn, + StackVariables & stack ) const + + { + // *********************************************** + // First, we call the base computeFlux to compute: + // 1) compFlux and its derivatives (including derivatives wrt temperature), + // 2) enthalpy part of energyFlux and its derivatives (including derivatives wrt temperature) + // + // Computing dFlux_dT and the enthalpy flux requires quantities already computed in the base computeFlux, + // such as potGrad, fluxVal, and the indices of the upwind cell + // We use the lambda below (called **inside** the phase loop of the base computeFlux) to access these variables + Base::computeFlux( iconn, stack, [&] ( localIndex const (&k)[2], + localIndex const (&seri)[2], + localIndex const (&sesri)[2], + localIndex const (&sei)[2], + localIndex const, + real64 const & alpha, + real64 const & mobility, + real64 const & potGrad, + real64 const & massFlux, + real64 const & dMassFlux_dTrans, + real64 const (&dMassFlux_dP)[2] ) + { + real64 trans[2] = {stack.transmissibility[0][0], stack.transmissibility[0][1]}; + real64 dMassFlux_dT[2]{}; + + singlePhaseFluxKernelsHelper::computeEnthalpyFlux( seri, sesri, sei, + trans, + m_enthalpy, + m_dEnthalpy, + m_gravCoef, + m_dDens, + m_dMob, + alpha, + mobility, + potGrad, + massFlux, + dMassFlux_dTrans, + dMassFlux_dP, + dMassFlux_dT, + stack.energyFlux, + stack.dEnergyFlux_dTrans, + stack.dEnergyFlux_dP, + stack.dEnergyFlux_dT ); + + // add dMassFlux_dT to localFluxJacobian + for( integer ke = 0; ke < 2; ++ke ) + { + localIndex const localDofIndexTemp = k[ke] * numDof + numDof - 1; + stack.localFluxJacobian[k[0]*numEqn][localDofIndexTemp] += m_dt * dMassFlux_dT[ke]; + stack.localFluxJacobian[k[1]*numEqn][localDofIndexTemp] -= m_dt * dMassFlux_dT[ke]; + } + } ); + + // ***************************************************** + // Computation of the conduction term in the energy flux + // Note that the enthalpy term in the energy was computed above + // Note that this term is computed using an explicit treatment of conductivity for now + + // Step 1: compute the thermal transmissibilities at this face + m_stencilWrapper.computeWeights( iconn, + m_thermalConductivity, + m_thermalConductivity, // we have to pass something here, so we just use thermal conductivity + stack.thermalTransmissibility, + stack.dTrans_dPres ); // again, we have to pass something here, but this is unused for now + + localIndex k[2]; + localIndex connectionIndex = 0; + + for( k[0] = 0; k[0] < stack.numFluxElems; ++k[0] ) + { + for( k[1] = k[0] + 1; k[1] < stack.numFluxElems; ++k[1] ) + { + real64 const thermalTrans[2] = { stack.thermalTransmissibility[connectionIndex][0], stack.thermalTransmissibility[connectionIndex][1] }; + + localIndex const seri[2] = {m_seri( iconn, k[0] ), m_seri( iconn, k[1] )}; + localIndex const sesri[2] = {m_sesri( iconn, k[0] ), m_sesri( iconn, k[1] )}; + localIndex const sei[2] = {m_sei( iconn, k[0] ), m_sei( iconn, k[1] )}; + + // Step 2: compute temperature difference at the interface + singlePhaseFluxKernelsHelper::computeConductiveFlux( seri, sesri, sei, m_temp, thermalTrans, stack.energyFlux, stack.dEnergyFlux_dT ); + + // add energyFlux and its derivatives to localFlux and localFluxJacobian + stack.localFlux[k[0]*numEqn + numEqn - 1] += m_dt * stack.energyFlux; + stack.localFlux[k[1]*numEqn + numEqn - 1] -= m_dt * stack.energyFlux; + + for( integer ke = 0; ke < 2; ++ke ) + { + integer const localDofIndexPres = k[ke] * numDof; + stack.localFluxJacobian[k[0]*numEqn + numEqn - 1][localDofIndexPres] = m_dt * stack.dEnergyFlux_dP[ke]; + stack.localFluxJacobian[k[1]*numEqn + numEqn - 1][localDofIndexPres] = -m_dt * stack.dEnergyFlux_dP[ke]; + integer const localDofIndexTemp = localDofIndexPres + 1; + stack.localFluxJacobian[k[0]*numEqn + numEqn - 1][localDofIndexTemp] = m_dt * stack.dEnergyFlux_dT[ke]; + stack.localFluxJacobian[k[1]*numEqn + numEqn - 1][localDofIndexTemp] = -m_dt * stack.dEnergyFlux_dT[ke]; + } + + connectionIndex++; + } + } + } + + /** + * @brief Performs the complete phase for the kernel. + * @param[in] iconn the connection index + * @param[inout] stack the stack variables + */ + GEOS_HOST_DEVICE + void complete( localIndex const iconn, + StackVariables & stack ) const + { + // Call Base::complete to assemble the mass balance equations + // In the lambda, add contribution to residual and jacobian into the energy balance equation + Base::complete( iconn, stack, [&] ( integer const i, + localIndex const localRow ) + { + // The no. of fluxes is equal to the no. of equations in m_localRhs and m_localMatrix + RAJA::atomicAdd( parallelDeviceAtomic{}, &SinglePhaseFVMAbstractBase::m_localRhs[localRow + numEqn-1], stack.localFlux[i * numEqn + numEqn-1] ); + + SinglePhaseFVMAbstractBase::m_localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( localRow + numEqn-1, + stack.dofColIndices.data(), + stack.localFluxJacobian[i * numEqn + numEqn-1].dataIfContiguous(), + stack.stencilSize * numDof ); + + } ); + } + + +private: + + /// Views on temperature + ElementViewConst< arrayView1d< real64 const > > const m_temp; + + /// Views on enthalpies + ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > > const m_enthalpy; + + /// Views on enthalpies + ElementViewConst< arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > > const m_dEnthalpy; + + /// View on thermal conductivity + ElementViewConst< arrayView3d< real64 const > > m_thermalConductivity; + +}; + + + +/** + * @class ConnectorBasedAssemblyKernelFactory + */ +class ConnectorBasedAssemblyKernelFactory +{ +public: + + /** + * @brief Create a new kernel and launch + * @tparam POLICY the policy used in the RAJA kernel + * @tparam STENCILWRAPPER the type of the stencil wrapper + * @param[in] rankOffset the offset of my MPI rank + * @param[in] dofKey string to get the element degrees of freedom numbers + * @param[in] solverName name of the solver (to name accessors) + * @param[in] elemManager reference to the element region manager + * @param[in] stencilWrapper reference to the stencil wrapper + * @param[in] dt time step size + * @param[inout] localMatrix the local CRS matrix + * @param[inout] localRhs the local right-hand side vector + */ + template< typename POLICY > + static void + createAndLaunch( globalIndex const rankOffset, + string const & flowDofKey, + string const & solverName, + ElementRegionManager const & elemManager, + SurfaceElementStencilWrapper const & stencilWrapper, + real64 const & dt, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) + { + integer constexpr NUM_DOF = 2; // pressure + temperature + integer constexpr NUM_EQN = 2; // mass balance + energy balance + + + ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > > flowDofNumberAccessor = + elemManager.constructArrayViewAccessor< globalIndex, 1 >( flowDofKey ); + flowDofNumberAccessor.setName( solverName + "/accessors/" + flowDofKey ); + + using kernelType = ConnectorBasedAssemblyKernel< NUM_EQN, NUM_DOF >; + typename kernelType::SinglePhaseFlowAccessors flowAccessors( elemManager, solverName ); + typename kernelType::ThermalSinglePhaseFlowAccessors thermalFlowAccessors( elemManager, solverName ); + + typename kernelType::SinglePhaseFluidAccessors fluidAccessors( elemManager, solverName ); + typename kernelType::ThermalSinglePhaseFluidAccessors thermalFluidAccessors( elemManager, solverName ); + + typename kernelType::PermeabilityAccessors permAccessors( elemManager, solverName ); + typename kernelType::FracturePermeabilityAccessors edfmPermAccessors( elemManager, solverName ); + typename kernelType::ThermalConductivityAccessors thermalConductivityAccessors( elemManager, solverName ); + + kernelType kernel( rankOffset, stencilWrapper, + flowDofNumberAccessor, + flowAccessors, thermalFlowAccessors, fluidAccessors, thermalFluidAccessors, + permAccessors, edfmPermAccessors, thermalConductivityAccessors, + dt, localMatrix, localRhs, dR_dAper ); + + kernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); + } +}; + + + +} // namespace SinglePhaseProppantFluxKernels + +} // namespace geos + +#endif // GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_THERMALSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp index 3faa6bdff59..c3384cb2bbd 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp @@ -27,6 +27,8 @@ #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp" #include "common/GEOS_RAJA_Interface.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" +#include "finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp" +#include "finiteElement/elementFormulations/H1_QuadrilateralFace_Lagrange1_GaussLegendre2.hpp" namespace geos { @@ -34,6 +36,8 @@ namespace geos using namespace dataRepository; using namespace constitutive; using namespace fields; +using namespace finiteElement; + ContactSolverBase::ContactSolverBase( const string & name, Group * const parent ): @@ -256,4 +260,328 @@ void ContactSolverBase::setConstitutiveNamesCallSuper( ElementSubRegionBase & su } } +void ContactSolverBase::computeFaceNodalArea( localIndex const kf0, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, + ArrayOfArraysView< localIndex const > const & faceToNodeMap, + ArrayOfArraysView< localIndex const > const & faceToEdgeMap, + arrayView2d< localIndex const > const & edgeToNodeMap, + arrayView2d< real64 const > const faceCenters, + arrayView2d< real64 const > const faceNormals, + arrayView1d< real64 const > const faceAreas, + stackArray1d< real64, FaceManager::maxFaceNodes() > & basisIntegrals ) const +{ + GEOS_MARK_FUNCTION; + localIndex const TriangularPermutation[3] = { 0, 1, 2 }; + localIndex const QuadrilateralPermutation[4] = { 0, 1, 3, 2 }; + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 ); + + basisIntegrals.resize( numNodesPerFace ); + for( localIndex a = 0; a < numNodesPerFace; ++a ) + { + basisIntegrals[a] = 0.0; + } + localIndex const * const permutation = ( numNodesPerFace == 3 ) ? TriangularPermutation : QuadrilateralPermutation; + if( numNodesPerFace == 3 ) + { + real64 xLocal[3][3]; + for( localIndex a = 0; a < numNodesPerFace; ++a ) + { + for( localIndex j = 0; j < 3; ++j ) + { + xLocal[a][j] = nodePosition[faceToNodeMap( kf0, permutation[a] )][j]; + } + } + real64 N[3]; + for( localIndex q=0; q 4 && numNodesPerFace <= m_maxFaceNodes ) + { + // we need to L2 projector based on VEM to approximate the quadrature weights + // we need to use extra geometry information to computing L2 projector + + localIndex const MFN = m_maxFaceNodes; // Max number of face vertices. + localIndex const faceIndex = kf0; + localIndex const numFaceNodes = faceToNodeMap[ faceIndex ].size(); + + // get the face center and normal. + real64 const faceArea = faceAreas[ faceIndex ]; + localIndex faceToNodes[ MFN ]; + localIndex faceToEdges[ MFN ]; + for( localIndex i = 0; i < numFaceNodes; ++i ) + { + faceToNodes[i] = faceToNodeMap[ faceIndex ][ i ]; + faceToEdges[i] = faceToEdgeMap[ faceIndex ][ i ]; + } + // - get outward face normal and center + real64 faceNormal[3] = { faceNormals[faceIndex][0], + faceNormals[faceIndex][1], + faceNormals[faceIndex][2] }; + real64 const faceCenter[3] { faceCenters[faceIndex][0], + faceCenters[faceIndex][1], + faceCenters[faceIndex][2] }; + // - compute integrals calling auxiliary method + real64 threeDMonomialIntegrals[3] = { 0.0 }; + real64 const invCellDiameter = 0.0; + real64 const cellCenter[3] { 0.0, 0.0, 0.0 }; + computeFaceIntegrals( nodePosition, + faceToNodes, + faceToEdges, + numFaceNodes, + faceArea, + faceCenter, + faceNormal, + edgeToNodeMap, + invCellDiameter, + cellCenter, + basisIntegrals, + threeDMonomialIntegrals ); + } + else + { + GEOS_ERROR( GEOS_FMT( "Face with {} nodes. Only triangles and quadrilaterals and PEBI prisms up to 11 sides are supported.", + numNodesPerFace ), + getDataContext() ); + } +} + +void ContactSolverBase::computeFaceIntegrals( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodesCoords, + localIndex const (&faceToNodes)[11], + localIndex const (&faceToEdges)[11], + localIndex const & numFaceVertices, + real64 const & faceArea, + real64 const (&faceCenter)[3], + real64 const (&faceNormal)[3], + arrayView2d< localIndex const > const & edgeToNodes, + real64 const & invCellDiameter, + real64 const (&cellCenter)[3], + stackArray1d< real64, FaceManager::maxFaceNodes() > & basisIntegrals, + real64 (& threeDMonomialIntegrals)[3] ) const +{ + GEOS_MARK_FUNCTION; + localIndex const MFN = m_maxFaceNodes; // Max number of face vertices. + basisIntegrals.resize( numFaceVertices ); + // Rotate the face. + // - compute rotation matrix. + real64 faceRotationMatrix[ 3 ][ 3 ]; + computationalGeometry::RotationMatrix_3D( faceNormal, faceRotationMatrix ); + // - below we compute the diameter, the rotated vertices and the rotated center. + real64 faceRotatedVertices[ MFN ][ 2 ]; + real64 faceDiameter = 0; + + for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) + { + // apply the transpose (that is the inverse) of the rotation matrix to face vertices. + // NOTE: + // the second and third rows of the transpose of the rotation matrix rotate on the 2D face. + faceRotatedVertices[numVertex][0] = + faceRotationMatrix[ 0 ][ 1 ]*nodesCoords( faceToNodes[ numVertex ], 0 ) + + faceRotationMatrix[ 1 ][ 1 ]*nodesCoords( faceToNodes[ numVertex ], 1 ) + + faceRotationMatrix[ 2 ][ 1 ]*nodesCoords( faceToNodes[ numVertex ], 2 ); + faceRotatedVertices[numVertex][1] = + faceRotationMatrix[ 0 ][ 2 ]*nodesCoords( faceToNodes[ numVertex ], 0 ) + + faceRotationMatrix[ 1 ][ 2 ]*nodesCoords( faceToNodes[ numVertex ], 1 ) + + faceRotationMatrix[ 2 ][ 2 ]*nodesCoords( faceToNodes[ numVertex ], 2 ); + } + + faceDiameter = computationalGeometry::computeDiameter< 2 >( faceRotatedVertices, + numFaceVertices ); + real64 const invFaceDiameter = 1.0/faceDiameter; + // - rotate the face centroid as done for the vertices. + real64 faceRotatedCentroid[2]; + faceRotatedCentroid[0] = + faceRotationMatrix[ 0 ][ 1 ]*faceCenter[0] + + faceRotationMatrix[ 1 ][ 1 ]*faceCenter[1] + + faceRotationMatrix[ 2 ][ 1 ]*faceCenter[2]; + faceRotatedCentroid[1] = + faceRotationMatrix[ 0 ][ 2 ]*faceCenter[0] + + faceRotationMatrix[ 1 ][ 2 ]*faceCenter[1] + + faceRotationMatrix[ 2 ][ 2 ]*faceCenter[2]; + // - compute edges' lengths, outward pointing normals and local edge-to-nodes map. + real64 edgeOutwardNormals[ MFN ][ 2 ]; + real64 edgeLengths[ MFN ]; + localIndex localEdgeToNodes[ MFN ][ 2 ]; + + for( localIndex numEdge = 0; numEdge < numFaceVertices; ++numEdge ) + { + if( edgeToNodes( faceToEdges[numEdge], 0 ) == faceToNodes[ numEdge ] ) + { + localEdgeToNodes[ numEdge ][ 0 ] = numEdge; + localEdgeToNodes[ numEdge ][ 1 ] = (numEdge+1)%numFaceVertices; + } + else + { + localEdgeToNodes[ numEdge ][ 0 ] = (numEdge+1)%numFaceVertices; + localEdgeToNodes[ numEdge ][ 1 ] = numEdge; + } + real64 edgeTangent[2]; + edgeTangent[0] = faceRotatedVertices[(numEdge+1)%numFaceVertices][0] - + faceRotatedVertices[numEdge][0]; + edgeTangent[1] = faceRotatedVertices[(numEdge+1)%numFaceVertices][1] - + faceRotatedVertices[numEdge][1]; + edgeOutwardNormals[numEdge][0] = edgeTangent[1]; + edgeOutwardNormals[numEdge][1] = -edgeTangent[0]; + real64 signTestVector[2]; + signTestVector[0] = faceRotatedVertices[numEdge][0] - faceRotatedCentroid[0]; + signTestVector[1] = faceRotatedVertices[numEdge][1] - faceRotatedCentroid[1]; + if( signTestVector[0]*edgeOutwardNormals[numEdge][0] + + signTestVector[1]*edgeOutwardNormals[numEdge][1] < 0 ) + { + edgeOutwardNormals[numEdge][0] = -edgeOutwardNormals[numEdge][0]; + edgeOutwardNormals[numEdge][1] = -edgeOutwardNormals[numEdge][1]; + } + edgeLengths[numEdge] = LvArray::math::sqrt< real64 >( edgeTangent[0]*edgeTangent[0] + + edgeTangent[1]*edgeTangent[1] ); + edgeOutwardNormals[numEdge][0] /= edgeLengths[numEdge]; + edgeOutwardNormals[numEdge][1] /= edgeLengths[numEdge]; + } + + // Compute boundary quadrature weights (also equal to the integrals of basis functions on the + // boundary). + real64 boundaryQuadratureWeights[ MFN ]; + for( localIndex numWeight = 0; numWeight < numFaceVertices; ++numWeight ) + boundaryQuadratureWeights[numWeight] = 0.0; + for( localIndex numEdge = 0; numEdge < numFaceVertices; ++numEdge ) + { + boundaryQuadratureWeights[ localEdgeToNodes[ numEdge ][ 0 ] ] += 0.5*edgeLengths[numEdge]; + boundaryQuadratureWeights[ localEdgeToNodes[ numEdge ][ 1 ] ] += 0.5*edgeLengths[numEdge]; + } + + // Compute scaled monomials' integrals on edges. + real64 monomBoundaryIntegrals[3] = { 0.0 }; + for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) + { + monomBoundaryIntegrals[0] += boundaryQuadratureWeights[ numVertex ]; + monomBoundaryIntegrals[1] += (faceRotatedVertices[ numVertex ][ 0 ] - faceRotatedCentroid[0]) * + invFaceDiameter*boundaryQuadratureWeights[ numVertex ]; + monomBoundaryIntegrals[2] += (faceRotatedVertices[ numVertex ][ 1 ] - faceRotatedCentroid[1]) * + invFaceDiameter*boundaryQuadratureWeights[ numVertex ]; + } + + // Compute non constant 2D and 3D scaled monomials' integrals on the face. + real64 monomInternalIntegrals[2] = { 0.0 }; + for( localIndex numSubTriangle = 0; numSubTriangle < numFaceVertices; ++numSubTriangle ) + { + localIndex const nextVertex = (numSubTriangle+1)%numFaceVertices; + // - compute value of 2D monomials at the quadrature point on the sub-triangle (the + // barycenter). + // The result is ((v(0)+v(1)+faceCenter)/3 - faceCenter) / faceDiameter = + // = (v(0) + v(1) - 2*faceCenter)/(3*faceDiameter). + real64 monomialValues[2]; + for( localIndex i = 0; i < 2; ++i ) + { + monomialValues[i] = (faceRotatedVertices[numSubTriangle][i] + + faceRotatedVertices[nextVertex][i] - + 2.0*faceRotatedCentroid[i]) / (3.0*faceDiameter); + } + // compute value of 3D monomials at the quadrature point on the sub-triangle (the + // barycenter). The result is + // ((v(0) + v(1) + faceCenter)/3 - cellCenter)/cellDiameter. + real64 threeDMonomialValues[3]; + for( localIndex i = 0; i < 3; ++i ) + { + threeDMonomialValues[i] = ( (faceCenter[i] + + nodesCoords[faceToNodes[ numSubTriangle ]][i] + + nodesCoords[faceToNodes[ nextVertex ]][i]) / 3.0 - + cellCenter[i] ) * invCellDiameter; + } + // compute quadrature weight associated to the quadrature point (the area of the + // sub-triangle). + real64 edgesTangents[2][2]; // used to compute the area of the sub-triangle + for( localIndex i = 0; i < 2; ++i ) + { + edgesTangents[0][i] = faceRotatedVertices[numSubTriangle][i] - faceRotatedCentroid[i]; + } + for( localIndex i = 0; i < 2; ++i ) + { + edgesTangents[1][i] = faceRotatedVertices[nextVertex][i] - faceRotatedCentroid[i]; + } + real64 subTriangleArea = 0.5*LvArray::math::abs + ( edgesTangents[0][0]*edgesTangents[1][1] - + edgesTangents[0][1]*edgesTangents[1][0] ); + // compute the integrals on the sub-triangle and add it to the global integrals + for( localIndex i = 0; i < 2; ++i ) + { + monomInternalIntegrals[ i ] += monomialValues[ i ]*subTriangleArea; + } + for( localIndex i = 0; i < 3; ++i ) + { + // threeDMonomialIntegrals is assumed to be initialized to 0 by the caller + threeDMonomialIntegrals[ i ] += threeDMonomialValues[ i ]*subTriangleArea; + } + } + + // Compute integral of basis functions times normal derivative of monomials on the boundary. + real64 basisTimesMonomNormalDerBoundaryInt[ MFN ][ 2 ]; + for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) + { + for( localIndex i = 0; i < 2; ++i ) + { + basisTimesMonomNormalDerBoundaryInt[ numVertex ][ i ] = 0.0; + } + } + for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) + { + for( localIndex i = 0; i < 2; ++i ) + { + real64 thisEdgeIntTimesNormal_i = edgeOutwardNormals[numVertex][i]*edgeLengths[numVertex]; + basisTimesMonomNormalDerBoundaryInt[ localEdgeToNodes[ numVertex ][ 0 ] ][i] += thisEdgeIntTimesNormal_i; + basisTimesMonomNormalDerBoundaryInt[ localEdgeToNodes[ numVertex ][ 1 ] ][i] += thisEdgeIntTimesNormal_i; + } + } + for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) + { + for( localIndex i = 0; i < 2; ++i ) + { + basisTimesMonomNormalDerBoundaryInt[ numVertex ][ i ] *= 0.5*invFaceDiameter; + } + } + + // Compute integral mean of basis functions on this face. + real64 const invFaceArea = 1.0/faceArea; + real64 const monomialDerivativeInverse = (faceDiameter*faceDiameter)*invFaceArea; + for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) + { + real64 piNablaDofs[ 3 ]; + piNablaDofs[ 1 ] = monomialDerivativeInverse * + basisTimesMonomNormalDerBoundaryInt[ numVertex ][ 0 ]; + piNablaDofs[ 2 ] = monomialDerivativeInverse * + basisTimesMonomNormalDerBoundaryInt[ numVertex ][ 1 ]; + piNablaDofs[ 0 ] = (boundaryQuadratureWeights[ numVertex ] - + piNablaDofs[ 1 ]*monomBoundaryIntegrals[ 1 ] - + piNablaDofs[ 2 ]*monomBoundaryIntegrals[ 2 ])/monomBoundaryIntegrals[ 0 ]; + basisIntegrals[ numVertex ] = piNablaDofs[ 0 ]*faceArea + + (piNablaDofs[ 1 ]*monomInternalIntegrals[ 0 ] + + piNablaDofs[ 2 ]*monomInternalIntegrals[ 1 ]); + } +} + } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp index ccfc68f6ea7..e9b642a69c9 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp @@ -30,6 +30,8 @@ namespace geos class ContactSolverBase : public SolidMechanicsLagrangianFEM { public: + static constexpr bool hasContactStabilization = false; + ContactSolverBase( const string & name, Group * const parent ); @@ -49,6 +51,29 @@ class ContactSolverBase : public SolidMechanicsLagrangianFEM void synchronizeFractureState( DomainPartition & domain ) const; + void computeFaceNodalArea( localIndex const kf0, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, + ArrayOfArraysView< localIndex const > const & faceToNodeMap, + ArrayOfArraysView< localIndex const > const & faceToEdgeMap, + arrayView2d< localIndex const > const & edgeToNodeMap, + arrayView2d< real64 const > const faceCenters, + arrayView2d< real64 const > const faceNormals, + arrayView1d< real64 const > const faceAreas, + stackArray1d< real64, FaceManager::maxFaceNodes() > & basisIntegrals ) const; + + void computeFaceIntegrals( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodesCoords, + localIndex const (&faceToNodes)[11], + localIndex const (&faceToEdges)[11], + localIndex const & numFaceVertices, + real64 const & faceArea, + real64 const (&faceCenter)[3], + real64 const (&faceNormal)[3], + arrayView2d< localIndex const > const & edgeToNodes, + real64 const & invCellDiameter, + real64 const (&cellCenter)[3], + stackArray1d< real64, FaceManager::maxFaceNodes() > & basisIntegrals, + real64 ( &threeDMonomialIntegrals )[3] ) const; + struct viewKeyStruct : SolidMechanicsLagrangianFEM::viewKeyStruct { constexpr static char const * fractureStateString() { return "fractureState"; } @@ -59,6 +84,9 @@ class ContactSolverBase : public SolidMechanicsLagrangianFEM }; + static constexpr localIndex m_maxFaceNodes = 11; // Maximum number of nodes on a contact face + + protected: virtual void postInputInitialization() override; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index f5ec003c05d..ea4f857724a 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -19,6 +19,7 @@ #include "SolidMechanicsAugmentedLagrangianContact.hpp" +#include "common/logger/Logger.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBase.hpp" #include "linearAlgebra/utilities/SparsityPatternUtilities.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" @@ -284,6 +285,15 @@ void SolidMechanicsAugmentedLagrangianContact::setupSystem( DomainPartition & do GEOS_MARK_FUNCTION; + updateFractureGeometry( domain ); + + PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity ); +} + +void SolidMechanicsAugmentedLagrangianContact::updateFractureGeometry( DomainPartition & domain ) +{ + GEOS_MARK_FUNCTION; + // Recompute geometric quantities (face normals, areas) after mesh topology changes. // This is critical for distorted/non-axis-aligned meshes and after fracture events. forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, @@ -327,8 +337,6 @@ void SolidMechanicsAugmentedLagrangianContact::setupSystem( DomainPartition & do // Create the list of cell elements that they are enriched with bubble functions. createBubbleCellList( domain ); - - PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity ); } void SolidMechanicsAugmentedLagrangianContact::postInputInitialization() @@ -553,12 +561,15 @@ void SolidMechanicsAugmentedLagrangianContact::assembleContact( real64 const tim NodeManager const & nodeManager = mesh.getNodeManager(); FaceManager const & faceManager = mesh.getFaceManager(); + GEOS_LOG_RANK_0("[debug] \t \t Inside contact reaching for dof keys"); string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); string const & bubbleDofKey = dofManager.getKey( contact::totalBubbleDisplacement::key() ); + GEOS_LOG_RANK_0("[debug] \t \t Inside contact reaching for dof num"); arrayView1d< globalIndex const > const dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + GEOS_LOG_RANK_0(GEOS_FMT("[debug] \t \t Inside contact reaching unique fract reg: {}", m_fractureRegionNames[0])); string const & fractureRegionName = getUniqueFractureRegionName(); forFiniteElementOnStickFractureSubRegions( meshName, [&] ( string const &, @@ -569,6 +580,7 @@ void SolidMechanicsAugmentedLagrangianContact::assembleContact( real64 const tim if( m_simultaneous ) { + GEOS_LOG_RANK_0("[debug] \t \t Inside contact reaching for simultaneous-stick kernels"); solidMechanicsALMKernels::ALMSimultaneousFactory kernelFactory( dispDofNumber, bubbleDofNumber, dofManager.rankOffset(), @@ -577,6 +589,7 @@ void SolidMechanicsAugmentedLagrangianContact::assembleContact( real64 const tim dt, faceElementList ); + GEOS_LOG_RANK_0("[debug] \t \t Inside contact reaching for Coulomb kernels"); real64 maxTraction = finiteElement::interfaceBasedKernelApplication< parallelDevicePolicy< >, CoulombFriction >( mesh, fractureRegionName, faceElementList, @@ -621,6 +634,7 @@ void SolidMechanicsAugmentedLagrangianContact::assembleContact( real64 const tim if( m_simultaneous ) { + GEOS_LOG_RANK_0("[debug] \t \t Inside contact reaching for simultaneous-slip kernels"); solidMechanicsALMKernels::ALMSimultaneousFactory kernelFactory( dispDofNumber, bubbleDofNumber, dofManager.rankOffset(), @@ -684,6 +698,7 @@ void SolidMechanicsAugmentedLagrangianContact::assembleContact( real64 const tim real64 const gravityVectorData[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( gravityVector() ); + GEOS_LOG_RANK_0("[debug] \t \t Filling in bubble matrices"); solidMechanicsConformingContactKernels::FaceBubbleFactory kernelFactory( dispDofNumber, bubbleDofNumber, dofManager.rankOffset(), diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 405b2004d7d..fa013a3fb73 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -21,6 +21,8 @@ #ifndef GEOS_PHYSICSSOLVERS_CONTACT_SOLIDMECHANICSAUGMENTEDLAGRANGIANCONTACT_HPP_ #define GEOS_PHYSICSSOLVERS_CONTACT_SOLIDMECHANICSAUGMENTEDLAGRANGIANCONTACT_HPP_ +#include "common/format/Format.hpp" +#include "common/logger/Logger.hpp" #include "physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp" namespace geos @@ -29,6 +31,8 @@ namespace geos class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase { public: + static constexpr bool hasContactStabilization = true; + SolidMechanicsAugmentedLagrangianContact( const string & name, Group * const parent ); @@ -114,6 +118,10 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase integer configurationLoopIter ) override final; + //TODO (jafranc) - bubble should be reframe as a stabilization + string getStabilizationName() const { return "TODO:bubble"; } + bool hasStabilization() const { return false;} + /** * @brief Loop over the finite element type on the fracture subregions of meshName and apply callback. * @tparam LAMBDA The callback function type @@ -150,6 +158,7 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase bool const isStickState = true; + GEOS_LOG_RANK_0(GEOS_FMT("[debug] \t \t fetching list stick at {}",meshName)); stdMap< string, array1d< localIndex > > const & faceTypesToFaceElements = m_faceTypesToFaceElementsStick.at( meshName ); @@ -209,6 +218,17 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase */ void createBubbleCellList( DomainPartition & domain ) const; + /** + * @brief Recompute face and element geometric quantities (normals, areas, centers, volumes) after mesh + * topology changes, and reorder kf1 fracture nodes to match kf0 for the conforming contact kernels. + * @param domain The physical domain object + * + * This is called from setupSystem() for the standalone contact solver. Coupled solvers built on top of this + * contact solver (e.g. poromechanics conforming fractures) that do not route through setupSystem() must call + * this explicitly before assembling contact-dependent sparsity/terms. + */ + void updateFractureGeometry( DomainPartition & domain ); + private: /** diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp index 4fb89a1fa6d..f7ff4449c08 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp @@ -55,7 +55,6 @@ using namespace constitutive; using namespace dataRepository; using namespace fields; using namespace finiteElement; -const localIndex geos::SolidMechanicsLagrangeContact::m_maxFaceNodes = 11; SolidMechanicsLagrangeContact::SolidMechanicsLagrangeContact( const string & name, Group * const parent ): @@ -644,7 +643,7 @@ void SolidMechanicsLagrangeContact::assembleSystem( real64 const time, localMatrix, localRhs ); - assembleContact( domain, dofManager, localMatrix, localRhs ); + assembleContact( time, dt, domain, dofManager, localMatrix, localRhs ); // for sequential: add (fixed) pressure force contribution into residual (no derivatives) if( m_isFixedStressPoromechanicsUpdate || m_performStressInitialization ) @@ -658,7 +657,9 @@ void SolidMechanicsLagrangeContact::assembleSystem( real64 const time, } } -void SolidMechanicsLagrangeContact::assembleContact( DomainPartition & domain, +void SolidMechanicsLagrangeContact::assembleContact( real64 const GEOS_UNUSED_PARAM(time), + real64 const GEOS_UNUSED_PARAM(dt), + DomainPartition & domain, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) @@ -724,9 +725,9 @@ void SolidMechanicsLagrangeContact:: Nbar[ 2 ] = faceNormal[elemsToFaces[kfe][0]][2] - faceNormal[elemsToFaces[kfe][1]][2]; LvArray::tensorOps::normalize< 3 >( Nbar ); - globalIndex rowDOF[3 * m_maxFaceNodes]; - real64 nodeRHS[3 * m_maxFaceNodes]; - stackArray1d< real64, 3 * m_maxFaceNodes > dRdP( 3*m_maxFaceNodes ); + globalIndex rowDOF[3 * ContactSolverBase::m_maxFaceNodes]; + real64 nodeRHS[3 * ContactSolverBase::m_maxFaceNodes]; + stackArray1d< real64, 3 * ContactSolverBase::m_maxFaceNodes > dRdP( 3*ContactSolverBase::m_maxFaceNodes ); for( localIndex kf=0; kf<2; ++kf ) { @@ -1000,329 +1001,7 @@ void SolidMechanicsLagrangeContact::computeRotationMatrices( DomainPartition & d } ); } -void SolidMechanicsLagrangeContact::computeFaceIntegrals( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodesCoords, - localIndex const (&faceToNodes)[11], - localIndex const (&faceToEdges)[11], - localIndex const & numFaceVertices, - real64 const & faceArea, - real64 const (&faceCenter)[3], - real64 const (&faceNormal)[3], - arrayView2d< localIndex const > const & edgeToNodes, - real64 const & invCellDiameter, - real64 const (&cellCenter)[3], - stackArray1d< real64, FaceManager::maxFaceNodes() > & basisIntegrals, - real64 (& threeDMonomialIntegrals)[3] ) const -{ - GEOS_MARK_FUNCTION; - localIndex const MFN = m_maxFaceNodes; // Max number of face vertices. - basisIntegrals.resize( numFaceVertices ); - // Rotate the face. - // - compute rotation matrix. - real64 faceRotationMatrix[ 3 ][ 3 ]; - computationalGeometry::RotationMatrix_3D( faceNormal, faceRotationMatrix ); - // - below we compute the diameter, the rotated vertices and the rotated center. - real64 faceRotatedVertices[ MFN ][ 2 ]; - real64 faceDiameter = 0; - - for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) - { - // apply the transpose (that is the inverse) of the rotation matrix to face vertices. - // NOTE: - // the second and third rows of the transpose of the rotation matrix rotate on the 2D face. - faceRotatedVertices[numVertex][0] = - faceRotationMatrix[ 0 ][ 1 ]*nodesCoords( faceToNodes[ numVertex ], 0 ) + - faceRotationMatrix[ 1 ][ 1 ]*nodesCoords( faceToNodes[ numVertex ], 1 ) + - faceRotationMatrix[ 2 ][ 1 ]*nodesCoords( faceToNodes[ numVertex ], 2 ); - faceRotatedVertices[numVertex][1] = - faceRotationMatrix[ 0 ][ 2 ]*nodesCoords( faceToNodes[ numVertex ], 0 ) + - faceRotationMatrix[ 1 ][ 2 ]*nodesCoords( faceToNodes[ numVertex ], 1 ) + - faceRotationMatrix[ 2 ][ 2 ]*nodesCoords( faceToNodes[ numVertex ], 2 ); - } - - faceDiameter = computationalGeometry::computeDiameter< 2 >( faceRotatedVertices, - numFaceVertices ); - real64 const invFaceDiameter = 1.0/faceDiameter; - // - rotate the face centroid as done for the vertices. - real64 faceRotatedCentroid[2]; - faceRotatedCentroid[0] = - faceRotationMatrix[ 0 ][ 1 ]*faceCenter[0] + - faceRotationMatrix[ 1 ][ 1 ]*faceCenter[1] + - faceRotationMatrix[ 2 ][ 1 ]*faceCenter[2]; - faceRotatedCentroid[1] = - faceRotationMatrix[ 0 ][ 2 ]*faceCenter[0] + - faceRotationMatrix[ 1 ][ 2 ]*faceCenter[1] + - faceRotationMatrix[ 2 ][ 2 ]*faceCenter[2]; - // - compute edges' lengths, outward pointing normals and local edge-to-nodes map. - real64 edgeOutwardNormals[ MFN ][ 2 ]; - real64 edgeLengths[ MFN ]; - localIndex localEdgeToNodes[ MFN ][ 2 ]; - - for( localIndex numEdge = 0; numEdge < numFaceVertices; ++numEdge ) - { - if( edgeToNodes( faceToEdges[numEdge], 0 ) == faceToNodes[ numEdge ] ) - { - localEdgeToNodes[ numEdge ][ 0 ] = numEdge; - localEdgeToNodes[ numEdge ][ 1 ] = (numEdge+1)%numFaceVertices; - } - else - { - localEdgeToNodes[ numEdge ][ 0 ] = (numEdge+1)%numFaceVertices; - localEdgeToNodes[ numEdge ][ 1 ] = numEdge; - } - real64 edgeTangent[2]; - edgeTangent[0] = faceRotatedVertices[(numEdge+1)%numFaceVertices][0] - - faceRotatedVertices[numEdge][0]; - edgeTangent[1] = faceRotatedVertices[(numEdge+1)%numFaceVertices][1] - - faceRotatedVertices[numEdge][1]; - edgeOutwardNormals[numEdge][0] = edgeTangent[1]; - edgeOutwardNormals[numEdge][1] = -edgeTangent[0]; - real64 signTestVector[2]; - signTestVector[0] = faceRotatedVertices[numEdge][0] - faceRotatedCentroid[0]; - signTestVector[1] = faceRotatedVertices[numEdge][1] - faceRotatedCentroid[1]; - if( signTestVector[0]*edgeOutwardNormals[numEdge][0] + - signTestVector[1]*edgeOutwardNormals[numEdge][1] < 0 ) - { - edgeOutwardNormals[numEdge][0] = -edgeOutwardNormals[numEdge][0]; - edgeOutwardNormals[numEdge][1] = -edgeOutwardNormals[numEdge][1]; - } - edgeLengths[numEdge] = LvArray::math::sqrt< real64 >( edgeTangent[0]*edgeTangent[0] + - edgeTangent[1]*edgeTangent[1] ); - edgeOutwardNormals[numEdge][0] /= edgeLengths[numEdge]; - edgeOutwardNormals[numEdge][1] /= edgeLengths[numEdge]; - } - - // Compute boundary quadrature weights (also equal to the integrals of basis functions on the - // boundary). - real64 boundaryQuadratureWeights[ MFN ]; - for( localIndex numWeight = 0; numWeight < numFaceVertices; ++numWeight ) - boundaryQuadratureWeights[numWeight] = 0.0; - for( localIndex numEdge = 0; numEdge < numFaceVertices; ++numEdge ) - { - boundaryQuadratureWeights[ localEdgeToNodes[ numEdge ][ 0 ] ] += 0.5*edgeLengths[numEdge]; - boundaryQuadratureWeights[ localEdgeToNodes[ numEdge ][ 1 ] ] += 0.5*edgeLengths[numEdge]; - } - // Compute scaled monomials' integrals on edges. - real64 monomBoundaryIntegrals[3] = { 0.0 }; - for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) - { - monomBoundaryIntegrals[0] += boundaryQuadratureWeights[ numVertex ]; - monomBoundaryIntegrals[1] += (faceRotatedVertices[ numVertex ][ 0 ] - faceRotatedCentroid[0]) * - invFaceDiameter*boundaryQuadratureWeights[ numVertex ]; - monomBoundaryIntegrals[2] += (faceRotatedVertices[ numVertex ][ 1 ] - faceRotatedCentroid[1]) * - invFaceDiameter*boundaryQuadratureWeights[ numVertex ]; - } - - // Compute non constant 2D and 3D scaled monomials' integrals on the face. - real64 monomInternalIntegrals[2] = { 0.0 }; - for( localIndex numSubTriangle = 0; numSubTriangle < numFaceVertices; ++numSubTriangle ) - { - localIndex const nextVertex = (numSubTriangle+1)%numFaceVertices; - // - compute value of 2D monomials at the quadrature point on the sub-triangle (the - // barycenter). - // The result is ((v(0)+v(1)+faceCenter)/3 - faceCenter) / faceDiameter = - // = (v(0) + v(1) - 2*faceCenter)/(3*faceDiameter). - real64 monomialValues[2]; - for( localIndex i = 0; i < 2; ++i ) - { - monomialValues[i] = (faceRotatedVertices[numSubTriangle][i] + - faceRotatedVertices[nextVertex][i] - - 2.0*faceRotatedCentroid[i]) / (3.0*faceDiameter); - } - // compute value of 3D monomials at the quadrature point on the sub-triangle (the - // barycenter). The result is - // ((v(0) + v(1) + faceCenter)/3 - cellCenter)/cellDiameter. - real64 threeDMonomialValues[3]; - for( localIndex i = 0; i < 3; ++i ) - { - threeDMonomialValues[i] = ( (faceCenter[i] + - nodesCoords[faceToNodes[ numSubTriangle ]][i] + - nodesCoords[faceToNodes[ nextVertex ]][i]) / 3.0 - - cellCenter[i] ) * invCellDiameter; - } - // compute quadrature weight associated to the quadrature point (the area of the - // sub-triangle). - real64 edgesTangents[2][2]; // used to compute the area of the sub-triangle - for( localIndex i = 0; i < 2; ++i ) - { - edgesTangents[0][i] = faceRotatedVertices[numSubTriangle][i] - faceRotatedCentroid[i]; - } - for( localIndex i = 0; i < 2; ++i ) - { - edgesTangents[1][i] = faceRotatedVertices[nextVertex][i] - faceRotatedCentroid[i]; - } - real64 subTriangleArea = 0.5*LvArray::math::abs - ( edgesTangents[0][0]*edgesTangents[1][1] - - edgesTangents[0][1]*edgesTangents[1][0] ); - // compute the integrals on the sub-triangle and add it to the global integrals - for( localIndex i = 0; i < 2; ++i ) - { - monomInternalIntegrals[ i ] += monomialValues[ i ]*subTriangleArea; - } - for( localIndex i = 0; i < 3; ++i ) - { - // threeDMonomialIntegrals is assumed to be initialized to 0 by the caller - threeDMonomialIntegrals[ i ] += threeDMonomialValues[ i ]*subTriangleArea; - } - } - - // Compute integral of basis functions times normal derivative of monomials on the boundary. - real64 basisTimesMonomNormalDerBoundaryInt[ MFN ][ 2 ]; - for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) - { - for( localIndex i = 0; i < 2; ++i ) - { - basisTimesMonomNormalDerBoundaryInt[ numVertex ][ i ] = 0.0; - } - } - for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) - { - for( localIndex i = 0; i < 2; ++i ) - { - real64 thisEdgeIntTimesNormal_i = edgeOutwardNormals[numVertex][i]*edgeLengths[numVertex]; - basisTimesMonomNormalDerBoundaryInt[ localEdgeToNodes[ numVertex ][ 0 ] ][i] += thisEdgeIntTimesNormal_i; - basisTimesMonomNormalDerBoundaryInt[ localEdgeToNodes[ numVertex ][ 1 ] ][i] += thisEdgeIntTimesNormal_i; - } - } - for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) - { - for( localIndex i = 0; i < 2; ++i ) - { - basisTimesMonomNormalDerBoundaryInt[ numVertex ][ i ] *= 0.5*invFaceDiameter; - } - } - - // Compute integral mean of basis functions on this face. - real64 const invFaceArea = 1.0/faceArea; - real64 const monomialDerivativeInverse = (faceDiameter*faceDiameter)*invFaceArea; - for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) - { - real64 piNablaDofs[ 3 ]; - piNablaDofs[ 1 ] = monomialDerivativeInverse * - basisTimesMonomNormalDerBoundaryInt[ numVertex ][ 0 ]; - piNablaDofs[ 2 ] = monomialDerivativeInverse * - basisTimesMonomNormalDerBoundaryInt[ numVertex ][ 1 ]; - piNablaDofs[ 0 ] = (boundaryQuadratureWeights[ numVertex ] - - piNablaDofs[ 1 ]*monomBoundaryIntegrals[ 1 ] - - piNablaDofs[ 2 ]*monomBoundaryIntegrals[ 2 ])/monomBoundaryIntegrals[ 0 ]; - basisIntegrals[ numVertex ] = piNablaDofs[ 0 ]*faceArea + - (piNablaDofs[ 1 ]*monomInternalIntegrals[ 0 ] + - piNablaDofs[ 2 ]*monomInternalIntegrals[ 1 ]); - } -} - -void SolidMechanicsLagrangeContact::computeFaceNodalArea( localIndex const kf0, - arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, - ArrayOfArraysView< localIndex const > const & faceToNodeMap, - ArrayOfArraysView< localIndex const > const & faceToEdgeMap, - arrayView2d< localIndex const > const & edgeToNodeMap, - arrayView2d< real64 const > const faceCenters, - arrayView2d< real64 const > const faceNormals, - arrayView1d< real64 const > const faceAreas, - stackArray1d< real64, FaceManager::maxFaceNodes() > & basisIntegrals ) const -{ - GEOS_MARK_FUNCTION; - localIndex const TriangularPermutation[3] = { 0, 1, 2 }; - localIndex const QuadrilateralPermutation[4] = { 0, 1, 3, 2 }; - localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 ); - - basisIntegrals.resize( numNodesPerFace ); - for( localIndex a = 0; a < numNodesPerFace; ++a ) - { - basisIntegrals[a] = 0.0; - } - localIndex const * const permutation = ( numNodesPerFace == 3 ) ? TriangularPermutation : QuadrilateralPermutation; - if( numNodesPerFace == 3 ) - { - real64 xLocal[3][3]; - for( localIndex a = 0; a < numNodesPerFace; ++a ) - { - for( localIndex j = 0; j < 3; ++j ) - { - xLocal[a][j] = nodePosition[faceToNodeMap( kf0, permutation[a] )][j]; - } - } - real64 N[3]; - for( localIndex q=0; q 4 && numNodesPerFace <= m_maxFaceNodes ) - { - // we need to L2 projector based on VEM to approximate the quadrature weights - // we need to use extra geometry information to computing L2 projector - - localIndex const MFN = m_maxFaceNodes; // Max number of face vertices. - localIndex const faceIndex = kf0; - localIndex const numFaceNodes = faceToNodeMap[ faceIndex ].size(); - - // get the face center and normal. - real64 const faceArea = faceAreas[ faceIndex ]; - localIndex faceToNodes[ MFN ]; - localIndex faceToEdges[ MFN ]; - for( localIndex i = 0; i < numFaceNodes; ++i ) - { - faceToNodes[i] = faceToNodeMap[ faceIndex ][ i ]; - faceToEdges[i] = faceToEdgeMap[ faceIndex ][ i ]; - } - // - get outward face normal and center - real64 faceNormal[3] = { faceNormals[faceIndex][0], - faceNormals[faceIndex][1], - faceNormals[faceIndex][2] }; - real64 const faceCenter[3] { faceCenters[faceIndex][0], - faceCenters[faceIndex][1], - faceCenters[faceIndex][2] }; - // - compute integrals calling auxiliary method - real64 threeDMonomialIntegrals[3] = { 0.0 }; - real64 const invCellDiameter = 0.0; - real64 const cellCenter[3] { 0.0, 0.0, 0.0 }; - computeFaceIntegrals( nodePosition, - faceToNodes, - faceToEdges, - numFaceNodes, - faceArea, - faceCenter, - faceNormal, - edgeToNodeMap, - invCellDiameter, - cellCenter, - basisIntegrals, - threeDMonomialIntegrals ); - } - else - { - GEOS_ERROR( GEOS_FMT( "Face with {} nodes. Only triangles and quadrilaterals and PEBI prisms up to 11 sides are supported.", - numNodesPerFace ), - getDataContext() ); - } -} void SolidMechanicsLagrangeContact:: assembleForceResidualDerivativeWrtTraction( MeshLevel const & mesh, @@ -1368,9 +1047,9 @@ void SolidMechanicsLagrangeContact:: localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elemsToFaces[kfe][0] ); - globalIndex rowDOF[3 * m_maxFaceNodes]; // this needs to be changed when dealing with arbitrary element types - real64 nodeRHS[3 * m_maxFaceNodes]; - stackArray2d< real64, 3 * m_maxFaceNodes * 3 > dRdT( 3 * m_maxFaceNodes, 3 ); + globalIndex rowDOF[3 * ContactSolverBase::m_maxFaceNodes]; // this needs to be changed when dealing with arbitrary element types + real64 nodeRHS[3 * ContactSolverBase::m_maxFaceNodes]; + stackArray2d< real64, 3 * ContactSolverBase::m_maxFaceNodes * 3 > dRdT( 3 * ContactSolverBase::m_maxFaceNodes, 3 ); globalIndex colDOF[3]; for( localIndex i = 0; i < 3; ++i ) { @@ -1498,7 +1177,7 @@ void SolidMechanicsLagrangeContact:: if( ghostRank[kfe] < 0 ) { localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elemsToFaces[kfe][0] ); - globalIndex nodeDOF[2 * 3 * m_maxFaceNodes]; + globalIndex nodeDOF[2 * 3 * ContactSolverBase::m_maxFaceNodes]; globalIndex elemDOF[3]; for( localIndex i = 0; i < 3; ++i ) { @@ -1508,7 +1187,7 @@ void SolidMechanicsLagrangeContact:: real64 elemRHS[3] = {0.0, 0.0, 0.0}; real64 const Ja = area[kfe]; - stackArray2d< real64, 2 * 3 * m_maxFaceNodes * 3 > dRdU( 3, 2 * 3 * m_maxFaceNodes ); + stackArray2d< real64, 2 * 3 * ContactSolverBase::m_maxFaceNodes * 3 > dRdU( 3, 2 * 3 * ContactSolverBase::m_maxFaceNodes ); stackArray2d< real64, 3 * 3 > dRdT( 3, 3 ); switch( fractureState[kfe] ) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp index 674a08876c1..8c20a6be0c5 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp @@ -102,7 +102,9 @@ class SolidMechanicsLagrangeContact : public ContactSolverBase void updateState( DomainPartition & domain ) override final; - void assembleContact( DomainPartition & domain, + void assembleContact( real64 const time, + real64 const dt, + DomainPartition & domain, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ); @@ -142,32 +144,10 @@ class SolidMechanicsLagrangeContact : public ContactSolverBase void computeTolerances( DomainPartition & domain ) const; - void computeFaceNodalArea( localIndex const kf0, - arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, - ArrayOfArraysView< localIndex const > const & faceToNodeMap, - ArrayOfArraysView< localIndex const > const & faceToEdgeMap, - arrayView2d< localIndex const > const & edgeToNodeMap, - arrayView2d< real64 const > const faceCenters, - arrayView2d< real64 const > const faceNormals, - arrayView1d< real64 const > const faceAreas, - stackArray1d< real64, FaceManager::maxFaceNodes() > & nodalArea ) const; - - void computeFaceIntegrals( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodesCoords, - localIndex const (&faceToNodes)[11], - localIndex const (&faceToEdges)[11], - localIndex const & numFaceVertices, - real64 const & faceArea, - real64 const (&faceCenter)[3], - real64 const (&faceNormal)[3], - arrayView2d< localIndex const > const & edgeToNodes, - real64 const & invCellDiameter, - real64 const (&cellCenter)[3], - stackArray1d< real64, FaceManager::maxFaceNodes() > & basisIntegrals, - real64 ( &threeDMonomialIntegrals )[3] ) const; - real64 const machinePrecision = std::numeric_limits< real64 >::epsilon(); string getStabilizationName() const { return m_stabilizationName; } + bool hasStabilization() const { return true;} protected: @@ -186,8 +166,6 @@ class SolidMechanicsLagrangeContact : public ContactSolverBase real64 m_stabilizationScalingCoefficient = 1.0; - static const localIndex m_maxFaceNodes; // Maximum number of nodes on a contact face - void computeFaceDisplacementJump( DomainPartition & domain ); struct viewKeyStruct : ContactSolverBase::viewKeyStruct