Merge pull request #9 from pnnl/separateSamplingFromRisk

Separate the GenerateSamples out of the EvaluateFail

SOSAT/risk_analysis/critically_oriented_fault.py

@@ -52,7 +52,30 @@ def __init__(self,
         else:
             self.mu_dist = mu_dist
 
-    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):
         """
         A method to evaluate the failure probability at pressures
         between the native pore pressure at self.dPmax.
@@ -66,6 +89,10 @@ def EvaluatePfail(self, Npressures=20, Nsamples=1e6):
         Nsamples : int
             The number of stress samples to use for the analysis
 
+        shmin, shmax, sv: arrays containing samples of the three
+        principal stress; default to be None; They can be calculated
+        using self.SampleStressPoints(Nsamples=1e6)
+
         Returns
         -------
         P, pfail : `numpy.ndarray`
@@ -74,10 +101,15 @@ def EvaluatePfail(self, Npressures=20, Nsamples=1e6):
         """
         Nsamples = int(Nsamples)
         Npressures = int(Npressures)
-        # generate samples of stress state
-        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)
+
         gamma = self.gamma_dist.rvs(Nsamples)
         mu = self.mu_dist.rvs(Nsamples)
 

SOSAT/risk_analysis/hydraulic_fracture.py

@@ -43,7 +43,30 @@ def __init__(self,
         self.T_dist = T_dist
         self.T_unit = T_unit
 
-    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):
         """
         A method to evaluate the failure probability at pressures
         between the native pore pressure at self.dPmax.
@@ -57,6 +80,10 @@ def EvaluatePfail(self, Npressures=20, Nsamples=1e6):
         Nsamples : int
             The number of stress samples to use for the analysis
 
+        shmin, shmax, sv: arrays containing samples of the three
+        principal stress; default to be None; They can be calculated
+        using self.SampleStressPoints(Nsamples=1e6)
+
         Returns
         -------
         P, pfail : `numpy.ndarray`
@@ -66,9 +93,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)

tests/risk_analysis/test_critical_fault_activation_risk.py

@@ -31,7 +31,9 @@ def test_cfa():
     dPmax = 4.0
     gamma_dist = uniform(0.4, (0.6 - 0.4))
     cfa = CriticalFaultActivation(ss, dPmax, gamma_dist)
-    pressures, Pfail = cfa.EvaluatePfail()
+    shmin, shmax, sv = cfa.SampleStressPoints()
+    pressures, Pfail = cfa.EvaluatePfail(
+        shmin=shmin, shmax=shmax, sv=sv)
     print("pressures= ", pressures)
     print("Pfail= ", Pfail)
     pressures[0] == pytest.approx(12.227)

tests/risk_analysis/test_hydraulic_fracturing_risk.py

@@ -30,8 +30,9 @@ def test_hf():
     T_dist = uniform(0.0, 5.0)
 
     hf = HydraulicFracturing(ss, dPmax, gamma_dist)
-
-    pressures, Pfail = hf.EvaluatePfail()
+    shmin, shmax, sv = hf.SampleStressPoints(Nsamples=1e5)
+    pressures, Pfail = hf.EvaluatePfail(
+        shmin=shmin, shmax=shmax, sv=sv)
 
     print("pressures= ", pressures)
     print("Pfail= ", Pfail)