WellControls.hpp

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 TotalEnergies
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
/*
 * @file WellControls.hpp
 */
 
#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_WELLCONTROLS_HPP
#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_WELLCONTROLS_HPP
 
#include "physicsSolvers/PhysicsSolverBase.hpp"
#include "common/format/EnumStrings.hpp"
#include "dataRepository/Group.hpp"
#include "functions/TableFunction.hpp"
#include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
 
#include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
#include "constitutive/fluid/singlefluid/SingleFluidBase.hpp"
#include "physicsSolvers/fluidFlow/wells/WellConstraintsBase.hpp"
#include "physicsSolvers/fluidFlow/wells/WellNewtonSolver.hpp"
 
namespace geos
{
namespace dataRepository
{
namespace keys
{
static constexpr auto wellControls = "WellControls";
}
}
 
class ElementsReporterBuffer;
 
 
class WellControls :  public dataRepository::Group
{
public:
 
  enum class Type : integer
  {
    PRODUCER,  
    INJECTOR   
  };
 
  enum class Status : integer
  {
    OPEN,  
    CLOSED   
  };
 
 
 
  explicit WellControls( string const & name, dataRepository::Group * const parent );
 
 
  ~WellControls() override;
 
  WellControls() = delete;
 
  WellControls( WellControls const & ) = delete;
 
  WellControls( WellControls && ) = delete;
 
  WellControls & operator=( WellControls const & ) = delete;
 
  WellControls & operator=( WellControls && ) = delete;
 
 
 
  static string coupledSolverAttributePrefix() { return "well"; }
 
  virtual Group * createChild( string const & childKey, string const & childName ) override;
 
  virtual void expandObjectCatalogs() override;
 
  virtual void registerWellDataOnMesh( WellElementSubRegion & subRegion );
  virtual void setConstitutiveNames( ElementSubRegionBase & subRegion ) const = 0;
  virtual void  createSeparator( WellElementSubRegion & subRegion ) = 0;
 
  virtual void validateWellConstraints( real64 const & time_n,
                                        real64 const & dt,
                                        WellElementSubRegion const & subRegion ) = 0;
 
  virtual bool isCompositional() const = 0;
  virtual void initializeWell( DomainPartition & domain, Group & meshBodies, string const & meshBodyName, MeshLevel & mesh, WellElementSubRegion & subRegion, real64 const & time_n ) = 0;
  real64 setNextDt( real64 const & currentTime,
                    real64 const & currentDt,
                    WellElementSubRegion & subRegion );
  // Bring the base class implicitStepSetup into scope to avoid hiding the overloaded virtual function
 
 
  virtual void implicitStepSetup( real64 const & time_n,
                                  real64 const & GEOS_UNUSED_PARAM( dt ),
                                  DomainPartition & domain,
                                  string const & meshBodyName,
                                  ElementRegionManager & elemManager,
                                  WellElementSubRegion & subRegion ) = 0;
 
  virtual void
  implicitStepComplete( real64 const & time,
                        real64 const & dt,
                        WellElementSubRegion const & subRegion ) = 0;
  virtual real64 updateSubRegionState( real64 const time_n,
                                       MeshBody const & meshBody,
                                       ElementRegionManager const & elemManager,
                                       WellElementSubRegion & subRegion ) = 0;
 
  bool evaluateConstraints( real64 const & time_n,
                            WellElementSubRegion & subRegion );
 
  bool evaluateConstraints( real64 const & time_n,
                            real64 const & dt,
                            integer const cycleNumber,
                            integer const coupledIterationNumber,
                            DomainPartition & domain,
                            MeshLevel & mesh,
                            ElementRegionManager & elemManager,
                            WellElementSubRegion & subRegion,
                            DofManager const & dofManager );
 
  void solveConstraint( WellConstraintBase *constraint,
                        real64 const & time_n,
                        real64 const & dt,
                        integer const cycleNumber,
                        integer const coupledIterationNumber,
                        DomainPartition & domain,
                        MeshLevel & mesh,
                        ElementRegionManager & elemManager,
                        WellElementSubRegion & subRegion,
                        DofManager const & dofManager );
 
