Skip to content

Commit

Permalink
feat: feature/byer3/thermal well (#3156)
Browse files Browse the repository at this point in the history
* Add support for mass constraits (inj wells only)

* Throw error if massinj constraint and useSurfaceConditions set to 0, fix error message for negative mass rate

* thermal devs- 1) energy eqn for perf in/out flows, 2) reservoir well coupling 3) residualnomr calcs, 4) some isothermal , 5) refactoring  6) next focus solution updating/scaling

* 1) fix isthermal well field runtime error 2) fix refactored isothermal flux index bug 3) 2 isothermal models tested with refactored code 4) thermal case crashes at 2nd linear solve well internal energy derivative = 0

* thermal well devs - 1) derivative bugfixes 2) added thermal to post linear solve well methods 3) next step verify simple test case results

* thermal devs - 1) fixes for > 1 well seg 2) added soln updating/scaling 3) removed debug std::couts 4) next steps a) address constant temp BC b) singlephase thermo dev c) deriv checker

* 1) Added simple temperature constraint 2) Fixes from running numerical derivative checker and treating warnings as errors

* 1) add iterative solver config for thermal res+well  2) skip well calcs for shutin well 3) call well initialization code for wells that come online at t> 0

* 1) add Aww numerical derivative unit tests, 2) thermal index fix

* 1) fix temp constraint eqn setup for well spanning > 1 cores , 2) well debug moved to well_debug branch, 3) stored global indices for well segs, 4) stored global reservoir indices for perfs, 5) added interfaces to access newton iteration number

* 1) remove schema addition for current nonlinear stats , 2) remove test code that zeroed energy derivatives if phase not present

* 1) for shutin wells zero completion perforation rates for clarity in reporting, 2) initialize seg temps with input inj temp, if user specifed surface+inj temp this caused differences with dev for isothermal runs

* Update src/coreComponents/unitTests/fluidFlowTests/testCompFlowUtils.hpp

* 1) dont allocate temperature_n for isothermal run, 2) uncrustify style

* perfs global res cell index set with wrong array (node instead of cell values),  quantites only used for reporting

* Update .integrated_tests.yaml
  • Loading branch information
tjb-ltk authored Oct 31, 2024
1 parent cd8175d commit 890c5e7
Show file tree
Hide file tree
Showing 52 changed files with 9,386 additions and 2,567 deletions.
2 changes: 1 addition & 1 deletion .integrated_tests.yaml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
baselines:
bucket: geosx
baseline: integratedTests/baseline_integratedTests-pr2878-8188-ed2ded2
baseline: integratedTests/baseline_integratedTests-pr3156-8366-aa71f2a

allow_fail:
all: ''
Expand Down
4 changes: 4 additions & 0 deletions BASELINE_NOTES.md
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines.
Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining.
These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD).

PR #3156 (2024-10-29)
====================
Restart check errors due to 1) schema node added to enable thermal option in well model and 2) arrays removed/added for option. Max difference errors due treatment of shutin wells. Previously non-zero rate value reported for shutin well, new code will set rate arrays to zero.

PR #2878 (2024-10-17)
=====================
Sorted region cellBlocks names alphabetically. Therefore affected ordering of: faceManager/elemSubRegionList, nodeManager/elemList, nodeManager/elemSubRegionList, SurfaceElementSubRegion::fractureElementsToCellSubRegions, field::perforation::reservoirElementSubregion.
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -386,7 +386,6 @@ CO2BrineFluid< PHASE1, PHASE2, FLASH >::KernelWrapper::

