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 {
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 
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  }
117 
118  virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override final
119  {
121 
122  this->template setConstitutiveName< constitutive::CoupledSolidBase >( subRegion,
123  viewKeyStruct::porousMaterialNamesString(), "coupled solid" );
124 
125  // This is needed by the way the surface generator currently does things.
126  this->template setConstitutiveName< constitutive::PorosityBase >( subRegion,
127  constitutive::CoupledSolidBase::viewKeyStruct::porosityModelNameString(), "porosity" );
128 
129  if( dynamic_cast< SurfaceElementSubRegion * >( &subRegion ) )
130  {
131  this->template setConstitutiveName< constitutive::HydraulicApertureBase >( subRegion,
132  viewKeyStruct::hydraulicApertureRelationNameString(), "hydraulic aperture" );
133  }
134  }
135 
136  virtual void initializePreSubGroups() override
137  {
139 
140  GEOS_THROW_IF( m_stabilizationType == stabilization::StabilizationType::Local,
142  ": Local stabilization has been temporarily disabled",
143  InputError );
144 
145  DomainPartition & domain = this->template getGroupByPath< DomainPartition >( "/Problem/domain" );
146 
147  this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&] ( string const &,
148  MeshLevel & mesh,
149  string_array const & regionNames )
150  {
151  ElementRegionManager & elementRegionManager = mesh.getElemManager();
152  elementRegionManager.forElementSubRegions< CellElementSubRegion >( regionNames,
153  [&]( localIndex const,
154  ElementSubRegionBase & subRegion )
155  {
156  if( subRegion.hasField< fields::poromechanics::bulkDensity >() )
157  {
158  // get the solid model to know the number of quadrature points and resize the bulk density
159  constitutive::CoupledSolidBase const & solid =
160  this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion,
161  subRegion.getReference< string >( viewKeyStruct::porousMaterialNamesString() ) );
162  subRegion.getField< fields::poromechanics::bulkDensity >().resizeDimension< 1 >( solid.getDensity().size( 1 ) );
163  }
164  } );
165  } );
166  }
167 
168  virtual void registerDataOnMesh( dataRepository::Group & meshBodies ) override
169  {
170  Base::registerDataOnMesh( meshBodies );
171 
172  if( this->getNonlinearSolverParameters().m_couplingType == NonlinearSolverParameters::CouplingType::Sequential )
173  {
174  // to let the solid mechanics solver that there is a pressure and temperature RHS in the mechanics solve
175  solidMechanicsSolver()->enableFixedStressPoromechanicsUpdate();
176  // to let the flow solver that saving pressure_k and temperature_k is necessary (for the fixed-stress porosity terms)
177  flowSolver()->enableFixedStressPoromechanicsUpdate();
178  }
179 
180  if( m_stabilizationType == stabilization::StabilizationType::Global || m_stabilizationType == stabilization::StabilizationType::Local )
181  {
182  flowSolver()->enableJumpStabilization();
183  }
184 
185  this->template forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &,
186  MeshLevel & mesh,
187  string_array const & regionNames )
188  {
189  ElementRegionManager & elemManager = mesh.getElemManager();
190 
191  elemManager.forElementSubRegions< CellElementSubRegion >( regionNames,
192  [&]( localIndex const,
193  ElementSubRegionBase & subRegion )
194  {
195  subRegion.registerWrapper< string >( viewKeyStruct::porousMaterialNamesString() ).
196  setPlotLevel( dataRepository::PlotLevel::NOPLOT ).
197  setRestartFlags( dataRepository::RestartFlags::NO_WRITE ).
198  setSizedFromParent( 0 );
199 
200  // This is needed by the way the surface generator currently does things.
201  subRegion.registerWrapper< string >( constitutive::CoupledSolidBase::viewKeyStruct::porosityModelNameString() ).
202  setPlotLevel( dataRepository::PlotLevel::NOPLOT ).
203  setRestartFlags( dataRepository::RestartFlags::NO_WRITE ).
204  setSizedFromParent( 0 );
205 
206  // register the bulk density for use in the solid mechanics solver
207  // ideally we would resize it here as well, but the solid model name is not available yet (see below)
208  subRegion.registerField< fields::poromechanics::bulkDensity >( this->getName() );
209 
210  if( m_stabilizationType == stabilization::StabilizationType::Global || m_stabilizationType == stabilization::StabilizationType::Local )
211  {
212  subRegion.registerField< fields::flow::macroElementIndex >( this->getName());
213  subRegion.registerField< fields::flow::elementStabConstant >( this->getName());
214  }
215  } );
216  } );
217  }
218 
219  virtual void implicitStepSetup( real64 const & time_n,
220  real64 const & dt,
221  DomainPartition & domain ) override
222  {
223  flowSolver()->setKeepVariablesConstantDuringInitStep( m_performStressInitialization );
224 
225  if( this->m_stabilizationType == stabilization::StabilizationType::Global || this->m_stabilizationType == stabilization::StabilizationType::Local )
226  {
227  this->updateStabilizationParameters( domain );
228  }
229 
230  Base::implicitStepSetup( time_n, dt, domain );
231  }
232 
233  virtual void setupDofs( DomainPartition const & domain,
234  DofManager & dofManager ) const override
235  {
236  // note that the order of operations matters a lot here (for instance for the MGR labels)
237  // we must set up dofs for solid mechanics first, and then for flow
238  // that's the reason why this function is here and not in CoupledSolvers.hpp
239  solidMechanicsSolver()->setupDofs( domain, dofManager );
240  flowSolver()->setupDofs( domain, dofManager );
241  this->setupCoupling( domain, dofManager );
242  }
243 
244  virtual bool checkSequentialConvergence( int const & iter,
245  real64 const & time_n,
246  real64 const & dt,
247  DomainPartition & domain ) override
248  {
249  // always force outer loop for initialization
250  auto & subcycling = this->getNonlinearSolverParameters().m_subcyclingOption;
251  auto const subcycling_orig = subcycling;
252  if( m_performStressInitialization )
253  subcycling = 1;
254 
255  bool isConverged = Base::checkSequentialConvergence( iter, time_n, dt, domain );
256 
257  // restore original
258  subcycling = subcycling_orig;
259 
260  return isConverged;
261  }
262 
267  MECHANICS_SOLVER * solidMechanicsSolver() const
268  {
269  return std::get< toUnderlying( SolverType::SolidMechanics ) >( m_solvers );
270  }
271 
276  FLOW_SOLVER * flowSolver() const
277  {
278  return std::get< toUnderlying( SolverType::Flow ) >( m_solvers );
279  }
280 
281  /*
282  * @brief Utility function to set the stress initialization flag
283  * @param[in] performStressInitialization true if the solver has to initialize stress, false otherwise
284  */
285  void setStressInitialization( bool const performStressInitialization )
286  {
287  m_performStressInitialization = performStressInitialization;
288  solidMechanicsSolver()->setStressInitialization( performStressInitialization );
289  }
290 
292  {
294  constexpr static char const * porousMaterialNamesString() { return "porousMaterialNames"; }
295 
297  constexpr static char const * isThermalString() { return "isThermal"; }
298 
300  constexpr static char const * stabilizationTypeString() {return "stabilizationType"; }
301 
303  constexpr static const char * stabilizationRegionNamesString() {return "stabilizationRegionNames"; }
304 
306  constexpr static const char * stabilizationMultiplierString() {return "stabilizationMultiplier"; }
307 
309  static constexpr char const * hydraulicApertureRelationNameString() {return "hydraulicApertureRelationName"; }
310 
311  };
312 
313  void updateStabilizationParameters( DomainPartition & domain ) const
314  {
315  // Step 1: we loop over the regions where stabilization is active and collect their name
316  set< string > regionFilter;
317  for( string const & regionName : m_stabilizationRegionNames )
318  {
319  regionFilter.insert( regionName );
320  }
321 
322  // Step 2: loop over target regions of solver, and tag the elements belonging to the stabilization regions
323  this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&] ( string const &,
324  MeshLevel & mesh,
325  string_array const & targetRegionNames )
326  {
327  //keep only target regions in filter
328  string_array filteredTargetRegionNames;
329  filteredTargetRegionNames.reserve( targetRegionNames.size() );
330 
331  for( string const & targetRegionName : targetRegionNames )
332  {
333  if( regionFilter.count( targetRegionName ) )
334  {
335  filteredTargetRegionNames.emplace_back( targetRegionName );
336  }
337  }
338 
339  // Loop over elements and update stabilization constant
340  mesh.getElemManager().forElementSubRegions( filteredTargetRegionNames, [&]( localIndex const,
341  ElementSubRegionBase & subRegion )
342  {
343  arrayView1d< integer > const macroElementIndex = subRegion.getField< fields::flow::macroElementIndex >();
344  arrayView1d< real64 > const elementStabConstant = subRegion.getField< fields::flow::elementStabConstant >();
345 
346  constitutive::CoupledSolidBase const & porousSolid =
347  this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion,
348  subRegion.getReference< string >( viewKeyStruct::porousMaterialNamesString() ) );
349 
350  arrayView1d< real64 const > const bulkModulus = porousSolid.getBulkModulus();
351  arrayView1d< real64 const > const shearModulus = porousSolid.getShearModulus();
352  arrayView1d< real64 const > const biotCoefficient = porousSolid.getBiotCoefficient();
353 
354  real64 const stabilizationMultiplier = m_stabilizationMultiplier;
355 
356  forAll< parallelDevicePolicy<> >( subRegion.size(), [bulkModulus,
357  shearModulus,
358  biotCoefficient,
359  stabilizationMultiplier,
360  macroElementIndex,
361  elementStabConstant] GEOS_HOST_DEVICE ( localIndex const ei )
362 
363  {
364  real64 const bM = bulkModulus[ei];
365  real64 const sM = shearModulus[ei];
366  real64 const bC = biotCoefficient[ei];
367 
368  macroElementIndex[ei] = 1;
369  elementStabConstant[ei] = stabilizationMultiplier * 9.0 * (bC * bC) / (32.0 * (10.0 * sM / 3.0 + bM));
370  } );
371  } );
372  } );
373 
374  }
375 
376  void updateBulkDensity( DomainPartition & domain )
377  {
378  this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&]( string const &,
379  MeshLevel & mesh,
380  string_array const & regionNames )
381  {
382  mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
383  auto & subRegion )
384  {
385  // update the bulk density
386  // TODO: ideally, we would not recompute the bulk density, but a more general "rhs" containing the body force and the
387  // pressure/temperature terms
388  updateBulkDensity( subRegion );
389  } );
390  } );
391  }
392 
393 protected:
394 
395  template< typename TYPE_LIST,
396  typename KERNEL_WRAPPER,
397  typename ... PARAMS >
398  real64 assemblyLaunch( MeshLevel & mesh,
399  DofManager const & dofManager,
400  string_array const & regionNames,
401  string const & materialNamesString,
402  CRSMatrixView< real64, globalIndex const > const & localMatrix,
403  arrayView1d< real64 > const & localRhs,
404  real64 const dt,
405  PARAMS && ... params )
406  {
408 
409  NodeManager const & nodeManager = mesh.getNodeManager();
410 
411  string const dofKey = dofManager.getKey( fields::solidMechanics::totalDisplacement::key() );
412  arrayView1d< globalIndex const > const & dofNumber = nodeManager.getReference< globalIndex_array >( dofKey );
413 
414  real64 const gravityVectorData[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( this->gravityVector() );
415 
416  KERNEL_WRAPPER kernelWrapper( dofNumber,
417  dofManager.rankOffset(),
418  localMatrix,
419  localRhs,
420  dt,
421  gravityVectorData,
422  std::forward< PARAMS >( params )... );
423 
424  return finiteElement::
425  regionBasedKernelApplication< parallelDevicePolicy< >,
426  TYPE_LIST >( mesh,
427  regionNames,
428  this->solidMechanicsSolver()->getDiscretizationName(),
429  materialNamesString,
430  kernelWrapper );
431  }
432 
433  /* Implementation of Nonlinear Acceleration (Aitken) of averageMeanTotalStressIncrement */
434 
435  void recordAverageMeanTotalStressIncrement( DomainPartition & domain,
436  array1d< real64 > & averageMeanTotalStressIncrement )
437  {
438  averageMeanTotalStressIncrement.resize( 0 );
439  this->template forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
440  MeshLevel & mesh,
441  string_array const & regionNames )
442  {
443  mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
444  auto & subRegion )
445  {
446  // get the solid model (to access stress increment)
447  string const & solidName =
448  subRegion.template getReference< string >( viewKeyStruct::porousMaterialNamesString() );
449  constitutive::CoupledSolidBase & solid =
450  this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, solidName );
451 
452  arrayView1d< const real64 > const & averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
453  for( localIndex k = 0; k < localIndex( averageMeanTotalStressIncrement_k.size()); k++ )
454  {
455  averageMeanTotalStressIncrement.emplace_back( averageMeanTotalStressIncrement_k[k] );
456  }
457  } );
458  } );
459  }
460 
461  void applyAcceleratedAverageMeanTotalStressIncrement( DomainPartition & domain,
462  array1d< real64 > & averageMeanTotalStressIncrement )
463  {
464  integer i = 0;
465  this->template forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
466  MeshLevel & mesh,
467  string_array const & regionNames )
468  {
469  mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
470  auto & subRegion )
471  {
472  // get the solid model (to access stress increment)
473  string const & solidName = subRegion.template getReference< string >( viewKeyStruct::porousMaterialNamesString() );
474  constitutive::CoupledSolidBase & solid =
475  this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, solidName );
476  auto & porosityModel = dynamic_cast< constitutive::BiotPorosity const & >( solid.getBasePorosityModel());
477  arrayView1d< real64 > const & averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
478  for( localIndex k = 0; k < localIndex( averageMeanTotalStressIncrement_k.size()); k++ )
479  {
480  porosityModel.updateAverageMeanTotalStressIncrement( k, averageMeanTotalStressIncrement[i] );
481  i++;
482  }
483  } );
484  } );
485  }
486 
487  real64 computeAitkenRelaxationFactor( array1d< real64 > const & s0,
488  array1d< real64 > const & s1,
489  array1d< real64 > const & s1_tilde,
490  array1d< real64 > const & s2_tilde,
491  real64 const omega0 )
492  {
493  array1d< real64 > r1 = axpy( s1_tilde, s0, -1.0 );
494  array1d< real64 > r2 = axpy( s2_tilde, s1, -1.0 );
495 
496  // diff = r2 - r1
497  array1d< real64 > diff = axpy( r2, r1, -1.0 );
498 
499  real64 const denom = dot( diff, diff );
500  real64 const numer = dot( r1, diff );
501 
502  real64 omega1 = 1.0;
503  if( !isZero( denom ))
504  {
505  omega1 = -1.0 * omega0 * numer / denom;
506  }
507  return omega1;
508  }
509 
510  array1d< real64 > computeUpdate( array1d< real64 > const & s1,
511  array1d< real64 > const & s2_tilde,
512  real64 const omega1 )
513  {
514  return axpy( scale( s1, 1.0 - omega1 ),
515  scale( s2_tilde, omega1 ),
516  1.0 );
517  }
518 
519  void startSequentialIteration( integer const & iter,
520  DomainPartition & domain ) override
521  {
522  if( this->getNonlinearSolverParameters().m_nonlinearAccelerationType == NonlinearSolverParameters::NonlinearAccelerationType::Aitken )
523  {
524  if( iter == 0 )
525  {
526  recordAverageMeanTotalStressIncrement( domain, m_s1 );
527  }
528  else
529  {
530  m_s0 = m_s1;
531  m_s1 = m_s2;
532  m_s1_tilde = m_s2_tilde;
533  m_omega0 = m_omega1;
534  }
535  }
536  }
537 
538  void finishSequentialIteration( integer const & iter,
539  DomainPartition & domain ) override
540  {
541  if( this->getNonlinearSolverParameters().m_nonlinearAccelerationType == NonlinearSolverParameters::NonlinearAccelerationType::Aitken )
542  {
543  if( iter == 0 )
544  {
545  m_s2 = m_s2_tilde;
546  m_omega1 = 1.0;
547  }
548  else
549  {
550  m_omega1 = computeAitkenRelaxationFactor( m_s0, m_s1, m_s1_tilde, m_s2_tilde, m_omega0 );
551  m_s2 = computeUpdate( m_s1, m_s2_tilde, m_omega1 );
552  applyAcceleratedAverageMeanTotalStressIncrement( domain, m_s2 );
553  }
554  }
555  }
556 
557  virtual void mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) override
558  {
560 
562  if( solverType == static_cast< integer >( SolverType::Flow ) )
563  {
564  updateBulkDensity( domain );
565  }
566 
568  if( solverType == static_cast< integer >( SolverType::SolidMechanics )
569  && !m_performStressInitialization ) // do not update during poromechanics initialization
570  {
571  // compute the average of the mean total stress increment over quadrature points
572  averageMeanTotalStressIncrement( domain );
573 
574  this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&]( string const &,
575  MeshLevel & mesh,
576  string_array const & regionNames )
577  {
578 
579  mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
580  auto & subRegion )
581  {
582  // update the porosity after a change in displacement (after mechanics solve)
583  // or a change in pressure/temperature (after a flow solve)
584  flowSolver()->updatePorosityAndPermeability( subRegion );
585  // update bulk density to reflect porosity change into mechanics
586  updateBulkDensity( subRegion );
587  } );
588  } );
589  }
590 
591  // needed to perform nonlinear acceleration
592  if( solverType == static_cast< integer >( SolverType::SolidMechanics ) &&
593  this->getNonlinearSolverParameters().m_nonlinearAccelerationType== NonlinearSolverParameters::NonlinearAccelerationType::Aitken )
594  {
595  recordAverageMeanTotalStressIncrement( domain, m_s2_tilde );
596  }
597  }
598 
604  {
605  this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&]( string const &,
606  MeshLevel & mesh,
607  string_array const & regionNames )
608  {
609  mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
610  auto & subRegion )
611  {
612  // get the solid model (to access stress increment)
613  string const & solidName = subRegion.template getReference< string >( viewKeyStruct::porousMaterialNamesString() );
614  constitutive::CoupledSolidBase & solid =
615  this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, solidName );
616 
617  arrayView2d< real64 const > const meanTotalStressIncrement_k = solid.getMeanTotalStressIncrement_k();
618  arrayView1d< real64 > const averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
619 
620  finiteElement::FiniteElementBase & subRegionFE =
621  subRegion.template getReference< finiteElement::FiniteElementBase >( solidMechanicsSolver()->getDiscretizationName() );
622 
623  // determine the finite element type
624  finiteElement::FiniteElementDispatchHandler< BASE_FE_TYPES >::
625  dispatch3D( subRegionFE, [&] ( auto const finiteElement )
626  {
627  using FE_TYPE = decltype( finiteElement );
628 
629  // call the factory and launch the kernel
630  AverageOverQuadraturePoints1DKernelFactory::
631  createAndLaunch< CellElementSubRegion,
632  FE_TYPE,
633  parallelDevicePolicy<> >( mesh.getNodeManager(),
634  mesh.getEdgeManager(),
635  mesh.getFaceManager(),
636  subRegion,
637  finiteElement,
638  meanTotalStressIncrement_k,
639  averageMeanTotalStressIncrement_k );
640  } );
641  } );
642  } );
643  }
644 
645  virtual void updateBulkDensity( ElementSubRegionBase & subRegion ) = 0;
646 
647  virtual void validateNonlinearAcceleration() override
648  {
649  if( MpiWrapper::commSize( MPI_COMM_GEOS ) > 1 )
650  {
651  GEOS_ERROR( "Nonlinear acceleration is not implemented for MPI runs" );
652  }
653  }
654 
657 
660 
662  stabilization::StabilizationType m_stabilizationType;
663 
666 
669 
671  array1d< real64 > m_s0; // Accelerated averageMeanTotalStresIncrement @ outer iteration v ( two iterations ago )
672  array1d< real64 > m_s1; // Accelerated averageMeanTotalStresIncrement @ outer iteration v + 1 ( previous iteration )
673  array1d< real64 > m_s1_tilde; // Unaccelerated averageMeanTotalStresIncrement @ outer iteration v + 1 ( previous iteration )
674  array1d< real64 > m_s2; // Accelerated averageMeanTotalStresIncrement @ outer iteration v + 2 ( current iteration )
675  array1d< real64 > m_s2_tilde; // Unaccelerated averageMeanTotalStresIncrement @ outer iteration v + 1 ( current iteration )
676  real64 m_omega0; // Old Aitken relaxation factor
677  real64 m_omega1; // New Aitken relaxation factor
678 
679 };
680 
681 } /* namespace geos */
682 
683 #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.
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:44
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.
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 initializePostInitialConditionsPreSubGroups() override
Called by InitializePostInitialConditions() prior to initializing sub-Groups.
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
virtual void initializePostInitialConditionsPreSubGroups()
Called by InitializePostInitialConditions() prior to initializing sub-Groups.
Definition: Group.hpp:1554
@ 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:401
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:410
ENUM_STRINGS(LinearSolverParameters::SolverType, "direct", "cg", "gmres", "fgmres", "bicgstab", "preconditioner")
Declare strings associated with enumeration values.
std::set< T > set
A set of local indices.
Definition: DataTypes.hpp:271
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
std::int32_t integer
Signed integer type.
Definition: DataTypes.hpp:81
@ None
No Schur complement - just block-GS/block-Jacobi preconditioner.
Array< T, 1 > array1d
Alias for 1D array.
Definition: DataTypes.hpp:184
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