From 17289f0fa68d2455aa4357fdf96b157b3d1171d3 Mon Sep 17 00:00:00 2001 From: Yannik Werner Date: Mon, 2 Feb 2026 17:41:29 +0100 Subject: [PATCH 1/4] new states for steady state --- src/weac/analysis/criteria_evaluator.py | 25 +++++- src/weac/components/layer.py | 3 + src/weac/core/slab_touchdown.py | 9 ++- src/weac/core/system_model.py | 3 +- tests/analysis/test_criteria_evaluator.py | 99 +++++++++++++++++++++++ 5 files changed, 129 insertions(+), 10 deletions(-) diff --git a/src/weac/analysis/criteria_evaluator.py b/src/weac/analysis/criteria_evaluator.py index 842160d..90de157 100644 --- a/src/weac/analysis/criteria_evaluator.py +++ b/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 diff --git a/src/weac/components/layer.py b/src/weac/components/layer.py index e944a94..4f47fd9 100644 --- a/src/weac/components/layer.py +++ b/src/weac/components/layer.py @@ -97,6 +97,9 @@ def _sigrist_tensile_strength(rho, unit: Literal["kPa", "MPa"] = "kPa"): # TODO: Compressive Strength from Schöttner # (11 +/- 7) * (rho/rho_0) ^ (5.4 +/- 0.5) +# TODO: tensile strength from Valle +# Extrapolated for higher densities (see Teams Chat) + class Layer(BaseModel): """ diff --git a/src/weac/core/slab_touchdown.py b/src/weac/core/slab_touchdown.py index 32aed9e..36fe09f 100644 --- a/src/weac/core/slab_touchdown.py +++ b/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 diff --git a/src/weac/core/system_model.py b/src/weac/core/system_model.py index 3aa1d37..c9a8426 100644 --- a/src/weac/core/system_model.py +++ b/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): diff --git a/tests/analysis/test_criteria_evaluator.py b/tests/analysis/test_criteria_evaluator.py index 7f1845c..89669a8 100644 --- a/tests/analysis/test_criteria_evaluator.py +++ b/tests/analysis/test_criteria_evaluator.py @@ -272,6 +272,105 @@ 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_mode_C_in_contact(self): + """Test evaluate_SteadyState with mode='C_in_contact'.""" + 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), + ) + results: SteadyStateResult = self.evaluator.evaluate_SteadyState( + system, mode="C_in_contact" + ) + self.assertTrue(results.converged) + self.assertEqual( + results.system.slab_touchdown.touchdown_mode, + "C_in_contact", + "Touchdown mode should match the mode parameter passed to evaluate_SteadyState", + ) + + def test_evaluate_SteadyState_mode_B_point_contact(self): + """Test evaluate_SteadyState with mode='B_point_contact'.""" + 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), + ) + results: SteadyStateResult = self.evaluator.evaluate_SteadyState( + system, mode="B_point_contact" + ) + self.assertTrue(results.converged) + self.assertEqual( + results.system.slab_touchdown.touchdown_mode, + "B_point_contact", + "Touchdown mode should match the mode parameter passed to evaluate_SteadyState", + ) + + def test_evaluate_SteadyState_mode_A_free_hanging(self): + """Test evaluate_SteadyState with mode='A_free_hanging'.""" + 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), + ) + results: SteadyStateResult = self.evaluator.evaluate_SteadyState( + system, mode="A_free_hanging" + ) + self.assertTrue(results.converged) + self.assertEqual( + results.system.slab_touchdown.touchdown_mode, + "A_free_hanging", + "Touchdown mode should match the mode parameter passed to evaluate_SteadyState", + ) + + def test_evaluate_SteadyState_default_mode(self): + """Test evaluate_SteadyState with default mode (should be 'C_in_contact').""" + 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), + ) + # Call without specifying mode - should default to 'C_in_contact' + results: SteadyStateResult = self.evaluator.evaluate_SteadyState(system) + self.assertTrue(results.converged) + self.assertEqual( + results.system.slab_touchdown.touchdown_mode, + "C_in_contact", + "Default touchdown mode should be 'C_in_contact'", + ) + if __name__ == "__main__": unittest.main() From d3b09ecfa7b9848ad892c9fc19d8f7304c904f0c Mon Sep 17 00:00:00 2001 From: Yannik Werner Date: Mon, 2 Feb 2026 17:53:46 +0100 Subject: [PATCH 2/4] Replace redundant code in steady state tests --- tests/analysis/test_criteria_evaluator.py | 129 ++++++---------------- 1 file changed, 35 insertions(+), 94 deletions(-) diff --git a/tests/analysis/test_criteria_evaluator.py b/tests/analysis/test_criteria_evaluator.py index 89669a8..b11d063 100644 --- a/tests/analysis/test_criteria_evaluator.py +++ b/tests/analysis/test_criteria_evaluator.py @@ -272,104 +272,45 @@ 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_mode_C_in_contact(self): - """Test evaluate_SteadyState with mode='C_in_contact'.""" - segments = [ - Segment(length=self.segments_length, has_foundation=True, m=0), - Segment(length=self.segments_length, has_foundation=True, m=0), + 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 ] - 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), - ) - results: SteadyStateResult = self.evaluator.evaluate_SteadyState( - system, mode="C_in_contact" - ) - self.assertTrue(results.converged) - self.assertEqual( - results.system.slab_touchdown.touchdown_mode, - "C_in_contact", - "Touchdown mode should match the mode parameter passed to evaluate_SteadyState", - ) - def test_evaluate_SteadyState_mode_B_point_contact(self): - """Test evaluate_SteadyState with mode='B_point_contact'.""" - 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), - ) - results: SteadyStateResult = self.evaluator.evaluate_SteadyState( - system, mode="B_point_contact" - ) - self.assertTrue(results.converged) - self.assertEqual( - results.system.slab_touchdown.touchdown_mode, - "B_point_contact", - "Touchdown mode should match the mode parameter passed to evaluate_SteadyState", - ) + 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), + ) - def test_evaluate_SteadyState_mode_A_free_hanging(self): - """Test evaluate_SteadyState with mode='A_free_hanging'.""" - 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), - ) - results: SteadyStateResult = self.evaluator.evaluate_SteadyState( - system, mode="A_free_hanging" - ) - self.assertTrue(results.converged) - self.assertEqual( - results.system.slab_touchdown.touchdown_mode, - "A_free_hanging", - "Touchdown mode should match the mode parameter passed to evaluate_SteadyState", - ) + if mode_param is None: + results: SteadyStateResult = self.evaluator.evaluate_SteadyState( + system + ) + else: + results: SteadyStateResult = self.evaluator.evaluate_SteadyState( + system, mode=mode_param + ) - def test_evaluate_SteadyState_default_mode(self): - """Test evaluate_SteadyState with default mode (should be 'C_in_contact').""" - 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), - ) - # Call without specifying mode - should default to 'C_in_contact' - results: SteadyStateResult = self.evaluator.evaluate_SteadyState(system) - self.assertTrue(results.converged) - self.assertEqual( - results.system.slab_touchdown.touchdown_mode, - "C_in_contact", - "Default touchdown mode should be 'C_in_contact'", - ) + self.assertTrue(results.converged) + self.assertEqual( + results.system.slab_touchdown.touchdown_mode, + expected_mode, + ) if __name__ == "__main__": From d75e6a6b37a6fb1256a0210631a9ec361b4504bd Mon Sep 17 00:00:00 2001 From: Yannik Werner Date: Mon, 2 Feb 2026 18:25:11 +0100 Subject: [PATCH 3/4] Tensile Strength Hybrid Approach according to Adam --- src/weac/components/layer.py | 76 +++++++++--- tests/components/test_layer.py | 220 +++++++++++++++++++++++++++++++++ 2 files changed, 279 insertions(+), 17 deletions(-) diff --git a/src/weac/components/layer.py b/src/weac/components/layer.py index 4f47fd9..cc335db 100644 --- a/src/weac/components/layer.py +++ b/src/weac/components/layer.py @@ -94,11 +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. -# TODO: tensile strength from Valle -# Extrapolated for higher densities (see Teams Chat) + 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): @@ -117,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 @@ -132,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( @@ -156,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): diff --git a/tests/components/test_layer.py b/tests/components/test_layer.py index c69243f..656f69f 100644 --- a/tests/components/test_layer.py +++ b/tests/components/test_layer.py @@ -12,9 +12,11 @@ from weac.components.layer import ( Layer, WeakLayer, + _adam_tensile_strength, _bergfeld_youngs_modulus, _gerling_youngs_modulus, _scapozza_youngs_modulus, + _sigrist_tensile_strength, ) from weac.constants import NU @@ -45,6 +47,224 @@ def test_gerling_calculation(self): self.assertGreater(E, 0, "Young's modulus should be positive") +class TestTensileStrengthCalculations(unittest.TestCase): + """Test tensile strength calculation functions.""" + + def test_sigrist_calculation_kPa(self): + """Test Sigrist tensile strength calculation in kPa.""" + # Test with typical snow density + ts = _sigrist_tensile_strength(rho=200.0, unit="kPa") + self.assertGreater(ts, 0, "Tensile strength should be positive") + self.assertTrue(np.isscalar(ts), "Result should be a scalar") + + # Test with different densities + ts_light = _sigrist_tensile_strength(rho=100.0, unit="kPa") + ts_heavy = _sigrist_tensile_strength(rho=400.0, unit="kPa") + self.assertLess(ts_light, ts_heavy, "Heavier snow should have higher strength") + + def test_sigrist_calculation_MPa(self): + """Test Sigrist tensile strength calculation in MPa.""" + ts_kPa = _sigrist_tensile_strength(rho=200.0, unit="kPa") + ts_MPa = _sigrist_tensile_strength(rho=200.0, unit="MPa") + self.assertAlmostEqual( + ts_kPa, ts_MPa * 1000, places=5, msg="Unit conversion should be correct" + ) + + def test_adam_calculation_kPa(self): + """Test Adam tensile strength calculation in kPa.""" + # Test with typical snow density + ts = _adam_tensile_strength(rho=300.0, unit="kPa") + self.assertGreater(ts, 0, "Tensile strength should be positive") + self.assertTrue(np.isscalar(ts), "Result should be a scalar") + + # Test with different densities + ts_light = _adam_tensile_strength(rho=150.0, unit="kPa") + ts_heavy = _adam_tensile_strength(rho=450.0, unit="kPa") + self.assertLess(ts_light, ts_heavy, "Heavier snow should have higher strength") + + def test_adam_calculation_MPa(self): + """Test Adam tensile strength calculation in MPa.""" + ts_kPa = _adam_tensile_strength(rho=300.0, unit="kPa") + ts_MPa = _adam_tensile_strength(rho=300.0, unit="MPa") + self.assertAlmostEqual( + ts_kPa, ts_MPa * 1000, places=5, msg="Unit conversion should be correct" + ) + + def test_sigrist_vs_adam_comparison(self): + """Compare Sigrist and Adam formulations at different densities.""" + # At low densities, compare the formulations + rho_low = 150.0 + ts_sigrist = _sigrist_tensile_strength(rho=rho_low, unit="kPa") + ts_adam = _adam_tensile_strength(rho=rho_low, unit="kPa") + # Both should give positive values + self.assertGreater(ts_sigrist, 0) + self.assertGreater(ts_adam, 0) + + # At high densities + rho_high = 400.0 + ts_sigrist_high = _sigrist_tensile_strength(rho=rho_high, unit="kPa") + ts_adam_high = _adam_tensile_strength(rho=rho_high, unit="kPa") + self.assertGreater(ts_sigrist_high, 0) + self.assertGreater(ts_adam_high, 0) + + +class TestLayerTensileStrength(unittest.TestCase): + """Test Layer class tensile strength functionality.""" + + def test_layer_default_tensile_strength_method(self): + """Test that default method is 'hybrid'.""" + layer = Layer(rho=200.0, h=100.0) + self.assertEqual( + layer.tensile_strength_method, + "hybrid", + "Default method should be 'hybrid'", + ) + self.assertGreater( + layer.tensile_strength, 0, "Tensile strength should be calculated" + ) + + def test_layer_sigrist_method(self): + """Test Layer with explicit Sigrist method.""" + layer = Layer(rho=200.0, h=100.0, tensile_strength_method="sigrist") + expected_ts = _sigrist_tensile_strength(rho=200.0, unit="kPa") + self.assertAlmostEqual( + layer.tensile_strength, + expected_ts, + places=5, + msg="Tensile strength should match Sigrist calculation", + ) + + def test_layer_adam_method(self): + """Test Layer with explicit Adam method.""" + layer = Layer(rho=300.0, h=100.0, tensile_strength_method="adam") + expected_ts = _adam_tensile_strength(rho=300.0, unit="kPa") + self.assertAlmostEqual( + layer.tensile_strength, + expected_ts, + places=5, + msg="Tensile strength should match Adam calculation", + ) + + def test_layer_hybrid_method_low_density(self): + """Test hybrid method uses Sigrist for density < 250.""" + rho = 200.0 # Below 250 threshold + layer = Layer(rho=rho, h=100.0, tensile_strength_method="hybrid") + expected_ts = _sigrist_tensile_strength(rho=rho, unit="kPa") + self.assertAlmostEqual( + layer.tensile_strength, + expected_ts, + places=5, + msg="Hybrid should use Sigrist for rho < 250", + ) + + def test_layer_hybrid_method_high_density(self): + """Test hybrid method uses Adam for density >= 250.""" + rho = 300.0 # Above 250 threshold + layer = Layer(rho=rho, h=100.0, tensile_strength_method="hybrid") + expected_ts = _adam_tensile_strength(rho=rho, unit="kPa") + self.assertAlmostEqual( + layer.tensile_strength, + expected_ts, + places=5, + msg="Hybrid should use Adam for rho >= 250", + ) + + def test_layer_hybrid_method_at_threshold(self): + """Test hybrid method behavior exactly at 250 kg/m³.""" + rho = 250.0 # Exactly at threshold + layer = Layer(rho=rho, h=100.0, tensile_strength_method="hybrid") + expected_ts = _adam_tensile_strength(rho=rho, unit="kPa") + self.assertAlmostEqual( + layer.tensile_strength, + expected_ts, + places=5, + msg="Hybrid should use Adam for rho = 250", + ) + + def test_layer_custom_tensile_strength(self): + """Test that custom tensile strength overrides calculation.""" + custom_ts = 50.0 + layer = Layer( + rho=200.0, + h=100.0, + tensile_strength=custom_ts, + tensile_strength_method="sigrist", + ) + self.assertEqual( + layer.tensile_strength, + custom_ts, + "Custom tensile strength should override calculation", + ) + + +class TestTensileStrengthPhysicalConsistency(unittest.TestCase): + """Test physical consistency of tensile strength calculations.""" + + def test_density_strength_relationship(self): + """Test that higher density leads to higher tensile strength.""" + layer_light = Layer(rho=150.0, h=100.0) + layer_heavy = Layer(rho=350.0, h=100.0) + + self.assertLess( + layer_light.tensile_strength, + layer_heavy.tensile_strength, + "Heavier snow should have higher tensile strength", + ) + + def test_hybrid_continuity_around_threshold(self): + """Test continuity of hybrid method around 250 kg/m³ threshold.""" + # Test just below threshold + layer_below = Layer(rho=249.0, h=100.0, tensile_strength_method="hybrid") + # Test just above threshold + layer_above = Layer(rho=251.0, h=100.0, tensile_strength_method="hybrid") + + # Both should have positive strength + self.assertGreater(layer_below.tensile_strength, 0) + self.assertGreater(layer_above.tensile_strength, 0) + + # Values should be reasonably close (within an order of magnitude) + # This is a loose check since the formulations differ + ratio = layer_above.tensile_strength / layer_below.tensile_strength + self.assertLess( + ratio, 10.0, "Strength shouldn't jump by more than 10x at threshold" + ) + self.assertGreater( + ratio, 0.1, "Strength shouldn't drop by more than 10x at threshold" + ) + + def test_all_methods_give_positive_strength(self): + """Test that all methods produce positive tensile strength.""" + rho_values = [100.0, 200.0, 300.0, 400.0] + methods = ["sigrist", "adam", "hybrid"] + + for rho in rho_values: + for method in methods: + layer = Layer(rho=rho, h=100.0, tensile_strength_method=method) + self.assertGreater( + layer.tensile_strength, + 0, + f"Method {method} with rho={rho} should give positive strength", + ) + + def test_tensile_strength_density_monotonicity(self): + """Test that tensile strength increases monotonically with density.""" + densities = [100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0] + methods = ["sigrist", "adam", "hybrid"] + + for method in methods: + strengths = [ + Layer(rho=rho, h=100.0, tensile_strength_method=method).tensile_strength + for rho in densities + ] + # Check that each strength is greater than the previous + for i in range(1, len(strengths)): + self.assertGreater( + strengths[i], + strengths[i - 1], + f"Strength should increase with density for {method} method", + ) + + class TestLayer(unittest.TestCase): """Test the Layer class functionality.""" From 1e25db691a8fac7bf96b6dd822ad18eafe98f97d Mon Sep 17 00:00:00 2001 From: Yannik Werner Date: Mon, 2 Feb 2026 18:31:14 +0100 Subject: [PATCH 4/4] v3.1.2 --- CITATION.cff | 2 +- demo/demo.ipynb | 2 +- pyproject.toml | 4 ++-- src/weac/__init__.py | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/CITATION.cff b/CITATION.cff index f9a2004..3dcf0f2 100644 --- a/CITATION.cff +++ b/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 diff --git a/demo/demo.ipynb b/demo/demo.ipynb index 2f55dff..38007da 100644 --- a/demo/demo.ipynb +++ b/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", diff --git a/pyproject.toml b/pyproject.toml index 32db2bd..cf35b60 100644 --- a/pyproject.toml +++ b/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" diff --git a/src/weac/__init__.py b/src/weac/__init__.py index de16c65..58ccd86 100644 --- a/src/weac/__init__.py +++ b/src/weac/__init__.py @@ -2,4 +2,4 @@ WEAC - Weak Layer Anticrack Nucleation Model """ -__version__ = "3.1.1" +__version__ = "3.1.2"