diff --git a/TODO.md b/TODO.md index e8fc2b3..333e2ca 100644 --- a/TODO.md +++ b/TODO.md @@ -2,12 +2,15 @@ ## Major +- [ ] Layer & Slab using pipelines from Mary-Kate + Attributes (value, calculated [bool], pipeline, uncertainty) +- [ ] Uncertainties propagation - [ ] Use Classes for Boundary Types - [ ] Automatically figure out type of system - [ ] Automatically set boundary conditions based on system ## Minor +- [ ] Swap to Pytest from Unittest - [ ] resolve fracture criterion also when lower than strength criterion - [ ] Florian CriterionEvaluator: clarify and fix damping behavior (find_minimum_force / evaluate_coupled_criterion) - Expected behavior diff --git a/src/weac/analysis/criteria_evaluator.py b/src/weac/analysis/criteria_evaluator.py index d271c17..de31697 100644 --- a/src/weac/analysis/criteria_evaluator.py +++ b/src/weac/analysis/criteria_evaluator.py @@ -721,6 +721,12 @@ def evaluate_SteadyState( ----------- system: SystemModel The system model. + mode: TouchdownMode, optional + Touchdown evaluation mode. The three supported modes are: + ``"C_in_contact"`` for a cut distance of ``2 * l_BC``, + ``"B_point_contact"`` for a cut distance just shorter than ``l_BC``, + and ``"A_free_hanging"`` for a cut distance just shorter than ``l_AB``. + Defaults to ``"C_in_contact"``. vertical: bool, optional Whether to evaluate the system in a vertical configuration. Defaults to False. @@ -1301,15 +1307,30 @@ def _calculate_maximal_stresses( ) max_principal_stress_norm = np.max(principal_stress_norm) max_Sxx_norm = np.max(Sxx_norm) - # Count height levels as failure-prone when tensile stress exceeds the layer - # strength or the slab density is below the configured weak-snow threshold. - # zmesh rho is t/mm^3, layer rho is kg/m^3 + # zmesh rows are ordered from slab top to bottom. Low-density levels only + # fail through downward growth from tensile failures above; when they do, + # they are excluded from the slab tensile criterion denominator. + # zmesh rho is t/mm^3, layer rho is kg/m^3. zmesh = analyzer.get_zmesh(dz=1) rho_kg_m3 = zmesh["rho"] * 1e12 tensile_exceeds = np.max(Sxx_norm, axis=1) > 1 low_density = rho_kg_m3 <= self.criteria_config.low_density_threshold_kg_m3 - height_level_prone_to_fail = tensile_exceeds | low_density - slab_tensile_criterion = np.mean(height_level_prone_to_fail) + height_level_prone_to_fail = np.zeros_like(tensile_exceeds, dtype=bool) + all_above_prone_to_fail = True + for index, is_low_density in enumerate(low_density): + if is_low_density: + height_level_prone_to_fail[index] = all_above_prone_to_fail + else: + height_level_prone_to_fail[index] = tensile_exceeds[index] + all_above_prone_to_fail &= height_level_prone_to_fail[index] + + low_density_prone_to_fail = low_density & height_level_prone_to_fail + load_bearing_levels = ~low_density_prone_to_fail + slab_tensile_criterion = ( + np.mean(height_level_prone_to_fail[load_bearing_levels]) + if np.any(load_bearing_levels) + else 1.0 + ) if print_call_stats: analyzer.print_call_stats( message="_calculate_maximal_stresses Call Statistics" diff --git a/src/weac/components/criteria_config.py b/src/weac/components/criteria_config.py index 3aa7374..d9cdfff 100644 --- a/src/weac/components/criteria_config.py +++ b/src/weac/components/criteria_config.py @@ -50,7 +50,8 @@ class CriteriaConfig(BaseModel): Order of magnitude for stress envelope. Default is 1.0. low_density_threshold_kg_m3 : float Slab density threshold in kg/m^3 below which a layer is treated as weak snow - and counted as prone to tensile failure in the slab tensile criterion. + and excluded from the slab tensile criterion percentage when broken through + directional growth from above. """ fn: float = Field( @@ -94,6 +95,6 @@ class CriteriaConfig(BaseModel): gt=0, description=( "Slab density threshold in kg/m^3 below which a layer is treated as weak " - "snow in the slab tensile criterion" + "snow and conditionally excluded from the slab tensile criterion percentage" ), ) diff --git a/src/weac/core/field_quantities.py b/src/weac/core/field_quantities.py index 15ed0b0..8203b76 100644 --- a/src/weac/core/field_quantities.py +++ b/src/weac/core/field_quantities.py @@ -107,13 +107,14 @@ def sig(self, Z: np.ndarray, unit: StressUnit = "MPa") -> float | np.ndarray: return -self._unit_factor(unit) * self.es.weak_layer.kn * self.w(Z) def tau(self, Z: np.ndarray, unit: StressUnit = "MPa") -> float | np.ndarray: - """Weak-layer shear stress `tau = -kt * (w' * h/2 - u(h=H/2))`""" + """Weak-layer shear stress `tau = kt * h * (w' / 2 - u(h=H/2) / h)`""" return ( - -self._unit_factor(unit) + self._unit_factor(unit) * self.es.weak_layer.kt + * self.es.weak_layer.h * ( - self.dw_dx(Z) * self.es.weak_layer.h / 2 - - self.u(Z, h0=self.es.slab.H / 2) + self.dw_dx(Z) / 2 + - self.u(Z, h0=self.es.slab.H / 2) / self.es.weak_layer.h ) ) @@ -122,7 +123,7 @@ def eps(self, Z: np.ndarray) -> float | np.ndarray: return -self.w(Z) / self.es.weak_layer.h def gamma(self, Z: np.ndarray) -> float | np.ndarray: - """Weak-layer shear strain `gamma = (w' * h/2 - u(h=H/2)) / h`""" + """Weak-layer shear strain `gamma = w' / 2 - u(h=H/2) / h`""" return ( self.dw_dx(Z) / 2 - self.u(Z, h0=self.es.slab.H / 2) / self.es.weak_layer.h ) diff --git a/tests/analysis/test_analyzer.py b/tests/analysis/test_analyzer.py index 8c039df..9c2014a 100644 --- a/tests/analysis/test_analyzer.py +++ b/tests/analysis/test_analyzer.py @@ -156,6 +156,37 @@ def test_energy_release_rates_shapes(self): self.assertEqual(Gdif.shape, (4,)) self.assertTrue(np.isfinite(Gdif).all()) + def test_energy_release_rate_integrands_non_negative(self): + """Test that ERR integrands are non-negative for matching stress/strain.""" + slope_angle = 20.0 + system = SystemModel( + model_input=ModelInput( + scenario_config=ScenarioConfig(phi=slope_angle, system_type="skier"), + layers=[Layer()], + weak_layer=WeakLayer(), + segments=[Segment(), Segment()], + ), + config=Config(), + ) + analyzer = Analyzer(system_model=system, printing_enabled=False) + + z_uncracked = np.array([[0.0], [0.0], [1.0], [0.2], [0.0], [0.0]]) + + def constant_solution(x): + return np.repeat(z_uncracked, np.size(np.atleast_1d(x)), axis=1) + + mode_i = analyzer._integrand_GI( # pylint: disable=protected-access + np.array([0.0, 1.0]), constant_solution, constant_solution + ) + mode_ii = analyzer._integrand_GII( # pylint: disable=protected-access + np.array([0.0, 1.0]), constant_solution, constant_solution + ) + + self.assertTrue(np.all(mode_i >= 0), "Mode I integrand should be non-negative") + self.assertTrue( + np.all(mode_ii >= 0), "Mode II integrand should be non-negative" + ) + def test_internal_and_external_potentials_pst(self): """Test internal and external potentials for PST.""" # Ensure PST-specific methods run diff --git a/tests/analysis/test_criteria_evaluator.py b/tests/analysis/test_criteria_evaluator.py index d0b5dd9..462b0a0 100644 --- a/tests/analysis/test_criteria_evaluator.py +++ b/tests/analysis/test_criteria_evaluator.py @@ -68,19 +68,19 @@ def test_stress_envelope_adam_unpublished(self): self.assertGreater(result[0], 0) @patch("weac.analysis.criteria_evaluator.Analyzer") - def test_calculate_maximal_stresses_uses_configured_low_density_threshold( + def test_calculate_maximal_stresses_applies_directional_low_density_exclusion( self, mock_analyzer_cls - ): - """Test that the slab tensile criterion uses the configured density cutoff.""" - sxx_kpa = np.zeros((3, 1)) - principal_stress_kPa = np.zeros((3, 1)) - sxx_norm = np.full((3, 1), 0.5) - principal_stress_norm = np.full((3, 1), 0.5) + ): # pylint: disable=protected-access + """Test that weak snow is excluded only after downward failure growth.""" + sxx_kpa = np.zeros((4, 1)) + principal_stress_kPa = np.zeros((4, 1)) + sxx_norm = np.array([[1.5], [0.5], [0.5], [0.5]]) + principal_stress_norm = np.full((4, 1), 0.5) mock_analyzer = mock_analyzer_cls.return_value mock_analyzer.rasterize_solution.return_value = ( None, - np.array([0, 1, 2]), + np.array([0, 1, 2, 3]), None, ) mock_analyzer.Sxx.side_effect = ( @@ -92,20 +92,22 @@ def test_calculate_maximal_stresses_uses_configured_low_density_threshold( ) ) mock_analyzer.get_zmesh.return_value = { - "rho": np.array([90.0, 110.0, 130.0]) * 1e-12 + "rho": np.array([130.0, 90.0, 90.0, 130.0]) * 1e-12 } system = SimpleNamespace(scenario=SimpleNamespace(phi=30.0)) # Access the helper directly so the test can isolate the density-threshold logic. - default_result = CriteriaEvaluator( + top_broken_result = CriteriaEvaluator( CriteriaConfig() )._calculate_maximal_stresses(system=system) # pylint: disable=protected-access - tuned_result = CriteriaEvaluator( - CriteriaConfig(low_density_threshold_kg_m3=120) + + sxx_norm = np.array([[0.5], [0.5], [0.5], [1.5]]) + top_unbroken_result = CriteriaEvaluator( + CriteriaConfig() )._calculate_maximal_stresses(system=system) # pylint: disable=protected-access - self.assertAlmostEqual(default_result.slab_tensile_criterion, 1 / 3) - self.assertAlmostEqual(tuned_result.slab_tensile_criterion, 2 / 3) + self.assertAlmostEqual(top_broken_result.slab_tensile_criterion, 1 / 2) + self.assertAlmostEqual(top_unbroken_result.slab_tensile_criterion, 1 / 4) def test_find_minimum_force_convergence(self): """Test the convergence of find_minimum_force.""" diff --git a/tests/analysis/test_slab_tensile_comparisons.py b/tests/analysis/test_slab_tensile_comparisons.py index 9d8282f..4188026 100644 --- a/tests/analysis/test_slab_tensile_comparisons.py +++ b/tests/analysis/test_slab_tensile_comparisons.py @@ -80,8 +80,8 @@ def _setup_from_cm(*layers: tuple[float, float]) -> SetupDefinition: ), ComparisonCase( name="case_2", - setup_a=_setup_from_cm((50, 75), (20, 225)), - setup_b=_setup_from_cm((30, 75), (20, 225)), + setup_a=_setup_from_cm((30, 75), (20, 225)), + setup_b=_setup_from_cm((50, 75), (20, 225)), ), ComparisonCase( name="case_3", @@ -115,8 +115,13 @@ def _setup_from_cm(*layers: tuple[float, float]) -> SetupDefinition: ), ComparisonCase( name="case_9", - setup_a=_setup_from_cm((15, 275), (40, 75)), - setup_b=_setup_from_cm((40, 75), (15, 275)), + setup_a=_setup_from_cm((40, 75), (15, 275)), + setup_b=_setup_from_cm((15, 275), (40, 75)), + ), + ComparisonCase( + name="case_10", + setup_a=_setup_from_cm((30, 75), (20, 275)), + setup_b=_setup_from_cm((50, 75), (20, 275)), ), ) diff --git a/tests/core/test_field_quantities.py b/tests/core/test_field_quantities.py index e2d168a..3363a89 100644 --- a/tests/core/test_field_quantities.py +++ b/tests/core/test_field_quantities.py @@ -272,7 +272,7 @@ def test_weak_layer_shear_stress(self): H = self.fq.es.slab.H u_surface = self.fq.u(self.Z, h0=H / 2) - expected = -self.fq.es.weak_layer.kt * (self.Z[3, :] * h / 2 - u_surface) + expected = self.fq.es.weak_layer.kt * (self.Z[3, :] * h / 2 - u_surface) np.testing.assert_array_almost_equal( tau, expected,