Pt/exp branch

inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_smoke.xml

@@ -49,7 +49,7 @@
       zCoords="{ 0, 1 }"
       nx="{ 5, 52, 5 }"
       ny="{ 5, 52, 5 }"
-      nz="{ 2 }"
+      nz="{ 3 }"
       cellBlockNames="{ cb1 }"/>
   </Mesh>
   <!-- SPHINX_MESH_END -->

src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.cpp

@@ -194,7 +194,7 @@ void CO2BrineFluid< PHASE1, PHASE2, FLASH >::checkTablesParameters( real64 const
 template< typename PHASE1, typename PHASE2, typename FLASH >
 void CO2BrineFluid< PHASE1, PHASE2, FLASH >::initializePreSubGroups()
 {
-  GEOS_THROW_IF( this->catalogName() == CO2BrineEzrokhiThermalFluid::catalogName(),
+  GEOS_THROW_IF( this->getCatalogName() == CO2BrineEzrokhiThermalFluid::catalogName(),
                  GEOS_FMT( "The `{}` model is disabled for now. Please use the other thermal CO2-brine model instead: `{}`",
                            CO2BrineEzrokhiThermalFluid::catalogName(),
                            CO2BrinePhillipsThermalFluid::catalogName() ),

src/coreComponents/constitutive/permeability/PermeabilityBase.cpp

@@ -68,6 +68,5 @@ void PermeabilityBase::allocateConstitutiveData( dataRepository::Group & parent,
   ConstitutiveBase::allocateConstitutiveData( parent, numConstitutivePointsPerParentIndex );
 }
 
-REGISTER_CATALOG_ENTRY( ConstitutiveBase, PermeabilityBase, string const &, Group * const )
 }
 } /* namespace geos */

src/coreComponents/constitutive/permeability/PermeabilityBase.hpp

@@ -113,10 +113,6 @@ class PermeabilityBase : public ConstitutiveBase
   virtual void allocateConstitutiveData( dataRepository::Group & parent,
                                          localIndex const numConstitutivePointsPerParentIndex ) override;
 
-  static string catalogName() { return "PermeabilityBase"; }
-
-  virtual string getCatalogName() const override { return catalogName(); }
-
   /**
    * @brief Const/non-mutable accessor for permeability.
    * @return Accessor

src/coreComponents/constitutive/solid/porosity/PorosityBase.hpp

@@ -137,10 +137,6 @@ class PorosityBase : public ConstitutiveBase
   virtual void allocateConstitutiveData( dataRepository::Group & parent,
                                          localIndex const numConstitutivePointsPerParentIndex ) override;
 
-  static string catalogName() { return "PorosityBase"; }
-
-  virtual string getCatalogName() const override { return catalogName(); }
-
   struct viewKeyStruct : public ConstitutiveBase::viewKeyStruct
   {
     static constexpr char const * defaultReferencePorosityString() { return "defaultReferencePorosity"; }

src/coreComponents/events/EventBase.hpp

@@ -44,12 +44,6 @@ class EventBase : public ExecutableGroup
   /// Destructor
   virtual ~EventBase() override;
 
-  /**
-   * @brief Catalog name interface.
-   * @return This type's catalog name.
-   **/
-  static string catalogName() { return "EventBase"; }
-
   /**
    * @brief If the event forecast is equal to 1, then signal the targets to prepare for execution
    *        during the next cycle.

src/coreComponents/events/tasks/TaskBase.hpp

@@ -39,12 +39,6 @@ class TaskBase : public ExecutableGroup
                      Group * const parent );
   virtual ~TaskBase( ) override;
 
-  /**
-   * @brief Catalog name interface
-   * @return This type's catalog name
-   */
-  static string catalogName() { return "TaskBase"; }
-
   /// The catalog interface type for TaskBase
   using CatalogInterface = dataRepository::CatalogInterface< TaskBase, string const &, Group * const >;
   /**

src/coreComponents/fileIO/Outputs/OutputBase.hpp

@@ -39,12 +39,6 @@ class OutputBase : public ExecutableGroup
   /// Destructor
   virtual ~OutputBase() override;
 
-  /**
-   * @brief Catalog name interface.
-   * @return This type's catalog name.
-   **/
-  static string catalogName() { return "OutputBase"; }
-
   /**
    * @brief Setter for the output directory
    * @param  outputDir The output directory

src/coreComponents/functions/FunctionBase.hpp

@@ -60,12 +60,6 @@ class FunctionBase : public dataRepository::Group
    */
   virtual ~FunctionBase() override = default;
 
-  /**
-   * @brief Static Factory Catalog Functions
-   * @return the catalog name
-   */
-  static string catalogName() { return "FunctionBase"; }
-
   /**
    * @brief Function initialization
    */

src/coreComponents/mesh/generators/MeshGeneratorBase.hpp

@@ -53,12 +53,6 @@ class MeshGeneratorBase : public dataRepository::Group
   explicit MeshGeneratorBase( string const & name,
                               Group * const parent );
 
-  /**
-   * @brief Return the name of the MeshGenerator in object catalog.
-   * @return string that contains the catalog name of the MeshGenerator
-   */
-  static string catalogName() { return "MeshGeneratorBase"; }
-
   /// This function is used to expand any catalogs in the data structure
   virtual void expandObjectCatalogs() override;
 

src/coreComponents/mesh/generators/WellGeneratorBase.hpp

@@ -46,12 +46,6 @@ class WellGeneratorBase : public WellGeneratorABC
   WellGeneratorBase( const string & name,
                      Group * const parent );
 
-  /**
-   * @brief Get the catalog name.
-   * @return the name of this type in the catalog
-   */
-  static string catalogName() { return "WellGeneratorBase"; }
-
   /// This function is used to expand any catalogs in the data structure
   virtual void expandObjectCatalogs() override;
 

src/coreComponents/physicsSolvers/NonlinearSolverParameters.cpp

@@ -60,6 +60,11 @@ NonlinearSolverParameters::NonlinearSolverParameters( string const & name,
     setInputFlag( InputFlags::OPTIONAL ).
     setDescription( "Iteration when line search starts." );
 
+  registerWrapper( viewKeysStruct::lineSearchResidualFactorString(), &m_lineSearchResidualFactor ).
+    setApplyDefaultValue( 1.0 ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setDescription( "Factor to determine residual increase (recommended values: 1.1 (conservative), 2.0 (relaxed), 10.0 (aggressive))." );
+
   registerWrapper( viewKeysStruct::normTypeString(), &m_normType ).
     setInputFlag( InputFlags::FALSE ).
     setApplyDefaultValue( solverBaseKernels::NormType::Linf ).
@@ -165,7 +170,7 @@ NonlinearSolverParameters::NonlinearSolverParameters( string const & name,
     setApplyDefaultValue( 0 ).
     setDescription( "Flag to decide whether to iterate between sequentially coupled solvers or not." );
 
-  this->registerWrapper( viewKeysStruct::nonlinearAccelerationTypeString(), &m_nonlinearAccelerationType ).
+  registerWrapper( viewKeysStruct::nonlinearAccelerationTypeString(), &m_nonlinearAccelerationType ).
     setApplyDefaultValue( NonlinearAccelerationType::None ).
     setInputFlag( dataRepository::InputFlags::OPTIONAL ).
     setDescription( "Nonlinear acceleration type for sequential solver." );
@@ -178,6 +183,9 @@ void NonlinearSolverParameters::postProcessInput()
                         getWrapperDataContext( viewKeysStruct::timeStepIncreaseIterLimString() ) <<
                         ": should be smaller than " << viewKeysStruct::timeStepDecreaseIterLimString() );
 
+  GEOS_ERROR_IF_LE_MSG( m_lineSearchResidualFactor, 0.0,
+                        getWrapperDataContext( viewKeysStruct::lineSearchResidualFactorString() ) << ": should be positive" );
+
   if( getLogLevel() > 0 )
   {
     GEOS_LOG_RANK_0( "Nonlinear solver parameters:" );
@@ -187,6 +195,8 @@ void NonlinearSolverParameters::postProcessInput()
       GEOS_LOG_RANK_0( GEOS_FMT( "  Line search interpolation type = {}", EnumStrings< LineSearchInterpolationType >::toString( m_lineSearchInterpType ) ) );
       GEOS_LOG_RANK_0( GEOS_FMT( "  Line search maximum number of cuts = {}", m_lineSearchMaxCuts ) );
       GEOS_LOG_RANK_0( GEOS_FMT( "  Line search cut factor = {}", m_lineSearchCutFactor ) );
+      GEOS_LOG_RANK_0( GEOS_FMT( "  Line starting iteration = {}", m_lineSearchStartingIteration ) );
+      GEOS_LOG_RANK_0( GEOS_FMT( "  Line residual increase factor = {}", m_lineSearchResidualFactor ) );
     }
     GEOS_LOG_RANK_0( GEOS_FMT( "  Norm type (flow solver) = {}", EnumStrings< solverBaseKernels::NormType >::toString( m_normType ) ) );
     GEOS_LOG_RANK_0( GEOS_FMT( "  Minimum residual normalizer = {}", m_minNormalizer ) );

src/coreComponents/physicsSolvers/NonlinearSolverParameters.hpp

@@ -64,6 +64,7 @@ class NonlinearSolverParameters : public dataRepository::Group
     m_lineSearchMaxCuts = params.m_lineSearchMaxCuts;
     m_lineSearchCutFactor = params.m_lineSearchCutFactor;
     m_lineSearchStartingIteration = params.m_lineSearchStartingIteration;
+    m_lineSearchResidualFactor = params.m_lineSearchResidualFactor;
 
     m_newtonTol = params.m_newtonTol;
     m_maxIterNewton = params.m_maxIterNewton;
@@ -104,6 +105,7 @@ class NonlinearSolverParameters : public dataRepository::Group
     static constexpr char const * lineSearchCutFactorString()     { return "lineSearchCutFactor"; }
     static constexpr char const * lineSearchInterpolationTypeString() { return "lineSearchInterpolationType"; }
     static constexpr char const * lineSearchStartingIterationString() { return "lineSearchStartingIteration"; }
+    static constexpr char const * lineSearchResidualFactorString() { return "lineSearchResidualFactor"; }
 
     static constexpr char const * normTypeString()                { return "normType"; }
     static constexpr char const * minNormalizerString()           { return "minNormalizer"; }
@@ -261,6 +263,9 @@ class NonlinearSolverParameters : public dataRepository::Group
   /// Iteration when line search starts
   integer m_lineSearchStartingIteration;
 
+  /// Factor to determine residual increase
+  real64 m_lineSearchResidualFactor;
+
   /// Norm used to check the nonlinear loop convergence
   solverBaseKernels::NormType m_normType;
 

src/coreComponents/physicsSolvers/SolverBase.cpp

@@ -896,7 +896,8 @@ bool SolverBase::solveNonlinearSystem( real64 const & time_n,
 
     // do line search in case residual has increased
     if( m_nonlinearSolverParameters.m_lineSearchAction != NonlinearSolverParameters::LineSearchAction::None
-        && residualNorm > lastResidual && newtonIter >= m_nonlinearSolverParameters.m_lineSearchStartingIteration )
+        && residualNorm > lastResidual * m_nonlinearSolverParameters.m_lineSearchResidualFactor
+        && newtonIter >= m_nonlinearSolverParameters.m_lineSearchStartingIteration )
     {
       bool lineSearchSuccess = false;
       if( m_nonlinearSolverParameters.m_lineSearchInterpType == NonlinearSolverParameters::LineSearchInterpolationType::Linear )

src/coreComponents/physicsSolvers/SolverBase.hpp

@@ -49,8 +49,6 @@ class SolverBase : public ExecutableGroup
   SolverBase & operator=( SolverBase const & ) = delete;
   SolverBase & operator=( SolverBase && ) = delete;
 
-  static string catalogName() { return "SolverBase"; }
-
   /**
    * @return Get the final class Catalog name
    */

src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp

@@ -537,83 +537,90 @@ real64 SolidMechanicsEmbeddedFractures::calculateResidualNorm( real64 const & ti
 
   if( !m_useStaticCondensation )
   {
+    real64 const fractureResidualNorm = calculateFractureResidualNorm( domain, dofManager, localRhs );
 
-    string const jumpDofKey = dofManager.getKey( contact::dispJump::key() );
+    return sqrt( solidResidualNorm * solidResidualNorm + fractureResidualNorm * fractureResidualNorm );
+  }
+  else
+  {
+    return solidResidualNorm;
+  }
+}
 
-    globalIndex const rankOffset = dofManager.rankOffset();
+real64 SolidMechanicsEmbeddedFractures::calculateFractureResidualNorm( DomainPartition const & domain,
+                                                                       DofManager const & dofManager,
+                                                                       arrayView1d< real64 const > const & localRhs ) const
+{
+  string const jumpDofKey = dofManager.getKey( contact::dispJump::key() );
+
+  globalIndex const rankOffset = dofManager.rankOffset();
 
-    RAJA::ReduceSum< parallelDeviceReduce, real64 > localSum( 0.0 );
+  RAJA::ReduceSum< parallelDeviceReduce, real64 > localSum( 0.0 );
 
-    // globalResidualNorm[0]: the sum of all the local sum(rhs^2).
-    // globalResidualNorm[1]: max of max force of each rank. Basically max force globally
-    real64 globalResidualNorm[2] = {0, 0};
+  // globalResidualNorm[0]: the sum of all the local sum(rhs^2).
+  // globalResidualNorm[1]: max of max force of each rank. Basically max force globally
+  real64 globalResidualNorm[2] = {0.0, 0.0};
 
-    // Fracture residual
-    forFractureRegionOnMeshTargets( domain.getMeshBodies(), [&] ( SurfaceElementRegion const & fractureRegion )
+  // Fracture residual
+  forFractureRegionOnMeshTargets( domain.getMeshBodies(), [&] ( SurfaceElementRegion const & fractureRegion )
+  {
+    fractureRegion.forElementSubRegions< SurfaceElementSubRegion >( [&]( SurfaceElementSubRegion const & subRegion )
     {
-      fractureRegion.forElementSubRegions< SurfaceElementSubRegion >( [&]( SurfaceElementSubRegion const & subRegion )
-      {
-        arrayView1d< globalIndex const > const &
-        dofNumber = subRegion.getReference< array1d< globalIndex > >( jumpDofKey );
-        arrayView1d< integer const > const & ghostRank = subRegion.ghostRank();
+      arrayView1d< globalIndex const > const &
+      dofNumber = subRegion.getReference< array1d< globalIndex > >( jumpDofKey );
+      arrayView1d< integer const > const & ghostRank = subRegion.ghostRank();
 
-        forAll< parallelDevicePolicy<> >( subRegion.size(),
-                                          [localRhs, localSum, dofNumber, rankOffset, ghostRank] GEOS_HOST_DEVICE ( localIndex const k )
+      forAll< parallelDevicePolicy<> >( subRegion.size(),
+                                        [localRhs, localSum, dofNumber, rankOffset, ghostRank] GEOS_HOST_DEVICE ( localIndex const k )
+      {
+        if( ghostRank[k] < 0 )
         {
-          if( ghostRank[k] < 0 )
+          localIndex const localRow = LvArray::integerConversion< localIndex >( dofNumber[k] - rankOffset );
+          for( int i = 0; i < 3; ++i )
           {
-            localIndex const localRow = LvArray::integerConversion< localIndex >( dofNumber[k] - rankOffset );
-            for( localIndex i = 0; i < 3; ++i )
-            {
-              localSum += localRhs[localRow + i] * localRhs[localRow + i];
-            }
+            localSum += localRhs[localRow + i] * localRhs[localRow + i];
           }
-        } );
-
+        }
       } );
 
-      real64 const localResidualNorm[2] = { localSum.get(), SolidMechanicsLagrangianFEM::getMaxForce() };
+    } );
 
+    real64 const localResidualNorm[2] = { localSum.get(), SolidMechanicsLagrangianFEM::getMaxForce() };
 
-      int const rank     = MpiWrapper::commRank( MPI_COMM_GEOSX );
-      int const numRanks = MpiWrapper::commSize( MPI_COMM_GEOSX );
-      array1d< real64 > globalValues( numRanks * 2 );
+    int const rank     = MpiWrapper::commRank( MPI_COMM_GEOSX );
+    int const numRanks = MpiWrapper::commSize( MPI_COMM_GEOSX );
+    array1d< real64 > globalValues( numRanks * 2 );
 
-      // Everything is done on rank 0
-      MpiWrapper::gather( localResidualNorm,
-                          2,
-                          globalValues.data(),
-                          2,
-                          0,
-                          MPI_COMM_GEOSX );
+    // Everything is done on rank 0
+    MpiWrapper::gather( localResidualNorm,
+                        2,
+                        globalValues.data(),
+                        2,
+                        0,
+                        MPI_COMM_GEOSX );
 
-      if( rank==0 )
+    if( rank==0 )
+    {
+      for( int r=0; r<numRanks; ++r )
       {
-        for( int r=0; r<numRanks; ++r )
-        {
-          // sum/max across all ranks
-          globalResidualNorm[0] += globalValues[r*2];
-          globalResidualNorm[1] = std::max( globalResidualNorm[1], globalValues[r*2+1] );
-        }
+        // sum/max across all ranks
+        globalResidualNorm[0] += globalValues[r*2];
+        globalResidualNorm[1] = std::max( globalResidualNorm[1], globalValues[r*2+1] );
       }
+    }
 
-      MpiWrapper::bcast( globalResidualNorm, 2, 0, MPI_COMM_GEOSX );
-    } );
-
-    real64 const fractureResidualNorm = sqrt( globalResidualNorm[0] )/(globalResidualNorm[1]+1);  // the + 1 is for the first
-    // time-step when maxForce = 0;
+    MpiWrapper::bcast( globalResidualNorm, 2, 0, MPI_COMM_GEOSX );
+  } );
 
-    if( getLogLevel() >= 1 && logger::internal::rank==0 )
-    {
-      std::cout << GEOS_FMT( "        ( RFracture ) = ( {:4.2e} )", fractureResidualNorm );
-    }
+  real64 const fractureResidualNorm = sqrt( globalResidualNorm[0] )/(globalResidualNorm[1]+1);  // the + 1 is for the first
+                                                                                                // time-step when maxForce = 0;
 
-    return sqrt( solidResidualNorm * solidResidualNorm + fractureResidualNorm * fractureResidualNorm );
-  }
-  else
+  if( getLogLevel() >= 1 && logger::internal::rank==0 )
   {
-    return solidResidualNorm;
+    std::cout << GEOS_FMT( "        ( RFracture ) = ( {:4.2e} )", fractureResidualNorm );
   }
+
+  return fractureResidualNorm;
 }
 
 void SolidMechanicsEmbeddedFractures::applySystemSolution( DofManager const & dofManager,

src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.hpp

@@ -77,6 +77,10 @@ class SolidMechanicsEmbeddedFractures : public ContactSolverBase
                          DofManager const & dofManager,
                          arrayView1d< real64 const > const & localRhs ) override;
 
+  real64 calculateFractureResidualNorm( DomainPartition const & domain,
+                                        DofManager const & dofManager,
+                                        arrayView1d< real64 const > const & localRhs ) const;
+
   virtual void
   applySystemSolution( DofManager const & dofManager,
                        arrayView1d< real64 const > const & localSolution,