new states for steady state

CITATION.cff

@@ -8,7 +8,7 @@ authors:
 - family-names: "Weissgraeber"
   given-names: "Philipp"
   orcid: "https://orcid.org/0000-0001-8320-8672"
-version: 3.1.1
+version: 3.1.2
 date-released: 2021-12-30
 identifiers:
   - description: Collection of archived snapshots of all versions of WEAC

demo/demo.ipynb

@@ -13,7 +13,7 @@
    "id": "695bafcb",
    "metadata": {},
    "source": [
-    "Note that instructions in this notebook refer to **release v3.1.1.** Please make sure you are running the latest version of weac using\n",
+    "Note that instructions in this notebook refer to **release v3.1.2.** Please make sure you are running the latest version of weac using\n",
     "\n",
     "```bash\n",
     "pip install -U weac\n",

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "weac"
-version = "3.1.1"
+version = "3.1.2"
 authors = [{ name = "2phi GbR", email = "mail@2phi.de" }]
 description = "Weak layer anticrack nucleation model"
 readme = "README.md"
@@ -123,7 +123,7 @@ ignore = [
 ]
 
 [tool.bumpversion]
-current_version = "3.1.1"
+current_version = "3.1.2"
 
 [[tool.bumpversion.files]]
 filename = "pyproject.toml"

src/weac/__init__.py

@@ -2,4 +2,4 @@
 WEAC - Weak Layer Anticrack Nucleation Model
 """
 
-__version__ = "3.1.1"
+__version__ = "3.1.2"

src/weac/analysis/criteria_evaluator.py

@@ -24,6 +24,7 @@
     WeakLayer,
 )
 from weac.constants import RHO_ICE
+from weac.core.slab_touchdown import TouchdownMode
 from weac.core.system_model import SystemModel
 
 logger = logging.getLogger(__name__)
@@ -683,6 +684,7 @@ def evaluate_coupled_criterion(
     def evaluate_SteadyState(
         self,
         system: SystemModel,
+        mode: TouchdownMode = "C_in_contact",
         vertical: bool = False,
         print_call_stats: bool = False,
     ) -> SteadyStateResult:
@@ -710,18 +712,33 @@ def evaluate_SteadyState(
                 UserWarning,
             )
         system_copy = copy.deepcopy(system)
+        # Evaluate touchdown distance for flat slab
         system_copy.toggle_touchdown(True)
-        system_copy.update_scenario(scenario_config=ScenarioConfig(phi=0.0))
-        l_BC = system_copy.slab_touchdown.l_BC
+        segments = [
+            Segment(length=5e3, has_foundation=True, m=0.0),
+            Segment(length=5e3, has_foundation=False, m=0.0),
+        ]
+        system_copy.update_scenario(
+            segments=segments, scenario_config=ScenarioConfig(phi=0.0)
+        )
+
+        cut_distance = 0
+        match mode:
+            case "C_in_contact":
+                cut_distance = 2 * system_copy.slab_touchdown.l_BC
+            case "B_point_contact":
+                cut_distance = system_copy.slab_touchdown.l_BC - 1e-3
+            case "A_free_hanging":
+                cut_distance = system_copy.slab_touchdown.l_AB - 1e-3
 
         segments = [
             Segment(length=5e3, has_foundation=True, m=0.0),
-            Segment(length=2 * l_BC, has_foundation=False, m=0.0),
+            Segment(length=cut_distance, has_foundation=False, m=0.0),
         ]
         scenario_config = ScenarioConfig(
             system_type="vpst-" if vertical else "pst-",
             phi=0.0,  # Slab Touchdown works only for flat slab
-            cut_length=2 * l_BC,
+            cut_length=cut_distance,
         )
         system_copy.update_scenario(segments=segments, scenario_config=scenario_config)
         touchdown_distance = system_copy.slab_touchdown.touchdown_distance

src/weac/components/layer.py

@@ -94,8 +94,43 @@ def _sigrist_tensile_strength(rho, unit: Literal["kPa", "MPa"] = "kPa"):
     return convert[unit] * 240 * (rho / RHO_ICE) ** 2.44
 
 
-# TODO: Compressive Strength from Schöttner
-# (11 +/- 7) * (rho/rho_0) ^ (5.4 +/- 0.5)
+def _adam_tensile_strength(rho, unit: Literal["kPa", "MPa"] = "kPa"):
+    """
+    Estimate the tensile strength of a slab layer from its density.
+
+    Uses the density parametrization of Adam (2025).
+
+    Arguments
+    ---------
+    rho : ndarray, float
+        Layer density (kg/m^3).
+    unit : str, optional
+        Desired output unit of the layer strength. Default is 'kPa'.
+
+    Returns
+    -------
+    ndarray
+        Tensile strength in specified unit.
+    """
+    convert = {"kPa": 1e3, "MPa": 1}
+    TS_0 = 1.0  # [MPa]
+    kappa = 3.45  # [-]
+    # Adam's equation is given in MPa
+    return TS_0 * (rho / RHO_ICE) ** kappa * convert[unit]
+
+
+# # TODO: Compressive Strength from Schöttner
+# def _schotter_compressive_strength(rho, unit: Literal["kPa", "MPa"] = "kPa"):
+#     """
+#     Estimate the compressive strength of a slab layer from its density.
+#     On the compressive strength of weak snow layers of depth hoar - Schöttner (2025).
+
+#     Uses the density parametrization of Schöttner (2025).
+#     """
+#     convert = {"kPa": 1e3, "MPa": 1}
+#     CS_0 = 11.0  # [MPa]
+#     CS_1 = 5.4  # [-]
+#     return CS_0 * (rho / RHO_ICE) ** CS_1 * convert[unit]
 
 
 class Layer(BaseModel):
@@ -114,6 +149,10 @@ class Layer(BaseModel):
         Young's modulus E [MPa].  If omitted it is derived from ``rho``.
     G : float, optional
         Shear modulus G [MPa].  If omitted it is derived from ``E`` and ``nu``.
+    tensile_strength: float
+        Tensile strength [kPa].
+    tensile_strength_method: Literal["sigrist", "adam", "hybrid"]
+        Method to calculate the tensile strength.
     """
 
     # has to be provided
@@ -129,8 +168,8 @@ class Layer(BaseModel):
     tensile_strength: float = Field(
         default=0.0, ge=0, description="Tensile strength [kPa]"
     )
-    tensile_strength_method: Literal["sigrist"] = Field(
-        default="sigrist",
+    tensile_strength_method: Literal["sigrist", "adam", "hybrid"] = Field(
+        default="hybrid",
         description="Method to calculate the tensile strength",
     )
     E_method: Literal["bergfeld", "scapazzo", "gerling"] = Field(
@@ -153,17 +192,23 @@ def model_post_init(self, _ctx):  # pylint: disable=arguments-differ
         else:
             raise ValueError(f"Invalid E_method: {self.E_method}")
         object.__setattr__(self, "G", self.G or self.E / (2 * (1 + self.nu)))
-        if self.tensile_strength_method == "sigrist":
-            object.__setattr__(
-                self,
-                "tensile_strength",
-                self.tensile_strength
-                or _sigrist_tensile_strength(self.rho, unit="kPa"),
-            )
-        else:
-            raise ValueError(
-                f"Invalid tensile_strength_method: {self.tensile_strength_method}"
-            )
+
+        if not self.tensile_strength:
+            if self.tensile_strength_method == "sigrist":
+                ts_value = _sigrist_tensile_strength(self.rho, unit="kPa")
+            elif self.tensile_strength_method == "adam":
+                ts_value = _adam_tensile_strength(self.rho, unit="kPa")
+            elif self.tensile_strength_method == "hybrid":
+                # Use Sigrist for rho < 250, Adam for rho >= 250
+                if self.rho < 250:
+                    ts_value = _sigrist_tensile_strength(self.rho, unit="kPa")
+                else:
+                    ts_value = _adam_tensile_strength(self.rho, unit="kPa")
+            else:
+                raise ValueError(
+                    f"Invalid tensile_strength_method: {self.tensile_strength_method}"
+                )
+            object.__setattr__(self, "tensile_strength", ts_value)
 
     @model_validator(mode="after")
     def validate_positive_E_G(self):

src/weac/core/slab_touchdown.py

@@ -20,6 +20,9 @@
 logger = logging.getLogger(__name__)
 
 
+TouchdownMode = Literal["A_free_hanging", "B_point_contact", "C_in_contact"]
+
+
 class SlabTouchdown:  # pylint: disable=too-many-instance-attributes,too-few-public-methods
     """
     Handling the touchdown situation in a PST.
@@ -56,7 +59,7 @@ class SlabTouchdown:  # pylint: disable=too-many-instance-attributes,too-few-pub
         Length of the crack for transition of stage A to stage B [mm]
     l_BC : float
         Length of the crack for transition of stage B to stage C [mm]
-    touchdown_mode : Literal["A_free_hanging", "B_point_contact", "C_in_contact"]
+    touchdown_mode : TouchdownMode
         Type of touchdown mode
     touchdown_distance : float
         Length of the touchdown segment [mm]
@@ -74,9 +77,7 @@ class SlabTouchdown:  # pylint: disable=too-many-instance-attributes,too-few-pub
     straight_scenario: Scenario
     l_AB: float
     l_BC: float
-    touchdown_mode: Literal[
-        "A_free_hanging", "B_point_contact", "C_in_contact"
-    ]  # Three types of contact with collapsed weak layer
+    touchdown_mode: TouchdownMode  # Three types of contact with collapsed weak layer
     touchdown_distance: float
     collapsed_weak_layer_kR: float | None = None
 

src/weac/core/system_model.py

@@ -335,8 +335,7 @@ def update_scenario(
             weak_layer=self.weak_layer,
             slab=self.slab,
         )
-        if self.config.touchdown:
-            self._invalidate_slab_touchdown()
+        self._invalidate_slab_touchdown()
         self._invalidate_constants()
 
     def toggle_touchdown(self, touchdown: bool):

tests/analysis/test_criteria_evaluator.py

@@ -272,6 +272,46 @@ def test_find_minimum_crack_length(self):
         self.assertIsInstance(new_segments, list)
         self.assertTrue(all(isinstance(s, Segment) for s in new_segments))
 
+    def test_evaluate_SteadyState_modes(self):
+        """Test evaluate_SteadyState with various modes."""
+        test_cases = [
+            ("C_in_contact", "C_in_contact"),
+            ("B_point_contact", "B_point_contact"),
+            ("A_free_hanging", "A_free_hanging"),
+            (None, "C_in_contact"),  # default mode
+        ]
+
+        for mode_param, expected_mode in test_cases:
+            with self.subTest(mode=mode_param):
+                segments = [
+                    Segment(length=self.segments_length, has_foundation=True, m=0),
+                    Segment(length=self.segments_length, has_foundation=True, m=0),
+                ]
+                system = SystemModel(
+                    model_input=ModelInput(
+                        layers=self.layers,
+                        weak_layer=self.weak_layer,
+                        segments=segments,
+                        scenario_config=ScenarioConfig(phi=self.phi),
+                    ),
+                    config=Config(touchdown=True),
+                )
+
+                if mode_param is None:
+                    results: SteadyStateResult = self.evaluator.evaluate_SteadyState(
+                        system
+                    )
+                else:
+                    results: SteadyStateResult = self.evaluator.evaluate_SteadyState(
+                        system, mode=mode_param
+                    )
+
+                self.assertTrue(results.converged)
+                self.assertEqual(
+                    results.system.slab_touchdown.touchdown_mode,
+                    expected_mode,
+                )
+
 
 if __name__ == "__main__":
     unittest.main()