Merge pull request #13 from pnnl/changeplot2dP

change function to read in presample stress and fault points

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()
 

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")