GEOS
CoupledReservoirAndWellsBase.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 
21 #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_COUPLEDRESERVOIRANDWELLSBASE_HPP_
22 #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_COUPLEDRESERVOIRANDWELLSBASE_HPP_
23 
25 
26 #include "common/TimingMacros.hpp"
27 #include "constitutive/permeability/PermeabilityFields.hpp"
28 #include "constitutive/permeability/PermeabilityBase.hpp"
30 #include "mesh/DomainPartition.hpp"
31 #include "physicsSolvers/fluidFlow/wells/WellControls.hpp"
33 
34 namespace geos
35 {
36 
37 namespace coupledReservoirAndWellsInternal
38 {
50 void
52  DomainPartition & domain,
53  DofManager & dofManager,
54  arrayView1d< localIndex > const & rowLengths,
55  integer const resNumDof,
56  integer const wellNumDof,
57  string const & resElemDofName,
58  string const & wellElemDofName );
59 
67 bool validateWellPerforations( PhysicsSolverBase const * const reservoirSolver,
68  WellSolverBase const * const wellSolver,
69  DomainPartition const & domain );
70 
71 }
72 
73 template< typename RESERVOIR_SOLVER, typename WELL_SOLVER >
74 class CoupledReservoirAndWellsBase : public CoupledSolver< RESERVOIR_SOLVER, WELL_SOLVER >
75 {
76 public:
77 
79  using Base::m_solvers;
80  using Base::m_names;
81  using Base::m_dofManager;
82  using Base::m_localMatrix;
83  using Base::m_rhs;
84  using Base::m_solution;
85 
86  enum class SolverType : integer
87  {
88  Reservoir = 0,
89  Well = 1
90  };
91 
93  static string coupledSolverAttributePrefix() { return "reservoirAndWells"; }
94 
100  CoupledReservoirAndWellsBase ( const string & name,
101  dataRepository::Group * const parent )
102  : Base( name, parent ),
104  {
105  this->template getWrapper< string >( Base::viewKeyStruct::discretizationString() ).
106  setInputFlag( dataRepository::InputFlags::FALSE );
107  }
108 
116  virtual void setSparsityPattern( DomainPartition & domain,
117  DofManager & dofManager,
118  CRSMatrix< real64, globalIndex > & localMatrix,
119  SparsityPattern< globalIndex > & pattern ) override
120  {
121  // Set the reservoir sparsity pattern without reservoir-well coupling
122  SparsityPattern< globalIndex > patternDiag;
123  reservoirSolver()->setSparsityPattern( domain, dofManager, localMatrix, patternDiag );
124 
125  // Get the original row lengths (diagonal blocks only)
126  array1d< localIndex > rowLengths( patternDiag.numRows());
127  for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow )
128  {
129  rowLengths[localRow] = patternDiag.numNonZeros( localRow );
130  }
131 
132  // Add the number of nonzeros induced by coupling on perforations
133  addCouplingNumNonzeros( domain, dofManager, rowLengths.toView());
134 
135  // Create a new pattern with enough capacity for coupled matrix
136  pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data());
137 
138  // Copy the original nonzeros
139  for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow )
140  {
141  globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous();
142  pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow ));
143  }
144 
145  // Add the nonzeros from coupling
146  addCouplingSparsityPattern( domain, dofManager, pattern.toView());
147  }
148 
155  RESERVOIR_SOLVER *
156  reservoirSolver() const { return std::get< toUnderlying( SolverType::Reservoir ) >( m_solvers ); }
157 
162  WELL_SOLVER *
163  wellSolver() const { return std::get< toUnderlying( SolverType::Well ) >( m_solvers ); }
164 
165  virtual void
167  {
169 
170  DomainPartition & domain = this->template getGroupByPath< DomainPartition >( "/Problem/domain" );
171 
172  // Validate well perforations: Ensure that each perforation is in a region targeted by the solver
173  if( !validateWellPerforations( domain ))
174  {
175  GEOS_ERROR( GEOS_FMT( "{}: well perforations validation failed, bad perforations found", this->getName()));
176  }
177  }
178 
179  virtual void
181  {
183 
184  setMGRStrategy();
185  }
186 
187  virtual void
188  implicitStepSetup( real64 const & time_n,
189  real64 const & dt,
190  DomainPartition & domain ) override
191  {
192  Base::implicitStepSetup( time_n, dt, domain );
193 
194  // we delay the computation of the transmissibility until the last minute
195  // because we want to make sure that the permeability has been updated (in the flow solver)
196  // this is necessary for some permeability models (like Karman-Kozeny) that do not use the imported permeability
197  // ultimately, we may want to use this mechanism to update the well transmissibility at each time step (if needed)
199  {
200  computeWellTransmissibility( domain );
202  }
203  }
204 
205  void initializeState( DomainPartition & domain ) const { return reservoirSolver()->initializeState( domain ); }
206 
207  void
208  assembleFluxTerms( real64 const dt,
209  DomainPartition const & domain,
210  DofManager const & dofManager,
211  CRSMatrixView< real64, globalIndex const > const & localMatrix,
212  arrayView1d< real64 > const & localRhs ) const
213  { reservoirSolver()->assembleFluxTerms( dt, domain, dofManager, localMatrix, localRhs ); }
214 
215  void
216  assembleStabilizedFluxTerms( real64 const dt,
217  DomainPartition const & domain,
218  DofManager const & dofManager,
219  CRSMatrixView< real64, globalIndex const > const & localMatrix,
220  arrayView1d< real64 > const & localRhs ) const
221  { reservoirSolver()->assembleStabilizedFluxTerms( dt, domain, dofManager, localMatrix, localRhs ); }
222 
223  real64 updateFluidState( ElementSubRegionBase & subRegion ) const
224  { return reservoirSolver()->updateFluidState( subRegion ); }
225  void updatePorosityAndPermeability( CellElementSubRegion & subRegion ) const
226  { reservoirSolver()->updatePorosityAndPermeability( subRegion ); }
227  void updateSolidInternalEnergyModel( ObjectManagerBase & dataGroup ) const
228  { reservoirSolver()->updateSolidInternalEnergyModel( dataGroup ); }
229 
230  integer & isThermal() { return reservoirSolver()->isThermal(); }
231 
232  void enableJumpStabilization()
233  { reservoirSolver()->enableJumpStabilization(); }
234 
235  void enableFixedStressPoromechanicsUpdate()
236  { reservoirSolver()->enableFixedStressPoromechanicsUpdate(); }
237 
238  void setKeepVariablesConstantDuringInitStep( bool const keepVariablesConstantDuringInitStep )
239  {
240  reservoirSolver()->setKeepVariablesConstantDuringInitStep( keepVariablesConstantDuringInitStep );
241  wellSolver()->setKeepVariablesConstantDuringInitStep( keepVariablesConstantDuringInitStep );
242  }
243 
244  virtual void saveSequentialIterationState( DomainPartition & domain ) override
245  { reservoirSolver()->saveSequentialIterationState( domain ); }
246 
247 protected:
248 
255  void
257  DofManager & dofManager,
258  arrayView1d< localIndex > const & rowLengths ) const
259  {
260  coupledReservoirAndWellsInternal::
261  addCouplingNumNonzeros( this,
262  domain,
263  dofManager,
264  rowLengths,
265  wellSolver()->numDofPerResElement(),
266  wellSolver()->numDofPerWellElement(),
267  wellSolver()->resElementDofName(),
268  wellSolver()->wellElementDofName() );
269  }
270 
277  virtual void
279  DofManager const & dofManager,
280  SparsityPatternView< globalIndex > const & pattern ) const = 0;
281 
282  virtual void setMGRStrategy()
283  {
285  GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}", this->getName(), this->getCatalogName()));
286  }
287 
290 
291 private:
292 
298  bool validateWellPerforations( DomainPartition const & domain ) const
299  {
300  return coupledReservoirAndWellsInternal::validateWellPerforations( reservoirSolver(), wellSolver(), domain );
301  }
302 
307  void computeWellTransmissibility( DomainPartition & domain ) const
308  {
309  this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&] ( string const &,
310  MeshLevel & meshLevel,
311  string_array const & regionNames )
312  {
313  ElementRegionManager & elemManager = meshLevel.getElemManager();
314 
315  ElementRegionManager::ElementViewAccessor< arrayView2d< real64 > > const elemCenter =
316  elemManager.constructViewAccessor< array2d< real64 >, arrayView2d< real64 > >( ElementSubRegionBase::viewKeyStruct::elementCenterString() );
317 
318  // loop over the wells
319  elemManager.forElementSubRegions< WellElementSubRegion >( regionNames, [&]( localIndex const,
320  WellElementSubRegion & subRegion )
321  {
322  array1d< array1d< arrayView3d< real64 const > > > const permeability =
323  elemManager.constructMaterialFieldAccessor< constitutive::PermeabilityBase,
324  fields::permeability::permeability >();
325 
326  PerforationData & perforationData = *subRegion.getPerforationData();
327  WellControls const & wellControls = wellSolver()->getWellControls( subRegion );
328 
329  // compute the Peaceman index (if not read from XML)
330  perforationData.computeWellTransmissibility( meshLevel, subRegion, permeability );
331 
332  // if the log level is 1, we output the value of the transmissibilities
333  if( wellControls.getLogLevel() >= 2 )
334  {
335  arrayView2d< real64 const > const perfLocation =
336  perforationData.getField< fields::perforation::location >();
337  arrayView1d< real64 const > const perfTrans =
338  perforationData.getField< fields::perforation::wellTransmissibility >();
339 
340  // get the element region, subregion, index
341  arrayView1d< localIndex const > const resElemRegion =
342  perforationData.getField< fields::perforation::reservoirElementRegion >();
343  arrayView1d< localIndex const > const resElemSubRegion =
344  perforationData.getField< fields::perforation::reservoirElementSubRegion >();
345  arrayView1d< localIndex const > const resElemIndex =
346  perforationData.getField< fields::perforation::reservoirElementIndex >();
347 
348  GEOS_UNUSED_VAR( perfLocation ); // unused if geos_error_if is nulld
349  GEOS_UNUSED_VAR( perfTrans ); // unused if geos_error_if is nulld
350  GEOS_UNUSED_VAR( resElemRegion ); // unused if geos_error_if is nulld
351  GEOS_UNUSED_VAR( resElemSubRegion ); // unused if geos_error_if is nulld
352  GEOS_UNUSED_VAR( resElemIndex ); // unused if geos_error_if is nulld
353 
354  forAll< serialPolicy >( perforationData.size(), [=] ( localIndex const iperf )
355  {
356  GEOS_UNUSED_VAR( iperf ); // unused if geos_error_if is nulld
357  GEOS_LOG_RANK( GEOS_FMT( "{}: perforation at ({},{},{}), perforated element center = ({},{},{}), transmissibility = {} [{}]",
358  this->getName(), perfLocation[iperf][0], perfLocation[iperf][1], perfLocation[iperf][2],
359  elemCenter[resElemRegion[iperf]][resElemSubRegion[iperf]][resElemIndex[iperf]][0],
360  elemCenter[resElemRegion[iperf]][resElemSubRegion[iperf]][resElemIndex[iperf]][1],
361  elemCenter[resElemRegion[iperf]][resElemSubRegion[iperf]][resElemIndex[iperf]][2],
362  perfTrans[iperf], getSymbol( units::Transmissibility ) ) );
363  } );
364  }
365  } );
366  } );
367  }
368 
369 };
370 
371 
372 } /* namespace geos */
373 
374 #endif /* GEOS_PHYSICSSOLVERS_MULTIPHYSICS_COUPLEDRESERVOIRANDWELLSBASE_HPP_ */
bool validateWellPerforations(PhysicsSolverBase const *const reservoirSolver, WellSolverBase const *const wellSolver, DomainPartition const &domain)
Validate the well perforations ensuring that each perforation is located in a reservoir region that i...
void addCouplingNumNonzeros(PhysicsSolverBase const *const solver, DomainPartition &domain, DofManager &dofManager, arrayView1d< localIndex > const &rowLengths, integer const resNumDof, integer const wellNumDof, string const &resElemDofName, string const &wellElemDofName)
Utility function for the implementation details of addCouplingNumZeros.
#define GEOS_ERROR(msg)
Raise a hard error and terminate the program.
Definition: Logger.hpp:157
bool m_isWellTransmissibilityComputed
Flag to determine whether the well transmissibility needs to be computed.
virtual void initializePostInitialConditionsPreSubGroups() override
Called by InitializePostInitialConditions() prior to initializing sub-Groups.
virtual void addCouplingSparsityPattern(DomainPartition const &domain, DofManager const &dofManager, SparsityPatternView< globalIndex > const &pattern) const =0
static string coupledSolverAttributePrefix()
String used to form the solverName used to register solvers in CoupledSolver.
CoupledReservoirAndWellsBase(const string &name, dataRepository::Group *const parent)
main constructor for ManagedGroup Objects
WELL_SOLVER * wellSolver() const
accessor for the pointer to the well solver
void addCouplingNumNonzeros(DomainPartition &domain, DofManager &dofManager, arrayView1d< localIndex > const &rowLengths) const
virtual void saveSequentialIterationState(DomainPartition &domain) override
Save the state of the solver for sequential iteration.
virtual void implicitStepSetup(real64 const &time_n, real64 const &dt, DomainPartition &domain) override
function to perform setup for implicit timestep
RESERVOIR_SOLVER * reservoirSolver() const
accessor for the pointer to the reservoir solver
std::array< string, sizeof...(SOLVERS) > m_names
Names of the single-physics solvers.
std::tuple< SOLVERS *... > m_solvers
Pointers of the single-physics solvers.
The DoFManager is responsible for allocating global dofs, constructing sparsity patterns,...
Definition: DofManager.hpp:45
Partition of the decomposed physical domain. It also manages the connexion information to its neighbo...
Group const & getMeshBodies() const
Get the mesh bodies, const version.
Class facilitating the representation of a multi-level discretization of a MeshBody.
Definition: MeshLevel.hpp:42
Base class for all physics solvers.
virtual string getCatalogName() const =0
CRSMatrix< real64, globalIndex > m_localMatrix
Local system matrix and rhs.
DofManager m_dofManager
Data structure to handle degrees of freedom.
ParallelVector m_solution
System solution vector.
ParallelVector m_rhs
System right-hand side vector.
LinearSolverParametersInput m_linearSolverParameters
Linear solver parameters.
string const & getName() const
Get group name.
Definition: Group.hpp:1331
virtual void initializePostInitialConditionsPreSubGroups()
Called by InitializePostInitialConditions() prior to initializing sub-Groups.
Definition: Group.hpp:1554
virtual void implicitStepSetup(real64 const &time_n, real64 const &dt, DomainPartition &domain) override
virtual void setSparsityPattern(DomainPartition &domain, DofManager &dofManager, CRSMatrix< real64, globalIndex > &localMatrix, SparsityPattern< globalIndex > &pattern) override
Set the sparsity pattern of the linear system matrix.
@ FALSE
Not read from input.
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:179
stdVector< string > string_array
A 1-dimensional array of geos::string types.
Definition: DataTypes.hpp:361
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
Definition: DataTypes.hpp:87
LvArray::CRSMatrix< T, COL_INDEX, INDEX_TYPE, LvArray::ChaiBuffer > CRSMatrix
Alias for CRS Matrix class.
Definition: DataTypes.hpp:305
LvArray::SparsityPatternView< COL_INDEX, INDEX_TYPE const, LvArray::ChaiBuffer > SparsityPatternView
Alias for Sparsity pattern View.
Definition: DataTypes.hpp:301
LvArray::SparsityPattern< COL_INDEX, INDEX_TYPE, LvArray::ChaiBuffer > SparsityPattern
Alias for Sparsity pattern class.
Definition: DataTypes.hpp:297
double real64
64-bit floating point type.
Definition: DataTypes.hpp:98
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:84
int integer
Signed integer type.
Definition: DataTypes.hpp:81
Array< T, 1 > array1d
Alias for 1D array.
Definition: DataTypes.hpp:175
@ mgr
Multigrid reduction (Hypre only)
PreconditionerType preconditionerType
Preconditioner type.