GEOS
HydrofractureSolverKernels.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_HYDROFRACTURESOLVERKERNELS_HPP_
22 #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_HYDROFRACTURESOLVERKERNELS_HPP_
23 
24 #include "physicsSolvers/solidMechanics/contact/FractureState.hpp"
25 namespace geos
26 {
27 
28 namespace hydrofractureSolverKernels
29 {
30 
31 struct DeformationUpdateKernel
32 {
33 
34  template< typename POLICY, typename HYDRAULICAPERTURE_WRAPPER, typename POROUS_WRAPPER >
35  static std::tuple< double, double, double, double, double, double >
36  launch( localIndex const size,
37  HYDRAULICAPERTURE_WRAPPER const & hydraulicApertureWrapper,
38  POROUS_WRAPPER const & porousMaterialWrapper,
39  arrayView2d< real64 const, nodes::TOTAL_DISPLACEMENT_USD > const & u,
40  arrayView2d< real64 const > const & faceNormal,
41  ArrayOfArraysView< localIndex const > const & faceToNodeMap,
42  arrayView2d< localIndex const > const & elemsToFaces,
43  arrayView1d< real64 const > const & area,
44  arrayView1d< real64 const > const & volume,
45  arrayView1d< real64 > const & deltaVolume,
46  arrayView1d< real64 > const & aperture,
47  arrayView1d< integer > const & fractureState,
48  arrayView1d< real64 > const & hydraulicAperture
49 #ifdef GEOS_USE_SEPARATION_COEFFICIENT
50  ,
51  arrayView1d< real64 const > const & apertureAtFailure,
52  arrayView1d< real64 > const & separationCoeff,
53  arrayView1d< real64 > const & dSeparationCoeff_dAper,
54  arrayView1d< real64 const > const & separationCoeff0
55 #endif
56  )
57  {
58 
59  RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxApertureChange( 0.0 );
60  RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxHydraulicApertureChange( 0.0 );
61  RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > minAperture( 1e10 );
62  RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxAperture( -1e10 );
63  RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > minHydraulicAperture( 1e10 );
64  RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxHydraulicAperture( -1e10 );
65 
66  forAll< POLICY >( size,
67  [=] GEOS_HOST_DEVICE ( localIndex const kfe ) mutable
68  {
69 
70  localIndex const kf0 = elemsToFaces[kfe][0];
71  localIndex const kf1 = elemsToFaces[kfe][1];
72  localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 );
73  real64 temp[ 3 ] = { 0 };
74  for( localIndex a=0; a<numNodesPerFace; ++a )
75  {
76  LvArray::tensorOps::add< 3 >( temp, u[ faceToNodeMap( kf0, a ) ] );
77  LvArray::tensorOps::subtract< 3 >( temp, u[ faceToNodeMap( kf1, a ) ] );
78  }
79 
80  // TODO this needs a proper contact based strategy for aperture
81  real64 const normalJump = -LvArray::tensorOps::AiBi< 3 >( temp, faceNormal[kf0] ) / numNodesPerFace;
82  maxApertureChange.max( std::fabs( normalJump - aperture[kfe] ));
83  aperture[kfe] = normalJump;
84  minAperture.min( aperture[kfe] );
85  maxAperture.max( aperture[kfe] );
86 
87  fractureState[kfe] = normalJump <= 0.0 ? fields::contact::FractureState::Stick : fields::contact::FractureState::Open;
88 
89  real64 normalTraction = 0.0;
90  real64 dHydraulicAperture_dNormalTraction = 0.0;
91  real64 dHydraulicAperture_dNormalJump = 0.0;
92  real64 const newHydraulicAperture = hydraulicApertureWrapper.computeHydraulicAperture( aperture[kfe],
93  normalTraction,
94  fractureState[kfe],
95  dHydraulicAperture_dNormalJump,
96  dHydraulicAperture_dNormalTraction );
97 
98  maxHydraulicApertureChange.max( std::fabs( newHydraulicAperture - hydraulicAperture[kfe] ));
99  real64 const oldHydraulicAperture = hydraulicAperture[kfe];
100  hydraulicAperture[kfe] = newHydraulicAperture;
101  minHydraulicAperture.min( hydraulicAperture[kfe] );
102  maxHydraulicAperture.max( hydraulicAperture[kfe] );
103 
104  real64 const jump[3] = { normalJump, 0.0, 0.0 };
105  real64 const traction[3] = {0.0, 0.0, 0.0};
106 
107  porousMaterialWrapper.updateStateFromPressureApertureJumpAndTraction( kfe, 0, 0.0,
108  oldHydraulicAperture, newHydraulicAperture,
109  dHydraulicAperture_dNormalJump,
110  jump, traction );
111 
112 #ifdef GEOS_USE_SEPARATION_COEFFICIENT
113  real64 const s = aperture[kfe] / apertureAtFailure[kfe];
114  if( separationCoeff0[kfe]<1.0 && s>separationCoeff0[kfe] )
115  {
116  if( s >= 1.0 )
117  {
118  separationCoeff[kfe] = 1.0;
119  dSeparationCoeff_dAper[kfe] = 0.0;
120  }
121  else
122  {
123  separationCoeff[kfe] = s;
124  dSeparationCoeff_dAper[kfe] = 1.0/apertureAtFailure[kfe];
125  }
126  }
127 #endif
128  deltaVolume[kfe] = hydraulicAperture[kfe] * area[kfe] - volume[kfe];
129  } );
130 
131  return std::make_tuple( maxApertureChange.get(), maxHydraulicApertureChange.get(), minAperture.get(), maxAperture.get(), minHydraulicAperture.get(), maxHydraulicAperture.get() );
132  }
133 };
134 
135 struct FluidMassResidualDerivativeAssemblyKernel
136 {
137  template< typename HYDRAULICAPERTURE_WRAPPER >
138  GEOS_HOST_DEVICE
139  inline
140  static void
141  computeAccumulationDerivative( HYDRAULICAPERTURE_WRAPPER const & hydraulicApertureWrapper,
142  localIndex const numNodesPerFace,
143  arraySlice1d< localIndex const > const elemsToFaces,
144  ArrayOfArraysView< localIndex const > const faceToNodeMap,
145  arrayView1d< globalIndex const > const dispDofNumber,
146  real64 const (&Nbar)[ 3 ],
147  real64 const & area,
148  real64 const & aperture,
149  real64 const & dens,
150  integer const & fractureState,
151  globalIndex (& nodeDOF)[8 * 3],
152  arraySlice1d< real64 > const dRdU )
153  {
154  real64 dHydraulicAperture_dNormalJump = 0.0;
155  real64 dHydraulicAperture_dTraction = 0.0;
156  real64 fractureTraction = 0.0;
157  real64 const hydraulicAperture = hydraulicApertureWrapper.computeHydraulicAperture( aperture,
158  fractureTraction,
159  fractureState,
160  dHydraulicAperture_dNormalJump,
161  dHydraulicAperture_dTraction );
162  GEOS_UNUSED_VAR( hydraulicAperture );
163 
164  constexpr integer kfSign[2] = { -1, 1 };
165  for( localIndex kf = 0; kf < 2; ++kf )
166  {
167  for( localIndex a = 0; a < numNodesPerFace; ++a )
168  {
169  for( int i = 0; i < 3; ++i )
170  {
171  nodeDOF[kf * 3 * numNodesPerFace + 3 * a + i] = dispDofNumber[faceToNodeMap( elemsToFaces[kf], a )] + i;
172 
173  real64 const dNormalJump_dDisplacement = kfSign[kf] * Nbar[i] / numNodesPerFace;
174  real64 const dHydraulicAperture_dDisplacement = dHydraulicAperture_dNormalJump * dNormalJump_dDisplacement;
175  real64 const dVolume_dDisplacement = area * dHydraulicAperture_dDisplacement;
176 
177  dRdU( kf * 3 * numNodesPerFace + 3 * a + i ) = dens * dVolume_dDisplacement;
178  }
179  }
180  }
181  }
182 
183  GEOS_HOST_DEVICE
184  inline
185  static void
186  computeFluxDerivative( localIndex const kfe2,
187  localIndex const numNodesPerFace,
188  arraySlice1d< localIndex const > const & columns,
189  arraySlice1d< real64 const > const & values,
190  arrayView2d< localIndex const > const elemsToFaces,
191  ArrayOfArraysView< localIndex const > const faceToNodeMap,
192  arrayView1d< globalIndex const > const dispDofNumber,
193  real64 const (&Nbar)[ 3 ],
194  globalIndex (& nodeDOF)[8 * 3],
195  arraySlice1d< real64 > const dRdU )
196  {
197  constexpr integer kfSign[2] = { -1, 1 };
198 
199  real64 const dR_dNormalJump = values[kfe2];
200  localIndex const ei2 = columns[kfe2];
201 
202  for( localIndex kf = 0; kf < 2; ++kf )
203  {
204  for( localIndex a = 0; a < numNodesPerFace; ++a )
205  {
206  for( localIndex i = 0; i < 3; ++i )
207  {
208  nodeDOF[kf * 3 * numNodesPerFace + 3 * a + i] = dispDofNumber[faceToNodeMap( elemsToFaces[ei2][kf], a )] + i;
209 
210  real64 const dNormalJump_dDisplacement = kfSign[kf] * Nbar[i] / numNodesPerFace;
211 
212  dRdU( kf * 3 * numNodesPerFace + 3 * a + i ) = dR_dNormalJump * dNormalJump_dDisplacement;
213  }
214  }
215  }
216  }
217 
218  template< typename POLICY, typename HYDRAULICAPERTURE_WRAPPER >
219  static void
220  launch( localIndex const size,
221  globalIndex const rankOffset,
222  HYDRAULICAPERTURE_WRAPPER const & hydraulicApertureWrapper,
223  integer const useQuasiNewton,
224  arrayView2d< localIndex const > const elemsToFaces,
225  ArrayOfArraysView< localIndex const > const faceToNodeMap,
226  arrayView2d< real64 const > const faceNormal,
227  arrayView1d< real64 const > const area,
228  arrayView1d< real64 const > const aperture,
229  arrayView1d< integer const > const fractureState,
230  arrayView1d< globalIndex const > const presDofNumber,
231  arrayView1d< globalIndex const > const dispDofNumber,
232  arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const dens,
233  CRSMatrixView< real64 const, localIndex const > const dFluxResidual_dNormalJump,
234  CRSMatrixView< real64, globalIndex const > const & localMatrix )
235  {
236  forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex ei )
237  {
238  localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elemsToFaces[ei][0] );
239 
240  real64 Nbar[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( faceNormal[elemsToFaces[ei][0]] );
241  LvArray::tensorOps::subtract< 3 >( Nbar, faceNormal[elemsToFaces[ei][1]] );
242  LvArray::tensorOps::normalize< 3 >( Nbar );
243 
244  globalIndex const rowNumber = presDofNumber[ei] - rankOffset;
245  globalIndex nodeDOF[8 * 3];
246  stackArray1d< real64, 24 > dRdU( 2 * numNodesPerFace * 3 );
247 //
248  computeAccumulationDerivative( hydraulicApertureWrapper,
249  numNodesPerFace,
250  elemsToFaces[ei],
251  faceToNodeMap,
252  dispDofNumber,
253  Nbar,
254  area[ei],
255  aperture[ei],
256  dens[ei][0],
257  fractureState[ei],
258  nodeDOF,
259  dRdU );
260 
261  if( rowNumber >= 0 && rowNumber < localMatrix.numRows() )
262  {
263  localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( rowNumber,
264  nodeDOF,
265  dRdU.data(),
266  2 * numNodesPerFace * 3 );
267  }
268 //
269  if( useQuasiNewton == 0 ) // when Quasi Newton is not enabled - add flux derivatives
270  {
271  localIndex const numColumns = dFluxResidual_dNormalJump.numNonZeros( ei );
272  arraySlice1d< localIndex const > const & columns = dFluxResidual_dNormalJump.getColumns( ei );
273  arraySlice1d< real64 const > const & values = dFluxResidual_dNormalJump.getEntries( ei );
274 
275  for( localIndex kfe2 = 0; kfe2 < numColumns; ++kfe2 )
276  {
277  computeFluxDerivative( kfe2,
278  numNodesPerFace,
279  columns,
280  values,
281  elemsToFaces,
282  faceToNodeMap,
283  dispDofNumber,
284  Nbar,
285  nodeDOF,
286  dRdU );
287 
288  if( rowNumber >= 0 && rowNumber < localMatrix.numRows() )
289  {
290  localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( rowNumber,
291  nodeDOF,
292  dRdU.data(),
293  2 * numNodesPerFace * 3 );
294  }
295  }
296  }
297  } );
298  }
299 };
300 
301 } /* namespace hydrofractureSolverKernels */
302 
303 } /* namespace geos */
304 
305 #endif // GEOS_PHYSICSSOLVERS_MULTIPHYSICS_HYDROFRACTURESOLVERKERNELS_HPP_
GEOS_HOST_DEVICE
#define GEOS_HOST_DEVICE
Marks a host-device function.
Definition: GeosxMacros.hpp:49
GEOS_UNUSED_VAR
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
Definition: GeosxMacros.hpp:84
geos::arrayView1d
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:188
geos::CRSMatrixView
LvArray::CRSMatrixView< T, COL_INDEX, localIndex const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
Definition: DataTypes.hpp:318
geos::globalIndex
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
Definition: DataTypes.hpp:87
geos::stackArray1d
StackArray< T, 1, MAXSIZE > stackArray1d
Alias for 1D stack array.
Definition: DataTypes.hpp:196
geos::ArrayOfArraysView
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:294
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::arraySlice1d
ArraySlice< T, 1, USD > arraySlice1d
Alias for 1D array slice.
Definition: DataTypes.hpp:192
geos::integer
std::int32_t integer
Signed integer type.
Definition: DataTypes.hpp:81
geos::arrayView2d
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:204
geos::fields::contact::FractureState::Open
@ Open
element is open: no constraints are imposed.
Definition: FractureState.hpp:46
geos::fields::contact::FractureState::Stick
@ Stick
element is closed: no jump across the discontinuity.
Definition: FractureState.hpp:43
geos::hydrofractureSolverKernels::DeformationUpdateKernel::launch
static std::tuple< double, double, double, double, double, double > launch(localIndex const size, HYDRAULICAPERTURE_WRAPPER const &hydraulicApertureWrapper, POROUS_WRAPPER const &porousMaterialWrapper, arrayView2d< real64 const, nodes::TOTAL_DISPLACEMENT_USD > const &u, arrayView2d< real64 const > const &faceNormal, ArrayOfArraysView< localIndex const > const &faceToNodeMap, arrayView2d< localIndex const > const &elemsToFaces, arrayView1d< real64 const > const &area, arrayView1d< real64 const > const &volume, arrayView1d< real64 > const &deltaVolume, arrayView1d< real64 > const &aperture, arrayView1d< integer > const &fractureState, arrayView1d< real64 > const &hydraulicAperture)
Definition: HydrofractureSolverKernels.hpp:36
Generated by doxygen 1.9.1