GEOS
HydrostaticPressureKernel.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 
20 #ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_HYDROSTATICPRESSUREKERNEL_HPP
21 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_HYDROSTATICPRESSUREKERNEL_HPP
22 
23 #include "common/DataTypes.hpp"
24 #include "common/GEOS_RAJA_Interface.hpp"
25 #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp"
26 
27 namespace geos
28 {
29 
30 namespace singlePhaseBaseKernels
31 {
32 
33 /******************************** HydrostaticPressureKernel ********************************/
34 
36 {
37 
38  template< typename FLUID_WRAPPER >
39  static bool
40  computeHydrostaticPressure( integer const maxNumEquilIterations,
41  real64 const & equilTolerance,
42  real64 const (&gravVector)[ 3 ],
43  FLUID_WRAPPER fluidWrapper,
44  real64 const & refElevation,
45  real64 const & refPres,
46  real64 const & refDens,
47  real64 const & newElevation,
48  real64 & newPres,
49  real64 & newDens )
50  {
51  // Step 1: guess the pressure with the refDensity
52 
53  real64 const gravCoef = gravVector[2] * ( refElevation - newElevation );
54  real64 pres0 = refPres - refDens * gravCoef;
55  real64 pres1 = 0.0;
56 
57  // Step 2: compute the mass density at this elevation using the guess, and update pressure
58 
59  real64 dens = 0.0;
60  real64 visc = 0.0;
61  constitutive::SingleFluidBaseUpdate::computeValues( fluidWrapper,
62  pres0,
63  dens,
64  visc );
65  pres1 = refPres - 0.5 * ( refDens + dens ) * gravCoef;
66 
67  // Step 3: fixed-point iteration until convergence
68 
69  bool equilHasConverged = false;
70  for( localIndex eqIter = 0; eqIter < maxNumEquilIterations; ++eqIter )
71  {
72 
73  // check convergence
74  equilHasConverged = ( LvArray::math::abs( pres0 - pres1 ) < equilTolerance );
75  pres0 = pres1;
76 
77  // if converged, move on
78  if( equilHasConverged )
79  {
80  break;
81  }
82 
83  // compute the density at this elevation using the previous pressure, and compute the new pressure
84  constitutive::SingleFluidBaseUpdate::computeValues( fluidWrapper,
85  pres0,
86  dens,
87  visc );
88  pres1 = refPres - 0.5 * ( refDens + dens ) * gravCoef;
89  }
90 
91  // Step 4: save the hydrostatic pressure and the corresponding density
92 
93  newPres = pres1;
94  newDens = dens;
95 
96  return equilHasConverged;
97  }
98 
99 
100  template< typename FLUID_WRAPPER >
101  static bool
102  launch( localIndex const size,
103  integer const maxNumEquilIterations,
104  real64 const equilTolerance,
105  real64 const (&gravVector)[ 3 ],
106  real64 const & minElevation,
107  real64 const & elevationIncrement,
108  real64 const & datumElevation,
109  real64 const & datumPres,
110  FLUID_WRAPPER fluidWrapper,
111  arrayView1d< arrayView1d< real64 > const > elevationValues,
112  arrayView1d< real64 > pressureValues )
113  {
114  bool hasConverged = true;
115 
116  // Step 1: compute the mass density at the datum elevation
117 
118  real64 datumDens = 0.0;
119  real64 datumVisc = 0.0;
120 
121  constitutive::SingleFluidBaseUpdate::computeValues( fluidWrapper,
122  datumPres,
123  datumDens,
124  datumVisc );
125 
126  // Step 2: find the closest elevation to datumElevation
127 
128  forAll< parallelHostPolicy >( size, [=] ( localIndex const i )
129  {
130  real64 const elevation = minElevation + i * elevationIncrement;
131  elevationValues[0][i] = elevation;
132  } );
133  integer const iRef = LvArray::sortedArrayManipulation::find( elevationValues[0].begin(),
134  elevationValues[0].size(),
135  datumElevation );
136 
137 
138  // Step 3: compute the mass density and pressure at the reference elevation
139 
140  array1d< real64 > dens( pressureValues.size() );
141 
142  bool const hasConvergedRef =
143  computeHydrostaticPressure( maxNumEquilIterations,
144  equilTolerance,
145  gravVector,
146  fluidWrapper,
147  datumElevation,
148  datumPres,
149  datumDens,
150  elevationValues[0][iRef],
151  pressureValues[iRef],
152  dens[iRef] );
153  if( !hasConvergedRef )
154  {
155  return false;
156  }
157 
158 
159  // Step 4: for each elevation above the reference elevation, compute the pressure
160 
161  localIndex const numEntriesAboveRef = size - iRef - 1;
162  forAll< serialPolicy >( numEntriesAboveRef, [=, &hasConverged] ( localIndex const i )
163  {
164  bool const hasConvergedAboveRef =
165  computeHydrostaticPressure( maxNumEquilIterations,
166  equilTolerance,
167  gravVector,
168  fluidWrapper,
169  elevationValues[0][iRef+i],
170  pressureValues[iRef+i],
171  dens[iRef+i],
172  elevationValues[0][iRef+i+1],
173  pressureValues[iRef+i+1],
174  dens[iRef+i+1] );
175  if( !hasConvergedAboveRef )
176  {
177  hasConverged = false;
178  }
179 
180 
181  } );
182 
183  // Step 5: for each elevation below the reference elevation, compute the pressure
184 
185  localIndex const numEntriesBelowRef = iRef;
186  forAll< serialPolicy >( numEntriesBelowRef, [=, &hasConverged] ( localIndex const i )
187  {
188  bool const hasConvergedBelowRef =
189  computeHydrostaticPressure( maxNumEquilIterations,
190  equilTolerance,
191  gravVector,
192  fluidWrapper,
193  elevationValues[0][iRef-i],
194  pressureValues[iRef-i],
195  dens[iRef-i],
196  elevationValues[0][iRef-i-1],
197  pressureValues[iRef-i-1],
198  dens[iRef-i-1] );
199  if( !hasConvergedBelowRef )
200  {
201  hasConverged = false;
202  }
203  } );
204 
205  return hasConverged;
206  }
207 };
208 
209 } // namespace singlePhaseBaseKernels
210 
211 } // namespace geos
212 
213 #endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_HYDROSTATICPRESSUREKERNEL_HPP
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:180
double real64
64-bit floating point type.
Definition: DataTypes.hpp:99
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:85
std::int32_t integer
Signed integer type.
Definition: DataTypes.hpp:82
Array< T, 1 > array1d
Alias for 1D array.
Definition: DataTypes.hpp:176