diff --git a/SOSAT/risk_analysis/user_prescribed_fault.py b/SOSAT/risk_analysis/user_prescribed_fault.py index 813f40b..dc5af0e 100644 --- a/SOSAT/risk_analysis/user_prescribed_fault.py +++ b/SOSAT/risk_analysis/user_prescribed_fault.py @@ -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. @@ -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` @@ -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) @@ -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) @@ -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() diff --git a/tests/risk_analysis/test_user_fault_activation_risk.py b/tests/risk_analysis/test_user_fault_activation_risk.py index e340ec0..1ad9a5c 100644 --- a/tests/risk_analysis/test_user_fault_activation_risk.py +++ b/tests/risk_analysis/test_user_fault_activation_risk.py @@ -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)) @@ -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")