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_COMPOSITIONAL_HYDROSTATICPRESSUREKERNEL_HPP
21 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_HYDROSTATICPRESSUREKERNEL_HPP
22 
23 #include "common/DataLayouts.hpp"
24 #include "common/DataTypes.hpp"
25 #include "common/GEOS_RAJA_Interface.hpp"
26 #include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
27 #include "functions/TableFunction.hpp"
28 
29 namespace geos
30 {
31 
32 namespace isothermalCompositionalMultiphaseBaseKernels
33 {
34 
35 /******************************** HydrostaticPressureKernel ********************************/
36 
37 struct HydrostaticPressureKernel
38 {
39  // TODO: this type of constants should be centralized somewhere or provided by fluid model
40  static real64 constexpr MIN_FOR_PHASE_PRESENCE = 1e-12;
41  static real64 constexpr MIN_FOR_COMP_PRESENCE = 1e-12;
42 
43  enum class ReturnType : integer
44  {
45  FAILED_TO_CONVERGE = 0,
46  DETECTED_MULTIPHASE_FLOW = 1,
47  SUCCESS = 2,
48  DETECTED_SINGLEPHASE_FLOW = 3,
49  PHASE_CORRECTION_NOT_NEEDED = 4
50  };
51 
52  template< typename FLUID_WRAPPER >
53  static ReturnType
54  computeHydrostaticPressure( integer const & numComps,
55  integer const & numPhases,
56  integer const & ipGas,
57  integer const & ipOil,
58  integer const & ipWater,
59  integer const & ipInit,
60  arrayView1d< real64 const > const & phaseContacts,
61  arrayView1d< real64 const > const & phaseMinVolumeFraction,
62  integer const maxNumEquilIterations,
63  real64 const & equilTolerance,
64  real64 const (&gravVector)[ 3 ],
65  FLUID_WRAPPER fluidWrapper,
66  arrayView1d< TableFunction::KernelWrapper const > compFracTableWrappers,
67  TableFunction::KernelWrapper tempTableWrapper,
68  real64 const & refElevation,
69  arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > const & refPres,
70  arraySlice1d< real64 const > const & refPhaseMassDens,
71  real64 const & newElevation,
72  arraySlice1d< real64, constitutive::multifluid::USD_PHASE - 2 > const & newPres,
73  arraySlice1d< real64 > const & newPhaseMassDens,
74  arraySlice1d< real64, constitutive::multifluid::USD_PHASE - 2 > const & newPhaseDens,
75  arraySlice2d< real64, constitutive::multifluid::USD_PHASE_COMP - 2 > const & newPhaseCompFrac )
76  {
77  // fluid properties at this elevation
78  StackArray< real64, 2, constitutive::MultiFluidBase::MAX_NUM_COMPONENTS, compflow::LAYOUT_COMP > inputComposition( 1, numComps );
79  StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > phaseFrac( 1, 1, numPhases );
80  StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > phaseDens( 1, 1, numPhases );
81  StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > phaseMassDens( 1, 1, numPhases );
82  StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > phaseVisc( 1, 1, numPhases );
83  StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > phaseEnthalpy( 1, 1, numPhases );
84  StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > phaseInternalEnergy( 1, 1, numPhases );
85  StackArray< real64, 4, constitutive::MultiFluidBase::MAX_NUM_PHASES *constitutive::MultiFluidBase::MAX_NUM_COMPONENTS,
86  constitutive::multifluid::LAYOUT_PHASE_COMP > phaseCompFrac( 1, 1, numPhases, numComps );
87  real64 totalDens = 0.0;
88 
89  bool const isSinglePhaseInitialisation = phaseContacts.size() == 0;
90 
91  bool isSinglePhaseFlow = true;
92 
93  // Step 1: compute the hydrostatic pressure at the current elevation
94 
95  real64 const gravCoef = gravVector[2] * ( refElevation - newElevation );
96  real64 const temperature = tempTableWrapper.compute( &newElevation );
97  for( integer ic = 0; ic < numComps; ++ic )
98  {
99  inputComposition[0][ic] = compFracTableWrappers[ic].compute( &newElevation );
100  }
101 
102  // Step 2: guess the pressure with the refPhaseMassDensity
103  StackArray< real64, 2, 2*constitutive::MultiFluidBase::MAX_NUM_PHASES > phasePressures( 2, numPhases );
104  StackArray< real64, 2, constitutive::MultiFluidBase::MAX_NUM_COMPONENTS, compflow::LAYOUT_COMP > compFrac( 1, numComps );
105  for( localIndex ip = 0; ip < numPhases; ++ip )
106  {
107  phasePressures[0][ip] = refPres[ip] - refPhaseMassDens[ip] * gravCoef;
108  }
109 
110  // Step 3: determine which phase pressure to use as the flash pressure
111  integer const flashPhase = isSinglePhaseInitialisation ? ipInit :
112  evaluateFlashPhaseIndex( numPhases,
113  ipGas,
114  ipOil,
115  ipWater,
116  ipInit,
117  newElevation,
118  phaseContacts );
119 
120  // Step 4: compute the mass density at this elevation using the guess, and update pressure
121  constitutive::MultiFluidBase::KernelWrapper::computeValues( fluidWrapper,
122  phasePressures[0][flashPhase], // flash pressure
123  temperature,
124  inputComposition[0],
125  phaseFrac[0][0],
126  phaseDens[0][0],
127  phaseMassDens[0][0],
128  phaseVisc[0][0],
129  phaseEnthalpy[0][0],
130  phaseInternalEnergy[0][0],
131  phaseCompFrac[0][0],
132  totalDens );
133 
134  // Step 5: Ensure the correct phases exist. If not, apply phase correction.
135  if( isSinglePhaseInitialisation )
136  {
137  for( integer ic = 0; ic < numComps; ++ic )
138  {
139  compFrac[0][ic] = inputComposition[0][ic];
140  }
141  }
142  else
143  {
144  phaseCorrection( numComps,
145  numPhases,
146  ipGas,
147  ipWater,
148  phaseMinVolumeFraction,
149  phasePressures[0][flashPhase],
150  temperature,
151  inputComposition[0],
152  phaseFrac[0][0],
153  compFrac[0],
154  fluidWrapper );
155  }
156  for( localIndex ip = 0; ip < numPhases; ++ip )
157  {
158  phasePressures[1][ip] = refPres[ip] - 0.5 * ( refPhaseMassDens[ip] + phaseMassDens[0][0][ip] ) * gravCoef;
159  }
160 
161  // Step 6: fixed-point iteration until convergence
162  bool equilHasConverged = false;
163  for( integer eqIter = 0; eqIter < maxNumEquilIterations; ++eqIter )
164  {
165  // check convergence
166  real64 pressureError = 0.0;
167  for( integer ip = 0; ip < numPhases; ++ip )
168  {
169  real64 const dp = LvArray::math::abs( phasePressures[0][ip] - phasePressures[1][ip] );
170  pressureError = LvArray::math::max( pressureError, dp );
171  phasePressures[0][ip] = phasePressures[1][ip];
172  }
173  equilHasConverged = ( pressureError < equilTolerance );
174 
175  // if converged, check number of phases and move on
176  if( equilHasConverged )
177  {
178  break;
179  }
180 
181  // compute the mass density at this elevation using the previous pressure, and compute the new pressure
182  constitutive::MultiFluidBase::KernelWrapper::computeValues( fluidWrapper,
183  phasePressures[0][flashPhase], // flash pressure
184  temperature,
185  compFrac[0],
186  phaseFrac[0][0],
187  phaseDens[0][0],
188  phaseMassDens[0][0],
189  phaseVisc[0][0],
190  phaseEnthalpy[0][0],
191  phaseInternalEnergy[0][0],
192  phaseCompFrac[0][0],
193  totalDens );
194  for( integer ip = 0; ip < numPhases; ++ip )
195  {
196  phasePressures[1][ip] = refPres[ip] - 0.5 * ( refPhaseMassDens[ip] + phaseMassDens[0][0][ip] ) * gravCoef;
197  }
198  }
199 
200  // Step 7: save the hydrostatic pressure and the corresponding density and composition
201  // Also count the number of phases on exit
202  integer numberOfPhases = 0;
203  for( integer ip = 0; ip < numPhases; ++ip )
204  {
205  newPres[ip] = phasePressures[1][ip];
206  newPhaseMassDens[ip] = phaseMassDens[0][0][ip];
207  newPhaseDens[ip] = phaseDens[0][0][ip];
208  for( integer ic = 0; ic < numComps; ++ic )
209  {
210  newPhaseCompFrac[ip][ic] = phaseCompFrac[0][0][ip][ic];
211  }
212  if( phaseFrac[0][0][ip] > MIN_FOR_PHASE_PRESENCE )
213  {
214  numberOfPhases++;
215  }
216  }
217  if( 1 < numberOfPhases )
218  {
219  isSinglePhaseFlow = false;
220  }
221 
222  if( !equilHasConverged )
223  {
224  return ReturnType::FAILED_TO_CONVERGE;
225  }
226  else if( !isSinglePhaseFlow )
227  {
228  return ReturnType::DETECTED_MULTIPHASE_FLOW;
229  }
230  else
231  {
232  return ReturnType::SUCCESS;
233  }
234  }
235 
236  template< typename FLUID_WRAPPER >
237  static ReturnType
238  computeHydrostaticPressureAtMultipleElevations( localIndex const & startElevationIndex,
239  localIndex const & endElevationIndex,
240  integer const & numComps,
241  integer const & numPhases,
242  integer const & ipGas,
243  integer const & ipOil,
244  integer const & ipWater,
245  integer const & ipInit,
246  arrayView1d< real64 const > const & phaseContacts,
247  arrayView1d< real64 const > const & phaseMinVolumeFraction,
248  integer const & maxNumEquilIterations,
249  real64 const & equilTolerance,
250  real64 const (&gravVector)[ 3 ],
251  FLUID_WRAPPER fluidWrapper,
252  arrayView1d< TableFunction::KernelWrapper const > compFracTableWrappers,
253  TableFunction::KernelWrapper tempTableWrapper,
254  arrayView1d< arrayView1d< real64 > const > elevationValues,
255  arrayView3d< real64, constitutive::multifluid::USD_PHASE > const & pressureValues,
256  arrayView2d< real64 > const & phaseMassDens,
257  arrayView3d< real64, constitutive::multifluid::USD_PHASE > const & phaseDens,
258  arrayView4d< real64, constitutive::multifluid::USD_PHASE_COMP > const & phaseCompFrac )
259  {
260  // startElevIndex is the reference point
261  localIndex const numEntries = LvArray::math::abs( startElevationIndex - endElevationIndex );
262  localIndex const step = ( endElevationIndex >= startElevationIndex ) ? 1 : -1;
263  ReturnType returnVal = ReturnType::SUCCESS;
264  forAll< serialPolicy >( numEntries, [=, &returnVal] ( localIndex const i )
265  {
266  localIndex const ref = startElevationIndex + i * step;
267  localIndex const next = ref + step;
268  ReturnType const iReturnVal =
269  computeHydrostaticPressure( numComps,
270  numPhases,
271  ipGas,
272  ipOil,
273  ipWater,
274  ipInit,
275  phaseContacts,
276  phaseMinVolumeFraction,
277  maxNumEquilIterations,
278  equilTolerance,
279  gravVector,
280  fluidWrapper,
281  compFracTableWrappers,
282  tempTableWrapper,
283  elevationValues[0][ref],
284  pressureValues[ref][0],
285  phaseMassDens[ref],
286  elevationValues[0][next],
287  pressureValues[next][0],
288  phaseMassDens[next],
289  phaseDens[next][0],
290  phaseCompFrac[next][0] );
291  if( iReturnVal == ReturnType::FAILED_TO_CONVERGE )
292  {
293  returnVal = ReturnType::FAILED_TO_CONVERGE;
294  }
295  else if( iReturnVal == ReturnType::DETECTED_MULTIPHASE_FLOW )
296  {
297  returnVal = ReturnType::DETECTED_MULTIPHASE_FLOW;
298  }
299  } );
300 
301  return returnVal;
302  }
303 
304  template< typename FLUID_WRAPPER >
305  static ReturnType
306  marchBetweenTwoElevations( real64 const & startElevation,
307  real64 const & endElevation,
308  integer const & numComps,
309  integer const & numPhases,
310  integer const & ipGas,
311  integer const & ipOil,
312  integer const & ipWater,
313  integer const & ipInit,
314  arrayView1d< real64 const > const & phaseContacts,
315  arrayView1d< real64 const > const & phaseMinVolumeFraction,
316  integer const maxNumEquilIterations,
317  real64 const & equilTolerance,
318  real64 const (&gravVector)[ 3 ],
319  FLUID_WRAPPER fluidWrapper,
320  arrayView1d< TableFunction::KernelWrapper const > compFracTableWrappers,
321  TableFunction::KernelWrapper tempTableWrapper,
322  arrayView1d< arrayView1d< real64 > const > elevationValues,
323  arrayView3d< real64, constitutive::multifluid::USD_PHASE > const & pressureValues,
324  arrayView2d< real64 > const & phaseMassDens,
325  arrayView3d< real64, constitutive::multifluid::USD_PHASE > const & phaseDens,
326  arrayView4d< real64, constitutive::multifluid::USD_PHASE_COMP > const & phaseCompFrac )
327  {
328  // Find the primary and contact phase indices
329  integer ipPP;
330  integer ipCP;
331  evaluatePrimaryAndContactPhaseIndices( numPhases,
332  ipGas,
333  ipOil,
334  ipWater,
335  startElevation,
336  endElevation,
337  phaseContacts,
338  ipPP,
339  ipCP );
340 
341  // Find index of the closest element in elevationValues to the start elevation
342  localIndex const iStart = LvArray::sortedArrayManipulation::find( elevationValues[0].begin(),
343  elevationValues[0].size(),
344  startElevation );
345  // Find index of the closest element in elevationValues to the end elevation
346  localIndex const iEnd = LvArray::sortedArrayManipulation::find( elevationValues[0].begin(),
347  elevationValues[0].size(),
348  endElevation );
349  // March from iStart to iEnd
350  ReturnType returnVal =
351  computeHydrostaticPressureAtMultipleElevations( iStart,
352  iEnd,
353  numComps,
354  numPhases,
355  ipGas,
356  ipOil,
357  ipWater,
358  ipInit,
359  phaseContacts,
360  phaseMinVolumeFraction,
361  maxNumEquilIterations,
362  equilTolerance,
363  gravVector,
364  fluidWrapper,
365  compFracTableWrappers,
366  tempTableWrapper,
367  elevationValues,
368  pressureValues,
369  phaseMassDens,
370  phaseDens,
371  phaseCompFrac );
372 
373  // Compute phase presssures and densities at end elevation using iEnd as the reference
374  array3d< real64, constitutive::multifluid::LAYOUT_PHASE > endPressure( 1, 1, numPhases );
375  array3d< real64 > endPhaseMassDens( 1, 1, numPhases );
376  array3d< real64, constitutive::multifluid::LAYOUT_PHASE > endPhaseDens( 1, 1, numPhases );
377  array4d< real64, constitutive::multifluid::LAYOUT_PHASE_COMP > endPhaseCompFrac( 1, 1, numPhases, numComps );
378  returnVal =
379  computeHydrostaticPressure( numComps,
380  numPhases,
381  ipGas,
382  ipOil,
383  ipWater,
384  ipInit,
385  phaseContacts,
386  phaseMinVolumeFraction,
387  maxNumEquilIterations,
388  equilTolerance,
389  gravVector,
390  fluidWrapper,
391  compFracTableWrappers,
392  tempTableWrapper,
393  elevationValues[0][iEnd],
394  pressureValues[iEnd][0],
395  phaseMassDens[iEnd],
396  endElevation,
397  endPressure[0][0],
398  endPhaseMassDens[0][0],
399  endPhaseDens[0][0],
400  endPhaseCompFrac[0][0] );
401  // Compute relative error defined as the relative difference between the phase pressures at end elevation
402  real64 err = LvArray::math::abs( endPressure[0][0][ipCP] - endPressure[0][0][ipPP] ) / endPressure[0][0][ipPP];
403  int constexpr maxMarchIterations = 10;
404  real64 constexpr pressureTolerance = 1.0e-5;
405 
406  // Marching Loop
407  for( int marchIter = 1; marchIter < maxMarchIterations; ++marchIter )
408  {
409  if( err < pressureTolerance )
410  {
411  break;
412  }
413  // equate the phase pressure at the end elevation
414  endPressure[0][0][ipCP] = endPressure[0][0][ipPP];
415  real64 const iStartPrimaryPressure = pressureValues[iStart][0][ipPP]; // saves the known phase pressure at iStart
416  // Estimate the unknown phase (contact) pressure and density at iStart using the updated pressures
417  returnVal =
418  computeHydrostaticPressure( numComps,
419  numPhases,
420  ipGas,
421  ipOil,
422  ipWater,
423  ipInit,
424  phaseContacts,
425  phaseMinVolumeFraction,
426  maxNumEquilIterations,
427  equilTolerance,
428  gravVector,
429  fluidWrapper,
430  compFracTableWrappers,
431  tempTableWrapper,
432  endElevation,
433  endPressure[0][0],
434  endPhaseMassDens[0][0],
435  elevationValues[0][iStart],
436  pressureValues[iStart][0],
437  phaseMassDens[iStart],
438  phaseDens[iStart][0],
439  phaseCompFrac[iStart][0] );
440  pressureValues[iStart][0][ipPP] = iStartPrimaryPressure;
441  // March from iStart to iEnd
442  returnVal =
443  computeHydrostaticPressureAtMultipleElevations( iStart,
444  iEnd,
445  numComps,
446  numPhases,
447  ipGas,
448  ipOil,
449  ipWater,
450  ipInit,
451  phaseContacts,
452  phaseMinVolumeFraction,
453  maxNumEquilIterations,
454  equilTolerance,
455  gravVector,
456  fluidWrapper,
457  compFracTableWrappers,
458  tempTableWrapper,
459  elevationValues,
460  pressureValues,
461  phaseMassDens,
462  phaseDens,
463  phaseCompFrac );
464  // Compute phase presssures and densities at the end elevation using iEnd as the reference
465  returnVal =
466  computeHydrostaticPressure( numComps,
467  numPhases,
468  ipGas,
469  ipOil,
470  ipWater,
471  ipInit,
472  phaseContacts,
473  phaseMinVolumeFraction,
474  maxNumEquilIterations,
475  equilTolerance,
476  gravVector,
477  fluidWrapper,
478  compFracTableWrappers,
479  tempTableWrapper,
480  elevationValues[0][iEnd],
481  pressureValues[iEnd][0],
482  phaseMassDens[iEnd],
483  endElevation,
484  endPressure[0][0],
485  endPhaseMassDens[0][0],
486  endPhaseDens[0][0],
487  endPhaseCompFrac[0][0] );
488  err = LvArray::math::abs( endPressure[0][0][ipCP] - endPressure[0][0][ipPP] ) / endPressure[0][0][ipPP];
489  }
490 
491  return returnVal;
492  }
493 
494  template< typename FLUID_WRAPPER >
495  static ReturnType
496  launch( localIndex const & size,
497  integer const & numComps,
498  integer const & numPhases,
499  integer const & ipGas,
500  integer const & ipOil,
501  integer const & ipWater,
502  integer const & ipInit,
503  integer const & maxNumEquilIterations,
504  arrayView1d< real64 const > const & phaseContacts,
505  arrayView1d< real64 const > const & phaseMinVolumeFraction,
506  real64 const equilTolerance,
507  real64 const (&gravVector)[ 3 ],
508  real64 const & datumElevation,
509  real64 const & datumPres,
510  FLUID_WRAPPER fluidWrapper,
511  arrayView1d< TableFunction::KernelWrapper const > compFracTableWrappers,
512  TableFunction::KernelWrapper tempTableWrapper,
513  arrayView1d< arrayView1d< real64 > const > elevationValues,
514  arrayView3d< real64, constitutive::multifluid::USD_PHASE > const & pressureValues,
515  arrayView3d< real64, constitutive::multifluid::USD_PHASE > const & phaseDens,
516  arrayView4d< real64, constitutive::multifluid::USD_PHASE_COMP > const & phaseCompFrac )
517  {
518  ReturnType returnVal = ReturnType::SUCCESS;
519 
520  // Contacts classified as "close" and "far"
521  int const numContacts = phaseContacts.size();
522  bool const isSinglePhase = (numContacts == 0);
523  // Find the index of the contact that is closest to the datum
524  integer iContactClose = 0;
525  integer iContactFar = -1;
526  if( numContacts > 1 )
527  {
528  iContactFar = 1;
529  // Assumes phaseContacts is of size 2
530  if( LvArray::math::abs( phaseContacts[1] - datumElevation ) <=
531  LvArray::math::abs( phaseContacts[0] - datumElevation ) )
532  {
533  iContactClose = 1;
534  iContactFar = 0;
535  }
536  }
537 
538  // Temporarily set all phase pressures at datum to input datum pressure.
539  array3d< real64, constitutive::multifluid::LAYOUT_PHASE > datumPresInput( 1, 1, numPhases );
540  for( localIndex ip = 0; ip < numPhases; ++ip )
541  {
542  datumPresInput[0][0][ip] = datumPres;
543  }
544 
545  // compute the phase mass densities at datum
546  real64 const datumTemp = tempTableWrapper.compute( &datumElevation );
547  array3d< real64, constitutive::multifluid::LAYOUT_PHASE > datumPhaseMassDens( 1, 1, numPhases );
548  computeDatumPhaseMassDens( numComps,
549  numPhases,
550  ipGas,
551  ipWater,
552  phaseMinVolumeFraction,
553  datumElevation,
554  datumPres,
555  datumTemp,
556  datumPhaseMassDens,
557  compFracTableWrappers,
558  isSinglePhase,
559  fluidWrapper );
560 
561 
562  // find the closest elevation to datumElevation. Its index is denoted as iDatum
563  localIndex const iDatum = LvArray::sortedArrayManipulation::find( elevationValues[0].begin(),
564  elevationValues[0].size(),
565  datumElevation );
566 
567  // compute the mass density and pressure at iDatum elevation
568  array2d< real64 > phaseMassDens( size, numPhases );
569  // temporary array without permutation to compile on GPU
570  array1d< real64 > datumPhaseMassDensTmp( numPhases );
571  for( integer ip = 0; ip < numPhases; ++ip )
572  {
573  datumPhaseMassDensTmp[ip] = datumPhaseMassDens[0][0][ip];
574  }
575 
576  ReturnType const iDatumReturnVal =
577  computeHydrostaticPressure( numComps,
578  numPhases,
579  ipGas,
580  ipOil,
581  ipWater,
582  ipInit,
583  phaseContacts,
584  phaseMinVolumeFraction,
585  maxNumEquilIterations,
586  equilTolerance,
587  gravVector,
588  fluidWrapper,
589  compFracTableWrappers,
590  tempTableWrapper,
591  datumElevation,
592  datumPresInput[0][0],
593  datumPhaseMassDensTmp,
594  elevationValues[0][iDatum],
595  pressureValues[iDatum][0],
596  phaseMassDens[iDatum],
597  phaseDens[iDatum][0],
598  phaseCompFrac[iDatum][0] );
599  if( iDatumReturnVal == ReturnType::FAILED_TO_CONVERGE )
600  {
601  return ReturnType::FAILED_TO_CONVERGE;
602  }
603  else if( iDatumReturnVal == ReturnType::DETECTED_MULTIPHASE_FLOW )
604  {
605  returnVal = ReturnType::DETECTED_MULTIPHASE_FLOW;
606  }
607 
608  if( isSinglePhase )
609  {
610  // compute hydrostatic pressure for each elevation above the reference elevation, compute the pressure
611  returnVal = computeHydrostaticPressureAtMultipleElevations( iDatum,
612  size - 1,
613  numComps,
614  numPhases,
615  ipGas,
616  ipOil,
617  ipWater,
618  ipInit,
619  phaseContacts,
620  phaseMinVolumeFraction,
621  maxNumEquilIterations,
622  equilTolerance,
623  gravVector,
624  fluidWrapper,
625  compFracTableWrappers,
626  tempTableWrapper,
627  elevationValues,
628  pressureValues,
629  phaseMassDens,
630  phaseDens,
631  phaseCompFrac );
632 
633  // compute hydrostatic pressure for each elevation above the reference elevation, compute the pressure
634  returnVal = computeHydrostaticPressureAtMultipleElevations( iDatum,
635  0,
636  numComps,
637  numPhases,
638  ipGas,
639  ipOil,
640  ipWater,
641  ipInit,
642  phaseContacts,
643  phaseMinVolumeFraction,
644  maxNumEquilIterations,
645  equilTolerance,
646  gravVector,
647  fluidWrapper,
648  compFracTableWrappers,
649  tempTableWrapper,
650  elevationValues,
651  pressureValues,
652  phaseMassDens,
653  phaseDens,
654  phaseCompFrac );
655  }
656  else
657  {
658  integer iContactCloseIndex = iDatum;
659  integer iContactTop = iDatum;
660  integer iContactBottom = iDatum;
661 
662  // March from datum to the close contact
663  if( LvArray::math::abs( datumElevation - phaseContacts[iContactClose] ) > 1e-12 )
664  {
665  // Find index of the closest element in elevationValues to the close contact, denoted as iContact
666  iContactCloseIndex = LvArray::sortedArrayManipulation::find( elevationValues[0].begin(),
667  elevationValues[0].size(),
668  phaseContacts[iContactClose] );
669 
670  // compute hydrostatic pressure for each elevation between the datum and the closest contact
671  returnVal = marchBetweenTwoElevations( datumElevation,
672  phaseContacts[iContactClose],
673  numComps,
674  numPhases,
675  ipGas,
676  ipOil,
677  ipWater,
678  ipInit,
679  phaseContacts,
680  phaseMinVolumeFraction,
681  maxNumEquilIterations,
682  equilTolerance,
683  gravVector,
684  fluidWrapper,
685  compFracTableWrappers,
686  tempTableWrapper,
687  elevationValues,
688  pressureValues,
689  phaseMassDens,
690  phaseDens,
691  phaseCompFrac );
692  }
693 
694  integer iContactFarIndex = iContactCloseIndex;
695  // March from close contact to far contact
696  if( iContactFar != -1 && LvArray::math::abs( phaseContacts[iContactClose] - phaseContacts[iContactFar] ) > 1e-12 )
697  {
698  // Find index of the closest element in elevationValues to the far contact, denoted as iContact
699  iContactFarIndex = LvArray::sortedArrayManipulation::find( elevationValues[0].begin(),
700  elevationValues[0].size(),
701  phaseContacts[iContactFar] );
702 
703  // compute hydrostatic pressure for each elevation between the closest and the farthest contacts
704  returnVal = marchBetweenTwoElevations( phaseContacts[iContactClose],
705  phaseContacts[iContactFar],
706  numComps,
707  numPhases,
708  ipGas,
709  ipOil,
710  ipWater,
711  ipInit,
712  phaseContacts,
713  phaseMinVolumeFraction,
714  maxNumEquilIterations,
715  equilTolerance,
716  gravVector,
717  fluidWrapper,
718  compFracTableWrappers,
719  tempTableWrapper,
720  elevationValues,
721  pressureValues,
722  phaseMassDens,
723  phaseDens,
724  phaseCompFrac );
725  }
726 
727  if( iContactFar == -1 )
728  {
729  iContactTop = iContactCloseIndex;
730  iContactBottom = iContactCloseIndex;
731  }
732  else if( phaseContacts[iContactClose] > phaseContacts[iContactFar] )
733  {
734  iContactTop = iContactCloseIndex;
735  iContactBottom = iContactFarIndex;
736  }
737  else
738  {
739  iContactTop = iContactFarIndex;
740  iContactBottom = iContactCloseIndex;
741  }
742 
743  // compute hydrostatic pressure for each elevation between the top contact and the topmost elevation
744  returnVal = computeHydrostaticPressureAtMultipleElevations( iContactTop,
745  size - 1,
746  numComps,
747  numPhases,
748  ipGas,
749  ipOil,
750  ipWater,
751  ipInit,
752  phaseContacts,
753  phaseMinVolumeFraction,
754  maxNumEquilIterations,
755  equilTolerance,
756  gravVector,
757  fluidWrapper,
758  compFracTableWrappers,
759  tempTableWrapper,
760  elevationValues,
761  pressureValues,
762  phaseMassDens,
763  phaseDens,
764  phaseCompFrac );
765 
766  // compute hydrostatic pressure for each elevation between the bottom contact and the bottom-most elevation
767  returnVal = computeHydrostaticPressureAtMultipleElevations( iContactBottom,
768  0,
769  numComps,
770  numPhases,
771  ipGas,
772  ipOil,
773  ipWater,
774  ipInit,
775  phaseContacts,
776  phaseMinVolumeFraction,
777  maxNumEquilIterations,
778  equilTolerance,
779  gravVector,
780  fluidWrapper,
781  compFracTableWrappers,
782  tempTableWrapper,
783  elevationValues,
784  pressureValues,
785  phaseMassDens,
786  phaseDens,
787  phaseCompFrac );
788  }
789  return returnVal;
790  }
791 
792  template< typename FLUID_WRAPPER >
793  static ReturnType
794  phaseCorrection( integer const & numComps,
795  integer const & numPhases,
796  integer const & ipGas,
797  integer const & ipWater,
798  arrayView1d< real64 const > const & phaseMinVolumeFraction,
799  real64 const & pres,
800  real64 const & temp,
801  arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & inputComposition,
802  arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > const & phaseFrac,
803  arraySlice1d< real64, compflow::USD_COMP - 1 > const & compFrac,
804  FLUID_WRAPPER fluidWrapper )
805  {
806  // TODO: should the if conditions be independent or else ifs after the first if
807  StackArray< real64, 2, constitutive::MultiFluidBase::MAX_NUM_COMPONENTS, compflow::LAYOUT_COMP > addedCompFrac( 1, numComps );
808  integer ip_phase = -1;
809  integer ip_otherPhase = -1;
810  bool phaseCorrectionNeeded = false;
811  if( ipGas >= 0 && LvArray::math::abs( phaseFrac[ipGas] ) < MIN_FOR_PHASE_PRESENCE
812  && LvArray::math::abs( phaseMinVolumeFraction[ipGas] ) > MIN_FOR_PHASE_PRESENCE )
813  {
814  addedCompFrac[0][0] = 1.0; // hard-coded (assumes co2 is at index 0)
815  ip_phase = ipGas;
816  ip_otherPhase = ipWater;
817  phaseCorrectionNeeded = true;
818  }
819  // if ( ipOil >= 0 && LvArray::math::abs( phaseFrac[ipOil] ) < MIN_FOR_PHASE_PRESENCE
820  // && LvArray::math::abs( phaseMinVolumeFraction[ipOil] ) > MIN_FOR_PHASE_PRESENCE )
821  // {
822  // addedCompFrac[0][0] = 1.0; // hard-coded (assumes some hydrocarbon is at index 0)
823  // ip_phase = ipOil;
824  // phaseCorrectionNeeded = true;
825  // }
826 
827  if( ipWater >= 0 && LvArray::math::abs( phaseFrac[ipWater] ) < MIN_FOR_PHASE_PRESENCE
828  && LvArray::math::abs( phaseMinVolumeFraction[ipWater] ) > MIN_FOR_PHASE_PRESENCE )
829  {
830  addedCompFrac[0][1] = 1.0; // hard-coded (assumes H2O is at index 1)
831  ip_phase = ipWater;
832  ip_otherPhase = ipGas;
833  phaseCorrectionNeeded = true;
834  }
835 
836  if( phaseCorrectionNeeded )
837  {
838  return applyPhaseCorrection( numComps,
839  numPhases,
840  ip_phase,
841  ip_otherPhase,
842  pres,
843  temp,
844  inputComposition,
845  addedCompFrac[0],
846  compFrac,
847  fluidWrapper );
848  }
849  else
850  {
851  for( localIndex ic = 0; ic < numComps; ++ic )
852  {
853  compFrac[ic] = inputComposition[ic];
854  }
855  return ReturnType::PHASE_CORRECTION_NOT_NEEDED;
856  }
857 
858  }
859 
860  template< typename FLUID_WRAPPER >
861  static ReturnType
862  applyPhaseCorrection( integer const & numComps,
863  integer const & numPhases,
864  integer const & ip_phase,
865  integer const & ip_otherPhase,
866  real64 const & pres,
867  real64 const & temp,
868  arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & inputComposition,
869  arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & addedCompFrac,
870  arraySlice1d< real64, compflow::USD_COMP - 1 > const & compFrac,
871  FLUID_WRAPPER fluidWrapper )
872  {
873  // flash inputs
874  array3d< real64, constitutive::multifluid::LAYOUT_PHASE > phaseFrac( 1, 1, numPhases );
875  array3d< real64, constitutive::multifluid::LAYOUT_PHASE > phaseDens( 1, 1, numPhases );
876  array3d< real64, constitutive::multifluid::LAYOUT_PHASE > phaseMassDens( 1, 1, numPhases );
877  array3d< real64, constitutive::multifluid::LAYOUT_PHASE > phaseVisc( 1, 1, numPhases );
878  array3d< real64, constitutive::multifluid::LAYOUT_PHASE > phaseEnthalpy( 1, 1, numPhases );
879  array3d< real64, constitutive::multifluid::LAYOUT_PHASE > phaseInternalEnergy( 1, 1, numPhases );
880  array4d< real64, constitutive::multifluid::LAYOUT_PHASE_COMP > phaseCompFrac( 1, 1, numPhases, numComps );
881  real64 totalDens = 0.0;
882 
883  real64 a_low = 0.0;
884  real64 a_high = 1.0;
885  real64 a = 0.0;
886  int maxIters = 100;
887  real64 targetPhaseFrac = 1e-4;
888  real64 err = 1e10;
889  for( int iter = 0; iter < maxIters; ++iter )
890  {
891  a = ( a_high + a_low ) * 0.5;
892  mixingStep( numComps,
893  a,
894  inputComposition,
895  addedCompFrac,
896  compFrac );
897  constitutive::MultiFluidBase::KernelWrapper::computeValues( fluidWrapper,
898  pres,
899  temp,
900  compFrac,
901  phaseFrac[0][0],
902  phaseDens[0][0],
903  phaseMassDens[0][0],
904  phaseVisc[0][0],
905  phaseEnthalpy[0][0],
906  phaseInternalEnergy[0][0],
907  phaseCompFrac[0][0],
908  totalDens );
909  err = ( phaseFrac[0][0][ip_phase] - targetPhaseFrac ) / targetPhaseFrac;
910  if( LvArray::math::abs( phaseFrac[0][0][ip_otherPhase] - 1.0 ) < MIN_FOR_PHASE_PRESENCE )
911  {
912  a_low = a;
913  }
914  else if( LvArray::math::abs( phaseFrac[0][0][ip_phase] - 1.0 ) < MIN_FOR_PHASE_PRESENCE )
915  {
916  a_high = a;
917  }
918  else
919  {
920  if( err > 0 )
921  {
922  a_high = a;
923  }
924  else
925  {
926  a_low = a;
927  }
928  if( LvArray::math::abs( err ) < 1e-5 )
929  {
930  break;
931  }
932  }
933  }
934 
935  return ReturnType::SUCCESS;
936  }
937 
938  static void mixingStep( integer const & numComps,
939  real64 const & a,
940  arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & inputComposition,
941  arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & addedCompFrac,
942  arraySlice1d< real64, compflow::USD_COMP - 1 > const & compFrac )
943  {
944  real64 tot = 0.0;
945  for( localIndex ic = 0; ic < numComps; ++ic )
946  {
947  compFrac[ic] = a * addedCompFrac[ic] + ( 1 - a ) * inputComposition[ic];
948  if( compFrac[ic] < MIN_FOR_COMP_PRESENCE )
949  {
950  compFrac[ic] = 0.0;
951  }
952  tot += compFrac[ic];
953  }
954  for( localIndex ic = 0; ic < numComps; ++ic )
955  {
956  compFrac[ic] /= tot;
957  }
958  }
959 
960  static integer evaluateFlashPhaseIndex( integer const & numPhases,
961  integer const & ipGas,
962  integer const & ipOil,
963  integer const & ipWater,
964  integer const & ipInit,
965  real64 const & elevation,
966  arrayView1d< real64 const > const & phaseContacts )
967  {
968  integer ipFP = -1;
969  if( numPhases == 1 )
970  {
971  if( ipInit != -1 )
972  {
973  ipFP = ipInit;
974  }
975  else
976  {
977  ipFP = ( ipGas >= 0 ) ? ipGas :
978  ( ipOil >= 0 ) ? ipOil :
979  ( ipWater >= 0 ) ? ipWater : 0;
980  }
981  }
982  else if( ipGas >= 0 && ipOil >= 0 && ipWater >= 0 )
983  {
984  // phases = gas + oil + water
985  real64 const goc = phaseContacts[1];
986  real64 const owc = phaseContacts[0];
987  ipFP = ipWater; // default
988  if( elevation > owc
989  && elevation < goc )
990  {
991  ipFP = ipOil;
992  }
993  else if( elevation > goc
994  || LvArray::math::abs( elevation - goc ) < 1e-12 )
995  {
996  ipFP = ipGas;
997  }
998  }
999  else if( ipGas >= 0 && ipWater >= 0 )
1000  {
1001  // phases = gas + water
1002  real64 const gwc = phaseContacts[0];
1003  ipFP = ipWater; // default
1004  if( elevation > gwc
1005  || LvArray::math::abs( elevation - gwc ) < 1e-12 )
1006  {
1007  ipFP = ipGas;
1008  }
1009  }
1010  else if( ipOil >= 0 && ipWater >= 0 )
1011  {
1012  // phases = oil + water
1013  real64 const owc = phaseContacts[0];
1014  ipFP = ipWater; // default
1015  if( elevation > owc
1016  || LvArray::math::abs( elevation - owc ) < 1e-12 )
1017  {
1018  ipFP = ipOil;
1019  }
1020  }
1021  else if( ipGas >= 0 && ipOil >= 0 )
1022  {
1023  // phases = gas + oil
1024  real64 const goc = phaseContacts[0];
1025  ipFP = ipOil; // default
1026  if( elevation > goc
1027  || LvArray::math::abs( elevation - goc ) < 1e-12 )
1028  {
1029  ipFP = ipGas;
1030  }
1031  }
1032  return ipFP;
1033  }
1034 
1035  static void evaluatePrimaryAndContactPhaseIndices( integer const & numPhases,
1036  integer const & ipGas,
1037  integer const & ipOil,
1038  integer const & ipWater,
1039  real64 const & startElevation,
1040  real64 const & endElevation,
1041  arrayView1d< real64 const > const & phaseContacts,
1042  integer & ipPP,
1043  integer & ipCP )
1044  {
1045  ipPP = -1;
1046  ipCP = -1;
1047  if( numPhases > 1 )
1048  {
1049  if( ipGas >= 0 && ipOil >= 0 && ipWater >= 0 )
1050  {
1051  // phases = gas + oil + water
1052  real64 const goc = phaseContacts[1];
1053  real64 const owc = phaseContacts[0];
1054  ipPP = ipWater; // default
1055  ipCP = ipOil; // default
1056  if( LvArray::math::abs( endElevation - owc ) < 1e-12
1057  && startElevation > owc )
1058  {
1059  ipPP = ipOil;
1060  ipCP = ipWater;
1061  }
1062  else if( LvArray::math::abs( endElevation - goc ) < 1e-12 )
1063  {
1064  if( startElevation < goc )
1065  {
1066  ipPP = ipOil;
1067  ipCP = ipGas;
1068  }
1069  else
1070  {
1071  ipPP = ipGas;
1072  ipCP = ipOil;
1073  }
1074  }
1075  }
1076  else if( ipGas >= 0 && ipWater >= 0 )
1077  {
1078  // phases = gas + water
1079  real64 const gwc = phaseContacts[0];
1080  ipPP = ipWater; // default
1081  ipCP = ipGas; // default
1082  if( startElevation > gwc )
1083  {
1084  ipPP = ipGas;
1085  ipCP = ipWater;
1086  }
1087  }
1088  else if( ipOil >= 0 && ipWater >= 0 )
1089  {
1090  // phases = oil + water
1091  real64 const owc = phaseContacts[0];
1092  ipPP = ipWater; // default
1093  ipCP = ipOil; // default
1094  if( startElevation > owc )
1095  {
1096  ipPP = ipOil;
1097  ipCP = ipWater;
1098  }
1099  }
1100  else if( ipGas >= 0 && ipOil >= 0 )
1101  {
1102  // phases = gas + oil
1103  real64 const goc = phaseContacts[0];
1104  ipPP = ipOil; // default
1105  ipCP = ipGas; // default
1106  if( startElevation > goc )
1107  {
1108  ipPP = ipGas;
1109  ipCP = ipOil;
1110  }
1111  }
1112  }
1113  }
1114 
1115 
1116  template< typename FLUID_WRAPPER >
1117  static void
1118  computeDatumPhaseMassDens( integer const & numComps,
1119  integer const & numPhases,
1120  integer const & ipGas,
1121  integer const & ipWater,
1122  arrayView1d< real64 const > const & phaseMinVolumeFraction,
1123  real64 const & datumElevation,
1124  real64 const & datumPres,
1125  real64 const & datumTemp,
1126  arrayView3d< real64, constitutive::multifluid::USD_PHASE > const & datumPhaseMassDens,
1127  arrayView1d< TableFunction::KernelWrapper const > compFracTableWrappers,
1128  bool singlePhase,
1129  FLUID_WRAPPER fluidWrapper )
1130  {
1131  // datum fluid properties
1132  array2d< real64, compflow::LAYOUT_COMP > datumInputComposition( 1, numComps );
1133  array2d< real64, compflow::LAYOUT_COMP > datumCompFrac( 1, numComps );
1134  array3d< real64, constitutive::multifluid::LAYOUT_PHASE > datumPhaseFrac( 1, 1, numPhases );
1135  array3d< real64, constitutive::multifluid::LAYOUT_PHASE > datumPhaseDens( 1, 1, numPhases );
1136  array3d< real64, constitutive::multifluid::LAYOUT_PHASE > datumPhaseVisc( 1, 1, numPhases );
1137  array3d< real64, constitutive::multifluid::LAYOUT_PHASE > datumPhaseEnthalpy( 1, 1, numPhases );
1138  array3d< real64, constitutive::multifluid::LAYOUT_PHASE > datumPhaseInternalEnergy( 1, 1, numPhases );
1139  array4d< real64, constitutive::multifluid::LAYOUT_PHASE_COMP > datumPhaseCompFrac( 1, 1, numPhases, numComps );
1140  real64 datumTotalDens = 0.0;
1141  for( integer ic = 0; ic < numComps; ++ic )
1142  {
1143  datumInputComposition[0][ic] = compFracTableWrappers[ic].compute( &datumElevation );
1144  }
1145  constitutive::MultiFluidBase::KernelWrapper::computeValues( fluidWrapper,
1146  datumPres,
1147  datumTemp,
1148  datumInputComposition[0],
1149  datumPhaseFrac[0][0],
1150  datumPhaseDens[0][0],
1151  datumPhaseMassDens[0][0],
1152  datumPhaseVisc[0][0],
1153  datumPhaseEnthalpy[0][0],
1154  datumPhaseInternalEnergy[0][0],
1155  datumPhaseCompFrac[0][0],
1156  datumTotalDens );
1157 
1158  if( singlePhase )
1159  {
1160  return;
1161  }
1162 
1163  ReturnType datumPhaseCorr = phaseCorrection( numComps,
1164  numPhases,
1165  ipGas,
1166  ipWater,
1167  phaseMinVolumeFraction,
1168  datumPres,
1169  datumTemp,
1170  datumInputComposition[0],
1171  datumPhaseFrac[0][0],
1172  datumCompFrac[0],
1173  fluidWrapper );
1174  if( datumPhaseCorr == ReturnType::SUCCESS )
1175  constitutive::MultiFluidBase::KernelWrapper::computeValues( fluidWrapper,
1176  datumPres,
1177  datumTemp,
1178  datumCompFrac[0],
1179  datumPhaseFrac[0][0],
1180  datumPhaseDens[0][0],
1181  datumPhaseMassDens[0][0],
1182  datumPhaseVisc[0][0],
1183  datumPhaseEnthalpy[0][0],
1184  datumPhaseInternalEnergy[0][0],
1185  datumPhaseCompFrac[0][0],
1186  datumTotalDens );
1187  }
1188 
1189 };
1190 
1191 } // namespace isothermalCompositionalMultiphaseBaseKernels
1192 
1193 } // namespace geos
1194 
1195 
1196 #endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_HYDROSTATICPRESSUREKERNEL_HPP
geos::TableFunction::KernelWrapper::compute
GEOS_HOST_DEVICE real64 compute(IN_ARRAY const &input) const
Interpolate in the table.
geos::arrayView1d
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:179
geos::array2d
Array< T, 2, PERMUTATION > array2d
Alias for 2D array.
Definition: DataTypes.hpp:191
geos::array3d
Array< T, 3, PERMUTATION > array3d
Alias for 3D array.
Definition: DataTypes.hpp:207
geos::StackArray
LvArray::StackArray< T, NDIM, PERMUTATION, localIndex, MAXSIZE > StackArray
Multidimensional stack-based array type. See LvArray:StackArray for details.
Definition: DataTypes.hpp:155
geos::arraySlice2d
ArraySlice< T, 2, USD > arraySlice2d
Alias for 2D array slice.
Definition: DataTypes.hpp:199
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:183
geos::array4d
Array< T, 4, PERMUTATION > array4d
Alias for 4D array.
Definition: DataTypes.hpp:223
geos::arrayView4d
ArrayView< T, 4, USD > arrayView4d
Alias for 4D array view.
Definition: DataTypes.hpp:227
geos::arrayView2d
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:195
geos::integer
int integer
Signed integer type.
Definition: DataTypes.hpp:81
geos::array1d
Array< T, 1 > array1d
Alias for 1D array.
Definition: DataTypes.hpp:175
geos::arrayView3d
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
Definition: DataTypes.hpp:211
Generated by doxygen 1.9.1