  void assembleSystem( real64 const & time_n,
                       real64 const & dt,
                       integer const cycleNumber,
                       ElementRegionManager & elemManager,
                       WellElementSubRegion & subRegion,
                       DofManager const & dofManager,
                       CRSMatrixView< real64, globalIndex const > const & localMatrix,
                       arrayView1d< real64 > const & localRhs );
  virtual void assembleWellAccumulationTerms( real64 const & time,
                                              real64 const & dt,
                                              WellElementSubRegion & subRegion,
                                              DofManager const & dofManager,
                                              CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                              arrayView1d< real64 > const & localRhs ) = 0;
  virtual void assembleWellPressureRelations( real64 const & time_n,
                                              real64 const & dt,
                                              WellElementSubRegion const & subRegion,
                                              DofManager const & dofManager,
                                              CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                              arrayView1d< real64 > const & localRhs ) = 0;
  virtual void assembleWellConstraintTerms( real64 const & time_n,
                                            real64 const & dt,
                                            WellElementSubRegion const & subRegion,
                                            DofManager const & dofManager,
                                            CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                            arrayView1d< real64 > const & localRhs ) = 0;
 
  virtual void computeWellPerforationRates( real64 const & time_n,
                                            real64 const & GEOS_UNUSED_PARAM( dt ),
                                            ElementRegionManager & elemManager,
                                            WellElementSubRegion & subRegion ) = 0;
 
  virtual void assembleWellFluxTerms( real64 const & time,
                                      real64 const & dt,
                                      WellElementSubRegion const & subRegion,
                                      DofManager const & dofManager,
                                      CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                      arrayView1d< real64 > const & localRhs ) = 0;
  virtual real64
  calculateWellResidualNorm( real64 const & time_n,
                             real64 const & dt,
                             NonlinearSolverParameters const & nonlinearSolverParameters,
                             WellElementSubRegion const & subRegion,
                             DofManager const & dofManager,
                             arrayView1d< real64 const > const & localRhs ) = 0;
 
  virtual array1d< real64 >
  calculateLocalWellResidualNorm( real64 const & time_n,
                                  real64 const & dt,
                                  NonlinearSolverParameters const & nonlinearSolverParameters,
                                  WellElementSubRegion const & subRegion,
                                  DofManager const & dofManager,
                                  arrayView1d< real64 const > const & localRhs ) = 0;
 
  virtual real64
  scalingForWellSystemSolution( WellElementSubRegion & subRegion,
                                DofManager const & dofManager,
                                arrayView1d< real64 const > const & localSolution ) = 0;
 
  virtual bool
  checkWellSystemSolution( WellElementSubRegion & subRegion,
                           DofManager const & dofManager,
                           arrayView1d< real64 const > const & localSolution,
                           real64 const scalingFactor,
                           real64 & minPressure,
                           real64 & minDensity,
                           real64 & minTotalDensity,
                           ElementsReporterBuffer & negPressureIds,
                           ElementsReporterBuffer & negDensityIds,
                           ElementsReporterBuffer & negTotalDensityIds ) = 0;
 
  virtual void
  applyWellSystemSolution( DofManager const & dofManager,
                           arrayView1d< real64 const > const & localSolution,
                           real64 const scalingFactor,
                           real64 const dt,
                           DomainPartition & domain,
                           MeshLevel & mesh,
                           WellElementSubRegion & subRegion ) = 0;
 
  virtual void applyWellBoundaryConditions( real64 const time_n,
                                            real64 const dt,
                                            ElementRegionManager & elemManager,
                                            WellElementSubRegion & subRegion,
                                            DofManager const & dofManager,
                                            arrayView1d< real64 > const & localRhs,
                                            CRSMatrixView< real64, globalIndex const > const & localMatrix ) = 0;
  virtual real64 updateWellState( MeshBody const & meshBody,
                                  ElementRegionManager const & elemManager,
                                  WellElementSubRegion & subRegion ) = 0;
 
  virtual void resetStateToBeginningOfStep( DomainPartition & domain,
                                            string const & meshBodyName, ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ) = 0;
 
  virtual void postInputInitialization() override;
 
  virtual void initializePreSubGroups() override;
 
  virtual void initializeWellPostInitialConditionsPreSubGroups( WellElementSubRegion & subRegion ) = 0;
  virtual void printRates( real64 const & time_n,
                           real64 const & dt,
                           WellElementSubRegion const & subRegion ) = 0;
 
 
  template< typename CONSTITUTIVE_BASE_TYPE >
  static string getConstitutiveName( ElementSubRegionBase const & subRegion );
 
  template< typename CONSTITUTIVE >
  void setConstitutiveName( ElementSubRegionBase & subRegion, string const & wrapperName, string const & constitutiveType ) const;
 
  string_array const & getTargetRegionNames() const {return m_targetRegionNames;}
  std::string getFlowSolverName() const { return m_flowSolverName; }
 
