Skip to content

Commit

Permalink
wip: Update xml files and python script to include more mesurement po…
Browse files Browse the repository at this point in the history 
…ints along the fault. Start debugging flow solver interaction with ExplicitQDRateAndState
  • Loading branch information
@VidarStiernstrom
VidarStiernstrom committed Jan 31, 2025
1 parent 2006498 commit 3fdd2a7
Show file tree
Hide file tree
Showing 4 changed files with 86 additions and 63 deletions.
79 changes: 62 additions & 17 deletions inputFiles/inducedSeismicity/SCEC_BP6_QD_S_base.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 }"/>
Expand All @@ -29,13 +29,58 @@
<Geometry>
<Box
name="faultPlane"
xMin="{ -0.01, -1.01e3, -10.0e3 }"
xMax="{ 0.01, 1.01e3, 10.0e3 }"/>
xMin="{ -0.01, -1.01e3, -20.0e3 }"
xMax="{ 0.01, 1.01e3, 20.0e3 }"/>

<!-- Boxes for single-cell sources/receivers on fault -->

<!-- Source at z = 0 m -->
<Box
name="source"
xMin="{ -0.01, -1.01e3, -50.0 }"
xMax="{ 0.01, 1.01e3, 50.0 }"/>
xMin="{ -0.01, -1.01e3, -0.1e3 }"
xMax="{ 0.01, 1.01e3, 0.1e3 }"/>

<!-- Receiver at z = 500 m-->
<Box
name="receiver1"
xMin="{ -0.01, -1.01e3, 0.4e3 }"
xMax="{ 0.01, 1.01e3, 0.6e3 }"/>

<!-- Receiver at z = 1500 m-->
<Box
name="receiver2"
xMin="{ -0.01, -1.01e3, 1.4e3 }"
xMax="{ 0.01, 1.01e3, 1.6e3 }"/>

<!-- Receiver at z = 2500 m-->
<Box
name="receiver3"
xMin="{ -0.01, -1.01e3, 2.4e3 }"
xMax="{ 0.01, 1.01e3, 2.6e3 }"/>

<!-- Receiver at z = 3500 m-->
<Box
name="receiver4"
xMin="{ -0.01, -1.01e3, 3.4e3 }"
xMax="{ 0.01, 1.01e3, 3.6e3 }"/>

<!-- Receiver at z = 5000 m-->
<Box
name="receiver5"
xMin="{ -0.01, -1.01e3, 4.9e3 }"
xMax="{ 0.01, 1.01e3, 5.1e3 }"/>

<!-- Receiver at z = 7500 m-->
<Box
name="receiver6"
xMin="{ -0.01, -1.01e3, 7.4e3 }"
xMax="{ 0.01, 1.01e3, 7.6e3 }"/>

<!-- Receiver at z = -1500 m-->
<Box
name="receiver7"
xMin="{ -0.01, -1.01e3, -1.6e3 }"
xMax="{ 0.01, 1.01e3, -1.4e3 }"/>
</Geometry>

<ElementRegions>
Expand All @@ -47,7 +92,7 @@
<SurfaceElementRegion
name="Fault"
materialList="{frictionLaw, faultMaterial, water}"
defaultAperture="1e-3"/>
defaultAperture="5e-3"/>
</ElementRegions>

<Constitutive>
Expand Down Expand Up @@ -83,7 +128,7 @@
name="rock"
defaultDensity="2670"
defaultBulkModulus="50.0e9"
defaultShearModulus="24.03e9"/>
defaultShearModulus="24.03e9"/> <!-- Note: Adjusted to translate antiplane shear to plane strain -->

<RateAndStateFriction
name="frictionLaw"
Expand All @@ -110,13 +155,13 @@
name="slipRateCollection"
objectPath="ElementRegions/Fault/FractureSubRegion"
fieldName="slipRate"
setNames="{source}"/>
setNames="{ receiver1, receiver2, receiver3, receiver4, receiver5, receiver6, receiver7, source}"/>

<PackCollection
name="pressureCollection"
objectPath="ElementRegions/Fault/FractureSubRegion"
fieldName="pressure"
setNames="{source}"/>
setNames="{ receiver1, receiver2, receiver3, receiver4, receiver5, receiver6, receiver7, source}"/>
</Tasks>

<ElementRegions>
Expand Down Expand Up @@ -169,7 +214,7 @@
objectPath="ElementRegions/Fault"
fieldName="stateVariable"
scale="0.68070857"
setNames="{ all }"
setNames="{ faultPlane }"
initialCondition="1"/>

<FieldSpecification
Expand All @@ -179,7 +224,7 @@
objectPath="ElementRegions/Fault/FractureSubRegion"
component="0"
scale="1.0e-12"
setNames="{all}"/>
setNames="{ faultPlane }"/>

