From c6e3171441295fd57fc103038a8b2db1fd237ffd Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Mon, 27 Jan 2025 20:24:33 -0800 Subject: [PATCH 1/9] iterative solver + reference solution in plotting script. --- .../SCEC_BP6_QD_S_implicit.xml | 9 +- .../inducedSeismicity/SCEC_BP7_QD_base.xml | 75 ----------------- .../SCEC_BP7_QD_implicit.xml | 84 +++++++++++++++++++ .../scripts/plotBP6_results.py | 25 +++++- 4 files changed, 112 insertions(+), 81 deletions(-) create mode 100644 inputFiles/inducedSeismicity/SCEC_BP7_QD_implicit.xml diff --git a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml index 87d8a3e8455..1e913da0c96 100644 --- a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml +++ b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml @@ -57,8 +57,8 @@ lineSearchMaxCuts="2" maxTimeStepCuts="1"/> @@ -78,8 +78,9 @@ newtonTol = "1.0e-6" newtonMaxIter ="8"/> + solverType="gmres" + preconditionerType="amg" + logLevel="0"/> diff --git a/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml b/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml index 8df6ff94139..6161fc7caf4 100644 --- a/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml +++ b/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml @@ -1,54 +1,5 @@ - - - - - - - - - - @@ -222,30 +173,4 @@ fieldName="pressure" setNames="{source}"/> --> - - - - - - - - - - - diff --git a/inputFiles/inducedSeismicity/SCEC_BP7_QD_implicit.xml b/inputFiles/inducedSeismicity/SCEC_BP7_QD_implicit.xml new file mode 100644 index 00000000000..9d1569167b6 --- /dev/null +++ b/inputFiles/inducedSeismicity/SCEC_BP7_QD_implicit.xml @@ -0,0 +1,84 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/inducedSeismicity/scripts/plotBP6_results.py b/inputFiles/inducedSeismicity/scripts/plotBP6_results.py index e869c3d5b5a..ee709bf9c53 100644 --- a/inputFiles/inducedSeismicity/scripts/plotBP6_results.py +++ b/inputFiles/inducedSeismicity/scripts/plotBP6_results.py @@ -3,6 +3,7 @@ from geos.hdf5_wrapper import hdf5_wrapper import matplotlib.pyplot as plt import os +from scipy.special import erfc def remove_padding(data): """Removes trailing zeros from a NumPy array.""" @@ -22,9 +23,23 @@ def getDataFromHDF5( hdf5FilePath, var_name ): # Remove padding var = remove_padding(var) time = time[:len(var)] # Match the length of the time array to the variable - return time, var + return time, var +def analytical_pressure( time, x ): + alpha = 0.1 # m^2 / s + phi = 0.1 + beta = 1e-8 # Pa^-1 + t_off = 100 * 24 * 3600 + q0 = 1.25 * 1e-6 #m/s + pressure = q0 / (beta * phi * np.sqrt(alpha)) * (G(time, x, alpha) - np.where(time > t_off, G(time - t_off, x, alpha), 0)) + return pressure + +def G( t, x, alpha ): + A = np.abs(x) / np.sqrt( 4*alpha*t ) + B = -x**2 / (4 * alpha * t) + return np.sqrt(t) * ( np.exp(B) / np.sqrt(np.pi) - A * erfc( A ) ) + if __name__ == "__main__": parser = argparse.ArgumentParser() parser.add_argument('-d', '--dir', type=str, help='Path to hdf5 file') @@ -43,6 +58,9 @@ def getDataFromHDF5( hdf5FilePath, var_name ): time, pressure = getDataFromHDF5( filePath, "pressure" ) time, slipRate = getDataFromHDF5( filePath, "slipRate" ) + pressure_analytical = analytical_pressure( time, 0.0 ) + print(pressure_analytical) + # Convert time to years time_in_years = time / (365 * 24 * 3600) # Assuming time is in seconds @@ -50,13 +68,16 @@ def getDataFromHDF5( hdf5FilePath, var_name ): fig, ax1 = plt.subplots() print(len(time_in_years)) - print(len(pressure)) + print(len(pressure)) + print(len(pressure_analytical)) # Plot pressure on the left y-axis ax1.set_xlabel('Time (years)') ax1.set_ylabel('Pressure (Pa)', color='tab:blue') ax1.plot(time_in_years, pressure, label="Pressure", color='tab:blue') + ax1.plot(time_in_years, pressure_analytical, label="Pressure Analytical", color='black', linestyle='--') ax1.tick_params(axis='y', labelcolor='tab:blue') + ax1.legend() # Plot slipRate on the right y-axis (log scale) From 3e2796d446ac515ef6e45e2489f169063e87df46 Mon Sep 17 00:00:00 2001 From: Matteo Cusini Date: Tue, 28 Jan 2025 16:43:22 -0800 Subject: [PATCH 2/9] remove separate displacement component filter. --- .../LagrangeContactBubbleStab_FixedSlip_smoke.xml | 7 ++++--- .../linearAlgebra/interfaces/hypre/HypreUtils.hpp | 3 ++- .../contact/SolidMechanicsLagrangeContactBubbleStab.cpp | 8 ++++++-- 3 files changed, 12 insertions(+), 6 deletions(-) diff --git a/inputFiles/lagrangianContactMechanics/LagrangeContactBubbleStab_FixedSlip_smoke.xml b/inputFiles/lagrangianContactMechanics/LagrangeContactBubbleStab_FixedSlip_smoke.xml index 5da14675d94..4c68e85d42b 100644 --- a/inputFiles/lagrangianContactMechanics/LagrangeContactBubbleStab_FixedSlip_smoke.xml +++ b/inputFiles/lagrangianContactMechanics/LagrangeContactBubbleStab_FixedSlip_smoke.xml @@ -22,9 +22,10 @@ lineSearchMaxCuts="2" maxTimeStepCuts="1"/> + logLevel="1"/> diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.hpp index 47f901330f0..51169a7fc99 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.hpp @@ -541,7 +541,8 @@ enum class MGRFRelaxationType : HYPRE_Int l1backwardGaussSeidel = 14, //!< \f$\ell_1\f$ Gauss-Seidel, backward solve l1jacobi = 18, //!< \f$\ell_1\f$-scaled Jacobi gsElimWPivoting = 99, //!< Gaussian Elimination with pivoting direct solver (for small systems) - gsElimWInverse = 199 //!< Direct Inversion with Gaussian Elimination (OK for larger systems) + gsElimWInverse = 199, //!< Direct Inversion with Gaussian Elimination (OK for larger systems) + blockJacobi = 0 }; /** diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContactBubbleStab.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContactBubbleStab.cpp index 7c2543a3a08..4e2663a6e2e 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContactBubbleStab.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContactBubbleStab.cpp @@ -48,7 +48,7 @@ SolidMechanicsLagrangeContactBubbleStab::SolidMechanicsLagrangeContactBubbleStab LinearSolverParameters & linSolParams = m_linearSolverParameters.get(); linSolParams.mgr.strategy = LinearSolverParameters::MGR::StrategyType::lagrangianContactMechanicsBubbleStab; - linSolParams.mgr.separateComponents = true; + linSolParams.mgr.separateComponents = false; linSolParams.dofsPerNode = 3; } @@ -896,6 +896,8 @@ void SolidMechanicsLagrangeContactBubbleStab::updateStickSlipList( DomainPartiti arrayView1d< integer const > const fractureState = subRegion.getField< contact::fractureState >(); + // std::cout << fractureState << std::endl; + forFiniteElementOnFractureSubRegions( meshName, [&] ( string const & finiteElementName, finiteElement::FiniteElementBase const &, arrayView1d< localIndex const > const & faceElementList ) @@ -910,6 +912,8 @@ void SolidMechanicsLagrangeContactBubbleStab::updateStickSlipList( DomainPartiti arrayView1d< localIndex > const keys_v = keys.toView(); arrayView1d< localIndex > const vals_v = vals.toView(); + + // std::cout<< faceElementList.size() << std::endl; forAll< parallelDevicePolicy<> >( faceElementList.size(), [ = ] GEOS_HOST_DEVICE ( localIndex const kfe ) @@ -965,7 +969,7 @@ void SolidMechanicsLagrangeContactBubbleStab::updateStickSlipList( DomainPartiti this->m_faceTypesToFaceElementsStick[meshName][finiteElementName] = stickList; this->m_faceTypesToFaceElementsSlip[meshName][finiteElementName] = slipList; - GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Configuration, GEOS_FMT( "# stick elements: {}, # slip elements: {}", nStick, nSlip )) + GEOS_LOG_LEVEL_INFO_BY_RANK( logInfo::Configuration, GEOS_FMT( "# stick elements: {}, # slip elements: {}", nStick, nSlip )) } ); } ); From 40c44dc8c9d03886feb75aeb0fd3da22313865d8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vidar=20Stiernstr=C3=B6m?= Date: Wed, 29 Jan 2025 09:06:55 -0800 Subject: [PATCH 3/9] Use iterative solvers --- inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml index ffab615c0f9..9b5049135d8 100644 --- a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml +++ b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml @@ -56,8 +56,8 @@ lineSearchMaxCuts="2" maxTimeStepCuts="1"/> @@ -77,8 +77,9 @@ newtonTol = "1.0e-6" newtonMaxIter ="8"/> + solverType="gmres" + preconditionerType="amg" + logLevel="0"/> From 20064980a3e01672132f6a2257ca526ec76eead1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vidar=20Stiernstr=C3=B6m?= Date: Wed, 29 Jan 2025 13:56:20 -0800 Subject: [PATCH 4/9] Adjust shear modulus in xml to translate from value used in antiplane shear to plane strain --- inputFiles/inducedSeismicity/SCEC_BP6_QD_S_base.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_base.xml b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_base.xml index 1257f1ae3ec..60a6832416d 100644 --- a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_base.xml +++ b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_base.xml @@ -83,7 +83,7 @@ name="rock" defaultDensity="2670" defaultBulkModulus="50.0e9" - defaultShearModulus="32.04e9"/> + defaultShearModulus="24.03e9"/> Date: Fri, 31 Jan 2025 13:21:18 -0800 Subject: [PATCH 5/9] wip: Update xml files and python script to include more mesurement points along the fault. Start debugging flow solver interaction with ExplicitQDRateAndState --- .../inducedSeismicity/SCEC_BP6_QD_S_base.xml | 79 +++++++++++++++---- .../SCEC_BP6_QD_S_explicit.xml | 19 +---- .../SCEC_BP6_QD_S_implicit.xml | 2 +- .../scripts/plotBP6_results.py | 49 +++++------- 4 files changed, 86 insertions(+), 63 deletions(-) diff --git a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_base.xml b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_base.xml index 60a6832416d..12e05d26fe6 100644 --- a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_base.xml +++ b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_base.xml @@ -18,9 +18,9 @@ name="mesh1" elementTypes="{ C3D8 }" xCoords="{ -10.0e3, 10.0e3 }" - yCoords="{ -0.5e3, 0.5e3 }" - zCoords="{ -20.0e3, -10.0e3, 10.0e3, 20.0e3 }" - nx="{ 100 }" + yCoords="{ -0.1e3, 0.1e3 }" + zCoords="{ -40.0e3, -10.0e3, 10.0e3, 40.0e3 }" + nx="{ 50 }" ny="{ 1 }" nz="{ 20, 201, 20 }" cellBlockNames="{ cb1 }"/> @@ -29,13 +29,58 @@ + xMin="{ -0.01, -1.01e3, -20.0e3 }" + xMax="{ 0.01, 1.01e3, 20.0e3 }"/> + + + + xMin="{ -0.01, -1.01e3, -0.1e3 }" + xMax="{ 0.01, 1.01e3, 0.1e3 }"/> + + + + + + + + + + + + + + + + + + + + + @@ -47,7 +92,7 @@ + defaultAperture="5e-3"/> @@ -83,7 +128,7 @@ name="rock" defaultDensity="2670" defaultBulkModulus="50.0e9" - defaultShearModulus="24.03e9"/> + defaultShearModulus="24.03e9"/> + setNames="{ receiver1, receiver2, receiver3, receiver4, receiver5, receiver6, receiver7, source}"/> + setNames="{ receiver1, receiver2, receiver3, receiver4, receiver5, receiver6, receiver7, source}"/> @@ -169,7 +214,7 @@ objectPath="ElementRegions/Fault" fieldName="stateVariable" scale="0.68070857" - setNames="{ all }" + setNames="{ faultPlane }" initialCondition="1"/> + setNames="{ faultPlane }"/> + setNames="{ faultPlane }"/> + setNames="{ faultPlane }"/> + setNames="{ faultPlane }"/> \ No newline at end of file diff --git a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml index 9b5049135d8..025017d8bff 100644 --- a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml +++ b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml @@ -83,9 +83,10 @@ - + + maxTime="1e6"> - - + maxTime="1e6"> Date: Fri, 31 Jan 2025 15:39:21 -0800 Subject: [PATCH 6/9] Try resetting variable_n after each step. --- .../inducedSeismicity/ExplicitQDRateAndState.cpp | 2 ++ .../inducedSeismicity/QuasiDynamicEarthQuake.cpp | 6 ++++++ .../inducedSeismicity/QuasiDynamicEarthQuake.hpp | 2 ++ 3 files changed, 10 insertions(+) diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp index d6fc3c74871..6494081aa52 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp @@ -91,6 +91,7 @@ real64 ExplicitQDRateAndState::solverStep( real64 const & time_n, real64 dtStage = m_butcherTable.c[1]*dtAdaptive; dtStress = updateStresses( time_n, dtStage, cycleNumber, domain ); updateSlipVelocity( time_n, dtStage, domain ); + resetStateToBeginningOfStep( domain ); // Remaining stages for( integer stageIndex = 1; stageIndex < m_butcherTable.numStages-1; stageIndex++ ) @@ -99,6 +100,7 @@ real64 ExplicitQDRateAndState::solverStep( real64 const & time_n, dtStage = m_butcherTable.c[stageIndex+1]*dtAdaptive; dtStress = updateStresses( time_n, dtStage, cycleNumber, domain ); updateSlipVelocity( time_n, dtStage, domain ); + resetStateToBeginningOfStep( domain ); } stepRateStateODEAndComputeError( dtAdaptive, domain ); diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.cpp index 7c45b81a7a1..3940d6e08ff 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.cpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.cpp @@ -107,6 +107,12 @@ real64 QuasiDynamicEarthQuake< RSSOLVER_TYPE >::updateStresses( real64 const & t return dtAccepted; } +template< typename RSSOLVER_TYPE > +void QuasiDynamicEarthQuake< RSSOLVER_TYPE >::resetStateToBeginningOfStep( DomainPartition & domain ) +{ + m_stressSolver->resetStateToBeginningOfStep( domain ); +} + template< typename RSSOLVER_TYPE > void QuasiDynamicEarthQuake< RSSOLVER_TYPE >::setTargetDispJump( DomainPartition & domain ) const { diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.hpp index 728ec97ca3d..aa870200a0d 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.hpp @@ -59,6 +59,8 @@ class QuasiDynamicEarthQuake : public RSSOLVER_TYPE const int cycleNumber, DomainPartition & domain ) const override final; + void resetStateToBeginningOfStep( DomainPartition & domain ) override final; + private: string m_stressSolverName; From c5eca39f039069d00d995da8f56f99a5de81a3fc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vidar=20Stiernstr=C3=B6m?= Date: Fri, 31 Jan 2025 16:20:04 -0800 Subject: [PATCH 7/9] Restructure reset calls to stress solver to work with the FSAL Runge-Kutta methods --- .../ExplicitQDRateAndState.cpp | 27 +++++++++++++++---- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp index 6494081aa52..a854b1c9009 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp @@ -86,23 +86,39 @@ real64 ExplicitQDRateAndState::solverStep( real64 const & time_n, { real64 dtStress; GEOS_UNUSED_VAR( dtStress ); + // // Initial Runge-Kutta stage + // + + // Evolve ODE:s for slip and state evolution stepRateStateODEInitialSubstage( dtAdaptive, domain ); - real64 dtStage = m_butcherTable.c[1]*dtAdaptive; + + // Compute stresses (linear mechanic + fluid solve) at next substage + real64 dtStage = m_butcherTable.c[1]*dtAdaptive; // Stage time step size dtStress = updateStresses( time_n, dtStage, cycleNumber, domain ); updateSlipVelocity( time_n, dtStage, domain ); - resetStateToBeginningOfStep( domain ); - // Remaining stages + // + // Remaining Runge-Kutta stages + // for( integer stageIndex = 1; stageIndex < m_butcherTable.numStages-1; stageIndex++ ) { + // Evolve ODE:s for slip and state evolution stepRateStateODESubstage( stageIndex, dtAdaptive, domain ); - dtStage = m_butcherTable.c[stageIndex+1]*dtAdaptive; + + // Compute stresses (linear mechanic + fluid solve) at next substage + // Need to reset stress solver to beginning of time step to not + // accumulate field in the stages. + resetStateToBeginningOfStep( domain ); + dtStage = m_butcherTable.c[stageIndex+1]*dtAdaptive; // Stage time step size dtStress = updateStresses( time_n, dtStage, cycleNumber, domain ); + + // Compute slip velocity using updated stresses and state updateSlipVelocity( time_n, dtStage, domain ); - resetStateToBeginningOfStep( domain ); } + // Evolve rate-and-state ODE:s to next time step and compute + // time step error stepRateStateODEAndComputeError( dtAdaptive, domain ); // Update timestep based on the time step error evalTimestep( domain ); @@ -111,6 +127,7 @@ real64 ExplicitQDRateAndState::solverStep( real64 const & time_n, // Compute stresses, and slip velocity and save results at time_n + dtAdapitve if( !m_butcherTable.FSAL ) { + resetStateToBeginningOfStep( domain ); // Reset stress fields dtStress = updateStresses( time_n, dtAdaptive, cycleNumber, domain ); updateSlipVelocity( time_n, dtAdaptive, domain ); } From 8dfaf678d9047002ee547ab62de783a787ff0821 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vidar=20Stiernstr=C3=B6m?= Date: Fri, 31 Jan 2025 16:21:05 -0800 Subject: [PATCH 8/9] Override resetStateToBeginningOfStep in SpringSlider --- .../physicsSolvers/inducedSeismicity/SpringSlider.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.hpp index 08c9a829c7e..c773b72aa33 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.hpp @@ -49,6 +49,12 @@ class SpringSlider : public RSSOLVER_TYPE real64 const & dt, const int cycleNumber, DomainPartition & domain ) const override final; + + void resetStateToBeginningOfStep( DomainPartition & domain ) override final + { + GEOS_UNUSED_VAR(domain); + return; + }; private: From a3e4fcaedf5e652c91366eb1cc7175749af955dd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vidar=20Stiernstr=C3=B6m?= Date: Fri, 31 Jan 2025 16:21:47 -0800 Subject: [PATCH 9/9] Add plot of slip rate to python script --- inputFiles/inducedSeismicity/scripts/plotBP6_results.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/inputFiles/inducedSeismicity/scripts/plotBP6_results.py b/inputFiles/inducedSeismicity/scripts/plotBP6_results.py index bd4ba0bea6e..79ea8fa99ab 100644 --- a/inputFiles/inducedSeismicity/scripts/plotBP6_results.py +++ b/inputFiles/inducedSeismicity/scripts/plotBP6_results.py @@ -57,7 +57,7 @@ def G( t, x, alpha ): # Plotting _, ax1 = plt.subplots() - # _, ax2 = plt.subplots() + _, ax2 = plt.subplots() positions_along_fault = [0., 500., 1500., 2500., 3500., 5000., 7500., -1500.] set_names = ["source", "receiver1", "receiver2", "receiver3", "receiver4", "receiver5", "receiver6", "receiver7"] @@ -73,10 +73,15 @@ def G( t, x, alpha ): ax1.tick_params(axis='y', labelcolor='tab:blue') ax1.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left') + _, slipRate = getDataFromHDF5( filePath, "slipRate", set_name) + ax2.plot(time_in_years, slipRate, label=f"Slip Rate z = {position} m") + ax2.set_yscale('log') + ax2.tick_params(axis='y', labelcolor='tab:red') + ax2.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left') # Set x-axis limits to 0 to 2 years ax1.set_xlim(0, np.max(time_in_years)) - # ax2.set_xlim(0, np.max(time_in_years)) + ax2.set_xlim(0, np.max(time_in_years)) # Add grid and title