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