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