GEOS
SurfaceGenerator.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 
19 #ifndef GEOS_PHYSICSSOLVERS_SURFACEGENERATION_SURFACEGENERATOR_HPP_
20 #define GEOS_PHYSICSSOLVERS_SURFACEGENERATION_SURFACEGENERATOR_HPP_
21 
22 #include "mesh/mpiCommunications/NeighborCommunicator.hpp"
23 #include "physicsSolvers/PhysicsSolverBase.hpp"
24 #include "mesh/DomainPartition.hpp"
25 #include "mesh/CellElementSubRegion.hpp"
26 #include "ParallelTopologyChange.hpp"
27 
28 namespace geos
29 {
30 
31 class SpatialPartition;
32 class NodeManager;
33 class EdgeManager;
34 class FaceManager;
35 class ElementRegionManager;
36 class ElementRegionBase;
37 
42 struct ElemLocComparator
43 {
44  bool operator()( std::pair< CellElementSubRegion const *, localIndex > const & a,
45  std::pair< CellElementSubRegion const *, localIndex > const & b ) const
46  {
47  if( a.first != b.first )
48  {
49  localIndex const regionA = a.first->getParent().getParent().getIndexInParent();
50  localIndex const regionB = b.first->getParent().getParent().getIndexInParent();
51  if( regionA != regionB )
52  return regionA < regionB;
53  return a.first->getIndexInParent() < b.first->getIndexInParent();
54  }
55  return a.second < b.second;
56  }
57 };
58 
60 using ElemLocMapType = stdMap< std::pair< CellElementSubRegion const *, localIndex >, int, ElemLocComparator >;
61 
68 class SurfaceGenerator : public PhysicsSolverBase
69 {
70 public:
71  SurfaceGenerator( const string & name,
72  Group * const parent );
73  ~SurfaceGenerator() override;
74 
75 
76  static string catalogName() { return "SurfaceGenerator"; }
80  string getCatalogName() const override { return catalogName(); }
81 
82  virtual void registerDataOnMesh( Group & MeshBody ) override final;
83 
91  virtual bool execute( real64 const time_n,
92  real64 const dt,
93  integer const cycleNumber,
94  integer const eventCounter,
95  real64 const eventProgress,
96  DomainPartition & domain ) override;
97 
98  virtual real64 solverStep( real64 const & time_n,
99  real64 const & dt,
100  integer const cycleNumber,
101  DomainPartition & domain ) override;
102 
105  inline string const getFractureRegionName() const { return m_fractureRegionName; }
106 
107  void postInputInitialization() override final;
108 
109 protected:
110 
111  virtual void initializePostInitialConditionsPreSubGroups() override final;
112  virtual void postRestartInitialization() override final;
113 
114 private:
115 
116  int separationDriver( DomainPartition & domain,
117  MeshLevel & mesh,
118  stdVector< NeighborCommunicator > & neighbors,
119  int const tileColor,
120  int const numTileColors,
121  const bool prefrac,
122  const real64 time_np1 );
123 
129  void assignNewGlobalIndicesSerial( ObjectManagerBase & object,
130  std::set< localIndex > const & indexList );
131 
132 
138  void
139  assignNewGlobalIndicesSerial( ElementRegionManager & elementManager,
140  map< std::pair< localIndex, localIndex >, std::set< localIndex > > const & indexList );
141 
151  void identifyRupturedFaces( DomainPartition const & domain,
152  NodeManager & nodeManager,
153  EdgeManager & edgeManager,
154  FaceManager & faceManager,
155  ElementRegionManager const & elementManager,
156  bool const prefrac );
157 
170  real64 calculateEdgeSif( DomainPartition const & domain,
171  localIndex const edgeID,
172  localIndex & trailFaceID,
173  NodeManager const & nodeManager,
174  EdgeManager & edgeManager,
175  FaceManager const & faceManager,
176  ElementRegionManager const & elementManager,
177  real64 ( &vecTipNorm )[3],
178  real64 ( &vecTip )[3] );
179 
188  void calculateNodeAndFaceSif( DomainPartition const & domain,
189  NodeManager & nodeManager,
190  EdgeManager const & edgeManager,
191  FaceManager & faceManager,
192  ElementRegionManager const & elementManager );
193 
208  int calculateElementForcesOnEdge( DomainPartition const & domain,
209  localIndex const edgeID,
210  real64 edgeLength,
211  localIndex_array & nodeIndices,
212  NodeManager const & nodeManager,
213  EdgeManager const & edgeManager,
214  ElementRegionManager const & elementManager,
215  real64 ( &vecTipNorm )[3],
216  real64 ( &fNode )[3],
217  real64 & GdivBeta,
218  bool threeNodesPinched,
219  bool calculatef_u );
220 
234  void markRuptureFaceFromEdge( localIndex const edgeID,
235  localIndex const & trailFaceID,
236  NodeManager const & nodeManager,
237  EdgeManager const & edgeManager,
238  FaceManager & faceManager,
239  ElementRegionManager const & elementManager,
240  real64 ( &vecTipNorm )[3],
241  real64 ( &vecTip )[3],
242  ModifiedObjectLists & modifiedObjects,
243  int const edgeMode );
244 
254  void markRuptureFaceFromNode( localIndex const nodeIndex,
255  NodeManager const & nodeManager,
256  EdgeManager const & edgeManager,
257  FaceManager & faceManager,
258  ElementRegionManager const & elementManager,
259  ModifiedObjectLists & modifiedObjects );
260 
270  void postUpdateRuptureStates( NodeManager const & nodeManager,
271  EdgeManager const & edgeManager,
272  FaceManager const & faceManager,
273  ElementRegionManager const & elementManager,
274  stdVector< std::set< localIndex > > & nodesToRupturedFaces,
275  stdVector< std::set< localIndex > > & edgesToRupturedFaces );
276 
292  stdVector< std::set< localIndex > > groupRupturedFacesIntoSets(
293  stdVector< std::set< localIndex > > const & nodesToRupturedFaces,
294  NodeManager const & nodeManager,
295  FaceManager const & faceManager ) const;
296 
308  void buildFilteredRuptureMaps(
309  std::set< localIndex > const & fractureSetFaces,
310  NodeManager const & nodeManager,
311  EdgeManager const & edgeManager,
312  FaceManager const & faceManager,
313  stdVector< std::set< localIndex > > & nodesToRupturedFaces,
314  stdVector< std::set< localIndex > > & edgesToRupturedFaces ) const;
315 
322  int checkOrphanElement( ElementRegionManager const & elementManager,
323  FaceManager const & faceManager,
324  localIndex iFace );
325 
335  int checkEdgeSplitability( localIndex const edgeID,
336  NodeManager const & nodeManager,
337  FaceManager const & faceManager,
338  EdgeManager const & edgeManager,
339  const bool prefrac );
340 
341 
342 // void UpdatePathCheckingArrays();
343 
358  bool processNode( const localIndex nodeID,
359  real64 const time,
360  NodeManager & nodeManager,
361  EdgeManager & edgeManager,
362  FaceManager & faceManager,
363  ElementRegionManager & elemManager,
364  stdVector< std::set< localIndex > > & nodesToRupturedFaces,
365  stdVector< std::set< localIndex > > & edgesToRupturedFaces,
366  ElementRegionManager & elementManager,
367  ModifiedObjectLists & modifiedObjects,
368  const bool prefrac );
369 
385  bool findFracturePlanes( localIndex const nodeID,
386  NodeManager const & nodeManager,
387  EdgeManager const & edgeManager,
388  FaceManager const & faceManager,
389  ElementRegionManager const & elemManager,
390  stdVector< std::set< localIndex > > const & nodesToRupturedFaces,
391  stdVector< std::set< localIndex > > const & edgesToRupturedFaces,
392  std::set< localIndex > & separationPathFaces,
393  map< localIndex, int > & edgeLocations,
394  map< localIndex, int > & faceLocations,
395  ElemLocMapType & elemLocations );
396 
397 
413  void performFracture( localIndex const nodeID,
414  real64 const time_np1,
415  NodeManager & nodeManager,
416  EdgeManager & edgeManager,
417  FaceManager & faceManager,
418  ElementRegionManager & elementManager,
419  ModifiedObjectLists & modifiedObjects,
420  stdVector< std::set< localIndex > > & nodesToRupturedFaces,
421  stdVector< std::set< localIndex > > & edgesToRupturedFaces,
422  std::set< localIndex > const & separationPathFaces,
423  map< localIndex, int > const & edgeLocations,
424  map< localIndex, int > const & faceLocations,
425  ElemLocMapType const & elemLocations );
426 
427  void mapConsistencyCheck( localIndex const nodeID,
428  NodeManager const & nodeManager,
429  EdgeManager const & edgeManager,
430  FaceManager const & faceManager,
431  ElementRegionManager const & elementManager,
432  ElemLocMapType const & elemLocations );
433 
465  bool assignLocationsBFS( std::set< localIndex > const & separationPathFaces,
466  ElementRegionManager const & elemManager,
467  FaceManager const & faceManager,
468  stdVector< std::pair< CellElementSubRegion const *, localIndex > > const & nodeToElementMaps,
469  map< localIndex, std::pair< localIndex, localIndex > > const & localFacesToEdges,
470  map< localIndex, int > & edgeLocations,
471  map< localIndex, int > & faceLocations,
472  ElemLocMapType & elemLocations );
473 
482  real64 calculateKinkAngle( localIndex const nodeID,
483  NodeManager const & nodeManager,
484  EdgeManager const & edgeManager,
485  FaceManager const & faceManager );
486 
495  void calculateKinkAngles( FaceManager const & faceManager,
496  EdgeManager & edgeManager,
497  NodeManager const & nodeManager,
498  ModifiedObjectLists const & modifiedObjects,
499  bool const prefrac );
500 
505  void synchronizeTipSets( FaceManager & faceManager,
506  EdgeManager & edgeManager,
507  NodeManager & nodeManager,
508  ModifiedObjectLists & receivedObjects );
509 
510 
511 // void setDegreeFromCrackTip( NodeManager & nodeManager,
512 // FaceManager & faceManager );
513 
522  real64 minimumToughnessOnEdge( localIndex const edgeID,
523  NodeManager const & nodeManager,
524  EdgeManager const & edgeManager,
525  FaceManager const & faceManager );
526 
535  real64 minimumToughnessOnNode( localIndex const nodeID,
536  NodeManager const & nodeManager,
537  EdgeManager const & edgeManager,
538  FaceManager const & faceManager );
539 
540 
541  real64 calculateRuptureRate( SurfaceElementRegion & faceElementRegion );
542 
557  bool isGhostBoundaryFace( localIndex const faceIndex ) const
558  {
559  return m_originalFacesToElemIndex( faceIndex, 1 ) == -1;
560  }
561 
562  real64 scalingToughness( R1Tensor const fractureOrigin,
563  real64 const (&faceCenter)[3],
564  real64 const initialRockToughness,
565  real64 const toughnessScalingFactor );
566 
570  struct viewKeyStruct : PhysicsSolverBase::viewKeyStruct
571  {
572  constexpr static char const * failCriterionString() { return "failCriterion"; }
573  constexpr static char const * solidMaterialNameString() { return "solidMaterialNames"; }
574  constexpr static char const * fExternalString() { return "fExternal"; }
575  constexpr static char const * tipNodesString() { return "tipNodes"; }
576  constexpr static char const * tipEdgesString() { return "tipEdges"; }
577  constexpr static char const * tipFacesString() { return "tipFaces"; }
578  constexpr static char const * trailingFacesString() { return "trailingFaces"; }
579  constexpr static char const * fractureRegionNameString() { return "fractureRegion"; }
580  constexpr static char const * mpiCommOrderString() { return "mpiCommOrder"; }
581  constexpr static char const * isPoroelasticString() {return "isPoroelastic";}
582 
583  //TODO: rock toughness should be a material parameter, and we need to make rock toughness to KIC a constitutive
584  // relation.
585  constexpr static char const * initialRockToughnessString() { return "initialRockToughness"; }
586  constexpr static char const * toughnessScalingFactorString() { return "toughnessScalingFactor"; }
587  //TODO: fracture origin can be obtained from the initial fracture geometry
588  constexpr static char const * fractureOriginString() { return "fractureOrigin"; }
589 
590 // //TODO: Once the node-based SIF criterion becomes mature and robust, remove the edge-based criterion.
591  constexpr static char const * nodeBasedSIFString() { return "nodeBasedSIF"; }
592  };
593 
594 
595 private:
596 
597  constexpr static real64 m_nonRuptureTime = 1e9;
598 
600  integer m_failCriterion=1;
601 
602  array1d< localIndex > m_solidMaterialFullIndex;
603 
604  int m_nodeBasedSIF;
605 
606  int m_isPoroelastic;
607 
608  real64 m_initialRockToughness;
609 
610  real64 m_toughnessScalingFactor;
611 
612  R1Tensor m_fractureOrigin;
613 
614  // Flag for consistent communication ordering
615  int m_mpiCommOrder;
616 
618  SortedArray< localIndex > m_separableFaceSet;
619 
621  ArrayOfSets< localIndex > m_originalNodetoFaces;
622 
624  ArrayOfSets< localIndex > m_originalNodetoEdges;
625 
627  ArrayOfArrays< localIndex > m_originalFaceToEdges;
628 
630  array1d< SortedArray< localIndex > > m_usedFacesForNode;
631 
633  array2d< localIndex > m_originalFacesToElemRegion;
634 
636  array2d< localIndex > m_originalFacesToElemSubRegion;
637 
639  array2d< localIndex > m_originalFacesToElemIndex;
640 
642  string m_fractureRegionName;
643 
644  SortedArray< localIndex > m_tipNodes;
645 
646  SortedArray< localIndex > m_tipEdges;
647 
648  SortedArray< localIndex > m_tipFaces;
649 
650  SortedArray< localIndex > m_trailingFaces;
651 
652  SortedArray< localIndex > m_faceElemsRupturedThisSolve;
653 
656  localIndex m_numFractureElementsBefore = 0;
657 
658 };
659 
660 } /* namespace geos */
661 
662 #endif /* GEOS_PHYSICSSOLVERS_SURFACEGENERATION_SURFACEGENERATOR_HPP_ */
geos::DomainPartition
Partition of the decomposed physical domain. It also manages the connexion information to its neighbo...
Definition: DomainPartition.hpp:53
geos::EdgeManager
This class provides an interface to ObjectManagerBase in order to manage edge data.
Definition: EdgeManager.hpp:43
geos::ElementRegionManager
The ElementRegionManager class provides an interface to ObjectManagerBase in order to manage ElementR...
Definition: ElementRegionManager.hpp:45
geos::FaceManager
The FaceManager class provides an interface to ObjectManagerBase in order to manage face data.
Definition: FaceManager.hpp:44
geos::MeshBody
The class is used to manage mesh body.
Definition: MeshBody.hpp:36
geos::MeshLevel
Class facilitating the representation of a multi-level discretization of a MeshBody.
Definition: MeshLevel.hpp:42
geos::NodeManager
The NodeManager class provides an interface to ObjectManagerBase in order to manage node data.
Definition: NodeManager.hpp:46
geos::ObjectManagerBase
The ObjectManagerBase is the base object of all object managers in the mesh data hierachy.
Definition: ObjectManagerBase.hpp:34
geos::PhysicsSolverBase
Base class for all physics solvers.
Definition: PhysicsSolverBase.hpp:52
geos::SurfaceGenerator::postInputInitialization
void postInputInitialization() override final
geos::SurfaceGenerator::initializePostInitialConditionsPreSubGroups
virtual void initializePostInitialConditionsPreSubGroups() override final
Called by InitializePostInitialConditions() prior to initializing sub-Groups.
geos::SurfaceGenerator::postRestartInitialization
virtual void postRestartInitialization() override final
Performs initialization required after reading from a restart file.
geos::SurfaceGenerator::getCatalogName
string getCatalogName() const override
Definition: SurfaceGenerator.hpp:80
geos::SurfaceGenerator::registerDataOnMesh
virtual void registerDataOnMesh(Group &MeshBody) override final
Register wrappers that contain data on the mesh objects.
geos::mapBase
Base template for ordered and unordered maps.
Definition: StdContainerWrappers.hpp:335
geos::SurfaceGenerator::solverStep
virtual real64 solverStep(real64 const &time_n, real64 const &dt, integer const cycleNumber, DomainPartition &domain) override
entry function to perform a solver step
geos::SurfaceGenerator::execute
virtual bool execute(real64 const time_n, real64 const dt, integer const cycleNumber, integer const eventCounter, real64 const eventProgress, DomainPartition &domain) override
Main extension point of executable targets.
geos::localIndex_array
array1d< localIndex > localIndex_array
A 1-dimensional array of geos::localIndex types.
Definition: DataTypes.hpp:367
geos::set
std::set< T > set
A set of local indices.
Definition: DataTypes.hpp:262
geos::real64
double real64
64-bit floating point type.
Definition: DataTypes.hpp:98
geos::localIndex
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:84
geos::stdMap
internal::StdMapWrapper< std::map< Key, T, Compare, Allocator >, USE_STD_CONTAINER_BOUNDS_CHECKING > stdMap
Definition: StdContainerWrappers.hpp:266
geos::R1Tensor
Tensor< real64, 3 > R1Tensor
Alias for a local (stack-based) rank-1 tensor type.
Definition: DataTypes.hpp:165
geos::ElemLocMapType
stdMap< std::pair< CellElementSubRegion const *, localIndex >, int, ElemLocComparator > ElemLocMapType
Type alias for the element location map with deterministic ordering.
Definition: SurfaceGenerator.hpp:60
geos::integer
int integer
Signed integer type.
Definition: DataTypes.hpp:81
geos::stdVector
internal::StdVectorWrapper< T, Allocator, USE_STD_CONTAINER_BOUNDS_CHECKING > stdVector
Definition: StdContainerWrappers.hpp:162
geos::ElemLocComparator
Comparator for (CellElementSubRegion const *, localIndex) pairs that provides deterministic ordering ...
Definition: SurfaceGenerator.hpp:43
geos::PhysicsSolverBase::viewKeyStruct
Structure to hold scoped key names.
Definition: PhysicsSolverBase.hpp:727
Generated by doxygen 1.9.1