Merge pull request #24 from TchilDill/extras-doc-unit-testing-and-more

0.5.0

CHANGELOG.md

@@ -7,8 +7,14 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/)
 and this project adheres to [Semantic Versioning](http://semver.org/),
 and [PEP 440](https://www.python.org/dev/peps/pep-0440/).
 
+## [0.5.1] - 2023-XX-XX (Note released yet)  
+- *blank*
 
-## [0.4.0] - 2023-07-30 (not released yet) 
+## [0.5.0] - 2023-08-02  
+- added soil spring `openpile.utils.py_curves.reese_weakrock`
+- added soil model `openpile.soilmodels.Reese_weakrock`
+
+## [0.4.0] - 2023-07-30 
 - Updates in documentation
 - new methods available for `openpile.analyze.Result` class:
   - new method `winkler()` replacing `simple_winkler_analysis()`, the use of the latter triggers a deprecation warning when used.

docs/source/API.rst

@@ -17,7 +17,7 @@ API
 
 .. automodule:: openpile.analyze
     :members: 
-    :exclude-members: simple_winkler_analysis, simple_beam_analysis, PydanticConfig, structural_forces_to_df, springs_mob_to_df, reaction_forces_to_df, disp_to_df, AnalyzeResult
+    :exclude-members: winkler, beam, PydanticConfig, structural_forces_to_df, springs_mob_to_df, reaction_forces_to_df, disp_to_df, AnalyzeResult
 
 
 `utils` module
@@ -43,6 +43,10 @@ API
     :exclude-members: random
 
 
+.. automodule:: openpile.utils.multipliers
+    :members:
+    :exclude-members: 
+
 The `Result` class
 ==================
 

src/openpile/globals.py

@@ -1,2 +1,2 @@
 # version of the package
-VERSION = "0.4.0"
+VERSION = "0.5.0"

src/openpile/soilmodels.py

@@ -793,3 +793,91 @@ def mt_spring_fct(
         m = 1 / 4 * np.pi * D**2 * tz_pos.reshape((1, -1)) * np.ones((output_length, 1))
 
         return t, m
+
+
+@dataclass(config=PydanticConfigFrozen)
+class Reese_weakrock(LateralModel):
+    """A class to establish the Reese weakrock model.
+
+    Parameters
+    ----------
+    Ei: float or list[top_value, bottom_value]
+        Initial modulus of rock [unit: kPa]
+    qu: float or list[top_value, bottom_value]
+        compressive strength of rock [unit: kPa]
+    RQD: float or list[top_value, bottom_value]
+        Rock Quality Designation [unit: %]
+    k: float
+        dimensional constant randing from 0.0005 to 0.00005, by default 0.0005
+    ztop: float
+        absolute depth of top layer elevation with respect to rock surface [m]
+    p_multiplier: float or function taking the depth as argument and returns the multiplier
+        multiplier for p-values
+    y_multiplier: float or function taking the depth as argument and returns the multiplier
+        multiplier for y-values
+
+    """
+
+    #: initial modulus of rock [kPa], if a variation in values, two values can be given.
+    Ei: Union[PositiveFloat, conlist(PositiveFloat, min_items=1, max_items=2)]
+    #: scompressive strength of rock [kPa], if a variation in values, two values can be given.
+    qu: Optional[Union[PositiveFloat, conlist(PositiveFloat, min_items=1, max_items=2)]]
+    #: Rock Quality Designation
+    RQD: confloat(ge=0.0, le=100.0)
+    #: dimnesional constant
+    k: confloat(ge=0.00005, le=0.0005)
+    #: absolute depth of top layer elevation with respect to rock surface [m]
+    ztop: confloat(ge=0.0)
+    #: p-multiplier
+    p_multiplier: Union[Callable[[float], float], confloat(ge=0.0)] = 1.0
+    #: y-multiplier
+    y_multiplier: Union[Callable[[float], float], confloat(gt=0.0)] = 1.0
+
+    # define class variables needed for all soil models
+    m_multiplier = 1.0
+    t_multiplier = 1.0
+
+    spring_signature = np.array([True, False, False, False], dtype=bool)
+
+    def __str__(self):
+        return f"\tReese weakrock\n\tEi = {var_to_str(self.Ei)}kPa\n\tqu = {var_to_str(self.qu)}kPa\n\tRQD = {var_to_str(self.RQD)}%"
+
+    def py_spring_fct(
+        self,
+        sig: float,
+        X: float,
+        layer_height: float,
+        depth_from_top_of_layer: float,
+        D: float,
+        L: float = None,
+        below_water_table: bool = True,
+        ymax: float = 0.0,
+        output_length: int = 15,
+    ):
+        # validation
+        if depth_from_top_of_layer > layer_height:
+            raise ValueError("Spring elevation outside layer")
+
+        # define Ei
+        Ei_t, Ei_b = from_list2x_parse_top_bottom(self.Ei)
+        Ei = Ei_t + (Ei_b - Ei_t) * depth_from_top_of_layer / layer_height
+
+        # define qu
+        qu_t, qu_b = from_list2x_parse_top_bottom(self.qu)
+        qu = qu_t + (qu_b - qu_t) * depth_from_top_of_layer / layer_height
+
+        y, p = py_curves.reese_weakrock(
+            Ei=Ei,
+            xr=(X + self.ztop),
+            RQD=self.RQD,
+            qu=qu,
+            D=D,
+            k=self.k,
+            output_length=output_length,
+        )
+
+        # parse multipliers and apply results
+        y_mult = self.y_multiplier if isinstance(self.y_multiplier, float) else self.y_multiplier(X)
+        p_mult = self.p_multiplier if isinstance(self.p_multiplier, float) else self.p_multiplier(X)
+
+        return y * y_mult, p * p_mult

src/openpile/utils/Hb_curves.py

@@ -47,10 +47,10 @@ def dunkirk_sand(
 
     Returns
     -------
-    1darray
-        Hb vector [unit: kN]
     1darray
         y vector [unit: m]
+    1darray
+        Hb vector [unit: kN]
     """
     # correct relative density for decimal value
     Dr = Dr / 100
@@ -109,17 +109,10 @@ def cowden_clay(
 
     Returns
     -------
-    1darray
-        Hb vector [unit: kN]
     1darray
         y vector [unit: m]
-
-    References
-    ----------
-    .. [1] Byrne, B. W., Houlsby, G. T., Burd, H. J., Gavin, K. G., Igoe, D. J. P., Jardine,
-           R. J., Martin, C. M., McAdam, R. A., Potts, D. M., Taborda, D. M. G. & Zdravkovic ́,
-           L. (2020). PISA design model for monopiles for offshore wind turbines: application
-           to a stiff glacial clay till. Géotechnique, https://doi.org/10.1680/ jgeot.18.P.255.
+    1darray
+        Hb vector [unit: kN]
 
     """
 

src/openpile/utils/Mb_curves.py

@@ -44,10 +44,10 @@ def cowden_clay(
 
     Returns
     -------
-    1darray
-        Mb vector [unit: kN]
     1darray
         t vector of length [unit: rad]
+    1darray
+        Mb vector [unit: kN]
 
     """
 
@@ -107,10 +107,10 @@ def dunkirk_sand(
 
     Returns
     -------
-    1darray
-        Mb vector [unit: kN]
     1darray
         t vector of length [unit: rad]
+    1darray
+        Mb vector [unit: kN]
     """
     # correct relative density for decimal value
     Dr = Dr / 100

src/openpile/utils/mt_curves.py

@@ -41,10 +41,10 @@ def cowden_clay(
 
     Returns
     -------
-    1darray
-        m vector [unit: kN]
     1darray
         t vector of length [unit: rad]
+    1darray
+        m vector [unit: kN]
 
     """
 
@@ -106,10 +106,10 @@ def dunkirk_sand(
 
     Returns
     -------
-    1darray
-        m vector [unit: kN]
     1darray
         t vector of length [unit: rad]
+    1darray
+        m vector [unit: kN]
     """
     # correct relative density for decimal value
     Dr = Dr / 100

src/openpile/utils/py_curves.py

@@ -41,10 +41,10 @@ def cowden_clay(
 
     Returns
     -------
-    1darray
-        p vector [unit: kN/m]
     1darray
         y vector [unit: m]
+    1darray
+        p vector [unit: kN/m]
 
     """
 
@@ -104,10 +104,10 @@ def dunkirk_sand(
 
     Returns
     -------
-    1darray
-        p vector [unit: kN/m]
     1darray
         y vector [unit: m]
+    1darray
+        p vector [unit: kN/m]
     """
     # correct relative density for decimal value
     Dr = Dr / 100
@@ -169,10 +169,10 @@ def api_sand(
 
     Returns
     -------
-    1darray
-        p vector [unit: kN/m]
     1darray
         y vector [unit: m]
+    1darray
+        p vector [unit: kN/m]
     """
     # A value - only thing that changes between cyclic or static
     if kind == "static":
@@ -273,11 +273,10 @@ def api_clay(
 
     Returns
     -------
-    1darray
-        p vector [unit: kN/m]
     1darray
         y vector [unit: m]
-    ---------
+    1darray
+        p vector [unit: kN/m]
     """
     # important variables
     y50 = 2.5 * eps50 * D
@@ -356,3 +355,74 @@ def api_clay(
             p[i] = 0.23 * Pmax
 
     return y, p
+
+
+@njit(parallel=True, cache=True)
+def reese_weakrock(
+    Ei: float,
+    qu: float,
+    RQD: float,
+    xr: float,
+    D: float,
+    k: float = 0.0005,
+    output_length=20,
+):
+    """creates the Reese weakrock p-y curve based on the work of Reese (1997) #TODO(ref).
+
+
+    Parameters
+    ----------
+    Ei : float
+        initial modulus of rock [kPa]
+    qu : float
+        compressive strength of rock [kPa]
+    RQD : float
+        Rock Quality Designation [%]
+    xr : float
+        depth from rock surface [m]
+    D : float
+        pile width [m]
+    k : float, optional
+        dimensional constant, randing from 0.0005 to 0.00005, by default 0.0005
+    output_length : int, optional
+        length of output arrays, by default 20
+
+    Returns
+    -------
+    1darray
+        y vector [unit: m]
+    1darray
+        p vector [unit: kN/m]
+
+    """
+
+    # Rqd forced to be within 0 and 100
+    rqd = max(min(100, RQD), 0)
+
+    # determine alpha
+    alpha = 1.0 - 2 / 3 * rqd / 100
+
+    # determine ultimate resistance of rock
+    Pmax = min(5.2 * alpha * qu * D, alpha * qu * D * (1 + 1.4 * xr / D))
+
+    # initial portion of p-y curve
+    Epy_i = Ei * min(500, 100 + 400 * xr / (3 * D))
+
+    # yA & yrm
+    yrm = k * D
+    yA = (Pmax / (2 * (yrm) ** 0.25 * Epy_i)) ** 1.333
+
+    # define y
+    ymax = max(1.05 * yA, (2 * yrm ** (0.25)) ** 4)
+    y = np.linspace(yA, ymax, output_length - 2)
+    y = np.concatenate((np.array([0.0]), y, np.array([1.2 * ymax])))
+
+    # define p
+    p = np.zeros(y.size)
+    for i in range(len(p)):
+        if y[i] <= yA:
+            p[i] = min(Pmax, Epy_i * y[i])
+        else:
+            p[i] = min(Pmax, Pmax / 2 * (y[i] / yrm) ** 0.25)
+
+    return y, p

src/openpile/utils/qz_curves.py

@@ -38,10 +38,10 @@ def api_sand(
 
     Returns
     -------
-    numpy 1darray
-        Q vector [unit: kPa]
     numpy 1darray
         z vector [unit: m]
+    numpy 1darray
+        Q vector [unit: kPa]
     """
     # cannot have less than 8
     if output_length < 8:
@@ -100,10 +100,10 @@ def api_clay(
 
     Returns
     -------
-    numpy 1darray
-        Q vector [unit: kPa]
     numpy 1darray
         z vector [unit: m]
+    numpy 1darray
+        Q vector [unit: kPa]
     """
     # cannot have less than 8
     if output_length < 8:

test/test_pycurves.py

@@ -114,3 +114,26 @@ def test_api_clay(make_pmax_api_clay,xsigma, xX, xSu, xe50, xD, xJ, xkind):
     assert m.isclose(p[-1],pres, rel_tol=0.01, abs_tol=0.1)
     assert m.isclose(np.max(p), pu, rel_tol=0.01, abs_tol=0.1)
 
+@pytest.mark.parametrize("xEi",[1,75,200])
+@pytest.mark.parametrize("xqu",[1,75,200])
+@pytest.mark.parametrize("xRQD",[0,25,50,75,100])
+@pytest.mark.parametrize("xxr",[0,5,20])
+@pytest.mark.parametrize("xD",[1,3,6])
+def test_reese_weakrock(xEi, xqu, xRQD, xxr, xD):
+    #create curve
+    y, p = py.reese_weakrock(Ei=xEi, qu=xqu, RQD=xRQD, xr=xxr, D=xD)
+    #helper fct
+    is_sorted = lambda a: np.all(a[:-1] <= a[1:])
+    #check if sorted
+    assert is_sorted(y)
+    #check if origin is (0,0)
+    assert p[0] == 0.0
+    assert y[0] == 0.0
+
+    alpha = 1 - 2/3 * xRQD/100
+    pu = min(alpha * xqu * xD * (1+1.4*xxr/xD), 5.2*alpha*xqu*xD)
+    assert m.isclose(np.max(p), pu, rel_tol=0.01, abs_tol=0.1)
+
+
+
+