  void setFlowSolverName( const std::string & flowSolverName )   { m_flowSolverName = flowSolverName;    }
 
  std::string getDiscretizationName() const { return m_discretizationName; }
 
  void setDiscretizationName( const std::string & discretizationName )   { m_discretizationName = discretizationName;    }
  ConstraintTypeId getControl() const { return m_currentControl; }
 
  void setControl( ConstraintTypeId const & newControl )  {  m_currentControl = newControl; }
 
  ConstraintTypeId getInputControl() const { return m_inputControl; }
 
  integer estimateSolution() const { return m_estimateSolution; }
 
 
  real64 getTargetBHP( real64 const & targetTime, const ConstraintSourceId source = ConstraintSourceId::USER ) const;
 
  real64 getInjectionTemperature() const;
 
  arrayView1d< real64 const > getInjectionStream() const;
 
  integer getConstraintPhaseIndex() const;
 
  real64 getReferenceElevation() const;
 
  integer useSurfaceConditions() const { return m_useSurfaceConditions; }
 
  string const & referenceReservoirRegion() const { return m_referenceReservoirRegion; }
 
  const real64 & getSurfacePressure() const { return m_surfacePres; }
 
  const real64 & getSurfaceTemperature() const { return m_surfaceTemp; }
 
  bool isInjector() const { return ( m_type == Type::INJECTOR ); }
 
  bool isProducer() const { return ( m_type == Type::PRODUCER ); }
 
  integer isThermal() const { return m_isThermal; }
 
  void setThermal( bool isThermal )   {  m_isThermal=isThermal; }
 
  void setUseMass( integer useMass )   {  m_useMass=useMass; }
 
  integer useMass( )   {  return m_useMass;}
 
  bool isWellOpen() const;
 
  void setWellState( bool open );
  bool getWellState() const;
 
  void setCurrentConstraint( WellConstraintBase * currentConstraint )
  {
    setControl( currentConstraint->getControl()  );
    m_currentConstraint = currentConstraint;
  }
  WellConstraintBase * getCurrentConstraint() { return m_currentConstraint; }
  WellConstraintBase const * getCurrentConstraint() const { return m_currentConstraint; }
 
  bool isCrossflowEnabled() const { return m_isCrossflowEnabled; }
 
  real64 getInitialPressureCoefficient() const { return m_initialPressureCoefficient; }
 
  void setNextDtFromTables( real64 const & currentTime, real64 & nextDt );
  void setKeepVariablesConstantDuringInitStep( bool const keepVariablesConstantDuringInitStep )
  { m_keepVariablesConstantDuringInitStep = keepVariablesConstantDuringInitStep; }
 
  void setFluidSeparator( std::unique_ptr< constitutive::ConstitutiveBase > fluidSeparatorPtr )  {  m_fluidSeparatorPtr = std::move( fluidSeparatorPtr );}
  constitutive::MultiFluidBase & getMultiFluidSeparator()  { return dynamicCast< constitutive::MultiFluidBase & >( *m_fluidSeparatorPtr ); }
 
  constitutive::SingleFluidBase & getSingleFluidSeparator()  { return dynamicCast< constitutive::SingleFluidBase & >( *m_fluidSeparatorPtr ); }
 
  real64 getRegionAveragePressure() const { return m_regionAveragePressure; }
 
  void setRegionAveragePressure( real64 regionAveragePressure ) { m_regionAveragePressure = regionAveragePressure; }
 
  real64 getRegionAverageTemperature() const { return m_regionAverageTemperature; }
 
  void setRegionAverageTemperature( real64 regionAverageTemperature ) { m_regionAverageTemperature = regionAverageTemperature; }
 
  void setWellStatus ( real64 const & currentTime, WellControls::Status status );
 
  WellControls::Status getWellStatus () const { return m_wellStatus; }
 
  virtual string wellElementDofName() const = 0;
 
  virtual string resElementDofName() const = 0;
 
  localIndex numDofPerWellElement() const { return m_numDofPerWellElement; }
 
  localIndex numDofPerResElement() const { return m_numDofPerResElement; }
 
 
  virtual localIndex numFluidComponents() const = 0;
 
  virtual localIndex numFluidPhases() const = 0;
 
  struct viewKeyStruct
  {
    static constexpr char const * fluidNamesString() { return "fluidNames"; }
    static constexpr char const * writeCSVFlagString() { return "writeCSV"; }
    static constexpr char const * timeStepFromTablesFlagString() { return "timeStepFromTables"; }
    static constexpr char const * targetRegionsString() { return "targetRegions"; }
 
