21 #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSCONFORMINGFRACTURES_HPP_
22 #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSCONFORMINGFRACTURES_HPP_
30 #include "constitutive/solid/CoupledSolidBase.hpp"
31 #include "constitutive/contact/HydraulicApertureBase.hpp"
32 #include "constitutive/contact/HydraulicApertureRelationSelector.hpp"
40 template<
template<
typename,
typename >
class POROMECHANICS_BASE,
typename FLOW_SOLVER >
44 using Base = POROMECHANICS_BASE< FLOW_SOLVER, SolidMechanicsLagrangeContact >;
48 : Base( name, parent )
56 Base::setupCoupling( domain, dofManager );
60 fields::contact::traction::key(),
71 this->flowSolver()->setSparsityPattern( domain, dofManager, localMatrix, patternOriginal );
75 for(
localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow )
77 rowLengths[localRow] = patternOriginal.numNonZeros( localRow );
84 pattern.resizeFromRowCapacities< parallelHostPolicy >( patternOriginal.numRows(),
85 patternOriginal.numColumns(),
89 for(
localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow )
91 globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous();
92 pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow ));
101 virtual void assembleSystem(
real64 const time_n,
103 DomainPartition & domain,
104 DofManager
const & dofManager,
105 CRSMatrixView< real64, globalIndex const >
const & localMatrix,
106 arrayView1d< real64 >
const & localRhs )
override
111 this->solidMechanicsSolver()->synchronizeFractureState( domain );
121 this->flowSolver()->assembleHydrofracFluxTerms( time_n,
127 getDerivativeFluxResidual_dNormalJump() );
138 virtual void updateState( DomainPartition & domain )
override
143 Base::updateState( domain );
145 this->solidMechanicsSolver()->updateState( domain );
148 this->flowSolver()->prepareStencilWeights( domain );
150 updateHydraulicApertureAndFracturePermeability( domain );
153 this->flowSolver()->updateStencilWeights( domain );
170 integer const numComp = numFluidComponents();
172 this->forDiscretizationOnMeshTargets( domain.
getMeshBodies(), [&] (
string const &,
176 ElementRegionManager const & elemManager = mesh.getElemManager();
178 string const flowDofKey = dofManager.getKey( getFlowDofKey() );
180 globalIndex const rankOffset = dofManager.rankOffset();
182 NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
183 FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager();
184 FluxApproximationBase const & stabilizationMethod = fvManager.getFluxApproximation( this->solidMechanicsSolver()->getStabilizationName() );
186 stabilizationMethod.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil )
188 for( localIndex iconn=0; iconn<stencil.size(); ++iconn )
190 localIndex const numFluxElems = stencil.stencilSize( iconn );
191 typename SurfaceElementStencil::IndexContainerViewConstType const & seri = stencil.getElementRegionIndices();
192 typename SurfaceElementStencil::IndexContainerViewConstType const & sesri = stencil.getElementSubRegionIndices();
193 typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices();
195 FaceElementSubRegion const & elementSubRegion =
196 elemManager.getRegion( seri[iconn][0] ).getSubRegion< FaceElementSubRegion >( sesri[iconn][0] );
198 ArrayOfArraysView< localIndex const > const elemsToNodes = elementSubRegion.nodeList().toViewConst();
200 arrayView1d< globalIndex const > const faceElementDofNumber =
201 elementSubRegion.getReference< array1d< globalIndex > >( flowDofKey );
203 for( localIndex k0=0; k0<numFluxElems; ++k0 )
205 globalIndex const activeFlowDOF = faceElementDofNumber[sei[iconn][k0]];
206 globalIndex const rowNumber = activeFlowDOF - rankOffset;
208 if( rowNumber >= 0 && rowNumber < rowLengths.size() )
210 for( localIndex k1=0; k1<numFluxElems; ++k1 )
216 localIndex const numNodesPerElement = elemsToNodes[sei[iconn][k1]].size();
217 for( integer ic = 0; ic < numComp; ic++ )
219 rowLengths[rowNumber + ic] += 3*numNodesPerElement;
243 integer const numComp = numFluidComponents();
245 this->forDiscretizationOnMeshTargets( domain.
getMeshBodies(), [&] (
string const &,
249 FaceManager const & faceManager = mesh.getFaceManager();
250 NodeManager const & nodeManager = mesh.getNodeManager();
251 ElementRegionManager const & elemManager = mesh.getElemManager();
253 string const dispDofKey = dofManager.getKey( fields::solidMechanics::totalDisplacement::key() );
254 string const flowDofKey = dofManager.getKey( getFlowDofKey() );
256 arrayView1d< globalIndex const > const &
257 dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey );
258 ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst();
261 NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
262 FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager();
263 FluxApproximationBase const & fvDiscretization = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() );
265 SurfaceElementRegion const & fractureRegion =
266 elemManager.getRegion< SurfaceElementRegion >( this->solidMechanicsSolver()->getUniqueFractureRegionName() );
267 FaceElementSubRegion const & fractureSubRegion =
268 fractureRegion.getUniqueSubRegion< FaceElementSubRegion >();
270 GEOS_ERROR_IF( !fractureSubRegion.hasWrapper( fields::flow::pressure::key() ),
271 this->getDataContext() <<
": The fracture subregion must contain pressure field." );
273 arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst();
275 arrayView1d< globalIndex const > const &
276 flowDofNumber = fractureSubRegion.getReference< globalIndex_array >( flowDofKey );
278 globalIndex const rankOffset = dofManager.rankOffset();
280 fvDiscretization.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil )
282 forAll< serialPolicy >( stencil.size(), [=] ( localIndex const iconn )
284 localIndex const numFluxElems = stencil.stencilSize( iconn );
287 if( numFluxElems == 2 )
289 typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices();
292 for( localIndex kf = 0; kf < 2; ++kf )
297 localIndex const rowIndex = flowDofNumber[sei[iconn][1-kf]] - rankOffset;
299 if( rowIndex >= 0 && rowIndex < pattern.numRows() )
303 localIndex const fractureIndex = sei[iconn][kf];
306 localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elem2dToFaces[fractureIndex][0] );
309 for( localIndex kf1 = 0; kf1 < 2; ++kf1 )
311 localIndex const faceIndex = elem2dToFaces[fractureIndex][kf1];
314 for( localIndex a=0; a<numNodesPerFace; ++a )
316 for( localIndex i = 0; i < 3; ++i )
318 globalIndex const colIndex = dispDofNumber[faceToNodeMap( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i );
319 for( integer ic = 0; ic < numComp; ic++ )
321 pattern.insertNonZero( rowIndex + ic, colIndex );
345 integer const numComp = numFluidComponents();
347 this->forDiscretizationOnMeshTargets( domain.
getMeshBodies(), [&] (
string const &,
351 std::unique_ptr< CRSMatrix< real64, localIndex > > & derivativeFluxResidual_dAperture = getRefDerivativeFluxResidual_dAperture();
355 localIndex numRows = 0;
356 mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames,
357 [&]( localIndex const, FaceElementSubRegion const & subRegion )
359 numRows += subRegion.size();
362 localIndex numCol = numRows;
366 derivativeFluxResidual_dAperture = std::make_unique< CRSMatrix< real64, localIndex > >( numRows, numCol );
367 derivativeFluxResidual_dAperture->setName( this->getName() +
"/derivativeFluxResidual_dAperture" );
369 derivativeFluxResidual_dAperture->reserveNonZeros( localMatrix.numNonZeros() );
370 localIndex maxRowSize = -1;
371 for( localIndex row = 0; row < localMatrix.numRows(); ++row )
373 localIndex const rowSize = localMatrix.numNonZeros( row );
374 maxRowSize = maxRowSize > rowSize ? maxRowSize : rowSize;
377 for( localIndex row = 0; row < numRows; ++row )
379 derivativeFluxResidual_dAperture->reserveNonZeros( row, maxRowSize );
385 FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() );
389 for( localIndex iconn = 0; iconn < stencil.size(); ++iconn )
391 localIndex const numFluxElems = stencil.stencilSize( iconn );
392 typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices();
394 for( localIndex k0 = 0; k0 < numFluxElems; ++k0 )
396 for( localIndex k1 = 0; k1 < numFluxElems; ++k1 )
398 for( integer ic = 0; ic < numComp; ic++ )
400 derivativeFluxResidual_dAperture->insertNonZero( sei[iconn][k0] * numComp + ic, sei[iconn][k1], 0.0 );
420 Base::assembleElementBasedTerms( time_n, dt, domain, dofManager, localMatrix, localRhs );
423 this->forDiscretizationOnMeshTargets( domain.
getMeshBodies(), [&] (
string const &,
427 mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const,
428 FaceElementSubRegion const & subRegion )
430 this->flowSolver()->accumulationAssemblyLaunch( dofManager, subRegion, localMatrix, localRhs );
434 this->solidMechanicsSolver()->assembleContact( domain, dofManager, localMatrix, localRhs );
446 this->forDiscretizationOnMeshTargets( domain.
getMeshBodies(), [&] (
string const &,
451 assembleForceResidualDerivativeWrtPressure( mesh, regionNames, dofManager, localMatrix, localRhs );
452 assembleFluidMassResidualDerivativeWrtDisplacement( mesh, regionNames, dofManager, localMatrix, localRhs );
456 void assembleForceResidualDerivativeWrtPressure(
MeshLevel const & mesh,
476 string const & dispDofKey = dofManager.
getKey( fields::solidMechanics::totalDisplacement::key() );
477 string const & flowDofKey = dofManager.
getKey( getFlowDofKey() );
495 forAll< serialPolicy >( subRegion.size(), [=](
localIndex const kfe )
497 localIndex const kf0 = elemsToFaces[kfe][0];
498 localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 );
501 Nbar[ 0 ] = faceNormal[elemsToFaces[kfe][0]][0] - faceNormal[elemsToFaces[kfe][1]][0];
502 Nbar[ 1 ] = faceNormal[elemsToFaces[kfe][0]][1] - faceNormal[elemsToFaces[kfe][1]][1];
503 Nbar[ 2 ] = faceNormal[elemsToFaces[kfe][0]][2] - faceNormal[elemsToFaces[kfe][1]][2];
504 LvArray::tensorOps::normalize< 3 >( Nbar );
505 globalIndex rowDOF[3 * m_maxFaceNodes];
506 real64 nodeRHS[3 * m_maxFaceNodes];
507 stackArray1d< real64, 3 * m_maxFaceNodes > dRdP( 3*m_maxFaceNodes );
508 globalIndex colDOF[1];
509 colDOF[0] = flowDofNumber[kfe];
511 for( localIndex kf=0; kf<2; ++kf )
513 localIndex const faceIndex = elemsToFaces[kfe][kf];
516 stackArray1d< real64, FaceManager::maxFaceNodes() > nodalArea;
517 this->solidMechanicsSolver()->computeFaceNodalArea( elemsToFaces[kfe][kf],
526 for( localIndex a=0; a<numNodesPerFace; ++a )
528 real64 const nodalForceMag = -( pressure[kfe] ) * nodalArea[a];
529 real64 globalNodalForce[ 3 ];
530 LvArray::tensorOps::scaledCopy< 3 >( globalNodalForce, Nbar, nodalForceMag );
532 for( localIndex i=0; i<3; ++i )
534 rowDOF[3*a+i] = dispDofNumber[faceToNodeMap( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i );
536 nodeRHS[3*a+i] = +globalNodalForce[i] * pow( -1, kf );
539 dRdP( 3*a+i ) = -nodalArea[a] * Nbar[i] * pow( -1, kf );
543 for( localIndex idof = 0; idof < numNodesPerFace * 3; ++idof )
545 localIndex const localRow = LvArray::integerConversion< localIndex >( rowDOF[idof] - rankOffset );
547 if( localRow >= 0 && localRow < localMatrix.numRows() )
549 localMatrix.addToRow< parallelHostAtomic >( localRow,
553 RAJA::atomicAdd( parallelHostAtomic{}, &localRhs[localRow], nodeRHS[idof] );
561 virtual void assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel
const & mesh,
563 DofManager
const & dofManager,
564 CRSMatrixView< real64, globalIndex const >
const & localMatrix,
565 arrayView1d< real64 >
const & localRhs ) = 0;
572 if( solverType ==
static_cast< integer >( Base::SolverType::SolidMechanics )
573 && !this->m_performStressInitialization )
576 this->flowSolver()->prepareStencilWeights( domain );
578 updateHydraulicApertureAndFracturePermeability( domain );
581 this->flowSolver()->updateStencilWeights( domain );
584 Base::mapSolutionBetweenSolvers( domain, solverType );
587 void updateHydraulicApertureAndFracturePermeability(
DomainPartition & domain )
589 using namespace constitutive;
591 this->forDiscretizationOnMeshTargets( domain.
getMeshBodies(), [&] (
string const &,
595 ElementRegionManager & elemManager = mesh.getElemManager();
597 elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames,
598 [&]( localIndex const,
599 FaceElementSubRegion & subRegion )
601 arrayView2d< real64 const > const dispJump = subRegion.getField< fields::contact::dispJump >();
602 arrayView1d< real64 const > const area = subRegion.getElementArea();
603 arrayView1d< real64 const > const volume = subRegion.getElementVolume();
604 arrayView2d< real64 const > const fractureTraction = subRegion.getField< fields::contact::traction >();
605 arrayView1d< real64 const > const pressure = subRegion.getField< fields::flow::pressure >();
606 arrayView1d< real64 const > const oldHydraulicAperture = subRegion.getField< fields::flow::aperture0 >();
608 arrayView1d< real64 > const aperture = subRegion.getElementAperture();
609 arrayView1d< real64 > const hydraulicAperture = subRegion.getField< fields::flow::hydraulicAperture >();
610 arrayView1d< real64 > const deltaVolume = subRegion.getField< fields::flow::deltaVolume >();
611 arrayView1d< integer > const & fractureState = subRegion.getField< fields::contact::fractureState >();
613 string const porousSolidName = subRegion.getReference< string >( FlowSolverBase::viewKeyStruct::solidNamesString() );
614 CoupledSolidBase & porousSolid = subRegion.getConstitutiveModel< CoupledSolidBase >( porousSolidName );
616 string const & hydraulicApertureRelationName = subRegion.template getReference< string >( viewKeyStruct::hydraulicApertureRelationNameString() );
617 HydraulicApertureBase const & hydraulicApertureModel = this->template getConstitutiveModel< HydraulicApertureBase >( subRegion, hydraulicApertureRelationName );
619 constitutiveUpdatePassThru( hydraulicApertureModel, [&] ( auto & castedHydraulicAperture )
621 using HydraulicApertureType = TYPEOFREF( castedHydraulicAperture );
622 typename HydraulicApertureType::KernelWrapper hydraulicApertureWrapper = castedHydraulicAperture.createKernelWrapper();
624 ConstitutivePassThru< CompressibleSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid )
626 typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousMaterialWrapper = castedPorousSolid.createKernelUpdates();
628 poromechanicsFracturesKernels::StateUpdateKernel::
629 launch< parallelDevicePolicy<> >( subRegion.size(),
630 porousMaterialWrapper,
631 hydraulicApertureWrapper,
638 oldHydraulicAperture,
649 std::unique_ptr< CRSMatrix< real64, localIndex > > & getRefDerivativeFluxResidual_dAperture()
651 return m_derivativeFluxResidual_dAperture;
654 CRSMatrixView< real64, localIndex const > getDerivativeFluxResidual_dNormalJump()
656 return m_derivativeFluxResidual_dAperture->toViewConstSizes();
659 CRSMatrixView< real64 const, localIndex const > getDerivativeFluxResidual_dNormalJump()
const
661 return m_derivativeFluxResidual_dAperture->toViewConst();
664 virtual string getFlowDofKey()
const = 0;
666 virtual integer numFluidComponents()
const = 0;
673 std::unique_ptr< CRSMatrix< real64, localIndex > > m_derivativeFluxResidual_dAperture;
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
#define GEOS_UNUSED_PARAM(X)
Mark an unused argument and silence compiler warnings.
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
The DoFManager is responsible for allocating global dofs, constructing sparsity patterns,...
globalIndex rankOffset(string const &fieldName) const
void addCoupling(string const &rowFieldName, string const &colFieldName, Connector connectivity, stdVector< FieldSupport > const ®ions={}, bool symmetric=true)
Add coupling between two fields.
@ Elem
connectivity is element (like in finite elements)
string const & getKey(string const &fieldName) const
Return the key used to record the field in the DofManager.
Partition of the decomposed physical domain. It also manages the connexion information to its neighbo...
Group const & getMeshBodies() const
Get the mesh bodies, const version.
NumericalMethodsManager const & getNumericalMethodManager() const
This class provides an interface to ObjectManagerBase in order to manage edge data.
NodeMapType & nodeList()
Get a node list.
The ElementRegionManager class provides an interface to ObjectManagerBase in order to manage ElementR...
void forElementSubRegions(LAMBDA &&lambda)
This function is used to launch kernel function over the element subregions of all the subregion type...
The FaceManager class provides an interface to ObjectManagerBase in order to manage face data.
array2d< real64 > & faceNormal()
Get a mutable accessor to an array containing all the face normals.
array1d< real64 > & faceArea()
Get a mutable accessor to an array containing all the face area.
array2d< real64 > & faceCenter()
Get a mutable accessor to an array containing all the face center.
NodeMapType & nodeList()
Get a mutable accessor to a map containing the list of each nodes for each faces.
EdgeMapType & edgeList()
Get a mutable accessor to a map containing the list of each edges for each faces.
void forStencils(MeshLevel const &mesh, LAMBDA &&lambda) const
Call a user-provided function for the each stencil according to the provided TYPE.
Class facilitating the representation of a multi-level discretization of a MeshBody.
NodeManager const & getNodeManager() const
Get the node manager.
FaceManager const & getFaceManager() const
Get the face manager.
ElementRegionManager const & getElemManager() const
Get the element region manager.
EdgeManager const & getEdgeManager() const
Get the edge manager.
The NodeManager class provides an interface to ObjectManagerBase in order to manage node data.
Provides management of the interior stencil points for a face elements when using Two-Point flux appr...
GEOS_DECLTYPE_AUTO_RETURN getReference(LOOKUP_TYPE const &lookup) const
Look up a wrapper and get reference to wrapped object.
array2d< real64, nodes::REFERENCE_POSITION_PERM > & referencePosition()
Get the mutable reference position array. This table will contain all the node coordinates.
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
stdVector< string > string_array
A 1-dimensional array of geos::string types.
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
array1d< globalIndex > globalIndex_array
A 1-dimensional array of geos::globalIndex types.
LvArray::CRSMatrix< T, COL_INDEX, INDEX_TYPE, LvArray::ChaiBuffer > CRSMatrix
Alias for CRS Matrix class.
LvArray::SparsityPatternView< COL_INDEX, INDEX_TYPE const, LvArray::ChaiBuffer > SparsityPatternView
Alias for Sparsity pattern View.
LvArray::SparsityPattern< COL_INDEX, INDEX_TYPE, LvArray::ChaiBuffer > SparsityPattern
Alias for Sparsity pattern class.
LvArray::ArrayOfArraysView< T, INDEX_TYPE const, CONST_SIZES, LvArray::ChaiBuffer > ArrayOfArraysView
View of array of variable-sized arrays. See LvArray::ArrayOfArraysView for details.
double real64
64-bit floating point type.
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
LvArray::CRSMatrixView< T, COL_INDEX, INDEX_TYPE const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
int integer
Signed integer type.
Array< T, 1 > array1d
Alias for 1D array.