Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

change function to read in presample stress and fault points #13

Merged
merged 1 commit into from
Jun 6, 2024
Merged

change function to read in presample stress and fault points #13

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
70 changes: 60 additions & 10 deletions SOSAT/risk_analysis/user_prescribed_fault.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -253,7 +253,32 @@ def CreateOrientations(self, Nsamples=1e6):

return n_shmax, n_fault, n_fault_p, ss_azi, strike_fault, dip_fault

def EvaluatePfail(self, Npressures=20, Nsamples=1e6):
def SampleStressPoints(self, Nsamples=1e6):
"""
A method to sample three principal stress points
from the joint stress posterior distribution

Parameters
----------
Nsamples : int
The number of stress samples to use for the analysis

Returns
-------
shmin, shmax, sv: arrays containing samples of the three
principal stress

"""
Nsamples = int(Nsamples)
# generate samples of stress state
stress_sampler = RejectionSampler(self.ss)
shmin, shmax, sv = stress_sampler.GenerateSamples(Nsamples)

return shmin, shmax, sv

def EvaluatePfail(self, Npressures=20, Nsamples=1e6,
shmin=None, shmax=None, sv=None,
n_fault_p=None):
"""
A method to evaluate the failure probability at pressures
between the native pore pressure at self.dPmax.
Expand All @@ -267,6 +292,16 @@ def EvaluatePfail(self, Npressures=20, Nsamples=1e6):
Nsamples : int
The number of stress samples to use for the analysis

shmin, shmax, sv: `numpy.ndarray`
arrays containing samples of the three
principal stress; default to be None; They can be calculated
using self.SampleStressPoints(Nsamples=1e6)

n_fault_p: `numpy.ndarray`
array containing the fault normal vectors in the
principal coordinate system; default to be None;
It can be calculated using self.CreateOrientations(Nsamples=1e6)

Returns
-------
P, pfail : `numpy.ndarray`
Expand All @@ -276,9 +311,13 @@ def EvaluatePfail(self, Npressures=20, Nsamples=1e6):
Nsamples = int(Nsamples)
Npressures = int(Npressures)

# Generate samples of stress state magnitudes
stress_sampler = RejectionSampler(self.ss)
shmin, shmax, sv = stress_sampler.GenerateSamples(Nsamples)
# generate samples of stress state if there is no stress inputs
if shmin is None or shmax is None or sv is None:
stress_sampler = RejectionSampler(self.ss)
shmin, shmax, sv = stress_sampler.GenerateSamples(Nsamples)
else:
# make sure the Nsamples equals to the given length of shmin
Nsamples = len(sv)

# Generate samples of stress path coefficients
gamma = self.gamma_dist.rvs(Nsamples)
Expand All @@ -287,8 +326,16 @@ def EvaluatePfail(self, Npressures=20, Nsamples=1e6):
mu = self.mu_dist.rvs(Nsamples)

# Create orientations for stress state and fault
n_shmax, n_fault, n_fault_p, ss_azi, strike_fault, dip_fault = \
self.CreateOrientations(Nsamples)
if n_fault_p is None:
_, _, n_fault_p, _, _, _ = \
self.CreateOrientations(Nsamples)
elif len(n_fault_p) != Nsamples:
error_message = 'The number of fault samples inputted should \
equal to the number of stress samples. \
Use the regenerated fault samples instead.'
print(error_message)
_, _, n_fault_p, _, _, _ = \
self.CreateOrientations(Nsamples)

Po = self.ss.pore_pressure
dP_array = np.linspace(0.0, self.dPmax, Npressures)
Expand Down Expand Up @@ -327,14 +374,17 @@ def EvaluatePfail(self, Npressures=20, Nsamples=1e6):
def PlotFailureProbability(self,
Npressures=20,
Nsamples=1e6,
figwidth=5.0):
figwidth=5.0,
n_fault_p=None):

P, Pfail = self.EvaluatePfail(Npressures, Nsamples)
P, Pfail = self.EvaluatePfail(Npressures, Nsamples,
n_fault_p=n_fault_p)
dP = P - self.ss.pore_pressure

fig = plt.figure(figsize=(figwidth, figwidth * 0.7))
ax = fig.add_subplot(111)
ax.plot(P, Pfail, "k")
ax.set_xlabel("Pore Pressure")
ax.plot(dP, Pfail, "k") # change to plot dP versus Pfail
ax.set_xlabel("Pore Pressure Change (" + self.ss.stress_unit + ")")
ax.set_ylabel("Probability of Fault Activation")
plt.tight_layout()

Expand Down
26 changes: 25 additions & 1 deletion tests/risk_analysis/test_user_fault_activation_risk.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,8 @@ def test_upfa():
fault_K=fault_K)
n_shmax, n_fault, n_fault_p, ss_azi, strike_fault, dip_fault =\
upfa.CreateOrientations()

# case 1: didn't pre-sample stress and fault points
pressures, Pfail = upfa.EvaluatePfail()

print("mean ss_azi = ", np.mean(ss_azi))
Expand All @@ -59,6 +61,28 @@ def test_upfa():
pressures[-1] == pytest.approx(ss.pore_pressure + dPmax)
Pfail[-1] == pytest.approx(0.5153, rel=0.5e-1)

fig = upfa.PlotFailureProbability()
# case 2: pre-sample stress points
shmin, shmax, sv = upfa.SampleStressPoints()
pressures, Pfail = upfa.EvaluatePfail(
shmin=shmin, shmax=shmax, sv=sv)

pressures[0] == pytest.approx(12.227, rel=1e-2)
Pfail[0] == pytest.approx(0.0533, rel=0.5e-1)
pressures[-1] == pytest.approx(ss.pore_pressure + dPmax)
Pfail[-1] == pytest.approx(0.5153, rel=0.5e-1)

# case 3: pre-sample stress points and fault points
_, _, n_fault_p, _, _, _ = upfa.CreateOrientations()
pressures, Pfail = upfa.EvaluatePfail(
shmin=shmin, shmax=shmax, sv=sv, n_fault_p=n_fault_p)

pressures[0] == pytest.approx(12.227, rel=1e-2)
Pfail[0] == pytest.approx(0.0533, rel=0.5e-1)
pressures[-1] == pytest.approx(ss.pore_pressure + dPmax)
Pfail[-1] == pytest.approx(0.5153, rel=0.5e-1)

fig = upfa.PlotFailureProbability()
fig.savefig("User_Fault_Activation_Probability.png")

fig = upfa.PlotFailureProbability(n_fault_p=n_fault_p)
fig.savefig("User_Fault_Activation_Probability_inputfault.png")
Loading