if( m_isThermal )
{

m_phase1.enthalpy.compute( pressure,
temperatureInCelsius,
phaseCompFraction.value[ip1].toSliceConst(), phaseCompFraction.derivs[ip1].toSliceConst(),
Expand Down
40 changes: 40 additions & 0 deletions src/coreComponents/constitutive/fluid/multifluid/Layouts.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -30,10 +30,23 @@ namespace geos
{
namespace constitutive
{

namespace singlefluid
{
struct DerivativeOffset
{
/// index of derivative wrt pressure
static integer constexpr dP = 0;
/// index of derivative wrt temperature
static integer constexpr dT = 1;

};
}
namespace multifluid
{

/// indices of pressure, temperature, and composition derivatives
// Fix me - if the order is changed the code crashes
struct DerivativeOffset
{
/// index of derivative wrt pressure
Expand All @@ -44,6 +57,33 @@ struct DerivativeOffset
static integer constexpr dC = 2;
};

/// indices of pressure, temperature, and composition derivatives
template< integer NC, integer IS_THERMAL >
struct DerivativeOffsetC {};

template< integer NC >
struct DerivativeOffsetC< NC, 1 >
{
/// index of derivative wrt pressure
static integer constexpr dP = 0;
/// index of derivative wrt temperature
static integer constexpr dT = dP + 1;
/// index of first derivative wrt compositions
static integer constexpr dC = dP+2;
/// number of derivatives
static integer constexpr nDer = NC + 2;
};
template< integer NC >
struct DerivativeOffsetC< NC, 0 >
{
/// index of derivative wrt pressure
static integer constexpr dP = 0;
/// index of first derivative wrt compositions
static integer constexpr dC = dP+1;
/// number of derivatives
static integer constexpr nDer = NC + 1;
};

#if defined( GEOS_USE_DEVICE )

/// Constitutive model phase property array layout
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -276,6 +276,14 @@ class SingleFluidBase : public ConstitutiveBase
virtual real64 defaultDensity() const = 0;
virtual real64 defaultViscosity() const = 0;

/**
* @brief Get the thermal flag.
* @return boolean value indicating whether the model can be used to assemble the energy balance equation or not
* @detail if isThermal is true, the constitutive model compute the enthalpy and internal energy of the phase.
* This can be used to check the compatibility of the constitutive model with the solver
*/
virtual bool isThermal() const { return false; }

protected:

virtual void postInputInitialization() override;
Expand Down
2 changes: 1 addition & 1 deletion src/coreComponents/dataRepository/Group.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -457,7 +457,6 @@ class Group
{
using T = std::conditional_t< std::is_const< CONTAINERTYPE >::value, CASTTYPE const, CASTTYPE >;
T * const castedContainer = dynamic_cast< T * >( &container );

if( castedContainer != nullptr )
{
lambda( *castedContainer );
Expand Down Expand Up @@ -598,6 +597,7 @@ class Group
void forSubGroups( LOOKUP_CONTAINER const & subGroupKeys, LAMBDA && lambda )
{
localIndex counter = 0;

for( auto const & subgroup : subGroupKeys )
{
applyLambdaToContainer< GROUPTYPE, GROUPTYPES... >( getGroup( subgroup ), [&]( auto & castedSubGroup )
Expand Down
1 change: 1 addition & 0 deletions src/coreComponents/linearAlgebra/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,7 @@ if( ENABLE_HYPRE )
interfaces/hypre/mgrStrategies/SinglePhaseReservoirHybridFVM.hpp
interfaces/hypre/mgrStrategies/SolidMechanicsEmbeddedFractures.hpp
interfaces/hypre/mgrStrategies/ThermalCompositionalMultiphaseFVM.hpp
interfaces/hypre/mgrStrategies/ThermalCompositionalMultiphaseReservoirFVM.hpp
interfaces/hypre/mgrStrategies/ThermalMultiphasePoromechanics.hpp
interfaces/hypre/mgrStrategies/ThermalSinglePhasePoromechanics.hpp )
list( APPEND linearAlgebra_sources
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@
#include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirFVM.hpp"
#include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirHybridFVM.hpp"
#include "linearAlgebra/interfaces/hypre/mgrStrategies/ThermalCompositionalMultiphaseFVM.hpp"
#include "linearAlgebra/interfaces/hypre/mgrStrategies/ThermalCompositionalMultiphaseReservoirFVM.hpp"
#include "linearAlgebra/interfaces/hypre/mgrStrategies/ThermalSinglePhasePoromechanics.hpp"
#include "linearAlgebra/interfaces/hypre/mgrStrategies/ThermalMultiphasePoromechanics.hpp"
#include "linearAlgebra/interfaces/hypre/mgrStrategies/SolidMechanicsEmbeddedFractures.hpp"
Expand Down Expand Up @@ -112,6 +113,11 @@ void hypre::mgr::createMGR( LinearSolverParameters const & params,
setStrategy< ThermalCompositionalMultiphaseFVM >( params.mgr, numComponentsPerField, precond, mgrData );
break;
}
case LinearSolverParameters::MGR::StrategyType::thermalCompositionalMultiphaseReservoirFVM:
{
setStrategy< ThermalCompositionalMultiphaseReservoirFVM >( params.mgr, numComponentsPerField, precond, mgrData );
break;
}
case LinearSolverParameters::MGR::StrategyType::hybridSinglePhasePoromechanics:
{
setStrategy< HybridSinglePhasePoromechanics >( params.mgr, numComponentsPerField, precond, mgrData );
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,131 @@
/*
* ------------------------------------------------------------------------------------------------------------
* SPDX-License-Identifier: LGPL-2.1-only
*
* Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
* Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
* Copyright (c) 2018-2020 TotalEnergies
* Copyright (c) 2019- GEOSX Contributors
* All rights reserved
*
* See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
* ------------------------------------------------------------------------------------------------------------
*/

/**
* @file ThermalCompositionalMultiphaseFVM.hpp
*/

#ifndef GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRTHERMALCOMPOSITIONALMULTIPHASERESERVOIRFVM_HPP_
#define GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRTHERMALCOMPOSITIONALMULTIPHASERESERVOIRFVM_HPP_

#include "linearAlgebra/interfaces/hypre/HypreMGR.hpp"

namespace geos
{

namespace hypre
{

namespace mgr
{

/**
* @brief ThermalCompositionalMultiphaseReservoirFVM strategy.
*
* Labels description stored in point_marker_array
* 0 = reservoir pressure
* 1 = reservoir density (component 1 )
* ... = reservoir densities (# components , NC)
* NC + 1 = reservoir temperature
*
* numResLabels - 1 = reservoir temperature
* numResLabels = well pressure
* numResLabels + 1 = well density
* ... = ... (well densities)
* numResLabels + numWellLabels - 3 = last well density
* numResLabels + numWellLabels - 2 = well rate
* numResLabels + numWellLabels - 1 = well temperature
*
* 3-level MGR reduction strategy
* - 1st level: eliminate the well block
* - 2nd level: eliminate the reservoir density associated with the volume constraint
* - 3rd level: eliminate the remaining the reservoir densities
* - The coarse grid (pressure and temperature system) is solved with BoomerAMG with numFunctions==2.
*
*/

class ThermalCompositionalMultiphaseReservoirFVM : public MGRStrategyBase< 3 >
{
public:
/**
* @brief Constructor.
* @param numComponentsPerField array with number of components for each field
*/
explicit ThermalCompositionalMultiphaseReservoirFVM( arrayView1d< int const > const & numComponentsPerField )
: MGRStrategyBase( LvArray::integerConversion< HYPRE_Int >( numComponentsPerField[0] + numComponentsPerField[1] ) )
{
HYPRE_Int const numResLabels = LvArray::integerConversion< HYPRE_Int >( numComponentsPerField[0] );

// Level 0: eliminate the well block
m_labels[0].resize( numResLabels );
std::iota( m_labels[0].begin(), m_labels[0].end(), 0 );

// Level 1: eliminate last density which corresponds to the volume constraint equation
m_labels[1].resize( numResLabels - 2 );
std::iota( m_labels[1].begin(), m_labels[1].end(), 0 );
m_labels[1].push_back( numResLabels-1 ); // keep temperature
// Level 2: eliminate the other densities
m_labels[2].push_back( 0 ); // keep pressure
m_labels[2].push_back( numResLabels-1 ); // keep temperature

setupLabels();

// level 0
m_levelFRelaxType[0] = MGRFRelaxationType::gsElimWInverse;
m_levelFRelaxIters[0] = 1;
m_levelInterpType[0] = MGRInterpolationType::blockJacobi;
m_levelRestrictType[0] = MGRRestrictionType::injection;
m_levelCoarseGridMethod[0] = MGRCoarseGridMethod::galerkin;
m_levelGlobalSmootherType[0] = MGRGlobalSmootherType::none;

m_levelFRelaxType[1] = MGRFRelaxationType::jacobi;
m_levelFRelaxIters[1] = 1;
m_levelInterpType[1] = MGRInterpolationType::jacobi; // Diagonal scaling (Jacobi)
m_levelRestrictType[1] = MGRRestrictionType::injection;
m_levelCoarseGridMethod[1] = MGRCoarseGridMethod::galerkin; // Standard Galerkin
m_levelGlobalSmootherType[1] = MGRGlobalSmootherType::blockGaussSeidel;
m_levelGlobalSmootherIters[1] = 1;

m_levelFRelaxType[2] = MGRFRelaxationType::jacobi;
m_levelFRelaxIters[2] = 1;
m_levelInterpType[2] = MGRInterpolationType::injection; // Injection
m_levelRestrictType[2] = MGRRestrictionType::injection;
m_levelCoarseGridMethod[2] = MGRCoarseGridMethod::cprLikeBlockDiag; // Non-Galerkin Quasi-IMPES CPR
m_levelGlobalSmootherType[2] = MGRGlobalSmootherType::ilu0;
m_levelGlobalSmootherIters[2] = 1;
}

/**
* @brief Setup the MGR strategy.
* @param precond preconditioner wrapper
* @param mgrData auxiliary MGR data
*/
void setup( LinearSolverParameters::MGR const &,
HyprePrecWrapper & precond,
HypreMGRData & mgrData )
{
setReduction( precond, mgrData );

// Configure the BoomerAMG solver used as mgr coarse solver for the pressure/temperature reduced system
setPressureTemperatureAMG( mgrData.coarseSolver );
}
};

} // namespace mgr

} // namespace hypre

} // namespace geos

#endif /*GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRCOMPOSITIONALMULTIPHASERESERVOIRFVM_HPP_*/
Original file line number Diff line number Diff line change
Expand Up @@ -271,28 +271,29 @@ struct LinearSolverParameters
*/
enum class StrategyType : integer
{
invalid, ///< default value, to ensure solver sets something
singlePhaseReservoirFVM, ///< finite volume single-phase flow with wells
singlePhaseHybridFVM, ///< hybrid finite volume single-phase flow
singlePhaseReservoirHybridFVM, ///< hybrid finite volume single-phase flow with wells
singlePhasePoromechanics, ///< single phase poromechanics with finite volume single phase flow
thermalSinglePhasePoromechanics, ///< thermal single phase poromechanics with finite volume single phase flow
hybridSinglePhasePoromechanics, ///< single phase poromechanics with hybrid finite volume single phase flow
singlePhasePoromechanicsEmbeddedFractures, ///< single phase poromechanics with FV embedded fractures
singlePhasePoromechanicsConformingFractures, ///< single phase poromechanics with FV conforming fractures
singlePhasePoromechanicsReservoirFVM, ///< single phase poromechanics with finite volume single phase flow with wells
compositionalMultiphaseFVM, ///< finite volume compositional multiphase flow
compositionalMultiphaseHybridFVM, ///< hybrid finite volume compositional multiphase flow
compositionalMultiphaseReservoirFVM, ///< finite volume compositional multiphase flow with wells
compositionalMultiphaseReservoirHybridFVM, ///< hybrid finite volume compositional multiphase flow with wells
reactiveCompositionalMultiphaseOBL, ///< finite volume reactive compositional flow with OBL
thermalCompositionalMultiphaseFVM, ///< finite volume thermal compositional multiphase flow
multiphasePoromechanics, ///< multiphase poromechanics with finite volume compositional multiphase flow
multiphasePoromechanicsReservoirFVM, ///< multiphase poromechanics with finite volume compositional multiphase flow with wells
thermalMultiphasePoromechanics, ///< thermal multiphase poromechanics with finite volume compositional multiphase flow
hydrofracture, ///< hydrofracture
lagrangianContactMechanics, ///< Lagrangian contact mechanics
solidMechanicsEmbeddedFractures ///< Embedded fractures mechanics
invalid, ///< default value, to ensure solver sets something
singlePhaseReservoirFVM, ///< finite volume single-phase flow with wells
singlePhaseHybridFVM, ///< hybrid finite volume single-phase flow
singlePhaseReservoirHybridFVM, ///< hybrid finite volume single-phase flow with wells
singlePhasePoromechanics, ///< single phase poromechanics with finite volume single phase flow
thermalSinglePhasePoromechanics, ///< thermal single phase poromechanics with finite volume single phase flow
hybridSinglePhasePoromechanics, ///< single phase poromechanics with hybrid finite volume single phase flow
singlePhasePoromechanicsEmbeddedFractures, ///< single phase poromechanics with FV embedded fractures
singlePhasePoromechanicsConformingFractures, ///< single phase poromechanics with FV embedded fractures
singlePhasePoromechanicsReservoirFVM, ///< single phase poromechanics with finite volume single phase flow with wells
compositionalMultiphaseFVM, ///< finite volume compositional multiphase flow
compositionalMultiphaseHybridFVM, ///< hybrid finite volume compositional multiphase flow
compositionalMultiphaseReservoirFVM, ///< finite volume compositional multiphase flow with wells
compositionalMultiphaseReservoirHybridFVM, ///< hybrid finite volume compositional multiphase flow with wells
reactiveCompositionalMultiphaseOBL, ///< finite volume reactive compositional flow with OBL
thermalCompositionalMultiphaseFVM, ///< finite volume thermal compositional multiphase flow
thermalCompositionalMultiphaseReservoirFVM,///< finite volume thermal compositional multiphase flow
multiphasePoromechanics, ///< multiphase poromechanics with finite volume compositional multiphase flow
multiphasePoromechanicsReservoirFVM, ///< multiphase poromechanics with finite volume compositional multiphase flow with wells
thermalMultiphasePoromechanics, ///< thermal multiphase poromechanics with finite volume compositional multiphase flow
hydrofracture, ///< hydrofracture
lagrangianContactMechanics, ///< Lagrangian contact mechanics
solidMechanicsEmbeddedFractures ///< Embedded fractures mechanics
};

StrategyType strategy = StrategyType::invalid; ///< Predefined MGR solution strategy (solver specific)
Expand Down Expand Up @@ -379,6 +380,7 @@ ENUM_STRINGS( LinearSolverParameters::MGR::StrategyType,
"compositionalMultiphaseReservoirHybridFVM",
"reactiveCompositionalMultiphaseOBL",
"thermalCompositionalMultiphaseFVM",
"thermalCompositionalMultiphaseReservoirFVM",
"multiphasePoromechanics",
"multiphasePoromechanicsReservoirFVM",
"thermalMultiphasePoromechanics",
Expand Down
2 changes: 2 additions & 0 deletions src/coreComponents/mesh/PerforationData.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,8 @@ PerforationData::PerforationData( string const & name, Group * const parent )
registerField( fields::perforation::reservoirElementRegion{}, &m_toMeshElements.m_toElementRegion );
registerField( fields::perforation::reservoirElementSubRegion{}, &m_toMeshElements.m_toElementSubRegion );
registerField( fields::perforation::reservoirElementIndex{}, &m_toMeshElements.m_toElementIndex );
registerField( fields::perforation::reservoirElementGlobalIndex{}, &m_reservoirElementGlobalIndex );

registerField( fields::perforation::wellElementIndex{}, &m_wellElementIndex );
registerField( fields::perforation::location{}, &m_location );
registerField( fields::perforation::wellTransmissibility{}, &m_wellTransmissibility );
Expand Down
17 changes: 17 additions & 0 deletions src/coreComponents/mesh/PerforationData.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -150,6 +150,20 @@ class PerforationData : public ObjectManagerBase
*/
arrayView1d< localIndex const > getWellElements() const { return m_wellElementIndex; }

/**
* @brief Get perforation-to-reservoir-element connectivity.
* @return list of global reservoir element index connected to each perforation
*/
arrayView1d< globalIndex > getReservoirElementGlobalIndex() { return m_reservoirElementGlobalIndex; }



/**
* @brief Provide an immutable accessor to a const perforation-to-reservoir-element connectivity.
* @return list of well element index connected to each perforation
*/
arrayView1d< globalIndex const > getReservoirElementGlobalIndex() const { return m_reservoirElementGlobalIndex; }


/**
* @brief Get perforation locations.
Expand Down Expand Up @@ -275,6 +289,9 @@ class PerforationData : public ObjectManagerBase
/// Indices of the well elements to which perforations are attached
array1d< localIndex > m_wellElementIndex;

/// Global indices of reservoir cell containing perforation
array1d< globalIndex > m_reservoirElementGlobalIndex;

/// Location of the perforations
array2d< real64 > m_location;

Expand Down
8 changes: 8 additions & 0 deletions src/coreComponents/mesh/PerforationFields.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,14 @@ DECLARE_FIELD( reservoirElementIndex,
WRITE_AND_READ,
"For each perforation, element index of the perforated element" );

DECLARE_FIELD( reservoirElementGlobalIndex,
"reservoirElementGlobalIndex",
array1d< globalIndex >,
0,
NOPLOT,
WRITE_AND_READ,
"For each perforation, global element index of the perforated element" );

DECLARE_FIELD( wellElementIndex,
"wellElementIndex",
array1d< localIndex >,
Expand Down
Loading

0 comments on commit 890c5e7

Please sign in to comment.