    static constexpr char const * refElevString() { return "referenceElevation"; }
    static constexpr char const * typeString() { return "type"; }
    static constexpr char const * currentControlString() { return "currentControl"; }
    static constexpr char const * inputControlString() { return "control"; }
    static constexpr char const * useSurfaceConditionsString() { return "useSurfaceConditions"; }
    static constexpr char const * referenceReservoirRegionString() { return "referenceReservoirRegion"; }
    static constexpr char const * surfacePressureString() { return "surfacePressure"; }
    static constexpr char const * surfaceTemperatureString() { return "surfaceTemperature"; }
 
    static constexpr char const * statusTableNameString() { return "statusTableName"; }
    static constexpr char const * perfStatusTableNameString() { return "perfStatusTableName"; }
    static constexpr char const * enableCrossflowString() { return "enableCrossflow"; }
    static constexpr char const * initialPressureCoefficientString() { return "initialPressureCoefficient"; }
 
    static constexpr char const * estimateWellSolutionString() { return "estimateWellSolution"; }
    static constexpr char const * enableIsoThermalEstimatorString() { return "enableIsoThermalEstimator"; }
 
    // control data (not registered on the mesh)
    static constexpr char const * massDensityString() { return "massDensity";}
 
    static constexpr char const * currentBHPString() { return "currentBHP"; }
 
    static constexpr char const * currentPhaseVolRateString() { return "currentPhaseVolumetricRate"; }
    static constexpr char const * currentVolRateString() { return "currentVolRate"; }
 
    static constexpr char const * currentTotalVolRateString() { return "currentTotalVolumetricRate"; }
 
    static constexpr char const * currentMassRateString() { return "currentMassRate"; }
  };
 
  struct groupKeyStruct
  {
    static constexpr char const * wellNewtonSolverString() { return "WellNewtonSolver"; }
  };
 
  void setPerforationStatus( real64 const & time_n, WellElementSubRegion & subRegion );
  void setGravCoef( WellElementSubRegion & subRegion, R1Tensor const & gravVector );
  static void setNextDtFromTable( TableFunction const * table, real64 const currentTime, real64 & nextDt );
 
  template< typename ConstraintType > void createConstraint ( string const & constraintName );
 
 
  WellConstraintBase const * getBHPConstraint( const ConstraintSourceId source = ConstraintSourceId::USER ) const;
  WellConstraintBase * getBHPConstraint( const ConstraintSourceId source = ConstraintSourceId::USER );
 
  stdVector< WellConstraintBase const * > getRateConstraints() const;
  stdVector< WellConstraintBase * > getRateConstraints();
 
  stdVector< WellConstraintBase const * > getAllConstraints() const;
  stdVector< WellConstraintBase * > getAllConstraints();
 
  void enableThermalEffects ( bool enable ) { m_thermalEffectsEnabled = enable; };
 
  bool thermalEffectsEnabled() const { return m_thermalEffectsEnabled; }
 
  bool isoThermalEstimatorEnabled() const { return m_enableIsoThermalEstimator; }
 
  void setupWellDofs( DomainPartition & domain, WellElementRegion & wellElementRegion, string const & meshBodyName, MeshLevel const & meshLevel );
  WellNewtonSolver & getWellNewtonSolver() { return m_wellNewtonSolver; }
 
  void selectWellConstraint( real64 const & time_n,
                             real64 const & dt,
                             integer const cycleNumber,
                             integer const coupledIterationNumber,
                             DomainPartition & domain,
                             string const & meshBodyName,
                             MeshLevel & mesh,
                             ElementRegionManager & elemManager,
                             WellElementSubRegion & subRegion,
                             DofManager const & dofManager );
 
protected:
  virtual void postRestartInitialization( )override;
 
  void logConstraint( WellConstraintBase const * constraint,
                      WellElementSubRegion const & region,
                      real64 time,
                      bool isLimiting = false ) const;
 
  bool validateReferenceRegion() const;
 
  //template< typename STATISTICS >
  //bool validateReferenceRegionStatistics( ElementRegionManager const & elementManager,
  //                                        real64 & averagePressure,
  //                                        real64 & averageTemperature ) const;
 
private:
  string_array m_targetRegionNames;
 
protected:
 
  Type m_type;
 
  std::string m_flowSolverName;
 
  std::string m_discretizationName;
 
  integer m_useMass;
 
  integer m_numPhases;
 
  integer m_numComponents;
 
