GEOS
PoromechanicsConformingFractures.hpp
Go to the documentation of this file.
1 /*
2  * ------------------------------------------------------------------------------------------------------------
3  * SPDX-License-Identifier: LGPL-2.1-only
4  *
5  * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
6  * Copyright (c) 2018-2024 TotalEnergies
7  * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
8  * Copyright (c) 2023-2024 Chevron
9  * Copyright (c) 2019- GEOS/GEOSX Contributors
10  * All rights reserved
11  *
12  * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
13  * ------------------------------------------------------------------------------------------------------------
14  */
15 
21 #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSCONFORMINGFRACTURES_HPP_
22 #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSCONFORMINGFRACTURES_HPP_
23 
30 #include "constitutive/solid/CoupledSolidBase.hpp"
31 #include "constitutive/contact/HydraulicApertureBase.hpp"
32 #include "constitutive/contact/HydraulicApertureRelationSelector.hpp"
34 #include "common/DataTypes.hpp"
35 #include "mesh/DomainPartition.hpp"
36 
37 namespace geos
38 {
39 
40 template< template< typename, typename > class POROMECHANICS_BASE, typename FLOW_SOLVER >
41 class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, SolidMechanicsLagrangeContact >
42 {
43 public:
44  using Base = POROMECHANICS_BASE< FLOW_SOLVER, SolidMechanicsLagrangeContact >;
45 
46  PoromechanicsConformingFractures( const string & name,
47  dataRepository::Group * const parent )
48  : Base( name, parent )
49  {}
50 
51  virtual void setupCoupling( DomainPartition const & domain,
52  DofManager & dofManager ) const override
53  {
55  // 1. Poromechanical coupling in the bulk
56  Base::setupCoupling( domain, dofManager );
57 
58  // 2. Traction - pressure coupling in the fracture
59  dofManager.addCoupling( getFlowDofKey(),
60  fields::contact::traction::key(),
62  }
63 
64  virtual void setSparsityPattern( DomainPartition & domain,
65  DofManager & dofManager,
67  SparsityPattern< globalIndex > & pattern ) override
68  {
69  // start with the flow solver sparsity pattern (it could be reservoir + wells)
70  SparsityPattern< globalIndex > patternOriginal;
71  this->flowSolver()->setSparsityPattern( domain, dofManager, localMatrix, patternOriginal );
72 
73  // Get the original row lengths (diagonal blocks only)
74  array1d< localIndex > rowLengths( patternOriginal.numRows());
75  for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow )
76  {
77  rowLengths[localRow] = patternOriginal.numNonZeros( localRow );
78  }
79 
80  // Add the number of nonzeros induced by coupling
81  addTransmissibilityCouplingNNZ( domain, dofManager, rowLengths.toView());
82 
83  // Create a new pattern with enough capacity for coupled matrix
84  pattern.resizeFromRowCapacities< parallelHostPolicy >( patternOriginal.numRows(),
85  patternOriginal.numColumns(),
86  rowLengths.data());
87 
88  // Copy the original nonzeros
89  for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow )
90  {
91  globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous();
92  pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow ));
93  }
94 
95  // Add the nonzeros from coupling
96  addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView());
97 
98  setUpDflux_dApertureMatrix( domain, dofManager, localMatrix );
99  }
100 
101  virtual void assembleSystem( real64 const time_n,
102  real64 const dt,
103  DomainPartition & domain,
104  DofManager const & dofManager,
105  CRSMatrixView< real64, globalIndex const > const & localMatrix,
106  arrayView1d< real64 > const & localRhs ) override
107  {
108 
110 
111  this->solidMechanicsSolver()->synchronizeFractureState( domain );
112 
114  dt,
115  domain,
116  dofManager,
117  localMatrix,
118  localRhs );
119 
120  // Assemble fluxes 3D/2D and get dFluidResidualDAperture
121  this->flowSolver()->assembleHydrofracFluxTerms( time_n,
122  dt,
123  domain,
124  dofManager,
125  localMatrix,
126  localRhs,
127  getDerivativeFluxResidual_dNormalJump() );
128 
129  // This step must occur after the fluxes are assembled because that's when DerivativeFluxResidual_dAperture is filled.
130  assembleCouplingTerms( time_n,
131  dt,
132  domain,
133  dofManager,
134  localMatrix,
135  localRhs );
136  }
137 
138  virtual void updateState( DomainPartition & domain ) override
139  {
141 
142  // call base poromechanics update
143  Base::updateState( domain );
144  // need to call solid mechanics update separately to compute face displacement jump
145  this->solidMechanicsSolver()->updateState( domain );
146 
147  // remove the contribution of the hydraulic aperture from the stencil weights
148  this->flowSolver()->prepareStencilWeights( domain );
149 
150  updateHydraulicApertureAndFracturePermeability( domain );
151 
152  // update the stencil weights using the updated hydraulic aperture
153  this->flowSolver()->updateStencilWeights( domain );
154  }
155 
156 protected:
157 
165  DofManager const & dofManager,
166  arrayView1d< localIndex > const & rowLengths ) const
167  {
169 
170  integer const numComp = numFluidComponents();
171 
172  this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, // meshBodyName,
173  MeshLevel const & mesh,
174  string_array const & ) // regionNames
175  {
176  ElementRegionManager const & elemManager = mesh.getElemManager();
177 
178  string const flowDofKey = dofManager.getKey( getFlowDofKey() );
179 
180  globalIndex const rankOffset = dofManager.rankOffset();
181 
182  NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
183  FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager();
184  FluxApproximationBase const & stabilizationMethod = fvManager.getFluxApproximation( this->solidMechanicsSolver()->getStabilizationName() );
185 
186  stabilizationMethod.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil )
187  {
188  for( localIndex iconn=0; iconn<stencil.size(); ++iconn )
189  {
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();
194 
195  FaceElementSubRegion const & elementSubRegion =
196  elemManager.getRegion( seri[iconn][0] ).getSubRegion< FaceElementSubRegion >( sesri[iconn][0] );
197 
198  ArrayOfArraysView< localIndex const > const elemsToNodes = elementSubRegion.nodeList().toViewConst();
199 
200  arrayView1d< globalIndex const > const faceElementDofNumber =
201  elementSubRegion.getReference< array1d< globalIndex > >( flowDofKey );
202 
203  for( localIndex k0=0; k0<numFluxElems; ++k0 )
204  {
205  globalIndex const activeFlowDOF = faceElementDofNumber[sei[iconn][k0]];
206  globalIndex const rowNumber = activeFlowDOF - rankOffset;
207 
208  if( rowNumber >= 0 && rowNumber < rowLengths.size() )
209  {
210  for( localIndex k1=0; k1<numFluxElems; ++k1 )
211  {
212  // The coupling with the nodal displacements of the cell itself has already been added by the dofManager
213  // so we only add the coupling with the nodal displacements of the neighbors.
214  if( k1 != k0 )
215  {
216  localIndex const numNodesPerElement = elemsToNodes[sei[iconn][k1]].size();
217  for( integer ic = 0; ic < numComp; ic++ )
218  {
219  rowLengths[rowNumber + ic] += 3*numNodesPerElement;
220  }
221  }
222  }
223  }
224  }
225  }
226  } );
227  } );
228  }
229 
238  DofManager const & dofManager,
239  SparsityPatternView< globalIndex > const & pattern ) const
240  {
242 
243  integer const numComp = numFluidComponents();
244 
245  this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
246  MeshLevel const & mesh,
247  string_array const & )
248  {
249  FaceManager const & faceManager = mesh.getFaceManager();
250  NodeManager const & nodeManager = mesh.getNodeManager();
251  ElementRegionManager const & elemManager = mesh.getElemManager();
252 
253  string const dispDofKey = dofManager.getKey( fields::solidMechanics::totalDisplacement::key() );
254  string const flowDofKey = dofManager.getKey( getFlowDofKey() );
255 
256  arrayView1d< globalIndex const > const &
257  dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey );
258  ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst();
259 
260  // Get the finite volume method used to compute the stabilization
261  NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
262  FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager();
263  FluxApproximationBase const & fvDiscretization = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() );
264 
265  SurfaceElementRegion const & fractureRegion =
266  elemManager.getRegion< SurfaceElementRegion >( this->solidMechanicsSolver()->getUniqueFractureRegionName() );
267  FaceElementSubRegion const & fractureSubRegion =
268  fractureRegion.getUniqueSubRegion< FaceElementSubRegion >();
269 
270  GEOS_ERROR_IF( !fractureSubRegion.hasWrapper( fields::flow::pressure::key() ),
271  this->getDataContext() << ": The fracture subregion must contain pressure field." );
272 
273  arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst();
274 
275  arrayView1d< globalIndex const > const &
276  flowDofNumber = fractureSubRegion.getReference< globalIndex_array >( flowDofKey );
277 
278  globalIndex const rankOffset = dofManager.rankOffset();
279 
280  fvDiscretization.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil )
281  {
282  forAll< serialPolicy >( stencil.size(), [=] ( localIndex const iconn )
283  {
284  localIndex const numFluxElems = stencil.stencilSize( iconn );
285 
286  // A fracture connector has to be an edge shared by two faces
287  if( numFluxElems == 2 )
288  {
289  typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices();
290 
291  // First index: face element. Second index: node
292  for( localIndex kf = 0; kf < 2; ++kf )
293  {
294  // Set row DOF index
295  // Note that the 1-kf index is intentional, as this is coupling the pressure of one face cell
296  // to the nodes of the adjacent cell
297  localIndex const rowIndex = flowDofNumber[sei[iconn][1-kf]] - rankOffset;
298 
299  if( rowIndex >= 0 && rowIndex < pattern.numRows() )
300  {
301 
302  // Get fracture, face and region/subregion/element indices (for elements on both sides)
303  localIndex const fractureIndex = sei[iconn][kf];
304 
305  // Get the number of nodes
306  localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elem2dToFaces[fractureIndex][0] );
307 
308  // Loop over the two sides of each fracture element
309  for( localIndex kf1 = 0; kf1 < 2; ++kf1 )
310  {
311  localIndex const faceIndex = elem2dToFaces[fractureIndex][kf1];
312 
313  // Save the list of DOF associated with nodes
314  for( localIndex a=0; a<numNodesPerFace; ++a )
315  {
316  for( localIndex i = 0; i < 3; ++i )
317  {
318  globalIndex const colIndex = dispDofNumber[faceToNodeMap( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i );
319  for( integer ic = 0; ic < numComp; ic++ )
320  {
321  pattern.insertNonZero( rowIndex + ic, colIndex );
322  }
323  }
324  }
325  }
326  }
327  }
328  }
329  } );
330  } );
331  } );
332  }
333 
342  DofManager const & GEOS_UNUSED_PARAM( dofManager ),
343  CRSMatrix< real64, globalIndex > & localMatrix )
344  {
345  integer const numComp = numFluidComponents();
346 
347  this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
348  MeshLevel const & mesh,
349  string_array const & regionNames )
350  {
351  std::unique_ptr< CRSMatrix< real64, localIndex > > & derivativeFluxResidual_dAperture = getRefDerivativeFluxResidual_dAperture();
352 
353  {
354  // calculate number of fracture elements
355  localIndex numRows = 0;
356  mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames,
357  [&]( localIndex const, FaceElementSubRegion const & subRegion )
358  {
359  numRows += subRegion.size();
360  } );
361  // number of columns (derivatives) = number of fracture elements
362  localIndex numCol = numRows;
363  // number of rows (equations) = number of fracture elements * number of components
364  numRows *= numComp;
365 
366  derivativeFluxResidual_dAperture = std::make_unique< CRSMatrix< real64, localIndex > >( numRows, numCol );
367  derivativeFluxResidual_dAperture->setName( this->getName() + "/derivativeFluxResidual_dAperture" );
368 
369  derivativeFluxResidual_dAperture->reserveNonZeros( localMatrix.numNonZeros() );
370  localIndex maxRowSize = -1;
371  for( localIndex row = 0; row < localMatrix.numRows(); ++row )
372  {
373  localIndex const rowSize = localMatrix.numNonZeros( row );
374  maxRowSize = maxRowSize > rowSize ? maxRowSize : rowSize;
375  }
376  // TODO This is way too much. The With the full system rowSize is not a good estimate for this.
377  for( localIndex row = 0; row < numRows; ++row )
378  {
379  derivativeFluxResidual_dAperture->reserveNonZeros( row, maxRowSize );
380  }
381  }
382 
383  NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
384  FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager();
385  FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() );
386 
387  fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil )
388  {
389  for( localIndex iconn = 0; iconn < stencil.size(); ++iconn )
390  {
391  localIndex const numFluxElems = stencil.stencilSize( iconn );
392  typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices();
393 
394  for( localIndex k0 = 0; k0 < numFluxElems; ++k0 )
395  {
396  for( localIndex k1 = 0; k1 < numFluxElems; ++k1 )
397  {
398  for( integer ic = 0; ic < numComp; ic++ )
399  {
400  derivativeFluxResidual_dAperture->insertNonZero( sei[iconn][k0] * numComp + ic, sei[iconn][k1], 0.0 );
401  }
402  }
403  }
404  }
405  } );
406  } );
407  }
408 
410  real64 const dt,
411  DomainPartition & domain,
412  DofManager const & dofManager,
413  CRSMatrixView< real64, globalIndex const > const & localMatrix,
414  arrayView1d< real64 > const & localRhs )
415  {
416  GEOS_UNUSED_VAR( time_n, dt );
417 
419 
420  Base::assembleElementBasedTerms( time_n, dt, domain, dofManager, localMatrix, localRhs );
421 
422  // Flow accumulation for fractures
423  this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
424  MeshLevel & mesh,
425  string_array const & regionNames )
426  {
427  mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const,
428  FaceElementSubRegion const & subRegion )
429  {
430  this->flowSolver()->accumulationAssemblyLaunch( dofManager, subRegion, localMatrix, localRhs );
431  } );
432  } );
433 
434  this->solidMechanicsSolver()->assembleContact( domain, dofManager, localMatrix, localRhs );
435  }
436 
437  virtual void assembleCouplingTerms( real64 const time_n,
438  real64 const dt,
439  DomainPartition const & domain,
440  DofManager const & dofManager,
441  CRSMatrixView< real64, globalIndex const > const & localMatrix,
442  arrayView1d< real64 > const & localRhs ) override
443  {
444  GEOS_UNUSED_VAR( time_n, dt );
445  // These 2 steps need to occur after the fluxes are assembled because that's when DerivativeFluxResidual_dAperture is filled.
446  this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
447  MeshLevel const & mesh,
448  string_array const & regionNames )
449  {
451  assembleForceResidualDerivativeWrtPressure( mesh, regionNames, dofManager, localMatrix, localRhs );
452  assembleFluidMassResidualDerivativeWrtDisplacement( mesh, regionNames, dofManager, localMatrix, localRhs );
453  } );
454  }
455 
456  void assembleForceResidualDerivativeWrtPressure( MeshLevel const & mesh,
457  string_array const & regionNames,
458  DofManager const & dofManager,
459  CRSMatrixView< real64, globalIndex const > const & localMatrix,
460  arrayView1d< real64 > const & localRhs )
461  {
463 
464  FaceManager const & faceManager = mesh.getFaceManager();
465  NodeManager const & nodeManager = mesh.getNodeManager();
466  EdgeManager const & edgeManager = mesh.getEdgeManager();
467  ElementRegionManager const & elemManager = mesh.getElemManager();
468 
469  ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst();
470  ArrayOfArraysView< localIndex const > const faceToEdgeMap = faceManager.edgeList().toViewConst();
471  arrayView2d< localIndex const > const & edgeToNodeMap = edgeManager.nodeList().toViewConst();
472  arrayView2d< real64 const > faceCenters = faceManager.faceCenter();
473  arrayView2d< real64 const > const & faceNormal = faceManager.faceNormal();
474  arrayView1d< real64 const > faceAreas = faceManager.faceArea();
475 
476  string const & dispDofKey = dofManager.getKey( fields::solidMechanics::totalDisplacement::key() );
477  string const & flowDofKey = dofManager.getKey( getFlowDofKey() );
478 
480  dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey );
481  globalIndex const rankOffset = dofManager.rankOffset();
482 
483  // Get the coordinates for all nodes
485 
486  elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames,
487  [&]( localIndex const,
488  FaceElementSubRegion const & subRegion )
489  {
491  flowDofNumber = subRegion.getReference< globalIndex_array >( flowDofKey );
492  arrayView1d< real64 const > const & pressure = subRegion.getReference< array1d< real64 > >( fields::flow::pressure::key() );
493  arrayView2d< localIndex const > const & elemsToFaces = subRegion.faceList().toViewConst();
494 
495  forAll< serialPolicy >( subRegion.size(), [=]( localIndex const kfe )
496  {
497  localIndex const kf0 = elemsToFaces[kfe][0];
498  localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 );
499 
500  real64 Nbar[3];
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]; // this needs to be changed when dealing with arbitrary element types
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]; // pressure is always first
510 
511  for( localIndex kf=0; kf<2; ++kf )
512  {
513  localIndex const faceIndex = elemsToFaces[kfe][kf];
514 
515  // Compute local area contribution for each node
516  stackArray1d< real64, FaceManager::maxFaceNodes() > nodalArea;
517  this->solidMechanicsSolver()->computeFaceNodalArea( elemsToFaces[kfe][kf],
518  nodePosition,
519  faceToNodeMap,
520  faceToEdgeMap,
521  edgeToNodeMap,
522  faceCenters,
523  faceNormal,
524  faceAreas,
525  nodalArea );
526  for( localIndex a=0; a<numNodesPerFace; ++a )
527  {
528  real64 const nodalForceMag = -( pressure[kfe] ) * nodalArea[a];
529  real64 globalNodalForce[ 3 ];
530  LvArray::tensorOps::scaledCopy< 3 >( globalNodalForce, Nbar, nodalForceMag );
531 
532  for( localIndex i=0; i<3; ++i )
533  {
534  rowDOF[3*a+i] = dispDofNumber[faceToNodeMap( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i );
535  // Opposite sign w.r.t. theory because of minus sign in stiffness matrix definition (K < 0)
536  nodeRHS[3*a+i] = +globalNodalForce[i] * pow( -1, kf );
537 
538  // Opposite sign w.r.t. theory because of minus sign in stiffness matrix definition (K < 0)
539  dRdP( 3*a+i ) = -nodalArea[a] * Nbar[i] * pow( -1, kf );
540  }
541  }
542 
543  for( localIndex idof = 0; idof < numNodesPerFace * 3; ++idof )
544  {
545  localIndex const localRow = LvArray::integerConversion< localIndex >( rowDOF[idof] - rankOffset );
546 
547  if( localRow >= 0 && localRow < localMatrix.numRows() )
548  {
549  localMatrix.addToRow< parallelHostAtomic >( localRow,
550  colDOF,
551  &dRdP[idof],
552  1 );
553  RAJA::atomicAdd( parallelHostAtomic{}, &localRhs[localRow], nodeRHS[idof] );
554  }
555  }
556  }
557  } );
558  } );
559  }
560 
561  virtual void assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh,
562  string_array const & regionNames,
563  DofManager const & dofManager,
564  CRSMatrixView< real64, globalIndex const > const & localMatrix,
565  arrayView1d< real64 > const & localRhs ) = 0;
566 
567  virtual void mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) override
568  {
570 
572  if( solverType == static_cast< integer >( Base::SolverType::SolidMechanics )
573  && !this->m_performStressInitialization ) // do not update during poromechanics initialization
574  {
575  // remove the contribution of the hydraulic aperture from the stencil weights
576  this->flowSolver()->prepareStencilWeights( domain );
577 
578  updateHydraulicApertureAndFracturePermeability( domain );
579 
580  // update the stencil weights using the updated hydraulic aperture
581  this->flowSolver()->updateStencilWeights( domain );
582  }
583 
584  Base::mapSolutionBetweenSolvers( domain, solverType );
585  }
586 
587  void updateHydraulicApertureAndFracturePermeability( DomainPartition & domain )
588  {
589  using namespace constitutive;
590 
591  this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
592  MeshLevel & mesh,
593  string_array const & regionNames )
594  {
595  ElementRegionManager & elemManager = mesh.getElemManager();
596 
597  elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames,
598  [&]( localIndex const,
599  FaceElementSubRegion & subRegion )
600  {
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 >();
607 
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 >();
612 
613  string const porousSolidName = subRegion.getReference< string >( FlowSolverBase::viewKeyStruct::solidNamesString() );
614  CoupledSolidBase & porousSolid = subRegion.getConstitutiveModel< CoupledSolidBase >( porousSolidName );
615 
616  string const & hydraulicApertureRelationName = subRegion.template getReference< string >( viewKeyStruct::hydraulicApertureRelationNameString() );
617  HydraulicApertureBase const & hydraulicApertureModel = this->template getConstitutiveModel< HydraulicApertureBase >( subRegion, hydraulicApertureRelationName );
618 
619  constitutiveUpdatePassThru( hydraulicApertureModel, [&] ( auto & castedHydraulicAperture )
620  {
621  using HydraulicApertureType = TYPEOFREF( castedHydraulicAperture );
622  typename HydraulicApertureType::KernelWrapper hydraulicApertureWrapper = castedHydraulicAperture.createKernelWrapper();
623 
624  ConstitutivePassThru< CompressibleSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid )
625  {
626  typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousMaterialWrapper = castedPorousSolid.createKernelUpdates();
627 
628  poromechanicsFracturesKernels::StateUpdateKernel::
629  launch< parallelDevicePolicy<> >( subRegion.size(),
630  porousMaterialWrapper,
631  hydraulicApertureWrapper,
632  dispJump,
633  pressure,
634  area,
635  volume,
636  deltaVolume,
637  aperture,
638  oldHydraulicAperture,
639  hydraulicAperture,
640  fractureTraction,
641  fractureState );
642 
643  } );
644  } );
645  } );
646  } );
647  }
648 
649  std::unique_ptr< CRSMatrix< real64, localIndex > > & getRefDerivativeFluxResidual_dAperture()
650  {
651  return m_derivativeFluxResidual_dAperture;
652  }
653 
654  CRSMatrixView< real64, localIndex const > getDerivativeFluxResidual_dNormalJump()
655  {
656  return m_derivativeFluxResidual_dAperture->toViewConstSizes();
657  }
658 
659  CRSMatrixView< real64 const, localIndex const > getDerivativeFluxResidual_dNormalJump() const
660  {
661  return m_derivativeFluxResidual_dAperture->toViewConst();
662  }
663 
664  virtual string getFlowDofKey() const = 0;
665 
666  virtual integer numFluidComponents() const = 0;
667 
668  struct viewKeyStruct : public Base::viewKeyStruct
669  {};
670 
671  static const localIndex m_maxFaceNodes = 11; // Maximum number of nodes on a contact face
672 
673  std::unique_ptr< CRSMatrix< real64, localIndex > > m_derivativeFluxResidual_dAperture;
674 
675 };
676 
677 }
678 
679 #endif //GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSCONFORMINGFRACTURES_HPP_
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
Definition: GeosxMacros.hpp:84
#define GEOS_UNUSED_PARAM(X)
Mark an unused argument and silence compiler warnings.
Definition: GeosxMacros.hpp:72
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
The DoFManager is responsible for allocating global dofs, constructing sparsity patterns,...
Definition: DofManager.hpp:45
globalIndex rankOffset(string const &fieldName) const
void addCoupling(string const &rowFieldName, string const &colFieldName, Connector connectivity, stdVector< FieldSupport > const &regions={}, 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.
Definition: EdgeManager.hpp:43
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.
Definition: FaceManager.hpp:44
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.
Definition: MeshLevel.hpp:42
NodeManager const & getNodeManager() const
Get the node manager.
Definition: MeshLevel.hpp:155
FaceManager const & getFaceManager() const
Get the face manager.
Definition: MeshLevel.hpp:194
ElementRegionManager const & getElemManager() const
Get the element region manager.
Definition: MeshLevel.hpp:207
EdgeManager const & getEdgeManager() const
Get the edge manager.
Definition: MeshLevel.hpp:181
The NodeManager class provides an interface to ObjectManagerBase in order to manage node data.
Definition: NodeManager.hpp:46
void setUpDflux_dApertureMatrix(DomainPartition &domain, DofManager const &GEOS_UNUSED_PARAM(dofManager), CRSMatrix< real64, globalIndex > &localMatrix)
Set up the Dflux_dApertureMatrix object.
void addTransmissibilityCouplingPattern(DomainPartition const &domain, DofManager const &dofManager, SparsityPatternView< globalIndex > const &pattern) const
Set up the Dflux_dApertureMatrix object.
void assembleElementBasedContributions(real64 const time_n, real64 const dt, DomainPartition &domain, DofManager const &dofManager, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs)
void addTransmissibilityCouplingNNZ(DomainPartition const &domain, DofManager const &dofManager, arrayView1d< localIndex > const &rowLengths) const
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
virtual void mapSolutionBetweenSolvers(DomainPartition &domain, integer const solverType) override
virtual void setupCoupling(DomainPartition const &domain, DofManager &dofManager) const override
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.
Definition: Group.hpp:1275
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.
Definition: DataTypes.hpp:179
stdVector< string > string_array
A 1-dimensional array of geos::string types.
Definition: DataTypes.hpp:361
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
Definition: DataTypes.hpp:87
array1d< globalIndex > globalIndex_array
A 1-dimensional array of geos::globalIndex types.
Definition: DataTypes.hpp:370
LvArray::CRSMatrix< T, COL_INDEX, INDEX_TYPE, LvArray::ChaiBuffer > CRSMatrix
Alias for CRS Matrix class.
Definition: DataTypes.hpp:305
LvArray::SparsityPatternView< COL_INDEX, INDEX_TYPE const, LvArray::ChaiBuffer > SparsityPatternView
Alias for Sparsity pattern View.
Definition: DataTypes.hpp:301
LvArray::SparsityPattern< COL_INDEX, INDEX_TYPE, LvArray::ChaiBuffer > SparsityPattern
Alias for Sparsity pattern class.
Definition: DataTypes.hpp:297
LvArray::ArrayOfArraysView< T, INDEX_TYPE const, CONST_SIZES, LvArray::ChaiBuffer > ArrayOfArraysView
View of array of variable-sized arrays. See LvArray::ArrayOfArraysView for details.
Definition: DataTypes.hpp:285
double real64
64-bit floating point type.
Definition: DataTypes.hpp:98
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:84
LvArray::CRSMatrixView< T, COL_INDEX, INDEX_TYPE const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
Definition: DataTypes.hpp:309
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:195
int integer
Signed integer type.
Definition: DataTypes.hpp:81
Array< T, 1 > array1d
Alias for 1D array.
Definition: DataTypes.hpp:175