diff --git a/README.md b/README.md index 6a09c00..4fc0825 100644 --- a/README.md +++ b/README.md @@ -7,6 +7,7 @@ A comprehensive Python package for slope stability analysis using the SHALSTAB ( [![CI](https://github.com/federicogmz/shalstab-python/actions/workflows/ci.yml/badge.svg)](https://github.com/federicogmz/shalstab-python/actions/workflows/ci.yml) [![codecov](https://codecov.io/gh/federicogmz/shalstab-python/branch/main/graph/badge.svg)](https://codecov.io/gh/federicogmz/shalstab-python) [![PyPI version](https://badge.fury.io/py/shalstab.svg)](https://pypi.org/project/shalstab/) +[![MCP Enabled](https://img.shields.io/badge/MCP-Enabled-blueviolet)](https://modelcontextprotocol.io) [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![Python 3.9+](https://img.shields.io/badge/python-3.9+-blue.svg)](https://www.python.org/downloads/) @@ -14,6 +15,7 @@ A comprehensive Python package for slope stability analysis using the SHALSTAB ( - [Overview](#overview) - [Installation](#installation) +- [Model Context Protocol (MCP) Integration](#model-context-protocol-mcp-integration) - [Quick Start](#quick-start) - [Scientific Background](#scientific-background) - [Mathematical Framework](#mathematical-framework) @@ -62,6 +64,76 @@ pip install -e . pip install -r requirements.txt ``` +## Model Context Protocol (MCP) Integration + +> **Note:** MCP support is included in the development branch (`v1.1.0.dev0`). For local development, install with `pip install -e ".[mcp]"`. + +`shalstab` natively supports the [Model Context Protocol (MCP)](https://modelcontextprotocol.io), allowing frontier AI agents (Claude Desktop, Claude Code, OpenAI Agents SDK) to execute deterministic slope stability physics directly from natural language. Note: MCP features require Python $\ge 3.10$. + +### Quickstart with Claude Desktop + +Add this to your `claude_desktop_config.json`: + +**Directly from GitHub (Recommended):** +```json +{ + "mcpServers": { + "shalstab": { + "command": "uvx", + "args": [ + "--python", "3.11", + "--from", "shalstab[mcp] @ git+https://github.com/federicogmz/shalstab-python.git", + "shalstab-mcp" + ] + } + } +} +``` + +**From local repository checkout:** +```json +{ + "mcpServers": { + "shalstab": { + "command": "uv", + "args": [ + "run", + "--python", "3.11", + "--directory", "/path/to/shalstab-python", + "--extra", "mcp", + "shalstab-mcp" + ] + } + } +} +``` + +### Quickstart with OpenAI Agents SDK + +*(Ensure `shalstab-mcp` is installed via `pip install -e ".[mcp]"` or from git `pip install "shalstab[mcp] @ git+https://github.com/federicogmz/shalstab-python.git"`, or run with `uvx`)*: + +```python +import asyncio +from agents import Agent, Runner +from agents.mcp import MCPServerStdio + +async def main(): + # Connect to the local SHALSTAB MCP server via standard I/O (requires Python >= 3.10) + async with MCPServerStdio(params={"command": "shalstab-mcp"}) as server: + agent = Agent(name="GeohazardsAgent", mcp_servers=[server]) + result = await Runner.run( + agent, "Evaluate slope stability for a 35 degree slope with 2 kPa cohesion and 32 degree friction angle" + ) + print(result.final_output) + +asyncio.run(main()) +``` + +### Exposed Agent Tools +- `compute_infinite_slope_stability`: Analytical Factor of Safety (FS) and 4-class stability classification under dry, steady-state, and fully saturated regimes. +- `diagnose_topographic_wetness`: Steady-state hydrologic saturation and critical recharge threshold $\log_{10}(q/T)$ (and natural $\ln(q/T)$) for shallow failure per Montgomery & Dietrich (1994). + *Note on formulation:* Analytical critical recharge strictly implements Montgomery & Dietrich (1994) Eq. 5, matching `Analyzer.calculate_log_qt`. + ## Training data The easiest way to test if SHALSTAB is installed correctly or to get familiar with its implementation is to use the provided training data: diff --git a/pyproject.toml b/pyproject.toml index 212dccd..43e523c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "shalstab" -version = "1.0.3" +version = "1.1.0.dev0" authors = [ {name = "Federico Gómez", email = "fjgomezc@eafit.edu.co"}, ] @@ -55,7 +55,12 @@ dependencies = [ ] [project.optional-dependencies] -test = ["pytest>=7.0.0", "pytest-cov>=4.0.0"] +test = [ + "pytest>=7.0.0", + "pytest-cov>=4.0.0", + "mcp>=1.2.0; python_version>='3.10'", +] +mcp = ["mcp>=1.2.0; python_version>='3.10'"] dev = [ "pytest>=7.0.0", "pytest-cov>=4.0.0", @@ -73,8 +78,12 @@ all = [ "sphinx-rtd-theme>=1.0.0", "jupyter>=1.0.0", "pre-commit>=2.20.0", + "mcp>=1.2.0; python_version>='3.10'", ] +[project.scripts] +shalstab-mcp = "shalstab.mcp_server:main" + [project.urls] Homepage = "https://github.com/federicogmz/shalstab-python" Repository = "https://github.com/federicogmz/shalstab-python.git" diff --git a/setup.py b/setup.py index 2e7e653..198909e 100644 --- a/setup.py +++ b/setup.py @@ -35,7 +35,7 @@ setup( name="shalstab", - version="1.0.0", + version="1.1.0.dev0", author="Federico Gómez", author_email="fjgomezc@eafit.edu.co", description="Shallow Landsliding STABility (SHALSTAB) model for slope stability analysis", diff --git a/shalstab/__init__.py b/shalstab/__init__.py index c4d04ca..0baa433 100644 --- a/shalstab/__init__.py +++ b/shalstab/__init__.py @@ -12,14 +12,39 @@ """ -from .slope_analyzer import Analyzer +import importlib.util +from typing import TYPE_CHECKING, Any + +if TYPE_CHECKING: + from .slope_analyzer import Analyzer + +__version__ = "1.1.0.dev0" +__all__ = ["Analyzer", "__version__"] -# Training data try: - from .training_data import training_dem, training_geology + if importlib.util.find_spec("shalstab.training_data") is not None: + __all__.extend(["training_dem", "training_geology"]) +except Exception: + pass + + +def __getattr__(name: str) -> Any: + """Lazy import heavy modules only when accessed to keep startup fast and lean.""" + if name == "Analyzer": + from .slope_analyzer import Analyzer + + return Analyzer + elif name in ("training_dem", "training_geology"): + try: + from . import training_data + + return getattr(training_data, name) + except ImportError: + raise AttributeError( + f"module {__name__!r} has no attribute {name!r}" + ) from None + raise AttributeError(f"module {__name__!r} has no attribute {name!r}") - __all__ = ["Analyzer", "training_dem", "training_geology"] -except ImportError: - __all__ = ["Analyzer"] -__version__ = "1.0.3" +def __dir__() -> list[str]: + return sorted(set(globals()) | set(__all__)) diff --git a/shalstab/mcp_server.py b/shalstab/mcp_server.py new file mode 100644 index 0000000..57d3b02 --- /dev/null +++ b/shalstab/mcp_server.py @@ -0,0 +1,394 @@ +""" +Model Context Protocol (MCP) server for SHALSTAB slope stability analysis. + +Exposes deterministic slope stability calculations and hydrological steady-state +diagnostic tools compatible with Claude Desktop, Claude Code, and OpenAI Agents SDK. +""" + +from __future__ import annotations + +import math +import sys +from typing import Any, Dict + +# Water unit weight constant matching Analyzer.GRAVITY_WATER +GRAVITY_WATER: float = 9.81 # kN/m³ + +# FastMCP / MCPServer import handled gracefully across mcp 1.x and 2.x +try: + from mcp.server.fastmcp import FastMCP +except ImportError: + try: + from mcp.server.mcpserver import MCPServer as FastMCP + except ImportError: + FastMCP = None + + +def compute_infinite_slope_stability( + slope_degrees: float, + cohesion_kpa: float, + friction_angle_degrees: float, + soil_unit_weight_kn_m3: float = 19.0, + soil_thickness_m: float = 1.5, + water_table_ratio: float = 0.0, +) -> Dict[str, Any]: + """ + Compute factor of safety and SHALSTAB stability classification. + + Evaluates Coulomb limit equilibrium combined with Montgomery & Dietrich (1994) + infinite slope mechanics across dry, saturated, and steady-state seepage. + + Note on stability codes: + - Code 1 (Unconditionally Stable): tan(theta) < stable threshold. The slope + cannot fail even under complete saturation (m=1.0). + - Code 2 (Unconditionally Unstable): tan(theta) >= dry strength threshold. + The slope is unstable even in dry conditions (m=0.0). + - Code 3 (Unstable): conditionally unstable under current water table ratio m. + - Code 4 (Stable): conditionally stable under current water table ratio m. + For unconditional categories (codes 1 and 2), the classification is an intrinsic + property of terrain slope and soil strength, independent of water table height. + + Parameters + ---------- + slope_degrees : float + Topographic slope angle theta in degrees (0 < slope_degrees < 90). + cohesion_kpa : float + Effective soil cohesion c' in kPa (cohesion_kpa >= 0). + friction_angle_degrees : float + Internal friction angle phi' in degrees (0 < friction_angle_degrees < 90). + soil_unit_weight_kn_m3 : float, default 19.0 + Total soil unit weight gamma in kN/m³ (must be > 9.81 kN/m³). + soil_thickness_m : float, default 1.5 + Vertical soil mantle thickness z in meters (soil_thickness_m > 0). + water_table_ratio : float, default 0.0 + Water table height ratio m = h/z (0.0 = dry, 1.0 = fully saturated). + + Returns + ------- + Dict[str, Any] + Dictionary containing factor_of_safety, stability_category, stability_code, + unconditional thresholds, physical diagnosis, and input parameters. + """ + if not (0.0 < slope_degrees < 90.0): + raise ValueError(f"slope_degrees must be in (0, 90), got {slope_degrees}") + if cohesion_kpa < 0.0: + raise ValueError(f"cohesion_kpa must be >= 0, got {cohesion_kpa}") + if not (0.0 < friction_angle_degrees < 90.0): + raise ValueError( + f"friction_angle_degrees must be in (0, 90), got {friction_angle_degrees}" + ) + if soil_unit_weight_kn_m3 <= GRAVITY_WATER: + raise ValueError( + f"soil_unit_weight_kn_m3 must be > water unit weight ({GRAVITY_WATER}), " + f"got {soil_unit_weight_kn_m3}" + ) + if soil_thickness_m <= 0.0: + raise ValueError(f"soil_thickness_m must be > 0, got {soil_thickness_m}") + if not (0.0 <= water_table_ratio <= 1.0): + raise ValueError( + f"water_table_ratio must be in [0, 1], got {water_table_ratio}" + ) + + theta = math.radians(slope_degrees) + phi = math.radians(friction_angle_degrees) + tan_slope = math.tan(theta) + tan_friction = math.tan(phi) + cos_slope = math.cos(theta) + sin_slope = math.sin(theta) + cos_slope_sq = cos_slope**2 + + # Cohesion resistance component + cohesion_term = cohesion_kpa / ( + soil_unit_weight_kn_m3 * soil_thickness_m * cos_slope_sq + ) + + # Unconditional stability thresholds matching Analyzer + unstable_threshold = tan_friction + cohesion_term + stable_threshold = ( + (1.0 - (GRAVITY_WATER / soil_unit_weight_kn_m3)) * tan_friction + ) + cohesion_term + + # Factor of Safety calculation (Skempton-Delory / Coulomb equation) + shear_strength = cohesion_kpa + ( + (soil_unit_weight_kn_m3 - water_table_ratio * GRAVITY_WATER) + * soil_thickness_m + * cos_slope_sq + * tan_friction + ) + shear_stress = soil_unit_weight_kn_m3 * soil_thickness_m * sin_slope * cos_slope + fs = shear_strength / shear_stress if shear_stress > 0.0 else float("inf") + + # SHALSTAB Classification + is_unconditionally_unstable = tan_slope >= unstable_threshold + is_unconditionally_stable = tan_slope < stable_threshold + + if is_unconditionally_unstable: + category = "Unconditionally Unstable" + code = 2 + dry_deg = math.degrees(math.atan(unstable_threshold)) + diagnosis = ( + f"Slope angle ({slope_degrees:.1f}°) exceeds dry failure threshold " + f"({dry_deg:.1f}°). Unstable even in dry conditions (FS = {fs:.2f})." + ) + elif is_unconditionally_stable: + category = "Unconditionally Stable" + code = 1 + sat_deg = math.degrees(math.atan(stable_threshold)) + diagnosis = ( + f"Slope angle ({slope_degrees:.1f}°) is gentler than saturated threshold " + f"({sat_deg:.1f}°). Stable even under full saturation (FS = {fs:.2f})." + ) + else: + if fs < 1.0: + category = "Unstable" + code = 3 + diagnosis = ( + f"Conditionally unstable at water table ratio " + f"m={water_table_ratio:.2f} (FS = {fs:.2f} < 1.0)." + ) + else: + category = "Stable" + code = 4 + diagnosis = ( + f"Conditionally stable at water table ratio " + f"m={water_table_ratio:.2f} (FS = {fs:.2f} >= 1.0)." + ) + + return { + "factor_of_safety": round(fs, 4), + "stability_category": category, + "stability_code": code, + "unconditional_stable_threshold_tan": round(stable_threshold, 4), + "unconditional_unstable_threshold_tan": round(unstable_threshold, 4), + "is_unconditionally_stable": is_unconditionally_stable, + "is_unconditionally_unstable": is_unconditionally_unstable, + "physical_diagnosis": diagnosis, + "parameters": { + "slope_degrees": slope_degrees, + "cohesion_kpa": cohesion_kpa, + "friction_angle_degrees": friction_angle_degrees, + "soil_unit_weight_kn_m3": soil_unit_weight_kn_m3, + "soil_thickness_m": soil_thickness_m, + "water_table_ratio": water_table_ratio, + }, + } + + +def diagnose_topographic_wetness( + specific_catchment_area_m: float, + slope_degrees: float, + transmissivity_m2_day: float, + rainfall_mm_day: float, + cohesion_kpa: float = 0.0, + friction_angle_degrees: float = 33.0, + soil_unit_weight_kn_m3: float = 19.0, + soil_thickness_m: float = 1.5, +) -> Dict[str, Any]: + """ + Diagnose topographic wetness index and critical rainfall failure threshold. + + Computes steady-state hydrologic saturation w = (q * a) / (T * sin(theta)) + and critical log(q/T) ratio per Montgomery & Dietrich (1994). + + Parameters + ---------- + specific_catchment_area_m : float + Specific contributing drainage area a = A / b in meters, where A is upslope + contributing area (m²) and b is contour width (m). + slope_degrees : float + Ground surface slope angle theta in degrees (0 < slope_degrees < 90). + transmissivity_m2_day : float + Soil profile transmissivity T = K * z * cos(theta) in m²/day. + rainfall_mm_day : float + Steady-state rainfall infiltration rate in mm/day. Converted internally + to recharge q in m/day (q = rainfall_mm_day / 1000.0). + cohesion_kpa : float, default 0.0 + Effective soil cohesion c' in kPa (cohesion_kpa >= 0.0). + friction_angle_degrees : float, default 33.0 + Internal friction angle phi' in degrees (0 < friction_angle_degrees < 90). + soil_unit_weight_kn_m3 : float, default 19.0 + Total soil unit weight gamma in kN/m³ (must be > 9.81 kN/m³). + soil_thickness_m : float, default 1.5 + Vertical soil depth z in meters (soil_thickness_m > 0). + + Returns + ------- + Dict[str, Any] + Hydrologic diagnosis including wetness ratio w, critical recharge ratios + ('log10_q_over_t_critical' per Montgomery & Dietrich 1994, + and 'ln_q_over_t_critical' matching Analyzer.calculate_log_qt), failure prediction. + """ + if specific_catchment_area_m <= 0.0: + raise ValueError( + f"specific_catchment_area_m must be > 0, got {specific_catchment_area_m}" + ) + if not (0.0 < slope_degrees < 90.0): + raise ValueError(f"slope_degrees must be in (0, 90), got {slope_degrees}") + if transmissivity_m2_day <= 0.0: + raise ValueError( + f"transmissivity_m2_day must be > 0, got {transmissivity_m2_day}" + ) + if rainfall_mm_day < 0.0: + raise ValueError(f"rainfall_mm_day must be >= 0, got {rainfall_mm_day}") + if cohesion_kpa < 0.0: + raise ValueError(f"cohesion_kpa must be >= 0, got {cohesion_kpa}") + if not (0.0 < friction_angle_degrees < 90.0): + raise ValueError( + f"friction_angle_degrees must be in (0, 90), got {friction_angle_degrees}" + ) + if soil_unit_weight_kn_m3 <= GRAVITY_WATER: + raise ValueError( + f"soil_unit_weight_kn_m3 must be > {GRAVITY_WATER}, " + f"got {soil_unit_weight_kn_m3}" + ) + if soil_thickness_m <= 0.0: + raise ValueError(f"soil_thickness_m must be > 0, got {soil_thickness_m}") + + theta = math.radians(slope_degrees) + phi = math.radians(friction_angle_degrees) + sin_slope = math.sin(theta) + cos_slope = math.cos(theta) + tan_slope = math.tan(theta) + tan_friction = math.tan(phi) + cos_slope_sq = cos_slope**2 + + # Recharge q in m/day + q_m_day = rainfall_mm_day / 1000.0 + + # Steady-state wetness ratio w = (q * a) / (T * sin(theta)) + w = (q_m_day * specific_catchment_area_m) / (transmissivity_m2_day * sin_slope) + w_effective = min(w, 1.0) + + # Check unconditional stability thresholds matching Analyzer + cohesion_term = cohesion_kpa / ( + soil_unit_weight_kn_m3 * soil_thickness_m * cos_slope_sq + ) + unstable_threshold = tan_friction + cohesion_term + stable_threshold = ( + (1.0 - (GRAVITY_WATER / soil_unit_weight_kn_m3)) * tan_friction + ) + cohesion_term + + is_unconditionally_unstable = tan_slope >= unstable_threshold + is_unconditionally_stable = tan_slope < stable_threshold + + # Critical recharge threshold (q/T)_crit (Montgomery & Dietrich, 1994) + term_friction = (soil_unit_weight_kn_m3 / GRAVITY_WATER) * ( + 1.0 - (tan_slope / tan_friction) + ) + term_cohesion = cohesion_kpa / ( + GRAVITY_WATER * soil_thickness_m * cos_slope_sq * tan_friction + ) + q_over_t_crit = (sin_slope / specific_catchment_area_m) * ( + term_friction + term_cohesion + ) + + if is_unconditionally_unstable: + log10_q_over_t_crit = None + ln_q_over_t_crit = None + failure_predicted = True + diagnosis = ( + "Slope is unconditionally unstable (tan(theta) exceeds dry strength). " + "Failure is predicted regardless of hydrological recharge." + ) + elif is_unconditionally_stable: + log10_q_over_t_crit = None + ln_q_over_t_crit = None + failure_predicted = False + diagnosis = ( + "Slope is unconditionally stable (tan(theta) is below threshold). " + "No rainfall amount can trigger shallow planar failure." + ) + else: + log10_q_over_t_crit = round(math.log10(q_over_t_crit), 4) + ln_q_over_t_crit = round(math.log(q_over_t_crit), 4) + q_over_t_current = q_m_day / transmissivity_m2_day + failure_predicted = q_over_t_current >= q_over_t_crit + if failure_predicted: + diagnosis = ( + f"Current recharge ratio q/T ({q_over_t_current:.6f}) exceeds " + f"critical threshold ({q_over_t_crit:.6f}). Failure predicted." + ) + else: + diagnosis = ( + f"Current recharge ratio q/T ({q_over_t_current:.6f}) is below " + f"critical threshold ({q_over_t_crit:.6f}). Stable under rainfall." + ) + + if w >= 1.0: + regime = "Overland flow / Full saturation (w >= 1.0)" + elif w > 0.5: + regime = "High pore pressure build-up (0.5 < w < 1.0)" + else: + regime = "Subsurface unsaturated regime (w <= 0.5)" + + return { + "wetness_ratio_w": round(w, 4), + "effective_saturation_ratio": round(w_effective, 4), + "saturation_percentage": round(w_effective * 100.0, 2), + "critical_recharge_ratio_q_over_t": ( + round(q_over_t_crit, 6) + if (not is_unconditionally_unstable and not is_unconditionally_stable) + else None + ), + "log10_q_over_t_critical": log10_q_over_t_crit, + "ln_q_over_t_critical": ln_q_over_t_crit, + "hydrologic_regime": regime, + "is_saturated": w >= 1.0, + "is_unconditionally_stable": is_unconditionally_stable, + "is_unconditionally_unstable": is_unconditionally_unstable, + "failure_predicted": failure_predicted, + "physical_diagnosis": diagnosis, + "parameters": { + "specific_catchment_area_m": specific_catchment_area_m, + "slope_degrees": slope_degrees, + "transmissivity_m2_day": transmissivity_m2_day, + "rainfall_mm_day": rainfall_mm_day, + "cohesion_kpa": cohesion_kpa, + "friction_angle_degrees": friction_angle_degrees, + "soil_unit_weight_kn_m3": soil_unit_weight_kn_m3, + "soil_thickness_m": soil_thickness_m, + }, + } + + +def create_mcp_server() -> Any: + """Create and configure the FastMCP server instance.""" + if FastMCP is None: + raise ImportError( + "The 'mcp' package is required to run the MCP server. " + "Install it with: pip install 'shalstab[mcp]'" + ) + + server = FastMCP( + "shalstab", + instructions=( + "Deterministic infinite slope stability and hydrological analysis " + "based on SHALSTAB (Montgomery & Dietrich, 1994). Evaluates " + "Factor of Safety and critical rainfall thresholds. " + "Requires Python >= 3.10." + ), + ) + + server.tool()(compute_infinite_slope_stability) + server.tool()(diagnose_topographic_wetness) + + return server + + +def main() -> None: + """Entry point for shalstab-mcp CLI command.""" + if FastMCP is None: + sys.stderr.write( + "Error: Model Context Protocol SDK ('mcp') is not installed " + "or current Python version is < 3.10.\n" + "The SHALSTAB MCP server requires Python >= 3.10. Install with:\n" + " pip install 'shalstab[mcp]'\n" + ) + sys.exit(1) + + server = create_mcp_server() + server.run() + + +if __name__ == "__main__": + main() diff --git a/tests/test_mcp.py b/tests/test_mcp.py new file mode 100644 index 0000000..5882833 --- /dev/null +++ b/tests/test_mcp.py @@ -0,0 +1,293 @@ +"""Tests for SHALSTAB Model Context Protocol (MCP) server analytical tools.""" + +import math + +import pytest + +from shalstab.mcp_server import ( + compute_infinite_slope_stability, + diagnose_topographic_wetness, +) + + +def test_dry_slope_infinite_factor_of_safety(): + """Verify dry slope factor of safety matches tan(phi)/tan(theta).""" + slope = 20.0 + phi = 35.0 + result = compute_infinite_slope_stability( + slope_degrees=slope, + cohesion_kpa=0.0, + friction_angle_degrees=phi, + water_table_ratio=0.0, + ) + expected_fs = math.tan(math.radians(phi)) / math.tan(math.radians(slope)) + assert math.isclose(result["factor_of_safety"], expected_fs, rel_tol=1e-3) + assert result["is_unconditionally_stable"] is False + assert result["is_unconditionally_unstable"] is False + assert result["stability_category"] == "Stable" + assert result["stability_code"] == 4 + + +def test_fully_saturated_factor_of_safety(): + """Verify saturated slope factor of safety matches analytical formula.""" + slope = 15.0 + phi = 32.0 + gamma = 20.0 + gamma_w = 9.81 + result = compute_infinite_slope_stability( + slope_degrees=slope, + cohesion_kpa=0.0, + friction_angle_degrees=phi, + soil_unit_weight_kn_m3=gamma, + water_table_ratio=1.0, + ) + expected_fs = (1.0 - (gamma_w / gamma)) * ( + math.tan(math.radians(phi)) / math.tan(math.radians(slope)) + ) + assert math.isclose(result["factor_of_safety"], expected_fs, rel_tol=1e-3) + + +def test_partially_saturated_intermediate_factor_of_safety(): + """Verify FS decreases monotonically and matches analytical equation at m=0.5.""" + slope = 18.0 + phi = 30.0 + gamma = 19.0 + cohesion = 2.0 + thickness = 1.5 + gamma_w = 9.81 + + res_dry = compute_infinite_slope_stability( + slope_degrees=slope, + cohesion_kpa=cohesion, + friction_angle_degrees=phi, + soil_unit_weight_kn_m3=gamma, + soil_thickness_m=thickness, + water_table_ratio=0.0, + ) + res_mid = compute_infinite_slope_stability( + slope_degrees=slope, + cohesion_kpa=cohesion, + friction_angle_degrees=phi, + soil_unit_weight_kn_m3=gamma, + soil_thickness_m=thickness, + water_table_ratio=0.5, + ) + res_sat = compute_infinite_slope_stability( + slope_degrees=slope, + cohesion_kpa=cohesion, + friction_angle_degrees=phi, + soil_unit_weight_kn_m3=gamma, + soil_thickness_m=thickness, + water_table_ratio=1.0, + ) + + assert ( + res_dry["factor_of_safety"] + > res_mid["factor_of_safety"] + > res_sat["factor_of_safety"] + ) + + # Exact analytical verification for m=0.5 + theta = math.radians(slope) + phi_rad = math.radians(phi) + cos_sq = math.cos(theta) ** 2 + shear_strength = cohesion + ( + (gamma - 0.5 * gamma_w) * thickness * cos_sq * math.tan(phi_rad) + ) + shear_stress = gamma * thickness * math.sin(theta) * math.cos(theta) + expected_mid_fs = shear_strength / shear_stress + assert math.isclose(res_mid["factor_of_safety"], expected_mid_fs, rel_tol=1e-3) + + +def test_unconditionally_unstable_slope(): + """Verify slope steeper than friction angle with zero cohesion is unstable.""" + result = compute_infinite_slope_stability( + slope_degrees=40.0, + cohesion_kpa=0.0, + friction_angle_degrees=30.0, + ) + assert result["is_unconditionally_unstable"] is True + assert result["stability_category"] == "Unconditionally Unstable" + assert result["stability_code"] == 2 + + +def test_unconditionally_stable_slope(): + """Verify gentle slope with high cohesion is stable even when saturated.""" + result = compute_infinite_slope_stability( + slope_degrees=5.0, + cohesion_kpa=25.0, + friction_angle_degrees=35.0, + water_table_ratio=1.0, + ) + assert result["is_unconditionally_stable"] is True + assert result["stability_category"] == "Unconditionally Stable" + assert result["stability_code"] == 1 + + +def test_invalid_parameters_raise_value_error(): + """Verify physical bounds validation and error messages.""" + with pytest.raises(ValueError, match="slope_degrees"): + compute_infinite_slope_stability( + slope_degrees=0.0, cohesion_kpa=0, friction_angle_degrees=30 + ) + with pytest.raises(ValueError, match="cohesion_kpa"): + compute_infinite_slope_stability( + slope_degrees=30.0, cohesion_kpa=-1, friction_angle_degrees=30 + ) + with pytest.raises(ValueError, match="friction_angle_degrees"): + invalid_phi = 95.0 + compute_infinite_slope_stability( + slope_degrees=30.0, cohesion_kpa=0, friction_angle_degrees=invalid_phi + ) + with pytest.raises(ValueError, match="soil_unit_weight_kn_m3"): + compute_infinite_slope_stability( + slope_degrees=30.0, + cohesion_kpa=0, + friction_angle_degrees=30, + soil_unit_weight_kn_m3=8.0, + ) + with pytest.raises(ValueError, match="water_table_ratio"): + compute_infinite_slope_stability( + slope_degrees=30.0, + cohesion_kpa=0, + friction_angle_degrees=30, + water_table_ratio=1.5, + ) + + +def test_topographic_wetness_diagnosis(): + """Verify steady-state hydrological diagnosis and critical threshold.""" + result = diagnose_topographic_wetness( + specific_catchment_area_m=25.0, + slope_degrees=25.0, + transmissivity_m2_day=10.0, + rainfall_mm_day=50.0, + cohesion_kpa=2.0, + friction_angle_degrees=32.0, + ) + assert result["wetness_ratio_w"] > 0 + assert 0 <= result["effective_saturation_ratio"] <= 1.0 + assert result["log10_q_over_t_critical"] is not None + assert result["ln_q_over_t_critical"] is not None + # Rigorous mathematical cross-validation between base-10 and natural log + assert math.isclose( + result["ln_q_over_t_critical"], + result["log10_q_over_t_critical"] * math.log(10), + rel_tol=1e-3, + ) + assert "cohesion_kpa" in result["parameters"] + assert "failure_predicted" in result + + +def test_conditional_slope_rain_induced_failure(): + """Verify conditional slope predicts failure when rainfall exceeds threshold.""" + result = diagnose_topographic_wetness( + specific_catchment_area_m=50.0, + slope_degrees=25.0, + transmissivity_m2_day=5.0, + rainfall_mm_day=150.0, + cohesion_kpa=1.0, + friction_angle_degrees=32.0, + ) + assert result["failure_predicted"] is True + assert result["log10_q_over_t_critical"] is not None + assert result["ln_q_over_t_critical"] is not None + + +def test_topographic_wetness_unconditional_instability(): + """Verify diagnostic when steep slope fails unconditionally without rain.""" + result = diagnose_topographic_wetness( + specific_catchment_area_m=10.0, + slope_degrees=45.0, + transmissivity_m2_day=10.0, + rainfall_mm_day=0.0, + cohesion_kpa=0.0, + friction_angle_degrees=30.0, + ) + assert result["is_unconditionally_unstable"] is True + assert result["critical_recharge_ratio_q_over_t"] is None + assert result["log10_q_over_t_critical"] is None + assert result["ln_q_over_t_critical"] is None + assert result["failure_predicted"] is True + + +def test_topographic_wetness_unconditional_stability(): + """Verify gentle slope cannot fail regardless of extreme rainfall rate.""" + result = diagnose_topographic_wetness( + specific_catchment_area_m=1000.0, + slope_degrees=5.0, + transmissivity_m2_day=1.0, + rainfall_mm_day=200.0, + cohesion_kpa=0.0, + friction_angle_degrees=33.0, + ) + assert result["is_unconditionally_stable"] is True + assert result["failure_predicted"] is False + assert result["log10_q_over_t_critical"] is None + assert result["ln_q_over_t_critical"] is None + assert "unconditionally stable" in result["physical_diagnosis"].lower() + + +def test_topographic_wetness_validation_errors(): + """Verify input validation on hydrological parameters.""" + with pytest.raises(ValueError, match="specific_catchment_area_m"): + diagnose_topographic_wetness( + specific_catchment_area_m=-5.0, + slope_degrees=20.0, + transmissivity_m2_day=10.0, + rainfall_mm_day=10.0, + ) + with pytest.raises(ValueError, match="transmissivity_m2_day"): + diagnose_topographic_wetness( + specific_catchment_area_m=10.0, + slope_degrees=20.0, + transmissivity_m2_day=-1.0, + rainfall_mm_day=10.0, + ) + + +def test_mcp_server_factory(): + """Verify FastMCP server instance creates tools properly if mcp is installed.""" + pytest.importorskip("mcp") + import asyncio + + from shalstab.mcp_server import create_mcp_server + + server = create_mcp_server() + assert server.name == "shalstab" + tools = asyncio.run(server.list_tools()) + tool_names = [t.name for t in tools] + assert "compute_infinite_slope_stability" in tool_names + assert "diagnose_topographic_wetness" in tool_names + + +def test_unconditional_thresholds_match_analytical_formulation(): + """Verify unconditional thresholds match Montgomery & Dietrich (1994) equations.""" + slope = 28.0 + c = 3.5 + phi = 34.0 + gamma = 19.5 + z = 1.2 + gw = 9.81 + + res = compute_infinite_slope_stability( + slope_degrees=slope, + cohesion_kpa=c, + friction_angle_degrees=phi, + soil_unit_weight_kn_m3=gamma, + soil_thickness_m=z, + ) + + theta = math.radians(slope) + phi_rad = math.radians(phi) + cos_sq = math.cos(theta) ** 2 + c_term = c / (gamma * z * cos_sq) + expected_unstable = math.tan(phi_rad) + c_term + expected_stable = ((1.0 - (gw / gamma)) * math.tan(phi_rad)) + c_term + + assert math.isclose( + res["unconditional_unstable_threshold_tan"], expected_unstable, rel_tol=1e-4 + ) + assert math.isclose( + res["unconditional_stable_threshold_tan"], expected_stable, rel_tol=1e-4 + )