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 Total, S.A
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  ArrayOfArraysView< 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  if( elemsToFaces.sizeOfArray( kfe ) != 2 )
70  { return; }
71 
72  localIndex const kf0 = elemsToFaces[kfe][0];
73  localIndex const kf1 = elemsToFaces[kfe][1];
74  localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 );
75  real64 temp[ 3 ] = { 0 };
76  for( localIndex a=0; a<numNodesPerFace; ++a )
77  {
78  LvArray::tensorOps::add< 3 >( temp, u[ faceToNodeMap( kf0, a ) ] );
79  LvArray::tensorOps::subtract< 3 >( temp, u[ faceToNodeMap( kf1, a ) ] );
80  }
81 
82  // TODO this needs a proper contact based strategy for aperture
83  real64 const normalJump = -LvArray::tensorOps::AiBi< 3 >( temp, faceNormal[kf0] ) / numNodesPerFace;
84  maxApertureChange.max( std::fabs( normalJump - aperture[kfe] ));
85  aperture[kfe] = normalJump;
86  minAperture.min( aperture[kfe] );
87  maxAperture.max( aperture[kfe] );
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  dHydraulicAperture_dNormalJump,
95  dHydraulicAperture_dNormalTraction );
96 
97  maxHydraulicApertureChange.max( std::fabs( newHydraulicAperture - hydraulicAperture[kfe] ));
98  real64 const oldHydraulicAperture = hydraulicAperture[kfe];
99  hydraulicAperture[kfe] = newHydraulicAperture;
100  minHydraulicAperture.min( hydraulicAperture[kfe] );
101  maxHydraulicAperture.max( hydraulicAperture[kfe] );
102 
103  real64 const jump[3] = { normalJump, 0.0, 0.0 };
104  real64 const traction[3] = {0.0, 0.0, 0.0};
105 
106  porousMaterialWrapper.updateStateFromPressureApertureJumpAndTraction( kfe, 0, 0.0,
107  oldHydraulicAperture, newHydraulicAperture,
108  dHydraulicAperture_dNormalJump,
109  jump, traction );
110 
111 #ifdef GEOS_USE_SEPARATION_COEFFICIENT
112  real64 const s = aperture[kfe] / apertureAtFailure[kfe];
113  if( separationCoeff0[kfe]<1.0 && s>separationCoeff0[kfe] )
114  {
115  if( s >= 1.0 )
116  {
117  separationCoeff[kfe] = 1.0;
118  dSeparationCoeff_dAper[kfe] = 0.0;
119  }
120  else
121  {
122  separationCoeff[kfe] = s;
123  dSeparationCoeff_dAper[kfe] = 1.0/apertureAtFailure[kfe];
124  }
125  }
126 #endif
127  deltaVolume[kfe] = hydraulicAperture[kfe] * area[kfe] - volume[kfe];
128  } );
129 
130  return std::make_tuple( maxApertureChange.get(), maxHydraulicApertureChange.get(), minAperture.get(), maxAperture.get(), minHydraulicAperture.get(), maxHydraulicAperture.get() );
131  }
132 };
133 
134 struct FluidMassResidualDerivativeAssemblyKernel
135 {
136  template< typename HYDRAULICAPERTURE_WRAPPER >
137  GEOS_HOST_DEVICE
138  inline
139  static void
140  computeAccumulationDerivative( HYDRAULICAPERTURE_WRAPPER const & hydraulicApertureWrapper,
141  localIndex const numNodesPerFace,
142  arraySlice1d< localIndex const > const elemsToFaces,
143  ArrayOfArraysView< localIndex const > const faceToNodeMap,
144  arrayView1d< globalIndex const > const dispDofNumber,
145  real64 const (&Nbar)[ 3 ],
146  real64 const & area,
147  real64 const & aperture,
148  real64 const & dens,
149  globalIndex (& nodeDOF)[8 * 3],
150  arraySlice1d< real64 > const dRdU )
151  {
152  real64 dHydraulicAperture_dNormalJump = 0.0;
153  real64 dHydraulicAperture_dTraction = 0.0;
154  real64 fractureTraction = 0.0;
155  real64 const hydraulicAperture = hydraulicApertureWrapper.computeHydraulicAperture( aperture,
156  fractureTraction,
157  dHydraulicAperture_dNormalJump,
158  dHydraulicAperture_dTraction );
159  GEOS_UNUSED_VAR( hydraulicAperture );
160 
161  constexpr integer kfSign[2] = { -1, 1 };
162  for( localIndex kf = 0; kf < 2; ++kf )
163  {
164  for( localIndex a = 0; a < numNodesPerFace; ++a )
165  {
166  for( int i = 0; i < 3; ++i )
167  {
168  nodeDOF[kf * 3 * numNodesPerFace + 3 * a + i] = dispDofNumber[faceToNodeMap( elemsToFaces[kf], a )] + i;
169 
170  real64 const dNormalJump_dDisplacement = kfSign[kf] * Nbar[i] / numNodesPerFace;
171  real64 const dHydraulicAperture_dDisplacement = dHydraulicAperture_dNormalJump * dNormalJump_dDisplacement;
172  real64 const dVolume_dDisplacement = area * dHydraulicAperture_dDisplacement;
173 
174  dRdU( kf * 3 * numNodesPerFace + 3 * a + i ) = dens * dVolume_dDisplacement;
175  }
176  }
177  }
178  }
179 
180  GEOS_HOST_DEVICE
181  inline
182  static void
183  computeFluxDerivative( localIndex const kfe2,
184  localIndex const numNodesPerFace,
185  arraySlice1d< localIndex const > const & columns,
186  arraySlice1d< real64 const > const & values,
187  ArrayOfArraysView< localIndex const > const elemsToFaces,
188  ArrayOfArraysView< localIndex const > const faceToNodeMap,
189  arrayView1d< globalIndex const > const dispDofNumber,
190  real64 const (&Nbar)[ 3 ],
191  globalIndex (& nodeDOF)[8 * 3],
192  arraySlice1d< real64 > const dRdU )
193  {
194  constexpr integer kfSign[2] = { -1, 1 };
195 
196  real64 const dR_dNormalJump = values[kfe2];
197  localIndex const ei2 = columns[kfe2];
198 
199  for( localIndex kf = 0; kf < 2; ++kf )
200  {
201  for( localIndex a = 0; a < numNodesPerFace; ++a )
202  {
203  for( localIndex i = 0; i < 3; ++i )
204  {
205  nodeDOF[kf * 3 * numNodesPerFace + 3 * a + i] = dispDofNumber[faceToNodeMap( elemsToFaces[ei2][kf], a )] + i;
206 
207  real64 const dNormalJump_dDisplacement = kfSign[kf] * Nbar[i] / numNodesPerFace;
208 
209  dRdU( kf * 3 * numNodesPerFace + 3 * a + i ) = dR_dNormalJump * dNormalJump_dDisplacement;
210  }
211  }
212  }
213  }
214 
215  template< typename POLICY, typename HYDRAULICAPERTURE_WRAPPER >
216  static void
217  launch( localIndex const size,
218  globalIndex const rankOffset,
219  HYDRAULICAPERTURE_WRAPPER const & hydraulicApertureWrapper,
220  integer const useQuasiNewton,
221  ArrayOfArraysView< localIndex const > const elemsToFaces,
222  ArrayOfArraysView< localIndex const > const faceToNodeMap,
223  arrayView2d< real64 const > const faceNormal,
224  arrayView1d< real64 const > const area,
225  arrayView1d< real64 const > const aperture,
226  arrayView1d< globalIndex const > const presDofNumber,
227  arrayView1d< globalIndex const > const dispDofNumber,
228  arrayView2d< real64 const > const dens,
229  CRSMatrixView< real64 const, localIndex const > const dFluxResidual_dNormalJump,
230  CRSMatrixView< real64, globalIndex const > const & localMatrix )
231  {
232  forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex ei )
233  {
234  localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elemsToFaces[ei][0] );
235 
236  real64 Nbar[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( faceNormal[elemsToFaces[ei][0]] );
237  LvArray::tensorOps::subtract< 3 >( Nbar, faceNormal[elemsToFaces[ei][1]] );
238  LvArray::tensorOps::normalize< 3 >( Nbar );
239 
240  globalIndex const rowNumber = presDofNumber[ei] - rankOffset;
241  globalIndex nodeDOF[8 * 3];
242  stackArray1d< real64, 24 > dRdU( 2 * numNodesPerFace * 3 );
243 //
244  computeAccumulationDerivative( hydraulicApertureWrapper,
245  numNodesPerFace,
246  elemsToFaces[ei],
247  faceToNodeMap,
248  dispDofNumber,
249  Nbar,
250  area[ei],
251  aperture[ei],
252  dens[ei][0],
253  nodeDOF,
254  dRdU );
255 
256  if( rowNumber >= 0 && rowNumber < localMatrix.numRows() )
257  {
258  localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( rowNumber,
259  nodeDOF,
260  dRdU.data(),
261  2 * numNodesPerFace * 3 );
262  }
263 //
264  if( useQuasiNewton == 0 ) // when Quasi Newton is not enabled - add flux derivatives
265  {
266  localIndex const numColumns = dFluxResidual_dNormalJump.numNonZeros( ei );
267  arraySlice1d< localIndex const > const & columns = dFluxResidual_dNormalJump.getColumns( ei );
268  arraySlice1d< real64 const > const & values = dFluxResidual_dNormalJump.getEntries( ei );
269 
270  for( localIndex kfe2 = 0; kfe2 < numColumns; ++kfe2 )
271  {
272  computeFluxDerivative( kfe2,
273  numNodesPerFace,
274  columns,
275  values,
276  elemsToFaces,
277  faceToNodeMap,
278  dispDofNumber,
279  Nbar,
280  nodeDOF,
281  dRdU );
282 
283  if( rowNumber >= 0 && rowNumber < localMatrix.numRows() )
284  {
285  localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( rowNumber,
286  nodeDOF,
287  dRdU.data(),
288  2 * numNodesPerFace * 3 );
289  }
290  }
291  }
292  } );
293  }
294 };
295 
296 } /* namespace hydrofractureSolverKernels */
297 
298 } /* namespace geos */
299 
300 #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, ArrayOfArraysView< 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