QuasiDynamicEQRK32.hpp

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 TotalEnergies
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQRK32_HPP
#define GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQRK32_HPP
 
#include "physicsSolvers/PhysicsSolverBase.hpp"
#include "kernels/RateAndStateKernels.hpp"
 
namespace geos
{
 
class QuasiDynamicEQRK32 : public PhysicsSolverBase
{
public:
  QuasiDynamicEQRK32() = delete;
 
  QuasiDynamicEQRK32( const string & name,
                      Group * const parent );
 
  virtual ~QuasiDynamicEQRK32() override;
 
  static string catalogName() { return "QuasiDynamicEQRK32"; }
 
  virtual string getCatalogName() const override { return catalogName(); }
 
  virtual void registerDataOnMesh( Group & meshBodies ) override;
 
  struct viewKeyStruct : public PhysicsSolverBase::viewKeyStruct
  {
    constexpr static char const * stressSolverNameString() { return "stressSolverName"; }
    constexpr static char const * frictionLawNameString() { return "frictionLawName"; }
    constexpr static char const * shearImpedanceString() { return "shearImpedance"; }
    constexpr static char const * timeStepTol() { return "timeStepTol"; }
  };
 
  virtual real64 solverStep( real64 const & time_n,
                             real64 const & dt,
                             integer const cycleNumber,
                             DomainPartition & domain ) override final;
 
 
  virtual real64 setNextDt( real64 const & currentDt,
                            DomainPartition & domain ) override final;
 
  void stepRateStateODEInitialSubstage( real64 const dt, DomainPartition & domain ) const;
 
  void stepRateStateODESubstage( integer const stageIndex,
                                 real64 const dt,
                                 DomainPartition & domain ) const;
 
  void stepRateStateODEAndComputeError( real64 const dt, DomainPartition & domain ) const;
 
  real64 updateStresses( real64 const & time_n,
                         real64 const & dt,
                         const int cycleNumber,
                         DomainPartition & domain ) const;
 
  void updateSlipVelocity( real64 const & time_n,
                           real64 const & dt,
                           DomainPartition & domain ) const;
 
  void saveState( DomainPartition & domain ) const;
 
private:
 
  virtual void postInputInitialization() override;
 
 
  PhysicsSolverBase * m_stressSolver;
 
  string m_stressSolverName;
 
  real64 m_shearImpedance;
 
  // TODO: The specific type should not be hardcoded!
  // Should be possible to change RK-method based on the table.
  rateAndStateKernels::BogackiShampine32Table m_butcherTable;
 
  bool m_successfulStep; // Flag indicating if the adative time step was accepted
 
  class PIDController
  {
public:
 
    GEOS_HOST_DEVICE
    PIDController( std::array< const real64, 3 > const & cparams,
                   const real64 atol,
                   const real64 rtol,
                   const real64 safety ):
      controlParameters{ cparams },
      absTol( atol ),
      relTol( rtol ),
      acceptSafety( safety ),
      errors{ {0.0, 0.0, 0.0} }
    {}
 
    PIDController( PIDController const & ) = default;
 
    PIDController( PIDController && ) = default;
 
    PIDController() = delete;
 
    PIDController & operator=( PIDController const & ) = delete;
 
    PIDController & operator=( PIDController && ) =  delete;
 
    const std::array< const real64, 3 > controlParameters; // Controller parameters
 
    real64 const absTol; // absolut tolerence
 
    real64 const relTol; // relative tolerence
 
    real64 const acceptSafety; // Acceptance safety
 
    std::array< real64, 3 > errors; // Errors for current and two previous updates
                                    // stored as [n+1, n, n-1]
 
    real64 computeUpdateFactor( integer const algHighOrder, integer const algLowOrder )
    {
      // PID error controller + limiter
      real64 const k = LvArray::math::min( algHighOrder, algLowOrder ) + 1.0;
      real64 const eps0 = 1.0/(errors[0] + std::numeric_limits< real64 >::epsilon()); // n + 1
      real64 const eps1 = 1.0/(errors[1] + std::numeric_limits< real64 >::epsilon()); // n
      real64 const eps2 = 1.0/(errors[2] + std::numeric_limits< real64 >::epsilon()); // n-1
      // Compute update factor eps0^(beta0/k)*eps1^(beta1/k)*eps2^(beta2/k) where
      // beta0 - beta2 are the control parameters. Also apply limiter to smoothen changes.
      // Limiter is 1.0 + atan(x - 1.0). Here use atan(x) = atan2(x, 1.0).
      return 1.0 + LvArray::math::atan2( pow( eps0, controlParameters[0] / k ) *
                                         pow( eps1, controlParameters[1] / k ) *
                                         pow( eps2, controlParameters[2] / k ) - 1.0, 1.0 );
    }
  };
 
  PIDController m_controller;
 
 
  class SpringSliderParameters
  {
public:
 
    GEOS_HOST_DEVICE
    SpringSliderParameters( real64 const normalTraction, real64 const a, real64 const b, real64 const Dc ):
      tauRate( 1e-4 ),
      springStiffness( 0.0 )
    {
      real64 const criticalStiffness = normalTraction * (b - a) / Dc;
      springStiffness = 0.9 * criticalStiffness;
    }
 
    SpringSliderParameters( SpringSliderParameters const & ) = default;
 
    SpringSliderParameters( SpringSliderParameters && ) = default;
 
    SpringSliderParameters() = delete;
 
    SpringSliderParameters & operator=( SpringSliderParameters const & ) = delete;
 
    SpringSliderParameters & operator=( SpringSliderParameters && ) =  delete;
 
    real64 tauRate;
 
    real64 springStiffness;
  };
};
 
} /* namespace geos */
 
#endif /* GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQRK32_HPP */