GEOS
PoromechanicsSolver.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_POROMECHANICSSOLVER_HPP_
22 #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSSOLVER_HPP_
23 
28 #include "constitutive/solid/CoupledSolidBase.hpp"
29 #include "constitutive/contact/HydraulicApertureBase.hpp"
30 #include "mesh/DomainPartition.hpp"
32 #include "codingUtilities/Utilities.hpp"
33 
34 namespace geos
35 {
36 
37 namespace stabilization
38 {
39 enum class StabilizationType : integer
40 {
41  None,
42  Global,
43  Local,
44 };
45 
46 ENUM_STRINGS( StabilizationType,
47  "None",
48  "Global",
49  "Local" );
50 }
51 
52 
53 template< typename FLOW_SOLVER, typename MECHANICS_SOLVER = SolidMechanicsLagrangianFEM >
54 class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER >
55 {
56 public:
57 
59  using Base::m_solvers;
60  using Base::m_dofManager;
61  using Base::m_localMatrix;
62  using Base::m_rhs;
63  using Base::m_solution;
64 
65  enum class SolverType : integer
66  {
67  Flow = 0,
68  SolidMechanics = 1
69  };
70 
72  static string coupledSolverAttributePrefix() { return "poromechanics"; }
73 
79  PoromechanicsSolver( const string & name,
80  dataRepository::Group * const parent )
81  : Base( name, parent ),
82  m_isThermal( 0 ),
84  {
86  setApplyDefaultValue( 0 ).
88  setDescription( "Flag indicating whether the problem is thermal or not. Set isThermal=\"1\" to enable the thermal coupling" );
89 
92  setDescription( "StabilizationType. Options are:\n" +
93  toString( stabilization::StabilizationType::None ) + "- Add no stabilization to mass equation \n" +
94  toString( stabilization::StabilizationType::Global ) + "- Add jump stabilization to all faces \n" +
95  toString( stabilization::StabilizationType::Local ) + "- Add jump stabilization on interior of macro elements" );
96 
98  setRTTypeName( rtTypes::CustomTypes::groupNameRefArray ).
100  setDescription( "Regions where stabilization is applied." );
101 
103  setApplyDefaultValue( 1.0 ).
104  setInputFlag( dataRepository::InputFlags::OPTIONAL ).
105  setDescription( "Constant multiplier of stabilization strength" );
106  }
107 
108  virtual void postInputInitialization() override
109  {
111 
112  GEOS_THROW_IF( this->m_isThermal && !this->flowSolver()->isThermal(),
113  GEOS_FMT( "{} {}: The attribute `{}` of the flow solver must be thermal since the poromechanics solver is thermal",
114  this->getCatalogName(), this->getName(), this->flowSolver()->getName() ),
115  InputError );
116 
118  GEOS_FMT( "{} {}: The attribute `{}` of solid mechanics solver `{}` must be `{}`",
119  this->getCatalogName(), this->getName(),
120  SolidMechanicsLagrangianFEM::viewKeyStruct::timeIntegrationOptionString(),
121  this->solidMechanicsSolver()->getName(),
123  InputError );
124  }
125 
126  virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override final
127  {
129 
130  this->template setConstitutiveName< constitutive::CoupledSolidBase >( subRegion,
131  viewKeyStruct::porousMaterialNamesString(), "coupled solid" );
132 
133  // This is needed by the way the surface generator currently does things.
134  this->template setConstitutiveName< constitutive::PorosityBase >( subRegion,
135  constitutive::CoupledSolidBase::viewKeyStruct::porosityModelNameString(), "porosity" );
136 
137  if( dynamic_cast< SurfaceElementSubRegion * >( &subRegion ) )
138  {
139  this->template setConstitutiveName< constitutive::HydraulicApertureBase >( subRegion,
140  viewKeyStruct::hydraulicApertureRelationNameString(), "hydraulic aperture" );
141  }
142  }
143 
144  virtual void initializePreSubGroups() override
145  {
147 
148  GEOS_THROW_IF( m_stabilizationType == stabilization::StabilizationType::Local,
150  ": Local stabilization has been temporarily disabled",
151  InputError );
152 
153  DomainPartition & domain = this->template getGroupByPath< DomainPartition >( "/Problem/domain" );
154 
155  this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&] ( string const &,
156  MeshLevel & mesh,
157  string_array const & regionNames )
158  {
159  ElementRegionManager & elementRegionManager = mesh.getElemManager();
160  elementRegionManager.forElementSubRegions< CellElementSubRegion >( regionNames,
161  [&]( localIndex const,
162  ElementSubRegionBase & subRegion )
163  {
164  if( subRegion.hasField< fields::poromechanics::bulkDensity >() )
165  {
166  // get the solid model to know the number of quadrature points and resize the bulk density
167  constitutive::CoupledSolidBase const & solid =
168  this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion,
169  subRegion.getReference< string >( viewKeyStruct::porousMaterialNamesString() ) );
170  subRegion.getField< fields::poromechanics::bulkDensity >().resizeDimension< 1 >( solid.getDensity().size( 1 ) );
171  }
172  } );
173  } );
174  }
175 
176  virtual void registerDataOnMesh( dataRepository::Group & meshBodies ) override
177  {
178  Base::registerDataOnMesh( meshBodies );
179 
180  if( this->getNonlinearSolverParameters().m_couplingType == NonlinearSolverParameters::CouplingType::Sequential )
181  {
182  // to let the solid mechanics solver that there is a pressure and temperature RHS in the mechanics solve
183  solidMechanicsSolver()->enableFixedStressPoromechanicsUpdate();
184  // to let the flow solver that saving pressure_k and temperature_k is necessary (for the fixed-stress porosity terms)
185  flowSolver()->enableFixedStressPoromechanicsUpdate();
186  }
187 
188  if( m_stabilizationType == stabilization::StabilizationType::Global || m_stabilizationType == stabilization::StabilizationType::Local )
189  {
190  flowSolver()->enableJumpStabilization();
191  }
192 
193  this->forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &,
194  MeshLevel & mesh,
195  string_array const & regionNames )
196  {
197  ElementRegionManager & elemManager = mesh.getElemManager();
198 
199  elemManager.forElementSubRegions< CellElementSubRegion >( regionNames,
200  [&]( localIndex const,
201  ElementSubRegionBase & subRegion )
202  {
203  subRegion.registerWrapper< string >( viewKeyStruct::porousMaterialNamesString() ).
204  setPlotLevel( dataRepository::PlotLevel::NOPLOT ).
205  setRestartFlags( dataRepository::RestartFlags::NO_WRITE ).
206  setSizedFromParent( 0 );
207 
208  // This is needed by the way the surface generator currently does things.
209  subRegion.registerWrapper< string >( constitutive::CoupledSolidBase::viewKeyStruct::porosityModelNameString() ).
210  setPlotLevel( dataRepository::PlotLevel::NOPLOT ).
211  setRestartFlags( dataRepository::RestartFlags::NO_WRITE ).
212  setSizedFromParent( 0 );
213 
214  // register the bulk density for use in the solid mechanics solver
215  // ideally we would resize it here as well, but the solid model name is not available yet (see below)
216  subRegion.registerField< fields::poromechanics::bulkDensity >( this->getName() );
217 
218  if( m_stabilizationType == stabilization::StabilizationType::Global || m_stabilizationType == stabilization::StabilizationType::Local )
219  {
220  subRegion.registerField< fields::flow::macroElementIndex >( this->getName());
221  subRegion.registerField< fields::flow::elementStabConstant >( this->getName());
222  }
223  } );
224  } );
225  }
226 
227  virtual void implicitStepSetup( real64 const & time_n,
228  real64 const & dt,
229  DomainPartition & domain ) override
230  {
231  flowSolver()->setKeepVariablesConstantDuringInitStep( m_performStressInitialization );
232 
233  if( this->m_stabilizationType == stabilization::StabilizationType::Global || this->m_stabilizationType == stabilization::StabilizationType::Local )
234  {
235  this->updateStabilizationParameters( domain );
236  }
237 
238  Base::implicitStepSetup( time_n, dt, domain );
239  }
240 
241  virtual void setupDofs( DomainPartition const & domain,
242  DofManager & dofManager ) const override
243  {
244  // note that the order of operations matters a lot here (for instance for the MGR labels)
245  // we must set up dofs for solid mechanics first, and then for flow
246  // that's the reason why this function is here and not in CoupledSolvers.hpp
247  solidMechanicsSolver()->setupDofs( domain, dofManager );
248  flowSolver()->setupDofs( domain, dofManager );
249  this->setupCoupling( domain, dofManager );
250  }
251 
252  virtual bool checkSequentialConvergence( int const & iter,
253  real64 const & time_n,
254  real64 const & dt,
255  DomainPartition & domain ) override
256  {
257  // always force outer loop for initialization
258  auto & subcycling = this->getNonlinearSolverParameters().m_subcyclingOption;
259  auto const subcycling_orig = subcycling;
260  if( m_performStressInitialization )
261  subcycling = 1;
262 
263  bool isConverged = Base::checkSequentialConvergence( iter, time_n, dt, domain );
264 
265  // restore original
266  subcycling = subcycling_orig;
267 
268  return isConverged;
269  }
270 
275  MECHANICS_SOLVER * solidMechanicsSolver() const
276  {
277  return std::get< toUnderlying( SolverType::SolidMechanics ) >( m_solvers );
278  }
279 
284  FLOW_SOLVER * flowSolver() const
285  {
286  return std::get< toUnderlying( SolverType::Flow ) >( m_solvers );
287  }
288 
289  /*
290  * @brief Utility function to set the stress initialization flag
291  * @param[in] performStressInitialization true if the solver has to initialize stress, false otherwise
292  */
293  void setStressInitialization( bool const performStressInitialization )
294  {
295  m_performStressInitialization = performStressInitialization;
296  solidMechanicsSolver()->setStressInitialization( performStressInitialization );
297  }
298 
300  {
302  constexpr static char const * porousMaterialNamesString() { return "porousMaterialNames"; }
303 
305  constexpr static char const * isThermalString() { return "isThermal"; }
306 
308  constexpr static char const * stabilizationTypeString() {return "stabilizationType"; }
309 
311  constexpr static const char * stabilizationRegionNamesString() {return "stabilizationRegionNames"; }
312 
314  constexpr static const char * stabilizationMultiplierString() {return "stabilizationMultiplier"; }
315 
317  static constexpr char const * hydraulicApertureRelationNameString() {return "hydraulicApertureRelationName"; }
318 
319  };
320 
321  void updateStabilizationParameters( DomainPartition & domain ) const
322  {
323  // Step 1: we loop over the regions where stabilization is active and collect their name
324  set< string > regionFilter;
325  for( string const & regionName : m_stabilizationRegionNames )
326  {
327  regionFilter.insert( regionName );
328  }
329 
330  // Step 2: loop over target regions of solver, and tag the elements belonging to the stabilization regions
331  this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&] ( string const &,
332  MeshLevel & mesh,
333  string_array const & targetRegionNames )
334  {
335  //keep only target regions in filter
336  string_array filteredTargetRegionNames;
337  filteredTargetRegionNames.reserve( targetRegionNames.size() );
338 
339  for( string const & targetRegionName : targetRegionNames )
340  {
341  if( regionFilter.count( targetRegionName ) )
342  {
343  filteredTargetRegionNames.emplace_back( targetRegionName );
344  }
345  }
346 
347  // Loop over elements and update stabilization constant
348  mesh.getElemManager().forElementSubRegions( filteredTargetRegionNames, [&]( localIndex const,
349  ElementSubRegionBase & subRegion )
350  {
351  arrayView1d< integer > const macroElementIndex = subRegion.getField< fields::flow::macroElementIndex >();
352  arrayView1d< real64 > const elementStabConstant = subRegion.getField< fields::flow::elementStabConstant >();
353 
354  constitutive::CoupledSolidBase const & porousSolid =
355  this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion,
356  subRegion.getReference< string >( viewKeyStruct::porousMaterialNamesString() ) );
357 
358  arrayView1d< real64 const > const bulkModulus = porousSolid.getBulkModulus();
359  arrayView1d< real64 const > const shearModulus = porousSolid.getShearModulus();
360  arrayView1d< real64 const > const biotCoefficient = porousSolid.getBiotCoefficient();
361 
362  real64 const stabilizationMultiplier = m_stabilizationMultiplier;
363 
364  forAll< parallelDevicePolicy<> >( subRegion.size(), [bulkModulus,
365  shearModulus,
366  biotCoefficient,
367  stabilizationMultiplier,
368  macroElementIndex,
369  elementStabConstant] GEOS_HOST_DEVICE ( localIndex const ei )
370 
371  {
372  real64 const bM = bulkModulus[ei];
373  real64 const sM = shearModulus[ei];
374  real64 const bC = biotCoefficient[ei];
375 
376  macroElementIndex[ei] = 1;
377  elementStabConstant[ei] = stabilizationMultiplier * 9.0 * (bC * bC) / (32.0 * (10.0 * sM / 3.0 + bM));
378  } );
379  } );
380  } );
381 
382  }
383 
384  void updateBulkDensity( DomainPartition & domain )
385  {
386  this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&]( string const &,
387  MeshLevel & mesh,
388  string_array const & regionNames )
389  {
390  mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
391  auto & subRegion )
392  {
393  // update the bulk density
394  // TODO: ideally, we would not recompute the bulk density, but a more general "rhs" containing the body force and the
395  // pressure/temperature terms
396  updateBulkDensity( subRegion );
397  } );
398  } );
399  }
400 
401 protected:
402 
403  template< typename TYPE_LIST,
404  typename KERNEL_WRAPPER,
405  typename ... PARAMS >
406  real64 assemblyLaunch( MeshLevel & mesh,
407  DofManager const & dofManager,
408  string_array const & regionNames,
409  string const & materialNamesString,
410  CRSMatrixView< real64, globalIndex const > const & localMatrix,
411  arrayView1d< real64 > const & localRhs,
412  real64 const dt,
413  PARAMS && ... params )
414  {
416 
417  NodeManager const & nodeManager = mesh.getNodeManager();
418 
419  string const dofKey = dofManager.getKey( fields::solidMechanics::totalDisplacement::key() );
420  arrayView1d< globalIndex const > const & dofNumber = nodeManager.getReference< globalIndex_array >( dofKey );
421 
422  real64 const gravityVectorData[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( this->gravityVector() );
423 
424  KERNEL_WRAPPER kernelWrapper( dofNumber,
425  dofManager.rankOffset(),
426  localMatrix,
427  localRhs,
428  dt,
429  gravityVectorData,
430  std::forward< PARAMS >( params )... );
431 
432  return finiteElement::
433  regionBasedKernelApplication< parallelDevicePolicy< >,
434  TYPE_LIST >( mesh,
435  regionNames,
436  this->solidMechanicsSolver()->getDiscretizationName(),
437  materialNamesString,
438  kernelWrapper );
439  }
440 
441  /* Implementation of Nonlinear Acceleration (Aitken) of averageMeanTotalStressIncrement */
442 
443  void recordAverageMeanTotalStressIncrement( DomainPartition & domain,
444  array1d< real64 > & averageMeanTotalStressIncrement )
445  {
446  averageMeanTotalStressIncrement.resize( 0 );
447  this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
448  MeshLevel & mesh,
449  string_array const & regionNames )
450  {
451  mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
452  auto & subRegion )
453  {
454  // get the solid model (to access stress increment)
455  string const & solidName =
456  subRegion.template getReference< string >( viewKeyStruct::porousMaterialNamesString() );
457  constitutive::CoupledSolidBase & solid =
458  this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, solidName );
459 
460  arrayView1d< const real64 > const & averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
461  for( localIndex k = 0; k < localIndex( averageMeanTotalStressIncrement_k.size()); k++ )
462  {
463  averageMeanTotalStressIncrement.emplace_back( averageMeanTotalStressIncrement_k[k] );
464  }
465  } );
466  } );
467  }
468 
469  void applyAcceleratedAverageMeanTotalStressIncrement( DomainPartition & domain,
470  array1d< real64 > & averageMeanTotalStressIncrement )
471  {
472  integer i = 0;
473  this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
474  MeshLevel & mesh,
475  string_array const & regionNames )
476  {
477  mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
478  auto & subRegion )
479  {
480  // get the solid model (to access stress increment)
481  string const & solidName = subRegion.template getReference< string >( viewKeyStruct::porousMaterialNamesString() );
482  constitutive::CoupledSolidBase & solid =
483  this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, solidName );
484  auto & porosityModel = dynamic_cast< constitutive::BiotPorosity const & >( solid.getBasePorosityModel());
485  arrayView1d< real64 > const & averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
486  for( localIndex k = 0; k < localIndex( averageMeanTotalStressIncrement_k.size()); k++ )
487  {
488  porosityModel.updateAverageMeanTotalStressIncrement( k, averageMeanTotalStressIncrement[i] );
489  i++;
490  }
491  } );
492  } );
493  }
494 
495  real64 computeAitkenRelaxationFactor( array1d< real64 > const & s0,
496  array1d< real64 > const & s1,
497  array1d< real64 > const & s1_tilde,
498  array1d< real64 > const & s2_tilde,
499  real64 const omega0 )
500  {
501  array1d< real64 > r1 = axpy( s1_tilde, s0, -1.0 );
502  array1d< real64 > r2 = axpy( s2_tilde, s1, -1.0 );
503 
504  // diff = r2 - r1
505  array1d< real64 > diff = axpy( r2, r1, -1.0 );
506 
507  real64 const denom = dot( diff, diff );
508  real64 const numer = dot( r1, diff );
509 
510  real64 omega1 = 1.0;
511  if( !isZero( denom ))
512  {
513  omega1 = -1.0 * omega0 * numer / denom;
514  }
515  return omega1;
516  }
517 
518  array1d< real64 > computeUpdate( array1d< real64 > const & s1,
519  array1d< real64 > const & s2_tilde,
520  real64 const omega1 )
521  {
522  return axpy( scale( s1, 1.0 - omega1 ),
523  scale( s2_tilde, omega1 ),
524  1.0 );
525  }
526 
527  void startSequentialIteration( integer const & iter,
528  DomainPartition & domain ) override
529  {
530  if( this->getNonlinearSolverParameters().m_nonlinearAccelerationType == NonlinearSolverParameters::NonlinearAccelerationType::Aitken )
531  {
532  if( iter == 0 )
533  {
534  recordAverageMeanTotalStressIncrement( domain, m_s1 );
535  }
536  else
537  {
538  m_s0 = m_s1;
539  m_s1 = m_s2;
540  m_s1_tilde = m_s2_tilde;
541  m_omega0 = m_omega1;
542  }
543  }
544  }
545 
546  void finishSequentialIteration( integer const & iter,
547  DomainPartition & domain ) override
548  {
549  if( this->getNonlinearSolverParameters().m_nonlinearAccelerationType == NonlinearSolverParameters::NonlinearAccelerationType::Aitken )
550  {
551  if( iter == 0 )
552  {
553  m_s2 = m_s2_tilde;
554  m_omega1 = 1.0;
555  }
556  else
557  {
558  m_omega1 = computeAitkenRelaxationFactor( m_s0, m_s1, m_s1_tilde, m_s2_tilde, m_omega0 );
559  m_s2 = computeUpdate( m_s1, m_s2_tilde, m_omega1 );
560  applyAcceleratedAverageMeanTotalStressIncrement( domain, m_s2 );
561  }
562  }
563  }
564 
565  virtual void mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) override
566  {
568 
570  if( solverType == static_cast< integer >( SolverType::Flow ) )
571  {
572  updateBulkDensity( domain );
573  }
574 
576  if( solverType == static_cast< integer >( SolverType::SolidMechanics )
577  && !m_performStressInitialization ) // do not update during poromechanics initialization
578  {
579  // compute the average of the mean total stress increment over quadrature points
580  averageMeanTotalStressIncrement( domain );
581 
582  this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&]( string const &,
583  MeshLevel & mesh,
584  string_array const & regionNames )
585  {
586 
587  mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
588  auto & subRegion )
589  {
590  // update the porosity after a change in displacement (after mechanics solve)
591  // or a change in pressure/temperature (after a flow solve)
592  flowSolver()->updatePorosityAndPermeability( subRegion );
593  // update bulk density to reflect porosity change into mechanics
594  updateBulkDensity( subRegion );
595  } );
596  } );
597  }
598 
599  // needed to perform nonlinear acceleration
600  if( solverType == static_cast< integer >( SolverType::SolidMechanics ) &&
601  this->getNonlinearSolverParameters().m_nonlinearAccelerationType== NonlinearSolverParameters::NonlinearAccelerationType::Aitken )
602  {
603  recordAverageMeanTotalStressIncrement( domain, m_s2_tilde );
604  }
605  }
606 
612  {
613  this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&]( string const &,
614  MeshLevel & mesh,
615  string_array const & regionNames )
616  {
617  mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
618  auto & subRegion )
619  {
620  // get the solid model (to access stress increment)
621  string const & solidName = subRegion.template getReference< string >( viewKeyStruct::porousMaterialNamesString() );
622  constitutive::CoupledSolidBase & solid =
623  this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, solidName );
624 
625  arrayView2d< real64 const > const meanTotalStressIncrement_k = solid.getMeanTotalStressIncrement_k();
626  arrayView1d< real64 > const averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
627 
628  finiteElement::FiniteElementBase & subRegionFE =
629  subRegion.template getReference< finiteElement::FiniteElementBase >( solidMechanicsSolver()->getDiscretizationName() );
630 
631  // determine the finite element type
632  finiteElement::FiniteElementDispatchHandler< BASE_FE_TYPES >::
633  dispatch3D( subRegionFE, [&] ( auto const finiteElement )
634  {
635  using FE_TYPE = decltype( finiteElement );
636 
637  // call the factory and launch the kernel
638  AverageOverQuadraturePoints1DKernelFactory::
639  createAndLaunch< CellElementSubRegion,
640  FE_TYPE,
641  parallelDevicePolicy<> >( mesh.getNodeManager(),
642  mesh.getEdgeManager(),
643  mesh.getFaceManager(),
644  subRegion,
645  finiteElement,
646  meanTotalStressIncrement_k,
647  averageMeanTotalStressIncrement_k );
648  } );
649  } );
650  } );
651  }
652 
653  virtual void updateBulkDensity( ElementSubRegionBase & subRegion ) = 0;
654 
655  virtual void validateNonlinearAcceleration() override
656  {
657  if( MpiWrapper::commSize( MPI_COMM_GEOS ) > 1 )
658  {
659  GEOS_ERROR( "Nonlinear acceleration is not implemented for MPI runs" );
660  }
661  }
662 
665 
668 
670  stabilization::StabilizationType m_stabilizationType;
671 
674 
677 
679  array1d< real64 > m_s0; // Accelerated averageMeanTotalStresIncrement @ outer iteration v ( two iterations ago )
680  array1d< real64 > m_s1; // Accelerated averageMeanTotalStresIncrement @ outer iteration v + 1 ( previous iteration )
681  array1d< real64 > m_s1_tilde; // Unaccelerated averageMeanTotalStresIncrement @ outer iteration v + 1 ( previous iteration )
682  array1d< real64 > m_s2; // Accelerated averageMeanTotalStresIncrement @ outer iteration v + 2 ( current iteration )
683  array1d< real64 > m_s2_tilde; // Unaccelerated averageMeanTotalStresIncrement @ outer iteration v + 1 ( current iteration )
684  real64 m_omega0; // Old Aitken relaxation factor
685  real64 m_omega1; // New Aitken relaxation factor
686 
687 };
688 
689 } /* namespace geos */
690 
691 #endif //GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSSOLVER_HPP_
#define GEOS_ERROR(msg)
Raise a hard error and terminate the program.
Definition: Logger.hpp:157
#define GEOS_THROW_IF(EXP, msg, TYPE)
Conditionally throw an exception.
Definition: Logger.hpp:151
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
virtual void postInputInitialization() override
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.
The ElementRegionManager class provides an interface to ObjectManagerBase in order to manage ElementR...
void forElementSubRegions(LAMBDA &&lambda)
This function is used to launch kernel function over the element subregions of all the subregion type...
Class facilitating the representation of a multi-level discretization of a MeshBody.
Definition: MeshLevel.hpp:42
ElementRegionManager const & getElemManager() const
Get the element region manager.
Definition: MeshLevel.hpp:207
virtual string getCatalogName() const =0
virtual void setConstitutiveNamesCallSuper(ElementSubRegionBase &subRegion) const
This function sets constitutive name fields on an ElementSubRegionBase, and calls the base function i...
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.
array1d< real64 > m_s0
Member variables needed for Nonlinear Acceleration ( Aitken ). Naming convention follows ( Jiang & Tc...
virtual void implicitStepSetup(real64 const &time_n, real64 const &dt, DomainPartition &domain) override
function to perform setup for implicit timestep
real64 m_stabilizationMultiplier
Stabilization Multiplier.
virtual void initializePreSubGroups() override
Called by Initialize() prior to initializing sub-Groups.
stabilization::StabilizationType m_stabilizationType
Type of stabilization used.
virtual void postInputInitialization() override
FLOW_SOLVER * flowSolver() const
accessor for the pointer to the flow solver
PoromechanicsSolver(const string &name, dataRepository::Group *const parent)
main constructor for CoupledSolver Objects
virtual void setConstitutiveNamesCallSuper(ElementSubRegionBase &subRegion) const override final
This function sets constitutive name fields on an ElementSubRegionBase, and calls the base function i...
MECHANICS_SOLVER * solidMechanicsSolver() const
accessor for the pointer to the solid mechanics solver
virtual void mapSolutionBetweenSolvers(DomainPartition &domain, integer const solverType) override
Maps the solution obtained from one solver to the fields used by the other solver(s)
integer m_isThermal
Flag to determine whether or not this is a thermal simulation.
bool m_performStressInitialization
Flag to indicate that the solver is going to perform stress initialization.
string_array m_stabilizationRegionNames
Names of regions where stabilization applied.
virtual void setupDofs(DomainPartition const &domain, DofManager &dofManager) const override
Populate degree-of-freedom manager with fields relevant to this solver.
virtual void registerDataOnMesh(dataRepository::Group &meshBodies) override
Register wrappers that contain data on the mesh objects.
void averageMeanTotalStressIncrement(DomainPartition &domain)
Helper function to average the mean total stress increment over quadrature points.
static string coupledSolverAttributePrefix()
String used to form the solverName used to register solvers in CoupledSolver.
virtual void initializePreSubGroups()
Called by Initialize() prior to initializing sub-Groups.
Definition: Group.hpp:1542
Wrapper< TBASE > & registerWrapper(string const &name, wrapperMap::KeyIndex::index_type *const rkey=nullptr)
Create and register a Wrapper around a new object.
string const & getName() const
Get group name.
Definition: Group.hpp:1331
DataContext const & getWrapperDataContext(KEY key) const
Definition: Group.hpp:1356
@ NOPLOT
Do not ever write to plot file.
@ OPTIONAL
Optional in input.
@ NO_WRITE
Do not write into restart.
stdVector< string > string_array
A 1-dimensional array of geos::string types.
Definition: DataTypes.hpp:361
int MPI_COMM_GEOS
Global MPI communicator used by GEOSX.
array1d< globalIndex > globalIndex_array
A 1-dimensional array of geos::globalIndex types.
Definition: DataTypes.hpp:370
std::set< T > set
A set of local indices.
Definition: DataTypes.hpp:262
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
ENUM_STRINGS(LinearSolverParameters::SolverType, "direct", "cg", "gmres", "fgmres", "bicgstab", "richardson", "preconditioner")
Declare strings associated with enumeration values.
Provides enum <-> string conversion facilities.
Exception class used to report errors in user input.
Definition: Logger.hpp:464
constexpr static char const * porousMaterialNamesString()
Names of the porous materials.
constexpr static const char * stabilizationRegionNamesString()
Name of regions on which stabilization applied.
static constexpr char const * hydraulicApertureRelationNameString()
Name of the hydraulicApertureRelationName.
constexpr static char const * isThermalString()
Flag to indicate that the simulation is thermal.
constexpr static char const * stabilizationTypeString()
Type of pressure stabilization.
constexpr static const char * stabilizationMultiplierString()
Multiplier on stabilization strength.
Structure to hold scoped key names.
Definition: Group.hpp:1444