  integer m_numDofPerWellElement;
 
  integer m_numDofPerResElement;
 
  integer m_isThermal;
  bool m_keepVariablesConstantDuringInitStep;
  integer m_writeCSV;
  string const m_ratesOutputDir;
 
  // flag to enable time step selection base on rates/bhp tables coordinates
  integer m_timeStepFromTables;
 
  real64 m_refElevation;
 
  real64 m_refGravCoef;
 
  ConstraintTypeId m_inputControl;
 
  ConstraintTypeId m_currentControl;
 
  integer m_useSurfaceConditions;
 
  // Fuild model to compute properties for constraint equation user specified conditions
  std::unique_ptr< constitutive::ConstitutiveBase >  m_fluidSeparatorPtr;
 
  string m_referenceFluidModelName;
 
  string m_referenceReservoirRegion;
 
  real64 m_surfacePres;
 
  real64 m_surfaceTemp;
 
  string m_perfStatusTableName;
 
  integer m_isCrossflowEnabled;
 
  real64 m_initialPressureCoefficient;
 
  // Current constrint
  WellConstraintBase * m_currentConstraint{};
 
  WellControls::Status m_wellStatus;
 
  bool m_wellOpen;
 
  string m_statusTableName;
 
  TableFunction const * m_statusTable;
 
  real64 m_regionAveragePressure;
 
  real64 m_regionAverageTemperature;
 
  integer m_estimateSolution;
  integer m_enableIsoThermalEstimator;
  bool m_thermalEffectsEnabled;
 
  WellNewtonSolver m_wellNewtonSolver;
 
 
  DofManager m_estimatorDoFManager;
  bool m_dofManagerInitialized;
};
 
 
// Use local aliases to avoid accidental macro expansion of the tokens 'Type' or 'Control'
using WellControls_Type = WellControls::Type;
ENUM_STRINGS( WellControls_Type,
              "producer",
              "injector" );
 
using WellControls_Control = ConstraintTypeId;
ENUM_STRINGS( WellControls_Control,
              "BHP",
              "phaseVolRate",
              "totalVolRate",
              "massRate",
              "uninitialized" );
 
 
template< typename CONSTITUTIVE >
void WellControls::setConstitutiveName( ElementSubRegionBase & subRegion, string const & wrapperName, string const & constitutiveType ) const
{
  subRegion.registerWrapper< string >( wrapperName ).
    setPlotLevel( dataRepository::PlotLevel::NOPLOT ).
    setRestartFlags( dataRepository::RestartFlags::NO_WRITE ).
    setSizedFromParent( 0 );
 
  string & constitutiveName = subRegion.getReference< string >( wrapperName );
  constitutiveName = getConstitutiveName< CONSTITUTIVE >( subRegion );
  GEOS_ERROR_IF( constitutiveName.empty(), GEOS_FMT( "{}: {} constitutive model not found on subregion {}",
                                                     getDataContext(), constitutiveType, subRegion.getName() ) );
}
template< typename CONSTITUTIVE_BASE_TYPE >
string WellControls::getConstitutiveName( ElementSubRegionBase const & subRegion )
{
  string validName;
  dataRepository::Group const & constitutiveModels = subRegion.getConstitutiveModels();
 
  constitutiveModels.forSubGroups< CONSTITUTIVE_BASE_TYPE >( [&]( dataRepository::Group const & model )
  {
    GEOS_ERROR_IF( !validName.empty(), "A valid constitutive model was already found." );
    validName = model.getName();
  } );
 
  return validName;
}
 
template< typename BASETYPE = constitutive::ConstitutiveBase, typename LOOKUP_TYPE >
static BASETYPE const & getConstitutiveModel( dataRepository::Group const & dataGroup, LOOKUP_TYPE const & key )
{
  dataRepository::Group const & constitutiveModels = dataGroup.getGroup( ElementSubRegionBase::groupKeyStruct::constitutiveModelsString() );
  return constitutiveModels.getGroup< BASETYPE >( key );
}
template< typename BASETYPE = constitutive::ConstitutiveBase, typename LOOKUP_TYPE >
static BASETYPE & getConstitutiveModel( dataRepository::Group & dataGroup, LOOKUP_TYPE const & key )
{
  dataRepository::Group & constitutiveModels = dataGroup.getGroup( ElementSubRegionBase::groupKeyStruct::constitutiveModelsString() );
  return constitutiveModels.getGroup< BASETYPE >( key );
}
} //namespace geos
 
#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_WELLCONTROLS_HPP