diff --git a/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml b/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml new file mode 100644 index 0000000000..08fc7cb993 --- /dev/null +++ b/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml @@ -0,0 +1,251 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/inputFiles/inducedSeismicity/SpringSliderExplicit_base.xml b/inputFiles/inducedSeismicity/SpringSliderExplicit_base.xml index a0c0381fb2..4e920ca8fe 100644 --- a/inputFiles/inducedSeismicity/SpringSliderExplicit_base.xml +++ b/inputFiles/inducedSeismicity/SpringSliderExplicit_base.xml @@ -1,14 +1,14 @@ - - + + filename="springSlider"/> diff --git a/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml b/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml index 2de09b79a0..bdcbf1280b 100644 --- a/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml +++ b/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml @@ -23,6 +23,32 @@ maxEventDt="1e4" target="/Solvers/SpringSlider"/> + + + + + + + + + + , 0, NOPLOT, @@ -98,6 +98,14 @@ DECLARE_FIELD( deltaSlip, WRITE_AND_READ, "Slip increment" ); +DECLARE_FIELD( deltaSlip_n, + "deltaSlip_n", + array2d< real64 >, + 0.0, + NOPLOT, + WRITE_AND_READ, + "Initial slip increment at this time step" ); + DECLARE_FIELD( deltaDispJump, "deltaDisplacementJump", array2d< real64 >, diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/CMakeLists.txt b/src/coreComponents/physicsSolvers/inducedSeismicity/CMakeLists.txt index a4be651a3f..5ff02a897a 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/CMakeLists.txt @@ -1,14 +1,16 @@ # Specify solver headers set( physicsSolvers_headers ${physicsSolvers_headers} + inducedSeismicity/ExplicitQDRateAndState.hpp + inducedSeismicity/ImplicitQDRateAndState.hpp inducedSeismicity/inducedSeismicityFields.hpp + inducedSeismicity/QuasiDynamicEarthQuake.hpp inducedSeismicity/rateAndStateFields.hpp - inducedSeismicity/QuasiDynamicEQ.hpp - inducedSeismicity/ImplicitQDRateAndState.hpp - inducedSeismicity/QuasiDynamicEQRK32.hpp inducedSeismicity/SeismicityRate.hpp inducedSeismicity/SpringSlider.hpp - inducedSeismicity/kernels/RateAndStateKernels.hpp + inducedSeismicity/kernels/ExplicitRateAndStateKernels.hpp + inducedSeismicity/kernels/ImplicitRateAndStateKernels.hpp + inducedSeismicity/kernels/RateAndStateKernelsBase.hpp inducedSeismicity/kernels/SeismicityRateKernels.hpp inducedSeismicity/tractionUpdateWrapper/FaultTractionUpdate.hpp inducedSeismicity/tractionUpdateWrapper/FaultTractionUpdateBase.hpp @@ -19,9 +21,9 @@ set( physicsSolvers_headers # Specify solver sources set( physicsSolvers_sources ${physicsSolvers_sources} - inducedSeismicity/QuasiDynamicEQ.cpp + inducedSeismicity/ExplicitQDRateAndState.cpp inducedSeismicity/ImplicitQDRateAndState.cpp - inducedSeismicity/QuasiDynamicEQRK32.cpp + inducedSeismicity/QuasiDynamicEarthQuake.cpp inducedSeismicity/SeismicityRate.cpp inducedSeismicity/SpringSlider.cpp inducedSeismicity/tractionUpdateWrapper/OneWayCoupledTractionUpdate.cpp diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp similarity index 71% rename from src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.cpp rename to src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp index 7b56911893..d1ab8f4bec 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.cpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp @@ -14,14 +14,14 @@ */ /** - * @file QuasiDynamicEQRK32.cpp + * @file ExplicitQDRateAndState.cpp */ -#include "QuasiDynamicEQRK32.hpp" +#include "ExplicitQDRateAndState.hpp" #include "dataRepository/InputFlags.hpp" #include "mesh/DomainPartition.hpp" -#include "kernels/RateAndStateKernels.hpp" +#include "kernels/ExplicitRateAndStateKernels.hpp" #include "rateAndStateFields.hpp" #include "physicsSolvers/contact/ContactFields.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" @@ -34,11 +34,9 @@ using namespace fields; using namespace constitutive; using namespace rateAndStateKernels; -QuasiDynamicEQRK32::QuasiDynamicEQRK32( const string & name, +ExplicitQDRateAndState::ExplicitQDRateAndState( const string & name, Group * const parent ): PhysicsSolverBase( name, parent ), - m_stressSolver( nullptr ), - m_stressSolverName( "SpringSlider" ), m_shearImpedance( 0.0 ), m_butcherTable( BogackiShampine32Table()), // TODO: The butcher table should be specified in the XML input. m_successfulStep( false ), @@ -48,31 +46,19 @@ QuasiDynamicEQRK32::QuasiDynamicEQRK32( const string & name, this->registerWrapper( viewKeyStruct::shearImpedanceString(), &m_shearImpedance ). setInputFlag( InputFlags::REQUIRED ). setDescription( "Shear impedance." ); - - this->registerWrapper( viewKeyStruct::stressSolverNameString(), &m_stressSolverName ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Name of solver for computing stress. If empty, the spring-slider model is run." ); } -void QuasiDynamicEQRK32::postInputInitialization() +void ExplicitQDRateAndState::postInputInitialization() { - - // Initialize member stress solver as specified in XML input - if( !m_stressSolverName.empty() ) - { - m_stressSolver = &this->getParent().getGroup< PhysicsSolverBase >( m_stressSolverName ); - } - PhysicsSolverBase::postInputInitialization(); } -QuasiDynamicEQRK32::~QuasiDynamicEQRK32() +ExplicitQDRateAndState::~ExplicitQDRateAndState() { // TODO Auto-generated destructor stub } - -void QuasiDynamicEQRK32::registerDataOnMesh( Group & meshBodies ) +void ExplicitQDRateAndState::registerDataOnMesh( Group & meshBodies ) { PhysicsSolverBase::registerDataOnMesh( meshBodies ); @@ -97,52 +83,19 @@ void QuasiDynamicEQRK32::registerDataOnMesh( Group & meshBodies ) setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); subRegion.registerField< rateAndState::slipVelocity_n >( getName() ). setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); - subRegion.registerField< rateAndState::deltaSlip >( getName() ). - setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); - subRegion.registerField< rateAndState::deltaSlip_n >( getName() ). - setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); // Runge-Kutta stage rates and error integer const numRKComponents = 3; subRegion.registerField< rateAndState::rungeKuttaStageRates >( getName() ).reference().resizeDimension< 1, 2 >( m_butcherTable.numStages, numRKComponents ); subRegion.registerField< rateAndState::error >( getName() ).reference().resizeDimension< 1 >( numRKComponents ); - - - if( !subRegion.hasWrapper( contact::dispJump::key() )) - { - // 3-component functions on fault - string const labels3Comp[3] = { "normal", "tangent1", "tangent2" }; - subRegion.registerField< contact::dispJump >( getName() ). - setDimLabels( 1, labels3Comp ). - reference().resizeDimension< 1 >( 3 ); - subRegion.registerField< contact::dispJump_n >( getName() ). - setDimLabels( 1, labels3Comp ). - reference().resizeDimension< 1 >( 3 ); - subRegion.registerField< contact::traction >( getName() ). - setDimLabels( 1, labels3Comp ). - reference().resizeDimension< 1 >( 3 ); - subRegion.registerField< contact::traction_n >( getName() ). - setDimLabels( 1, labels3Comp ). - reference().resizeDimension< 1 >( 3 ); - - subRegion.registerWrapper< string >( viewKeyStruct::frictionLawNameString() ). - setPlotLevel( PlotLevel::NOPLOT ). - setRestartFlags( RestartFlags::NO_WRITE ). - setSizedFromParent( 0 ); - - string & frictionLawName = subRegion.getReference< string >( viewKeyStruct::frictionLawNameString() ); - frictionLawName = PhysicsSolverBase::getConstitutiveName< FrictionBase >( subRegion ); - GEOS_ERROR_IF( frictionLawName.empty(), GEOS_FMT( "{}: FrictionBase model not found on subregion {}", - getDataContext(), subRegion.getDataContext() ) ); - } } ); } ); } -real64 QuasiDynamicEQRK32::solverStep( real64 const & time_n, - real64 const & dt, - int const cycleNumber, - DomainPartition & domain ) +real64 ExplicitQDRateAndState::solverStep( real64 const & time_n, + real64 const & dt, + int const cycleNumber, + DomainPartition & domain ) { if( cycleNumber == 0 ) { @@ -208,7 +161,7 @@ real64 QuasiDynamicEQRK32::solverStep( real64 const & time_n, return dtAdaptive; } -void QuasiDynamicEQRK32::stepRateStateODEInitialSubstage( real64 const dt, DomainPartition & domain ) const +void ExplicitQDRateAndState::stepRateStateODEInitialSubstage( real64 const dt, DomainPartition & domain ) const { forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, @@ -247,7 +200,7 @@ void QuasiDynamicEQRK32::stepRateStateODEInitialSubstage( real64 const dt, Domai } ); } -void QuasiDynamicEQRK32::stepRateStateODESubstage( integer const stageIndex, +void ExplicitQDRateAndState::stepRateStateODESubstage( integer const stageIndex, real64 const dt, DomainPartition & domain ) const { @@ -276,7 +229,7 @@ void QuasiDynamicEQRK32::stepRateStateODESubstage( integer const stageIndex, } ); } -void QuasiDynamicEQRK32::stepRateStateODEAndComputeError( real64 const dt, DomainPartition & domain ) const +void ExplicitQDRateAndState::stepRateStateODEAndComputeError( real64 const dt, DomainPartition & domain ) const { forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, @@ -316,62 +269,7 @@ void QuasiDynamicEQRK32::stepRateStateODEAndComputeError( real64 const dt, Domai } ); } -real64 QuasiDynamicEQRK32::updateStresses( real64 const & time_n, - real64 const & dt, - const int cycleNumber, - DomainPartition & domain ) const -{ - GEOS_LOG_LEVEL_RANK_0( 1, "Stress solver" ); - // Call member variable stress solver to update the stress state - if( m_stressSolver ) - { - // 1. Solve the momentum balance - real64 const dtStress = m_stressSolver->solverStep( time_n, dt, cycleNumber, domain ); - - return dtStress; - } - else - { - // Spring-slider shear traction computation - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) - - { - mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, - [&]( localIndex const, - SurfaceElementSubRegion & subRegion ) - { - - arrayView2d< real64 const > const deltaSlip = subRegion.getField< rateAndState::deltaSlip >(); - arrayView2d< real64 > const traction = subRegion.getField< fields::contact::traction >(); - arrayView2d< real64 const > const traction_n = subRegion.getField< fields::contact::traction_n >(); - - string const & fricitonLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() ); - RateAndStateFriction const & frictionLaw = getConstitutiveModel< RateAndStateFriction >( subRegion, fricitonLawName ); - - RateAndStateFriction::KernelWrapper frictionKernelWrapper = frictionLaw.createKernelUpdates(); - - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) - { - SpringSliderParameters springSliderParameters = SpringSliderParameters( traction[k][0], - frictionKernelWrapper.getACoefficient( k ), - frictionKernelWrapper.getBCoefficient( k ), - frictionKernelWrapper.getDcCoefficient( k ) ); - - - traction[k][1] = traction_n[k][1] + springSliderParameters.tauRate * dt - - springSliderParameters.springStiffness * deltaSlip[k][0]; - traction[k][2] = traction_n[k][2] + springSliderParameters.tauRate * dt - - springSliderParameters.springStiffness * deltaSlip[k][1]; - } ); - } ); - } ); - return dt; - } -} - -void QuasiDynamicEQRK32::updateSlipVelocity( real64 const & time_n, +void ExplicitQDRateAndState::updateSlipVelocity( real64 const & time_n, real64 const & dt, DomainPartition & domain ) const { @@ -394,7 +292,7 @@ void QuasiDynamicEQRK32::updateSlipVelocity( real64 const & time_n, } ); } -void QuasiDynamicEQRK32::saveState( DomainPartition & domain ) const +void ExplicitQDRateAndState::saveState( DomainPartition & domain ) const { forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, @@ -407,13 +305,13 @@ void QuasiDynamicEQRK32::saveState( DomainPartition & domain ) const { arrayView1d< real64 const > const stateVariable = subRegion.getField< rateAndState::stateVariable >(); arrayView2d< real64 const > const slipVelocity = subRegion.getField< rateAndState::slipVelocity >(); - arrayView2d< real64 const > const deltaSlip = subRegion.getField< rateAndState::deltaSlip >(); + arrayView2d< real64 const > const deltaSlip = subRegion.getField< contact::deltaSlip >(); arrayView2d< real64 const > const dispJump = subRegion.getField< contact::dispJump >(); arrayView2d< real64 const > const traction = subRegion.getField< contact::traction >(); arrayView1d< real64 > const stateVariable_n = subRegion.getField< rateAndState::stateVariable_n >(); arrayView2d< real64 > const slipVelocity_n = subRegion.getField< rateAndState::slipVelocity_n >(); - arrayView2d< real64 > const deltaSlip_n = subRegion.getField< rateAndState::deltaSlip >(); + arrayView2d< real64 > const deltaSlip_n = subRegion.getField< contact::deltaSlip >(); arrayView2d< real64 > const dispJump_n = subRegion.getField< contact::dispJump_n >(); arrayView2d< real64 > const traction_n = subRegion.getField< contact::traction_n >(); @@ -429,7 +327,7 @@ void QuasiDynamicEQRK32::saveState( DomainPartition & domain ) const } ); } -real64 QuasiDynamicEQRK32::setNextDt( real64 const & currentDt, DomainPartition & domain ) +real64 ExplicitQDRateAndState::setNextDt( real64 const & currentDt, DomainPartition & domain ) { // Spring-slider shear traction computation @@ -473,6 +371,4 @@ real64 QuasiDynamicEQRK32::setNextDt( real64 const & currentDt, DomainPartition return nextDt; } -REGISTER_CATALOG_ENTRY( PhysicsSolverBase, QuasiDynamicEQRK32, string const &, dataRepository::Group * const ) - } // namespace geos diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.hpp similarity index 91% rename from src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.hpp rename to src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.hpp index 0979fd22dd..f8715ac050 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.hpp @@ -17,31 +17,26 @@ #define GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQRK32_HPP #include "physicsSolvers/PhysicsSolverBase.hpp" -#include "kernels/RateAndStateKernels.hpp" +#include "kernels/ExplicitRateAndStateKernels.hpp" namespace geos { -class QuasiDynamicEQRK32 : public PhysicsSolverBase +class ExplicitQDRateAndState : public PhysicsSolverBase { public: /// The default nullary constructor is disabled to avoid compiler auto-generation: - QuasiDynamicEQRK32() = delete; + ExplicitQDRateAndState() = delete; /// The constructor needs a user-defined "name" and a parent Group (to place this instance in the tree structure of classes) - QuasiDynamicEQRK32( const string & name, + ExplicitQDRateAndState( const string & name, Group * const parent ); /// Destructor - virtual ~QuasiDynamicEQRK32() override; - - static string catalogName() { return "QuasiDynamicEQRK32"; } - - /** - * @return Get the final class Catalog name - */ - virtual string getCatalogName() const override { return catalogName(); } + virtual ~ExplicitQDRateAndState() override; + static string derivedSolverPrefix() { return "Explicit";}; + /// This method ties properties with their supporting mesh virtual void registerDataOnMesh( Group & meshBodies ) override; @@ -90,11 +85,6 @@ class QuasiDynamicEQRK32 : public PhysicsSolverBase */ void stepRateStateODEAndComputeError( real64 const dt, DomainPartition & domain ) const; - real64 updateStresses( real64 const & time_n, - real64 const & dt, - const int cycleNumber, - DomainPartition & domain ) const; - /** * @brief Updates rate-and-state slip velocity * @param domain @@ -109,17 +99,14 @@ class QuasiDynamicEQRK32 : public PhysicsSolverBase */ void saveState( DomainPartition & domain ) const; -private: +protected: virtual void postInputInitialization() override; - - /// pointer to stress solver - PhysicsSolverBase * m_stressSolver; - - /// stress solver name - string m_stressSolverName; - + virtual real64 updateStresses( real64 const & time_n, + real64 const & dt, + const int cycleNumber, + DomainPartition & domain ) const = 0; /// shear impedance real64 m_shearImpedance; diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.cpp index 3fe2f0334a..620485771f 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.cpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.cpp @@ -21,7 +21,7 @@ #include "dataRepository/InputFlags.hpp" #include "mesh/DomainPartition.hpp" -#include "kernels/RateAndStateKernels.hpp" +#include "kernels/ImplicitRateAndStateKernels.hpp" #include "rateAndStateFields.hpp" #include "physicsSolvers/contact/ContactFields.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" @@ -80,6 +80,8 @@ void ImplicitQDRateAndState::registerDataOnMesh( Group & meshBodies ) string const labels2Comp[2] = {"tangent1", "tangent2" }; subRegion.registerField< rateAndState::slipVelocity >( getName() ). setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); + subRegion.registerField< rateAndState::slipVelocity_n >( getName() ). + setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); subRegion.registerField< rateAndState::shearTraction >( getName() ). setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); } ); @@ -106,13 +108,16 @@ void ImplicitQDRateAndState::applyInitialConditionsToFault( int const cycleNumbe { arrayView1d< real64 const > const slipRate = subRegion.getField< rateAndState::slipRate >(); arrayView1d< real64 const > const stateVariable = subRegion.getField< rateAndState::stateVariable >(); + arrayView2d< real64 const > const traction = subRegion.getField< contact::traction >(); arrayView1d< real64 > const stateVariable_n = subRegion.getField< rateAndState::stateVariable_n >(); arrayView1d< real64 > const slipRate_n = subRegion.getField< rateAndState::slipRate_n >(); + arrayView2d< real64 > const traction_n = subRegion.getField< contact::traction_n >(); forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) { slipRate_n [k] = slipRate[k]; stateVariable_n[k] = stateVariable[k]; + LvArray::tensorOps::copy< 3 >( traction_n[k], traction[k] ); } ); } ); } ); @@ -135,10 +140,10 @@ void ImplicitQDRateAndState::solveRateAndStateEquations( real64 const time_n, SurfaceElementSubRegion & subRegion ) { // solve rate and state equations. - createAndLaunch< ImplicitFixedStressRateAndStateKernel, parallelDevicePolicy<> >( subRegion, + rateAndStateKernels::createAndLaunch< rateAndStateKernels::ImplicitFixedStressRateAndStateKernel, parallelDevicePolicy<> >( subRegion, viewKeyStruct::frictionLawNameString(), m_shearImpedance, - maxIterNewton, newtonTol, + maxNewtonIter, newtonTol, time_n, dt ); // save old state updateSlip( subRegion, dt ); @@ -152,8 +157,9 @@ real64 ImplicitQDRateAndState::solverStep( real64 const & time_n, DomainPartition & domain ) { applyInitialConditionsToFault( cycleNumber, domain ); - updateStresses( dt, domain ); + updateStresses( time_n, dt, cycleNumber, domain ); solveRateAndStateEquations( time_n, dt, domain ); + saveState( domain ); return dt; } @@ -169,6 +175,41 @@ void ImplicitQDRateAndState::updateSlip( ElementSubRegionBase & subRegion, real6 } ); } +void ImplicitQDRateAndState::saveState( DomainPartition & domain ) const +{ + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + + { + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + arrayView1d< real64 const > const stateVariable = subRegion.getField< rateAndState::stateVariable >(); + arrayView2d< real64 const > const slipVelocity = subRegion.getField< rateAndState::slipVelocity >(); + arrayView2d< real64 const > const deltaSlip = subRegion.getField< contact::deltaSlip >(); + arrayView2d< real64 const > const dispJump = subRegion.getField< contact::dispJump >(); + arrayView2d< real64 const > const traction = subRegion.getField< contact::traction >(); + + arrayView1d< real64 > const stateVariable_n = subRegion.getField< rateAndState::stateVariable_n >(); + arrayView2d< real64 > const slipVelocity_n = subRegion.getField< rateAndState::slipVelocity_n >(); + arrayView2d< real64 > const deltaSlip_n = subRegion.getField< contact::deltaSlip >(); + arrayView2d< real64 > const dispJump_n = subRegion.getField< contact::dispJump_n >(); + arrayView2d< real64 > const traction_n = subRegion.getField< contact::traction_n >(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + stateVariable_n[k] = stateVariable[k]; + LvArray::tensorOps::copy< 2 >( deltaSlip_n[k], deltaSlip[k] ); + LvArray::tensorOps::copy< 2 >( slipVelocity_n[k], slipVelocity[k] ); + LvArray::tensorOps::copy< 3 >( dispJump_n[k], dispJump[k] ); + LvArray::tensorOps::copy< 3 >( traction_n[k], traction[k] ); + } ); + } ); + } ); +} + real64 ImplicitQDRateAndState::setNextDt( real64 const & currentDt, DomainPartition & domain ) { GEOS_UNUSED_VAR( currentDt ); diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.hpp index 8a79114068..b75e506566 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.hpp @@ -63,10 +63,17 @@ class ImplicitQDRateAndState : public PhysicsSolverBase integer const cycleNumber, DomainPartition & domain ) override final; + /** + * @brief save the current state + * @param domain + */ + void saveState( DomainPartition & domain ) const; protected: - virtual real64 updateStresses( real64 const dt, + virtual real64 updateStresses( real64 const & time_n, + real64 const & dt, + const int cycleNumber, DomainPartition & domain ) const = 0; virtual void solveRateAndStateEquations( real64 const time_n, diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.cpp index b81c868ff3..6b6d27c30a 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.cpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.cpp @@ -21,11 +21,12 @@ #include "dataRepository/InputFlags.hpp" #include "mesh/DomainPartition.hpp" -#include "kernels/RateAndStateKernels.hpp" #include "rateAndStateFields.hpp" #include "physicsSolvers/contact/ContactFields.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" +#include "ExplicitQDRateAndState.hpp" + namespace geos { @@ -34,57 +35,71 @@ using namespace fields; using namespace constitutive; using namespace rateAndStateKernels; -QuasiDynamicEarthQuake::QuasiDynamicEarthQuake( const string & name, - Group * const parent ): - ImplicitQDRateAndState( name, parent ), - m_stressSolver( nullptr ), - m_stressSolverName() +template< typename RSSOLVER_TYPE > +QuasiDynamicEarthQuake< RSSOLVER_TYPE >::QuasiDynamicEarthQuake( const string & name, + Group * const parent ): + RSSOLVER_TYPE( name, parent ), + m_stressSolverName(), + m_stressSolver( nullptr ) { this->registerWrapper( viewKeyStruct::stressSolverNameString(), &m_stressSolverName ). setInputFlag( InputFlags::REQUIRED ). setDescription( "Name of solver for computing stress." ); } -void QuasiDynamicEarthQuake::postInputInitialization() +template< typename RSSOLVER_TYPE > +void QuasiDynamicEarthQuake< RSSOLVER_TYPE >::postInputInitialization() { // Initialize member stress solver as specified in XML input - m_stressSolver = &this->getParent().getGroup< SolverBase >( m_stressSolverName ); + m_stressSolver = &this->getParent().template getGroup< PhysicsSolverBase >( m_stressSolverName ); PhysicsSolverBase::postInputInitialization(); } -QuasiDynamicEarthQuake::~QuasiDynamicEarthQuake() +template< typename RSSOLVER_TYPE > +QuasiDynamicEarthQuake< RSSOLVER_TYPE >::~QuasiDynamicEarthQuake() { // TODO Auto-generated destructor stub } -real64 QuasiDynamicEarthQuake::updateStresses( real64 const & time_n, - real64 const & dt, - const int cycleNumber, - DomainPartition & domain ) const +template< typename RSSOLVER_TYPE > +real64 QuasiDynamicEarthQuake< RSSOLVER_TYPE >::updateStresses( real64 const & time_n, + real64 const & dt, + const int cycleNumber, + DomainPartition & domain ) const { // 1. Solve the momentum balance - return m_stressSolver->solverStep( time_n, dt, cycleNumber, domain ); + real64 const dtAccepted = m_stressSolver->solverStep( time_n, dt, cycleNumber, domain ); + // 2. Add background stress and possible forcing. + return dtAccepted; } -void QuasiDynamicEarthQuake::saveOldStateAndUpdateSlip( ElementSubRegionBase & subRegion, real64 const dt ) const +// void QuasiDynamicEarthQuake::updateSlip( ElementSubRegionBase & subRegion, real64 const dt ) const +// { +// arrayView2d< real64 const > const slipVelocity = subRegion.getField< rateAndState::slipVelocity >(); +// arrayView2d< real64 > const deltaSlip = subRegion.getField< contact::deltaSlip >(); + +// arrayView2d< real64 > const dispJump = subRegion.getField< contact::targetIncrementalJump >(); + +// forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) +// { +// deltaSlip[k][0] = slipVelocity[k][0] * dt; +// deltaSlip[k][1] = slipVelocity[k][1] * dt; +// // Update tangential components of the displacement jump +// dispJump[k][1] = slipVelocity[k][0] * dt; +// dispJump[k][2] = slipVelocity[k][1] * dt; +// } ); +// } +template class QuasiDynamicEarthQuake< ImplicitQDRateAndState >; +template class QuasiDynamicEarthQuake< ExplicitQDRateAndState >; + +namespace { - arrayView2d< real64 const > const slipVelocity = subRegion.getField< rateAndState::slipVelocity >(); - arrayView2d< real64 > const deltaSlip = subRegion.getField< contact::deltaSlip >(); - - arrayView2d< real64 > const dispJump = subRegion.getField< contact::targetIncrementalJump >(); - - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) - { - deltaSlip[k][0] = slipVelocity[k][0] * dt; - deltaSlip[k][1] = slipVelocity[k][1] * dt; - // Update tangential components of the displacement jump - dispJump[k][1] = slipVelocity[k][0] * dt; - dispJump[k][2] = slipVelocity[k][1] * dt; - } ); +typedef QuasiDynamicEarthQuake< ImplicitQDRateAndState > ImplicitQuasiDynamicEarthQuake; +typedef QuasiDynamicEarthQuake< ExplicitQDRateAndState > ExplicitQuasiDynamicEarthQuake; +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, ImplicitQuasiDynamicEarthQuake, string const &, dataRepository::Group * const ) +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, ExplicitQuasiDynamicEarthQuake, string const &, dataRepository::Group * const ) } -REGISTER_CATALOG_ENTRY( SolverBase, QuasiDynamicEarthQuake, string const &, dataRepository::Group * const ) - } // namespace geos \ No newline at end of file diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.hpp index 7bb80f2ed3..e11ec9f972 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.hpp @@ -18,13 +18,13 @@ #ifndef GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEARTHQUAKE_HPP #define GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEARTHQUAKE_HPP -#include "physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.hpp" +#include "ImplicitQDRateAndState.hpp" namespace geos { -template< typename QDEQBASE > -class QuasiDynamicEarthQuake : public QDEQBASE +template< typename RSSOLVER_TYPE = ImplicitQDRateAndState > +class QuasiDynamicEarthQuake : public RSSOLVER_TYPE { public: @@ -33,26 +33,35 @@ class QuasiDynamicEarthQuake : public QDEQBASE /// The constructor needs a user-defined "name" and a parent Group (to place this instance in the tree structure of classes) QuasiDynamicEarthQuake( const string & name, - Group * const parent ); + dataRepository::Group * const parent ); /// Destructor - virtual ~QausiDynamicEarthQuake() override; + virtual ~QuasiDynamicEarthQuake() override; - static string catalogName() { return "QuasiDynamicEarthQuake"; } + static string catalogName() { return RSSOLVER_TYPE::derivedSolverPrefix() + "QuasiDynamicEQ"; } /** * @return Get the final class Catalog name */ virtual string getCatalogName() const override { return catalogName(); } - struct viewKeyStruct : public ImplicitQDRateAndState::viewKeyStruct + struct viewKeyStruct : public RSSOLVER_TYPE::viewKeyStruct { /// stress solver name - static constexpr char const * contactSolverNameString() { return "stressSolverName"; } + static constexpr char const * stressSolverNameString() { return "stressSolverName"; } }; + void postInputInitialization() override final; + private: +virtual real64 updateStresses( real64 const & time_n, + real64 const & dt, + const int cycleNumber, + DomainPartition & domain ) const override final; + +string m_stressSolverName; + PhysicsSolverBase * m_stressSolver; }; diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.cpp index 8f9241c3c7..1c0235a1a0 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.cpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.cpp @@ -21,10 +21,11 @@ #include "dataRepository/InputFlags.hpp" #include "mesh/DomainPartition.hpp" -#include "kernels/RateAndStateKernels.hpp" #include "rateAndStateFields.hpp" #include "physicsSolvers/contact/ContactFields.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" +#include "constitutive/contact/RateAndStateFriction.hpp" +#include "ExplicitQDRateAndState.hpp" namespace geos { @@ -65,12 +66,23 @@ void SpringSlider< RSSOLVER_TYPE >::registerDataOnMesh( Group & meshBodies ) subRegion.registerField< contact::dispJump >( this->getName() ). setDimLabels( 1, labels3Comp ). reference().template resizeDimension< 1 >( 3 ); + subRegion.registerField< contact::dispJump_n >( this->getName() ). + setDimLabels( 1, labels3Comp ). + reference().template resizeDimension< 1 >( 3 ); subRegion.registerField< contact::traction >( this->getName() ). setDimLabels( 1, labels3Comp ). reference().template resizeDimension< 1 >( 3 ); + subRegion.registerField< contact::traction_n >( this->getName() ). + setDimLabels( 1, labels3Comp ). + reference().template resizeDimension< 1 >( 3 ); + subRegion.registerField< contact::deltaSlip >( this->getName() ). setDimLabels( 1, labels3Comp ). - reference().template resizeDimension< 1 >( 3 ); + reference().template resizeDimension< 1 >( 3 ); + + subRegion.registerField< contact::deltaSlip_n >( this->getName() ). + setDimLabels( 1, labels3Comp ). + reference().template resizeDimension< 1 >( 3 ); subRegion.registerWrapper< string >( viewKeyStruct::frictionLawNameString() ). setPlotLevel( PlotLevel::NOPLOT ). @@ -86,9 +98,13 @@ void SpringSlider< RSSOLVER_TYPE >::registerDataOnMesh( Group & meshBodies ) } template< typename RSSOLVER_TYPE > -real64 SpringSlider< RSSOLVER_TYPE >::updateStresses( real64 const dt, - DomainPartition & domain ) const +real64 SpringSlider< RSSOLVER_TYPE >::updateStresses( real64 const & time_n, + real64 const & dt, + const int cycleNumber, + DomainPartition & domain ) const + { + GEOS_UNUSED_VAR( cycleNumber, time_n); // Spring-slider shear traction computation this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, @@ -102,6 +118,8 @@ real64 SpringSlider< RSSOLVER_TYPE >::updateStresses( real64 const dt, arrayView2d< real64 const > const deltaSlip = subRegion.getField< fields::contact::deltaSlip >(); arrayView2d< real64 > const traction = subRegion.getField< fields::contact::traction >(); + arrayView2d< real64 > const traction_n = subRegion.getField< fields::contact::traction_n >(); + string const & fricitonLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() ); RateAndStateFriction const & frictionLaw = this->template getConstitutiveModel< RateAndStateFriction >( subRegion, fricitonLawName ); @@ -116,10 +134,13 @@ real64 SpringSlider< RSSOLVER_TYPE >::updateStresses( real64 const dt, frictionKernelWrapper.getDcCoefficient( k ) ); - traction[k][1] = traction[k][1] + springSliderParameters.tauRate * dt - - springSliderParameters.springStiffness * deltaSlip[k][0]; - traction[k][2] = traction[k][2] + springSliderParameters.tauRate * dt - - springSliderParameters.springStiffness * deltaSlip[k][1]; + + + + traction[k][1] = traction_n[k][1] + springSliderParameters.tauRate * dt + - springSliderParameters.springStiffness * deltaSlip[k][0]; + traction[k][2] = traction_n[k][2] + springSliderParameters.tauRate * dt + - springSliderParameters.springStiffness * deltaSlip[k][1]; } ); } ); } ); @@ -127,11 +148,14 @@ real64 SpringSlider< RSSOLVER_TYPE >::updateStresses( real64 const dt, } template class SpringSlider< ImplicitQDRateAndState >; +template class SpringSlider< ExplicitQDRateAndState >; namespace { typedef SpringSlider< ImplicitQDRateAndState > ImplicitSpringSlider; +typedef SpringSlider< ExplicitQDRateAndState > ExplicitSpringSlider; REGISTER_CATALOG_ENTRY( PhysicsSolverBase, ImplicitSpringSlider, string const &, dataRepository::Group * const ) +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, ExplicitSpringSlider, string const &, dataRepository::Group * const ) } } // namespace geos diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.hpp index c8edd92c98..d2d5ea5ba0 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.hpp @@ -47,8 +47,10 @@ class SpringSlider : public RSSOLVER_TYPE private: - real64 updateStresses( real64 const dt, - DomainPartition & domain ) const override final; + virtual real64 updateStresses( real64 const & time_n, + real64 const & dt, + const int cycleNumber, + DomainPartition & domain ) const override final; class SpringSliderParameters { diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernels.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/ExplicitRateAndStateKernels.hpp similarity index 61% rename from src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernels.hpp rename to src/coreComponents/physicsSolvers/inducedSeismicity/kernels/ExplicitRateAndStateKernels.hpp index fdd58d595e..c336238f29 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernels.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/ExplicitRateAndStateKernels.hpp @@ -13,13 +13,10 @@ * ------------------------------------------------------------------------------------------------------------ */ -#ifndef GEOS_PHYSICSSOLVERS_RATEANDSTATEKERNELS_HPP_ -#define GEOS_PHYSICSSOLVERS_RATEANDSTATEKERNELS_HPP_ +#ifndef GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_EXPLICITRATEANDSTATEKERNELS_HPP_ +#define GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_EXPLICITRATEANDSTATEKERNELS_HPP_ -#include "common/DataTypes.hpp" -#include "common/GEOS_RAJA_Interface.hpp" -#include "constitutive/contact/RateAndStateFriction.hpp" -#include "physicsSolvers/inducedSeismicity/rateAndStateFields.hpp" +#include "RateAndStateKernelsBase.hpp" #include "denseLinearAlgebra/denseLASolvers.hpp" namespace geos @@ -28,175 +25,6 @@ namespace geos namespace rateAndStateKernels { -// TBD: Pass the kernel and add getters for relevant fields to make this function general purpose and avoid -// wrappers? -GEOS_HOST_DEVICE -static void projectSlipRateBase( localIndex const k, - real64 const frictionCoefficient, - real64 const shearImpedance, - arrayView2d< real64 const > const traction, - arrayView1d< real64 const > const slipRate, - arrayView2d< real64 > const slipVelocity ) -{ - // Project slip rate onto shear traction to get slip velocity components - real64 const frictionForce = traction[k][0] * frictionCoefficient; - real64 const projectionScaling = 1.0 / ( shearImpedance + frictionForce / slipRate[k] ); - slipVelocity[k][0] = projectionScaling * traction[k][1]; - slipVelocity[k][1] = projectionScaling * traction[k][2]; -} - -/** - * @class ImplicitFixedStressRateAndStateKernel - * - * @brief - * - * @details - */ -class ImplicitFixedStressRateAndStateKernel -{ -public: - - ImplicitFixedStressRateAndStateKernel( SurfaceElementSubRegion & subRegion, - constitutive::RateAndStateFriction const & frictionLaw, - real64 const shearImpedance ): - m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ), - m_stateVariable( subRegion.getField< fields::rateAndState::stateVariable >() ), - m_stateVariable_n( subRegion.getField< fields::rateAndState::stateVariable_n >() ), - m_slipRate_n( subRegion.getField< fields::rateAndState::slipRate_n >() ), - m_traction( subRegion.getField< fields::contact::traction >() ), - m_slipVelocity( subRegion.getField< fields::rateAndState::slipVelocity >() ), - m_shearImpedance( shearImpedance ), - m_frictionLaw( frictionLaw.createKernelUpdates() ) - {} - - /** - * @struct StackVariables - * @brief Kernel variables located on the stack - */ - struct StackVariables - { -public: - - GEOS_HOST_DEVICE - StackVariables( ) - {} - - real64 jacobian[2][2]{}; - - real64 rhs[2]{}; - - }; - - GEOS_HOST_DEVICE - void setup( localIndex const k, - real64 const dt, - StackVariables & stack ) const - { - real64 const normalTraction = m_traction[k][0]; - real64 const shearTractionMagnitude = LvArray::math::sqrt( m_traction[k][1] * m_traction[k][1] + m_traction[k][2] * m_traction[k][2] ); - - // Eq 1: Scalar force balance for slipRate and shear traction magnitude - stack.rhs[0] = shearTractionMagnitude - m_shearImpedance * m_slipRate[k] - - normalTraction * m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] ); - real64 const dFriction[2] = { -normalTraction * m_frictionLaw.dFrictionCoefficient_dStateVariable( k, m_slipRate[k], m_stateVariable[k] ), - -m_shearImpedance - normalTraction * m_frictionLaw.dFrictionCoefficient_dSlipRate( k, m_slipRate[k], m_stateVariable[k] ) }; - - // Eq 2: slip law - stack.rhs[1] = (m_stateVariable[k] - m_stateVariable_n[k]) / dt - m_frictionLaw.stateEvolution( k, m_slipRate[k], m_stateVariable[k] ); - real64 const dStateEvolutionLaw[2] = { 1.0 / dt - m_frictionLaw.dStateEvolution_dStateVariable( k, m_slipRate[k], m_stateVariable[k] ), - -m_frictionLaw.dStateEvolution_dSlipRate( k, m_slipRate[k], m_stateVariable[k] ) }; - - // Assemble Jacobian matrix - stack.jacobian[0][0] = dFriction[0]; // derivative of Eq 1 w.r.t. stateVariable - stack.jacobian[0][1] = dFriction[1]; // derivative of Eq 1 w.r.t. slipRate - stack.jacobian[1][0] = dStateEvolutionLaw[0]; // derivative of Eq 2 w.r.t. stateVariable - stack.jacobian[1][1] = dStateEvolutionLaw[1]; // derivative of Eq 2 w.r.t. slipRate - } - - GEOS_HOST_DEVICE - void solve( localIndex const k, - StackVariables & stack ) const - { - /// Solve 2x2 system - real64 solution[2] = {0.0, 0.0}; - LvArray::tensorOps::scale< 2 >( stack.rhs, -1.0 ); - denseLinearAlgebra::solve< 2 >( stack.jacobian, stack.rhs, solution ); - - // Slip rate is bracketed between [0, shear traction magnitude / shear impedance] - // Check that the update did not end outside of the bracket. - real64 const shearTractionMagnitude = LvArray::math::sqrt( m_traction[k][1] * m_traction[k][1] + m_traction[k][2] * m_traction[k][2] ); - real64 const upperBound = shearTractionMagnitude / m_shearImpedance - m_slipRate[k]; - real64 const lowerBound = -m_slipRate[k]; - - real64 scalingFactor = 1.0; - if( solution[1] > upperBound ) - { - scalingFactor = 0.5 * upperBound / solution[1]; - } - else if( solution[1] < lowerBound ) - { - scalingFactor = 0.5 * lowerBound / solution[1]; - } - - LvArray::tensorOps::scale< 2 >( solution, scalingFactor ); - - // Update variables - m_stateVariable[k] += solution[0]; - m_slipRate[k] += solution[1]; - } - - GEOS_HOST_DEVICE - void projectSlipRate( localIndex const k ) const - { - real64 const frictionCoefficient = m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] ); - projectSlipRateBase( k, frictionCoefficient, m_shearImpedance, m_traction, m_slipRate, m_slipVelocity ); - } - - GEOS_HOST_DEVICE - void udpateVariables( localIndex const k ) const - { - projectSlipRate( k ); - m_stateVariable_n[k] = m_stateVariable[k]; - m_slipRate_n[k] = m_slipRate[k]; - } - - GEOS_HOST_DEVICE - camp::tuple< int, real64 > checkConvergence( StackVariables const & stack, - real64 const tol ) const - { - real64 const residualNorm = LvArray::tensorOps::l2Norm< 2 >( stack.rhs ); - int const converged = residualNorm < tol ? 1 : 0; - camp::tuple< int, real64 > result { converged, residualNorm }; - return result; - } - - GEOS_HOST_DEVICE - void resetState( localIndex const k ) const - { - m_stateVariable[k] = m_stateVariable_n[k]; - m_slipRate[k] = m_slipRate_n[k]; - } - -private: - - arrayView1d< real64 > const m_slipRate; - - arrayView1d< real64 > const m_stateVariable; - - arrayView1d< real64 > const m_stateVariable_n; - - arrayView1d< real64 > const m_slipRate_n; - - arrayView2d< real64 const > const m_traction; - - arrayView2d< real64 > const m_slipVelocity; - - real64 const m_shearImpedance; - - constitutive::RateAndStateFriction::KernelWrapper m_frictionLaw; - -}; - /** * @class ExplicitRateAndStateKernel * @@ -286,132 +114,57 @@ class ExplicitRateAndStateKernel projectSlipRateBase( k, frictionCoefficient, m_shearImpedance, m_traction, m_slipRate, m_slipVelocity ); } -private: - - arrayView1d< real64 > const m_slipRate; - - arrayView1d< real64 > const m_stateVariable; - - arrayView2d< real64 const > const m_traction; - - arrayView2d< real64 > const m_slipVelocity; - - real64 const m_shearImpedance; - - constitutive::RateAndStateFriction::KernelWrapper m_frictionLaw; - -}; - -template< typename POLICY > -static bool newtonSolve( SurfaceElementSubRegion & subRegion, - RateAndStateKernel & kernel, - real64 const dt, - integer const maxNewtonIter - real64 const newtonTol) -{ - bool allConverged = false; - for( integer iter = 0; iter < maxIterNewton; iter++ ) + GEOS_HOST_DEVICE + void udpateVariables( localIndex const k ) const { - RAJA::ReduceMin< parallelDeviceReduce, int > converged( 1 ); - RAJA::ReduceMax< parallelDeviceReduce, real64 > residualNorm( 0.0 ); - forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) - { - typename KernelType::StackVariables stack; - kernel.setup( k, dt, stack ); - kernel.solve( k, stack ); - auto const [elementConverged, elementResidualNorm] = kernel.checkConvergence( stack, newtonTol ); - converged.min( elementConverged ); - residualNorm.max( elementResidualNorm ); - } ); - - real64 const maxResidualNorm = MpiWrapper::max( residualNorm.get() ); - GEOS_LOG_RANK_0( GEOS_FMT( " Newton iter {} : residual = {:.10e} ", iter, maxResidualNorm ) ); - - if( converged.get() ) - { - allConverged = true; - break; - } + projectSlipRate( k ); } - return allConverged; -} -template< typename POLICY > -static real64 solveRateAndStateEquation( SurfaceElementSubRegion & subRegion, - RateAndStateKernel & kernel, - real64 dt, - integer const maxNewtonIter - real64 const newtonTol ) -{ - bool converged = false; - for( integer attempt = 0; attempt < 5; attempt++ ) - { - if( attempt > 0 ) - { - forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) - { - kernel.resetState( k ); - } ); - } - GEOS_LOG_RANK_0( GEOS_FMT( " Attempt {} ", attempt ) ); - converged = newtonSolve< POLICY >( subRegion, kernel, dt, maxNewtonIter, newtonTol ); - if( converged ) - { - forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) - { - kernel.udpateVariables( k ); - } ); - return dt; - } - else - { - GEOS_LOG_RANK_0( GEOS_FMT( " Attempt {} failed. Halving dt and retrying.", attempt ) ); - dt *= 0.5; - } - } - if( !converged ) + GEOS_HOST_DEVICE + void resetState( localIndex const k ) const { - GEOS_ERROR( "Maximum number of attempts reached without convergence." ); + GEOS_UNUSED_VAR( k ); } - return dt; -} -/** + /** * @brief Performs the kernel launch + * @tparam KernelType The Rate-and-state kernel to launch * @tparam POLICY the policy used in the RAJA kernels */ template< typename POLICY > -static void -createAndLaunch( SurfaceElementSubRegion & subRegion, - string const & frictionLawNameKey, - real64 const shearImpedance, - integer const maxNewtonIter, - real64 const newtonTol, - real64 const time_n, - real64 const totalDt ) +static real64 +solveRateAndStateEquation( SurfaceElementSubRegion & subRegion, + ExplicitRateAndStateKernel & kernel, + real64 dt, + integer const maxNewtonIter, + real64 const newtonTol ) { GEOS_MARK_FUNCTION; + + newtonSolve< POLICY >( subRegion, kernel, dt, maxNewtonIter, newtonTol ); - GEOS_UNUSED_VAR( time_n ); + forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + kernel.projectSlipRate( k ); + } ); + return dt; +} - string const & frictionaLawName = subRegion.getReference< string >( frictionLawNameKey ); - constitutive::RateAndStateFriction const & frictionLaw = subRegion.getConstitutiveModel< constitutive::RateAndStateFriction >( frictionaLawName ); - RateAndStateKernel kernel( subRegion, frictionLaw, shearImpedance ); +private: - real64 dtRemaining = totalDt; - real64 dt = totalDt; - for( integer subStep = 0; subStep < 5 && dtRemaining > 0.0; ++subStep ) - { - real64 dtAccepted = solveRateAndStateEquation< POLICY >( subRegion, kernel, dt, maxNewtonIter, newtonTol ); - dtRemaining -= dtAccepted; + arrayView1d< real64 > const m_slipRate; - if( dtRemaining > 0.0 ) - { - dt = dtAccepted; - } - GEOS_LOG_RANK_0( GEOS_FMT( " sub-step = {} completed, dt = {}, remaining dt = {}", subStep, dt, dtRemaining ) ); - } -} + arrayView1d< real64 > const m_stateVariable; + + arrayView2d< real64 const > const m_traction; + + arrayView2d< real64 > const m_slipVelocity; + + real64 const m_shearImpedance; + + constitutive::RateAndStateFriction::KernelWrapper m_frictionLaw; + +}; /** * @brief Butcher table for embedded RK3(2) method using Kuttas third order @@ -469,8 +222,8 @@ class EmbeddedRungeKuttaKernel m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ), m_slipVelocity( subRegion.getField< fields::rateAndState::slipVelocity >() ), m_slipVelocity_n( subRegion.getField< fields::rateAndState::slipVelocity_n >() ), - m_deltaSlip( subRegion.getField< fields::rateAndState::deltaSlip >() ), - m_deltaSlip_n( subRegion.getField< fields::rateAndState::deltaSlip_n >() ), + m_deltaSlip( subRegion.getField< fields::contact::deltaSlip >() ), + m_deltaSlip_n( subRegion.getField< fields::contact::deltaSlip_n >() ), m_dispJump( subRegion.getField< fields::contact::dispJump >() ), m_dispJump_n( subRegion.getField< fields::contact::dispJump_n >() ), m_error( subRegion.getField< fields::rateAndState::error >() ), @@ -670,4 +423,4 @@ class EmbeddedRungeKuttaKernel } /* namespace geos */ -#endif /* GEOS_PHYSICSSOLVERS_RATEANDSTATEKERNELS_HPP_ */ +#endif /* GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_EXPLICITRATEANDSTATEKERNELS_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/ImplicitRateAndStateKernels.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/ImplicitRateAndStateKernels.hpp new file mode 100644 index 0000000000..cb1922def2 --- /dev/null +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/ImplicitRateAndStateKernels.hpp @@ -0,0 +1,224 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-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_INDUCEDSEISMICITY_KERNELS_IMPLICITRATEANDSTATEKERNELS_HPP_ +#define GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_IMPLICITRATEANDSTATEKERNELS_HPP_ + +#include "RateAndStateKernelsBase.hpp" +#include "denseLinearAlgebra/denseLASolvers.hpp" + +namespace geos +{ + +namespace rateAndStateKernels +{ + +/** + * @class ImplicitFixedStressRateAndStateKernel + * + * @brief + * + * @details + */ +class ImplicitFixedStressRateAndStateKernel +{ +public: + + ImplicitFixedStressRateAndStateKernel( SurfaceElementSubRegion & subRegion, + constitutive::RateAndStateFriction const & frictionLaw, + real64 const shearImpedance ): + m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ), + m_stateVariable( subRegion.getField< fields::rateAndState::stateVariable >() ), + m_stateVariable_n( subRegion.getField< fields::rateAndState::stateVariable_n >() ), + m_slipRate_n( subRegion.getField< fields::rateAndState::slipRate_n >() ), + m_traction( subRegion.getField< fields::contact::traction >() ), + m_slipVelocity( subRegion.getField< fields::rateAndState::slipVelocity >() ), + m_shearImpedance( shearImpedance ), + m_frictionLaw( frictionLaw.createKernelUpdates() ) + {} + + /** + * @struct StackVariables + * @brief Kernel variables located on the stack + */ + struct StackVariables + { +public: + + GEOS_HOST_DEVICE + StackVariables( ) + {} + + real64 jacobian[2][2]{}; + + real64 rhs[2]{}; + + }; + + GEOS_HOST_DEVICE + void setup( localIndex const k, + real64 const dt, + StackVariables & stack ) const + { + real64 const normalTraction = m_traction[k][0]; + real64 const shearTractionMagnitude = LvArray::math::sqrt( m_traction[k][1] * m_traction[k][1] + m_traction[k][2] * m_traction[k][2] ); + + // Eq 1: Scalar force balance for slipRate and shear traction magnitude + stack.rhs[0] = shearTractionMagnitude - m_shearImpedance * m_slipRate[k] + - normalTraction * m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] ); + real64 const dFriction[2] = { -normalTraction * m_frictionLaw.dFrictionCoefficient_dStateVariable( k, m_slipRate[k], m_stateVariable[k] ), + -m_shearImpedance - normalTraction * m_frictionLaw.dFrictionCoefficient_dSlipRate( k, m_slipRate[k], m_stateVariable[k] ) }; + + // Eq 2: slip law + stack.rhs[1] = (m_stateVariable[k] - m_stateVariable_n[k]) / dt - m_frictionLaw.stateEvolution( k, m_slipRate[k], m_stateVariable[k] ); + real64 const dStateEvolutionLaw[2] = { 1.0 / dt - m_frictionLaw.dStateEvolution_dStateVariable( k, m_slipRate[k], m_stateVariable[k] ), + -m_frictionLaw.dStateEvolution_dSlipRate( k, m_slipRate[k], m_stateVariable[k] ) }; + + // Assemble Jacobian matrix + stack.jacobian[0][0] = dFriction[0]; // derivative of Eq 1 w.r.t. stateVariable + stack.jacobian[0][1] = dFriction[1]; // derivative of Eq 1 w.r.t. slipRate + stack.jacobian[1][0] = dStateEvolutionLaw[0]; // derivative of Eq 2 w.r.t. stateVariable + stack.jacobian[1][1] = dStateEvolutionLaw[1]; // derivative of Eq 2 w.r.t. slipRate + } + + GEOS_HOST_DEVICE + void solve( localIndex const k, + StackVariables & stack ) const + { + /// Solve 2x2 system + real64 solution[2] = {0.0, 0.0}; + LvArray::tensorOps::scale< 2 >( stack.rhs, -1.0 ); + denseLinearAlgebra::solve< 2 >( stack.jacobian, stack.rhs, solution ); + + // Slip rate is bracketed between [0, shear traction magnitude / shear impedance] + // Check that the update did not end outside of the bracket. + real64 const shearTractionMagnitude = LvArray::math::sqrt( m_traction[k][1] * m_traction[k][1] + m_traction[k][2] * m_traction[k][2] ); + real64 const upperBound = shearTractionMagnitude / m_shearImpedance - m_slipRate[k]; + real64 const lowerBound = -m_slipRate[k]; + + real64 scalingFactor = 1.0; + if( solution[1] > upperBound ) + { + scalingFactor = 0.5 * upperBound / solution[1]; + } + else if( solution[1] < lowerBound ) + { + scalingFactor = 0.5 * lowerBound / solution[1]; + } + + LvArray::tensorOps::scale< 2 >( solution, scalingFactor ); + + // Update variables + m_stateVariable[k] += solution[0]; + m_slipRate[k] += solution[1]; + } + + GEOS_HOST_DEVICE + void projectSlipRate( localIndex const k ) const + { + real64 const frictionCoefficient = m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] ); + projectSlipRateBase( k, frictionCoefficient, m_shearImpedance, m_traction, m_slipRate, m_slipVelocity ); + } + + GEOS_HOST_DEVICE + void udpateVariables( localIndex const k ) const + { + projectSlipRate( k ); + m_stateVariable_n[k] = m_stateVariable[k]; + m_slipRate_n[k] = m_slipRate[k]; + } + + GEOS_HOST_DEVICE + camp::tuple< int, real64 > checkConvergence( StackVariables const & stack, + real64 const tol ) const + { + real64 const residualNorm = LvArray::tensorOps::l2Norm< 2 >( stack.rhs ); + int const converged = residualNorm < tol ? 1 : 0; + camp::tuple< int, real64 > result { converged, residualNorm }; + return result; + } + + GEOS_HOST_DEVICE + void resetState( localIndex const k ) const + { + m_stateVariable[k] = m_stateVariable_n[k]; + m_slipRate[k] = m_slipRate_n[k]; + } + +template< typename POLICY > +static real64 solveRateAndStateEquation( SurfaceElementSubRegion & subRegion, + ImplicitFixedStressRateAndStateKernel & kernel, + real64 dt, + integer const maxNewtonIter, + real64 const newtonTol ) +{ + bool converged = false; + for( integer attempt = 0; attempt < 5; attempt++ ) + { + if( attempt > 0 ) + { + forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + kernel.resetState( k ); + } ); + } + GEOS_LOG_RANK_0( GEOS_FMT( " Attempt {} ", attempt ) ); + converged = newtonSolve< POLICY >( subRegion, kernel, dt, maxNewtonIter, newtonTol ); + if( converged ) + { + forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + kernel.udpateVariables( k ); + } ); + return dt; + } + else + { + GEOS_LOG_RANK_0( GEOS_FMT( " Attempt {} failed. Halving dt and retrying.", attempt ) ); + dt *= 0.5; + } + } + if( !converged ) + { + GEOS_ERROR( "Maximum number of attempts reached without convergence." ); + } + return dt; +} + +private: + + arrayView1d< real64 > const m_slipRate; + + arrayView1d< real64 > const m_stateVariable; + + arrayView1d< real64 > const m_stateVariable_n; + + arrayView1d< real64 > const m_slipRate_n; + + arrayView2d< real64 const > const m_traction; + + arrayView2d< real64 > const m_slipVelocity; + + real64 const m_shearImpedance; + + constitutive::RateAndStateFriction::KernelWrapper m_frictionLaw; + +}; + +} /* namespace rateAndStateKernels */ + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_IMPLICITRATEANDSTATEKERNELS_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernelsBase.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernelsBase.hpp new file mode 100644 index 0000000000..62ef835346 --- /dev/null +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernelsBase.hpp @@ -0,0 +1,162 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-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_INDUCEDSEISMICITY_KERNELS_RATEANDSTATEKERNELSBASE_HPP_ +#define GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_RATEANDSTATEKERNELSBASE_HPP_ + +#include "common/DataTypes.hpp" +#include "common/GEOS_RAJA_Interface.hpp" +#include "constitutive/contact/RateAndStateFriction.hpp" +#include "physicsSolvers/inducedSeismicity/rateAndStateFields.hpp" + +namespace geos +{ + +namespace rateAndStateKernels +{ + +// TBD: Pass the kernel and add getters for relevant fields to make this function general purpose and avoid +// wrappers? +GEOS_HOST_DEVICE +static void projectSlipRateBase( localIndex const k, + real64 const frictionCoefficient, + real64 const shearImpedance, + arrayView2d< real64 const > const traction, + arrayView1d< real64 const > const slipRate, + arrayView2d< real64 > const slipVelocity ) +{ + // Project slip rate onto shear traction to get slip velocity components + real64 const frictionForce = traction[k][0] * frictionCoefficient; + real64 const projectionScaling = 1.0 / ( shearImpedance + frictionForce / slipRate[k] ); + slipVelocity[k][0] = projectionScaling * traction[k][1]; + slipVelocity[k][1] = projectionScaling * traction[k][2]; +} + +template< typename POLICY, typename KERNEL_TYPE > +static bool newtonSolve( SurfaceElementSubRegion & subRegion, + KERNEL_TYPE & kernel, + real64 const dt, + integer const maxNewtonIter, + real64 const newtonTol ) +{ + bool allConverged = false; + for( integer iter = 0; iter < maxNewtonIter; iter++ ) + { + RAJA::ReduceMin< parallelDeviceReduce, int > converged( 1 ); + RAJA::ReduceMax< parallelDeviceReduce, real64 > residualNorm( 0.0 ); + forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + typename KERNEL_TYPE::StackVariables stack; + kernel.setup( k, dt, stack ); + kernel.solve( k, stack ); + auto const [elementConverged, elementResidualNorm] = kernel.checkConvergence( stack, newtonTol ); + converged.min( elementConverged ); + residualNorm.max( elementResidualNorm ); + } ); + + real64 const maxResidualNorm = MpiWrapper::max( residualNorm.get() ); + GEOS_LOG_RANK_0( GEOS_FMT( " Newton iter {} : residual = {:.10e} ", iter, maxResidualNorm ) ); + + if( converged.get() ) + { + allConverged = true; + break; + } + } + return allConverged; +} + +// template< typename POLICY, typename KERNEL_TYPE > +// static real64 solveRateAndStateEquation( SurfaceElementSubRegion & subRegion, +// KERNEL_TYPE & kernel, +// real64 dt, +// integer const maxNewtonIter, +// real64 const newtonTol ) +// { +// bool converged = false; +// for( integer attempt = 0; attempt < 5; attempt++ ) +// { +// if( attempt > 0 ) +// { +// forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) +// { +// kernel.resetState( k ); +// } ); +// } +// GEOS_LOG_RANK_0( GEOS_FMT( " Attempt {} ", attempt ) ); +// converged = newtonSolve< POLICY >( subRegion, kernel, dt, maxNewtonIter, newtonTol ); +// if( converged ) +// { +// forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) +// { +// kernel.udpateVariables( k ); +// } ); +// return dt; +// } +// else +// { +// GEOS_LOG_RANK_0( GEOS_FMT( " Attempt {} failed. Halving dt and retrying.", attempt ) ); +// dt *= 0.5; +// } +// } +// if( !converged ) +// { +// GEOS_ERROR( "Maximum number of attempts reached without convergence." ); +// } +// return dt; +// } + +/** + * @brief Performs the kernel launch + * @tparam POLICY the policy used in the RAJA kernels + */ +template< typename KERNEL_TYPE, typename POLICY > +static void +createAndLaunch( SurfaceElementSubRegion & subRegion, + string const & frictionLawNameKey, + real64 const shearImpedance, + integer const maxNewtonIter, + real64 const newtonTol, + real64 const time_n, + real64 const totalDt ) +{ + GEOS_MARK_FUNCTION; + + GEOS_UNUSED_VAR( time_n ); + + string const & frictionaLawName = subRegion.getReference< string >( frictionLawNameKey ); + constitutive::RateAndStateFriction const & frictionLaw = subRegion.getConstitutiveModel< constitutive::RateAndStateFriction >( frictionaLawName ); + KERNEL_TYPE kernel( subRegion, frictionLaw, shearImpedance ); + + real64 dtRemaining = totalDt; + real64 dt = totalDt; + for( integer subStep = 0; subStep < 5 && dtRemaining > 0.0; ++subStep ) + { + real64 dtAccepted = KERNEL_TYPE::template solveRateAndStateEquation< POLICY >( subRegion, kernel, dt, maxNewtonIter, newtonTol ); + dtRemaining -= dtAccepted; + + if( dtRemaining > 0.0 ) + { + dt = dtAccepted; + } + GEOS_LOG_RANK_0( GEOS_FMT( " sub-step = {} completed, dt = {}, remaining dt = {}", subStep, dt, dtRemaining ) ); + } +} + +} /* namespace rateAndStateKernels */ + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_RATEANDSTATEKERNELSBASE_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp index 3a16ff2b72..d271c851ca 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp @@ -86,16 +86,7 @@ DECLARE_FIELD( shearTraction, 29.2e6, LEVEL_0, WRITE_AND_READ, - "shear Traction" ); - -DECLARE_FIELD( deltaSlip_n, - "deltaSlip_n", - array2d< real64 >, - 0.0, - NOPLOT, - WRITE_AND_READ, - "Initial slip increment at this time step" ); - + "shear Traction" ); DECLARE_FIELD( rungeKuttaStageRates, "rungeKuttaStageRates", diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index eed28a4ccd..1c9a232d20 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -403,8 +403,8 @@ - - + + @@ -2296,7 +2296,7 @@ Level 0 outputs no specific information for this solver. Higher levels require m - + @@ -3623,7 +3623,7 @@ Level 0 outputs no specific information for this solver. Higher levels require m ======= - + diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index 0fb2b8d913..7a366ce326 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -539,7 +539,7 @@ - + @@ -1044,7 +1044,7 @@ <<<<<<< HEAD ======= - +