<FieldSpecification
name="slipVelocity2"
Expand All @@ -188,7 +233,7 @@
objectPath="ElementRegions/Fault/FractureSubRegion"
component="1"
scale="0."
setNames="{all}"/>
setNames="{ faultPlane }"/>

<FieldSpecification
name="backgroundShearStress"
Expand All @@ -197,21 +242,21 @@
objectPath="ElementRegions/Fault/FractureSubRegion"
component="0"
scale="29.2e6"
setNames="{all}"/>
setNames="{ faultPlane }"/>

<FieldSpecification
name="backgroundNormalStress"
fieldName="backgroundNormalStress"
initialCondition="1"
objectPath="ElementRegions/Fault/FractureSubRegion"
scale="50.0e6"
setNames="{all}"/>
setNames="{ faultPlane }"/>

<SourceFlux
name="sourceTerm"
objectPath="ElementRegions/Fault"
scale="-1.25e-3"
setNames="{ source }"
setNames="{source}"
endTime="8.64e6"/>
</FieldSpecifications>
</Problem>
19 changes: 4 additions & 15 deletions inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,9 +83,10 @@
</SinglePhaseFVM>
</Solvers>


<!-- <Events
maxTime="6.312e7"> -->
<Events
maxTime="6.312e7">
maxTime="1e6">

<!-- fault generation -->
<SoloEvent
Expand All @@ -101,21 +102,9 @@

<!-- solver application -->
<PeriodicEvent
name="solverApplications1"
name="solverApplications"
maxEventDt="1e4"
endTime="8.64e6"
target="/Solvers/QDSolver"/>
<PeriodicEvent
name="solverApplications2"
maxEventDt="10"
beginTime="8.64e6"
endTime="8.6499e6"
target="/Solvers/QDSolver"/>
<PeriodicEvent
name="solverApplications3"
maxEventDt="1e3"
beginTime="8.6499e6"
target="/Solvers/QDSolver"/>

<!-- Time history collection -->
<PeriodicEvent
Expand Down
2 changes: 1 addition & 1 deletion inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,7 @@


<Events
maxTime="6.312e7">
maxTime="1e6">

<!-- fault generation -->
<SoloEvent
Expand Down
49 changes: 19 additions & 30 deletions inputFiles/inducedSeismicity/scripts/plotBP6_results.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,10 +13,10 @@ def remove_padding(data):
last_nonzero = nonzero_indices[-1]
return data[:last_nonzero + 1]

def getDataFromHDF5( hdf5FilePath, var_name ):
def getDataFromHDF5( hdf5FilePath, var_name, set_name):
# Read HDF5
data = hdf5_wrapper(f'{hdf5FilePath}').get_copy()
var = data[f'{var_name} source']
var = data[f'{var_name} {set_name}']
var = np.asarray(var)
time = data[f'{var_name} Time']

Expand Down Expand Up @@ -55,45 +55,34 @@ def G( t, x, alpha ):
print(f"Error: {filePath} not found.")
exit(1)

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

# Plotting
fig, ax1 = plt.subplots()
_, ax1 = plt.subplots()
# _, ax2 = plt.subplots()

print(len(time_in_years))
print(len(pressure))
print(len(pressure_analytical))

positions_along_fault = [0., 500., 1500., 2500., 3500., 5000., 7500., -1500.]
set_names = ["source", "receiver1", "receiver2", "receiver3", "receiver4", "receiver5", "receiver6", "receiver7"]
# 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)
ax2 = ax1.twinx()
ax2.set_ylabel('Slip Rate (m/s)', color='tab:red')
ax2.plot(time_in_years, slipRate, label="Slip Rate", color='tab:red')
ax2.set_yscale('log')
ax2.tick_params(axis='y', labelcolor='tab:red')
for position, set_name in zip(positions_along_fault, set_names):
time, pressure = getDataFromHDF5( filePath, "pressure" , set_name)
time_in_years = time / (365 * 24 * 3600) # Convert time to years, assuming time is in seconds
pressure_analytical = analytical_pressure( time, position )
ax1.plot(time_in_years, pressure, label=f"Pressure z = {position} m")
ax1.plot(time_in_years, pressure_analytical, label=f"Pressure Analytical z = {position} m", linestyle='--')
ax1.set_xlabel('Time (years)')
ax1.set_ylabel('Pressure (Pa)', color='tab:blue')
ax1.tick_params(axis='y', labelcolor='tab:blue')
ax1.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))


# Add grid and title
plt.title("Pressure and Slip Rate vs Time")
plt.grid()
plt.tight_layout()

# Show plot
plt.show()

0 comments on commit 3fdd2a7

Please sign in to comment.